#ifndef lint
static char *RCSid() { return RCSid("$Id: plot3d.c,v 1.20.2.6 2000/10/24 18:58:12 broeker Exp $"); }
#endif
/* GNUPLOT - plot3d.c */
/*[
* Copyright 1986 - 1993, 1998 Thomas Williams, Colin Kelley
*
* Permission to use, copy, and distribute this software and its
* documentation for any purpose with or without fee is hereby granted,
* provided that the above copyright notice appear in all copies and
* that both that copyright notice and this permission notice appear
* in supporting documentation.
*
* Permission to modify the software is granted, but not the right to
* distribute the complete modified source code. Modifications are to
* be distributed as patches to the released version. Permission to
* distribute binaries produced by compiling modified sources is granted,
* provided you
* 1. distribute the corresponding source modifications from the
* released version in the form of a patch file along with the binaries,
* 2. add special version identification to distinguish your version
* in addition to the base release version number,
* 3. provide your name and address as the primary contact for the
* support of your modified version, and
* 4. retain our contact information in regard to use of the base
* software.
* Permission to distribute the released version of the source code along
* with corresponding source modifications in the form of a patch file is
* granted with same provisions 2 through 4 for binary distributions.
*
* This software is provided "as is" without express or implied warranty
* to the extent permitted by applicable law.
]*/
#include "plot3d.h"
#include "gp_types.h"
#include "alloc.h"
#include "axis.h"
#include "binary.h"
#include "command.h"
#include "contour.h"
#include "datafile.h"
#include "graph3d.h"
#include "internal.h"
#include "misc.h"
#include "parse.h"
/* #include "setshow.h" */
#include "term_api.h"
#include "util.h"
#ifdef THIN_PLATE_SPLINES_GRID
#include "matrix.h"
#endif
#ifndef _Windows
# include "help.h"
#endif
/* global variables exported by this module */
t_data_mapping mapping3d = MAP3D_CARTESIAN;
int dgrid3d_row_fineness = 10;
int dgrid3d_col_fineness = 10;
int dgrid3d_norm_value = 1;
TBOOLEAN dgrid3d = FALSE;
/* static prototypes */
static void calculate_set_of_isolines __PROTO((AXIS_INDEX value_axis, TBOOLEAN cross, struct iso_curve **this_iso,
AXIS_INDEX iso_axis, double iso_min, double iso_step, int num_iso_to_use,
AXIS_INDEX sam_axis, double sam_min, double sam_step, int num_sam_to_use));
static void get_3ddata __PROTO((struct surface_points * this_plot));
static void print_3dtable __PROTO((int pcount));
static void eval_3dplots __PROTO((void));
static void grid_nongrid_data __PROTO((struct surface_points * this_plot));
static void parametric_3dfixup __PROTO((struct surface_points * start_plot, int *plot_num));
static struct surface_points * sp_alloc __PROTO((int num_samp_1, int num_iso_1, int num_samp_2, int num_iso_2));
static void sp_replace __PROTO((struct surface_points *sp, int num_samp_1, int num_iso_1, int num_samp_2, int num_iso_2));
/* the curves/surfaces of the plot */
struct surface_points *first_3dplot = NULL;
static struct udft_entry plot_func;
int plot3d_num=0;
/* HBB 20000508: moved these functions to the only module that uses them
* so they can be turned 'static' */
/*
* sp_alloc() allocates a surface_points structure that can hold 'num_iso_1'
* iso-curves with 'num_samp_2' samples and 'num_iso_2' iso-curves with
* 'num_samp_1' samples.
* If, however num_iso_2 or num_samp_1 is zero no iso curves are allocated.
*/
static struct surface_points *
sp_alloc(num_samp_1, num_iso_1, num_samp_2, num_iso_2)
int num_samp_1, num_iso_1, num_samp_2, num_iso_2;
{
struct surface_points *sp;
sp = (struct surface_points *) gp_alloc(sizeof(struct surface_points),
"surface");
sp->next_sp = NULL;
sp->title = NULL;
sp->contours = NULL;
sp->iso_crvs = NULL;
sp->num_iso_read = 0;
if (num_iso_2 > 0 && num_samp_1 > 0) {
int i;
struct iso_curve *icrv;
for (i = 0; i < num_iso_1; i++) {
icrv = iso_alloc(num_samp_2);
icrv->next = sp->iso_crvs;
sp->iso_crvs = icrv;
}
for (i = 0; i < num_iso_2; i++) {
icrv = iso_alloc(num_samp_1);
icrv->next = sp->iso_crvs;
sp->iso_crvs = icrv;
}
} else
sp->iso_crvs = (struct iso_curve *) NULL;
return (sp);
}
/*
* sp_replace() updates a surface_points structure so it can hold 'num_iso_1'
* iso-curves with 'num_samp_2' samples and 'num_iso_2' iso-curves with
* 'num_samp_1' samples.
* If, however num_iso_2 or num_samp_1 is zero no iso curves are allocated.
*/
static void
sp_replace(sp, num_samp_1, num_iso_1, num_samp_2, num_iso_2)
struct surface_points *sp;
int num_samp_1, num_iso_1, num_samp_2, num_iso_2;
{
int i;
struct iso_curve *icrv, *icrvs = sp->iso_crvs;
while (icrvs) {
icrv = icrvs;
icrvs = icrvs->next;
iso_free(icrv);
}
sp->iso_crvs = NULL;
if (num_iso_2 > 0 && num_samp_1 > 0) {
for (i = 0; i < num_iso_1; i++) {
icrv = iso_alloc(num_samp_2);
icrv->next = sp->iso_crvs;
sp->iso_crvs = icrv;
}
for (i = 0; i < num_iso_2; i++) {
icrv = iso_alloc(num_samp_1);
icrv->next = sp->iso_crvs;
sp->iso_crvs = icrv;
}
} else
sp->iso_crvs = (struct iso_curve *) NULL;
}
/*
* sp_free() releases any memory which was previously malloc()'d to hold
* surface points.
*/
/* HBB 20000506: don't risk stack havoc by recursion, use iterative list
* cleanup unstead */
void
sp_free(sp)
struct surface_points *sp;
{
while (sp) {
struct surface_points *next = sp->next_sp;
if (sp->title)
free(sp->title);
while (sp->contours) {
struct gnuplot_contours *next_cntrs = sp->contours->next;
free(sp->contours->coords);
free(sp->contours);
sp->contours = next_cntrs;
}
while (sp->iso_crvs) {
struct iso_curve *next_icrvs = sp->iso_crvs->next;
iso_free(sp->iso_crvs);
sp->iso_crvs = next_icrvs;
}
free(sp);
sp = next;
}
}
/* support for dynamic size of input line */
void
plot3drequest()
/*
* in the parametric case we would say splot [u= -Pi:Pi] [v= 0:2*Pi] [-1:1]
* [-1:1] [-1:1] sin(v)*cos(u),sin(v)*cos(u),sin(u) in the non-parametric
* case we would say only splot [x= -2:2] [y= -5:5] sin(x)*cos(y)
*
*/
{
int dummy_token0 = -1, dummy_token1 = -1;
AXIS_INDEX u_axis, v_axis;
is_3d_plot = TRUE;
if (parametric && strcmp(set_dummy_var[0], "t") == 0) {
strcpy(set_dummy_var[0], "u");
strcpy(set_dummy_var[1], "v");
}
/* put stuff into arrays to simplify access */
AXIS_INIT3D(FIRST_X_AXIS, 0, 0);
AXIS_INIT3D(FIRST_Y_AXIS, 0, 0);
AXIS_INIT3D(FIRST_Z_AXIS, 0, 1);
AXIS_INIT3D(U_AXIS, 1, 0);
AXIS_INIT3D(V_AXIS, 1, 0);
if (!term) /* unknown */
int_error(c_token, "use 'set term' to set terminal type first");
u_axis = (parametric ? U_AXIS : FIRST_X_AXIS);
v_axis = (parametric ? V_AXIS : FIRST_Y_AXIS);
PARSE_NAMED_RANGE(u_axis, dummy_token0);
PARSE_NAMED_RANGE(v_axis, dummy_token1);
if (parametric) {
PARSE_RANGE(FIRST_X_AXIS);
PARSE_RANGE(FIRST_Y_AXIS);
} /* parametric */
PARSE_RANGE(FIRST_Z_AXIS);
CHECK_REVERSE(FIRST_X_AXIS);
CHECK_REVERSE(FIRST_Y_AXIS);
CHECK_REVERSE(FIRST_Z_AXIS);
/* use the default dummy variable unless changed */
if (dummy_token0 >= 0)
copy_str(c_dummy_var[0], dummy_token0, MAX_ID_LEN);
else
(void) strcpy(c_dummy_var[0], set_dummy_var[0]);
if (dummy_token1 >= 0)
copy_str(c_dummy_var[1], dummy_token1, MAX_ID_LEN);
else
(void) strcpy(c_dummy_var[1], set_dummy_var[1]);
eval_3dplots();
}
#ifdef THIN_PLATE_SPLINES_GRID
static double splines_kernel __PROTO((double h));
/* HBB 991025 FIXME: these don't belong in here --> move to 'matrix' */
static void lu_decomp __PROTO((double **, int, int *, double *));
static void lu_backsubst __PROTO((double **, int n, int *, double *));
static double
splines_kernel(h)
double h;
{
/* this is normaly not usefull ... */
h = fabs(h);
if (h != 0.0) {
return h * h * log(h);
} else {
return 0;
}
}
#define Swap(a,b) {double tmp; tmp=a; a=b; b=tmp;}
static void
lu_decomp(a, n, indx, d)
double **a;
int n;
int *indx;
double *d;
{
int i, imax = -1, j, k; /* HBB: added initial value, to shut up gcc -Wall */
double large, dummy, temp, **ar, **lim, *limc, *ac, *dp, *vscal;
dp = vscal = vec(n);
*d = 1.0;
for (ar = a, lim = &(a[n]); ar < lim; ar++) {
large = 0.0;
for (ac = *ar, limc = &(ac[n]); ac < limc;)
if ((temp = fabs(*ac++)) > large)
large = temp;
if (large == 0.0)
int_error(NO_CARET, "Singular matrix in LU-DECOMP");
*dp++ = 1 / large;
}
ar = a;
for (j = 0; j < n; j++, ar++) {
for (i = 0; i < j; i++) {
ac = &(a[i][j]);
for (k = 0; k < i; k++)
*ac -= a[i][k] * a[k][j];
}
large = 0.0;
dp = &(vscal[j]);
for (i = j; i < n; i++) {
ac = &(a[i][j]);
for (k = 0; k < j; k++)
*ac -= a[i][k] * a[k][j];
if ((dummy = *dp++ * fabs(*ac)) >= large) {
large = dummy;
imax = i;
}
}
if (j != imax) {
ac = a[imax];
dp = *ar;
for (k = 0; k < n; k++, ac++, dp++)
Swap(*ac, *dp);
*d = -(*d);
vscal[imax] = vscal[j];
}
indx[j] = imax;
if (*(dp = &(*ar)[j]) == 0)
*dp = 1e-30;
if (j != n - 1) {
dummy = 1 / (*ar)[j];
for (i = j + 1; i < n; i++)
a[i][j] *= dummy;
}
}
free(vscal);
}
static void
lu_backsubst(a, n, indx, b)
double **a;
int n;
int *indx;
double *b;
{
int i, memi = -1, ip, j;
double sum, *bp, *bip, **ar, *ac;
ar = a;
for (i = 0; i < n; i++, ar++) {
ip = indx[i];
sum = b[ip];
b[ip] = b[i];
if (memi >= 0) {
ac = &((*ar)[memi]);
bp = &(b[memi]);
for (j = memi; j = 0; i--) {
ac = &(*ar)[i + 1];
bp = &(b[i + 1]);
bip = &(b[i]);
for (j = i + 1; j < n; j++)
*bip -= *ac++ * *bp++;
*bip /= (*ar--)[i];
}
}
#endif
static void
grid_nongrid_data(this_plot)
struct surface_points *this_plot;
{
int i, j, k;
double x, y, z, w, dx, dy, xmin, xmax, ymin, ymax;
struct iso_curve *old_iso_crvs = this_plot->iso_crvs;
struct iso_curve *icrv, *oicrv, *oicrvs;
#ifdef THIN_PLATE_SPLINES_GRID
double *b, **K, *xx, *yy, *zz, d;
int *indx, numpoints;
#endif
/* Compute XY bounding box on the original data. */
xmin = xmax = old_iso_crvs->points[0].x;
ymin = ymax = old_iso_crvs->points[0].y;
for (icrv = old_iso_crvs; icrv != NULL; icrv = icrv->next) {
struct coordinate GPHUGE *points = icrv->points;
for (i = 0; i < icrv->p_count; i++, points++) {
if (xmin > points->x)
xmin = points->x;
if (xmax < points->x)
xmax = points->x;
if (ymin > points->y)
ymin = points->y;
if (ymax < points->y)
ymax = points->y;
}
}
dx = (xmax - xmin) / (dgrid3d_col_fineness - 1);
dy = (ymax - ymin) / (dgrid3d_row_fineness - 1);
/* Create the new grid structure, and compute the low pass filtering from
* non grid to grid structure.
*/
this_plot->iso_crvs = NULL;
this_plot->num_iso_read = dgrid3d_col_fineness;
this_plot->has_grid_topology = TRUE;
#ifdef THIN_PLATE_SPLINES_GRID
numpoints = 0;
for (oicrv = old_iso_crvs; oicrv != NULL; oicrv = oicrv->next) {
numpoints += oicrv->p_count;
}
xx = (double *) gp_alloc(sizeof(double) * (numpoints + 3) * (numpoints + 8),
"thin plate splines in dgrid3d");
/* the memory needed is not really (n+3)*(n+8) for now,
but might be if I take into account errors ... */
K = (double **) gp_alloc(sizeof(double *) * (numpoints + 3),
"matrix : thin plate splines 2d");
yy = xx + numpoints;
zz = yy + numpoints;
b = zz + numpoints;
i = 0;
for (oicrv = old_iso_crvs; oicrv != NULL; oicrv = oicrv->next) {
struct coordinate GPHUGE *opoints = oicrv->points;
for (k = 0; k < oicrv->p_count; k++, opoints++) {
xx[i] = opoints->x;
yy[i] = opoints->y;
zz[i] = opoints->z;
i++;
}
}
for (i = 0; i < numpoints + 3; i++) {
K[i] = b + (numpoints + 3) * (i + 1);
}
for (i = 0; i < numpoints; i++) {
for (j = i + 1; j < numpoints; j++) {
double dx = xx[i] - xx[j], dy = yy[i] - yy[j];
K[i][j] = K[j][i] = -splines_kernel(sqrt(dx * dx + dy * dy));
}
K[i][i] = 0.0; /* here will come the weights for errors */
b[i] = zz[i];
}
for (i = 0; i < numpoints; i++) {
K[i][numpoints] = K[numpoints][i] = 1.0;
K[i][numpoints + 1] = K[numpoints + 1][i] = xx[i];
K[i][numpoints + 2] = K[numpoints + 2][i] = yy[i];
}
b[numpoints] = 0.0;
b[numpoints + 1] = 0.0;
b[numpoints + 2] = 0.0;
K[numpoints][numpoints] = 0.0;
K[numpoints][numpoints + 1] = 0.0;
K[numpoints][numpoints + 2] = 0.0;
K[numpoints + 1][numpoints] = 0.0;
K[numpoints + 1][numpoints + 1] = 0.0;
K[numpoints + 1][numpoints + 2] = 0.0;
K[numpoints + 2][numpoints] = 0.0;
K[numpoints + 2][numpoints + 1] = 0.0;
K[numpoints + 2][numpoints + 2] = 0.0;
indx = (int *) gp_alloc(sizeof(int) * (numpoints + 3), "indexes lu");
/* actually, K is *not* positive definite, but
has only non zero real eigenvalues ->
we can use an lu_decomp safely */
lu_decomp(K, numpoints + 3, indx, &d);
lu_backsubst(K, numpoints + 3, indx, b);
#endif /* THIN_PLATE_SPLINES_GRID */
for (i = 0, x = xmin; i < dgrid3d_col_fineness; i++, x += dx) {
struct coordinate GPHUGE *points;
icrv = iso_alloc(dgrid3d_row_fineness + 1);
icrv->p_count = dgrid3d_row_fineness;
icrv->next = this_plot->iso_crvs;
this_plot->iso_crvs = icrv;
points = icrv->points;
for (j = 0, y = ymin; j < dgrid3d_row_fineness; j++, y += dy, points++) {
z = w = 0.0;
#ifndef BUGGY_DGRID_RANGING /* HBB 981209 */
/* as soon as ->type is changed to UNDEFINED, break out of
* two inner loops! */
points->type = INRANGE;
#endif
#ifdef THIN_PLATE_SPLINES_GRID
z = b[numpoints];
for (k = 0; k < numpoints; k++) {
double dx = xx[k] - x, dy = yy[k] - y;
z = z - b[k] * splines_kernel(sqrt(dx * dx + dy * dy));
}
z = z + b[numpoints + 1] * x + b[numpoints + 2] * y;
#else
for (oicrv = old_iso_crvs; oicrv != NULL; oicrv = oicrv->next) {
struct coordinate GPHUGE *opoints = oicrv->points;
for (k = 0; k < oicrv->p_count; k++, opoints++) {
double dist, dist_x = fabs(opoints->x - x), dist_y = fabs(opoints->y - y);
switch (dgrid3d_norm_value) {
case 1:
dist = dist_x + dist_y;
break;
case 2:
dist = dist_x * dist_x + dist_y * dist_y;
break;
case 4:
dist = dist_x * dist_x + dist_y * dist_y;
dist *= dist;
break;
case 8:
dist = dist_x * dist_x + dist_y * dist_y;
dist *= dist;
dist *= dist;
break;
case 16:
dist = dist_x * dist_x + dist_y * dist_y;
dist *= dist;
dist *= dist;
dist *= dist;
break;
default:
dist = pow(dist_x, (double) dgrid3d_norm_value) +
pow(dist_y, (double) dgrid3d_norm_value);
break;
}
/* The weight of this point is inverse proportional
* to the distance.
*/
if (dist == 0.0) {
#ifndef BUGGY_DGRID_RANGING
/* HBB 981209: revised flagging as undefined */
/* Supporting all those infinities on various
* platforms becomes tiresome, to say the least :-(
* Let's just return the first z where this happens,
* unchanged, and be done with this, period. */
points->type = UNDEFINED;
z = opoints->z;
w = 1.0;
break; /* out of for (k...) loop */
#else
#if !defined(AMIGA_SC_6_1) && !defined(__PUREC__)
dist = VERYLARGE;
#else /* !AMIGA_SC_6_1 && !__PUREC__ */
/* Multiplying VERYLARGE by opoints->z below
* might yield Inf (i.e. a number that can't
* be represented on the machine). This will
* result in points->z being set to NaN. It's
* better to have a pretty large number that is
* also on the safe side... The numbers that are
* read by gnuplot are float values anyway, so
* they can't be bigger than FLT_MAX. So setting
* dist to FLT_MAX^2 will make dist pretty large
* with respect to any value that has been read. */
dist = ((double) FLT_MAX) * ((double) FLT_MAX);
#endif /* !AMIGA_SC_6_1 && !__PUREC__ */
#endif /* BUGGY_DGRID_RANGING */
} else
dist = 1.0 / dist;
z += opoints->z * dist;
w += dist;
}
#ifndef BUGGY_DGRID_RANGING
if (points->type != INRANGE)
break; /* out of the second-inner loop as well ... */
#endif
}
#endif /* THIN_PLATE_SPLINES_GRID */
#ifndef BUGGY_DGRID_RANGING
/* Now that we've escaped the loops safely, we know that we
* do have a good value in z and w, so we can proceed just as
* if nothing had happened at all. Nice, isn't it? */
points->type = INRANGE;
STORE_WITH_LOG_AND_UPDATE_RANGE(points->x, x, points->type, x_axis, NOOP, continue);
STORE_WITH_LOG_AND_UPDATE_RANGE(points->y, y, points->type, y_axis, NOOP, continue);
#ifndef THIN_PLATE_SPLINES_GRID
STORE_WITH_LOG_AND_UPDATE_RANGE(points->z, z / w, points->type, z_axis, NOOP, continue);
#else
STORE_WITH_LOG_AND_UPDATE_RANGE(points->z, z, points->type, z_axis, NOOP, continue);
#endif
#else
/* HBB 981026: original, short version of this code */
points->x = x;
points->y = y;
points->z = z / w;
points->type = INRANGE;
#endif
}
}
#ifdef THIN_PLATE_SPLINES_GRID
free(K);
free(xx);
free(indx);
#endif
/* Delete the old non grid data. */
for (oicrvs = old_iso_crvs; oicrvs != NULL;) {
oicrv = oicrvs;
oicrvs = oicrvs->next;
iso_free(oicrv);
}
}
static void
get_3ddata(this_plot)
struct surface_points *this_plot;
/* this_plot->token is end of datafile spec, before title etc
* will be moved passed title etc after we return
*/
{
int xdatum = 0;
int ydatum = 0;
int i, j;
double v[3];
int pt_in_iso_crv = 0;
struct iso_curve *this_iso;
if (mapping3d == MAP3D_CARTESIAN) {
if (df_no_use_specs == 2)
int_error(this_plot->token, "Need 1 or 3 columns for cartesian data");
} else {
if (df_no_use_specs == 1)
int_error(this_plot->token, "Need 2 or 3 columns for polar data");
}
this_plot->num_iso_read = 0;
this_plot->has_grid_topology = TRUE;
/* we ought to keep old memory - most likely case
* is a replot, so it will probably exactly fit into
* memory already allocated ?
*/
if (this_plot->iso_crvs != NULL) {
struct iso_curve *icrv, *icrvs = this_plot->iso_crvs;
while (icrvs) {
icrv = icrvs;
icrvs = icrvs->next;
iso_free(icrv);
}
this_plot->iso_crvs = NULL;
}
/* data file is already open */
if (df_matrix)
xdatum = df_3dmatrix(this_plot);
else {
/*{{{ read surface from text file */
struct iso_curve *local_this_iso = iso_alloc(samples_1);
struct coordinate GPHUGE *cp;
double x, y, z;
while ((j = df_readline(v, 3)) != DF_EOF) {
if (j == DF_SECOND_BLANK)
break; /* two blank lines */
if (j == DF_FIRST_BLANK) {
/* one blank line */
if (pt_in_iso_crv == 0) {
if (xdatum == 0)
continue;
pt_in_iso_crv = xdatum;
}
if (xdatum > 0) {
local_this_iso->p_count = xdatum;
local_this_iso->next = this_plot->iso_crvs;
this_plot->iso_crvs = local_this_iso;
this_plot->num_iso_read++;
if (xdatum != pt_in_iso_crv)
this_plot->has_grid_topology = FALSE;
local_this_iso = iso_alloc(pt_in_iso_crv);
xdatum = 0;
ydatum++;
}
continue;
}
/* its a data point or undefined */
if (xdatum >= local_this_iso->p_max) {
/*
* overflow about to occur. Extend size of points[] array. We
* either double the size, or add 1000 points, whichever is a
* smaller increment. Note i = p_max.
*/
iso_extend(local_this_iso,
xdatum + (xdatum < 1000 ? xdatum : 1000));
}
cp = local_this_iso->points + xdatum;
if (j == DF_UNDEFINED) {
cp->type = UNDEFINED;
continue;
}
cp->type = INRANGE; /* unless we find out different */
switch (mapping3d) {
case MAP3D_CARTESIAN:
switch (j) {
case 1:
x = xdatum;
y = ydatum;
z = v[0];
break;
case 3:
x = v[0];
y = v[1];
z = v[2];
break;
default:
{
int_error(this_plot->token,
"Need 1 or 3 columns - line %d",
df_line_number);
return; /* avoid gcc -Wuninitialised for x,y,z */
}
}
break;
case MAP3D_SPHERICAL:
if (j < 2)
int_error(this_plot->token, "Need 2 or 3 columns");
if (j < 3)
v[2] = 1; /* default radius */
/* Convert to radians. */
v[0] *= ang2rad;
v[1] *= ang2rad;
x = v[2] * cos(v[0]) * cos(v[1]);
y = v[2] * sin(v[0]) * cos(v[1]);
z = v[2] * sin(v[1]);
break;
case MAP3D_CYLINDRICAL:
if (j < 2)
int_error(this_plot->token, "Need 2 or 3 columns");
if (j < 3)
v[2] = 1; /* default radius */
/* Convert to radians. */
v[0] *= ang2rad;
x = v[2] * cos(v[0]);
y = v[2] * sin(v[0]);
z = v[1];
break;
default:
int_error(NO_CARET, "Internal error : Unknown mapping type");
return;
}
/* adjust for logscales. Set min/max and point types.
* store in cp
*/
cp->type = INRANGE;
/* cannot use continue, as macro is wrapped in a loop.
* I regard this as correct goto use
*/
STORE_WITH_LOG_AND_UPDATE_RANGE(cp->x, x, cp->type, x_axis, NOOP, goto come_here_if_undefined);
STORE_WITH_LOG_AND_UPDATE_RANGE(cp->y, y, cp->type, y_axis, NOOP, goto come_here_if_undefined);
STORE_WITH_LOG_AND_UPDATE_RANGE(cp->z, z, cp->type, z_axis, NOOP, goto come_here_if_undefined);
/* some may complain, but I regard this as the correct use
* of goto
*/
come_here_if_undefined:
++xdatum;
} /* end of whileloop - end of surface */
if (xdatum > 0) {
this_plot->num_iso_read++; /* Update last iso. */
local_this_iso->p_count = xdatum;
local_this_iso->next = this_plot->iso_crvs;
this_plot->iso_crvs = local_this_iso;
if (xdatum != pt_in_iso_crv)
this_plot->has_grid_topology = FALSE;
} else {
iso_free(local_this_iso); /* Free last allocation. */
}
if (dgrid3d && this_plot->num_iso_read > 0)
grid_nongrid_data(this_plot);
/*}}} */
}
if (this_plot->num_iso_read has_grid_topology = FALSE;
if (this_plot->has_grid_topology && !hidden3d) {
struct iso_curve *new_icrvs = NULL;
int num_new_iso = this_plot->iso_crvs->p_count, len_new_iso = this_plot->num_iso_read;
/* Now we need to set the other direction (pseudo) isolines. */
for (i = 0; i < num_new_iso; i++) {
struct iso_curve *new_icrv = iso_alloc(len_new_iso);
new_icrv->p_count = len_new_iso;
for (j = 0, this_iso = this_plot->iso_crvs;
this_iso != NULL;
j++, this_iso = this_iso->next) {
/* copy whole point struct to get type too.
* wasteful for windows, with padding */
/* more efficient would be extra pointer to same struct */
new_icrv->points[j] = this_iso->points[i];
}
new_icrv->next = new_icrvs;
new_icrvs = new_icrv;
}
/* Append the new iso curves after the read ones. */
for (this_iso = this_plot->iso_crvs;
this_iso->next != NULL;
this_iso = this_iso->next);
this_iso->next = new_icrvs;
}
}
static void
print_3dtable(pcount)
int pcount;
{
register struct surface_points *this_plot;
int i, curve, surface;
struct iso_curve *icrvs;
struct coordinate GPHUGE *points;
char *table_format = NULL;
char *pcat;
table_format = gp_alloc(strlen(axis_array[FIRST_X_AXIS].formatstring)
+strlen(axis_array[FIRST_Y_AXIS].formatstring)
+strlen(axis_array[FIRST_Z_AXIS].formatstring)
+6,
"table format");
strcpy(table_format, axis_array[FIRST_X_AXIS].formatstring);
strcat(table_format, " ");
strcat(table_format, axis_array[FIRST_Y_AXIS].formatstring);
strcat(table_format, " ");
strcat(table_format, axis_array[FIRST_Z_AXIS].formatstring);
pcat = &table_format[strlen(table_format)];
for (surface = 0, this_plot = first_3dplot; surface < pcount;
this_plot = this_plot->next_sp, surface++) {
fprintf(gpoutfile, "\n#Surface %d of %d surfaces\n", surface, pcount);
icrvs = this_plot->iso_crvs;
curve = 0;
if (draw_surface) {
strcpy(pcat," %c\n");
/* only the curves in one direction */
while (icrvs && curve < this_plot->num_iso_read) {
fprintf(gpoutfile, "\n#IsoCurve %d, %d points\n#x y z type\n",
curve, icrvs->p_count);
for (i = 0, points = icrvs->points; i < icrvs->p_count; i++) {
fprintf(gpoutfile, table_format,
points[i].x,
points[i].y,
points[i].z,
points[i].type == INRANGE ? 'i'
: points[i].type == OUTRANGE ? 'o'
: 'u');
}
icrvs = icrvs->next;
curve++;
}
putc('\n', gpoutfile);
}
if (draw_contour) {
int number = 0;
struct gnuplot_contours *c = this_plot->contours;
strcpy(pcat,"\n");
while (c) {
int count = c->num_pts;
struct coordinate GPHUGE *p = c->coords;
if (c->isNewLevel)
/* dont display count - contour split across chunks */
/* put # in case user wants to use it for a plot */
/* double blank line to allow plot ... index ... */
fprintf(gpoutfile, "\n# Contour %d, label: %s\n", number++, c->label);
for (; --count >= 0; ++p)
fprintf(gpoutfile, table_format, p->x, p->y, p->z);
/* blank line between segments of same contour */
putc('\n', gpoutfile);
c = c->next;
}
}
}
fflush(gpoutfile);
free(table_format);
}
/* HBB 20000501: code isolated from eval_3dplots(), where practically
* identical code occured twice, for direct and crossing isolines,
* respectively. The latter only are done for in non-hidden3d
* mode. */
static void
calculate_set_of_isolines(value_axis, cross, this_iso,
iso_axis, iso_min, iso_step, num_iso_to_use,
sam_axis, sam_min, sam_step, num_sam_to_use
)
AXIS_INDEX iso_axis, sam_axis, value_axis;
struct iso_curve **this_iso;
TBOOLEAN cross;
double iso_min, iso_step, sam_min, sam_step;
int num_iso_to_use, num_sam_to_use;
{
int i, j;
struct coordinate GPHUGE *points = (*this_iso)->points;
for (j = 0; j < num_iso_to_use; j++) {
double iso = iso_min + j * iso_step;
/* HBB 20000501: with the new code, it should
* be safe to rely on the actual 'v' axis not
* to be improperly logscaled... */
(void) Gcomplex(&plot_func.dummy_values[cross ? 0 : 1],
AXIS_DE_LOG_VALUE(iso_axis, iso), 0.0);
for (i = 0; i < num_sam_to_use; i++) {
double sam = sam_min + i * sam_step;
struct value a;
double temp;
(void) Gcomplex(&plot_func.dummy_values[cross ? 1 : 0],
AXIS_DE_LOG_VALUE(sam_axis, sam), 0.0);
if (cross) {
points[i].x = iso;
points[i].y = sam;
} else {
points[i].x = sam;
points[i].y = iso;
}
evaluate_at(plot_func.at, &a);
if (undefined || (fabs(imag(&a)) > zero)) {
points[i].type = UNDEFINED;
continue;
}
temp = real(&a);
points[i].type = INRANGE;
STORE_WITH_LOG_AND_UPDATE_RANGE(points[i].z, temp, points[i].type,
value_axis, NOOP, NOOP);
}
(*this_iso)->p_count = num_sam_to_use;
*this_iso = (*this_iso)->next;
points = (*this_iso) ? (*this_iso)->points : NULL;
}
}
/*
* This parses the splot command after any range specifications. To support
* autoscaling on the x/z axis, we want any data files to define the x/y
* range, then to plot any functions using that range. We thus parse the
* input twice, once to pick up the data files, and again to pick up the
* functions. Definitions are processed twice, but that won't hurt.
* div - okay, it doesn't hurt, but every time an option as added for
* datafiles, code to parse it has to be added here. Change so that
* we store starting-token in the plot structure.
*/
static void
eval_3dplots()
{
int i;
struct surface_points **tp_3d_ptr;
int start_token, end_token;
int begin_token;
TBOOLEAN some_data_files = FALSE, some_functions = FALSE;
int plot_num, line_num, point_num;
/* part number of parametric function triplet: 0 = z, 1 = y, 2 = x */
int crnt_param = 0;
char *xtitle;
char *ytitle;
/* Reset first_3dplot. This is usually done at the end of this function.
* If there is an error within this function, the memory is left allocated,
* since we cannot call sp_free if the list is incomplete
*/
if (first_3dplot && plot3d_num>0)
sp_free(first_3dplot);
plot3d_num=0;
first_3dplot = NULL;
x_axis = FIRST_X_AXIS;
y_axis = FIRST_Y_AXIS;
z_axis = FIRST_Z_AXIS;
tp_3d_ptr = &(first_3dplot);
plot_num = 0;
line_num = 0; /* default line type */
point_num = 0; /* default point type */
xtitle = NULL;
ytitle = NULL;
begin_token = c_token;
/*** First Pass: Read through data files ***/
/*
* This pass serves to set the x/yranges and to parse the command, as
* well as filling in every thing except the function data. That is done
* after the x/yrange is defined.
*/
while (TRUE) {
if (END_OF_COMMAND)
int_error(c_token, "function to plt3d expected");
start_token = c_token;
if (is_definition(c_token)) {
define();
} else {
int specs;
struct surface_points *this_plot;
if (isstring(c_token)) { /* data file to plot */
/*{{{ data file */
if (parametric && crnt_param != 0)
int_error(c_token, "previous parametric function not fully specified");
if (!some_data_files) {
if (axis_array[FIRST_X_AXIS].autoscale & 1) {
axis_array[FIRST_X_AXIS].min = VERYLARGE;
}
if (axis_array[FIRST_X_AXIS].autoscale & 2) {
axis_array[FIRST_X_AXIS].max = -VERYLARGE;
}
if (axis_array[FIRST_Y_AXIS].autoscale & 1) {
axis_array[FIRST_Y_AXIS].min = VERYLARGE;
}
if (axis_array[FIRST_Y_AXIS].autoscale & 2) {
axis_array[FIRST_Y_AXIS].max = -VERYLARGE;
}
some_data_files = TRUE;
}
if (*tp_3d_ptr)
this_plot = *tp_3d_ptr;
else { /* no memory malloc()'d there yet */
/* Allocate enough isosamples and samples */
this_plot = sp_alloc(0, 0, 0, 0);
*tp_3d_ptr = this_plot;
}
this_plot->plot_type = DATA3D;
this_plot->plot_style = data_style;
specs = df_open(3);
/* parses all datafile-specific modifiers */
/* we will load the data after parsing title,with,... */
/* for capture to key */
this_plot->token = end_token = c_token - 1;
/* this_plot->token is temporary, for errors in get_3ddata() */
if (specs < 3) {
if (axis_array[FIRST_X_AXIS].is_timedata) {
int_error(c_token, "Need full using spec for x time data");
}
if (axis_array[FIRST_Y_AXIS].is_timedata) {
int_error(c_token, "Need full using spec for y time data");
}
/* df_axis[0] = FIRST_Z_AXIS; */
} /* else */ { /* HBB 20000725: testestest */
df_axis[0] = FIRST_X_AXIS;
df_axis[1] = FIRST_Y_AXIS;
df_axis[2] = FIRST_Z_AXIS;
}
/*}}} */
} else { /* function to plot */
/*{{{ function */
++plot_num;
if (parametric) {
/* Rotate between x/y/z axes */
/* +2 same as -1, but beats -ve problem */
crnt_param = (crnt_param + 2) % 3;
}
if (*tp_3d_ptr) {
this_plot = *tp_3d_ptr;
if (!hidden3d)
sp_replace(this_plot, samples_1, iso_samples_1,
samples_2, iso_samples_2);
else
sp_replace(this_plot, iso_samples_1, 0,
0, iso_samples_2);
} else { /* no memory malloc()'d there yet */
/* Allocate enough isosamples and samples */
if (!hidden3d)
this_plot = sp_alloc(samples_1, iso_samples_1,
samples_2, iso_samples_2);
else
this_plot = sp_alloc(iso_samples_1, 0,
0, iso_samples_2);
*tp_3d_ptr = this_plot;
}
this_plot->plot_type = FUNC3D;
this_plot->has_grid_topology = TRUE;
this_plot->plot_style = func_style;
this_plot->num_iso_read = iso_samples_2;
dummy_func = &plot_func;
plot_func.at = temp_at();
dummy_func = NULL;
/* ignore it for now */
some_functions = TRUE;
end_token = c_token - 1;
/*}}} */
} /* end of IS THIS A FILE OR A FUNC block */
/*{{{ title */
if (this_plot->title) {
free(this_plot->title);
this_plot->title = NULL;
}
if (almost_equals(c_token, "t$itle")) {
this_plot->title_no_enhanced = 0; /* can be enhanced */
if (parametric) {
if (crnt_param != 0)
int_error(c_token, "\"title\" allowed only after parametric function fully specified");
else {
if (xtitle != NULL)
xtitle[0] = NUL; /* Remove default title . */
if (ytitle != NULL)
ytitle[0] = NUL; /* Remove default title . */
}
}
if (isstring(++c_token))
m_quote_capture(&(this_plot->title), c_token, c_token);
else
int_error(c_token, "expecting \"title\" for plot");
/* end of new method */
++c_token;
} else if (almost_equals(c_token, "not$itle")) {
if (xtitle != NULL)
xtitle[0] = '\0';
if (ytitle != NULL)
ytitle[0] = '\0';
/* this_plot->title = NULL; */
++c_token;
} else {
this_plot->title_no_enhanced = 1; /* filename or function cannot be enhanced */
m_capture(&(this_plot->title), start_token, end_token);
if (crnt_param == 2)
xtitle = this_plot->title;
else if (crnt_param == 1)
ytitle = this_plot->title;
}
/*}}} */
/*{{{ line types, widths, ... */
this_plot->lp_properties.l_type = line_num;
this_plot->lp_properties.p_type = point_num;
if (almost_equals(c_token, "w$ith")) {
this_plot->plot_style = get_style();
}
/* pick up line/point specs
* - point spec allowed if style uses points, ie style&2 != 0
* - keywords are optional
*/
lp_parse(&this_plot->lp_properties, 1, this_plot->plot_style & 2,
line_num, point_num);
/* allow old-style syntax too - ignore case lt 3 4 for example */
if (!equals(c_token, ",") && !END_OF_COMMAND) {
struct value t;
this_plot->lp_properties.l_type =
this_plot->lp_properties.p_type = (int) real(const_express(&t)) - 1;
if (!equals(c_token, ",") && !END_OF_COMMAND)
this_plot->lp_properties.p_type = (int) real(const_express(&t)) - 1;
}
if (this_plot->plot_style & 2) /* lines, linesp, ... */
if (crnt_param == 0)
point_num +=
1 + (draw_contour != 0)
+ (hidden3d != 0);
if (crnt_param == 0)
line_num += 1 + (draw_contour != 0)
+ (hidden3d != 0);
/*}}} */
/* now get the data... having to think hard here...
* first time through, we fill in this_plot. For second
* surface in file, we have to allocate another surface
* struct. BUT we may allocate this store only to
* find that it is merely some blank lines at end of file
* tp_3d_ptr is still pointing at next field of prev. plot,
* before : prev_or_first -> this_plot -> possible_preallocated_store
* tp_3d_ptr--^
* after : prev_or_first -> first -> second -> last -> possibly_more_store
* tp_3d_ptr ----^
* if file is empty, tp_3d_ptr is not moved. this_plot continues
* to point at allocated storage, but that will be reused later
*/
assert(this_plot == *tp_3d_ptr);
if (this_plot->plot_type == DATA3D) {
/*{{{ read data */
/* remember settings for second surface in file */
struct lp_style_type *these_props = &(this_plot->lp_properties);
enum PLOT_STYLE this_style = this_plot->plot_style;
int this_token = this_plot->token;
while (!df_eof) {
this_plot = *tp_3d_ptr;
assert(this_plot != NULL);
/* dont move tp_3d_ptr until we are sure we
* have read a surface
*/
/* used by get_3ddata() */
this_plot->token = this_token;
get_3ddata(this_plot);
/* for second pass */
this_plot->token = c_token;
if (this_plot->num_iso_read == 0)
/* probably df_eof, in which case we
* will leave loop. if not eof, then
* how come we got no surface ? - retry
* in neither case do we update tp_3d_ptr
*/
continue;
/* okay, we have read a surface */
++plot_num;
tp_3d_ptr = &(this_plot->next_sp);
if (df_eof)
break;
/* there might be another surface so allocate
* and prepare another surface structure
* This does no harm if in fact there are
* no more surfaces to read
*/
if ((this_plot = *tp_3d_ptr) != NULL) {
if (this_plot->title) {
free(this_plot->title);
this_plot->title = NULL;
}
} else {
/* Allocate enough isosamples and samples */
this_plot = *tp_3d_ptr = sp_alloc(0, 0, 0, 0);
}
this_plot->plot_type = DATA3D;
this_plot->plot_style = this_style;
/* Struct copy */
this_plot->lp_properties = *these_props;
}
df_close();
/*}}} */
} else { /* not a data file */
tp_3d_ptr = &(this_plot->next_sp);
this_plot->token = c_token; /* store for second pass */
}
} /* !is_definition() : end of scope of this_plot */
if (equals(c_token, ","))
c_token++;
else
break;
} /* while(TRUE), ie first pass */
if (parametric && crnt_param != 0)
int_error(NO_CARET, "parametric function not fully specified");
/*** Second Pass: Evaluate the functions ***/
/*
* Everything is defined now, except the function data. We expect no
* syntax errors, etc, since the above parsed it all. This makes the code
* below simpler. If axis_array[FIRST_Y_AXIS].autoscale, the yrange may still change.
* - eh ? - z can still change. x/y/z can change if we are parametric ??
*/
if (some_functions) {
/* I've changed the controlled variable in fn plots to u_min etc since
* it's easier for me to think parametric - 'normal' plot is after all
* a special case. I was confused about x_min being both minimum of
* x values found, and starting value for fn plots.
*/
register double u_min, u_max, u_step, v_min, v_max, v_step;
double u_isostep, v_isostep;
AXIS_INDEX u_axis, v_axis;
struct surface_points *this_plot;
/* Make these point out the right 'u' and 'v' axis. In
* non-parametric mode, x is used as u, and y as v */
u_axis = parametric ? U_AXIS : FIRST_X_AXIS;
v_axis = parametric ? V_AXIS : FIRST_Y_AXIS;
if (!parametric) {
/*{{{ check ranges */
/* give error if xrange badly set from missing datafile error
* parametric fn can still set ranges
* if there are no fns, we'll report it later as 'nothing to plot'
*/
/* check that xmin -> xmax is not too small */
axis_checked_extend_empty_range(FIRST_X_AXIS, "x range is invalid");
axis_checked_extend_empty_range(FIRST_Y_AXIS, "y range is invalid");
/*}}} */
}
if (parametric && !some_data_files) {
/*{{{ set ranges */
/* parametric fn can still change x/y range */
if (axis_array[FIRST_X_AXIS].autoscale & 1)
axis_array[FIRST_X_AXIS].min = VERYLARGE;
if (axis_array[FIRST_X_AXIS].autoscale & 2)
axis_array[FIRST_X_AXIS].max = -VERYLARGE;
if (axis_array[FIRST_Y_AXIS].autoscale & 1)
axis_array[FIRST_Y_AXIS].min = VERYLARGE;
if (axis_array[FIRST_Y_AXIS].autoscale & 2)
axis_array[FIRST_Y_AXIS].max = -VERYLARGE;
/*}}} */
}
/*{{{ figure ranges, taking logs etc into account */
u_min = axis_log_value_checked(u_axis, axis_array[u_axis].min, "x range");
u_max = axis_log_value_checked(u_axis, axis_array[u_axis].max, "x range");
v_min = axis_log_value_checked(v_axis, axis_array[v_axis].min, "y range");
v_max = axis_log_value_checked(v_axis, axis_array[v_axis].max, "y range");
/*}}} */
if (samples_1 < 2 || samples_2 < 2 || iso_samples_1 < 2 ||
iso_samples_2 < 2) {
int_error(NO_CARET, "samples or iso_samples < 2. Must be at least 2.");
}
/* start over */
this_plot = first_3dplot;
c_token = begin_token;
/* why do attributes of this_plot matter ? */
/* FIXME HBB 20000501: I think they don't, actually. I'm
* taking out references to has_grid_topology in this part of
* the code. It only deals with function, which always is
* gridded */
if (hidden3d) {
u_step = (u_max - u_min) / (iso_samples_1 - 1);
v_step = (v_max - v_min) / (iso_samples_2 - 1);
} else {
u_step = (u_max - u_min) / (samples_1 - 1);
v_step = (v_max - v_min) / (samples_2 - 1);
}
u_isostep = (u_max - u_min) / (iso_samples_1 - 1);
v_isostep = (v_max - v_min) / (iso_samples_2 - 1);
/* Read through functions */
while (TRUE) {
if (is_definition(c_token)) {
define();
} else {
if (!isstring(c_token)) { /* func to plot */
/*{{{ evaluate function */
struct iso_curve *this_iso = this_plot->iso_crvs;
int num_sam_to_use, num_iso_to_use;
/* crnt_param is used as the axis number. As the
* axis array indices are ordered z, y, x, we have
* to count *backwards*, starting starting at 2,
* to properly store away contents to x, y and
* z. The following little gimmick does that. */
if (parametric)
crnt_param = (crnt_param + 2) % 3;
dummy_func = &plot_func;
plot_func.at = temp_at(); /* reparse function */
dummy_func = NULL;
num_iso_to_use = iso_samples_2;
num_sam_to_use = hidden3d ? iso_samples_1 : samples_1;
calculate_set_of_isolines(crnt_param, FALSE, &this_iso,
v_axis, v_min, v_isostep,
num_iso_to_use,
u_axis, u_min, u_step,
num_sam_to_use);
if (!hidden3d) {
num_iso_to_use = iso_samples_1;
num_sam_to_use = samples_2;
calculate_set_of_isolines(crnt_param, TRUE, &this_iso,
u_axis, u_min, u_isostep,
num_iso_to_use,
v_axis, v_min, v_step,
num_sam_to_use);
}
/*}}} */
} /* end of ITS A FUNCTION TO PLOT */
/* we saved it from first pass */
c_token = this_plot->token;
/* one data file can make several plots */
do
this_plot = this_plot->next_sp;
while (this_plot && this_plot->token == c_token);
} /* !is_definition */
if (equals(c_token, ","))
c_token++;
else
break;
} /* while(TRUE) */
if (parametric) {
/* Now actually fix the plot triplets to be single plots. */
parametric_3dfixup(first_3dplot, &plot_num);
}
} /* some functions */
/* if first_3dplot is NULL, we have no functions or data at all.
* This can happen, if you type "splot x=5", since x=5 is a
* variable assignment
*/
if (plot_num == 0 || first_3dplot == NULL) {
int_error(c_token, "no functions or data to plot");
}
axis_checked_extend_empty_range(FIRST_X_AXIS, "All points x value undefined");
axis_checked_extend_empty_range(FIRST_Y_AXIS, "All points y value undefined");
axis_checked_extend_empty_range(FIRST_Z_AXIS, "All points z value undefined");
axis_revert_and_unlog_range(FIRST_X_AXIS);
axis_revert_and_unlog_range(FIRST_Y_AXIS);
axis_revert_and_unlog_range(FIRST_Z_AXIS);
setup_tics(FIRST_X_AXIS, 20);
setup_tics(FIRST_Y_AXIS, 20);
setup_tics(FIRST_Z_AXIS, 20);
AXIS_WRITEBACK(FIRST_X_AXIS);
if (plot_num == 0 || first_3dplot == NULL) {
int_error(c_token, "no functions or data to plot");
}
/* Creates contours if contours are to be plotted as well. */
if (draw_contour) {
struct surface_points *this_plot;
for (this_plot = first_3dplot, i = 0;
i < plot_num;
this_plot = this_plot->next_sp, i++) {
if (this_plot->contours) {
struct gnuplot_contours *cntrs = this_plot->contours;
while (cntrs) {
struct gnuplot_contours *cntr = cntrs;
cntrs = cntrs->next;
free(cntr->coords);
free(cntr);
}
}
/* Make sure this one can be contoured. */
if (!this_plot->has_grid_topology) {
this_plot->contours = NULL;
fputs("Notice: cannot contour non grid data!\n", stderr);
/* changed from int_error by recommendation of
* rkc@xn.ll.mit.edu
*/
} else if (this_plot->plot_type == DATA3D) {
this_plot->contours = contour(this_plot->num_iso_read,
this_plot->iso_crvs);
} else {
this_plot->contours = contour(iso_samples_2,
this_plot->iso_crvs);
}
}
} /* draw_contour */
/* perform the plot */
if (strcmp(term->name, "table") == 0)
print_3dtable(plot_num);
else {
START_LEAK_CHECK(); /* assert no memory leaks here ! */
do_3dplot(first_3dplot, plot_num, 0);
END_LEAK_CHECK();
/* after do_3dplot(), axis_array[] and max_array[].min
* contain the plotting range actually used (rounded
* to tic marks, not only the min/max data values)
* --> save them now for writeback if requested
*/
SAVE_WRITEBACK_ALL_AXES;
}
/* if we get here, all went well, so record the line for replot */
if (plot_token != -1) {
/* note that m_capture also frees the old replot_line */
m_capture(&replot_line, plot_token, c_token - 1);
plot_token = -1;
}
if (first_3dplot)
sp_free(first_3dplot);
first_3dplot = NULL;
/* record that all went well */
plot3d_num=plot_num;
}
static void
parametric_3dfixup(start_plot, plot_num)
struct surface_points *start_plot;
int *plot_num;
/*
* The hardest part of this routine is collapsing the FUNC plot types in the
* list (which are gauranteed to occur in (x,y,z) triplets while preserving
* the non-FUNC type plots intact. This means we have to work our way
* through various lists. Examples (hand checked):
* start_plot:F1->F2->F3->NULL ==> F3->NULL
* start_plot:F1->F2->F3->F4->F5->F6->NULL ==> F3->F6->NULL
* start_plot:F1->F2->F3->D1->D2->F4->F5->F6->D3->NULL ==>
* F3->D1->D2->F6->D3->NULL
*/
{
/*
* I initialized *free_list with NULL, because my compiler warns some lines
* later that it might be uninited. The code however seems to not access that
* line in that case, but if I'm right, my change is OK and if not, this is a
* serious bug in the code.
*
* x and y ranges now fixed in eval_3dplots
*/
struct surface_points *xp, *new_list, *free_list = NULL;
struct surface_points **last_pointer = &new_list;
size_t tlen;
int i, surface;
char *new_title;
/*
* Ok, go through all the plots and move FUNC3D types together. Note:
* this originally was written to look for a NULL next pointer, but
* gnuplot wants to be sticky in grabbing memory and the right number of
* items in the plot list is controlled by the plot_num variable.
*
* Since gnuplot wants to do this sticky business, a free_list of
* surface_points is kept and then tagged onto the end of the plot list
* as this seems more in the spirit of the original memory behavior than
* simply freeing the memory. I'm personally not convinced this sort of
* concern is worth it since the time spent computing points seems to
* dominate any garbage collecting that might be saved here...
*/
new_list = xp = start_plot;
for (surface = 0; surface < *plot_num; surface++) {
if (xp->plot_type == FUNC3D) {
struct surface_points *yp = xp->next_sp;
struct surface_points *zp = yp->next_sp;
/* Here's a FUNC3D parametric function defined as three parts.
* Go through all the points and assign the x's and y's from xp and
* yp to zp. min/max already done
*/
struct iso_curve *xicrvs = xp->iso_crvs;
struct iso_curve *yicrvs = yp->iso_crvs;
struct iso_curve *zicrvs = zp->iso_crvs;
(*plot_num) -= 2;
assert(INRANGE < OUTRANGE && OUTRANGE < UNDEFINED);
while (zicrvs) {
struct coordinate GPHUGE *xpoints = xicrvs->points, GPHUGE * ypoints =
yicrvs->points, GPHUGE * zpoints = zicrvs->points;
for (i = 0; i < zicrvs->p_count; ++i) {
zpoints[i].x = xpoints[i].z;
zpoints[i].y = ypoints[i].z;
if (zpoints[i].type < xpoints[i].type)
zpoints[i].type = xpoints[i].type;
if (zpoints[i].type < ypoints[i].type)
zpoints[i].type = ypoints[i].type;
}
xicrvs = xicrvs->next;
yicrvs = yicrvs->next;
zicrvs = zicrvs->next;
}
/* Ok, fix up the title to include xp and yp plots. */
if (((xp->title && xp->title[0] != '\0') ||
(yp->title && yp->title[0] != '\0')) && zp->title) {
tlen = (xp->title ? strlen(xp->title) : 0) +
(yp->title ? strlen(yp->title) : 0) +
(zp->title ? strlen(zp->title) : 0) + 5;
new_title = gp_alloc(tlen, "string");
new_title[0] = 0;
if (xp->title && xp->title[0] != '\0') {
strcat(new_title, xp->title);
strcat(new_title, ", "); /* + 2 */
}
if (yp->title && yp->title[0] != '\0') {
strcat(new_title, yp->title);
strcat(new_title, ", "); /* + 2 */
}
strcat(new_title, zp->title);
free(zp->title);
zp->title = new_title;
}
/* add xp and yp to head of free list */
assert(xp->next_sp == yp);
yp->next_sp = free_list;
free_list = xp;
/* add zp to tail of new_list */
*last_pointer = zp;
last_pointer = &(zp->next_sp);
xp = zp->next_sp;
} else { /* its a data plot */
assert(*last_pointer == xp); /* think this is true ! */
last_pointer = &(xp->next_sp);
xp = xp->next_sp;
}
}
/* Ok, append free list and write first_plot */
*last_pointer = free_list;
first_3dplot = new_list;
}