/* 

   part of the shansyn spherical harmonics package, see COPYRIGHT for license 

   $Id: read_she_model.c,v 1.9 2001/07/14 19:51:26 becker Exp becker $ 

*/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include <limits.h>
#include "function_macros.h"
#include "trig_constants.h"
#include "precision.h"
#include "spherical_harmonics_functions.h"
#include "legendre_macros.h"
#include "myio.h"


/*
  read in spherical harmonic expansion model with depth layers
  
  the first integer, the layer number n, determines the type of 
  model. 

  0<n<1000: discrete layer, will need to read depths
  n < 0   : Chebyshev expansion of order -n
  1000<=n : splines interpolation of order n-1000
  
  fname: filename
  lmax: 
         if the input value lmax is set to <=0, no change in l_max
         if > 0, the model resolution will be changed to l_max,
         adding zeroes if necessary

  nexp:  number of expansions given in a column (1,2,3..; normally 1)

*/
void read_she_model(char *fname,struct mod *model,
		    int lmax,int nexp)
{
  FILE *in;
  int j,lms,l,m,i;
  lms = 0;			/* for compiler */
  model->dmin=FLT_MAX;model->dmax=-FLT_MAX;
  in=fopen(fname,"r");
  if(!in){
    fprintf(stderr,"read_she_model: error opening spherical harmonics model file \"%s\"\n",
	    fname);
    exit(-1);
  }
  
  if(fscanf(in,"%i",&model[0].n)!=1){
    fprintf(stderr,"read_she_model: read error model file %s, could not read nl\n",fname);
    exit(-1);
  }
  if(model[0].n < 0){
    model[0].n *= -1;
    model[0].type=CHEBYSHEV;
  }else{
    if(model[0].n > 1000){
      model[0].type = SPLINES;
      model[0].n -= 1000.0;
    }else
      model[0].type = DISCRETE;
  }
  if(model[0].n < 1){
    fprintf(stderr,"read_she_model: read error model file %s, number of layers is %i\n",
	    fname,model[0].n);
    exit(-1);
  }
  model[0].lmax = 0;
  for(i=1;i<nexp;i++){
    model[i].n = model[0].n;
    model[i].type = model[0].type;
    model[i].lmax = model[0].lmax;
  }
  for(i=0;i<nexp;i++)
    if(!(model[i].d = (COMP_PRECISION *)
	 malloc(sizeof(COMP_PRECISION)*model[i].n)))
    MEMERROR;
  for(j=0;j<model[0].n;j++){
    /* 
       depth loop
    */
    if(model[0].type == DISCRETE){// for discrete model, need to read in depth 
      if(fscanf(in,DATA_FSCAN_FORMAT,&model[0].d[j])!=1){
	fprintf(stderr,"read_she_model: read error model file %s, layer %i, could not read d, nl=%i, lmax=%i\n",
		fname,j,model[0].n,model[0].lmax);exit(-1);}
      if(model[0].d[j] < 0){
	fprintf(stderr,"read_she_model: depth argument has to be positive, %g\n",
		model[0].d[j]);
	exit(-1);
      }
      if(j && (model[0].d[j]>model[0].d[j-1])){
	fprintf(stderr,
		"read_she_model: model %s\nread_she_model: discrete layers have to sorted be from bottom to top,\nread_she_model: largest z (positive) first, %i: %g %i: %g\n",
		fname,j,model[0].d[j],j-1,model[0].d[j-1]);
	exit(-1);
      }
      for(i=1;i<nexp;i++)
	model[i].d[j] = model[0].d[j];
    }else{// Chebyshev or splines model, read dmin and dmax only once
      if(j == 0){
	if(fscanf(in,TWO_DATA_FSCAN_FORMAT,&model[0].dmin,
		  &model[0].dmax)!=2){
	  fprintf(stderr,"read_she_model: read error dmin / dmax Chebyshev model\n");
	}
	for(i=1;i<nexp;i++){
	  model[i].dmin = model[0].dmin;
	  model[i].dmax = model[0].dmax;
	}
      }
    }
    if(model[0].type == DISCRETE){
      if(model[0].d[j] > model[0].dmax)
	model[0].dmax = model[0].d[j];
      if(model[0].d[j] < model[0].dmin)
	model[0].dmin = model[0].d[j];
      for(i=1;i<nexp;i++){
	model[i].dmin = model[0].dmin;
	model[i].dmax = model[0].dmax;
      }
    }
    if(fscanf(in,"%i",&model[0].lmax)!=1){
      fprintf(stderr,"read_she_model: read error file %s, layer %i could not read lmax\n",
	      fname,j);exit(-1);}
    for(i=1;i<nexp;i++){
      model[i].lmax = model[0].lmax;
    }
    for(i=0;i<nexp;i++){
      if(j==0){// first layer, allocate pointers for coefficients
	lms= (int)((((COMP_PRECISION)model[0].lmax)+1.0)*(((COMP_PRECISION)model[0].lmax)+2)/2.0);
	if(lms==0){
	  fprintf(stderr,"read_she_model: model %s, lms is %i\n",
		  fname,lms);
	  exit(-1);
	}
	model[i].a=(COMP_PRECISION **)malloc(sizeof(COMP_PRECISION *)*model[i].n);
	model[i].b=(COMP_PRECISION **)malloc(sizeof(COMP_PRECISION *)*model[i].n);
	if(!model[i].a || !model[i].b)MEMERROR;
      }
      model[i].a[j]=(COMP_PRECISION *)malloc(sizeof(COMP_PRECISION)*lms*model[i].n);
      model[i].b[j]=(COMP_PRECISION *)malloc(sizeof(COMP_PRECISION)*lms*model[i].n);
      if(!model[i].a[j] || !model[i].b[j])MEMERROR;
    }
    for(l=0;l<=model[0].lmax;l++){
      for(m=0;m<=l;m++){
	for(i=0;i<nexp;i++){
	  if((fscanf(in,TWO_DATA_FSCAN_FORMAT,(model[i].a[j]+POSLM(l, m)),(model[i].b[j]+POSLM(l, m))))!=2){
	    fprintf(stderr,"read_she_model: read error file %s, layer %i, l %i, m %i, model %i out of %i\n",
		    fname,j,l,m,i,nexp);exit(-1);
	  }
	}
      }
    }
  }
  fclose(in);
  for(i=0;i<nexp;i++){
    // modify l_max
    if(lmax > 0){
      lms= (int)((((COMP_PRECISION)lmax)+1.0)*(((COMP_PRECISION)lmax)+2)/2.0);
      for(j=0;j<model[i].n;j++){
	model[i].a[j]=(COMP_PRECISION *)
	  realloc(model[i].a[j],sizeof(COMP_PRECISION)*lms*model[i].n);
	model[i].b[j]=(COMP_PRECISION *)
	  realloc(model[i].b[j],sizeof(COMP_PRECISION)*lms*model[i].n);
	if(!model[i].a[j] || !model[i].b[j])MEMERROR;
	if(lmax>model[i].lmax)
	  for(l=model[i].lmax+1;l<=lmax;l++)
	    for(m=0;m<=l;m++)
	      *(model[i].a[j]+POSLM(l, m))= *(model[i].b[j]+POSLM(l, m))=0.0;
	if((j==0)&&(i==0))
	  fprintf(stderr,"read_she_model: changed l_max from %i to %i\n",model[i].lmax,lmax);
      }
      model[i].lmax=lmax;
    }
    if(i==0)
      switch(model[i].type){
      case DISCRETE:{
	fprintf(stderr,"read_she_model: read in discrete layer model %s with %i layers, lmax is %i, dmin/max : %g %g\n",
		fname,model[i].n,model[i].lmax,model[i].dmin,model[i].dmax);
	break;
      }
      case CHEBYSHEV:{
	fprintf(stderr,"read_she_model: read in Chebyshev model %s with order %i, lmax is %i, dmin/max : %g %g\n",
	      fname,model[i].n,model[i].lmax,model[i].dmin,model[i].dmax);
	break;
      }
      case SPLINES:{
	fprintf(stderr,"read_she_model: read in spline model %s with order %i, lmax is %i, dmin/max : %g %g\n",
		fname,model[i].n,model[i].lmax,model[i].dmin,model[i].dmax);
	break;
      }
      default:{
	fprintf(stderr,"read_she_model: ended up with wrong model code: %i\n",
		model[i].type);
	exit(-1);
      }}
  }

    
}

