/* part of the shansyn spherical harmonics package, see COPYRIGHT for license */
/* $Id: model2scatter.c,v 1.2 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"

//
// reads two spherical harmonics models  and writes all of their
// coefficients as x y scatter data to stdout (linear regression, e.g.)
//
int main(int argc, char **argv)
{
 
  int i,lmax,lmsize,l,m,weighted=0;
  COMP_PRECISION dz,zmin,zmax,z,*a,*b,r,w;
  struct mod model[2];
  // boundaries
  zmin=50.0;
  zmax=2850;
  dz=50.0;

  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],"%i",&weighted);
    break;
  }
  default:{
    fprintf(stderr,"%s file1 file2 [zmin(%g) zmax(%g) dz(%g)] [weighted, 0]\n",
	    argv[0],zmin,zmax,dz);
    fprintf(stderr,"writes all non zero coefficients as x y data for both models to stdout\n");
    fprintf(stderr,"(can be used for linear regression, say). if weighted is set, will weigh by r^2\n");
    exit(-1);
    break;
  }}
  fprintf(stderr,"%s: weighted: %i\n",argv[0],weighted);
  // 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;
  for(z=zmin;z<=zmax+EPS_COMP_PREC;z+=dz){
    r=6371.0-z;
    /* don't allow extrapolation */
    interpolate_she_model(a,b,model,z,lmax,FALSE);
    interpolate_she_model((a+lmsize),(b+lmsize),(model+1),z,lmax,FALSE);
    /*  */
    w=(weighted)?(r*r/40589641.):(1.0);
    for(l=0;l<=lmax;l++)
      for(m=0;m<=l;m++){
	fprintf(stdout,"%12.5e %12.5e\n",w*a[POSLM(l,m)],w*a[lmsize+POSLM(l,m)]);
	if(m!=0)
	  fprintf(stdout,"%12.5e %12.5e\n",w*b[POSLM(l,m)],w*b[lmsize+POSLM(l,m)]);
      }
  }
  return 0;
}


