/* part of the shansyn spherical harmonics package, see COPYRIGHT for license */
/* $Id: cradialcorr.c,v 1.6 2009/04/16 00:55:29 becker Exp $ */
#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 the radial correlation function
//
// there is also a newer program called calc_radial_corr
//
// $Id: cradialcorr.c,v 1.6 2009/04/16 00:55:29 becker Exp $
//

int main(int argc, char **argv)
{
 
  int mode,i,j,lmsize,cmode;
  COMP_PRECISION *a1,*a2,*b1,*b2,zmin,zmax,dz,z1,z2,tmp[3],rdist,dmode,z,r,dr,a;
  struct mod model;
  zmin=0;
  zmax=2850;
  dz=50.0;
  dmode = 0;			/* will be mode */
  cmode = 1;
  switch(argc){
  case 2:{
    break;
  }
  case 3:{
    sscanf(argv[2],"%lf",&dmode);
    break;
  }
  case 6:{
    sscanf(argv[2],"%lf",&dmode);
    sscanf(argv[3],DATA_FSCAN_FORMAT,&zmin);
    sscanf(argv[4],DATA_FSCAN_FORMAT,&zmax);
    sscanf(argv[5],DATA_FSCAN_FORMAT,&dz);
    break;
  }
  default:{
    fprintf(stderr,"%s file1 [mode, %i] [zmin(%g) zmax(%g) dz(%g)]\n",
	    argv[0],(int)dmode,zmin,zmax,dz);
    fprintf(stderr,"calculates the radial correlation function of model file1\n");
    fprintf(stderr,"output is:\nz_1 z_2 r_8 r_20 r_total\n");
    fprintf(stderr,"mode 0: uses only the input model depths\n");
    fprintf(stderr,"mode 1: uses z values from zmin to zmax in dz steps\n");
    fprintf(stderr,"mode<0: find distance diagonal to correlation of -mode\n"); /* this is actually mode 2 internally */
    exit(-1);
    break;
  }}
  /* deal with mode */
  if(dmode < 0){
    rdist = -dmode;
    fprintf(stderr,"%s: finding distance from diagonal to correlation %g\n",argv[0],rdist);
    mode = 2;
  }else{
    mode=(int)dmode;
  }
  read_she_model(argv[1],&model,-1,1);
  switch(mode){
  case 0:{			/* compute matrix at model depth levels */
    if(model.type != DISCRETE){
      fprintf(stderr,"%s: mode 0 only works for discrete input level models\n",
	      argv[0]);
      exit(-1);
    }
    for(i=0;i<model.n;i++)
      for(j=i;j<model.n;j++){
	if(i==j){
	  fprintf(stdout,"%g %g %g %g %g\n",model.d[i],model.d[j],1.0,1.0,1.0);
	}else{
	  tmp[0]=correlation(model.a[i],model.b[i],model.a[j],model.b[j],-((model.lmax>=8)?(8):(model.lmax)),0,1,cmode);
	  tmp[1]=correlation(model.a[i],model.b[i],model.a[j],model.b[j],-((model.lmax>=20)?(20):(model.lmax)),0,1,cmode);
	  tmp[2]=correlation(model.a[i],model.b[i],model.a[j],model.b[j],-model.lmax,0,1,cmode);
	  fprintf(stdout,"%g %g %g %g %g\n",model.d[i],model.d[j],tmp[0],tmp[1],tmp[2]);
	  fprintf(stdout,"%g %g %g %g %g\n",model.d[j],model.d[i],tmp[0],tmp[1],tmp[2]);
	}
      }
    break;
  }
  case 1:{
    /* compute matrix at interpolated values */
    lmsize=(int)((((COMP_PRECISION)model.lmax)+1.0)*(((COMP_PRECISION)model.lmax)+2)/2.0);
    a1=(COMP_PRECISION *)malloc(sizeof(COMP_PRECISION)*lmsize);
    b1=(COMP_PRECISION *)malloc(sizeof(COMP_PRECISION)*lmsize);
    a2=(COMP_PRECISION *)malloc(sizeof(COMP_PRECISION)*lmsize);
    b2=(COMP_PRECISION *)malloc(sizeof(COMP_PRECISION)*lmsize);
    if(!a1 || !a2 || !b1 || !b2)MEMERROR;
    for(z1=zmin;z1<=zmax+EPS_COMP_PREC;z1+=dz)
      for(z2=z1;z2<=zmax+EPS_COMP_PREC;z2+=dz)
	if(z1==z2){
	  fprintf(stdout,"%g %g %g %g %g\n",z1,z2,1.0,1.0,1.0);
	}else{
	  interpolate_she_model(a1,b1,&model,z1,model.lmax,FALSE);
	  interpolate_she_model(a2,b2,&model,z2,model.lmax,FALSE);
	  tmp[0]=correlation(a1,b1,a2,b2,-((model.lmax>=8)?(8):(model.lmax)),0,1,cmode);
	  tmp[1]=correlation(a1,b1,a2,b2,-((model.lmax>=20)?(20):(model.lmax)),0,1,cmode);
	  tmp[2]=correlation(a1,b1,a2,b2,-model.lmax,0,1,cmode);
	  fprintf(stdout,"%g %g %g %g %g\n",z1,z2,tmp[0],tmp[1],tmp[2]);
	  fprintf(stdout,"%g %g %g %g %g\n",z2,z1,tmp[0],tmp[1],tmp[2]);
	}
    break;
  } 
  case 2:{
    /* find distance to a cretain correlation */
    lmsize=(int)((((COMP_PRECISION)model.lmax)+1.0)*(((COMP_PRECISION)model.lmax)+2)/2.0);
    a1=(COMP_PRECISION *)malloc(sizeof(COMP_PRECISION)*lmsize);
    b1=(COMP_PRECISION *)malloc(sizeof(COMP_PRECISION)*lmsize);
    a2=(COMP_PRECISION *)malloc(sizeof(COMP_PRECISION)*lmsize);
    b2=(COMP_PRECISION *)malloc(sizeof(COMP_PRECISION)*lmsize);
    if(!a1 || !a2 || !b1 || !b2)MEMERROR;
    dr=0.1;
    for(z=zmin;z<=zmax+EPS_COMP_PREC;z+=dz){
      r = 0;tmp[0]=1;z1=zmin;z2=zmin;
      while((tmp[0] > rdist)&&(z1>=zmin)&&(z1<=zmax)&& 
	    (z2>=zmin)&&(z2<=zmax)){
	r += dr;
	a = r * HALF_SQRT_TWO;
	z1 = z - a;
	z2 = z + a;
	interpolate_she_model(a1,b1,&model,z1,model.lmax,FALSE);
	interpolate_she_model(a2,b2,&model,z2,model.lmax,FALSE);
	tmp[0]=correlation(a1,b1,a2,b2,-model.lmax,0,1,cmode);
      }
      if(tmp[0] <= rdist)
	fprintf(stdout,"%11g %11g\n",z,r);
      else
	fprintf(stdout,"%11g NaN\n",z,r);
    }
    break;
  }
  default:{
    fprintf(stderr,"%s: mode %i is undefined\n",argv[0],mode);
    exit(-1);
  }}
  return 0;
}


