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

#define DPI     600
#define DPCM    (DPI/2.54)
#define DP300um 7           /* = (int)(DPCM*0.03) */
#define SCL     "0.12"      /* = (72./DPI) */
#define MAXPEN  8
#define RAD     (180./3.14159)
#define MAXLVL  50

/***************************************************************/
/*   psopn / pscls / plots / plote / plot / clspath / scisor   */
/*     / factor / penatr / newpen / qpatr / qpsfp / where      */
/***************************************************************/

static FILE *fp;
static double xorg, yorg, px, py, fs, xr, yr, xa, xb, ya, yb, penw;
static double wpen[MAXPEN] = { 0.01, 0.03, 0.06, 0.1, 0.01, 0.03, 0.06, 0.1 };
static int lpent[MAXPEN] = { 0, 0, 0, 0, 1, 1, 1, 1 };
static int lpenc[MAXPEN] = { 0, 0, 0, 0, 0, 0, 0, 0 };
static int lpc, ixs, iys;
static int kplt=0, keps;
static int ipen, kpatr, lpt, kclp, npage, nseq;
static int initx, inity, irot;
static int bbx0=0, bby0=0, bbx1=0, bby1=0, kbox=0;

static void setpatr()
{   int ipw, ids, r, g, b;
    if (lpc == 0) fprintf(fp, "0 T ");
    else if (lpc < 0) fprintf(fp, "%.3f T ", -lpc/255.);
    else {
        r = lpc/0x010000; g = (lpc/256)%256; b = lpc%256;
        fprintf(fp,"%.3f %.3f %.3f RGB ",r/255.,g/255.,b/255.);
    }
    ipw = (int)(penw*DPCM);
    fprintf(fp, "%d W ", ipw);
    if (ipw < DP300um) ids = DP300um*2+2; else ids = ipw*2+2;
    switch (lpt) {
    case 0:
        fprintf(fp, "[] 0 D\n"); break;
    case 1:
        fprintf(fp, "[%d %d] %d D\n", ids*2, ids, ids); break;
    case 2:
        fprintf(fp, "[1 %d] 0 D\n", ipw*3+2); break;
    case 3:
        fprintf(fp, "[%d %d 1 %d] %d D\n", ids*2-1, ids, ids, ids-1);
    }
    kpatr = 1;
}

int epsbox(double obbx, double obby, double sbbx, double sbby)
{   if ((obbx<0.) || (obby<0.) || (sbbx<0.) || (sbby<0.))
        { fprintf(stderr, "epsbox: illegal bounding box\n"); return(-1); }
    if ((sbbx>19.7) || (sbby>28.0))
        { fprintf(stderr, "epsbox: too large bounding box\n"); return(-1); }
    bbx0 = (int)(obbx*72./2.54) + 20;
    bby0 = (int)(obby*72./2.54) + 25;
    bbx1 = (int)((obbx+sbbx)*72./2.54) + 21;
    bby1 = (int)((obby+sbby)*72./2.54) + 26;
    kbox = 1; return(0);
}

int psopn(char *pfile, char *sheet)
{   char sht[4]="---";
    int h[2][6] = { 3168,2232,1584,1116,792,558, 3888,2736,1944,1368,972,684 };
    int k, i, wh, ww=0;
    time_t dtm; char *adtm, *user;
    time(&dtm); adtm = asctime(localtime(&dtm));
    *(adtm+strlen(adtm)-1) = '\0';
    if ((user=getenv("LOGNAME")) == NULL) user = getenv("USERNAME");
    if (kplt != 0) {
        fprintf(stderr, "psopn: psplot already open\n");
        return (-1);
    }
    if ((pfile==NULL) || (*pfile=='\0')) fp = stdout;
    else {
        if ((fp=fopen(pfile,"w")) == NULL) {
            fprintf(stderr, "psopn: unable to open psfile (%s)\n", pfile);
            return (-1);
        }
    }
    initx = 20; inity = 25; irot = 0; keps = 0;
    if (sheet != NULL) {
        if ((*sheet=='a') || (*sheet=='A')) { sht[0] = 'A'; k = 0; }
        if ((*sheet=='b') || (*sheet=='B')) { sht[0] = 'B'; k = 1; }
        if ((sht[0]=='A') || (sht[0]=='B')) {
            sheet++; i = *sheet-'0';
            if ((i>=0) && (i<=4)) {
                sht[1] = *sheet; sheet++; wh = h[k][i]; ww = h[k][i+1];
                if ((*sheet=='p') || (*sheet=='P')) sht[2] = 'P';
                else if ((*sheet=='l') || (*sheet=='L'))
                    { sht[2] = 'L'; initx = ww+20; irot = 90; }
                if ((sht[2]=='P') || (sht[2]=='L')) {
                    sheet++;
                    if ((*sheet=='e') || (*sheet=='E')) keps = 1;
                }
            }
        }
    }
    fprintf(fp, "%%!PS-Adobe-3.0\n");
    fprintf(fp, "%%%%Title: (%s) psplot 2009.06 coded by T. Nakatsuka\n", sht);
    if (user != NULL) fprintf(fp, "%%%%Creator: (%s)\n", user);
    fprintf(fp, "%%%%CreationDate: (%s)\n", adtm);
    if (keps == 0) fprintf(fp, "%%%%Pages: (atend)\n");
    if (kbox == 1)
        fprintf(fp, "%%%%BoundingBox: %d %d %d %d\n", bbx0,bby0,bbx1,bby1);
    if ((kbox==0) && (ww != 0))
        fprintf(fp, "%%%%BoundingBox: 15 20 %d %d\n", ww+25,wh+30);
    fprintf(fp, "%%%%EndComments\n\n");
    fprintf(fp, "%%%%BeginProlog\n");
    fprintf(fp, "22 dict begin\n");
    fprintf(fp, "/W /setlinewidth load def\n");
    fprintf(fp, "/D /setdash load def\n");
    fprintf(fp, "/T /setgray load def\n");
    fprintf(fp, "/RGB /setrgbcolor load def\n");
    fprintf(fp, "/S /stroke load def\n");
    fprintf(fp, "/N /newpath load def\n");
    fprintf(fp, "/M /moveto load def\n");
    fprintf(fp, "/L /lineto load def\n");
    fprintf(fp, "/J /setlinejoin load def\n");
    fprintf(fp, "/C /closepath load def\n");
    fprintf(fp, "/A /arc load def\n");
    fprintf(fp, "/Count 0 def\n");
    fprintf(fp, "/Byte 1 string def\n");
    fprintf(fp, "/Color 3 string def\n");
    fprintf(fp, "/Pixels 768 string def\n");
    fprintf(fp, "/Get1B {currentfile Byte readhexstring pop} bind def\n");
    fprintf(fp, "/GetCol {currentfile Color readhexstring pop} bind def\n");
    fprintf(fp, "/GrayPacket {/Count Get1B 0 get def Get1B pop\n");
    fprintf(fp, "   0 1 Count {Pixels exch Byte putinterval} for\n");
    fprintf(fp, "   Pixels 0 Count 1 add getinterval} bind def\n");
    fprintf(fp, "/ColorPacket {/Count Get1B 0 get 3 mul def GetCol pop\n");
    fprintf(fp, "   0 3 Count {Pixels exch Color putinterval} for\n");
    fprintf(fp, "   Pixels 0 Count 3 add getinterval} bind def\n");
    fprintf(fp, "/Cmap %d string def\n", (MAXLVL+3)*3);
    fprintf(fp, "/CmapPacket {/Count Get1B 0 get 3 mul def\n");
    fprintf(fp, "   0 3 Count {Cmap exch GetCol putinterval} for} bind def\n");
    fprintf(fp, "/CcodePacket {/Count Get1B 0 get 3 mul def  Color 0\n");
    fprintf(fp, "   Cmap Get1B 0 get 3 mul 3 getinterval putinterval\n");
    fprintf(fp, "   0 3 Count {Pixels exch Color putinterval} for\n");
    fprintf(fp, "   Pixels 0 Count 3 add getinterval} bind def\n");
    fprintf(fp, "%%%%EndProlog\n\n");
    ipen = 0; lpc = 0; lpt = 0; penw = 0.; fs = DPCM; npage = 0;
    kplt = 1; return(0);
}

int plots(double xo, double yo)
{   if (kplt != 1) {
        if ((kplt==2) || (kplt==3))
             fprintf(stderr, "plots: illegal call sequence\n");
        else fprintf(stderr, "plots: psplot not open\n");
        return(-1);
    }
    xorg = xo*fs; yorg = yo*fs; px = 0.; py = 0.;
    xr = xorg; yr = yorg; npage++;
    if (keps == 0) { fprintf(fp, "%%%%Page: %d %d\n", npage, npage); }
    else if (keps == 1) keps = 2;
    else { fprintf(stderr, "plots: multi-page on EPS\n"); return(-1); }
    fprintf(fp, "gsave %d %d translate %d rotate", initx,inity,irot);
    fprintf(fp, " %s %s scale 1 setlinecap 1 J\n", SCL, SCL);
    kpatr = 0; ipen = 1; lpc = lpenc[0]; lpt = lpent[0]; penw = wpen[0];
    kclp = 0; kplt = 2; return(0);
}

int plote(void)
{   if (kplt == 3) { fprintf(fp, " S\n"); kplt = 2; }
    if (kplt != 2) {
        if (kplt == 1) fprintf(stderr, "plote: illegal call sequence\n");
        else           fprintf(stderr, "plote: psplot not open\n");
        return(-1);
    }
    fprintf(fp, "grestore showpage\n\n");
    kplt = 1; ipen = 0; lpc = 0; lpt = 0; penw = 0.;
    return(0);
}

int pscls(void)
{   if ((kplt<1) || (kplt>3)) {
        fprintf(stderr, "pscls: psplot not open\n");
        return(-1);
    }
    if (kplt >= 2) plote();
    fprintf(fp, "%%%%Trailer\n");
    fprintf(fp, "end\n");
    if (keps == 0) { fprintf(fp, "%%%%Pages: %d\n", npage); }
    fprintf(fp, "%%%%EOF\n");
    if (fp != stdout) fclose(fp);
    kplt = 0; keps = 0; return(0);
}

int plot(double xp, double yp, int md)
{   double x0, y0, x1, y1;
    int ix, iy;
    if ((kplt<2) || (kplt>3)) {
        if (kplt == 1) fprintf(stderr, "plot: illegal call sequence\n");
        else           fprintf(stderr, "plot: psplot not open\n");
        return(-1);
    }
    px = xp; x0 = xr; xr = xorg + px*fs;
    py = yp; y0 = yr; yr = yorg + py*fs;
    if (md < 0) { xorg = xr; yorg = yr; px = 0.; py = 0.; md = -md; }
    if (md != 2) {
        if (kplt == 3) { fprintf(fp, " S\n"); kplt = 2; }
        return(0);
    }
    if (kclp == 0) {
        ix = (int)xr; iy = (int)yr;
        if (kplt == 2) {
            if (kpatr == 0) setpatr();
            fprintf(fp, "N %d %d M %d %d L", (int)x0, (int)y0, ix, iy);
            nseq = 2; kplt = 3; ixs = ix; iys = iy;
        }
        else {
            if ((ix!=ixs) || (iy!=iys)) {
                if (nseq > 5) { fputc('\n', fp); nseq = 0; }
                fprintf(fp, " %d %d L", ix, iy); nseq++; ixs = ix; iys = iy;
            }
        }
        return(0);
    }
    if (kplt == 2) {
        if (x0 < xa) {
            if (xr < xa) return(0);
            else { y0 += (yr-y0)*(xa-x0)/(xr-x0); x0 = xa; }
        }
        else if (x0 > xb) {
            if (xr > xb) return(0);
            else { y0 += (yr-y0)*(x0-xb)/(x0-xr); x0 = xb; }
        }
        if (y0 < ya) {
            if (yr < ya) return(0);
            else { x0 += (xr-x0)*(ya-y0)/(yr-y0); y0 = ya; }
        }
        else if (y0 > yb) {
            if (yr > yb) return(0);
            else { x0 += (xr-x0)*(y0-yb)/(y0-yr); y0 = yb; }
        }
        if ((x0 < xa) && (xr < xa)) return(0);
        if ((x0 > xb) && (xr > xb)) return(0);
    }
    x1 = xr; y1 = yr; kclp = 1;
    if (x1 < xa) { y1 += (y0-y1)*(xa-x1)/(x0-x1); x1 = xa; kclp = 2; }
    if (x1 > xb) { y1 += (y0-y1)*(x1-xb)/(x1-x0); x1 = xb; kclp = 2; }
    if (y1 < ya) { x1 += (x0-x1)*(ya-y1)/(y0-y1); y1 = ya; kclp = 2; }
    if (y1 > yb) { x1 += (x0-x1)*(y1-yb)/(y1-y0); y1 = yb; kclp = 2; }
    ix = (int)x1; iy = (int)y1;
    if (kplt == 2) {
        if (kpatr == 0) setpatr();
        fprintf(fp, "N %d %d M %d %d L", (int)x0, (int)y0, ix, iy);
        nseq = 2; kplt = 3; ixs = ix; iys = iy;
    }
    else {
        if ((ix!=ixs) || (iy!=iys)) {
            if (nseq > 5) { fputc('\n', fp); nseq = 0; }
            fprintf(fp, " %d %d L", ix, iy); nseq++; ixs = ix; iys = iy;
        }
    }
    if (kclp == 2) { fprintf(fp, " S\n"); kplt = 2; }
    return(0);
}

int clspath(int join)
{   if (kplt != 3) {
        if (kplt == 2) fprintf(stderr, "clspath: current path not open\n");
        if (kplt == 1) fprintf(stderr, "clspath: illegal call sequence\n");
        else           fprintf(stderr, "clspath: psplot not open\n");
        return(-1);
    }
    if ((join==0) || (join==2)) fprintf(fp, " C %d J S 1 J\n", join);
    else fprintf(fp, " C S\n");
    kplt = 2; return(0);
}

int scisor(double xlm, double ylm, double wlm, double hlm)
{   if ((kplt<2) || (kplt>3)) {
        if (kplt == 1) fprintf(stderr, "scisor: illegal call sequence\n");
        else           fprintf(stderr, "scisor: psplot not open\n");
        return(-1);
    }
    if ((wlm==0.) && (hlm==0.)) { kclp = 0; return(0); }
    else if ((wlm<=0.) || (hlm<=0.))
        { fprintf(stderr, "scisor: not accepted\n"); return(1); }
    else {
        xa = xorg + xlm*fs; xb = xorg + (xlm+wlm)*fs;
        ya = yorg + ylm*fs; yb = yorg + (ylm+hlm)*fs;
        if (kplt == 3) { fprintf(fp, " S\n"); kplt = 2; }
        kclp = 1; return(0);
    }
}

int factor(double fac)
{   if ((kplt<1) || (kplt>3))
        { fprintf(stderr, "factor: psplot not open\n"); return(-1); }
    fs = fac * DPCM; return(0);
}

int penatr(int jpen, int icol, int ityp, double thick)
{   if ((jpen<0) || (jpen>MAXPEN))
        { fprintf(stderr, "penatr: illegal pen-number\n"); return(1); }
    if ((icol<-255) || (icol>0xffffff))
        { fprintf(stderr, "penatr: illegal line-color\n"); return(1); }
    if ((ityp<0) || (ityp>3))
        { fprintf(stderr, "penatr: illegal line-type\n"); return(1); }
    if (thick < 0.)
        { fprintf(stderr, "penatr: illegal thickness\n"); return(1); }
    if (thick < 0.01) thick = 0.01;
    if (jpen == 0) ipen = 0;
    else { lpenc[jpen-1] = icol; lpent[jpen-1] = ityp; wpen[jpen-1] = thick; }
    if (jpen == ipen) {
        if (kplt == 3) { fprintf(fp, " S\n"); kplt = 2; }
        if (kplt != 2) {
            if (kplt == 1) fprintf(stderr, "penatr: illegal call sequence\n");
            else           fprintf(stderr, "penatr: psplot not open\n");
            return(-1);
        }
        kpatr = 0; lpc = icol; lpt = ityp; penw = thick;
    }
    return(0);
}

int newpen(int npen)
{   if ((npen<1) || (npen>MAXPEN))
        { fprintf(stderr, "newpen: illegal pen-number\n"); return(1); }
    ipen = npen;
    if (kplt == 3) { fprintf(fp, " S\n"); kplt = 2; }
    if (kplt != 2) {
        if (kplt == 1) fprintf(stderr, "newpen: illegal call sequence\n");
        else           fprintf(stderr, "newpen: psplot not open\n");
        return(-1);
    }
    kpatr = 0; lpc = lpenc[ipen-1]; lpt = lpent[ipen-1]; penw = wpen[ipen-1];
    return(0);
}

int qpatr(int mpen, int *lcol, int *ltyp, double *wthick)
{   if ((kplt<2) || (kplt>3)) {
        if (kplt == 1) fprintf(stderr, "qpatr: illegal call sequence\n");
        else           fprintf(stderr, "qpatr: psplot not open\n");
        return(-1);
    }
    if (kplt == 3) { fprintf(fp, " S\n"); kplt = 2; }
    if ((mpen<1) || (mpen>MAXPEN))
        { *lcol = lpc; *ltyp = lpt; *wthick = penw; }
    else {
        *lcol = lpenc[mpen-1]; *ltyp = lpent[mpen-1];
        *wthick = wpen[mpen-1];
    }
    return(ipen);
}

int qpsfp(FILE **fps)
{   if ((kplt<1) || (kplt>3))
        { fprintf(stderr, "qpsfp: psplot not open\n"); return(-1); }
    *fps = fp; return(0);
}

int qapos(int *ixp, int *iyp)
{   if ((kplt<2) || (kplt>3)) {
        if (kplt == 1) fprintf(stderr, "qapos: illegal call sequence\n");
        else fprintf(stderr, "qapos: psplot not open\n");
        return(-1);
    }
    if (kplt == 3) { fprintf(fp, " S\n"); kplt = 2; }
    *ixp = (int)xr; *iyp = (int)yr; return(0);
}

int where(double *wx, double *wy, double *wfac)
{   if ((kplt<2) || (kplt>3)) {
        if (kplt == 1) fprintf(stderr, "where: illegal call sequence\n");
        else           fprintf(stderr, "where: psplot not open\n");
        return(-1);
    }
    *wx = px; *wy = py; *wfac = fs / DPCM; return(0);
}

/*********************/
/*   wrect / wcirc   */
/*********************/

int wrect(double xl, double yl, double xsz, double ysz)
{   int ixl, iyl, ixh, iyh;
    if (kplt == 3) { fprintf(fp, " S\n"); kplt = 2; }
    if (kplt != 2) {
        if (kplt == 1) fprintf(stderr, "wrect: illegal call sequence\n");
        else           fprintf(stderr, "wrect: psplot not open\n");
        return(-1);
    }
    if (kpatr == 0) setpatr();
    ixl = (int)(xorg + xl*fs); iyl = (int)(yorg + yl*fs);
    ixh = ixl + (int)(xsz*fs); iyh = iyl + (int)(ysz*fs);
    fprintf(fp, "N %d %d M %d %d L ", ixl, iyl, ixh, iyl);
    fprintf(fp, "%d %d L %d %d L C 0 J S 1 J\n", ixh, iyh, ixl, iyh);
    return(0);
}

int wcirc(double xc, double yc, double radius, double ang1, double ang2)
{   int ixc, iyc, ira;
    int ia1, ia2;
    if (kplt == 3) { fprintf(fp, " S\n"); kplt = 2; }
    if (kplt != 2) {
        if (kplt == 1) fprintf(stderr, "wcirc: illegal call sequence\n");
        else           fprintf(stderr, "wcirc: psplot not open\n");
        return(-1);
    }
    if (kpatr == 0) setpatr();
    ixc = (int)(xorg + xc*fs); iyc = (int)(yorg + yc*fs);
    ira = (int)(radius * fs); ia1 = (int)ang1; ia2 = (int)ang2;
    if (ang1 == (double)ia1) {
        if (ang2 == (double)ia2)
             fprintf(fp, "N %d %d %d %d %d A S\n", ixc, iyc, ira, ia1, ia2);
        else fprintf(fp, "N %d %d %d %d %f A S\n", ixc, iyc, ira, ia1, ang2);
    }
    else {
        if (ang2 == (double)ia2)
             fprintf(fp, "N %d %d %d %f %d A S\n", ixc, iyc, ira, ang1, ia2);
        else fprintf(fp, "N %d %d %d %f %f A S\n", ixc, iyc, ira, ang1, ang2);
    }
    return(0);
}

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

#define W       double

/******************************************************/
/*   psopn_ / pscls_ / plots_ / plote_ / plot_ /      */
/*   scisor_ / factor_ / where_ / penatr_ / newpen_   */
/******************************************************/

void epsbox_(float *obbx, float *obby, float *sbbx, float *sbby)
{ if (epsbox((W)*obbx,(W)*obby,(W)*sbbx,(W)*sbby) < 0) exit(1); }

void psopn_(char *fnam, char *psize, size_t lf, size_t lp)
{   char pfile[256], sheet[5];
    while ((lf>0) && (*(fnam+lf-1)==' ')) lf--;
    if (lf > 255) lf=255; pfile[0]='\0'; if (lf > 0) strncat(pfile,fnam,lf);
    if (lp > 4) lp=4; sheet[0]='\0'; if (lp > 0) strncat(sheet,psize,lp);
    if (psopn(pfile,sheet) < 0) exit(1);
}

void pscls_(void)
{ if (pscls() < 0) exit(1); }

void plots_(float *xo, float *yo)
{ if (plots((W)*xo,(W)*yo) < 0) exit(1); }

void plote_(void)
{ if (plote() < 0) exit(1); }

void plot_(float *xp, float *yp, int *md)
{ if (plot((W)*xp,(W)*yp,*md) < 0) exit(1); }

void scisor_(float *xlm, float *ylm, float *wlm, float *hlm)
{ if (scisor((W)*xlm,(W)*ylm,(W)*wlm,(W)*hlm) < 0) exit(1); }

void factor_(float *fac)
{ if (factor((W)*fac) < 0) exit(1); }

void penatr_(int *jpen, int *icol, int *ityp, float *thick)
{ if (penatr(*jpen,*icol,*ityp,(W)*thick) < 0) exit(1); }

void newpen_(int *npen)
{ if (newpen(*npen) < 0) exit(1); }

void where_(float *wx, float *wy, float *wfac)
{   double wwx, wwy, wwfac;
    if (where(&wwx,&wwy,&wwfac) < 0) exit(1);
    *wx = wwx; *wy = wwy; *wfac = wwfac;
}

/*********************************/
/*   wrect_ / wpolyg_ / wcirc_   */
/*********************************/

void wrect_(float *xl, float *yl, float *xsz, float *ysz)
{ if (wrect((W)*xl,(W)*yl,(W)*xsz,(W)*ysz) < 0) exit(1); }

void wpolyg_(float *xm, float *ym, int *npt, int *inc)
{   int n;
    if ((*npt<3) || (*inc<1))
        fprintf(stderr, "wpolyg: invalid parameter\n");
    else {
        plot((W)*xm,(W)*ym,3);
        for (n=1; n<*npt; n++)
            { xm += *inc; ym += *inc; plot((W)*xm,(W)*ym,2); }
        clspath(0);
    }
}

void wcirc_(float *xc, float *yc, float *radius, float *ang1, float *ang2)
{ if (wcirc((W)*xc,(W)*yc,(W)*radius,(W)*ang1,(W)*ang2) < 0) exit(1); }