/* part of the shansyn spherical harmonics package, see COPYRIGHT for license */
/* $Id: cmodelmeancorr.c,v 1.5 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, the fields are however 
// first averaged over a certain depth interval
//

int main(int argc, char **argv)
{
 
  int lmax,lmsize,l,steps;
  COMP_PRECISION dz,zmin,zmax,z,z1,z2,
    *am,*bm,rms[2],deltaz;
  struct mod model[2];
  zmin=0;
  zmax=2900;
  deltaz=500;
  dz=50;
  switch(argc){
  case 3:{
    break;
  }
  case 4:{
    sscanf(argv[3],DATA_FSCAN_FORMAT,&zmin);
    break;
  }
  case 5:{
    sscanf(argv[3],DATA_FSCAN_FORMAT,&zmin);
    sscanf(argv[4],DATA_FSCAN_FORMAT,&zmax);
    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,&deltaz);
    break;
  }
  default:{
    fprintf(stderr,"%s file1 file2 [zmin, %g] [zmax, %g] [dz, %g] [delta z, %g]\n",
	    argv[0],zmin,zmax,dz,deltaz);
    fprintf(stderr,"calculates correlation between fields\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 r_{l_max=8} r_{l_max=20} r_tot^w r_{l=1} r_{l=2} ... r_{l=lmax} \n\n");
    fprintf(stderr,"the depth will move from zmin (%g) to zmax (%g) in dz (%g) steps,\n",
	    zmin,zmax,dz);
    fprintf(stderr,"while fields will be average from z-deltaz/2 to z+deltaz/2, deltaz/2: %g\n",
	    deltaz/2.0);
    exit(-1);
    break;
  }}
  fprintf(stderr,"%s: %g - %g - %g, using delta z: %g\n",
	  argv[0],zmin,dz,zmax,deltaz);
  if(deltaz<0){
    fprintf(stderr,"%s: delta z should be >= 0, %g\n",
	    argv[0],deltaz);
    exit(-1);
  }
  // read in models
  read_she_model(argv[1],(model),-1,1);
  read_she_model(argv[2],(model+1),-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);
  }

  steps = (int)deltaz/20;
  deltaz /= 2.0;
  // depth loop
  for(z=zmin;z<=zmax+1.0e-5;z+=dz){
    z1=z-deltaz;
    if(z1<zmin)
      z1=zmin;
    z2=z+deltaz;
    if(z2>zmax)
      z2=zmax;

    mean_expansions(lmax,&am,&bm,2,model,z1, z2,steps,TRUE);
    rms[0]=calc_rms(am,bm,lmax);
    rms[1]=calc_rms((am+lmsize),(bm+lmsize),lmax);
    fprintf(stdout,"%12f %12.5e %12.5e %12.5e %12.5e %12.5e %12.5e ",
	    z1+deltaz/2.,rms[0],rms[1],
	    correlation(am,bm,(am+lmsize),(bm+lmsize),-lmax,0,1),
	    correlation(am,bm,(am+lmsize),(bm+lmsize),-((lmax>=8)?(8):(lmax)),0,1),
	    correlation(am,bm,(am+lmsize),(bm+lmsize),-((lmax>=20)?(20):(lmax)),0,1),
	    weighted_correlation(am,bm,(am+lmsize),(bm+lmsize),-lmax,1));
    for(l=1;l<=lmax;l++)
      fprintf(stdout,"%12.5e ",
	      correlation(am,bm,(am+lmsize),(bm+lmsize),l,0,1));
    fprintf(stdout,"\n");
  }
  return 0;
}
