#include <math.h>
#include <stdio.h>
#include <string.h>
#include "calc_mag.h"
#include "calc_mb.h"

double mb_table[43][29];
double mb_0_table[500];

double calc_mb (double amplitude, double frequency, Station_Struct sta_info, double lat, double lon, double depth) {
  static int first=1;
  int test=-1;
  double period=1/frequency;
  double correction=log10(amplitude/period);
  double delta=calc_dist(sta_info.lat,sta_info.lon,lat,lon);
  char mb_table_file[]="/usr/contrib/data/bdsn/mb.chart";
  char mb_0_table_file[]="/usr/contrib/data/bdsn/mb_0.chart";

  if (first) {
    first=0;
    if (read_in_mb_table(mb_table_file)==-1) {
      return 0;
    }
     if (read_in_mb_0_table(mb_0_table_file)==-1) { 
       return 0; 
     } 
  }

  if (depth!=0) {
    return q_intrp(delta,depth)+correction;
  } else {
    return mb_0_table[integer(delta)]+correction;
  }
}

int read_in_mb_table(char *file) {
  FILE *in_file=fopen(file,"r");
 char in_line[500],*word;
 int i,j;

 if (in_file==NULL) {
   fprintf(stderr,"Can not open %s\n",file);
   return -1;
 }

 for(i=0;i<43;i++) {
   fgets(in_line,500,in_file);
   word=in_line;

   if (!isdigit(in_line[0])) {
     i-=1;
     continue;
   }

   for(j=0;j<29;j++) {
     mb_table[i][j]=atof(word);
     word+=6;
   }
 }
 return 0;
}

int read_in_mb_0_table(char *file) {
  FILE *in_file=fopen(file,"r");
  char in_line[500],*word;
  double delta,corr;
  char *token;
  char line[200];
  int i;

  if (in_file==NULL) {
    fprintf(stderr,"Can not open %s\n",file);
    return -1;
  }

  for (i=0;i<200;i++) {
    mb_0_table[i]=-1;
  }

  while(!feof(in_file)) {
    if (fgets(line,200,in_file)==NULL) {break;};
    token=strtok(line," ");
    delta=atoi(token);
    token=strtok(NULL," ");
    corr=atof(token);
    
    mb_0_table[integer(delta)]=corr;
  }
  
  return 0;
}

int integer(double number) {
  int temp;
  temp=floor(number+.5);
  return number;
}

double l4pi( double fm1, double f0, double fp1, double fp2, double p){
  float one=1.,two=2.,six=6.;
  double am1,a0,ap1,ap2,f;
  
  am1=-p*(p-one)*(p-two)/six;
  a0 = (p*p-one)*(p-two)/two;
  ap1=-p*(p+one)*(p-two)/two;
  ap2=       p*(p*p-one)/six;
  f=am1*fm1+a0*f0+ap1*fp1+ap2*fp2;
  
  return(f);
}

double q_intrp(float fx, float fy) {
  double f1,f2,f3,f4,ph=0.,pd=0.,qval;
  double temp;
  int x,y;
  
  y=integer(fy/25);
  if(y<1) {
    y=1;
  }
  else
    if(y>26){
      y=26;
    }
  ph=(fy-y*25)/25;
  
  fx-=5;
  x=integer(fx/2.5);
  if(x<1){
    x=1;
  }
  else
    if(x>40){
      x=40;
    }
  pd=(fx-2.5*x)/2.5;	
  
  if(fy<=2.5) qval=l4pi(mb_table[x-1][1],mb_table[x][1],mb_table[x+1][1],
			 mb_table[x+2][1],pd);
  else {
    f1=l4pi(mb_table[x-1][y-1],mb_table[x][y-1],
            mb_table[x+1][y-1],mb_table[x+2][y-1],pd);
    f2=l4pi(mb_table[x-1][y],mb_table[x][y],
            mb_table[x+1][y],mb_table[x+2][y],pd);
    f3=l4pi(mb_table[x-1][y+1],mb_table[x][y+1],
            mb_table[x+1][y+1],mb_table[x+2][y+1],pd);
    f4=l4pi(mb_table[x-1][y+2],mb_table[x][y+2],
            mb_table[x+1][y+2],mb_table[x+2][y+2],pd);
    qval=l4pi(f1,f2,f3,f4,ph);
  }
  
  return(qval);
}