/* 
   modified from Toh's code rec_3D_4.c 
   read A3 model and output ascii model in grid (grid_A3d.out)

   new format A3d base file, allowing different level of knots
   
   gung 02/2004  

   * discontinuity to be added 

*/

#include <stdio.h>
#include <stdlib.h>
#include <sys/file.h>
#include <math.h>
#include <malloc.h>
#include <fcntl.h>

#define PARAM 5
#define NVKNOT 30
#define NSTRING 100
#define NHKNOT 1500
#define MXCOEF NVKNOT*NHKNOT*PARAM

extern float spl(/*ord,nknots,knot,xi*/);
extern float *farray1();
float spbsp(float hdel, float ahdel);
int loadknots(float *knot);
float getdel(float a0, float o0, float a, float o);

main(ac, av)
  int ac;
  char **av;

{ 
  FILE *knotfp,*fcoefp2,*recfp;
  char knotfname[100],string[NSTRING];
  float oknot[NHKNOT],aknot[NHKNOT],*kntrad;
  int nknotA2[PARAM],nknotA1[PARAM];
  char block[PARAM],par;
  float rr[100];
  int i,ii,j,jd,kk,k,nn,J,ind,jump,index[NHKNOT];
  int io,ia,level[NHKNOT];
  int nhknot,nvknot,ndata,nhknot1,ndisc,nanpar;
  float dv,del,tmp;
  float bh[NHKNOT],bv,b;
  float x[MXCOEF],r,a,o,rec,step;
  float adel[6]={0, 63.4, 31.7, 15.8, 7.9, 3.95};
  char outf[80];

  if (ac != 2) stop("Usage: %s [spacing of grid]\n", av[0]);
  /*  adel[0]=adel[atoi(av[1])];*/

  step=atof(av[1]);
  
  if ((knotfp=fopen("lu_hknots","r"))==NULL) /* base file of model */
    stop("grid_mA3d_anis: Cannot open knotfp\n");
  
  fgets(string,NSTRING,knotfp);
  sscanf(string,"%d",&nhknot);

  i=0;
  while (fgets(string,NSTRING,knotfp)!=NULL) {
    sscanf(string,"%f%f%d",&oknot[i],&aknot[i],&level[i]);
    i++;
  }
  
  if (i!=nhknot) stop("grid_mA3d_anis: inconsistent in nhknot - 1");
  
  printf("Number of spherical splines %d\n",nhknot);
  fclose(knotfp);

  if ((fcoefp2=fopen("lu_coef","r"))==0)    /* new model */
    stop("grid_mA3d_anis: Cannot open new model file");
  
  fscanf(fcoefp2,"%d",&nanpar);
  if (nanpar>PARAM) stop("grid_mA3d_anis: Too many parametes");
  for (k=0;k<nanpar;k++) {
     fscanf(fcoefp2,"%d",&nknotA2[k]);
     fscanf(fcoefp2,"%d",&nknotA1[k]); 
     fscanf(fcoefp2,"%s",&block[k]); 
  }
  fscanf(fcoefp2,"%d",&ndisc);
  for (k=nanpar;k<ndisc+nanpar;k++) {
     fscanf(fcoefp2,"%d",&nknotA2[k]);
     fscanf(fcoefp2,"%s",&block[k]);
     printf("knots %d for %c\n",nknotA2[k],block[k]);
  }

  printf("Number of radial splines %d\n",nknotA1[0]);

  if (nanpar!=0) {
     kntrad=farray1(0,nknotA1[0]-1);
     for (k=0;k<nknotA1[0];k++)
       fscanf(fcoefp2,"%f",&kntrad[k]);
  }

  i=0;
  for (ii=0;ii<nanpar;ii++) {
    for (nn=0;nn<nknotA1[ii];nn++) {    
      for (j=0;j<nknotA2[ii];j++) { 
        fscanf(fcoefp2,"%g",&a);
        x[i++]=a;	   
      }
    }
  }
  
  if (ndisc!=0) {
    for (ii=nanpar;ii<ndisc+nanpar;ii++) {
      for (j=0;j<nknotA2[ii];j++) {
        fscanf(fcoefp2,"%g",&a);
	x[i++]=a;
      }
    }
  }
  
  fclose(fcoefp2);
  
  if ((recfp=fopen("grid_A3d_anis.out","w"))==0)
    stop("grid_A3d_anis_na: Cannot open recover.txt");
  
  printf("Isotropic S velocity (S), X, Topography (B) or Moho (M)?\n");
  scanf("%s",string);
  sscanf(string,"%s",&par);
  
  if (par=='S') {
    jump=0;
    ind=0;
  } else if (par=='X') {
    jump=nknotA2[0]*nknotA1[0];
    ind=1;
  } else if (par=='M') {
    jump=0;
    ind=3;
    for (ii=0;ii<nanpar;ii++) jump+=nknotA2[ii]*nknotA1[ii];
    jump+=nknotA2[2];
  } else stop("grid_A3d_anis: parameter not recognized!");  
  
  printf("Parameter: %c\n", par);
  
  if (par!='M') {
    kk=0;
    do {
      printf("Depth of horizontal slice or type -1 to break\n");
      scanf("%s",string);
      sscanf(string,"%f",&r); 
      if (r<0) break;
      if (r>(6371.0-kntrad[0]) || r<(6371.0-kntrad[nknotA1[0]-1])) {
         printf("\nr %f rad[0] %f rad[last] %f\n",r,kntrad[0],kntrad[nknotA1[0]-1]);
         printf("\nDepth out of range!! \n\n ");
      } else {     
         rr[kk++]=6371.0-r;
         fprintf(recfp," %10.2f ",r);
      }
    } while(r>0);
  
    fprintf(recfp,"\n");
  }
  
  for (io=-45;io<=85;io=io+step) {
    for (ia=5;ia<=60;ia=ia+step) {
      a=ia*1.0; 
      o=io*1.0;
      J=0;
      for (j=0;j<nknotA2[ind];j++) {
        del=getdel(a,o,aknot[j],oknot[j]);
        if (del<=adel[level[j]]*2) {
	  index[J]=j;
	  bh[J++]= spbsp(del,adel[level[j]]);
	}
      }
      
      fprintf(recfp,"%7.3f %7.3f ",o,a);
      
      if (par=='S' || par=='X') {
        for (k=0;k<kk;k++) {	  	
	  rec=0.;       
	  for (jd=0;jd<nknotA1[ind];jd++) {
	    bv=spl(jd,nknotA1[ind],kntrad,rr[k]);
	    for(j=0;j<J;j++){
	      b=bh[j]*bv;
	      i=index[j]+nknotA2[ind]*jd+jump;
	      rec=rec+x[i]*b;	    
	    }
	  }
	  fprintf(recfp,"%10.5f\t",rec*100.0);  /* IN percentage */
        }
      } else if (par=='M') {
        rec=0.;
        for (j=0;j<J;j++) {
	  i=index[j]+jump;
	  rec=rec+x[i]*bh[j];
	}  
        fprintf(recfp,"%10.5f\t",-1*(rec*6346.63-24.37));  /* IN km */
      } else stop("grid_A3d_anis_na: parameter not recognized!");
        	
      fprintf(recfp,"\n");
    }
   }
   
   fclose(recfp);
}

/*****************************************************************************/

float spbsp(float hdel, float ahdel)
{
  float x;
    if(hdel < ahdel) 
      x=0.75*(hdel/ahdel)*(hdel/ahdel)*(hdel/ahdel)-1.5*(hdel/ahdel)*(hdel/ahdel)+1.;
    else if(hdel <= ahdel*2.) {
      x=0.25*(2.-hdel/ahdel)*(2.-hdel/ahdel)*(2.-hdel/ahdel);
      if(x < 0)x=0.; 
      }
    else
	x=0.;

    return x ;  
}



float getdel(float a0, float o0, float a, float o) 
{
  double d2r,r2d;
  float del;
  double q0,sq0,cq0,q,sq,cq,ff,sff,cff,arg;
  

  d2r=3.1415926535/180.;
  r2d=180./3.1415926535;

  q0=(90.-a0)*d2r;
  sq0=sin(q0);
  cq0=cos(q0);

  q=(90-a)*d2r;
  sq=sin(q);
  cq=cos(q);

  ff=(o-o0)*d2r;
  sff=sin(ff);
  cff=cos(ff);

  arg=cq*cq0+sq*sq0*cff;
  if(arg > 1.) arg= 1.;
  if(arg < -1.) arg=-1.;
  del=r2d*acos(arg);

      return del;

}