/*************/
/* drange2.c */
/**************************************************/
/* Get range info. (min/max/mean) of a map data   */
/*   along a specified senction location.         */
/* Normally Returns message starting with "OK"    */
/*   otherwise, returns error info without "OK"   */
/**************************************************/
/*  to compile:  gcc drange2.c -o drange2 -lm     */
/**************************************************/
/* Usage:   drange2  path  parm                   */
/*            path :  filename of input data      */
/*            parm :  m:x1;y1:x2:y2:n             */
/*              (x1,y1): Start point loc          */
/*              (x2,y2): End point loc            */
/*              m: by 0:Coord.val 1:Dist.from LL  */
/*              n: No. of points to form section  */
/**************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#define  LMX   6001
#define  LMY   6001

void jcross(int mk, float pv[7], float xpa, float ypa, float xpb, float ypb,
            int *jj, float *ds1, float *ds2, float *hi1, float *hi2)
/*************************************************************************
*--<input data>                                                          *
*   mk:  (model-code)  1 (Block), 2 (hRect), 3 (vLine), 4 (Point)        *
*   pv[0-6]: [xc, xr, yc, yr, za, zr, vol]  source model locaton         *
*   (xpa, ypa) - (xpb, ypb) :                                            *
*     section location of (A)start and (B)end points                     *
*--<output data>                                                         *
*   jj=0: section line is apart from source model                        *
*   jj=1: section line intersects to source model (Block/hRect)          *
*          along line-segment (from ds1 to ds2) in plan view             *
*   jj=2: section line is close to source model (vLine/Point) at (ds1)   *
*          (within small clearance of 1m) in plan view                   *
*   ds1: distance(1) from (xpa,ypa) (km)                                 *
*   ds2: distance(2) from (xpa,ypa) (km) [0. if jj==2]                   *
*   hi1: source top altitude (-za) (m)                                   *
*-- hi2: source bottom altitude (-za-zr) (m) [0. if vLine/Point model]   *
*************************************************************************/
{ int kx;
  float dx, dy, ds, rd1, rd2, xp1, yp1, xp2, yp2;
  float sxl, sxh, syl, syh, x0, y0, dist, d0pa, d0pb;

  dx = xpb-xpa; dy = ypb-ypa; ds = sqrtf(dx*dx + dy*dy);
  kx = 0;
  if ((mk==1) || (mk==3)) {                     // Block/hRect model
    sxl = pv[0] - pv[1]/2.; sxh = pv[0] + pv[1]/2.;
    syl = pv[2] - pv[3]/2.; syh = pv[2] + pv[3]/2.;
    if (fabsf(dx) > fabsf(dy)) {                    //-- section nearer to N-S
      rd1 = (sxl-xpa)/dx; rd2 = (sxh-xpa)/dx;
      if ((rd1>1.) && (rd2>1.)) {
      } else if ((rd1<0.) && (rd2<0.)) {
      } else {
        yp1 = ypa + dy*rd1; yp2 = ypa + dy*rd2;
        if ((yp1>syh) && (yp2>syh)) {
        } else if ((yp1<syl) && (yp2<syl)) {
        } else {
          if (yp1 > syh)
            { rd1 = (syh-ypa)/dy; xp1 = xpa + dx*rd1; yp1 = ypa + dy*rd1; }
          if (yp2 > syh)
            { rd2 = (syh-ypa)/dy; xp2 = xpa + dx*rd2; yp2 = ypa + dy*rd2; }
          if (yp1 < syl)
            { rd1 = (syl-ypa)/dy; xp1 = xpa + dx*rd1; yp1 = ypa + dy*rd1; }
          if (yp2 < syl)
            { rd2 = (syl-ypa)/dy; xp2 = xpa + dx*rd2; yp2 = ypa + dy*rd2; }
        }
        kx = 1;
      }
    } else {                                        //-- section nearer to E-W
      rd1 = (syl-ypa)/dy; rd2 = (syh-ypa)/dy;
      if ((rd1>1.) && (rd2>1.)) {
      } else if ((rd1<0.) && (rd2<0.)) {
      } else {
        xp1 = xpa + dx*rd1; xp2 = xpa + dx*rd2;
        if ((xp1>sxh) && (xp2>sxh)) {
        } else if ((xp1<sxl) && (xp2<sxl)) {
        } else {
          if (xp1 > sxh)
            { rd1 = (sxh-xpa)/dx; yp1 = ypa + dy*rd1; xp1 = xpa + dx*rd1; }
          if (xp2 > sxh)
            { rd2 = (sxh-xpa)/dx; yp2 = ypa + dy*rd2; xp2 = xpa + dx*rd2; }
          if (xp1 < sxl)
            { rd1 = (sxl-xpa)/dx; yp1 = ypa + dy*rd1; xp1 = xpa + dx*rd1; }
          if (xp2 < sxl)
            { rd2 = (sxl-xpa)/dx; yp2 = ypa + dy*rd2; xp2 = xpa + dx*rd2; }
        }
      }
      kx = 1;
    }
    if (kx == 1) {
      *jj = 1; *hi1 = -pv[4]*1000.; *hi2 = -(pv[4]+pv[5])*1000.;
      if (rd1 < rd2) {
        *ds1 = ds*rd1; *ds2 = ds*rd2;
      } else {
        *ds1 = ds*rd2; *ds2 = ds*rd1;
      }
    } else {
      *jj = 0;
    }
  } else {                                        // vLine/Point model
/****************************************************************
*  (xpa,ypa), (xpb,ypb) を通る直線の一般式は                    *
*    A = (ypa - ypb), B = (xpb - xpa), C = (xpa*ypb - xpb*ypa)  *
*  として  Ax + By + C = 0 で表わされ，点(x0,y0) からの距離は   *
*    |A x0＋B y0＋C|／sqrt(A*A + B*B)  となる.                  *
****************************************************************/
    x0 = pv[0];  y0 = pv[2];
    dist = fabsf((xpa*ypb-xpb*ypa) - dy*x0 + dx*y0) / ds;
    if (dist < 0.001) {
      d0pa = sqrtf((x0-xpa)*(x0-xpa) + (y0-ypa)*(y0-ypa) - dist*dist);
      d0pb = sqrtf((x0-xpb)*(x0-xpb) + (y0-ypb)*(y0-ypb) - dist*dist);
      if (d0pa < d0pb) {
        if (d0pb <= ds) { kx = 2; *ds1 = d0pa; }
      } else {
        if (d0pa <= ds) { kx = 2; *ds1 = ds - d0pb; }
      }
    }
    if (kx == 2) {
      *jj = 2; *ds2 = 0.; *hi1 = -pv[4]*1000.; *hi2 = 0.;
    } else {
      *jj = 0;
    }
  }
  return;
}

int main(int argc, char *argv[])
{ int mkp[20], jj[20];
  float ds1[20], ds2[20], hi1[20], hi2[20], pam[20], pth[20], pph[20];
  float pv[20][7];
  int m, ix1, iy1, ix2, iy2, n, i, j, k, nv;
  int mdl=0, nmd=0, mk, ks;
  int nc, ncc, iorg, korg, ispa, ispb;
  int ixs, iys, mszx, mszy, nmx, nmy;
  int nc2, ncc2, iorg2, korg2, ispa2, ispb2;
  int ixs2, iys2, mszx2, mszy2, nmx2, nmy2;
  float am, the, phi, dip, dec, v0, v1, v2, v3, v4, v5, v6;
  float vnul, alt, hnul, xs, ys;
  float smx, smy, xp, yp, xinc, yinc, fx, fy, dx, dy;
  float xpa, ypa, xpb, ypb;
  float vv, vv1, vv2, vsum, vmin, vmax, hsum, hmin, hmax, deep;
  float *vm, *vh, *vs, *pm, *ph, *ps;
  char buf[81], area[9], area2[9], bnc1[5], bnc2[5];
  FILE *fin;

  strcpy(buf, "/pub/AMXS/data/");
  strcat(buf, argv[1]);
  sscanf(argv[2], "%d:%d:%d:%d:%d:%d", &m,&ix1,&iy1,&ix2,&iy2,&n);
  xpa = (float)ix1/1000.; ypa = (float)iy1/1000.;
  xpb = (float)ix2/1000.; ypb = (float)iy2/1000.;

  if ((fin=fopen(buf,"r")) == NULL)
    { printf("drange2: input-file open fail\n"); return(1); }

  while (1) {                           // read model parameters 
    if (fgets(buf,81,fin) == NULL)
      { fclose(fin); printf("drange2: no data found\n"); return(1); }
    if (buf[0] != '#') break;
    if (mdl == 0) {
      if (strncmp(buf, "### CALMAS model ###", 20) == 0) mdl = 1;
    } else if (mdl == 1) {
      if ((strncmp(buf, "##Block  1  ", 12) == 0) ||
          (strncmp(buf, "##hRect  2  ", 12) == 0) ||
          (strncmp(buf, "##vLine  3  ", 12) == 0) ||
          (strncmp(buf, "##Point  4  ", 12) == 0)) {
        if (sscanf(&buf[8], "%d%f%f%f", &mk,&am,&the,&phi) == 4) {
          mdl = 2;
        } else { mdl = -1; }
      } else if (strncmp(buf, "##AmbMF  0  ", 12) == 0) {
        if (sscanf(&buf[8], "%d%f%f%f", &mk,&am,&the,&phi) == 4) {
          dip = the; dec = phi; mdl = 3;        // model parm normal end
        } else { mdl = -1; }
      } else { mdl = -1; }
    } else if (mdl == 2) {
      if (strncmp(buf, "#   ", 4) == 0) {
        if (sscanf(&buf[2], "%f%f%f%f%f%f%f",
            &v0,&v1,&v2,&v3,&v4,&v5,&v6) == 7) {
          if (nmd < 20) {
            mkp[nmd] = mk; pam[nmd] = am; pth[nmd] = the; pph[nmd] = phi;
            pv[nmd][0] = v0; pv[nmd][1] = v1; pv[nmd][2] = v2;
            pv[nmd][3] = v3; pv[nmd][4] = v4; pv[nmd][5] = v5;
            pv[nmd][6] = v6; nmd++; mdl = 1;
          } else { mdl = 4; }                   // too many models truncated
        } else { mdl = -1; }
      }
    }
  }
  if ((mdl != 0) && (mdl != 3) && (mdl != 4)) {
    mdl = -1;                                   // invalid CALMAS inf neglected
  }                                               

  strncpy(area, buf, 8); area[8] = '\0';
  strncpy(bnc1, &buf[8], 4); bnc1[4] = '\0';
  strncpy(bnc2, &buf[12], 4); bnc2[4] = '\0';
  sscanf(&buf[16], "%d %d %d %d", &iorg,&korg,&ispa,&ispb);
  if (strcmp(bnc2,"    ") == 0) {
    ncc = 0; sscanf(bnc1, "%d", &ncc);
    if ((ncc<0) || (ncc>=400)) ncc = -1;
    else if (ncc<200) ncc += 800;
    else  ncc -= 200;
  } else { ncc = 0; sscanf(bnc2, "%d", &ncc); }
  nc = ((ncc<800)? ncc: (ncc-800));
  if ((nc>=1) && (nc<=60)) { iorg = 0; korg = (nc-30)*360 - 180; }
  else if (nc == 61) { iorg = +5400; korg = 0; }
  else if (nc == 62) { iorg = -5400; korg = 0; }
  fgets(buf,81,fin);
  sscanf(buf, "%d %d %d %d %d %d %f %f",
         &ixs,&iys,&mszx,&mszy,&nmx,&nmy,&vnul,&alt);
  if ((nmx>LMX) || (nmy>LMY))
    { fclose(fin); printf("drange2: too large array size\n"); return(1); }
  if ((nmx<2) || (nmy<2))
    { fclose(fin); printf("drange2: illegal array size\n"); return(1); }
  if ((vm=(float *)malloc(sizeof(float)*nmx*nmy)) == NULL)
    { fclose(fin); printf("drange2: malloc (indata) fail\n"); return(1); }
  if ((vh=(float *)malloc(sizeof(float)*nmx*nmy)) == NULL)
    { fclose(fin); printf("drange2: malloc (indata2) fail\n"); return(1); }
  for (pm=vm,j=0; j<nmy; j++) {
    for (i=0; i<nmx; pm++,i++) {
      fscanf(fin, "%f", &vv); *pm = vv;
    }
  }
  if (fgets(buf,81,fin) == NULL)
    { fclose(fin); printf("drange2: abnormal end of data\n"); return(1); }
  if (alt == 0.) {
    while (1) {
      if (fgets(buf,81,fin) == NULL)
        { fclose(fin); printf("drange2: 2nd set data lost\n"); return(1); }
      if (buf[0] != '#') break;
    }
    strncpy(area2, buf, 8); area2[8] = '\0';
    strncpy(bnc1, &buf[8], 4); bnc1[4] = '\0';
    strncpy(bnc2, &buf[12], 4); bnc2[4] = '\0';
    sscanf(&buf[16], "%d %d %d %d", &iorg2,&korg2,&ispa2,&ispb2);
    if (strcmp(bnc2,"    ") == 0) {
      ncc2 = 0; sscanf(bnc1, "%d", &ncc2);
      if ((ncc2<0) || (ncc2>=400)) ncc2 = -1;
      else if (ncc2<200) ncc2 += 800;
      else  ncc2 -= 200;
    } else { ncc2 = 0; sscanf(bnc2, "%d", &ncc2); }
    nc2 = ((ncc2<800)? ncc2: (ncc2-800));
    if ((nc2>=1) && (nc2<=60)) { iorg2 = 0; korg2 = (nc2-30)*360 - 180; }
    else if (nc2 == 61) { iorg2 = +5400; korg2 = 0; }
    else if (nc2 == 62) { iorg2 = -5400; korg2 = 0; }
    if ((strcmp(area,area2) != 0) || (ncc != ncc2) || (iorg != iorg2)
         || (korg != korg2) || (ispa != ispa2) || (ispb != ispb2))
      { printf("drange2: Header1 inconsistent\n"); return(1); }
    fgets(buf,81,fin);
    sscanf(buf, "%d %d %d %d %d %d %f",
           &ixs2,&iys2,&mszx2,&mszy2,&nmx2,&nmy2,&hnul);
    if ((ixs != ixs2) || (iys != iys2) || (mszx != mszx2)
         || (mszy != mszy2) || (nmx != nmx2) || (nmy != nmy2))
      { fclose(fin); printf("drange2: Header2 inconsistent\n"); return(1); }
    for (ph=vh,j=0; j<nmy; j++) {
      for (i=0; i<nmx; ph++,i++) {
        fscanf(fin, "%f", &vv); *ph = vv;
      }
    }
  } else {
    hnul = 9999.9;
    for (pm=vm,ph=vh,j=0; j<nmy; j++) {
      for (i=0; i<nmx; pm++,ph++,i++) {
        if (*pm == vnul) *ph = hnul; else *ph = ((alt<0.)? hnul: alt);
      }
    }
  }
  fclose(fin); xs = (float)ixs/1000.; ys = (float)iys/1000.;

//-- Generate section data
  if ((vs=(float *)malloc(sizeof(float)*n*2)) == NULL)
    { fclose(fin); printf("drange2: malloc (secdata) fail\n"); return(1); }
  if (m == 0) {
    ix1 -= ixs; iy1 -=iys; ix2 -= ixs; iy2 -= iys;
  }
  smx = (float)mszx; smy = (float)mszy;
  xp = (float)ix1; xinc = (float)(ix2-ix1)/(n-1);
  yp = (float)iy1; yinc = (float)(iy2-iy1)/(n-1);
  for (ps=vs,k=0; k<n; ps+=2,k++) {
    fx = xp/smx;  i = (int)fx; dx = fx -(float)i;
    if ((i==(nmx-1)) && (dx==0.)) { i--; dx = 1.; }
    fy = yp/smy;  j = (int)fy; dy = fy -(float)j;
    if ((j==(nmy-1)) && (dy==0.)) { j--; dy = 1.; }

    if ((i<0) || (dx<0.) || (j<0) || (dy<0.)) { *ps = vnul; *(ps+1) = hnul; }
    else if ((i>=(nmx-1)) || (j>=(nmy-1))) { *ps = vnul; *(ps+1) = hnul; }
    else if (*(vm+j*nmx+i) == vnul) { *ps = vnul; *(ps+1) = hnul; }
    else if (dy == 0.) {
      if (dx == 0.) { *ps = *(vm+j*nmx+i); *(ps+1) = *(vh+j*nmx+i); }
      else if (*(vm+j*nmx+i+1) == vnul) { *ps = vnul; *(ps+1) = hnul; }
      else {
        *ps = *(vm+j*nmx+i)*(1.-dx) + *(vm+j*nmx+i+1)*dx;
        if ((*(vh+j*nmx+i)==hnul) || (*(vh+j*nmx+i+1)==vnul)) {
          *(ps+1) = hnul;
        } else {
          *(ps+1) = *(vh+j*nmx+i)*(1.-dx) + *(vh+j*nmx+i+1)*dx;
        }
      }
    } else if (*(vm+j*nmx+i+1) == vnul) {
      *ps = vnul; *(ps+1) = hnul;
    } else {
      if (dx == 0.) {
        *ps = *(vm+j*nmx+i)*(1.-dy) + *(vm+(j+1)*nmx+i)*dy;
        if ((*(vh+j*nmx+i)==hnul) || (*(vh+(j+1)*nmx+i)==vnul)) {
          *(ps+1) = hnul;
        } else {
          *(ps+1) = *(vh+j*nmx+i)*(1.-dy) + *(vh+(j+1)*nmx+i)*dy;
        }
      } else if ((*(vm+(j+1)*nmx+i)==vnul) || (*(vm+(j+1)*nmx+i+1)==vnul)) {
        *ps = vnul; *(ps+1) = hnul;
      } else {
        vv1 = *(vm+j*nmx+i)*(1.-dx) + *(vm+j*nmx+i+1)*dx;
        vv2 = *(vm+(j+1)*nmx+i)*(1.-dx) + *(vm+(j+1)*nmx+i+1)*dx;
        *ps = vv1*(1-dy) + vv2*dy;
        if ((*(vh+j*nmx+i)==hnul) || (*(vh+j*nmx+i+1)==vnul) ||
            (*(vh+(j+1)*nmx+i)==hnul) || (*(vh+(j+1)*nmx+i+1)==vnul)) {
          *(ps+1) = hnul;
        } else {
          vv1 = *(vh+j*nmx+i)*(1.-dx) + *(vh+j*nmx+i+1)*dx;
          vv2 = *(vh+(j+1)*nmx+i)*(1.-dx) + *(vh+(j+1)*nmx+i+1)*dx;
          *(ps+1) = vv1*(1-dy) + vv2*dy;
        }
      }
    }
    xp += xinc; yp +=yinc;
  }

  for (nv=0,ps=vs,k=0; k<n; k++,ps+=2) {
    if (*ps == vnul) continue;
    vsum = *ps; vmin = *ps; vmax = *ps; nv = 1; break;
  }
  if (nv == 0) { printf("drange2: no valid section values\n"); return(1); }
  for (k++,ps+=2; k<n; k++,ps+=2) {
    if (*ps == vnul) continue;
    nv++; vsum += *ps;
    if (*ps < vmin) vmin = *ps;
    if (*ps > vmax) vmax = *ps;
  }
  vsum = vsum / (float)nv;
  for (nv=0,ps=vs+1,k=0; k<n; k++,ps+=2) {
    if (*ps == hnul) continue;
    hsum = *ps; hmin = *ps; hmax = *ps; nv = 1; break;
  }
  if (nv == 1) {
    for (k++,ps+=2; k<n; k++,ps+=2) {
      if (*ps == hnul) continue;
      nv++; hsum += *ps;
      if (*ps < hmin) hmin = *ps;
      if (*ps > hmax) hmax = *ps;
    }
    hsum = hsum / (float)nv;
  }

  for (ks=0,i=0; i<nmd; i++) {
    jcross(mkp[i], pv[i], xpa, ypa, xpb, ypb,
           &jj[i], &ds1[i], &ds2[i], &hi1[i], &hi2[i]);
    if (jj[i] != 0) {
      if ((mkp[i] == 1) || (mkp[i] == 3)) {
        deep = hi2[i];
      } else {
        deep = hi1[i];
      }
      ks = 1; break;
    }
  }
  if (ks == 1) {
    for (i=0; i<nmd; i++) {
      jcross(mkp[i], pv[i], xpa, ypa, xpb, ypb,
             &jj[i], &ds1[i], &ds2[i], &hi1[i], &hi2[i]);
      if (jj[i] != 0) {
        if ((mkp[i] == 1) || (mkp[i] == 3)) {
          if (hi2[i] < deep) deep = hi2[i];
        } else {
          if (hi1[i] < deep) deep = hi1[i];
        }
      }
      ks = 2;
    }
  }

  printf("OK\n#--- Ncoord =%4d\n", ncc);
  printf("%8s    %4d %7d %7d %7d %7d\n", area, ncc, iorg,korg,ispa,ispb);
  printf("            xs        ys     msz-x   msz-y    mx    my    vnul\n");
  printf("       %9.2f %9.2f %7d %7d %5d %5d %7.1f\n",
         xs, ys, mszx, mszy, nmx, nmy, vnul);
  printf("\nData Range : [%7.1f ,%7.1f ]   Mean :%7.1f\n", vmin, vmax, vsum);
  if (nv == 0) {
    printf("\nAlt. Range :   (no valid value is assigned)\n");
  } else {
    printf("\nAlt. Range : [%7.1f ,%7.1f ]   Mean :%7.1f\n", hmin, hmax, hsum);
  }
  if (ks == 2) {
    printf("    Deepest source bottom (along profile) : %.1f\n", deep);
  } else {
    printf("    There is no source section along the profile.\n");
  }

  free(vm); free(vh); free(vs);
  return(0);
}