/* sample program for getting P and S wave velocities at lat, lon and dep 

          usage: sample lat. lon. dep.

                                  Programed by Masaki Nakamura (29/Sep./2008) */
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include <stdarg.h>
#include <float.h>
#include <ctype.h>
#include <string.h>
#include <malloc.h>
#include <locale.h>
#include <unistd.h>
#include <limits.h>
#include <sys/types.h>
#include <sys/stat.h>

/* some defines */
#define MIN_NRAY 10
#define WEIGHT_THRESHOLD 0.025 /* 1.0 / 8 * 0.2 */
#define MODEL_P "./model.p"
#define MODEL_S "./model.s"
#define MODEL_P_PARAM "./model.p.parm"
#define MODEL_S_PARAM "./model.s.parm"
#define STANDARD_P "./standard.p"
#define STANDARD_S "./standard.s"
#define DISCONTINUITY "./discontinuity"
#define LIMIT(a,b,z)    (((a) > (z)) ? (a) : ((b) < (z)) ? (b) : (z))

typedef struct{
    int ndiscon,nlat,nlon;
    double *lat,*lon; /* lat ... geocentoric latitude */
    double ***dep; /* [n][lat][lon] */
    char ***flag; /* 'N' ... NG   ' ' ... OK */
} Discontinuity;

typedef struct _link_info{
    short nmodel,nlat,nlon,ndep;
    /* in some cases, for vel_data, ndep >= 0
                      for vel_sub_data, ndep = -(nsub_data+1) */
} Link_info;

typedef struct _vel_data{
    double pb; /* velocity perturbation */
    double err; /* standard error */
    int nray; /* the number of rays
                 If nray < 0, abs(nray) is subscript of link_info */
    int idx; /* >-2 ... normal grid
                       (-1 ... winnowed data   else ... index to unknowns)
                 -2 ... fix value at the grid */
} Vel_data;

typedef struct _vel_model{
    short nlat,nlon,ndep;
    double *lat,*lon,*dep; /* lat ... geocentoric latitude */
    Vel_data ***vel_data; /* [lat.][lon.][dep.] */
    Vel_data ***vel_sub_data; /* [sub_data][lat.][lon.]
                                 [2*n-1] ... just above discontinuity
                                 [2*n] ... just bellow discontinuity */
    short **no_pb_max_ndep[3]; /* [3][lat.][lon.]
          [0][][] ... ndep value of vel_data for not checking pb_max 
                     (target <= ndep)
                      -1  ... check pb_max for all depth
                      >-1 ... not check it for the above mensioned depth 
          [1|2][][] ... nsub_data value of vel_sub_data for not checking pb_max
                      -1 ... no target                                        */
} Vel_model;

typedef struct{
    int nlink_info;
    short nmodel,nsub_data;
    Vel_model *vel_model;
    Link_info *link_info; /* [0] ... dummy */
} Model_space;

typedef struct{
    int ndata;
    double interval; /* unit km */
    double *vel;
} Standard_vel;

/* function proto type */
static short check_endian_3dmodel(int *nlink_info,
                                  short *nmodel,short *nsub_data);
static short check_nray(Model_space model_space,Link_info link_info[][2][2],
                        double grid_w[][2][2]);
static void cnv_big_endian2little_endian(unsigned char dat[],int ndat,
                                         int nbyte);
static double discontinuity_dep(const double coordinate[],
                                Discontinuity discontinuity,int n);
static short discontinuity_kind(double lat,double lon,short ndiscon,
                                Discontinuity discontinuity);
static void error_exit(char string[]);
static void mk_3dmodel(Model_space *model_space,Discontinuity discon,FILE *fp);
static double pb3d(Model_space model_space,short nmodel,short nlat,short nlon,
                   short ndep,Link_info *link_info);
static double pb3d_sub(Model_space model_space,short nmodel,short nlat,
                       short nlon,
                       short ndep,short layer,short nsub_flag,
                       Discontinuity discon,
                       Link_info *link_info,short *flag);
static short read_3dmodel(Model_space *model_space, FILE *fp);
static void read_discon(Discontinuity *discon,char filename[]);
static void read_standard_vel(Standard_vel *standard_vel,char *fname);
static double vel1d(double coordinate[],Model_space dummy1,
                    Standard_vel standard_vel,Discontinuity dummy2,
                    Link_info dummy3[][2][2],
                    double dummy4[][2][2],double dummy5[][2][2]);
static double vel3d(double coordinate[],Model_space model_space,
                    Standard_vel standard_vel,Discontinuity discontinuity,
                    Link_info link_info[][2][2],
                    double grid_vel[][2][2],double weight[][2][2]);
static short which_layer_only(double coordinate[],Discontinuity discontinuity);

main(int argc,char *argv[])
{
     Discontinuity discontinuity;
     Model_space model_space[2];
     Standard_vel standard_vel[2];
     Link_info link_info[2][2][2];
     double grid_w[2][2][2],grid_vel[2][2][2];
     double lat,lon,dep;
     double coordinate[3],vel[2];
     double rad = 180/M_PI;
     FILE *fp;

     lat = atof(argv[1]);
     lon = atof(argv[2]);
     dep = atof(argv[3]);

     read_discon(&discontinuity,DISCONTINUITY);
     if (!(fp = fopen(MODEL_P_PARAM,"r"))){
        fprintf(stderr,"Open File Error: File Name ... %s \n",MODEL_P_PARAM);
        exit(1);
     }
     mk_3dmodel(&model_space[0],discontinuity,fp);
     fclose(fp);
     if (!(fp = fopen(MODEL_P,"r"))){
        fprintf(stderr,"Open File Error: File Name ... %s \n",MODEL_P);
        exit(1);
     }
     read_3dmodel(&model_space[0],fp);
     fclose(fp);
     read_standard_vel(&standard_vel[0],STANDARD_P);
     if (!(fp = fopen(MODEL_S_PARAM,"r"))){
        fprintf(stderr,"Open File Error: File Name ... %s \n",MODEL_S_PARAM);
        exit(1);
     }
     mk_3dmodel(&model_space[1],discontinuity,fp);
     fclose(fp);
     if (!(fp = fopen(MODEL_S,"r"))){
        fprintf(stderr,"Open File Error: File Name ... %s \n",MODEL_S);
        exit(1);
     }
     read_3dmodel(&model_space[1],fp);
     fclose(fp);
     read_standard_vel(&standard_vel[1],STANDARD_S);

     coordinate[0] = (90.0 - lat) / rad;
     coordinate[1] = lon / rad;
     coordinate[2] = dep;
     vel[0] = vel3d(coordinate,model_space[0],standard_vel[0],
                    discontinuity,link_info,grid_vel,grid_w);
     if (check_nray(model_space[0],link_info,grid_w) == -1) vel[0] = -9.99999;
     vel[1] = vel3d(coordinate,model_space[1],standard_vel[1],
                    discontinuity,link_info,grid_vel,grid_w);
     if (check_nray(model_space[1],link_info,grid_w) == -1) vel[1] = -9.99999;

     printf("Lat. = %12.6f   Lon. = %12.6f   Dep. = %12.6f\n",lat,lon,dep);
     printf("P ... %12.6f km/s   S ... %12.6f km/s \n", vel[0],vel[1]);
}

static short check_endian_3dmodel(int *nlink_info,
                                  short *nmodel,short *nsub_data)
/* check the endian type and convert the type if necessary
   this routine is not perfect, but usually work well      */
{
   short flag;
   short tmp[2];

   flag = 0;
   if (*nmodel < 0 || *nmodel > 30 || *nsub_data < 0 || *nsub_data > 30){
      tmp[0] = *nmodel;
      tmp[1] = *nsub_data;
      cnv_big_endian2little_endian((unsigned char *)nmodel,1,sizeof(short));
      cnv_big_endian2little_endian((unsigned char *)nsub_data,1,sizeof(short));
   }
   else return(flag);

   if (*nmodel < 0 || *nmodel > 30 || *nsub_data < 0 || *nsub_data > 30){
      *nmodel = tmp[0];
      *nsub_data = tmp[1];
   }
   else{
      flag = 1;
      cnv_big_endian2little_endian((unsigned char *)nlink_info,1,sizeof(int));
   }
   return(flag);
}

static short check_nray(Model_space model_space,Link_info link_info[][2][2],
                        double grid_w[][2][2])
/* return value ... -1 ... NG   0 ... OK */
{
   int i,j,k;
   short nmodel,nlat,nlon,ndep,nsub;

   for(i=0;i<2;i++){
     for(j=0;j<2;j++){
       for(k=0;k<2;k++){
          if (grid_w[i][j][k] < WEIGHT_THRESHOLD) continue;
          nmodel = link_info[i][j][k].nmodel;
          nlat = link_info[i][j][k].nlat;
          nlon = link_info[i][j][k].nlon;
          ndep = link_info[i][j][k].ndep;
          if (ndep >= 0){
             if (model_space.vel_model[nmodel].vel_data[nlat][nlon][ndep].nray < MIN_NRAY) return(-1);
          }
          else{
             nsub = -ndep-1;
             if (model_space.vel_model[nmodel].vel_sub_data[nsub][nlat][nlon].nray < MIN_NRAY) return(-1);
          }
       }
     }
   }
   return(0);
}

static void cnv_big_endian2little_endian(unsigned char dat[],int ndat,int nbyte)
/* Convert a big-endian data array to a little-endian data array */
/* nbyte ... byte size of dat[0] */
{
    int i,j,l,m,n;
    unsigned char dummy;

    for(i=0;i<ndat;i++){
       n = i * nbyte;
       m = (i + 1) * nbyte-1;
       l = nbyte / 2;
       for(j=0;j<l;j++){
          dummy = dat[n+j];
          dat[n+j] = dat[m-j];
          dat[m-j] = dummy;
       }
    }
}

static double discontinuity_dep(const double coordinate[],
                                Discontinuity discontinuity,int n)
/* return value ... depth of n-th discontinuity (km) 
                    < -9999.0 ... error              
   coordinate ... [0] ... colatitude   [1] ... longitude   (unit: radian) */
{
    double lat,lon,dep;
    double rad = 180/M_PI;
    double pf,pf1,rf,rf1;
    double wv[4];
    int i,i1,j,j1;

    if (discontinuity.ndiscon <= n) return(-9999.9);

    lat = 90.0 - coordinate[0] * rad;
    lon = coordinate[1] * rad;
    lat = LIMIT(discontinuity.lat[0],
                discontinuity.lat[discontinuity.nlat-1],lat);
    lon = LIMIT(discontinuity.lon[0],
                discontinuity.lon[discontinuity.nlon-1],lon);

    for(i=0;i<(discontinuity.nlat-1);i++){
       i1 = i+1;
/*     if (lat >= discontinuity.lat[i] && lat < discontinuity.lat[i1]) break; */
       if (lat < discontinuity.lat[i1]) break; /* for speeding up */
    }
    for(j=0;j<(discontinuity.nlon-1);j++){
       j1 = j+1;
/*     if (lon >= discontinuity.lon[j] && lon < discontinuity.lon[j1]) break; */
       if (lon < discontinuity.lon[j1]) break; /* for speeding up */
    }

    if (i==(discontinuity.nlat-1)){
       if (lat < discontinuity.lat[0]){
          i=0;
          i1=0;
       }
       else{
          i1 = i;
       }
       pf = 1.0;
       pf1 = 0.0;
    }
    else{
       pf = (lat-discontinuity.lat[i])
          / (discontinuity.lat[i1]-discontinuity.lat[i]);
       pf1 = 1.0 - pf;
    }
    if (j==(discontinuity.nlon-1)){
       if (lon < discontinuity.lon[0]){
          j=0;
          j1=0;
       }
       else{
          j1 = j;
       }
       rf = 1.0;
       rf1 = 0.0;
    }
    else{
       rf = (lon-discontinuity.lon[j])
          / (discontinuity.lon[j1]-discontinuity.lon[j]);
       rf1 = 1.0 - rf;
    }

    wv[0] = pf1 * rf1;
    wv[1] = pf * rf1;
    wv[2] = pf1 * rf;
    wv[3] = pf * rf;
    dep = wv[0]*discontinuity.dep[n][i][j]+wv[1]*discontinuity.dep[n][i1][j]
        + wv[2]*discontinuity.dep[n][i][j1]+wv[3]*discontinuity.dep[n][i1][j1];

    return(dep);
}

static short discontinuity_kind(double lat,double lon,short ndiscon,
                                Discontinuity discontinuity)
/* return value ... -1 ... include NG discontinuity data
                     0 ... not include NG discontinuity data */
{
    int i,i1,j,j1;

/* new version but slow. why ?
    for(i1=0;i1<discontinuity.nlat;i1++){
       if (lat < discontinuity.lat[i1]) break;
    }
    if      (i1==0) i = 0;
    else if (i1==discontinuity.nlat){
       i = discontinuity.nlat-1;
       i1 = discontinuity.nlat-1;
    }
    else i = i1 - 1;

    for(j1=0;j1<discontinuity.nlon;j1++){
       if (lon < discontinuity.lon[j1]) break;
    }
    if      (j1==0) j = 0;
    else if (j1==discontinuity.nlon){
       j = discontinuity.nlon-1;
       j1 = discontinuity.nlon-1;
    }
    else j = j1 - 1;
*/

    for(i=0;i<(discontinuity.nlat-1);i++){
       i1 = i+1;
       if (lat >= discontinuity.lat[i] && lat < discontinuity.lat[i1]) break;
    }
    if (i==(discontinuity.nlat-1)){
       if (lat < discontinuity.lat[0]){
          i = 0;
          i1 = 0;
       }
       else i1 = i;
    }
    for(j=0;j<(discontinuity.nlon-1);j++){
       j1 = j+1;
       if (lon >= discontinuity.lon[j] && lon < discontinuity.lon[j1]) break;
    }
    if (j==(discontinuity.nlon-1)){
       if (lon < discontinuity.lon[0]){
          j = 0;
          j1 = 0;
       }
       else j1 = j;
    }

    if (discontinuity.flag[ndiscon][i ][j ] == 'N' ||
        discontinuity.flag[ndiscon][i ][j1] == 'N' ||
        discontinuity.flag[ndiscon][i1][j ] == 'N' ||
        discontinuity.flag[ndiscon][i1][j1] == 'N') return(-1);
    else return(0);
}

static void error_exit(char string[])
{
    fprintf(stderr,string);
    exit(1);
}

static void mk_3dmodel(Model_space *model_space,Discontinuity discon,FILE *fp)
{
     char string[82];
     int i,j,k,l;

     model_space->nsub_data = 2*(discon.ndiscon+1)-2;

     while(fgets(string,82,fp)) if (string[0] != '/') break;
     sscanf(string,"%hd",&(model_space->nmodel));

     model_space->vel_model = 
                 (struct _vel_model *)malloc(sizeof(struct _vel_model) 
                * model_space->nmodel);
     if (model_space->vel_model == NULL){
        error_exit("Memory Allocate Error: vel_model in mk_3dmodel");
     }

     for(i=0;i<model_space->nmodel;i++){
        while(fgets(string,82,fp)) if (string[0] != '/') break;
        sscanf(string,"%hd%hd%hd",&(model_space->vel_model[i].nlat),
                                  &(model_space->vel_model[i].nlon),
                                  &(model_space->vel_model[i].ndep));

        model_space->vel_model[i].lat = (double *)malloc(sizeof(double) 
                                      * model_space->vel_model[i].nlat);
        model_space->vel_model[i].lon = (double *)malloc(sizeof(double) 
                                      * model_space->vel_model[i].nlon);
        model_space->vel_model[i].dep = (double *)malloc(sizeof(double) 
                                      * model_space->vel_model[i].ndep);
        if (model_space->vel_model[i].lat == NULL ||
            model_space->vel_model[i].lon == NULL ||
            model_space->vel_model[i].dep == NULL){
           error_exit("Memory Alloc. Err.: lat, lon, dep in mk_3dmodel");
        }
        for(j=0;j<model_space->vel_model[i].nlat;j=j+10){
           if (fgets(string,82,fp)==NULL){
              fprintf(stderr,"Strange record in model parameter file \n");
              fprintf(stderr,"%s",string);
              exit(1);
           }
           for(k=0;k<10;k++){
              l=j+k;
              if (l >= model_space->vel_model[i].nlat) break;
              if (sscanf(&string[7*k],"%lf",
                         &model_space->vel_model[i].lat[l]) != 1){
                 fprintf(stderr,"Strange record in model parameter file \n");
                 fprintf(stderr,"%s",string);
                 exit(1);
              }
           }
        }
        for(j=0;j<model_space->vel_model[i].nlon;j=j+10){
           if (fgets(string,82,fp)==NULL){
              fprintf(stderr,"Strange record in model parameter file \n");
              fprintf(stderr,"%s",string);
              exit(1);
           }
           for(k=0;k<10;k++){
              l=j+k;
              if (l >= model_space->vel_model[i].nlon) break;
              if (sscanf(&string[7*k],"%lf",
                         &model_space->vel_model[i].lon[l]) != 1){
                 fprintf(stderr,"Strange record in model parameter file \n");
                 fprintf(stderr,"%s",string);
                 exit(1);
              }
           }
        }
        for(j=0;j<model_space->vel_model[i].ndep;j=j+10){
           if (fgets(string,82,fp)==NULL){
              fprintf(stderr,"Strange record in model parameter file \n");
              fprintf(stderr,"%s",string);
              exit(1);
           }
           for(k=0;k<10;k++){
              l=j+k;
              if (l >= model_space->vel_model[i].ndep) break;
              if (sscanf(&string[7*k],"%lf",
                         &model_space->vel_model[i].dep[l]) != 1){
                 fprintf(stderr,"Strange record in model parameter file \n");
                 fprintf(stderr,"%s",string);
                 exit(1);
              }
           }
        }

        model_space->vel_model[i].vel_sub_data = 
                    (struct _vel_data ***)malloc(sizeof(struct _vel_data **)
                   * model_space->nsub_data);
        if (model_space->vel_model[i].vel_sub_data == NULL){
           error_exit("Memory Alloc. Err.: vel_sub_data in mk_3dmodel");
        }
        model_space->vel_model[i].vel_data = 
                    (struct _vel_data ***)malloc(sizeof(struct _vel_data **) 
                   * model_space->vel_model[i].nlat);
        for(j=0;j<model_space->nsub_data;j++){
           model_space->vel_model[i].vel_sub_data[j] = 
                       (struct _vel_data **)malloc(sizeof(struct _vel_data *) 
                      * model_space->vel_model[i].nlat);
        }
        if (model_space->vel_model[i].vel_data == NULL ||
            model_space->vel_model[i].vel_sub_data[0] == NULL ||
            model_space->vel_model[i].vel_sub_data[1] == NULL){
           error_exit("Memory Alloc. Err.: vel_data in mk_3dmodel");
        }
        for(j=0;j<model_space->vel_model[i].nlat;j++){
           model_space->vel_model[i].vel_data[j] = 
                       (struct _vel_data **)malloc(sizeof(struct _vel_data *) 
                      * model_space->vel_model[i].nlon);
           for(k=0;k<model_space->nsub_data;k++){
              model_space->vel_model[i].vel_sub_data[k][j] = 
                   (struct _vel_data *)malloc(sizeof(struct _vel_data) 
                  * model_space->vel_model[i].nlon);
           }
           if (model_space->vel_model[i].vel_data[j] == NULL ||
               model_space->vel_model[i].vel_sub_data[0][j] == NULL ||
               model_space->vel_model[i].vel_sub_data[1][j] == NULL){
              error_exit("Memory Alloc. Err.: vel_data in mk_3dmodel");
           }
           for(k=0;k<model_space->vel_model[i].nlon;k++){
              model_space->vel_model[i].vel_data[j][k] = 
                          (struct _vel_data *)malloc(sizeof(struct _vel_data) 
                         * model_space->vel_model[i].ndep);
              if (model_space->vel_model[i].vel_data[j][k] == NULL){
                 error_exit("Memory Alloc. Err.: vel_data in mk_3dmodel");
              }
           }
        }

        for(k=0;k<3;k++){
        model_space->vel_model[i].no_pb_max_ndep[k] =
             (short **)malloc(sizeof(short *) * model_space->vel_model[i].nlat);
        if (model_space->vel_model[i].no_pb_max_ndep[k] == NULL){
           error_exit("Memory Alloc. Err.: no_pb_max_ndep in mk_3dmodel");
        }
        for(j=0;j<model_space->vel_model[i].nlat;j++){
           model_space->vel_model[i].no_pb_max_ndep[k][j] = 
                (short *)malloc(sizeof(short) * model_space->vel_model[i].nlon);
           if (model_space->vel_model[i].no_pb_max_ndep[k][j] == NULL){
              error_exit("Memory Alloc. Err.: no_pb_max_ndep in mk_3dmodel");
           }
        }
        }

     }
}

static double pb3d(Model_space model_space,short nmodel,short nlat,short nlon,
                   short ndep,Link_info *link_info)
{
     int index;
     
     index=model_space.vel_model[nmodel].vel_data[nlat][nlon][ndep].nray;
     if (index < 0){
        index = -index;
        *link_info = model_space.link_info[index];
        return(model_space.vel_model[link_info->nmodel].vel_data[link_info->nlat][link_info->nlon][link_info->ndep].pb);
     }
     else{
        link_info->nmodel = nmodel;
        link_info->nlat = nlat;
        link_info->nlon = nlon;
        link_info->ndep = ndep;
        return(model_space.vel_model[nmodel].vel_data[nlat][nlon][ndep].pb);
     }
}

static double pb3d_sub(Model_space model_space,short nmodel,short nlat,
                       short nlon,
                       short ndep,short layer,short nsub_flag,
                       Discontinuity discon,
                       Link_info *link_info,short *flag)
/* nsub_flag ... < 0 ... vel_sub_data above the boundary
                 > 0 ... vel_sub_data under the boundary
   flag ... 0 ... link_info is for vel_sub_data
            1 ... link_info is for vel_data     */
{
     int index;
     short ndiscon,nsub;

     if (nsub_flag < 0){
        ndiscon=layer-1;
        nsub=2*layer-1;
     }
     else if (nsub_flag > 0){
        ndiscon=layer;
        nsub=2*layer;
     }
     else{
        *flag = 1;
        return(pb3d(model_space,nmodel,nlat,nlon,ndep,link_info));
     }

     if (discontinuity_kind(model_space.vel_model[nmodel].lat[nlat],
                            model_space.vel_model[nmodel].lon[nlon],
                            ndiscon,discon) < 0){
        *flag = 1;
        return(pb3d(model_space,nmodel,nlat,nlon,ndep,link_info));
     }
     else{
        *flag = 0;
        index=model_space.vel_model[nmodel].vel_sub_data[nsub][nlat][nlon].nray;
        if (index < 0){
           index = -index;
           *link_info = model_space.link_info[index];
           link_info->ndep = -(nsub+1);
           return(model_space.vel_model[link_info->nmodel].vel_sub_data[nsub][link_info->nlat][link_info->nlon].pb);
        }
        else{
           link_info->nmodel = nmodel;
           link_info->nlat = nlat;
           link_info->nlon = nlon;
           link_info->ndep = -(nsub+1);
           return(model_space.vel_model[nmodel].vel_sub_data[nsub][nlat][nlon].pb);
        }
     }
}

static short read_3dmodel(Model_space *model_space, FILE *fp)
{
    int i,j,k,l;
    short endian_flag; /* 0 ... not convert endian   1 ... convert endian */

    if (fread(&(model_space->nlink_info),sizeof(int),1,fp) != 1){
       fprintf(stderr,"Data Read Error: nlink_info in read_3dmodel \n");
       return(-1);
    }
    if (fread(&(model_space->nmodel),sizeof(short),1,fp) != 1){
       fprintf(stderr,"Data Read Error: nmodel in read_3dmodel \n");
       return(-1);
    }
    if (fread(&(model_space->nsub_data),sizeof(short),1,fp) != 1){
       fprintf(stderr,"Data Read Error: nsub_data in read_3dmodel \n");
       return(-1);
    }
    endian_flag = check_endian_3dmodel(&(model_space->nlink_info),
                &(model_space->nmodel),&(model_space->nsub_data));
    for(i=0;i<model_space->nmodel;i++){
       if (fread(&(model_space->vel_model[i].nlat),sizeof(short),1,fp) != 1){
          fprintf(stderr,"Data Read Error: nlat in read_3dmodel \n");
          return(-1);
       }
       if (fread(&(model_space->vel_model[i].nlon),sizeof(short),1,fp) != 1){
          fprintf(stderr,"Data Read Error: nlon in read_3dmodel \n");
          return(-1);
       }
       if (fread(&(model_space->vel_model[i].ndep),sizeof(short),1,fp) != 1){
          fprintf(stderr,"Data Read Error: ndep in read_3dmodel \n");
          return(-1);
       }
       if (endian_flag != 0){ 
          cnv_big_endian2little_endian((unsigned char *)(&(model_space->vel_model[i].nlat)),1,sizeof(short));
          cnv_big_endian2little_endian((unsigned char *)(&(model_space->vel_model[i].nlon)), 1,sizeof(short));
          cnv_big_endian2little_endian((unsigned char *)(&(model_space->vel_model[i].ndep)),1,sizeof(short));
       }
       if (fread(model_space->vel_model[i].lat,sizeof(double),model_space->vel_model[i].nlat,fp)!=model_space->vel_model[i].nlat){
          fprintf(stderr,"Data Read Error: lat in read_3dmodel \n");
          return(-1);
       }
       if (fread(model_space->vel_model[i].lon,sizeof(double),model_space->vel_model[i].nlon,fp)!=model_space->vel_model[i].nlon){
          fprintf(stderr,"Data Read Error: lon in read_3dmodel \n");
          return(-1);
       }
       if (fread(model_space->vel_model[i].dep,sizeof(double),model_space->vel_model[i].ndep,fp)!=model_space->vel_model[i].ndep){
          fprintf(stderr,"Data Read Error: dep in read_3dmodel \n");
          return(-1);
       }
       if (endian_flag != 0){
          cnv_big_endian2little_endian(
                  (unsigned char *)(model_space->vel_model[i].lat),
                  model_space->vel_model[i].nlat,sizeof(double));
          cnv_big_endian2little_endian(
                  (unsigned char *)(model_space->vel_model[i].lon),
                  model_space->vel_model[i].nlon,sizeof(double));
          cnv_big_endian2little_endian(
                  (unsigned char *)(model_space->vel_model[i].dep),
                  model_space->vel_model[i].ndep,sizeof(double));
       }
       for(j=0;j<model_space->vel_model[i].nlat;j++){
          for(k=0;k<model_space->vel_model[i].nlon;k++){
             if (fread(model_space->vel_model[i].vel_data[j][k],sizeof(struct _vel_data),model_space->vel_model[i].ndep,fp)!=model_space->vel_model[i].ndep){
                fprintf(stderr,"Data Read Error: vel_data in read_3dmodel\n");
                return(-1);
             }
             if (endian_flag != 0){
                for(l=0;l<model_space->vel_model[i].ndep;l++){
                   cnv_big_endian2little_endian((unsigned char *)(&(model_space->vel_model[i].vel_data[j][k][l].pb)),1,sizeof(double));
                   cnv_big_endian2little_endian((unsigned char *)(&(model_space->vel_model[i].vel_data[j][k][l].err)),1,sizeof(double));
                   cnv_big_endian2little_endian((unsigned char *)(&(model_space->vel_model[i].vel_data[j][k][l].nray)),1,sizeof(int));
                   cnv_big_endian2little_endian((unsigned char *)(&(model_space->vel_model[i].vel_data[j][k][l].idx)),1,sizeof(int));
                }
             }
          }
       }
       for(j=0;j<model_space->nsub_data;j++){
          for(k=0;k<model_space->vel_model[i].nlat;k++){
            if (fread(model_space->vel_model[i].vel_sub_data[j][k],sizeof(struct _vel_data),model_space->vel_model[i].nlon,fp)!=model_space->vel_model[i].nlon){
                fprintf(stderr,"Data Read Error: vel_data in read_3dmodel\n");
                return(-1);
            }
            if (endian_flag != 0){
                for(l=0;l<model_space->vel_model[i].nlon;l++){
                   cnv_big_endian2little_endian((unsigned char *)(&(model_space->vel_model[i].vel_sub_data[j][k][l].pb)),1,sizeof(double));
                   cnv_big_endian2little_endian((unsigned char *)(&(model_space->vel_model[i].vel_sub_data[j][k][l].err)),1,sizeof(double));
                   cnv_big_endian2little_endian((unsigned char *)(&(model_space->vel_model[i].vel_sub_data[j][k][l].nray)),1,sizeof(int));
                   cnv_big_endian2little_endian((unsigned char *)(&(model_space->vel_model[i].vel_sub_data[j][k][l].idx)),1,sizeof(int));
                }
            }
          }
       }

       for(k=0;k<3;k++){
       for(j=0;j<model_space->vel_model[i].nlat;j++){
          if (fread(model_space->vel_model[i].no_pb_max_ndep[k][j],sizeof(short),model_space->vel_model[i].nlon,fp)!=model_space->vel_model[i].nlon){
             fprintf(stderr,"Data Read Error: vel_model in read_3dmodel\n");
             return(-1);
          }
          if (endian_flag != 0){
             cnv_big_endian2little_endian((unsigned char *)(model_space->vel_model[i].no_pb_max_ndep[k][j]),model_space->vel_model[i].nlon,sizeof(short));
          }
       }
       }

    }
    model_space->link_info = (struct _link_info *)malloc(sizeof(struct _link_info) * model_space->nlink_info);
    if (model_space->link_info == NULL){
       fprintf(stderr,"Memory Allocate Error: link_info in read_3dmodel\n");
       return(-1);
    }
    if (fread(model_space->link_info,sizeof(struct _link_info),model_space->nlink_info,fp)!=model_space->nlink_info){
       fprintf(stderr,"Data Read Error: link_info in read_3dmodel \n");
       return(-1);
    }
    if (endian_flag != 0){
       for(i=0;i<model_space->nlink_info;i++){
          cnv_big_endian2little_endian((unsigned char *)(&(model_space->link_info[i].nmodel)),1,sizeof(short));
          cnv_big_endian2little_endian((unsigned char *)(&(model_space->link_info[i].nlat)),1,sizeof(short));
          cnv_big_endian2little_endian((unsigned char *)(&(model_space->link_info[i].nlon)),1,sizeof(short));
          cnv_big_endian2little_endian((unsigned char *)(&(model_space->link_info[i].ndep)),1,sizeof(short));
       }  
    }
    return(0);
}

static void read_discon(Discontinuity *discon,char filename[])
{
    FILE *fp;
    char string[72];
    int i,j,k,l,m;

    if (!(fp = fopen(filename,"r"))){
       fprintf(stderr,"Open File Error: File Name ... %s \n",DISCONTINUITY);
       exit(1);
    }

    /* read header record */
    if (fgets(string,72,fp)==NULL){
       fprintf(stderr,"Strange record in %s \n",DISCONTINUITY);
       fprintf(stderr,"%s",string);
       exit(1);
    }
    if (sscanf(string,"%d%d%d",
               &(discon->ndiscon),&(discon->nlat),&(discon->nlon))!=3){
       fprintf(stderr,"Strange record in %s \n",DISCONTINUITY);
       fprintf(stderr,"%s",string);
       exit(1);
    }

    /* memory allocate for discontinuity data */
    discon->lat = (double *)malloc(sizeof(double) * discon->nlat);
    discon->lon = (double *)malloc(sizeof(double) * discon->nlon);
    if (discon->lat == NULL || discon->lon == NULL){
       error_exit("Memory Allocate Error: discon->lat,lon in read_discon\n");
    }
    discon->dep = (double ***)malloc(sizeof(double **) * discon->ndiscon);
    discon->flag = (char ***)malloc(sizeof(char **) * discon->ndiscon);
    if (discon->dep == NULL || discon->flag == NULL){
       error_exit("Memory Allocate Error: discon->dep,flag in read_discon\n");
    }
    for(i=0;i<discon->ndiscon;i++){
       discon->dep[i] = (double **)malloc(sizeof(double *) * discon->nlat);
       discon->flag[i] = (char **)malloc(sizeof(char *) * discon->nlat);
       if (discon->dep[i] == NULL || discon->flag[i] == NULL){
          error_exit("Memory Alloc. Err.: discon->dep,flag in read_discon\n");
       }
       for(j=0;j<discon->nlat;j++){
          discon->dep[i][j] = (double *)malloc(sizeof(double) * discon->nlon);
          discon->flag[i][j] = (char *)malloc(sizeof(char) * discon->nlon);
       }
    }
   
    /* read lat. data */
    for(i=0;i<discon->nlat;i=i+10){
       if (fgets(string,72,fp)==NULL){
          fprintf(stderr,"Strange record in %s \n",DISCONTINUITY);
          fprintf(stderr,"%s",string);
          exit(1);
       }
       for(j=0;j<10;j++){
          k=i+j;
          if (k >= discon->nlat) break;
          if (sscanf(&string[7*j],"%lf",&discon->lat[k]) != 1){
             fprintf(stderr,"Strange record in %s \n",DISCONTINUITY);
             fprintf(stderr,"%s",string);
             exit(1);
          }
       }
    }

    /* read lon. data */
    for(i=0;i<discon->nlon;i=i+10){
       if (fgets(string,72,fp)==NULL){
          fprintf(stderr,"Strange record in %s \n",DISCONTINUITY);
          fprintf(stderr,"%s",string);
          exit(1);
       }
       for(j=0;j<10;j++){
          k=i+j;
          if (k >= discon->nlon) break;
          if (sscanf(&string[7*j],"%lf",&discon->lon[k]) != 1){
             fprintf(stderr,"Strange record in %s \n",DISCONTINUITY);
             fprintf(stderr,"%s",string);
             exit(1);
          }
       }
    }

    /* read dep. data */
    for(i=0;i<discon->ndiscon;i++){
       for(j=0;j<discon->nlat;j++){
          for(k=0;k<discon->nlon;k=k+10){
             if (fgets(string,72,fp)==NULL){
                fprintf(stderr,"Strange record in %s \n",DISCONTINUITY);
                fprintf(stderr,"%s",string);
                exit(1);
             }
             for(l=0;l<10;l++){
                m=k+l;
                if (m >= discon->nlon) break;
                discon->flag[i][j][m] = string[7*l];
                if (sscanf(&string[7*l+1],"%lf",&discon->dep[i][j][m]) != 1){
                   fprintf(stderr,"Strange record in %s \n",DISCONTINUITY);
                   fprintf(stderr,"%s",string);
                   exit(1);
                }
             }
          }
       }
    }

    fclose(fp);
}

static void read_standard_vel(Standard_vel *standard_vel,char *fname)
{
    FILE *fp;
    int i,j,k;
    char string[82];

    if (!(fp = fopen(fname,"r"))){
       fprintf(stderr,"Open File Error: File Name ... %s \n",fname);
       exit(1);
    }

    if (fgets(string,82,fp)==NULL){
       fprintf(stderr,"Strange record in file:%s \n",fname);
       exit(1);
    }
    if (sscanf(string,"%lf%d",
               &(standard_vel->interval),&(standard_vel->ndata)) != 2){
       fprintf(stderr,"Strange record in file:%s \n",fname);
       fprintf(stderr,"%s",string);
       exit(1);
    }
    standard_vel->vel = (double *)malloc(sizeof(double)*standard_vel->ndata);
    if (standard_vel->vel == NULL){
       error_exit("malloc error: standard_vel->vel in read_standard_vel\n");
    }
    
    for(i=0;i<standard_vel->ndata;i=i+5){
       if (fgets(string,82,fp)==NULL){
          fprintf(stderr,"Strange record in file:%s \n",fname);
          fprintf(stderr,"%s",string);
          exit(1);
       }
       for(j=0;j<5;j++){
          k=i+j;
          if (k >= standard_vel->ndata) break;
          if (sscanf(&string[7*j],"%lf",&(standard_vel->vel[k])) != 1){
             fprintf(stderr,"Strange record in file:%s \n",fname);
             fprintf(stderr,"%s",string);
             exit(1);
          }
       }
    }
}

static double vel1d(double coordinate[],Model_space dummy1,
                    Standard_vel standard_vel,Discontinuity dummy2,
                    Link_info dummy3[][2][2],
                    double dummy4[][2][2],double dummy5[][2][2])
{
    double dep,vel,weight;
    int i;

    dep = coordinate[2];
    i = dep / standard_vel.interval + 0.5;
/*  this version does not work well, because velocity gradient is zero, 
    pseudo bending method says that adequate ray-path is straight line.
    if (i>(standard_vel.ndata-2)) vel = standard_vel.vel[standard_vel.ndata-1];
    else if (i<0) vel = standard_vel.vel[0];
    else{
       weight = (dep - i * standard_vel.interval) / standard_vel.interval;
       vel = standard_vel.vel[i] * (1.0 - weight)
           + standard_vel.vel[i+1] * weight;
    }
*/
/* new version */
    if (i<0){
       i = 0;
       weight = (dep - i * standard_vel.interval) / standard_vel.interval;
       vel = standard_vel.vel[i] * (1.0 - weight)
           + standard_vel.vel[i+1] * weight;
       if (vel < standard_vel.vel[0]*0.5) vel = standard_vel.vel[0]*0.5;
    }
    else{
       if (i>(standard_vel.ndata-2)) i = standard_vel.ndata-2;
       weight = (dep - i * standard_vel.interval) / standard_vel.interval;
       vel = standard_vel.vel[i] * (1.0 - weight)
           + standard_vel.vel[i+1] * weight;
    }

    return(vel);
}

static double vel3d(double coordinate[],Model_space model_space,
                    Standard_vel standard_vel,Discontinuity discontinuity,
                    Link_info link_info[][2][2],
                    double grid_vel[][2][2],double weight[][2][2])
/* link_info, weight ... [lat.][lon.][dep.] */
{
    double pb,vel,tmp;
    double lat,lon,dep;
    double rad = 180/M_PI;
    double grid[3];
    double coordinate_tmp[3];
    int nmodel,nlat,nlon,ndep,i,j,k;
    short layer[2]; /* [0] ... for target point   [1] ... for grid point */
    short dummy;

    lat = 90.0 - coordinate[0] * rad;
    lon = coordinate[1] * rad;
    dep = coordinate[2];
    for(i=0;i<2;i++) for(j=0;j<2;j++) for(k=0;k<2;k++) weight[i][j][k]=1.0;
    lat = LIMIT(model_space.vel_model[0].lat[0],model_space.vel_model[0].lat[model_space.vel_model[0].nlat-1],lat);
    lon = LIMIT(model_space.vel_model[0].lon[0],model_space.vel_model[0].lon[model_space.vel_model[0].nlon-1],lon);
    dep = LIMIT(model_space.vel_model[0].dep[0],model_space.vel_model[0].dep[model_space.vel_model[0].ndep-1],dep);

    /* which model */
    for(nmodel=(model_space.nmodel-1);nmodel>=0;nmodel--){
       if (lat < model_space.vel_model[nmodel].lat[0] ||
           lat > model_space.vel_model[nmodel].lat[model_space.vel_model[nmodel].nlat-1] ||
           lon < model_space.vel_model[nmodel].lon[0] ||
           lon > model_space.vel_model[nmodel].lon[model_space.vel_model[nmodel].nlon-1] ||
           dep < model_space.vel_model[nmodel].dep[0] ||
           dep > model_space.vel_model[nmodel].dep[model_space.vel_model[nmodel].ndep-1]) continue;
       break;
    }

    /* search lat. direction */
    for(nlat=0;nlat<(model_space.vel_model[nmodel].nlat-1);nlat++){
       if (lat <= model_space.vel_model[nmodel].lat[nlat+1]){
          tmp = (lat - model_space.vel_model[nmodel].lat[nlat])
              / (model_space.vel_model[nmodel].lat[nlat+1]
              - model_space.vel_model[nmodel].lat[nlat]);
          for(j=0;j<2;j++) for(k=0;k<2;k++){
             weight[1][j][k] *= tmp;
             weight[0][j][k] *= (1.0 - tmp);
          }
          break;
       }
    }
    /* search lon. direction */
    for(nlon=0;nlon<(model_space.vel_model[nmodel].nlon-1);nlon++){
       if (lon <= model_space.vel_model[nmodel].lon[nlon+1]){
          tmp = (lon - model_space.vel_model[nmodel].lon[nlon])
              / (model_space.vel_model[nmodel].lon[nlon+1]
              - model_space.vel_model[nmodel].lon[nlon]);
          for(i=0;i<2;i++) for(k=0;k<2;k++){
             weight[i][1][k] *= tmp;
             weight[i][0][k] *= (1.0 - tmp);
          }
          break;
       }
    }
    /* search dep. direction */
    for(ndep=0;ndep<(model_space.vel_model[nmodel].ndep-1);ndep++){
       if (dep <= model_space.vel_model[nmodel].dep[ndep+1]){
          tmp = (dep - model_space.vel_model[nmodel].dep[ndep])
              / (model_space.vel_model[nmodel].dep[ndep+1]
              - model_space.vel_model[nmodel].dep[ndep]);
          for(i=0;i<2;i++) for(j=0;j<2;j++){
             weight[i][j][1] *= tmp;
             weight[i][j][0] *= (1.0 - tmp);
          }
          break;
       }
    }

    pb = 0.0;
    layer[0] = which_layer_only(coordinate,discontinuity); 
    for(i=0;i<2;i++) for(j=0;j<2;j++) for(k=0;k<2;k++){
       grid[0] = (90.0 - model_space.vel_model[nmodel].lat[nlat+i]) / rad;
       grid[1] = model_space.vel_model[nmodel].lon[nlon+j] / rad;
       grid[2] = model_space.vel_model[nmodel].dep[ndep+k];
       layer[1] = which_layer_only(grid,discontinuity); 
       if (layer[0] != layer[1]){
          tmp = pb3d_sub(model_space,nmodel,(nlat+i),(nlon+j),(ndep+k),
                         layer[0],(layer[1]-layer[0]),discontinuity,
                         &(link_info[i][j][k]),&dummy);
          grid_vel[i][j][k] = vel1d(grid,model_space,standard_vel,discontinuity,
                                    NULL,NULL,NULL) * (1.0 + tmp * 0.01);
          tmp *= weight[i][j][k];
          pb += tmp;
       }
       else{
          tmp = pb3d(model_space,nmodel,(nlat+i),(nlon+j),(ndep+k),
                     &(link_info[i][j][k]));
          grid_vel[i][j][k] = vel1d(grid,model_space,standard_vel,discontinuity,
                                    NULL,NULL,NULL) * (1.0 + tmp * 0.01);
          tmp *= weight[i][j][k];
          pb += tmp;
       }
    }

    coordinate_tmp[2] = dep;
    vel = vel1d(coordinate_tmp,model_space,standard_vel,discontinuity,
                NULL,NULL,NULL) * (1.0 + pb*0.01);
    return(vel);
}

static short which_layer_only(double coordinate[],Discontinuity discontinuity)
/* return value ... the number of layer
   coordinate ... [0] ... colatitude   [1] ... longitude   (unit: radian)
                  [2] ... depth (unit: km)                                */
{
    double discon_dep;
    int i;
    short layer;

    layer = 0;
    for(i=0;i<(discontinuity.ndiscon);i++){
       discon_dep = discontinuity_dep(coordinate,discontinuity,i);
       if (coordinate[2] <= discon_dep) break;
       layer = layer + 1;
    }

    return(layer);
}