/* part of the shansyn spherical harmonics package, see COPYRIGHT for license */
/* $Id: cmodelcorr.c,v 1.14 2005/04/05 19:06:14 becker Exp becker $ */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include "function_macros.h"
#include "trig_constants.h"
#include "precision.h"
#include "spherical_harmonics_functions.h"
#include "legendre_macros.h"
#include "myio.h"

//
// calculates correlation coefficients for each l
// for two different spherical harmonic expansion
// model files at various depths
//
// NOTE: there is also cmodelmeancorr for boxcar mean correlations
//
int main(int argc, char **argv)
{
 
  int i,lmax,lmsize,l,llim;
  COMP_PRECISION dz,zmin,zmax,z,*a,*b,rms[2],
    shift,stretch;
  struct mod model[2];
  // boundaries
  zmin=50.0;
  zmax=2850;
  dz=50.0;
  // shift of model 2 with depth
  shift=0.0;
  // stretch of depth layers of model 2 
  stretch=1.0;
  // llim, if set to -1, will use all coefficients, ie. l_max
  llim=-1;
  switch(argc){
  case 3:{
    break;
  }
  case 6:{
    sscanf(argv[3],DATA_FSCAN_FORMAT,&zmin);
    sscanf(argv[4],DATA_FSCAN_FORMAT,&zmax);
    sscanf(argv[5],DATA_FSCAN_FORMAT,&dz);
    break;
  }
  case 7:{
    sscanf(argv[3],DATA_FSCAN_FORMAT,&zmin);
    sscanf(argv[4],DATA_FSCAN_FORMAT,&zmax);
    sscanf(argv[5],DATA_FSCAN_FORMAT,&dz);
    sscanf(argv[6],DATA_FSCAN_FORMAT,&shift);
    break;
  }
  case 8:{
    sscanf(argv[3],DATA_FSCAN_FORMAT,&zmin);
    sscanf(argv[4],DATA_FSCAN_FORMAT,&zmax);
    sscanf(argv[5],DATA_FSCAN_FORMAT,&dz);
    sscanf(argv[6],DATA_FSCAN_FORMAT,&shift);
    sscanf(argv[7],DATA_FSCAN_FORMAT,&stretch);
    break;
  }
  case 9:{
    sscanf(argv[3],DATA_FSCAN_FORMAT,&zmin);
    sscanf(argv[4],DATA_FSCAN_FORMAT,&zmax);
    sscanf(argv[5],DATA_FSCAN_FORMAT,&dz);
    sscanf(argv[6],DATA_FSCAN_FORMAT,&shift);
    sscanf(argv[7],DATA_FSCAN_FORMAT,&stretch);
    sscanf(argv[8],"%i",&llim);
    break;
  }
  default:{
    fprintf(stderr,"%s file1 file2 [zmin(%g) zmax(%g) dz(%g)] [shift, %g] [stretch, %g] [L, %i]\n",
	    argv[0],zmin,zmax,dz,shift,stretch,llim);
    fprintf(stderr,"calculates correlation at depth intervals dz from zmin to zmax\n");
    fprintf(stderr,"based on two spherical harmonic model sets in file1 and file2\n");
    fprintf(stderr,"output is:\n\n depth RMS_mod_1 RMS_mod_2 r_{tot/L} r_{l_max=8} r_{l_max=20} r_tot^w r_{l=1} r_{l=2} ... r_{l=lmax} \n\n");
    fprintf(stderr,"if shift is set to non-zero, will shift all z's for model 2 by shift in depth\n");
    fprintf(stderr,"if stretch is set to non unity, will stretch depth levels of model 2 by strech\n");
    fprintf(stderr,"if L is set to a positive number, will limit r_{tot} to r_L, else r_{tot} will be up to the l_max of the models\n");
    exit(-1);
    break;
  }}
  if(dz<=0){ 
    dz=50;
    fprintf(stderr,"%s: adjusting dz to %g\n",argv[0],dz);
  }
  // read in models
  for(i=0;i<2;i++)
    read_she_model(argv[1+i],(model+i),-1,1);
  // pick smaller lmax 
  lmax=(model[0].lmax > model[1].lmax)?
    (model[1].lmax):(model[0].lmax);
  lmsize=(int)((((COMP_PRECISION)lmax)+1.0)*
	       (((COMP_PRECISION)lmax)+2)/2.0);
  fprintf(stderr,"%s: lmax models: %i/%i, using min lmax=%i\n",
	  argv[0],model[0].lmax,model[1].lmax,lmax);
  if(lmsize<1){
    fprintf(stderr,"%s: error: lmsize is %i\n",
	    argv[0],lmsize);
    exit(-1);
  }
  // interpolated coefficients for both models
  a=(COMP_PRECISION *)malloc(sizeof(COMP_PRECISION)*lmsize*2);
  b=(COMP_PRECISION *)malloc(sizeof(COMP_PRECISION)*lmsize*2);
  if(!a || !b)MEMERROR;

  // loop through all depths
  if(shift != 0.0){
    fprintf(stderr,"%s: WARNING: shift is non-zero, but %g\n",
	    argv[0],shift);
    
  }
  // stretch depth layers
  if(stretch != 1.0){
    fprintf(stderr,"%s: WARNING: stretch is non-unity, but %g\n",
	    argv[0],stretch);
    if(stretch < 1.0){
      zmax = model[1].d[0] * stretch;
      fprintf(stderr,"%s: limiting zmax to %g\n",
	      argv[0],zmax);
    }
    for(i=0;i<model[1].n;i++)
      model[1].d[i] *= stretch;
  }
  if(llim>=1)
    fprintf(stderr,"%s: using l_max=%i for r_{tot}, not %i\n",
	    argv[0],llim,lmax);
  for(z=zmin;z<=zmax+EPS_COMP_PREC;z+=dz){
    /* no extrapolation */
    // model A
    interpolate_she_model(a,b,model,z,lmax,FALSE);
    rms[0]=calc_rms(a,b,lmax);
    // model B
    interpolate_she_model((a+lmsize),(b+lmsize),(model+1),z+shift,lmax,FALSE);
    rms[1]=calc_rms((a+lmsize),(b+lmsize),lmax);
    fprintf(stdout,"%12f %12.5e %12.5e %12.5e %12.5e %12.5e %12.5e ",
	    z,rms[0],rms[1],
	    correlation(a,b,(a+lmsize),(b+lmsize),-((llim >=   1)?(llim):(lmax)),0,1),
	    correlation(a,b,(a+lmsize),(b+lmsize),-((lmax >=   8)?(8):(lmax)),0,1),
	    correlation(a,b,(a+lmsize),(b+lmsize),-((lmax >=  20)?(20):(lmax)),0,1),
	    weighted_correlation(a,b,(a+lmsize),(b+lmsize),-lmax,1));
    for(l=1;l<=lmax;l++)
      fprintf(stdout,"%12.5e ",
	      correlation(a,b,(a+lmsize),(b+lmsize),l,0,1));
    fprintf(stdout,"\n");
  }
  return 0;
}



