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

#define DIRSHORE "/home/SHARE/data/shore"       /* 40 bytes max */

int plot(double xp, double yp, int md);
int qpatr(int mpen, int *lcol, int *ltyp, double *wthick);
int penatr(int jpen, int icol, int ityp, double thick);
int newpen(int npen);

static int il=20, ih=50, kl=120, kh=155;
static int mg=1;

void rshore(int m, int is, int in, int kw, int ke)
{ if (m == 1) mg = 1;
  else if (m == 0) mg = 0;
  else { fprintf(stderr, "rshore: geodetic must be 1 or 0\n"); exit(1); }
  il = is; ih = in; kl = kw; kh = ke;
  if (il < -90) il = -90;
  if (ih > +90) ih = +90;
  if (kl < -180) kl = -180;
  if (kh > +180) kh = +180;
}

void pshore(double xorg, double yorg, int mpen, double scl,
            void (*conv)(float *alat, float *alon, float *xe, float *yn))
{ int nn, jj, i, ii, k, kk, m, n, ipen, lat, lon, lpen, icol0, ityp0;
  char cgj, cgi, cgk, nam[60], cs[3][6]={"coast", "river", "prefb"};
  float alat, alon, x[200], y[200];
  double thick0;
  FILE *fp;
  mpen = abs(mpen);
  lpen = qpatr(0, &icol0, &ityp0, &thick0);
  for (nn=0; nn<3; nn++) {
    ipen = mpen%10; mpen /= 10;
    if (ipen == 0) continue;
    newpen(ipen);
    for (i=il; i<ih; i++) {
      if (i < 0) { ii = -i; cgi = 'S'; }
      else { ii = i; cgi = 'N'; }
      cgj = ((mg==1)? 'Z': 't');
      for (k=kl; k<kh; k++) {
        jj = (k+180)/6 + 1;
        if (k < 0) { kk = -k; cgk = 'W'; }
        else { kk = k; cgk = 'E'; }
        sprintf(nam, "%s/%c%02d/%c%02i%c%03i.%s",
                DIRSHORE, cgj,jj, cgi,ii, cgk,kk, cs[nn]);
        if ((fp=fopen(nam,"r")) == NULL) continue;
        while (fscanf(fp, "%d", &n) != EOF) {
          for (m=0; m<n; m++) {
            fscanf(fp,"%d%d",&lat,&lon);
            alat = (i*6000+lat)/100.F; alon = (k*6000+lon)/100.F;
            conv(&alat, &alon, &x[m], &y[m]);
          }
          plot(xorg+x[0]/scl, yorg+y[0]/scl, 3);
          for (m=1; m<n; m++) plot(xorg+x[m]/scl, yorg+y[m]/scl, 2);
        }
        fclose(fp);
      }
    }
  }
  if (lpen == 0) penatr(0, icol0, ityp0, thick0); else newpen(lpen);
}

void wshore(double xorg, double yorg, int mpen, double scl,
            void (*conv)(float *alat, float *alon, float *xe, float *yn))
{ int nn, m, n, ipen, lat, lon, lpen, icol0, ityp0;
  char bf[82], nam[40], ew[2], cs[2][7]={"wcoast", "wriver"};
  float alat, alon, x[200], y[200];
  double thick0;
  FILE *fp;
  mpen = abs(mpen);
  lpen = qpatr(0, &icol0, &ityp0, &thick0);
  for (nn=0; nn<2; nn++) {
    if ((mpen==99) || (mpen==999))
      { nn++; sprintf(nam, "%s/world.cst", DIRSHORE); }
    else {
      ipen = mpen%10; mpen /= 10;
      if (ipen == 0) continue;
      sprintf(nam, "%s/%s.jpn", DIRSHORE,cs[nn]);
      newpen(ipen);
    }
    if ((fp=fopen(nam,"r")) == NULL) continue;
    while (fgets(bf,82,fp) != NULL) {
      if (bf[0] == '#') continue;
      if (sscanf(bf, "%s%d%1s", nam,&n,ew) < 2) continue;
      for (m=0; m<n; m++) {
        fscanf(fp,"%d%d",&lat,&lon);
        alat = (float)lat; alon = (float)lon;
        if ((mpen==999) && (ew[0]=='W')) alon-=21600.;
        conv(&alat, &alon, &x[m], &y[m]);
      }
      plot(xorg+x[0]/scl, yorg+y[0]/scl, 3);
      for (m=1; m<n; m++) plot(xorg+x[m]/scl, yorg+y[m]/scl, 2);
    }
    fclose(fp);
  }
  if (lpen == 0) penatr(0, icol0, ityp0, thick0); else newpen(lpen);
}

void rshore_(int *m, int *is, int *in, int *kw, int *ke)
{ if (*m == 0) rshore(0, *is, *in, *kw, *ke);
  else if (*m == 1) rshore(1, *is, *in, *kw, *ke);
  else rshore(0, *m, *is, *in, *kw);
}

void pshore_(float *xorg, float *yorg, int *mpen, float *scl,
            void (*conv)(float *alat, float *alon, float *xe, float *yn))
{ pshore((double)*xorg, (double)*yorg, *mpen, (double)*scl, conv); }

void wshore_(float *xorg, float *yorg, int *mpen, float *scl,
            void (*conv)(float *alat, float *alon, float *xe, float *yn))
{ wshore((double)*xorg, (double)*yorg, *mpen, (double)*scl, conv); }