Usenet 1994 January/usenetsourcesnewsgroupsinfomagicjanuary1994.iso/sources/misc/volume24/gnuplot3/part13

Text File | 1991-10-27 | 50KB | 1,758 lines

Newsgroups: comp.sources.misc
From: [e-mail] (Elber Gershon)
Subject:  v24i035:  gnuplot3 - interactive function plotting utility, Part13/26
Message-ID: <[e-mail]>
X-Md4-Signature: 93f50cac37b5b1fd89d0424c73f82ec2
Date: Mon, 28 Oct 1991 00:22:21 GMT
Approved: [e-mail]

Submitted-by: [e-mail] (Elber Gershon)
Posting-number: Volume 24, Issue 35
Archive-name: gnuplot3/part13
Environment: UNIX, MS-DOS, VMS
Supersedes: gnuplot2: Volume 11, Issue 65-79

#!/bin/sh
# this is Part.13 (part 13 of a multipart archive)
# do not concatenate these parts, unpack them in order with /bin/sh
# file gnuplot/command.c continued
#
if test ! -r _shar_seq_.tmp; then
    echo 'Please unpack part 1 first!'
    exit 1
fi
(read Scheck
 if test "$Scheck" != 13; then
    echo Please unpack part "$Scheck" next!
    exit 1
 else
    exit 0
 fi
) < _shar_seq_.tmp || exit 1
if test ! -f _shar_wnt_.tmp; then
    echo 'x - still skipping gnuplot/command.c'
else
echo 'x - continuing file gnuplot/command.c'
sed 's/^X//' << 'SHAR_EOF' >> 'gnuplot/command.c' &&
X        /* function ! */
X        start_token = c_token;
X        copy_str(c_dummy_var[0], c_token + 2);
X        if(equals(c_token+3,",")) {
X            copy_str(c_dummy_var[1], c_token + 4);
X            c_token += 2;  /* skip the  , dummy2 */
X        }
X        c_token += 5; /* skip (, dummy, ) and = */
X        if (END_OF_COMMAND)
X            int_error("function definition expected",c_token);
X        udf = dummy_func = add_udf(start_token);
X        if (udf->at)                /* already a dynamic a.t. there */
X            free((char *)udf->at);    /* so free it first */
X        if ((udf->at = perm_at()) == (struct at_type *)NULL)
X            int_error("not enough memory for function",start_token);
X        m_capture(&(udf->definition),start_token,c_token-1);
X    }
X    else {
X        /* variable ! */
X        start_token = c_token;
X        c_token +=2;
X        udv = add_udv(start_token);
X        (void) const_express(&(udv->udv_value));
X        udv->udv_undef = FALSE;
X    }
}
X
X
get_data(this_plot)
struct curve_points *this_plot;
{
register int i, l_num, datum;
register FILE *fp;
float x, y;
float ylow, yhigh;            /* for error bars */
float temp;
BOOLEAN yfirst;
char format[MAX_LINE_LEN+1], data_file[MAX_LINE_LEN+1], line[MAX_LINE_LEN+1];
char *float_format = "%f", *float_skip = "%*f";
int xcol = 1, ycol = 2, yemin = 3, yemax = 4;
struct value colvalue;
#ifdef AMIGA_LC_5_1
int num_perc_ast;
char *start_search;
#endif /* AMIGA_LC_5_1 */
X
X    quote_str(data_file, c_token);
X    this_plot->plot_type = DATA;
X    if ((fp = fopen(data_file, "r")) == (FILE *)NULL)
X        os_error("can't open data file", c_token);
X
X    format[0] = '\0';
X    yfirst = FALSE;
X    c_token++;    /* skip data file name */
X    if (almost_equals(c_token,"u$sing")) {
X        c_token++;      /* skip "using" */
X        if (!END_OF_COMMAND && !isstring(c_token)) {
X            struct value a;
X            ycol = (int)magnitude(const_express(&a));
X            if (equals(c_token,":")) {
X                c_token++;      /* skip ":" */
X                xcol = ycol;
X                ycol = (int)magnitude(const_express(&a));
X                if (equals(c_token,":")) {
X                    c_token++;      /* skip ":" */
X                    yemin = (int)magnitude(const_express(&a));
X                    if (equals(c_token,":")) {
X                        c_token++;      /* skip ":" */
X                        yemax = (int)magnitude(const_express(&a));
X                    }
X                    else
X                            yemax = -1;
X                }
X                else {
X                        yemin = -1;
X                        yemax = -1;
X                }
X            }
X            else {
X                    yemin = -1;
X                    yemax = -1;
X            }
X            if (xcol > ycol) yfirst = TRUE;
X        }
X        if (!END_OF_COMMAND && isstring(c_token)) {
X            quotel_str(format, c_token);
X            c_token++;    /* skip format */
X        }
X    }
X    if (strlen(format) == 0) {
X        for(i = 1; i <= ((xcol > ycol) ? xcol : ycol); i++)
X            if ((i == xcol) || (i == ycol))
X                (void) strcat(format,float_format);
X            else
X                (void) strcat(format,float_skip);
X
X        if (yemin > 0) {
X            /* We have error bars - handle them. */
X            yemin -= (xcol > ycol) ? xcol : ycol;
X            yemax -= (xcol > ycol) ? xcol : ycol;
X            if (yemin > yemax && yemax > 0) {
X                i = yemin;
X                yemin = yemax;
X                yemax = i;
X            }
X
X            if (yemin == yemax)
X                int_error("Two error bar columns are the same",
X                      NO_CARET);
X
X            if (yemin <= 0)
X                int_error("Error bar columns must follow data columns",
X                      NO_CARET);
X
X            for (i = 1; i < yemin; i++)
X                (void) strcat(format,float_skip);
X            (void) strcat(format,float_format);
X
X            if (yemax > 0) {
X                for (i = 1; i < yemax - yemin; i++)
X                    (void) strcat(format,float_skip);
X                (void) strcat(format,float_format);
X            }
X        }
X    }
X
X    l_num = 0;
X     datum = 0;
X    i = 0;
#ifdef AMIGA_LC_5_1
X    num_perc_ast = 0;
X    start_search = format;
X    while (*start_search != '\0') {
X        if (start_search[0] == '%')
X            if (start_search[1] == '*') num_perc_ast++;
X        start_search++;
X    }
#endif /* AMIGA_LC_5_1 */
X    while ( fgets(line, MAX_LINE_LEN, fp) != (char *)NULL ) {
X        l_num++;
X        if (is_comment(line[0]))
X            continue;        /* ignore comments */
X        if (i >= this_plot->p_max) {
X            /* overflow about to occur. Extend size of points[]
X            * array. We either double the size, or add 1000 points,
X            * whichever is a smaller increment. Note i=p_max.
X            */
X            cp_extend(this_plot, i + (i < 1000 ? i : 1000));
X        }
X        if (!line[1]) { /* is it blank line ? */
X            /* break in data, make next point undefined */
X            this_plot->points[i].type = UNDEFINED;
X            i++;
X            continue;
X        }
X
#ifdef AMIGA_LC_5_1
X        switch (sscanf(line, format, &x, &y, &ylow, &yhigh) -
X             num_perc_ast) {
#else /* AMIGA_LC_5_1 */
X        switch (sscanf(line, format, &x, &y, &ylow, &yhigh)) {
#endif /* AMIGA_LC_5_1 */
X            case 1: {        /* only one number on the line */
X               y = x;        /* assign that number to y */
X               x = datum;    /* and make the index into x */
X               /* no break; !!! */
X            }
X            case 2: {
X               if (yfirst) { /* exchange x and y */
X                  temp = y;
X                  y = x;
X                  x = temp;
X               }
X               datum++;
X
X               /* ylow and yhigh are same as y */
X               store2d_point(this_plot, i++, x, y, y, y);
X               break;
X            }
X            case 3: {        /* x, y, ydelta */
X               if (yfirst) { /* exchange x and y */
X                  temp = y;
X                  y = x;
X                  x = temp;
X               }
X               datum++;
X
X               /* ydelta is in the ylow variable */
X               store2d_point(this_plot, i++, x, y, y-ylow, y+ylow);
X               break;
X            }
X            case 4: {        /* x, y, ylow, yhigh */
X               if (yfirst) { /* exchange x and y */
X                  temp = y;
X                  y = x;
X                  x = temp;
X               }
X               datum++;
X
X               store2d_point(this_plot, i++, x, y, ylow, yhigh);
X               break;
X            }
X            default: {
X               (void) sprintf(line, "bad data on line %d", l_num);
X               int_error(line,c_token);
X            }
X        }
X    }
X    this_plot->p_count = i;
X    cp_extend(this_plot, i);    /* shrink to fit */
X    (void) fclose(fp);
}
X
/* called by get_data for each point */
store2d_point(this_plot, i, x, y, ylow, yhigh)
struct curve_points *this_plot;
int i;                    /* point number */
#ifdef AMIGA_LC_5_1
double x, y;
double ylow, yhigh;
#else
float x, y;
float ylow, yhigh;
#endif
{
X    struct coordinate *cp = &(this_plot->points[i]);
X
X    /* the easy part: */
X    cp->type = INRANGE;
X    cp->x = x;
X    cp->y = y;
X    cp->ylow = ylow;
X    cp->yhigh = yhigh;
X
X    /* Adjust for log scale. */
X    if (log_x)
X     cp->type = adjustlog(cp->type, &(cp->x));
X    if (log_y) {
X       cp->type = adjustlog(cp->type, &(cp->y));
X       /* Note ylow,yhigh can't affect cp->type. */
X       (void)     adjustlog(cp->type, &(cp->ylow));
X       (void)     adjustlog(cp->type, &(cp->yhigh));
X    }
X
X    /* Now adjust the xrange, or declare the point out of range */
X    /* The yrange is handled later, once we know whether to
X    * include ylow, yhigh in the calculation. See adjust_yrange()
X    */
X    if (cp->type == INRANGE)
X     if (autoscale_lx || inrange(x,xmin,xmax)) {
X        if (autoscale_lx) {
X            if (x < xmin) xmin = x;
X            if (x > xmax) xmax = x;
X        }
X     } else {
X        cp->type = OUTRANGE;
X     }
}
X
/* Adjust for log scale:
X * take the log of the second parameter, in place, if possible.
X * If not possible, return new type for point; if possible, then
X * return old type for point.
X */
static enum coord_type
adjustlog(type, val)
X    enum coord_type type;
X    coordval *val;
{
X    if (*val < 0.0) {
X       return(UNDEFINED);
X    } else if (*val == 0.0) {
X       *val = -VERYLARGE;
X       return(OUTRANGE);
X    } else {
X       *val = log10(*val);
X       return(type);
X    }
}
X
X
/* now adjust the yrange, or declare the point out of range */
/* this does all points in a curve */
adjust_yrange(curve)
X    struct curve_points *curve;
{
X    BOOLEAN ebars = (curve->plot_style == ERRORBARS);
X    int npoints = curve->p_count; /* number of points */
X    coordval y, ylow, yhigh;    /* one point value */
X    struct coordinate *cp;    /* one coordinate */
X    int i;                /* index into points */
X
X    for (i = 0; i < npoints; i++) {
X       cp = &(curve->points[i]);
X       if (cp->type == INRANGE) {
X          y = cp->y;
X          if ((autoscale_ly ||
X               inrange(log_y ? pow(10.0, y) : y, ymin, ymax) ||
X               polar)) {
X             if (autoscale_ly) {
X                if (y < ymin) ymin = y;
X                if (y > ymax) ymax = y;
X                if (ebars) {
X                    ylow = cp->ylow;
X                    yhigh = cp->yhigh;
X                    if (ylow < ymin) ymin = ylow;
X                    if (ylow > ymax) ymax = ylow;
X                    if (yhigh < ymin) ymin = yhigh;
X                    if (yhigh > ymax) ymax = yhigh;
X                }
X             }
X          } else {
X             cp->type = OUTRANGE;
X          }
X       }
X    }
}
X
X
get_3ddata(this_plot)
struct surface_points *this_plot;
{
register int i, j, l_num, xdatum, ydatum;
register FILE *fp;
float x, y, z;
BOOLEAN only_z = FALSE;
char format[MAX_LINE_LEN+1], data_file[MAX_LINE_LEN+1], line[MAX_LINE_LEN+1];
char *float_format = "%f", *float_skip = "%*f";
int xcol = 1, ycol = 2, zcol = 3, pt_in_iso_crv = 0, maxcol;
enum XYZ_order_type {XYZ, YXZ, ZXY, XZY, ZYX, YZX, XY, YX} xyz_order;
struct iso_curve *this_iso;
#ifdef AMIGA_LC_5_1
int num_perc_ast;
char *start_search;
#endif /* AMIGA_LC_5_1 */
X
X    quote_str(data_file, c_token);
X    this_plot->plot_type = DATA3D;
X    this_plot->has_grid_topology = TRUE;
X    if ((fp = fopen(data_file, "r")) == (FILE *)NULL)
X        os_error("can't open data file", c_token);
X
X    format[0] = '\0';
X    c_token++;    /* skip data file name */
X    if (almost_equals(c_token,"u$sing")) {
X        c_token++;      /* skip "using" */
X        if (!END_OF_COMMAND && !isstring(c_token)) {
X            struct value a;
X            zcol = (int)magnitude(const_express(&a));
X            only_z = TRUE;
X            if (equals(c_token,":")) {
X                c_token++;    /* skip ":" */
X                only_z = FALSE;
X                ycol = zcol;
X                zcol = (int)magnitude(const_express(&a));
X                if (equals(c_token,":")) {
X                    c_token++;    /* skip ":" */
X                    xcol = ycol;
X                    ycol = zcol;
X                    zcol = (int)magnitude(const_express(&a));
X                }
X                else {
X                    if (mapping3d == MAP3D_CARTESIAN)
X                        int_error("Must specify 1 or 3 columns",c_token);
X                    xcol = ycol;
X                    ycol = zcol;
X                }
X            }
X            if (!only_z)
X                if ( (xcol == ycol) || (ycol == zcol) || (xcol == zcol))
X                    int_error("Columns must be distinct",c_token);
X        }
X        if (!END_OF_COMMAND && isstring(c_token)) {
X            quotel_str(format, c_token);
X            c_token++;    /* skip format */
X        }
X    }
X
X    switch (mapping3d) {
X        case MAP3D_CARTESIAN:
X        maxcol = (xcol > ycol) ? xcol : ycol;
X        maxcol = (maxcol > zcol) ? maxcol : zcol;
X        if (!only_z) {    /* Determine ordering of input columns */
X             if (zcol == maxcol) {
X                 if (xcol < ycol)
X                     xyz_order = XYZ;  /* scanf(x,y,z) */
X                 else
X                     xyz_order = YXZ;  /* scanf(y,x,z) */
X             }
X             else if (ycol == maxcol) {
X                 if (xcol < zcol)
X                     xyz_order = XZY;  /* scanf(x,z,y) */
X                 else
X                     xyz_order = ZXY;  /* scanf(z,x,y) */
X             }
X             else {
X                 if (ycol < zcol)
X                     xyz_order = YZX;  /* scanf(y,z,x) */
X                 else
X                     xyz_order = ZYX;  /* scanf(z,y,x) */
X             }
X        }
X        if (strlen(format) == 0) {
X            if (only_z) {
X                for(i = 1; i <= zcol; i++)
X                    if (i == zcol)
X                        (void) strcat(format,float_format);
X                    else
X                        (void) strcat(format,float_skip);
X            }
X            else {
X                for(i = 1; i <= maxcol; i++)
X                    if ((i == xcol) || (i == ycol) || (i == zcol))
X                        (void) strcat(format,float_format);
X                    else
X                        (void) strcat(format,float_skip);
X            }
X        }
X            break;
X        case MAP3D_SPHERICAL:
X        case MAP3D_CYLINDRICAL:
X        if (only_z)
X            int_error("Two columns for spherical/cylindrical coords.",c_token);
X        maxcol = (xcol > ycol) ? xcol : ycol;
X        xyz_order = (xcol < ycol) ? XY : YX;
X        for(i = 1; i <= maxcol; i++)
X            if ((i == xcol) || (i == ycol))
X                (void) strcat(format,float_format);
X            else
X                (void) strcat(format,float_skip);
X    }
#ifdef AMIGA_LC_5_1
X    num_perc_ast = 0;
X    start_search = format;
X    while (*start_search != '\0') {
X        if (start_search[0] == '%')
X            if (start_search[1] == '*') num_perc_ast++;
X        start_search++;
X    }
#endif /* AMIGA_LC_5_1 */
X
X    l_num = 0;
X    xdatum = 0;
X    ydatum = 0;
X    this_plot->num_iso_read = 0;
X    this_plot->has_grid_topology = TRUE;
X    if ( this_plot->iso_crvs != NULL ) {
X        struct iso_curve *icrv, *icrvs = this_plot->iso_crvs;
X
X        while ( icrvs ) {
X        icrv = icrvs;
X        icrvs = icrvs->next;
X        iso_free( icrv );
X        }
X        this_plot->iso_crvs = NULL;
X    }
X    this_iso = iso_alloc( samples );
X
X    while ( fgets(line, MAX_LINE_LEN, fp) != (char *)NULL ) {
X        l_num++;
X        if (is_comment(line[0]))
X            continue;        /* ignore comments */
X        if (!line[1]) {            /* is it blank line ? */
X                if (pt_in_iso_crv == 0) {
X                    if (xdatum == 0)
X                    continue;
X                    pt_in_iso_crv = xdatum;
X            }
X
X            if (xdatum > 0) {
X                this_iso->p_count = xdatum;
X                this_iso->next = this_plot->iso_crvs;
X                this_plot->iso_crvs = this_iso;
X                this_plot->num_iso_read++;
X
X                if (xdatum != pt_in_iso_crv)
X                    this_plot->has_grid_topology = FALSE;
X
X                this_iso = iso_alloc(pt_in_iso_crv);
X                xdatum = 0;
X                ydatum++;
X            }
X            continue;
X        }
X
X        if (xdatum >= this_iso->p_max)
X        {
X            /* overflow about to occur. Extend size of points[]
X            * array. We either double the size, or add 1000 points,
X            * whichever is a smaller increment. Note i=p_max.
X            */
X            iso_extend(this_iso,
X                   xdatum + (xdatum < 1000 ? xdatum : 1000));
X        }
X
#ifdef AMIGA_LC_5_1
X        switch (sscanf(line, format, &x, &y, &z) - num_perc_ast) {
#else
X        switch (sscanf(line, format, &x, &y, &z)) {
#endif
X            case 3:         /* All parameter are specified. */
X               if (!only_z) {
X                switch (xyz_order) {
X                   case XYZ:    /* scanf(x,y,z) */
X                    this_iso->points[xdatum].x = x;
X                    this_iso->points[xdatum].y = y;
X                    this_iso->points[xdatum].z = z;
X                    break;
X                   case XZY:    /* scanf(x,z,y) */
X                    this_iso->points[xdatum].x = x;
X                    this_iso->points[xdatum].y = z;
X                    this_iso->points[xdatum].z = y;
X                    break;
X                   case YXZ:     /* scanf(y,x,z) */
X                    this_iso->points[xdatum].x = y;
X                    this_iso->points[xdatum].y = x;
X                    this_iso->points[xdatum].z = z;
X                    break;
X                   case ZXY:    /* scanf(z,x,y) */
X                    this_iso->points[xdatum].x = y;
X                    this_iso->points[xdatum].y = z;
X                    this_iso->points[xdatum].z = x;
X                    break;
X                   case YZX:    /* scanf(y,z,x) */
X                    this_iso->points[xdatum].x = z;
X                    this_iso->points[xdatum].y = x;
X                    this_iso->points[xdatum].z = y;
X                    break;
X                   case ZYX:    /* scanf(z,y,x) */
X                    this_iso->points[xdatum].x = z;
X                    this_iso->points[xdatum].y = y;
X                    this_iso->points[xdatum].z = x;
X                    break;
X                }
X                if (xyz_order != XYZ) {
X                    x = this_iso->points[xdatum].x;
X                    y = this_iso->points[xdatum].y;
X                    z = this_iso->points[xdatum].z;
X                }
X                if (!parametric)
X                    int_error("Must be in parametric mode.",
X                          NO_CARET);
X                   break;
X               }
X            case 1:         /* only one number on the line */
X               if (!only_z)
X                int_error("3 columns expected, only 1 found", c_token);
X               /* assign that number to z */
X               this_iso->points[xdatum].z = x;
X               z = x;
X               this_iso->points[xdatum].x = xdatum;
X               x = this_iso->points[xdatum].x;
X               this_iso->points[xdatum].y = ydatum;
X               y = this_iso->points[xdatum].y;
X               if (parametric)
X                int_error("Must be in non parametric mode.",
X                      NO_CARET);
X               break;
X            case 2:
X               switch (xyz_order) {
X                   case YX:
X                   z = x;    /* Use z as temp */
X                   x = y;
X                   y = z;
X                   break;
X                   default:
X                   break;
X               }
X               switch (mapping3d) {
X                   case MAP3D_CARTESIAN:
X                       int_error("2 columns found, 3 expected",
X                         c_token);
X                   break;
X                   case MAP3D_SPHERICAL:
X                   if (angles_format == ANGLES_DEGREES) {
X                       x *= DEG2RAD; /* Convert to radians. */
X                       y *= DEG2RAD;
X                   }
X                   this_iso->points[xdatum].x = cos(x) * cos(y);
X                   this_iso->points[xdatum].y = sin(x) * cos(y);
X                   this_iso->points[xdatum].z = sin(y);
X                   break;
X                   case MAP3D_CYLINDRICAL:
X                   if (angles_format == ANGLES_DEGREES)
X                       x *= DEG2RAD; /* Convert to radians. */
X                   this_iso->points[xdatum].x = cos(x);
X                   this_iso->points[xdatum].y = sin(x);
X                   this_iso->points[xdatum].z = y;
X                   break;
X               }
X               x = this_iso->points[xdatum].x;
X               y = this_iso->points[xdatum].y;
X               z = this_iso->points[xdatum].z;
X               break;
X            default:
X               (void) sprintf(line, "bad data on line %d", l_num);
X               int_error(line,c_token);
X        }
X
X        if (log_x) {
X            if (x < 0.0)
X            int_error("X value must be above 0 for log scale!",
X                  NO_CARET);
X            else
X            this_iso->points[xdatum].x =
X                log10(this_iso->points[xdatum].x);
X        }
X        if (log_y) {
X            if (y < 0.0)
X            int_error("Y value must be above 0 for log scale!",
X                  NO_CARET);
X            else
X            this_iso->points[xdatum].y =
X                log10(this_iso->points[xdatum].y);
X        }
X        if (log_z) {
X            if (z < 0.0)
X            int_error("Z value must be above 0 for log scale!",
X                  NO_CARET);
X            else
X            this_iso->points[xdatum].z =
X                log10(this_iso->points[xdatum].z);
X        }
X
X        if (autoscale_lx) {
X            if (x < xmin) xmin = x;
X            if (x > xmax) xmax = x;
X        }
X        if (autoscale_ly) {
X            if (y < ymin) ymin = y;
X            if (y > ymax) ymax = y;
X        }
X        if (autoscale_lz) {
X            if (z < zmin) zmin = z;
X            if (z > zmax) zmax = z;
X        }
X
X        xdatum++;
X    }
X
X    if (xdatum > 0) {
X        this_plot->num_iso_read++; /* Update last iso. */
X        this_iso->p_count = xdatum;
X
X        this_iso->next = this_plot->iso_crvs;
X        this_plot->iso_crvs = this_iso;
X
X        if (xdatum != pt_in_iso_crv)
X            this_plot->has_grid_topology = FALSE;
X
X    }
X    else {
X        iso_free(this_iso);/* Free last allocation. */
X    }
X
X    (void) fclose(fp);
X
X    if (this_plot->has_grid_topology) {
X            struct iso_curve *new_icrvs = NULL;
X        int num_new_iso = this_plot->iso_crvs->p_count,
X            len_new_iso = this_plot->num_iso_read;
X
X        /* Now we need to set the other direction (pseudo) isolines. */
X            for (i = 0; i < num_new_iso; i++) {
X            struct iso_curve *new_icrv = iso_alloc(len_new_iso);
X
X            new_icrv->p_count = len_new_iso;
X
X            for (j = 0, this_iso = this_plot->iso_crvs;
X                 this_iso != NULL;
X             j++, this_iso = this_iso->next) {
X            new_icrv->points[j].x = this_iso->points[i].x;
X                 new_icrv->points[j].y = this_iso->points[i].y;
X                 new_icrv->points[j].z = this_iso->points[i].z;
X            }
X
X            new_icrv->next = new_icrvs;
X            new_icrvs = new_icrv;
X        }
X
X        /* Append the new iso curves after the read ones. */
X        for (this_iso = this_plot->iso_crvs;
X             this_iso->next != NULL;
X             this_iso = this_iso->next);
X        this_iso->next = new_icrvs;
X    }
}
X
/* print_points:
X * a debugging routine to print out the points of a curve,
X * and the curve structure. If curve<0, then we print the
X * list of curves.
X */
static char *plot_type_names[4] = {
X    "Function", "Data", "3D Function", "3d data"
};
static char *plot_style_names[6] = {
X    "Lines", "Points", "Impulses", "LinesPoints", "Dots", "Errorbars"
};
X
print_points(curve)
X    int curve;            /* which curve to print */
{
X    register struct curve_points *this_plot;
X    int i;
X
X    if (curve < 0) {
X       for (this_plot = first_plot, i=0;
X           this_plot != NULL;
X           i++, this_plot = this_plot->next_cp) {
X          printf("Curve %d:\n", i);
X          if ((int)this_plot->plot_type >= 0 && (int)(this_plot->plot_type) < 4)
X            printf("Plot type %d: %s\n", (int)(this_plot->plot_type),
X                 plot_type_names[(int)(this_plot->plot_type)]);
X          else
X            printf("Plot type %d: BAD\n", (int)(this_plot->plot_type));
X          if ((int)this_plot->plot_style >= 0 && (int)(this_plot->plot_style) < 6)
X            printf("Plot style %d: %s\n", (int)(this_plot->plot_style),
X                 plot_style_names[(int)(this_plot->plot_style)]);
X          else
X            printf("Plot style %d: BAD\n", (int)(this_plot->plot_style));
X          printf("Plot title: '%s'\n", this_plot->title);
X          printf("Line type %d\n", this_plot->line_type);
X          printf("Point type %d\n", this_plot->point_type);
X          printf("max points %d\n", this_plot->p_max);
X          printf("current points %d\n", this_plot->p_count);
X          printf("\n");
X       }
X    } else {
X       for (this_plot = first_plot, i = 0;
X           i < curve && this_plot != NULL;
X           i++, this_plot = this_plot->next_cp)
X        ;
X       if (this_plot == NULL)
X        printf("Curve %d does not exist; list has %d curves\n", curve, i);
X       else {
X          printf ("Curve %d, %d points\n", curve, this_plot->p_count);
X          for (i = 0; i < this_plot->p_count; i++) {
X             printf("%c x=%g y=%g z=%g ylow=%g yhigh=%g\n",
X                   this_plot->points[i].type == INRANGE ? 'i'
X                   : this_plot->points[i].type == OUTRANGE ? 'o'
X                   : 'u',
X                   this_plot->points[i].x,
X                   this_plot->points[i].y,
X                   this_plot->points[i].z,
X                   this_plot->points[i].ylow,
X                   this_plot->points[i].yhigh);
X          }
X          printf("\n");
X       }
X    }
}
X
X
/* This parses the plot command after any range specifications.
X * To support autoscaling on the x axis, we want any data files to
X * define the x range, then to plot any functions using that range.
X * We thus parse the input twice, once to pick up the data files,
X * and again to pick up the functions. Definitions are processed
X * twice, but that won't hurt.
X */
eval_plots()
{
register int i;
register struct curve_points *this_plot, **tp_ptr;
register int start_token, end_token;
register int begin_token;
double x_min, x_max, y_min, y_max;
register double x, xdiff, temp;
static struct value a;
BOOLEAN ltmp, some_data_files = FALSE;
int plot_num, line_num, point_num, xparam=0;
char *xtitle;
void parametric_fixup();
X
X    if (autoscale_ly) {
X        ymin = VERYLARGE;
X        ymax = -VERYLARGE;
X    } else if (log_y && (ymin <= 0.0 || ymax <= 0.0))
X            int_error("y range must be above 0 for log scale!",
X                NO_CARET);
X
X    tp_ptr = &(first_plot);
X    plot_num = 0;
X    line_num = 0;     /* default line type */
X    point_num = 0;    /* default point type */
X
X    xtitle = NULL;
X
X    begin_token = c_token;
X
/*** First Pass: Read through data files ***/
/* This pass serves to set the xrange and to parse the command, as well
X * as filling in every thing except the function data. That is done after
X * the xrange is defined.
X */
X    while (TRUE) {
X        if (END_OF_COMMAND)
X            int_error("function to plot expected",c_token);
X
X        start_token = c_token;
X
X        if (is_definition(c_token)) {
X            define();
X        } else {
X            plot_num++;
X
X            if (isstring(c_token)) {            /* data file to plot */
X                if (parametric && xparam)
X                    int_error("previous parametric function not fully specified",
X                                                                    c_token);
X
X                if (!some_data_files && autoscale_lx) {
X                    xmin = VERYLARGE;
X                    xmax = -VERYLARGE;
X                }
X                some_data_files = TRUE;
X
X                if (*tp_ptr)
X                  this_plot = *tp_ptr;
X                else {        /* no memory malloc()'d there yet */
X                    this_plot = cp_alloc(MIN_CRV_POINTS);
X                    *tp_ptr = this_plot;
X                }
X                this_plot->plot_type = DATA;
X                this_plot->plot_style = data_style;
X                end_token = c_token;
X                get_data(this_plot); /* this also parses the using option */
X            }
X            else {                            /* function to plot */
X                if (parametric)            /* working on x parametric function */
X                    xparam = 1 - xparam;
X                if (*tp_ptr) {
X                    this_plot = *tp_ptr;
X                    cp_extend(this_plot, samples+1);
X                } else {        /* no memory malloc()'d there yet */
X                    this_plot = cp_alloc(samples+1);
X                    *tp_ptr = this_plot;
X                }
X                this_plot->plot_type = FUNC;
X                this_plot->plot_style = func_style;
X                dummy_func = &plot_func;
X                plot_func.at = temp_at();
X                /* ignore it for now */
X                end_token = c_token-1;
X            }
X
X            if (almost_equals(c_token,"t$itle")) {
X                if (parametric) {
X                    if (xparam)
X                        int_error(
X        "\"title\" allowed only after parametric function fully specified",
X                                                                    c_token);
X                    else if (xtitle != NULL)
X                        xtitle[0] = '\0';  /* Remove default title .*/
X                }
X                c_token++;
X                if ( isstring( c_token ) ) {
X                    m_quote_capture(&(this_plot->title),c_token,c_token);
X                }
X                else {
X                    int_error("expecting \"title\" for plot",c_token);
X                }
X                c_token++;
X            }
X              else {
X                  m_capture(&(this_plot->title),start_token,end_token);
X                 if (xparam) xtitle = this_plot->title;
X              }
X
X              this_plot->line_type = line_num;
X            this_plot->point_type = point_num;
X
X            if (almost_equals(c_token,"w$ith")) {
X                if (parametric && xparam)
X                    int_error("\"with\" allowed only after parametric function fully specified",
X                                    c_token);
X                this_plot->plot_style = get_style();
X            }
X
X            if ( !equals(c_token,",") && !END_OF_COMMAND ) {
X                struct value t;
X                this_plot->line_type = (int)real(const_express(&t))-1;
X            }
X            if ( !equals(c_token,",") && !END_OF_COMMAND ) {
X                struct value t;
X                this_plot->point_type = (int)real(const_express(&t))-1;
X            }
X            if ( (this_plot->plot_style == POINTS) ||
X                 (this_plot->plot_style == LINESPOINTS) ||
X                 (this_plot->plot_style == ERRORBARS) )
X                    if (!xparam) point_num++;
X            if (!xparam) line_num++;
X
X            if (this_plot->plot_type == DATA)
X              /* now that we know the plot style, adjust the yrange */
X              adjust_yrange(this_plot);
X
X            tp_ptr = &(this_plot->next_cp);
X        }
X
X        if (equals(c_token,","))
X            c_token++;
X        else
X            break;
X    }
X
X    if (parametric && xparam)
X        int_error("parametric function not fully specified", NO_CARET);
X
X    if (parametric) {
X    /* Swap t and x ranges for duration of these eval_plot computations. */
X        ltmp = autoscale_lx; autoscale_lx = autoscale_lt; autoscale_lt = ltmp;
X        temp = xmin; xmin = tmin; tmin = temp;
X        temp = xmax; xmax = tmax; tmax = temp;
X    }
X
/*** Second Pass: Evaluate the functions ***/
/* Everything is defined now, except the function data. We expect
X * no syntax errors, etc, since the above parsed it all. This makes
X * the code below simpler. If autoscale_ly, the yrange may still change.
X */
X     if (fabs(xmax-xmin) < zero)
X      if (autoscale_lx) {
X         fprintf(stderr, "Warning: empty %c range [%g:%g], ",
X            parametric ? 't' : 'x', xmin,xmax);
X         if (fabs(xmin) < zero) {
X            /* completely arbitary */
X            xmin = -1.;
X            xmax = 1.;
X         } else {
X            /* expand range by 10% in either direction */
X            xmin = xmin * 0.9;
X            xmax = xmax * 1.1;
X         }
X         fprintf(stderr, "adjusting to [%g:%g]\n", xmin,xmax);
X      } else {
X         int_error("x range is less than `zero`", c_token);
X      }
X
X    /* give error if xrange badly set from missing datafile error */
X    if (xmin == VERYLARGE || xmax == -VERYLARGE) {
X        int_error("x range is invalid", c_token);
X    }
X
X    if (log_x) {
X       if (xmin <= 0.0 || xmax <= 0.0)
X        int_error("x range must be greater than 0 for log scale!",NO_CARET);
X       x_min = log10(xmin);
X       x_max = log10(xmax);
X    } else {
X       x_min = xmin;
X       x_max = xmax;
X    }
X
X    xdiff = (x_max - x_min) / (samples - 1);
X
X    tp_ptr = &(first_plot);
X    plot_num = 0;
X    this_plot = first_plot;
X    c_token = begin_token;    /* start over */
X
X    /* Read through functions */
X    while (TRUE) {
X        if (is_definition(c_token)) {
X            define();
X        } else {
X            plot_num++;
X            if (isstring(c_token)) {            /* data file to plot */
X                /* ignore this now */
X                c_token++;
X                if (almost_equals(c_token,"u$sing")) {
X                    c_token++;      /* skip "using" */
X                        if (!isstring(c_token)) {
X                        struct value a;
X                        (void)magnitude(const_express(&a)); /* skip xcol */
X                        if (equals(c_token,":")) {
X                            c_token++;      /* skip ":" */
X                            (void)magnitude(const_express(&a)); /* skip ycol */
X                        }
X                        if (equals(c_token,":")) {
X                            c_token++;      /* skip ":" */
X                            (void)magnitude(const_express(&a)); /* skip yemin */
X                        }
X                        if (equals(c_token,":")) {
X                            c_token++;      /* skip ":" */
X                            (void)magnitude(const_express(&a)); /* skip yemax */
X                        }
X                    }
X                    if (isstring(c_token))
X                        c_token++;      /* skip format string */
X                }
X            }
X            else {                    /* function to plot */
X                if (parametric)            /* working on x parametric function */
X                    xparam = 1 - xparam;
X                dummy_func = &plot_func;
X                plot_func.at = temp_at(); /* reparse function */
X
X                for (i = 0; i < samples; i++) {
X                    x = x_min + i*xdiff;
X                    /* if (log_x) PEM fix logscale x axis */
X                    /* x = pow(10.0,x); 26-Sep-89 */
X                    (void) complex(&plot_func.dummy_values[0],
X                                log_x ? pow(10.0,x) : x,
X                                0.0);
X
X                    evaluate_at(plot_func.at,&a);
X
X                    if (undefined || (fabs(imag(&a)) > zero)) {
X                       this_plot->points[i].type = UNDEFINED;
X                       continue;
X                    }
X
X                    temp = real(&a);
X
X                    if (log_y && temp < 0.0) {
X                       this_plot->points[i].type = UNDEFINED;
X                       continue;
X                    }
X
X                    this_plot->points[i].x = x;
X                    if (log_y) {
X                       if (temp == 0.0) {
X                          this_plot->points[i].type = OUTRANGE;
X                          this_plot->points[i].y = -VERYLARGE;
X                          continue;
X                       } else {
X                          this_plot->points[i].y = log10(temp);
X                       }
X                    } else
X                     this_plot->points[i].y = temp;
X
X                    if (autoscale_ly || polar
X                       || inrange(temp, ymin, ymax)) {
X                       this_plot->points[i].type = INRANGE;
X                    /* When xparam is 1 we are not really computing y's! */
X                        if (!xparam && autoscale_ly) {
X                           if (temp < ymin) ymin = temp;
X                           if (temp > ymax) ymax = temp;
X                        }
X                    } else
X                     this_plot->points[i].type = OUTRANGE;
X                }
X                this_plot->p_count = i; /* samples */
X             }
X
X            /* title was handled above */
X            if (almost_equals(c_token,"t$itle")) {
X                c_token++;
X                c_token++;
X            }
X
X            /* style was handled above */
X            if (almost_equals(c_token,"w$ith")) {
X                c_token++;
X                c_token++;
X            }
X
X            /* line and point types were handled above */
X            if ( !equals(c_token,",") && !END_OF_COMMAND ) {
X                struct value t;
X                (void)real(const_express(&t));
X            }
X            if ( !equals(c_token,",") && !END_OF_COMMAND ) {
X                struct value t;
X                (void)real(const_express(&t));
X            }
X
X             tp_ptr = &(this_plot->next_cp); /* used below */
X            this_plot = this_plot->next_cp;
X         }
X
X        if (equals(c_token,","))
X          c_token++;
X        else
X          break;
X     }
X
X    /* throw out all curve_points at end of list, that we don't need  */
X    cp_free(*tp_ptr);
X    *tp_ptr = NULL;
X
X    if (fabs(ymax - ymin) < zero)
X     /* if autoscale, widen range */
X     if (autoscale_ly) {
X        fprintf(stderr, "Warning: empty y range [%g:%g], ", ymin, ymax);
X        if (fabs(ymin) < zero) {
X            ymin = -1.;
X            ymax = 1.;
X        } else {
X            /* expand range by 10% in either direction */
X            ymin = ymin * 0.9;
X            ymax = ymax * 1.1;
X        }
X        fprintf(stderr, "adjusting to [%g:%g]\n", ymin, ymax);
X     } else {
X        int_error("y range is less than `zero`", c_token);
X     }
X
/* Now we finally know the real ymin and ymax */
X    if (log_y) {
X        y_min = log10(ymin);
X        y_max = log10(ymax);
X    } else {
X        y_min = ymin;
X        y_max = ymax;
X    }
X    capture(replot_line,plot_token,c_token);
X
X    if (parametric) {
X    /* Now put t and x ranges back before we actually plot anything. */
X        ltmp = autoscale_lx; autoscale_lx = autoscale_lt; autoscale_lt = ltmp;
X        temp = xmin; xmin = tmin; tmin = temp;
X        temp = xmax; xmax = tmax; tmax = temp;
X        if (some_data_files && autoscale_lx) {
X        /*
X            Stop any further autoscaling in this case (may be a mistake, have
X              to consider what is really wanted some day in the future--jdc).
X        */
X            autoscale_lx = 0;
X        }
X    /* Now actually fix the plot pairs to be single plots. */
X        parametric_fixup (first_plot, &plot_num, &x_min, &x_max);
X    }
X
X    do_plot(first_plot,plot_num,x_min,x_max,y_min,y_max);
X    cp_free(first_plot);
X    first_plot = NULL;
}
X
/* This parses the splot command after any range specifications.
X * To support autoscaling on the x/z axis, we want any data files to
X * define the x/y range, then to plot any functions using that range.
X * We thus parse the input twice, once to pick up the data files,
X * and again to pick up the functions. Definitions are processed
X * twice, but that won't hurt.
X */
eval_3dplots()
{
register int i,j,k;
register struct surface_points *this_plot, **tp_3d_ptr;
register int start_token, end_token;
register int begin_token;
double x_min, x_max, y_min, y_max, z_min, z_max;
register double x, xdiff, xisodiff, y, ydiff, yisodiff, temp;
static struct value a;
BOOLEAN ltmp, some_data_files = FALSE;
int plot_num, line_num, point_num,
X    crnt_param = 0; /* 0=x, 1=y, 2=z */
char *xtitle;
char *ytitle;
void parametric_3dfixup();
X
X    if (autoscale_lz) {
X        zmin = VERYLARGE;
X        zmax = -VERYLARGE;
X    } else if (log_z && (zmin <= 0.0 || zmax <= 0.0))
X            int_error("z range must be above 0 for log scale!",
X                NO_CARET);
X
X    tp_3d_ptr = &(first_3dplot);
X    plot_num = 0;
X    line_num = 0;     /* default line type */
X    point_num = 0;    /* default point type */
X
X    xtitle = NULL;
X    ytitle = NULL;
X
X    begin_token = c_token;
X
/*** First Pass: Read through data files ***/
/* This pass serves to set the x/yranges and to parse the command, as well
X * as filling in every thing except the function data. That is done after
X * the x/yrange is defined.
X */
X    while (TRUE) {
X        if (END_OF_COMMAND)
X            int_error("function to plt3d expected",c_token);
X
X        start_token = c_token;
X
X        if (is_definition(c_token)) {
X            define();
X        } else {
X            plot_num++;
X
X            if (isstring(c_token)) {            /* data file to plot */
X                if (parametric && crnt_param != 0)
X                    int_error("previous parametric function not fully specified",
X                                                                    c_token);
X
X                if (!some_data_files) {
X                    if (autoscale_lx) {
X                        xmin = VERYLARGE;
X                        xmax = -VERYLARGE;
X                    }
X                    if (autoscale_ly) {
X                        ymin = VERYLARGE;
X                        ymax = -VERYLARGE;
X                    }
X                }
X
X                some_data_files = TRUE;
X
X                if (*tp_3d_ptr)
X                    this_plot = *tp_3d_ptr;
X                else {    /* no memory malloc()'d there yet */
X                    /* Allocate samples * iso_samples twice for */
X                    /* Each of the isoparametric direction. */
X                    this_plot = sp_alloc(0,0);
X                    *tp_3d_ptr = this_plot;
X                }
X
X                this_plot->plot_type = DATA3D;
X                this_plot->plot_style = data_style;
X                end_token = c_token;
X                get_3ddata(this_plot); /* this also parses the using option */
X            }
X            else {                        /* function to plot */
X                if (parametric) /* Rotate between x/y/z axes */
X                    crnt_param = (crnt_param+1) % 3;
X                if (*tp_3d_ptr) {
X                    this_plot = *tp_3d_ptr;
X                    sp_replace(this_plot,samples,2*iso_samples);
X                }
X                else {    /* no memory malloc()'d there yet */
X                    /* Allocate samples * iso_samples twice for */
X                    /* Each of the isoparametric direction. */
X                    this_plot = sp_alloc(samples,2*iso_samples);
X                    *tp_3d_ptr = this_plot;
X                }
X
X                this_plot->plot_type = FUNC3D;
X                this_plot->has_grid_topology = TRUE;
X                this_plot->plot_style = func_style;
X                dummy_func = &plot_func;
X                plot_func.at = temp_at();
X                /* ignore it for now */
X                end_token = c_token-1;
X            }
X
X            if (almost_equals(c_token,"t$itle")) {
X                if (parametric) {
X                        if (crnt_param)
X                        int_error(
X        "\"title\" allowed only after parametric function fully specified",
X                                                                    c_token);
X                    else {
X                        /* Remove default title */
X                        if (xtitle != NULL)
X                            xtitle[0] = '\0';
X                        if (ytitle != NULL)
X                            ytitle[0] = '\0';
X                    }
X                    }
X                c_token++;
X                if ( isstring( c_token ) ) {
X                    m_quote_capture(&(this_plot->title),c_token,c_token);
X                }
X                else {
X                    int_error("expecting \"title\" for plot",c_token);
X                }
X                c_token++;
X            }
X              else {
X                  m_capture(&(this_plot->title),start_token,end_token);
X                if (crnt_param == 1) xtitle = this_plot->title;
X                if (crnt_param == 2) ytitle = this_plot->title;
X              }
X
X              this_plot->line_type = line_num;
X            this_plot->point_type = point_num;
X
X            if (almost_equals(c_token,"w$ith")) {
X                this_plot->plot_style = get_style();
X            }
X
X            if ( !equals(c_token,",") && !END_OF_COMMAND ) {
X                struct value t;
X                this_plot->line_type = (int)real(const_express(&t))-1;
X            }
X            if ( !equals(c_token,",") && !END_OF_COMMAND ) {
X                struct value t;
X                this_plot->point_type = (int)real(const_express(&t))-1;
X            }
X            if ( (this_plot->plot_style == POINTS) ||
X                 (this_plot->plot_style == LINESPOINTS) ||
X                 (this_plot->plot_style == ERRORBARS) )
X                    if (crnt_param == 0)
X                        point_num +=
X                            1 + (draw_contour != 0);
X            if (crnt_param == 0)
X                line_num += 1 + (draw_contour != 0);
X
X            tp_3d_ptr = &(this_plot->next_sp);
X        }
X
X        if (equals(c_token,","))
X            c_token++;
X        else
X            break;
X    }
X
X    if (parametric && crnt_param != 0)
X        int_error("parametric function not fully specified", NO_CARET);
X
X    if (parametric) {
X    /* Swap u/v and x/y ranges for duration of these eval_plot computations. */
X        ltmp = autoscale_lx; autoscale_lx = autoscale_lu; autoscale_lu = ltmp;
X        ltmp = autoscale_ly; autoscale_ly = autoscale_lv; autoscale_lv = ltmp;
X        temp = xmin; xmin = umin; umin = temp;
X        temp = xmax; xmax = umax; umax = temp;
X        temp = ymin; ymin = vmin; vmin = temp;
X        temp = ymax; ymax = vmax; vmax = temp;
X    }
/*** Second Pass: Evaluate the functions ***/
/* Everything is defined now, except the function data. We expect
X * no syntax errors, etc, since the above parsed it all. This makes
X * the code below simpler. If autoscale_ly, the yrange may still change.
X */
X     if (xmin == xmax)
X      if (autoscale_lx) {
X         fprintf(stderr, "Warning: empty x range [%g:%g], ",
X             xmin,xmax);
X         if (xmin == 0.0) {
X            /* completely arbitary */
X            xmin = -1.;
X            xmax = 1.;
X         } else {
X            /* expand range by 10% in either direction */
X            xmin = xmin * 0.9;
X            xmax = xmax * 1.1;
X         }
X         fprintf(stderr, "adjusting to [%g:%g]\n", xmin,xmax);
X      } else {
X         int_error("x range is empty", c_token);
X      }
X
X     if (ymin == ymax)
X      if (autoscale_ly) {
X         fprintf(stderr, "Warning: empty y range [%g:%g], ",
X            ymin,ymax);
X         if (ymin == 0.0) {
X            /* completely arbitary */
X            ymin = -1.;
X            ymax = 1.;
X         } else {
X            /* expand range by 10% in either direction */
X            ymin = ymin * 0.9;
X            ymax = ymax * 1.1;
X         }
X         fprintf(stderr, "adjusting to [%g:%g]\n", ymin,ymax);
X      } else {
X         int_error("y range is empty", c_token);
X      }
X
X    /* give error if xrange badly set from missing datafile error */
X    if (xmin == VERYLARGE || xmax == -VERYLARGE) {
X    int_error("x range is invalid", c_token);
X    }
X
X    if (log_x) {
X       if (xmin <= 0.0 || xmax <= 0.0)
X        int_error("x range must be greater than 0 for log scale!",NO_CARET);
X       x_min = log10(xmin);
X       x_max = log10(xmax);
X    } else {
X       x_min = xmin;
X       x_max = xmax;
X    }
X
X    if (log_y) {
X       if (ymin <= 0.0 || ymax <= 0.0)
X        int_error("y range must be greater than 0 for log scale!",NO_CARET);
X       y_min = log10(ymin);
X       y_max = log10(ymax);
X    } else {
X       y_min = ymin;
X       y_max = ymax;
X    }
X
X    if (samples < 2 || iso_samples < 2)
X    int_error("samples or iso_samples < 2. Must be at least 2.\n");
X
X    xdiff = (x_max - x_min) / (samples - 1);
X    ydiff = (y_max - y_min) / (samples - 1);
X    xisodiff = (x_max - x_min) / (iso_samples - 1);
X    yisodiff = (y_max - y_min) / (iso_samples - 1);
X
X    plot_num = 0;
X    this_plot = first_3dplot;
X    c_token = begin_token;    /* start over */
X
X    /* Read through functions */
X    while (TRUE) {
X        if (is_definition(c_token)) {
X            define();
X        } else {
X            plot_num++;
X            if (isstring(c_token)) {            /* data file to plot */
X                /* ignore this now */
X                c_token++;
X                if (almost_equals(c_token,"u$sing")) {
X                    c_token++;      /* skip "using" */
X                                if (!isstring(c_token)) {
X                        struct value a;
X                        (void)magnitude(const_express(&a));    /* skip xcol */
X                        if (equals(c_token,":")) {
X                            c_token++;    /* skip ":" */
X                            (void)magnitude(const_express(&a));    /* skip ycol */
X                            if (equals(c_token,":")) {
X                                c_token++;    /* skip ":" */
X                                (void)magnitude(const_express(&a));    /* skip zcol */
X                            }
X                        }
X                    }
X                                if (isstring(c_token))
X                        c_token++;    /* skip format string */
X                }
X            }
X            else {                    /* function to plot */
X                struct iso_curve *this_iso = this_plot->iso_crvs;
X                struct coordinate *points = this_iso->points;
X
X                if (parametric)
X                    crnt_param = (crnt_param+1) % 3;
X                dummy_func = &plot_func;
X                plot_func.at = temp_at(); /* reparse function */
X
X                for (j = 0; j < iso_samples; j++) {
X                    y = y_min + j*yisodiff;
X                    /* if (log_y) PEM fix logscale y axis */
X                    /* y = pow(10.0,y); 26-Sep-89 */
X                    (void) complex(&plot_func.dummy_values[1],
X                                log_y ? pow(10.0,y) : y,
X                                0.0);
X
X                    for (i = 0; i < samples; i++) {
X                        x = x_min + i*xdiff;
X                        /* if (log_x) PEM fix logscale x axis */
X                        /* x = pow(10.0,x); 26-Sep-89 */
X                        (void) complex(&plot_func.dummy_values[0],
X                                log_x ? pow(10.0,x) : x,
X                                0.0);
X
X                        points[i].x = x;
X                        points[i].y = y;
X
X                        evaluate_at(plot_func.at,&a);
X
X                        if (undefined || (fabs(imag(&a)) > zero)) {
X                       points[i].type = UNDEFINED;
X                       continue;
X                        }
X
X                        temp = real(&a);
X
X                        if (log_z && temp < 0.0) {
X                       points[i].type = UNDEFINED;
X                       continue;
X                        }
X
X                        if (log_z) {
X                       if (temp == 0.0) {
X                          points[i].type = OUTRANGE;
X                          points[i].z = -VERYLARGE;
X                          continue;
X                       } else {
X                          points[i].z = log10(temp);
X                       }
X                        } else
X                       points[i].z = temp;
X
X                        if (autoscale_lz
X                       || inrange(temp, zmin, zmax)) {
X                       points[i].type = INRANGE;
X                       if (autoscale_lz) {
X                          if (temp < zmin) zmin = temp;
X                          if (temp > zmax) zmax = temp;
X                       }
X                        } else
X                       points[i].type = OUTRANGE;
X                    }
X                    this_iso->p_count = samples;
X                    this_iso = this_iso->next;
X                    points = this_iso->points;
X                }
X
X                for (i = 0; i < iso_samples; i++) {
X                    x = x_min + i*xisodiff;
X                    /* if (log_x) PEM fix logscale x axis */
X                    /* x = pow(10.0,x); 26-Sep-89 */
X                    (void) complex(&plot_func.dummy_values[0],
X                                log_x ? pow(10.0,x) : x,
X                                0.0);
X
X                    for (j = 0; j < samples; j++) {
X                        y = y_min + j*ydiff;
X                        /* if (log_y) PEM fix logscale y axis */
X                        /* y = pow(10.0,y); 26-Sep-89 */
X                        (void) complex(&plot_func.dummy_values[1],
X                                log_y ? pow(10.0,y) : y,
X                                0.0);
X
X                        points[j].x = x;
X                        points[j].y = y;
X
X                        evaluate_at(plot_func.at,&a);
X
X                        if (undefined || (fabs(imag(&a)) > zero)) {
X                       points[j].type = UNDEFINED;
X                       continue;
X                        }
X
X                        temp = real(&a);
X
X                        if (log_z && temp < 0.0) {
X                       points[j].type = UNDEFINED;
X                       continue;
X                        }
X
X                        if (log_z) {
X                       if (temp == 0.0) {
X                          points[j].type = OUTRANGE;
X                          points[j].z = -VERYLARGE;
X                          continue;
X                       } else {
X                          points[j].z = log10(temp);
X                       }
X                        } else
X                       points[j].z = temp;
X
X                        if (autoscale_lz
X                       || inrange(temp, zmin, zmax)) {
X                       points[j].type = INRANGE;
X                       if (autoscale_lz) {
X                          if (temp < zmin) zmin = temp;
X                          if (temp > zmax) zmax = temp;
X                       }
X                        } else
X                       points[j].type = OUTRANGE;
X                    }
X                        this_iso->p_count = samples;
X                    this_iso = this_iso->next;
X                    points = this_iso ? this_iso->points : NULL;
X                }
X             }
X
X            /* title was handled above */
X            if (almost_equals(c_token,"t$itle")) {
X                c_token++;
X                c_token++;
X            }
X
X            /* style was handled above */
X            if (almost_equals(c_token,"w$ith")) {
X                c_token++;
X                c_token++;
X            }
X
X            /* line and point types were handled above */
X            if ( !equals(c_token,",") && !END_OF_COMMAND ) {
X                struct value t;
X                (void)real(const_express(&t));
X            }
X            if ( !equals(c_token,",") && !END_OF_COMMAND ) {
X                struct value t;
X                (void)real(const_express(&t));
X            }
X
X            this_plot = this_plot->next_sp;
X         }
X
X        if (equals(c_token,","))
X          c_token++;
X        else
X          break;
X     }
X
X    if (fabs(zmax - zmin) < zero)
X     /* if autoscale, widen range */
X     if (autoscale_lz) {
X        fprintf(stderr, "Warning: empty z range [%g:%g], ", zmin, zmax);
X        if (fabs(zmin) < zero ) {
X            zmin = -1.;
X            zmax = 1.;
X        } else {
X            /* expand range by 10% in either direction */
X            zmin = zmin * 0.9;
X            zmax = zmax * 1.1;
X        }
X        fprintf(stderr, "adjusting to [%g:%g]\n", zmin, zmax);
X     } else {
X        int_error("z range is less than `zero`", c_token);
X     }
X
/* Now we finally know the real zmin and zmax */
X    if (log_z) {
X        if (zmin <= 0.0 || zmax <= 0.0)
X            int_error("z range must be greater than 0 for log scale!",NO_CARET);
X        z_min = log10(zmin);
X        z_max = log10(zmax);
X    } else {
X        z_min = zmin;
X        z_max = zmax;
X    }
X    capture(replot_line,plot_token,c_token);
X
X    if (parametric) {
X    /* Now put u/v and x/y ranges back before we actually plot anything. */
X        ltmp = autoscale_lx; autoscale_lx = autoscale_lu; autoscale_lu = ltmp;
X        ltmp = autoscale_ly; autoscale_ly = autoscale_lv; autoscale_lv = ltmp;
X        temp = xmin; xmin = umin; umin = temp;
X        temp = xmax; xmax = umax; umax = temp;
X        temp = ymin; ymin = vmin; vmin = temp;
X        temp = ymax; ymax = vmax; vmax = temp;
X
X    /* Now actually fix the plot triplets to be single plots. */
X        parametric_3dfixup(first_3dplot, &plot_num,
X                   &x_min, &x_max, &y_min, &y_max,
X                   &z_min, &z_max);
X        if (log_x) {
X            if (x_min <= 0.0 || x_max <= 0.0)
X                int_error("x range must be greater than 0 for log scale!",NO_CARET);
X            x_min = log10(x_min);
X            x_max = log10(x_max);
X        }
X
X        if (log_y) {
X            if (y_min <= 0.0 || y_max <= 0.0)
X                int_error("y range must be greater than 0 for log scale!",NO_CARET);
X            y_min = log10(y_min);
X            y_max = log10(y_max);
X        }
X
X        if (log_z) {
X            if (z_min <= 0.0 || z_max <= 0.0)
X                int_error("z range must be greater than 0 for log scale!",NO_CARET);
X            z_min = log10(z_min);
X            z_max = log10(z_max);
X        }
X    }
X
X
X    /* Creates contours if contours are to be plotted as well. */
X    if (draw_contour) {
X        for (this_plot=first_3dplot, i=0;
X             i < plot_num;
X             this_plot=this_plot->next_sp, i++) {
X            if (this_plot->contours) {
X                struct gnuplot_contours *cntr, *cntrs = this_plot->contours;
X
X                while (cntrs) {
X                    cntr = cntrs;
X                    cntrs = cntrs->next;
X                    free(cntr->coords);
X                    free(cntr);
X                }
X            }
X            /* Make sure this one can be contoured. */
X            if (!this_plot->has_grid_topology) {
X                this_plot->contours = NULL;
X                int_error("Can not contour non grid data!",NO_CARET);
X            }
X            if (this_plot->plot_type == DATA3D)
X                this_plot->contours = contour(
X                           this_plot->num_iso_read,
X                           this_plot->iso_crvs,
X                           contour_levels, contour_pts,
X                           contour_kind, contour_order);
X            else
X                this_plot->contours = contour(iso_samples,
X                           this_plot->iso_crvs,
X                           contour_levels, contour_pts,
X                           contour_kind, contour_order);
X        }
X    }
X
X    do_3dplot(first_3dplot,plot_num,x_min,x_max,y_min,y_max,z_min,z_max);
X    sp_free(first_3dplot);
X    first_3dplot = NULL;
}
X
done(status)
int status;
{
X    if (term && term_init)
X        (*term_tbl[term].reset)();
#ifdef vms
X    vms_reset();
#endif
X    exit(status);
}
X
void parametric_fixup (start_plot, plot_num, x_min, x_max)
struct curve_points *start_plot;
int *plot_num;
double *x_min, *x_max;
/*
X    The hardest part of this routine is collapsing the FUNC plot types
X   in the list (which are gauranteed to occur in (x,y) pairs while
X    preserving the non-FUNC type plots intact.  This means we have to
X    work our way through various lists.  Examples (hand checked):
X        start_plot:F1->F2->NULL ==> F2->NULL
X        start_plot:F1->F2->F3->F4->F5->F6->NULL ==> F2->F4->F6->NULL
X        start_plot:F1->F2->D1->D2->F3->F4->D3->NULL ==> F2->D1->D2->F4->D3->NULL
SHAR_EOF
true || echo 'restore of gnuplot/command.c failed'
fi
echo 'End of  part 13'
echo 'File gnuplot/command.c is continued in part 14'
echo 14 > _shar_seq_.tmp
exit 0

exit 0 # Just in case...
--
Kent Landfield                   INTERNET: [e-mail]
Sterling Software, IMD           UUCP:     uunet!sparky!kent
Phone:    [phone]         FAX:      [phone]
Please send comp.sources.misc-related mail to [e-mail].