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

#define MAXLPCM 20

int plot(double xp, double yp, int md);
int penatr(int mpen, int icol, int ityp, double thick);
int newpen(int npen);
int qpatr(int mpen, int *lcol, int *ltyp, double *wthick);
int pcstr(double xp, double yp, double chi, char *text, double ang, int ns);

/***************************************************/
/*   static variables  and  internal subroutines   */
/***************************************************/

static double sxo, syo, sxi, syi, sxj, syj, szi, szj;
static double scs, ss[5], wd0, wd1, wd2, dci, dcj;
static unsigned char cp[8] = { 0x01,0x02,0x04,0x08,0x10,0x20,0x40,0x80 };
static unsigned char *mm1, *mm2;
static char fval[17], fmt[8];
static float vc;
static int kdir=0, icol0, ltyp0, icol1, icol2;
static int imx, jmx, kc, kpen, npen, ihsz, mbuf;

static void moveto(double x, double y)
{   double px, py;
    px = sxi*x + sxj*y; py = syi*x + syj*y;
    if (kdir == 2) plot(syo+py, sxo+px, 3);
    else plot(sxo+px, syo+py, 3);
}

static void drawto(double x, double y)
{   double px, py;
    px = sxi*x + sxj*y; py = syi*x + syj*y;
    if (kdir == 2) plot(syo+py, sxo+px, 2);
    else plot(sxo+px, syo+py, 2);
}

static void drw1cont(float *v, int ih, int il, int j, int m)
{   static int mn[16] = { 8,4,0,12, 12,0,4,8, 0,12,8,4, 4,8,12,0 };
    unsigned char *pm1, *pm2;
    int i, mp, md, khv;
    float *vp, v0, v1;
    double dv, frac, px, py, x, y, xr, yr, r, rr;
    i = ih*8+il-1; vp = v+j*imx+i; khv = (m%8)/4+1; mp = 0; pm1=mm1; pm2=mm2;
    if (npen != kpen) {
        if (kpen < 4) { npen = kpen; penatr(0, icol1, npen, wd1); }
        else { npen = kpen; penatr(0, icol2, npen-4, wd2); }
    }
    while (khv != 0) {
        if (khv == 1) { v0 = *vp; v1 = *(vp+1); dv = (v1-v0)/dci; }
        else        { v0 = *vp; v1 = *(vp+imx); dv = (v1-v0)/dcj; }
        if ((kc<=0) || ((int)dv==0)) {
            frac = (vc-v0)/(v1-v0);
            if (khv == 1) { x = (double)i + frac; y = (double)j; }
            else          { x = (double)i; y = (double)j + frac; }
            if (mp == 0) { xr = x; yr = y; mp = 1; }
            else {
                if (mp == 1) { moveto(xr, yr); mp = 2; }
                drawto(x, y);
            }
        }
        else mp = 0;
        while (khv != 0) {
            switch (++m) {
            case 1: case 11:
                khv = 2; break;
            case 7: case 13:
                khv = 1; break;
            case 6:
                j--; pm1 -= ihsz; pm2 -= ihsz; vp -=imx; khv = 2; break;
            case 2:
                j++; pm1 += ihsz; pm2 += ihsz; vp +=imx; khv = 1; break;
            case 3: case 9:
                j--; pm1 -= ihsz; pm2 -= ihsz; vp -=imx;
            case 14:
                i++; if (++il > 7) { il = 0; pm1++; pm2++; }
                vp++; khv = 2; break;
            case 5: case 15:
                j++; pm1 += ihsz; pm2 += ihsz; vp +=imx;
            case 10:
                i--; if (--il < 0) { il = 7; pm1--; pm2--; }
                vp--; khv = 1; break;
            default:
                khv = 0;
            }
            if ((khv==1) && (*pm1 & cp[il]))
                { *pm1 = *pm1 ^ cp[il]; m = mn[m]; break; }
            if ((khv==2) && (*pm2 & cp[il]))
                { *pm2 = *pm2 ^ cp[il]; m = mn[m]; break; }
        }
    }
    if ((kc==0) && (scs>0.)) {
        px = sxi*x + sxj*y; py = syi*x + syj*y; m = m/4;
        if (kdir == 1) {
            if (szj>0.) m = 4-m; if (szi>0.) m = 6-m;
            if (m>3) m -= 4; px = px+ss[3-m]; py = py+ss[4-m];
            pcstr(sxo+px, syo+py, scs, &fval[4], m*90., -12);
        }
        else {
            m = ((szi>0.)? (m-1): (5-m)); if (szj>0.) m = 6-m;
            if (m>3) m -= 4; px = px+ss[m]; py = py+ss[m+1];
            pcstr(syo+py, sxo+px, scs, &fval[4], m*90., -12);
        }
    }
}

static int contrx(float *v, float undef, float step, int lcapt,
          float vlml, float vlmh, int kln)
{   int j, i, ih, il, ihmx;
    int ispen, klnt, klnp, klnc, kv, ke;
    float vmin, vmax, *vp, *vpe, *vp1, *vp2;
    unsigned char *pbf;
    if (kdir == 0)
        { fprintf(stderr, "contrx: invalid array size\n"); return(-1); }
    if ((pbf=calloc(mbuf,1)) == NULL)
        { fprintf(stderr, "contrx: memory allocation fail\n"); return(-1); }
    ihmx = (imx+9)/8; vpe = v+(imx*jmx);
    for (vp=v; vp<vpe; vp++)
        if (*vp != undef) { vmax = vmin = *vp; break; }
    if (vp >= vpe)
        { fprintf(stderr, "contrx: no valid data exist\n"); return(-1); }
    while (++vp < vpe) {
        if (*vp != undef) {
            if (*vp < vmin) vmin = *vp;
            else if (*vp > vmax) vmax = *vp;
        }
    }
    kv = (int)(vmin/step); if (vmin > 0.) kv++;
    ke = (int)(vmax/step); if (vmax > 0.) ke++;
    dci = szi*MAXLPCM*step; dcj = szj*MAXLPCM*step;
    ispen = qpatr(0,&icol0,&ltyp0,&wd0);
    qpatr(1,&icol1,&klnt,&wd1); qpatr(2,&icol2,&klnt,&wd2);
    klnt = (kln % 10); klnp = (klnt-1)%3+1;
    klnc = ((klnt<4)? klnp: ((klnt<7)? 0: -1));
    if (kln >= 10) klnc += 4; npen = -1;
    while (kv < ke) {
        vc = step*(float)kv; kc = (abs(kv) % lcapt); kv++;
        if ((vc<=vlml) || (vc>=vlmh)) continue;
        if ((kc==0) && (klnc<0)) continue;
        if (lcapt < 0) kc = -kc;
        if (kc == 0) { kpen = klnc; sprintf(fval, fmt, vc); }
        else kpen = klnp;
        for (vp=v,vp1=v+1,vp2=v+imx,j=0; j<jmx; j++) {
            mm1 = pbf+ihsz*(j+1); mm2 = mm1+ihsz*(jmx+1);
            for (i=0,il=0; i<imx; i++,vp++,vp1++,vp2++) {
                if (++il == 8) { il = 0; mm1++; mm2++; }
                if (*vp != undef) {
                    if (((i+1) != imx) && (*vp1 != undef)) {
                        if ((*vp>=vc) && (*vp1<vc) || (*vp<vc) && (*vp1>=vc))
                            *mm1 = *mm1 | cp[il];
                    }
                    if (((j+1) != jmx) && (*vp2 != undef)) {
                        if ((*vp>=vc) && (*vp2<vc) || (*vp<vc) && (*vp2>=vc))
                            *mm2 = *mm2 | cp[il];
                    }
                }
            }
        }
        for (j=0; j<jmx; j++) {
            mm1 = pbf+ihsz*(j+1); mm2 = mm1+ihsz*(jmx+1);
            for (ih=0; ih<ihmx; ih++,mm1++,mm2++) {
                if (*mm1 == 0) continue;
                for (il=0; il<8; il++) {
                    if (*mm1 & cp[il]) drw1cont(v, ih, il, j, 8);
                    if (*mm1 & cp[il])
                        { *mm1 = *mm1 ^ cp[il]; drw1cont(v, ih, il, j, 0); }
                }
            }
            if ((j+1) == jmx) continue;
            mm1 = pbf+ihsz*(j+1); mm2 = mm1+ihsz*(jmx+1);
            for (ih=0; ih<ihmx; ih++,mm1++,mm2++) {
                if (*mm2 == 0) continue;
                for (il=0; il<8; il++) {
                    if (*mm2 & cp[il]) drw1cont(v, ih, il, j, 4);
                    if (*mm2 & cp[il])
                        { *mm2 = *mm2 ^ cp[il]; drw1cont(v, ih, il, j, 12); }
                }
            }
        }
    }
    free(pbf); if (ispen == 0) penatr(0,icol0,ltyp0,wd0); else newpen(ispen); 
    return(0);
}

/*********************/
/*   conts / contx   */
/*********************/

int conts(double xorg, double yorg, double wide, double high,
          int imax, int jmax, double csize)
{   if ((wide==0.) || (high==0.))
        { fprintf(stderr, "conts: area size invalid\n"); return(-1); }
    if ((imax>1) && (jmax>1)) { kdir = 1; imx = imax; jmx = jmax; }
    else if ((imax<-1) && (jmax<-1)) { kdir = 2; imx = -jmax; jmx = -imax; }
    if (kdir == 0)
        { fprintf(stderr, "conts: array size rejected\n"); return(-1); }
    ihsz = (imx+9)/8; mbuf = ihsz*(jmx*2+2);
    if (mbuf > 2048000)
        { fprintf(stderr, "conts: huge array (%d)*(%d)\n", imx, jmx); }
    if (kdir == 1)
        { sxo = xorg; syo = yorg; szi = wide/(imx-1); szj = high/(jmx-1); }
    else
        { sxo = yorg; syo = xorg; szi = high/(imx-1); szj = wide/(jmx-1); }
    sxi = szi; syj = szj; syi = sxj = 0.;
    scs = csize; ss[1] = ((kdir==1)? -(scs/2.): (scs/2.));
    ss[2] = scs/12.; ss[3] = -ss[1]; ss[4] = -ss[2]; ss[0] = ss[4];
    return(0);
}

int contso(double x0, double y0, double x1, double y1,
           double x2, double y2, int imax, int jmax)
{   if ((imax<2) || (jmax<2))
        { fprintf(stderr, "contso: array size rejected\n"); return(-1); }
    kdir = 1; imx = imax; jmx = jmax; ihsz = (imx+9)/8; mbuf = ihsz*(jmx*2+2);
    if (mbuf > 2048000)
        { fprintf(stderr, "contso: huge array (%d)*(%d)\n", imx, jmx); }
    sxo = x0; sxi = (x1-x0)/(imx-1); sxj = (x2-x0)/(jmx-1);
    syo = y0; syi = (y1-y0)/(imx-1); syj = (y2-y0)/(jmx-1);
    szi = sqrt(sxi*sxi + syi*syi); szj = sqrt(sxj*sxj + syj*syj);
    scs = 0.; ss[4] = ss[3] = ss[2] = ss[1] = ss[0] = 0.;
    return(0);
}

int contx(float *v, float undef, int istep, int lcapt,
          int lml, int lmh, int kln)
{   int iret; float step, vlml, vlmh;
    if ((istep<=0) || (lcapt==0) || (lml>lmh) || (kln<0) || (kln>16))
        { fprintf(stderr, "contx: parameter value invalid\n"); return(-1); }
    strcpy(fmt,"%15.0f."); step = (float)istep;
    vlml = (float)lml - 0.5F; vlmh = (float)lmh + 0.5F;
    iret = contrx(v, undef, step, lcapt, vlml, vlmh, kln);
    return(iret);
}

int contr(float *v, float undef, float step, int lcapt,
          float vlml, float vlmh, int kln)
{   int n, iret;
    if ((step<=0.F) || (lcapt==0) || (vlml>=vlmh) || (kln<0) || (kln>16))
        { fprintf(stderr, "contr: parameter value invalid\n"); return(-1); }
    strcpy(fmt,"%15.5f."); sprintf(fval, fmt, step);
    for (n=5; ((n>0) && (fval[n+9]=='0')); n--); fmt[4] = '0' + n;
    iret = contrx(v, undef, step, lcapt, vlml, vlmh, kln);
    return(iret);
}

/*******************************************/
/*   FORTRAN Language Interface Routines   */
/*******************************************/

void conts_(float *xorg, float *yorg, float *wide, float *high,
            int *imax, int *jmax, float *csize)
{ if (conts((double)*xorg, (double)*yorg, (double)*wide,
            (double)*high, *imax, *jmax, (double)*csize) < 0) exit(1); }

void contso_(float *x0, float *y0, float *x1, float *y1,
             float *x2, float *y2, int *imax, int *jmax)
{ if (contso((double)*x0, (double)*y0, (double)*x1, (double)*y1,
             (double)*x2, (double)*y2, *imax, *jmax) < 0) exit(1); }

void contx_(float *v, float *undef, int *istep, int *lcapt,
            int *lml, int *lmh, int *kln)
{ if (contx(v, *undef, *istep, *lcapt, *lml, *lmh, *kln) < 0) exit(1); }

void contr_(float *v, float *undef, float *step, int *lcapt,
             float *vlml, float *vlmh, int *kln)
{ if (contr(v, *undef, *step, *lcapt, *vlml, *vlmh, *kln) < 0) exit(1); }