#include    <stdio.h>
#include    <string.h>
#include    <stdlib.h>
#include    <dirent.h>

#define DIRgeoid "/home/SHARE/data/geoid"       /* 40 bytes max */

static char model[21] = "NGA";                  /* default model name */
static char mdl[2][21] = { "", "" }, llc[2][12] = { "", "" };
static int ker, kr, ks=1, kw=0, krw=0, km=0, kmsg=3;
static int mxn[2], myn[2], kalc[2] = { 0, 0 }, knvf[2] = { 0, 0 };
static double prms[2][4];
static float vnul[2];
static float *hg[2], *hgw, *hgt;

int sgeoid(char *mname)
{ DIR *pdir; FILE *fp; char bf[82], fnam[62];
  int ncc, n0, n1, n2, n3, n4, n5, i, k; float vnvd;

  if (strlen(mname) > 20)
    { fprintf(stderr, "sgeoid: too long model name\n"); return(-1); }
  if ((strcmp(mname,"")==0) || (strcmp(mname,"/")==0))
    { fprintf(stderr, "sgeoid: model '' not found\n"); return(-1); }
  if ((strcmp(mname,"..")==0) || (strcmp(mname,"../")==0))
    { fprintf(stderr, "sgeoid: model '..' not found\n"); return(-1); }
  sprintf(fnam, "%s/%s", DIRgeoid, mname);
  if ((pdir=opendir(fnam)) == NULL)
    { fprintf(stderr, "sgeoid: model '%s' not found\n", mname); return(-1); }
  closedir(pdir);
  if ((strcmp(mname,".")==0) || (strcmp(mname,"./")==0)) {
    kw = 1;
    if (krw == 0) { krw = 1;
      sprintf(fnam, "%s/world.hgeoid", DIRgeoid);
      if ((fp=fopen(fnam,"r")) == NULL) {
        fprintf(stderr, "sgeoid: world.hgeoid file not found\n");
        return(-1);
      }
      bf[0] = '#';
      while (bf[0] == '#') {
        if (fgets(bf,82,fp) == NULL) {
          fprintf(stderr, "sgeoid: world.hgeoid empty file\n"); return(-1);
        }
      }
      if (strncmp(&bf[12], "    ", 4) == 0) {       /*-- v2005  --*/
        bf[12] = '\0';
        if ((sscanf(&bf[8],"%d",&ncc) != 1) || (ncc != 399)) {
          fprintf(stderr, "sgeoid: world.hgeoid header1 error\n"); return(-1);
        }
      } else {                                      /*-- v2018 --*/
        bf[16] = '\0';
        if ((sscanf(&bf[12],"%d",&ncc) != 1) || (ncc != 199)) {
          fprintf(stderr, "sgeoid: world.hgeoid header1 error\n"); return(-1);
        }
      }
      if (fgets(bf,82,fp) == NULL) {
        fprintf(stderr, "sgeoid: world.hgeoid empty file\n"); return(-1);
      }
      if (sscanf(bf,"%d%d%d%d%d%d%f", &n0,&n1,&n2,&n3,&n4,&n5,&vnvd) != 7) {
        fprintf(stderr, "sgeoid: world.hgeoid header2 broken\n"); return(-1);
      }
      if ((n0 != -5400000) || (n1 != -10800000) ||
          (n2 != 10000) || (n3 != 15000) || (n4 != 1081) ||
          (n5 != 1441) || (vnvd != 9999.9F)) {
        fprintf(stderr, "sgeoid: world.hgeoid header2 error"); return(-1);
      }
      if ((hgw=malloc(1081*1441*sizeof(float))) == NULL)
        { fprintf(stderr, "sgeoid: memory allocation fail\n"); return(-1); }
      for (hgt=hgw,k=0; k<1441; k++) {
        for (i=0; i<1081; i++,hgt++) { fscanf(fp,"%f",hgt); }
      }
      if (fscanf(fp,"%*s") != EOF) {
        fprintf(stderr, "sgeoid: world.hgeoid format error\n"); return(-1);
      }
      fclose(fp);
    }
  } else {
    strcpy(model, mname); kw = 0;
  }
  return(0);
}
int hgeoid(double flat, double flon, double *h)
{ int li, lj, lk, lin, lke, mi, mk, i, k, j; double fx, fy;
  int ncc, iorg, korg, n0, n1, n2, n3, n4, n5;
  FILE *fp; char cns, cew, llcd[12], fnam[81];
  char bf[82]; float vnvd, h0, h1, h2, h3;

  if ((flat < -5400.) || (flat > +5400.)) {
    fprintf(stderr, "hgeoid: Lat.=%.5lf out of reasonable range\n", flat);
    return(2);
  }
  if ((flon < -10800.) || (flon > +21600.)) {
    fprintf(stderr, "hgeoid: Lon.=%.5lf out of reasonable range\n", flon);
    return(2);
  }
  if (flon > +10800.) flon -= 21600.;   /* allow +180 to +360 deg.Long. */
  if (kw == 1) {
    fx = (flat+5400.)/10.; fy = (flon+10800.)/15.;
    mi = (int)fx; fx -= (double)mi;
    mk = (int)fy; fy -= (double)mk;
    j = 1081; hgt = hgw + j*mk + mi;
    if (fx == 0.) {
      if (fy == 0.) *h = (double)*hgt;
      else          *h = ((double)*hgt)*(1.-fy) + ((double)*(hgt+j))*fy;
    } else {
      if (fy == 0.) *h = ((double)*hgt)*(1.-fx) + ((double)*(hgt+1))*fx;
      else {
        *h =  (((double)*hgt)*(1.-fy) + ((double)*(hgt+j))*fy) * (1.-fx)
            + (((double)*(hgt+1))*(1.-fy) + ((double)*(hgt+j+1))*fy) * fx;
      }
    }
    return(0);
  }

  if (flat < 0.) {
    li = 90-(int)(flat/60.+90.); if (li == 0) li = 1;
    cns = 'S'; lin = -li*60;
  } else {
    li = (int)(flat/60.); if (li == 90) li = 89;
    cns = 'N'; lin = li*60;
  }
  if (flon < 0.) {
    lk = 180-(int)(flon/60.+180.); if (lk == 0) lk = 1;
    cew = 'W'; lke = -lk*60; lj = 30 - (lk-1)/6;
  } else {
    lk = (int)(flon/60.); if (lk == 180) lk = 179;
    cew = 'E'; lke = lk*60; lj = lk/6 + 31;
  }
  sprintf(llcd,"Z%02d/%c%02d%c%03d", lj, cns,li, cew,lk);
  if ((strcmp(mdl[0],model)==0) && (strcmp(llc[0],llcd)==0)) {
    ks = 0;
  } else if ((strcmp(mdl[1],model)==0) && (strcmp(llc[1],llcd)==0)) {
    ks = 1;
  } else {
    ks = 1 - ks;
    strcpy(mdl[ks], model); strcpy(llc[ks], llcd);
    while (1) {
      ker = 0; bf[0] = '#';
      sprintf(fnam, "%s/%s/%s.hgeoid", DIRgeoid,model,llcd);
      if ((fp=fopen(fnam,"r")) == NULL) {
        if (km < 3) { fprintf(stderr, "hgeoid: %s/%s.hgeoid file not found\n",
                              model,llcd); }
        ker = 1; break;
      }
      while (bf[0] == '#') {
        if (fgets(bf,82,fp) == NULL) {
          fprintf(stderr, "hgeoid: %s/%s.hgeoid empty file\n", model,llcd);
          ker = 1; break;
        }
      }
      if (ker == 1) break;
      bf[32] = '\0';
      if (strncmp(&bf[12], "    ", 4) == 0) {   //-- v2005
        bf[12] = '\0';
        if ((sscanf(&bf[8],"%d",&ncc) != 1) || (ncc != 399)) {
          fprintf(stderr, "hgeoid: %s/%s.hgeoid header1 error\n",model,llcd);
          ker = 1; break;
        }
        if (sscanf(&bf[16],"%d%d",&iorg,&korg) != 2) {
          fprintf(stderr, "hgeoid: %s/%s.hgeoid header1 broken\n",model,llcd);
          ker = 1; break;
        }
      } else {                                  //-- v2018
        if ((sscanf(&bf[12],"%d%d%d",&ncc,&iorg,&korg) != 3) || (ncc != 199)) {
          fprintf(stderr, "hgeoid: %s/%s.hgeoid header1 error\n",model,llcd);
          ker = 1; break;
        }
      }
      if (fgets(bf,82,fp) == NULL) {
        fprintf(stderr, "hgeoid: %s/%s.hgeoid empty file\n", model,llcd);
        ker = 1; break;
      }
      if (sscanf(bf,"%d%d%d%d%d%d%f", &n0,&n1,&n2,&n3,&n4,&n5,&vnvd) != 7) {
        fprintf(stderr, "hgeoid: %s/%s.hgeoid header2 broken\n", model,llcd);
        ker = 1; break;
      }
      if ( (n2*(n4-1) != 60000) || (n3*(n5-1) != 60000) ||
           ((iorg-lin)*1000 != n0) || ((korg-lke)*1000 != n1) ) {
        fprintf(stderr, "hgeoid: %s/%s.hgeoid header2 error\n", model,llcd);
        ker = 1; break;
      }
      if (kalc[ks] == 0) {
        kalc[ks] = n4*n5;
        if ((hg[ks]=malloc(n4*n5*sizeof(float))) == NULL)
          { fprintf(stderr, "hgeoid: memory alloc fail\n"); return(-1); }
      } else if (kalc[ks] < n4*n5) {
        if ((hg[ks]=realloc(hg[ks],n4*n5*sizeof(float))) == NULL)
          { fprintf(stderr, "hgeoid: memory realloc fail\n"); return(-1); }
      }
      for (hgt=hg[ks],k=0; k<n5; k++) {
        for (i=0; i<n4; i++,hgt++) { fscanf(fp,"%f",hgt); }
      }
      if (fscanf(fp,"%*s") != EOF) {
        fprintf(stderr, "hgeoid: %s/%s.hgeoid format error\n", model,llcd);
        ker = 1; break;
      }
      fclose(fp);
      prms[ks][0] = (double)iorg + (double)n0/1000.;
      prms[ks][1] = (double)korg + (double)n1/1000.;
      prms[ks][2] = (double)n2/1000.;
      prms[ks][3] = (double)n3/1000.;
      mxn[ks] = n4; myn[ks] = n5; vnul[ks] = vnvd; break;
    }
    knvf[ks] = ker;
  }
  if (knvf[ks] == 0) {
    fx = (flat-prms[ks][0])/prms[ks][2]; mi = (int)fx; fx -= (double)mi;
    fy = (flon-prms[ks][1])/prms[ks][3]; mk = (int)fy; fy -= (double)mk;
    j = mxn[ks]; hgt = hg[ks] + j*mk + mi; kr = 0;
    if ((h0=*hgt) == vnul[ks]) kr = 2;
    else if (fx == 0.) {
      if (fy == 0.) *h = (double)h0;
      else if ((h2=*(hgt+j)) == vnul[ks]) kr = 2;
      else  *h = ((double)h0)*(1.-fy) + ((double)h2)*fy;
    }
    else {
      if ((h1=*(hgt+1)) == vnul[ks]) kr = 2;
      else if (fy == 0.) *h = ((double)h0)*(1.-fx) + ((double)h1)*fx;
      else {
        if (((h2=*(hgt+j)) == vnul[ks]) ||
            ((h3=*(hgt+j+1)) == vnul[ks])) kr = 2;
        else  *h = (((double)h0)*(1.-fy) + ((double)h2)*fy) * (1.-fx)
                  + (((double)h1)*(1.-fy) + ((double)h3)*fy) * fx;
      }
    }
  } else {
    kr = 2;
  }
  if (kr != 0) {
    if (km < 3) {
      fprintf(stderr, "hgeoid: valid data not avalable");
      fprintf(stderr, " at Lat.=%.5lf' Lon.=%.5f'\n", flat, flon);
    }
    else if (km == kmsg) {
      fprintf(stderr, "hgeoid: [valid data not avalable]");
      fprintf(stderr, " condition exceeded %d times.\n", kmsg);
      kmsg = ((kmsg%3)==0 ? (kmsg/3*10): kmsg*3);
    }
    km++;
  }
  return(kr);
}

void sgeoid_(char *name, size_t lm)
{ char mname[21]; int i, j;
  for (i=0; i<lm; i++) { if ((name[i] != ' ') && (name[i] != '/')) break; }
  for (j=0; (i+j)<lm; j++) {
    if ((name[i+j] == ' ') || (name[i+j] == '/')) break;
    else mname[j] = name[i+j];
    if (j == 20)
      { fprintf(stderr, "sgeoid: too long model name\n"); exit(1); }
  }
  mname[j] ='\0'; 
  if (sgeoid(mname) < 0) exit(1);
}
float hgeoid_(float *flat, float *flon)
{ double hh;
  if (hgeoid((double)*flat,(double)*flon,&hh) != 0) exit(1);
  return((float)hh);
}
int jhgeoid_(float *flat, float *flon, float *hg)
{ double hh; int j;
  j = hgeoid((double)*flat,(double)*flon,&hh);
  if (j == 0) { *hg = (float)hh; return(0); }
  else if (j > 0) { *hg = 0.; return(2); }
  else exit(1);
}