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

#define MAXLVL  50
#define NLVLDF  20
#define NRPS    8       /* number of rasters per sample (both axis)  */
                        /* raster size is defined by DPI in psplot.c */
#define RAD     (180./3.14159)
#define S60     (1.7320508/2.)
#define C60     (1./2.)

int plot(double xp, double yp, int md);
int penatr(int jpen, int icol, int ityp, double thick);
int newpen(int npen);
int qpatr(int mpen, int *lcol, int *ltyp, double *wthick);
int qpsfp(FILE **fps);
int qapos(int *ixp, int *iyp);
int where(double *wx, double *wy, double *wfac);

/************************************************************/
/*   <<<  PSPAINT  >>>   dftone / dfrgbt / dfcols           */
/*                       dframe / paintm / paintw           */
/*                       dfpcol / paintc / paintr / paintp  */
/************************************************************/

static int tgry[NLVLDF]
        = {  28, 36, 48, 60, 72, 84, 96,108,120,132,
            144,156,168,180,192,204,216,228,240,248 };
static int tone[MAXLVL+3] = { 0,
             28, 36, 48, 60, 72, 84, 96,108,120,132,
            144,156,168,180,192,204,216,228,240,248, 252,255 };
static int tcol[MAXLVL+3] = { 0x005576,
            0x00676a, 0x007b5e, 0x009351, 0x00b33f, 0x00e422,
            0x1aff07, 0x46fe00, 0x74ff01, 0xa3ff0d, 0xd0ff16,
            0xffff00, 0xffd016, 0xffa30d, 0xff7401, 0xfe4600,
            0xff1a07, 0xe40022, 0xb3003f, 0x930051, 0x7b005e,
            0x67006a, 0xffffff };
static int tcw[MAXLVL+2][20] = {
          { 0x024059, 0x024a66, 0x025474, 0x025e82, 0x03678f,
            0x03719d, 0x037bab, 0x0485b8, 0x048fc6, 0x0499d4,
            0x0ca1dc, 0x1ca7de, 0x2cade0, 0x3bb3e3, 0x4bb9e5,
            0x5abfe7, 0x6ac5e9, 0x7acbeb, 0x89d1ee, 0x99d7f0 },
          { 0x004c4e, 0x00585a, 0x006466, 0x007072, 0x007c7e,
            0x00878a, 0x009397, 0x009fa3, 0x00abaf, 0x00b7bb,
            0x07bfc3, 0x17c3c7, 0x27c7cb, 0x37cbce, 0x47cfd2,
            0x57d3d6, 0x67d7da, 0x77dbde, 0x87e0e2, 0x97e4e5 },
          { 0x015742, 0x01644d, 0x017157, 0x017f61, 0x028c6c,
            0x029a76, 0x02a780, 0x02b58a, 0x02c295, 0x03cf9f,
            0x0bd7a7, 0x1adaad, 0x2adcb2, 0x3adfb8, 0x4ae1be,
            0x59e4c3, 0x69e6c9, 0x79e9cf, 0x89ecd4, 0x98eeda },
          { 0x055f36, 0x066d3f, 0x077c47, 0x078a50, 0x089958,
            0x09a861, 0x0ab669, 0x0bc572, 0x0cd37a, 0x0de282,
            0x15ea8a, 0x24eb92, 0x33ed99, 0x42eea1, 0x51efa8,
            0x60f1b0, 0x6ff2b7, 0x7ef3bf, 0x8df5c6, 0x9cf6ce },
          { 0x0b642b, 0x0d7331, 0x0f8338, 0x11923f, 0x13a245,
            0x14b14c, 0x16c053, 0x18d059, 0x1adf60, 0x1cef67,
            0x24f76f, 0x32f778, 0x40f881, 0x4ef88a, 0x5cf994,
            0x6af99d, 0x78faa6, 0x87fab0, 0x95fbb9, 0xa3fbc2 },
          { 0x136720, 0x167625, 0x19862a, 0x1d962f, 0x20a634,
            0x23b639, 0x26c63e, 0x29d643, 0x2ce548, 0x2ff54d,
            0x37fd55, 0x44fd60, 0x51fd6b, 0x5dfe76, 0x6afe81,
            0x77fe8c, 0x84fe97, 0x91fea1, 0x9efeac, 0xabfeb7 },
          { 0x1d6716, 0x21771a, 0x25871d, 0x2a9721, 0x2ea724,
            0x33b728, 0x37c62b, 0x3cd62f, 0x40e633, 0x45f636,
            0x4dfe3e, 0x58fe4a, 0x64fe57, 0x6ffe63, 0x7bfe70,
            0x86fe7c, 0x92fe89, 0x9dfe95, 0xa8fea1, 0xb4feae },
          { 0x26650e, 0x2c7511, 0x328413, 0x389415, 0x3ea417,
            0x44b31a, 0x4ac31c, 0x50d31e, 0x56e220, 0x5cf223,
            0x64fa2b, 0x6efa38, 0x78fb46, 0x82fb54, 0x8cfb61,
            0x96fb6f, 0xa0fc7d, 0xaafc8a, 0xb4fc98, 0xbefda6 },
          { 0x316108, 0x387109, 0x40800a, 0x478f0c, 0x4f9e0d,
            0x56ad0e, 0x5ebc10, 0x65cb11, 0x6dda12, 0x74e913,
            0x7cf11b, 0x85f22a, 0x8df339, 0x96f447, 0x9ef556,
            0xa6f565, 0xaff673, 0xb7f782, 0xc0f891, 0xc8f99f },
          { 0x3b5c03, 0x446a04, 0x4d7904, 0x568705, 0x5f9506,
            0x68a306, 0x71b207, 0x7ac007, 0x83ce08, 0x8cdc08,
            0x94e410, 0x9be620, 0xa2e82f, 0xa9e93e, 0xb0eb4e,
            0xb7ed5d, 0xbdee6d, 0xc4f07c, 0xcbf28b, 0xd2f39b },
          { 0x445500, 0x4f6201, 0x597001, 0x647d01, 0x6f8a01,
            0x799701, 0x84a401, 0x8eb201, 0x99bf02, 0xa3cc02,
            0xabd40a, 0xb1d719, 0xb6d929, 0xbbdc39, 0xc1df49,
            0xc6e259, 0xcce468, 0xd1e778, 0xd6ea88, 0xdced98 },
          { 0x4d4d00, 0x595900, 0x656500, 0x717100, 0x7d7d00,
            0x898900, 0x959500, 0xa1a100, 0xadad00, 0xb9b900,
            0xc1c107, 0xc5c517, 0xc9c927, 0xcdcd37, 0xd1d147,
            0xd5d557, 0xd9d967, 0xdddd77, 0xe1e187, 0xe5e597 },
          { 0x554400, 0x624f01, 0x705901, 0x7d6401, 0x8a6f01,
            0x977901, 0xa48401, 0xb28e01, 0xbf9902, 0xcca302,
            0xd4ab0a, 0xd7b119, 0xd9b629, 0xdcbb39, 0xdfc149,
            0xe2c659, 0xe4cc68, 0xe7d178, 0xead688, 0xeddc98 },
          { 0x5c3b03, 0x6a4404, 0x794d04, 0x875605, 0x955f06,
            0xa36806, 0xb27107, 0xc07a07, 0xce8308, 0xdc8c08,
            0xe49410, 0xe69b20, 0xe8a22f, 0xe9a93e, 0xebb04e,
            0xedb75d, 0xeebd6d, 0xf0c47c, 0xf2cb8b, 0xf3d29b },
          { 0x613108, 0x713809, 0x80400a, 0x8f470c, 0x9e4f0d,
            0xad560e, 0xbc5e10, 0xcb6511, 0xda6d12, 0xe97413,
            0xf17c1b, 0xf2852a, 0xf38d39, 0xf49647, 0xf59e56,
            0xf5a665, 0xf6af73, 0xf7b782, 0xf8c091, 0xf9c89f },
          { 0x65260e, 0x752c11, 0x843213, 0x943815, 0xa43e17,
            0xb3441a, 0xc34a1c, 0xd3501e, 0xe25620, 0xf25c23,
            0xfa642b, 0xfa6e38, 0xfb7846, 0xfb8254, 0xfb8c61,
            0xfb966f, 0xfca07d, 0xfcaa8a, 0xfcb498, 0xfdbea6 },
          { 0x671d16, 0x77211a, 0x87251d, 0x972a21, 0xa72e24,
            0xb73328, 0xc6372b, 0xd63c2f, 0xe64033, 0xf64536,
            0xfe4d3e, 0xfe584a, 0xfe6457, 0xfe6f63, 0xfe7b70,
            0xfe867c, 0xfe9289, 0xfe9d95, 0xfea8a1, 0xfeb4ae },
          { 0x671320, 0x761625, 0x86192a, 0x961d2f, 0xa62034,
            0xb62339, 0xc6263e, 0xd62943, 0xe52c48, 0xf52f4d,
            0xfd3755, 0xfd4460, 0xfd516b, 0xfe5d76, 0xfe6a81,
            0xfe778c, 0xfe8497, 0xfe91a1, 0xfe9eac, 0xfeabb7 },
          { 0x640b2b, 0x730d31, 0x830f38, 0x92113f, 0xa21345,
            0xb1144c, 0xc01653, 0xd01859, 0xdf1a60, 0xef1c67,
            0xf7246f, 0xf73278, 0xf84081, 0xf84e8a, 0xf95c94,
            0xf96a9d, 0xfa78a6, 0xfa87b0, 0xfb95b9, 0xfba3c2 },
          { 0x5f0536, 0x6d063f, 0x7c0747, 0x8a0750, 0x990858,
            0xa80961, 0xb60a69, 0xc50b72, 0xd30c7a, 0xe20d82,
            0xea158a, 0xeb2492, 0xed3399, 0xee42a1, 0xef51a8,
            0xf160b0, 0xf26fb7, 0xf37ebf, 0xf58dc6, 0xf69cce },
          { 0x570142, 0x64014d, 0x710157, 0x7f0161, 0x8c026c,
            0x9a0276, 0xa70280, 0xb5028a, 0xc20295, 0xcf039f,
            0xd70ba7, 0xda1aad, 0xdc2ab2, 0xdf3ab8, 0xe14abe,
            0xe459c3, 0xe669c9, 0xe979cf, 0xec89d4, 0xee98da },
          { 0x4c004e, 0x58005a, 0x640066, 0x700072, 0x7c007e,
            0x87008a, 0x930097, 0x9f00a3, 0xab00af, 0xb700bb,
            0xbf07c3, 0xc317c7, 0xc727cb, 0xcb37ce, 0xcf47d2,
            0xd357d6, 0xd767da, 0xdb77de, 0xe087e2, 0xe497e5 } } ;
static int nlvl=NLVLDF;
static int pcol=0;
static int kcol=1, kpcol=0, nrps=NRPS, kdir=0, imx, jmx, nx, ny;
static double t1, t2, t3, t4, t5, t6, dx, dy;
static FILE *fp;

int dftone(int nlevel, int mtone[], int ktnlo, int ktnhi, int ktnvd)
{   int i, j, k=0;
    if (nlevel > MAXLVL)
        { fprintf(stderr, "dftone: too many levels\n"); return(-1); }
    kcol = 0;
    if (nlevel < 0) { nlvl = 256; return(0); }
    if (nlevel == 0) {
        nlvl = NLVLDF;
        tone[0] = 0;
        for (i=1; i<=nlvl; i++) tone[i] = tgry[i-1];
        tone[nlvl+1] = 252;
        tone[nlvl+2] = 255;
    }
    else {
        nlvl = nlevel;
        for (i=1; i<=nlvl; i++) {
            if (((j=mtone[i-1])>=0) && (j<=255)) tone[i] = j;
            else { tone[i] = 0; k++; }
        }
        if ((ktnlo>=0) && (ktnlo<=255)) tone[0] = ktnlo;
        else { tone[0] = 0; k++; }
        if ((ktnhi>=0) && (ktnhi<=255)) tone[nlvl+1] = ktnhi;
        else { tone[nlvl+1] = 0; k++; }
        if ((ktnvd>=0) && (ktnvd<=255)) tone[nlvl+2] = ktnvd;
        else tone[nlvl+2] = 255;
    }
    if (k != 0) fprintf(stderr, "dftone: out-of-range tone forced to black\n");
    return(k);
}

int dfcols(int nclev, double fh, double fd, int kcnv)
{   int k, i, j, r, g, b;
    double d, x, a, t, cosa, sina, cr, cg, cb, c1, c2, ss;

    if (nclev == 0) { nclev = NLVLDF; fh = 1., fd = 1., kcnv = -1; }
    if (nclev < 0)
        { fprintf(stderr, "dfcols: no level specified\n"); return(-1); }
    if (nclev > MAXLVL)
        { fprintf(stderr, "dfcols: too many levels\n"); return(-1); }
    nlvl = nclev; kcol = 1; k = 0;
    if ((fh>1.2) || (fh<0.2)) {
        fprintf(stderr, "dfcols: out-of-range hue factor forced to 1.\n");
        fh = 1.; k++;
    }
    if ((fd>1.6) || (fd<0.)) {
        fprintf(stderr, "dfcols: out-of-range darkness factor forced to 1.\n");
        fd = 1.; k++;
    }
    d = 1. / (double)(nlvl/2 + 1);
    for (i=0; i<(nlvl+2); i++) {
        x = (d*i - 1.) * fh; a = (16.*x*x + 120.) * x;
        cosa = cos(a/RAD); sina = sin(a/RAD);
        cr = cosa*C60 + sina*S60; cg = cosa*C60 - sina*S60; cb = -cosa;
        t = 1. - 1./(9.*x*x + 1.); t = 170. - (340./3.)*fd*t;
        c1 = cr; if (cg>c1) c1 = cg; if (cb>c1) c1 = cb;
        c2 = cr; if (cg<c2) c2 = cg; if (cb<c2) c2 = cb;
        ss = (255.-t) / c1; if (ss > -t/c2) ss = -t/c2;
        r = t + cr*ss; g = t + cg*ss; b = t + cb*ss;
        tcol[i] = r*0x010000 + g*0x0100 + b;
        for (j=0; j<20; j++) {
            t = ((double)(j+j-19)/64.*fd + 0.5) * 255.;
            if (t < 127.5) ss = t; else ss = 255. - t;
            r = t + cr*ss; g = t + cg*ss; b = t + cb*ss;
            tcw[i][j] = r*0x010000 + g*0x0100 + b;
        }
    }
    if ((kcnv>=0) && (kcnv<=0xffffff)) tcol[nlvl+2] = kcnv;
    else    tcol[nlvl+2] = 0xffffff;
    return(k);
}

int dfrgbt(int nclev, int mrgb[], int kclo, int kchi, int kcnv)
{   int i, k=0;
    int kc;
    if (nclev < 0)
        { fprintf(stderr, "dfrgbt: no level specified\n"); return(-1); }
    if (nclev > MAXLVL)
        { fprintf(stderr, "dfrgbt: too many levels\n"); return(-1); }
    kcol = 1;
    if (nclev == 0) dfcols(NLVLDF, 1., 1., -1);
    else {
        nlvl = nclev;
        for (i=1; i<=nlvl; i++) {
            if (((kc=mrgb[i-1])>=0) && (kc<=0xffffff)) tcol[i] = kc;
            else { tcol[i] = 0; k++; }
        }
        if ((kclo>=0) && (kclo<=0xffffff)) tcol[0] = kclo;
        else { tcol[0] = 0; k++; }
        if ((kchi>=0) && (kchi<=0xffffff)) tcol[nlvl+1] = kchi;
        else { tcol[nlvl+1] = 0; k++; }
        if ((kcnv>=0) && (kcnv<=0xffffff)) tcol[nlvl+2] = kcnv;
        else tcol[nlvl+2] = 0xffffff;
    }
    if (k != 0)
        fprintf(stderr, "dfrgbt: out-of-range color forced to black\n");
    return(k);
}

int dfc40s(void)
{   int kclo=0x1e1eff, kchi=0x6918ff, kcnv=0xffffff;
    int mrgb[40]={ 0x4141ff, 0x3b58ff, 0x2f6aff, 0x217bff, 0x108fff, 0x00a5fc,
         0x00b8ec, 0x00cbdc, 0x00decc, 0x00f1bc, 0x08ffa9, 0x25ff92, 0x3dff80,
         0x51ff72, 0x65ff64, 0x78ff5a, 0x8cff4f, 0xa2ff45, 0xbaff39, 0xd7ff25,
         0xffff00, 0xffe11b, 0xffc92a, 0xffb532, 0xffa239, 0xff9042, 0xff7e4b,
         0xff6e55, 0xff5f5e, 0xff4e69, 0xff3d74, 0xff2c82, 0xff1893, 0xff01a7,
         0xec00b8, 0xd900c8, 0xc700d7, 0xb400e7, 0xa100f7, 0x860dff };
    return(dfrgbt(40,mrgb,kclo,kchi,kcnv));
}

int dfpcol(int key, int mcol)
{   int k=0;
    if ((key==2) && (nlvl==256)) key = 0;
    switch (key) {
    case 0:
        if ((mcol<0) || (mcol>255)) {
            fprintf(stderr, "dfpcol: out-of-range color forced to black\n");
            kpcol = 0; pcol = 0; k = 1;
        }
        else { kpcol = 0; pcol = mcol; }
        break;
    case 1:
        if ((mcol<-255) || (mcol>0xffffff)) {
            fprintf(stderr, "dfpcol: out-of-range color forced to black\n");
            kpcol = 0; pcol = 0; k = 1;
        }
        else if (mcol < 0) { kpcol = 0; pcol = -mcol; }
        else { kpcol = 1; pcol = mcol; }
        break;
    case 2:
        if ((mcol<0) || (mcol>(nlvl+2))) {
            fprintf(stderr, "dfpcol: out-of-range tone assumed black\n");
            kpcol = 0; pcol = 0; k = 1;
        }
        else {
            if (kcol) { kpcol = 1; pcol = tcol[mcol]; }
            else { kpcol = 0; pcol = tone[mcol]; }
        }
        break;
    default:
        fprintf(stderr, "dfpcol: illegal key value\n"); k = -1;
    }
    return(k);
}

int dresol(int nras)
{   if (nras < 1) { nrps = NRPS; return(1); }
    else { nrps = nras; return(0); }
}

int dframe(double xo, double yo, double wd, double hi, int imax, int jmax)
{   int ixl, iyl, ixh, iyh, iwd, ihi;
    double wx, wy, wf;
    if ((imax>1) && (jmax>1)) { kdir = 1; imx = imax-1; jmx = jmax-1; }
    else if ((imax<-1) && (jmax<-1))
        { kdir = 2; imx = -jmax-1; jmx = -imax-1; }
    else { fprintf(stderr, "dframe: array size rejected\n"); return(-1); }
    if (where(&wx,&wy,&wf) < 0) {
        fprintf(stderr, "dframe: error at calling where\n");
        kdir = 0; return(-1);
    }
    plot(xo,yo,3); qapos(&ixl,&iyl);
    plot(xo+wd,yo+hi,3); qapos(&ixh,&iyh); plot(wx,wy,3);
    iwd = abs(ixh-ixl); ihi = abs(iyh-iyl);
    if (iwd < nrps)
        { fprintf(stderr, "dframe: too narrow\n"); kdir = 0; return(-1); }
    if (ihi < nrps)
        { fprintf(stderr, "dframe: too short\n"); kdir = 0; return(-1); }
    if (kdir == 1) {
        nx = (iwd-1)/nrps + 1; ny = (ihi-1)/nrps + 1;
        t1 = (double)nx/(ixh-ixl); t4 = (double)ny/(iyh-iyl);
        t2 = 0.; t3 = 0.; t5 = -t1*ixl; t6 = -t4*iyl;
    }
    else {
        nx = (ihi-1)/nrps + 1; ny = (iwd-1)/nrps + 1;
        t2 = (double)ny/(ixh-ixl); t3 = (double)nx/(iyh-iyl);
        t1 = 0.; t4 = 0.; t5 = -t3*iyl; t6 = -t2*ixl;
    }
    dx = (double)imx/nx; dy = (double)jmx/ny;
    return(0); 
}

int dframo(double x0, double y0, double x1, double y1,
           double x2, double y2, int imax, int jmax)
{   int ix0, iy0, ix1, iy1, ix2, iy2;
    double wx, wy, wf, dxi, dyi, dxj, dyj, wd, hi, dds;
    kdir = 0; imx = imax-1; jmx = jmax-1;
    if ((imax<2) || (jmax<2))
        { fprintf(stderr, "dframo: array size rejected\n"); return(-1); }
    if (where(&wx,&wy,&wf) < 0)
        { fprintf(stderr, "dframo: error at calling where\n"); return(-1); }
    plot(x0,y0,3); qapos(&ix0,&iy0); plot(x1,y1,3); qapos(&ix1,&iy1);
    plot(x2,y2,3); qapos(&ix2,&iy2); plot(wx,wy,3);
    dxi = (double)(ix1-ix0); dyi = (double)(iy1-iy0);
    dxj = (double)(ix2-ix0); dyj = (double)(iy2-iy0);
    wd = sqrt(dxi*dxi + dyi*dyi); hi = sqrt(dxj*dxj + dyj*dyj);
    nx = (int)wd/nrps + 1; ny = (int)hi/nrps + 1;
    if ((nx<2) || (ny<2))
        { fprintf(stderr, "dframo: too small picture\n"); return(-1); }
    dds = dxi*dyj - dxj*dyi;
    if (dds == 0.)
        { fprintf(stderr, "dframo: invalid coord.\n"); return(-1); }
    t1 =  dyj*nx/dds; t2 = -dyi*ny/dds; t3 = -dxj*nx/dds; t4 =  dxi*ny/dds;
    t5 = - (ix0*t1 + iy0*t3); t6 = - (ix0*t2 + iy0*t4); kdir = 3;
    dx = (double)imx/nx; dy = (double)jmx/ny;
    return(0); 
}

int paintm(float *v, float vlo, float vhi, float vnvd)
{   double wt, xb, yb, x, y;
    double vstp, vmin, va, vb, vc;
    float *pv, *pv1, *pv2, *pv3, *pv4;
    int pen, col, typ;
    int n, r, g, b, ix, iy, l, ls, ct;
    if (kdir == 0)
        { fprintf(stderr, "paintm: frame undefined\n"); return(-1); }
    if (vhi <= vlo)
        { fprintf(stderr, "paintm: invalid parameter\n"); return(-1); }
    vstp = (vhi-vlo)/nlvl; vmin = vlo-vstp;
    if ((pen=qpatr(0,&col,&typ,&wt)) < 0)
        { fprintf(stderr, "paintm: error at calling qpatr\n"); return(-1); }
    qpsfp(&fp);
    if (kcol) {
        fprintf(fp, "CmapPacket\n");
        fprintf(fp, " %02x", nlvl+2);
        for (n=0; n<(nlvl+3); n++) {
            if ((n%10) == 0) fputc('\n',fp);
            r = tcol[n]/0x010000; g = (tcol[n]/256)%256; b = tcol[n]%256;
            fprintf(fp, " %02x%02x%02x", r,g,b);
        }
        fputc('\n',fp);
    }
    fprintf(fp, "%d %d 8 ", nx,ny);
    fprintf(fp, "[%.5f %.5f %.5f %.5f %.2f %.2f]", t1,t2,t3,t4,t5,t6);
    if (kcol) fprintf(fp,"\n {CcodePacket} false 3 colorimage\n ");
    else fprintf(fp," {GrayPacket} image\n ");
    n = 0; ls = -1; pv2 = v; yb = 0.;
    for (y=dy/2,iy=0; iy<ny; iy++,y+=dy) {
        while (y > yb) { pv1 = pv2; pv2+=(imx+1); yb+=1.; }
        if ((*pv1==vnvd) || (*pv2==vnvd)) vb = vnvd;
        else vb = *pv1*(yb-y) + *pv2*(y+1.-yb);
        pv3 = pv1; pv4 = pv2; xb = 0.;
        for (x=dx/2,ix=0; ix<nx; ix++,x+=dx) {
            while (x > xb) {
                pv3++; pv4++; xb+=1.; va = vb;
                if ((*pv3==vnvd) || (*pv4==vnvd)) vb = vnvd;
                else vb = *pv3*(yb-y) + *pv4*(y+1.-yb);
            }
            if ((va==vnvd) || (vb==vnvd)) l = nlvl+2;
            else {
                vc = va*(xb-x) + vb*(x+1.-xb);
                l = (int)((vc-vmin)/vstp);
                if (l < 0) l = 0;
                else if (l > nlvl) l = nlvl+1;
            }
            if (ls == -1) { ls = l; ct = 0; }
            else if ((ls==l) && (ct<255)) ct++;
            else {
                if (!kcol) {
                    if (nlvl == 256) {
                        if (ls != 0) ls--;
                        if (ls > 255) ls = 255;
                    }
                    else ls = tone[ls];
                }
                fprintf(fp, "%02x%02x", ct,ls); n++;
                ls = l; ct = 0;
            }
            if (n == 18) { fprintf(fp,"\n "); n = 0; }
        }
        if (!kcol) {
            if (nlvl == 256) {
                if (ls != 0) ls--;
                if (ls > 255) ls = 255;
            }
            else ls = tone[ls];
        }
        fprintf(fp, "%02x%02x\n ", ct,ls); n = 0; ls = -1;
    }
    fputc('\n',fp);
    if (pen == 0) penatr(0, col, typ, wt);
    else newpen(pen);
    return(0);
}

int paintw(float *v, float vlo, float vhi, float vnvd, float *u)
{   double wt, xb, yb, x, y;
    double vstp, vmin, va, vb, vc, ua, ub, uc;
    float *pv, *pv1, *pv2, *pv3, *pv4;
    float *pu, *pu1, *pu2, *pu3, *pu4;
    int pen, col, typ, wcol;
    int n, r, g, b, ix, iy, l, ls, m, ms, ct;
    if (kdir == 0)
        { fprintf(stderr, "paintw: frame undefined\n"); return(-1); }
    if (vhi <= vlo)
        { fprintf(stderr, "paintw: invalid parameter\n"); return(-1); }
    vstp = (vhi-vlo)/nlvl; vmin = vlo-vstp;
    if ((pen=qpatr(0,&col,&typ,&wt)) < 0)
        { fprintf(stderr, "paintw: error at calling qpatr\n"); return(-1); }
    qpsfp(&fp);
    fprintf(fp, "%d %d 8 ", nx,ny);
    fprintf(fp, "[%.5f %.5f %.5f %.5f %.2f %.2f]\n", t1,t2,t3,t4,t5,t6);
    fprintf(fp," {ColorPacket} false 3 colorimage\n ");
    n = 0; ls = -1; pv2 = v; pu2 = u; yb = 0.;
    for (y=dy/2,iy=0; iy<ny; iy++,y+=dy) {
        while (y > yb)
            { pv1 = pv2; pv2+=(imx+1); pu1 = pu2; pu2+=(imx+1); yb+=1.; }
        if ((*pv1==vnvd) || (*pv2==vnvd)) vb = vnvd;
        else { vb = *pv1*(yb-y) + *pv2*(y+1.-yb);
               ub = *pu1*(yb-y) + *pu2*(y+1.-yb); }
        pv3 = pv1; pv4 = pv2; pu3 = pu1; pu4 = pu2; xb = 0.;
        for (x=dx/2,ix=0; ix<nx; ix++,x+=dx) {
            while (x > xb) {
                pv3++; pv4++; pu3++; pu4++; xb+=1.; va = vb; ua = ub;
                if ((*pv3==vnvd) || (*pv4==vnvd)) vb = vnvd;
                else { vb = *pv3*(yb-y) + *pv4*(y+1.-yb);
                       ub = *pu3*(yb-y) + *pu4*(y+1.-yb); }
            }
            if ((va==vnvd) || (vb==vnvd)) l = nlvl+2;
            else {
                vc = va*(xb-x) + vb*(x+1.-xb);
                uc = ua*(xb-x) + ub*(x+1.-xb);
                l = (int)((vc-vmin)/vstp); m = (int)(uc*20.);
                if (l < 0) l = 0; else if (l > nlvl) l = nlvl+1;
                if (m < 0) m = 0; else if (m > 19) m = 19;
            }
            if (ls == -1) { ls = l; ms = m; ct = 0; }
            else if ((ls==l) && (ms==m) && (ct<255)) ct++;
            else {
                if (ls == (nlvl+2)) wcol = tcol[ls]; else wcol = tcw[ls][ms];
                fprintf(fp, "%02x%06x", ct,wcol); n++;
                ls = l; ms = m; ct = 0;
            }
            if (n == 9) { fprintf(fp,"\n "); n = 0; }
        }
        if (ls == (nlvl+2)) wcol = tcol[ls]; else wcol = tcw[ls][ms];
        fprintf(fp, "%02x%06x\n ", ct,wcol); n = 0; ls = -1;
    }
    fputc('\n',fp);
    if (pen == 0) penatr(0, col, typ, wt);
    else newpen(pen);
    return(0);
}

int paintc(double xc, double yc, double radius)
{   int pen, col, typ, ix0, iy0, ird, r, g, b;
    double wx, wy, wt;
    if (where(&wx,&wy,&wt) < 0)
        { fprintf(stderr, "paintc: error at calling where\n"); return(-1); }
    pen=qpatr(0,&col,&typ,&wt); qpsfp(&fp);
    plot(xc+radius,yc,3); qapos(&ird,&iy0);
    plot(xc,yc,3); qapos(&ix0,&iy0); ird -= ix0;
    if (kpcol) {
        r = pcol/0x010000; g = (pcol/256)%256; b = pcol%256;
        fprintf(fp,"%.3f %.3f %.3f RGB\n",r/255.,g/255.,b/255.);
    }
    else fprintf(fp,"%.3f T\n",pcol/255.);
    fprintf(fp,"N %d %d %d 0 360 A C fill\n",ix0,iy0,ird);
    if (pen == 0) penatr(0, col, typ, wt);
    else newpen(pen);
    return(plot(wx,wy,3));
}

int paintr(double xorg, double yorg, double wide, double high)
{   int pen, col, typ, ix1, iy1, ix2, iy2, r, g, b;
    double wx, wy, wt;
    if (where(&wx,&wy,&wt) < 0)
        { fprintf(stderr, "paintr: error at calling where\n"); return(-1); }
    pen=qpatr(0,&col,&typ,&wt); qpsfp(&fp);
    plot(xorg,yorg,3); qapos(&ix1,&iy1);
    plot(xorg+wide,yorg+high,3); qapos(&ix2,&iy2);
    if (kpcol) {
        r = pcol/0x010000; g = (pcol/256)%256; b = pcol%256;
        fprintf(fp,"%.3f %.3f %.3f RGB\n",r/255.,g/255.,b/255.);
    }
    else fprintf(fp,"%.3f T\n",pcol/255.);
    fprintf(fp,"N %d %d M %d %d L",ix1,iy1,ix2,iy1);
    fprintf(fp," %d %d L %d %d L C fill\n",ix2,iy2,ix1,iy2);
    if (pen == 0) penatr(0, col, typ, wt);
    else newpen(pen);
    return(plot(wx,wy,3));
}

int paintp(float xa[], float ya[], int npt, int inc)
{   int pen, col, typ, i, ix, iy, r, g, b, nseq;
    double wx, wy, wt;
    if ((npt<3) || (inc<1))
        { fprintf(stderr, "paintp: invalid parameter\n"); return(-1); }
    if (where(&wx,&wy,&wt) < 0)
        { fprintf(stderr, "paintp: error at calling where\n"); return(-1); }
    pen=qpatr(0,&col,&typ,&wt); qpsfp(&fp);
    if (kpcol) {
        r = pcol/0x010000; g = (pcol/256)%256; b = pcol%256;
        fprintf(fp, "%.3f %.3f %.3f RGB\n", r/255.,g/255.,b/255.);
    }
    else fprintf(fp, "%.3f T\n", pcol/255.);
    plot((double)xa[0],(double)ya[0],3); qapos(&ix,&iy);
    fprintf(fp, "N %d %d M", ix,iy); nseq = 1;
    for (i=1; i<npt; i++) {
        if (nseq > 5) { fputc('\n',fp); nseq = 0; }
        plot((double)xa[i*inc],(double)ya[i*inc],3); qapos(&ix,&iy);
        fprintf(fp, " %d %d L", ix,iy); nseq++;
    }
    fprintf(fp, " C fill\n");
    if (pen == 0) penatr(0, col, typ, wt);
    else newpen(pen);
    return(plot(wx,wy,3));
}

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

#define W   double

void dftone_(int *nlevel, int mtone[], int *ktnlo, int *ktnhi, int *ktnvd)
{   int rc;
    if (*nlevel < 0) rc = dftone(-1,tgry,0,0,0);
    else if (*nlevel == 0) rc = dftone(0,tgry,0,0,0);
    else rc = dftone(*nlevel,mtone,*ktnlo,*ktnhi,*ktnvd);
    if (rc < 0) exit(1);
}

void dfcols_(int *nclev, float *fh, float *fd, int *kcnv)
{   if (*nclev == 0) dfcols(0,0,0,0);
    else { if (dfcols(*nclev,*fh,*fd,*kcnv) < 0) exit(1); }
}

void dfrgbt_(int *nclev, int mrgb[], int *kclo, int *kchi, int *kcnv)
{   if (*nclev == 0) dfcols(0,0,0,0);
    else { if (dfrgbt(*nclev,mrgb,*kclo,*kchi,*kcnv) < 0) exit(1); }
}

void dfpcol_(int *key, int *mcol)
{ if (dfpcol(*key,*mcol) < 0) exit(1); }

void dresol_(int *nras) { dresol(*nras); }

void dframe_(float *xo, float *yo, float *wd, float *hi, int *imax, int *jmax)
{ if (dframe((W)*xo,(W)*yo,(W)*wd,(W)*hi,*imax,*jmax) < 0) exit(1); }

void dframo_(float *x0, float *y0, float *x1, float *y1,
             float *x2, float *y2, int *imax, int *jmax)
{ if (dframo((W)*x0,(W)*y0,(W)*x1,(W)*y1,
             (W)*x2,(W)*y2,*imax,*jmax) < 0) exit(1); }

void paintm_(float *v, float *vlo, float *vhi, float *vnvd)
{ if (paintm(v,*vlo,*vhi,*vnvd) < 0) exit(1); }

void paintw_(float *v, float *vlo, float *vhi, float *vnvd, float *u)
{ if (paintw(v,*vlo,*vhi,*vnvd,u) < 0) exit(1); }

void paintc_(float *xc, float *yc, float *radius)
{ if (paintc((W)*xc,(W)*yc,(W)*radius) < 0) exit(1); }

void paintr_(float *xorg, float *yorg, float *wide, float *high)
{ if (paintr((W)*xorg,(W)*yorg,(W)*wide,(W)*high) < 0) exit(1); }

void paintp_(float xa[], float ya[], int *npt, int *inc)
{ if (paintp(xa,ya,*npt,*inc) < 0) exit(1); }

void dfc40s_(void) { if (dfc40s() < 0) exit(1); }