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

int plot(double xp, double yp, int md);
int where(double *wx, double *wy, double *wfac);

#define RAD     (180./3.14159)

/**************************************************/
/*   p1char / pchars / chtype / psmset / pscdef   */
/**************************************************/

static double ax=0.1, by=0.1, bx=0., ay=0.;
static char csub='\0', csup='\0';
static int ksm=0;

int p1char(char c)
{   static short int bst[95] = {
          0,  8, 16, 24, 36, 52, 61, 65, 69, 73, 79, 83, 89, 91, 96, 98,
        107,112,120,133,137,146,156,161,177,189,199,210,213,217,220,232,
        248,256,268,276,283,290,295,306,312,318,325,331,334,339,343,352,
        359,370,379,391,395,401,404,409,413,418,422,426,428,432,435,437,
        441,455,465,473,483,493,499,511,518,528,539,544,549,559,566,575,
        585,595,601,611,617,624,627,632,636,642,646,653,655,662,668 };
    static unsigned char buf[668] = {
        0xc4,0x59,0x49,0x44, 0xc3,0x32,0x42,0x43, 0xb7,0x39,0x49,0x37,
        0xd7,0x59,0x69,0x57, 0xb2,0x39,0xa7,0x67, 0xd9,0x52,0xe4,0x24,
        0xa4,0x33,0x53,0x64, 0x65,0x26,0x27,0x38, 0x58,0x67,0xc9,0x42,
        0xa2,0x69,0xc8,0x38, 0x27,0x26,0x36,0x47, 0x48,0xe5,0x55,0x44,
        0x43,0x53,0x64,0x65, 0xe2,0x28,0x39,0x48, 0x25,0x23,0x32,0x42,
        0x64,0xc6,0x39,0x59, 0x46,0xd9,0x48,0x43, 0x52,0xb9,0x48,0x43,
        0x32,0xa4,0x66,0xe4, 0x26,0xc7,0x43,0xa5, 0x65,0xc7,0x43,0xc3,
        0x33,0x34,0x44,0x43, 0x32,0xa5,0x65,0xb2, 0x33,0x43,0x42,0x32,
        0xa2,0x69,0xd2,0x42, 0x33,0x37,0x48,0x58, 0x67,0x63,0x52,0xb7,
        0x48,0x42,0xb2,0x52, 0xa7,0x38,0x58,0x67, 0x66,0x23,0x22,0x62,
        0xa7,0x38,0x58,0x67, 0x66,0x55,0x35,0xd5, 0x64,0x63,0x52,0x32,
        0x23,0xd2,0x58,0x24, 0x64,0xa3,0x32,0x52, 0x63,0x64,0x55,0x25,
        0x28,0x68,0xe8,0x37, 0x25,0x23,0x32,0x52, 0x63,0x64,0x55,0x25,
        0xa7,0x28,0x68,0x43, 0x42,0xb5,0x26,0x27, 0x38,0x58,0x67,0x66,
        0x55,0x35,0x24,0x23, 0x32,0x52,0x63,0x64, 0x55,0xa3,0x32,0x52,
        0x63,0x67,0x58,0x38, 0x27,0x26,0x35,0x55, 0x66,0xb6,0x37,0x47,
        0x46,0x36,0xb3,0x34, 0x44,0x43,0x33,0xc6, 0x36,0x37,0x47,0x46,
        0xc3,0x33,0x34,0x44, 0x43,0x32,0xe7,0x25, 0x63,0xa6,0x66,0xe4,
        0x24,0xa7,0x65,0x23, 0xa7,0x28,0x39,0x59, 0x68,0x67,0x45,0x44,
        0xc3,0x42,0x32,0x43, 0xd5,0x57,0x47,0x36, 0x34,0x44,0x55,0x64,
        0x67,0x58,0x38,0x27, 0x24,0x33,0x53,0x64, 0xa2,0x28,0x39,0x59,
        0x68,0x62,0xa5,0x65, 0xa2,0x29,0x59,0x68, 0x67,0x56,0x26,0xd6,
        0x65,0x63,0x52,0x22, 0xe8,0x59,0x39,0x28, 0x23,0x32,0x52,0x63,
        0xa2,0x29,0x59,0x68, 0x63,0x52,0x22,0xa2, 0x29,0x69,0xd6,0x26,
        0xa2,0x62,0xa2,0x29, 0x69,0xd6,0x26,0xc5, 0x75,0xe5,0x63,0x52,
        0x32,0x23,0x28,0x39, 0x59,0x68,0xa2,0x29, 0xa6,0x66,0xe9,0x62,
        0xb9,0x59,0xc9,0x42, 0xb2,0x52,0xa3,0x32, 0x42,0x53,0x59,0xc9,
        0x69,0xa2,0x29,0xe9, 0x25,0xb6,0x62,0xa9, 0x22,0x62,0xa2,0x29,
        0x45,0x69,0x62,0xa2, 0x29,0x62,0x69,0xb2, 0x23,0x28,0x39,0x59,
        0x68,0x63,0x52,0x32, 0xa2,0x29,0x59,0x68, 0x67,0x56,0x26,0xb2,
        0x23,0x28,0x39,0x59, 0x68,0x63,0x52,0x32, 0xc4,0x62,0xa2,0x29,
        0x59,0x68,0x67,0x56, 0x26,0xc6,0x62,0xa3, 0x32,0x52,0x63,0x65,
        0x56,0x36,0x27,0x28, 0x39,0x59,0x68,0xc2, 0x49,0xa9,0x69,0xa9,
        0x23,0x32,0x52,0x63, 0x69,0xa9,0x42,0x69, 0xa9,0x32,0x46,0x52,
        0x69,0xa2,0x69,0xa9, 0x62,0xc2,0x46,0x29, 0xc6,0x69,0xa9,0x69,
        0x22,0x62,0xe9,0x39, 0x32,0x62,0xa9,0x62, 0xa2,0x52,0x59,0x29,
        0xa8,0x4a,0x68,0x91, 0x71,0xc9,0x39,0x58, 0x49,0xa6,0x37,0x57,
        0x66,0x64,0x62,0xe4, 0x55,0x35,0x24,0x23, 0x32,0x52,0x63,0xa2,
        0x29,0xa5,0x36,0x56, 0x65,0x63,0x52,0x32, 0x23,0xe6,0x57,0x37,
        0x26,0x23,0x32,0x52, 0x63,0xe9,0x62,0xe3, 0x52,0x32,0x23,0x25,
        0x36,0x56,0x65,0xa4, 0x65,0x66,0x57,0x37, 0x26,0x23,0x32,0x52,
        0x63,0xc2,0x48,0x59, 0x69,0xa6,0x66,0xa1, 0x30,0x50,0x61,0x66,
        0x57,0x37,0x26,0x24, 0x33,0x53,0x64,0xa2, 0x29,0xa5,0x36,0x56,
        0x65,0x62,0xb8,0x39, 0x49,0x48,0x38,0xb6, 0x46,0x42,0xb2,0x52,
        0xc8,0x49,0x59,0x58, 0x48,0xc6,0x56,0x51, 0x40,0x30,0x21,0xa9,
        0x22,0xe6,0x24,0x62, 0xb9,0x49,0x42,0xb2, 0x52,0xa2,0x27,0xa6,
        0x37,0x46,0x57,0x66, 0x62,0xc6,0x42,0xa2, 0x27,0xa6,0x37,0x57,
        0x66,0x62,0xb2,0x23, 0x26,0x37,0x57,0x66, 0x63,0x52,0x32,0xa0,
        0x27,0xa6,0x37,0x57, 0x66,0x64,0x53,0x33, 0x24,0xe0,0x67,0xe6,
        0x57,0x37,0x26,0x24, 0x33,0x53,0x64,0xa2, 0x27,0xa6,0x37,0x57,
        0x66,0xa3,0x32,0x52, 0x63,0x64,0x25,0x26, 0x37,0x57,0x66,0xa6,
        0x66,0xc9,0x43,0x52, 0x62,0xa7,0x23,0x32, 0x52,0x63,0xe2,0x67,
        0xa7,0x42,0x67,0xa7, 0x32,0x45,0x52,0x67, 0xa2,0x67,0xa7,0x62,
        0x90,0x20,0x31,0x67, 0xa7,0x43,0xa7,0x67, 0x22,0x62,0xda,0x49,
        0x47,0x36,0x45,0x43, 0x52,0xc9,0x42,0xba, 0x49,0x47,0x56,0x45,
        0x43,0x32,0x98,0x2a, 0x3a,0x58,0x68,0x7a };
    int k, i, ip;
    double x0, y0, x1, y1, an, bn, w;
    if ((c<' ') || (c>'~')) return(1);
    where(&x0, &y0, &w);
    if ((k=(c-' ')) != 0) {
        if ((c==csub) && (ksm!=2)) { ksm = 1 - ksm; return(0); }
        if ((c==csup) && (ksm!=1)) { ksm = 2 - ksm; return(0); }
        for (i=bst[k-1]; i<bst[k]; i++) {
            ip = ((buf[i]>>7) & 1) + 2;
            an = (double)(int)(((buf[i]>>4) & 7) - 2);
            bn = (double)(int)((buf[i] & 0x0f) - 2);
            if (ksm == 1) { an *= 0.7; bn = bn*0.7 - 0.6; }
            if (ksm == 2) { an *= 0.7; bn = bn*0.7 + 4.3; }
            x1 = x0 + an*ax + bn*bx;
            y1 = y0 + an*ay + bn*by;
            if (plot(x1,y1,ip) < 0) return(-1);
        }
    }
    x1 = x0 + ax*6.; y1 = y0 + ay*6.;
    if (plot(x1,y1,3) < 0) return(-1);
    return(0);
}

int pchars(char *bcd)
{   char c;
    while ((c=*bcd) != '\0')
        { bcd++; if (p1char(c) < 0) return(-1); }
    return(0);
}

int chtype(double ang, double slope, double high, double ratio)
{   double a, s, ar, br;
    a = ang/RAD; s = slope/RAD;
    ar = (high/7.) * ratio;
    br = (high/7.) / cos(s);
    ax = ar * cos(a); by = br * cos(s-a);
    ay = ar * sin(a); bx = br * sin(s-a);
    return(0);
}

int psmset(int ksms)
{   if ((ksms<0) || (ksms>2))
        { fprintf(stderr, "psmset: illegal pchar-mode\n"); return(1); }
    else { ksm = ksms; return(0); }
}

int pscdef(char csb, char csp)
{   csub = csb; csup = csp; return(0); }

/**************/
/*   pcsymb   */
/**************/

int pcsymb(int ics, double size, double ang)
{   static short int bst[40] = {
           6,  13,  24,   0,  29,  35,  41,  47,  52,  59,
          65,  79,  91,  96,  99, 106, 108, 110, 116, 122,
          -6, -13, -24, 128, 137, -35, 146, 151, 156, 161,
         -65, 166, 183, 192, -99, 195, 201,-110,-116,-122 };
    static unsigned char buf[207] = {
        0x35,0x31,0x33,0x13, 0x53,0x80,0x35,0x15, 0x11,0x51,0x55,0x35,
        0x80,0x35,0x25,0x14, 0x12,0x21,0x41,0x52, 0x54,0x45,0x35,0x80,
        0x35,0x12,0x52,0x35, 0x80,0x55,0x11,0x33, 0x15,0x51,0x80,0x35,
        0x13,0x31,0x53,0x35, 0x80,0x35,0x13,0x53, 0x35,0x31,0x80,0x11,
        0x55,0x15,0x51,0x80, 0x55,0x15,0x55,0x11, 0x51,0x11,0x80,0x55,
        0x33,0x15,0x33,0x31, 0x80,0x44,0x55,0x44, 0x24,0x15,0x24,0x22,
        0x11,0x22,0x42,0x51, 0x42,0x44,0x80,0x35, 0x31,0x33,0x15,0x51,
        0x33,0x13,0x53,0x33, 0x11,0x55,0x80,0x55, 0x15,0x51,0x11,0x80,
        0x35,0x31,0x80,0x36, 0x11,0x64,0x04,0x51, 0x36,0x80,0x63,0x80,
        0x36,0x80,0x35,0x23, 0x31,0x43,0x35,0x80, 0x34,0x13,0x32,0x53,
        0x34,0x80,0x34,0x14, 0x31,0x54,0x34,0x80, 0x3c,0x35,0x2b,0x13,
        0x3a,0x31,0x4b,0x53, 0x88,0x4c,0x55,0x2c, 0x15,0x2a,0x11,0x4a,
        0x51,0x88,0x30,0x33, 0x22,0x42,0x80,0x36, 0x33,0x44,0x24,0x80,
        0x63,0x33,0x42,0x44, 0x80,0x03,0x33,0x24, 0x22,0x80,0x3c,0x35,
        0x2c,0x15,0x2b,0x13, 0x2a,0x11,0x3a,0x31, 0x4a,0x51,0x4b,0x53,
        0x4c,0x55,0x88,0x4c, 0x55,0x15,0x24,0x2a, 0x11,0x51,0x42,0x88,
        0x53,0x13,0x80,0x63, 0x54,0x63,0x52,0x63, 0x80,0x36,0x25,0x36,
        0x45,0x36,0x80 };
    int ip, n, na;
    double x0, y0, x1, y1, dx, dy, fa, fb, w;
    if ((ics<0) || (ics>39)) return(1);
    where(&x0, &y0, &w);
    dx = (size/4.) * cos(ang/RAD);
    dy = (size/4.) * sin(ang/RAD);
    n = (na=bst[ics]); ip = 2; if (na < 0) { n = -na; ip = 3; }
    while ((buf[n] & 0x80) == 0) {
        if (buf[n] & 8) ip = 3;
        fa = (double)(int)(((buf[n]>>4) & 7) - 3);
        fb = (double)(int)((buf[n] & 7) - 3);
        x1 = x0 + fa*dx - fb*dy;
        y1 = y0 + fa*dy + fb*dx;
        if (plot(x1,y1,ip) < 0) return(-1);
        ip = 2; n++;
    }
    if ((buf[n] & 8) || (na < 0)) ip = 3;
    return( plot(x0,y0,ip) );
}

/***********************/
/*   symbol / number   */
/***********************/

int symbol(double xp, double yp, double high, char *text, double ang, int nc)
{   double w, wf;
    int k;
    if (xp == 999.) { where(&xp, &w, &wf); if (yp == 999.) yp = w; }
    else if (yp == 999.) where(&w, &yp, &wf);
    if (nc <= 0) {
        if (nc < -1) k = plot(xp,yp,2); else k = plot(xp,yp,3);
        if (k < 0) return(-1);
        pcsymb((int)*text, high, ang);
    }
    else {
        if (plot(xp,yp,3) < 0) return(-1);
        chtype(ang, 0., high, 7./6.);
        while(nc > 0) { p1char(*text); text++; nc--; }
    }
    return(0);
}

int number(double xp, double yp, double high, double fv, double ang, int ndec)
{   double w, wf;
    char str[32], fmt[7], *s;
    if (ndec < -9) ndec = -9;
    if (ndec > +9) ndec = +9;
    if (ndec > 0) sprintf(fmt, "%%%2d.%df", 22+ndec, ndec);
    else if (ndec == 0) sprintf(fmt, "%%21.0f.");
    else { while (++ndec < 0) fv /= 10.; sprintf(fmt, "%%21.0f"); }
    if (sprintf(str,fmt,fv) < 0)
        { fprintf(stderr, "number: too long expression\n"); return(2); }
    for (s=str; *s==' '; s++) if (*s == '\0') return(1);
    if (xp == 999.) { where(&xp, &w, &wf); if (yp == 999.) yp = w; }
    else if (yp == 999.) where(&w, &yp, &wf);
    if (plot(xp,yp,3) < 0) return(-1);
    chtype(ang, 0., high*7./6., 1.); pchars(s);
    return(0);
}

/********************/
/*   scale / line   */
/********************/

int scale(float va[], double axlen, int npt, int inc)
{   int n, ipt, ke, kf;
    double vl, vh, vd, ul, uh, uf;
    if ((npt<=1) || (inc<=0) || (axlen<2.))
        { fprintf(stderr, "scale: parameter value invalid\n"); return(-1); }
    vl = vh = va[0];
    for (n=inc,ipt=1; ipt<npt; ipt++,n+=inc)
        { if (va[n] < vl) vl = va[n]; if (va[n] > vh) vh = va[n]; }
    if (vl == vh)
        { fprintf(stderr, "scale: identical data values\n"); return(-1); }
    vd = (vh-vl)/axlen; ke = (int)log10(vd); if (vd < 1.) ke--;
    kf = (int)(vd/pow(10.,(double)ke));
    do {
        kf++;
        if (kf > 8) { kf = 1; ke++; }
        else if (kf > 5) kf = 8;
        else if (kf == 3) kf = 4;
        vd = (double)kf*pow(10.,(double)ke);
        ul = vl/vd; uh = vh/vd;
        uf = (double)(int)ul; if (ul < uf) uf = uf-1.;
    } while ((uh-uf) > axlen);
    va[n] = uf*vd; va[n+inc] = vd;
    return(0);
}

int line(float xa[], float ya[], int npt, int inc, int ltype, int mark)
{   int nf, n, ipt;
    double xf, yf, xd, yd, x, y;
    if ((npt<=0) || (inc<=0))
        { fprintf(stderr, "line: parameter value invalid\n"); return(-1); }
    nf = npt*inc; xf = xa[nf]; yf = ya[nf]; xd = xa[nf+inc]; yd = ya[nf+inc];
    if (xd == 0.) xd = 1.;
    if (yd == 0.) yd = 1.;
    x = (xa[0]-xf)/xd; y = (ya[0]-yf)/yd; if (plot(x,y,3) < 0) return(-1);
    if (ltype != 0) pcsymb(mark, 0.21, 0.);
    for (n=inc,ipt=1; ipt<npt; ipt++,n+=inc) {
        x = (xa[n]-xf)/xd; y = (ya[n]-yf)/yd;
        if (ltype < 0) plot(x, y, 3); else plot(x, y, 2);
        if (ltype != 0) if ((ipt%ltype) == 0) pcsymb(mark, 0.21, 0.);
    }
    return(0);
}

/************/
/*   axis   */
/************/

int axis(double xp, double yp, char *text, int ltx,
         double axlen, double angle, double firstv, double deltav)
{   int ke, i;
    double bt, bv, bs, ca, sa, fe, as, xs, ys, xt, yt, vf;
    if ((ltx==0) || (axlen<2.))
        { fprintf(stderr, "axis: parameter value invalid\n"); return(-1); }
    if (ltx > 0)     { bt =  0.735; bv =  0.298; bs =  0.178; }
    else { ltx = -ltx; bt = -1.015; bv = -0.508; bs = -0.178; }
    ca = cos(angle/RAD); sa = sin(angle/RAD);
    if (deltav == 0.) deltav = 1.;
    ke = (int)log10(fabs(deltav*1.00001));
    if (ke > 0) ke--;
    else if (ke < 0) ke++;
    fe = pow(10.,(double)ke); firstv /= fe; deltav /= fe;
    as = axlen/2.-0.14*(double)ltx; if (ke != 0) as -= 0.88;
    xs = xp+ca*as-sa*bt; ys = yp+sa*as+ca*bt;
    if (plot(xs,ys,3) < 0) return(-1);
    chtype(angle, 0., 0.28, 7./6.);
    for(i=0; i<ltx; i++,text++) p1char(*text);
    if (ke != 0) {
        pchars("  *10"); as = (ltx+5)*0.28;
        xs = xs+ca*as-sa*0.18; ys = ys+sa*as+ca*0.18;
        number(xs, ys, 0.18, (double)ke, angle, -1);
    }
    xt = xp; yt = yp; plot(xt, yt, 3);
    xs = xt-sa*bs; ys = yt+ca*bs; plot(xs, ys, 2);
    xs = xt-ca*0.28-sa*bv; ys = yt-sa*0.28+ca*bv; vf = firstv;
    number(xs, ys, 0.21, vf, angle, 2);
    for (i=1; i<=(int)axlen; i++) {
        xt += ca; yt += sa; vf += deltav; plot(xt, yt, 3);
        xs = xt-sa*bs; ys = yt+ca*bs; plot(xs, ys, 2);
        if ((i%2) == 0) {
            xs = xt-ca*0.49-sa*bv; ys = yt-sa*0.49+ca*bv;
            number(xs, ys, 0.21, vf, angle, 2);
        }
    }
    xs = xp+ca*axlen; ys = yp+sa*axlen; plot(xs, ys, 3); plot(xp, yp, 2);
    return(0);
}

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

#define W       double

/*********************************************/
/*   csymbl_ / gsymbl_ / chtype_ / pchar_    */
/*    pnorm_ / psubsc_ / psupsc_ / pshift_   */
/*********************************************/

void csymbl_(float *xs, float *ys, float *sz, int *ics, float *ang, int *mp)
{   int ip;
    ip = abs(*mp)%2 + 2;
    if (plot((W)*xs,(W)*ys,ip) < 0) exit(1);
    if (pcsymb(*ics,(W)*sz,(W)*ang) < 0) exit(1);
}

void gsymbl_(float *xs, float *ys, float *high, char *s, float *ang, int *nc)
{   int i;
    plot((W)*xs, (W)*ys, 3); chtype((W)*ang, 0., (W)*high, 1.);
    for (i=0; i<(*nc); i++,s++)
        { if (p1char(*s) < 0) exit(1); }
}

void chtype_(float *ang, float *slope, float *high, float *ratio)
{   chtype((W)*ang, (W)*slope, (W)*high, (W)*ratio); }

void pchar_(char *s, int *nc, float *xs, float *ys, int *kk)
{   int i;
    if ((*kk) != 0) plot((W)*xs, (W)*ys, 3);
    for (i=0; i<(*nc); i++,s++)
        { if (p1char(*s) < 0) exit(1); }
}

void pnorm_(void)
{   psmset(0); }

void psubsc_(void)
{   psmset(1); }

void psupsc_(void)
{   psmset(2); }

void pshift_(char cs[2])
{   pscdef(cs[0], cs[1]); }

/**************************************************/
/*   symbol_ / number_ / scale_ / line_ / axis_   */
/**************************************************/

void symbol_(float *xp, float *yp, float *high, char *text,
             float *ang, int *nc)
{ if (symbol((W)*xp,(W)*yp,(W)*high,text,(W)*ang,*nc) < 0) exit(1); }

void number_(float *xp, float *yp, float *high, float *fv,
             float *ang, int *ndec)
{ if (number((W)*xp,(W)*yp,(W)*high,(W)*fv,(W)*ang,*ndec) < 0) exit(1); }

void scale_(float va[], float *axlen, int *npt, int *inc)
{ if (scale(va,(W)*axlen,*npt,*inc) < 0) exit(1); }

void line_(float xa[], float ya[], int *npt, int *inc, int *ltype, int *mark)
{ if (line(xa,ya,*npt,*inc,*ltype,*mark) < 0) exit(1); }

void axis_(float *xp, float *yp, char *text, int *ltx,
           float *axlen, float *angle, float *firstv, float *deltav)
{ if (axis((W)*xp,(W)*yp,text,*ltx,(W)*axlen,
           (W)*angle,(W)*firstv,(W)*deltav) < 0) exit(1); }