/* part of the shansyn spherical harmonics package, see COPYRIGHT for license */
/* $Id: cradialcorr.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 the radial correlation function
//
// $Id: cradialcorr.c,v 1.5 2005/04/05 19:06:14 becker Exp becker $
//

int main(int argc, char **argv)
{
 
  int mode,i,j,lmsize;
  COMP_PRECISION *a1,*a2,*b1,*b2,zmin,zmax,dz,z1,z2,tmp[3];
  struct mod model;
  zmin=0;
  zmax=2850;
  dz=50.0;
  switch(argc){
  case 2:{
    mode=0;
    break;
  }
  case 3:{
    sscanf(argv[2],"%i",&mode);
    break;
  }
  case 6:{
    sscanf(argv[2],"%i",&mode);
    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, 0] [zmin(%g) zmax(%g) dz(%g)]\n",
	    argv[0],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");
    exit(-1);
    break;
  }}
  read_she_model(argv[1],&model,-1,1);
  switch(mode){
  case 0:{
    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);
	  tmp[1]=correlation(model.a[i],model.b[i],model.a[j],model.b[j],-((model.lmax>=20)?(20):(model.lmax)),0,1);
	  tmp[2]=correlation(model.a[i],model.b[i],model.a[j],model.b[j],-model.lmax,0,1);
	  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:{
    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);
	  tmp[1]=correlation(a1,b1,a2,b2,-((model.lmax>=20)?(20):(model.lmax)),0,1);
	  tmp[2]=correlation(a1,b1,a2,b2,-model.lmax,0,1);
	  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;
  }
  default:{
    fprintf(stderr,"%s: mode %i is undefined\n",argv[0],mode);
    exit(-1);
  }}
  return 0;
}


