/* part of the shansyn spherical harmonics package, see COPYRIGHT for license */
/* $Id: interpolate_she_model.c,v 1.10 2002/02/09 20:46:08 becker Exp becker $ */
#include <stdio.h>
#include <stdlib.h>
#include <malloc.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"

void polint(COMP_PRECISION *,COMP_PRECISION *,int ,COMP_PRECISION ,
	    COMP_PRECISION *,COMP_PRECISION *);
COMP_PRECISION my_make_nan(void);
/*

  interpolate a SHE at depth z

*/
#define LINEAR_INTERPOLATION

int interpolate_she_model(COMP_PRECISION *a, 
			   COMP_PRECISION *b,
			   struct mod *model, 
			   COMP_PRECISION z, 
			   int lmax,
			   int allow_lin_model_extrapolation) /*  */
{
  int j,k,l,m,i,interpolate;
  COMP_PRECISION fac,fac2,*ca,*cb;
  
  switch(model->type){
  case DISCRETE:{
    // discrete layers, interpolate layers somehow
    
#ifdef LINEAR_INTERPOLATION
    // check for levels
    for(i=1;i<model->n;i++)
    if(model->d[i]>=model->d[i-1]){
      fprintf(stderr,"interpolate_she_model: depth levels in discrete model not ordered in descending order\n");
      exit(-1);
    }

    if((z > model->d[0]) || (z < model->d[model->n-1])){ /* outside bounds */
      if(allow_lin_model_extrapolation){
	interpolate = TRUE;
	fprintf(stderr,"interpolate_she_model: WARNING: extrapolating at z: %g\n",z);
      }else
	interpolate = FALSE;
    }else{
      interpolate = TRUE;
    }

    if(interpolate){		/* regular operation */
      i=model->n-1;
      while((i>0)&&(model->d[i]<z))
	i--;
      if(i == (model->n - 1))
	i=model->n-2;
      j=i+1;
      fac=(z - model->d[j])/(model->d[i]-model->d[j]);
      fac2=1.0-fac;
      for(l=0;l<=lmax;l++)
	for(m=0;m<=l;m++){
	  k=POSLM(l,m);
	  a[k] =
	    fac  * *(model->a[i]+k) + fac2 * *(model->a[j]+k);
	  b[k] =
	    fac  * *(model->b[i]+k) + fac2 * *(model->b[j]+k);
	}
    }else{
      /* return NaN */
      for(l=0;l<=lmax;l++)
	for(m=0;m<=l;m++){
	  k=POSLM(l,m);
	  a[k] = my_make_nan();
	  b[k] = my_make_nan();
	}
      return 1;
    }
#endif
    break;
  }
  case CHEBYSHEV:{
    if((ca=(COMP_PRECISION *)malloc(sizeof(COMP_PRECISION)*
				   model->n))==NULL)MEMERROR;
    if((cb=(COMP_PRECISION *)malloc(sizeof(COMP_PRECISION)*
				   model->n))==NULL)MEMERROR;
    for(l=0;l<=lmax;l++)
      for(m=0;m<=l;m++){
	k=POSLM(l,m);
	for(i=0;i<model->n;i++){
	  ca[i]= *(model->a[i]+k);
	  cb[i]= *(model->b[i]+k);
	}
	a[k]=CHEBEV_FUNC(model->dmin,model->dmax,ca,model->n,z);
	b[k]=CHEBEV_FUNC(model->dmin,model->dmax,cb,model->n,z);
      }
    free(ca);free(cb);
    break;
  }
  case SPLINES:{
    if((ca=(COMP_PRECISION *)malloc(sizeof(COMP_PRECISION)*
				    model->n))==NULL)MEMERROR;
    if((cb=(COMP_PRECISION *)malloc(sizeof(COMP_PRECISION)*
				    model->n))==NULL)MEMERROR;
    for(l=0;l<=lmax;l++)
      for(m=0;m<=l;m++){
	k=POSLM(l,m);
	for(i=0;i<model->n;i++){
	  ca[i]= *(model->a[i]+k);
	  cb[i]= *(model->b[i]+k);
	}
	a[k]=spline_base(model->dmin,model->dmax,ca,model->n,z);
	b[k]=spline_base(model->dmin,model->dmax,cb,model->n,z);
      }
    free(ca);free(cb);
    break;

  }
  default:{
    fprintf(stderr,"interpolate_she_model: can not deal with model type %i\n",
	    model->type);
    exit(-1);
  }}
  return 0;
} 

/*

  obtain nmodel mean expansions by averaging 
  nmodel spherical harmonic models 
  from z1 to z2 in steps steps
  

  as opposed to interpolate_she_model, arrays are allocated here



*/

void mean_expansions(int lmax, COMP_PRECISION **am, 
		     COMP_PRECISION **bm,
		     int nmodel, struct mod *model, 
		     COMP_PRECISION z1, COMP_PRECISION z2,
		     int steps,int use_r2_weights)
{
  COMP_PRECISION *a,*b,dz,z,r,w,ws=0.;
  int i,l,m,lmsize,j=0;
  if(z2 < z1){
    fprintf(stderr,"mean_expansions: z2 should be bigger than z1, %g %g\n",
	    z1,z2);
    exit(-1);
  }
  lmsize=(int)((((COMP_PRECISION)lmax)+1.0)*
	       (((COMP_PRECISION)lmax)+2)/2.0);
  // temporary, one set of coefficients
  a=(COMP_PRECISION *)malloc(sizeof(COMP_PRECISION)*lmsize);
  b=(COMP_PRECISION *)malloc(sizeof(COMP_PRECISION)*lmsize);
  fprintf(stderr,"ok0\n");
  if(!a || !b)
    MEMERROR;
  //
  // output values for nmodel models
  //
  *am = (COMP_PRECISION *)malloc(sizeof(COMP_PRECISION)*
				 lmsize*nmodel);
  *bm = (COMP_PRECISION *)malloc(sizeof(COMP_PRECISION)*
				 lmsize*nmodel);
  if(!*am || !*bm)
    MEMERROR;
  // init mean expansions with zeroes
  for(i=0;i < nmodel;i++)
    for(l=0;l <= lmax;l++)
      for(m=0;m <= l;m++){
	*(*am+i*lmsize+POSLM(l,m))= 
	  *(*bm+i*lmsize+POSLM(l,m))=0.0;
      }
  //
  // delta z
  //
  dz=(z2-z1)/(COMP_PRECISION)(steps-1);
  if(dz>0)
    for(i=0;i<nmodel;i++){// loop over models
      ws=0.0;// sum of weights    
      // average from z1 to z2
      for(j=0,z=z1;z <= z2+1e-5;z += dz,j++){
	interpolate_she_model(a,b,(model+i),z,lmax,TRUE); /* allow extrapolation */
	if(use_r2_weights == 0){// simple mean
	  for(l=0;l<=lmax;l++)// add to mean
	    for(m=0;m<=l;m++){
	      *(*am+POSLM(l,m)+i*lmsize) += a[POSLM(l,m)];
	      *(*bm+POSLM(l,m)+i*lmsize) += b[POSLM(l,m)];
	    }
	  ws += 1.0;
	}else{// weight by radius^2, assuming constant Delta z
	  r=(REARTH-z)/REARTH;// radius
	  if((r>1) || (r<0)){
	    fprintf(stderr,"mean_models: error, r (%g) out of range\n",r);
	    exit(-1);
	  }
	  w=r*r;// weight
	  for(l=0;l<=lmax;l++)// add to mean
	    for(m=0;m<=l;m++){
	      *(*am+POSLM(l,m)+i*lmsize) += a[POSLM(l,m)] * w;
	      *(*bm+POSLM(l,m)+i*lmsize) += b[POSLM(l,m)] * w;
	    }
	  ws += w;
	}
      }
    }
  else
    ws=0.0;
  if(ws == 0.0){
    fprintf(stderr,"mean_models: error, no averaging, j: %i ws: %g\n",j,ws);
    exit(-1);
  }
  for(i=0;i<nmodel;i++)
    for(l=0;l<=lmax;l++)
      for(m=0;m<=l;m++){
	*(*am+POSLM(l,m)+i*lmsize) /= ws;
	*(*bm+POSLM(l,m)+i*lmsize) /= ws;
      }
  free(a);free(b);
}

#define GMT_make_fnan(x) (((unsigned int *) &x)[0] = 0x7fffffff)
#define GMT_make_dnan(x) (((unsigned int *) &x)[0] = 0xffffffff, ((unsigned int *) &x)[1] = 0x7fffffff)

COMP_PRECISION my_make_nan(void)
{
  static int been_here = 0;
  static COMP_PRECISION val;
  if(!been_here){
    //GMT_make_fnan(val);
    val = sqrt(-1.0);
  }
  return val;
}
