--- /dev/null
+#ifndef lint
+static char *RCSid() { return RCSid("$Id: plot3d.c,v 1.131.2.6 2008/03/20 09:39:46 sfeam Exp $"); }
+#endif
+
+/* GNUPLOT - plot3d.c */
+
+/*[
+ * Copyright 1986 - 1993, 1998, 2004 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 "eval.h"
+#include "graph3d.h"
+#include "misc.h"
+#include "parse.h"
+#include "setshow.h"
+#include "term_api.h"
+#include "util.h"
+#include "pm3d.h"
+#ifdef BINARY_DATA_FILE
+#include "plot.h"
+#endif
+
+#ifdef EAM_DATASTRINGS
+#include "plot2d.h" /* Only for store_label() */
+#endif
+
+#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,
+ TBOOLEAN need_palette));
+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));
+#ifdef THIN_PLATE_SPLINES_GRID
+static double splines_kernel __PROTO((double h));
+#endif
+
+/* 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(int num_samp_1, int num_iso_1, int num_samp_2, int num_iso_2)
+{
+ struct lp_style_type default_lp_properties = DEFAULT_LP_STYLE_TYPE;
+
+ struct surface_points *sp = gp_alloc(sizeof(*sp), "surface");
+ memset(sp,0,sizeof(struct surface_points));
+
+ /* Initialize various fields */
+ sp->lp_properties = default_lp_properties;
+ default_arrow_style(&(sp->arrow_properties));
+
+ 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;
+ }
+ }
+
+ 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(
+ struct surface_points *sp,
+ int num_samp_1, int num_iso_1, int num_samp_2, int 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 = 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(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;
+ }
+
+#ifdef EAM_DATASTRINGS
+ if (sp->labels) {
+ free_labels(sp->labels);
+ sp->labels = (struct text_label *)NULL;
+ }
+#endif
+
+ 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;
+
+ /* change view to become map if requested by 'set view map' */
+ if (splot_map == TRUE)
+ splot_map_activate();
+
+ 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);
+ AXIS_INIT3D(COLOR_AXIS, 0, 1);
+
+ 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);
+ if (splot_map == TRUE && !parametric) /* v_axis==FIRST_Y_AXIS */
+ splot_map_deactivate();
+ PARSE_NAMED_RANGE(v_axis, dummy_token1);
+ if (splot_map == TRUE && !parametric) /* v_axis==FIRST_Y_AXIS */
+ splot_map_activate();
+
+ if (parametric) {
+ PARSE_RANGE(FIRST_X_AXIS);
+ if (splot_map == TRUE)
+ splot_map_deactivate();
+ PARSE_RANGE(FIRST_Y_AXIS);
+ if (splot_map == TRUE)
+ splot_map_activate();
+ } /* 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(double h)
+{
+#if 0
+ /* this is normaly not useful ... */
+ h = fabs(h);
+
+ if (h != 0.0) {
+#else
+ if (h > 0) {
+#endif
+ return h * h * log(h);
+ } else {
+ return 0;
+ }
+}
+
+#endif
+
+static void
+grid_nongrid_data(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. */
+ /* FIXME HBB 20010424: Does this make any sense? Shouldn't we just
+ * use whatever the x and y ranges have been found to be, and
+ * that's that? The largest difference this is going to make is if
+ * we plot a datafile that doesn't span the whole x/y range
+ * used. Do we want a dgrid3d over the actual data rectangle, or
+ * over the xrange/yrange area? */
+ 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++) {
+ /* HBB 20010424: avoid crashing for undefined input */
+ if (points->type == UNDEFINED)
+ continue;
+ 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 = gp_alloc(sizeof(xx[0]) * (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 = gp_alloc(sizeof(K[0]) * (numpoints + 3),
+ "matrix : thin plate splines 2d");
+ yy = xx + numpoints;
+ zz = yy + numpoints;
+ b = zz + numpoints;
+
+ /* HBB 20010424: Count actual input points without the UNDEFINED
+ * ones, as we copy them */
+ numpoints = 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++) {
+ /* HBB 20010424: avoid crashing for undefined input */
+ if (opoints->type == UNDEFINED)
+ continue;
+ xx[numpoints] = opoints->x;
+ yy[numpoints] = opoints->y;
+ zz[numpoints] = opoints->z;
+ numpoints++;
+ }
+ }
+
+ 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 = gp_alloc(sizeof(indx[0]) * (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;
+
+ /* as soon as ->type is changed to UNDEFINED, break out of
+ * two inner loops! */
+ points->type = INRANGE;
+
+#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) {
+ /* 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 inner loop */
+ } else {
+ dist = 1.0 / dist;
+ z += opoints->z * dist;
+ w += dist;
+ }
+ }
+ if (points->type != INRANGE)
+ break; /* out of the second-inner loop as well ... */
+ }
+#endif /* THIN_PLATE_SPLINES_GRID */
+
+ /* 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;
+
+ /* HBB 20010424: if log x or log y axis, we don't want to
+ * log() the value again --> just store it, and trust that
+ * it's always inrange */
+ points->x = x;
+ points->y = y;
+
+#ifndef THIN_PLATE_SPLINES_GRID
+ z = z / w;
+#endif
+
+ STORE_WITH_LOG_AND_UPDATE_RANGE(points->z, z, points->type, z_axis, NOOP, continue);
+
+ if (this_plot->pm3d_color_from_column)
+ int_error(NO_CARET, "Gridding of the color column is not implemented");
+ else {
+ COLOR_STORE_WITH_LOG_AND_UPDATE_RANGE(points->CRD_COLOR, z, points->type, COLOR_AXIS, NOOP, continue);
+ }
+ }
+ }
+
+#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);
+ }
+}
+
+/* Get 3D data from file, and store into this_plot data
+ * structure. Takes care of 'set mapping' and 'set dgrid3d'.
+ *
+ * Notice: this_plot->token is end of datafile spec, before title etc
+ * will be moved past title etc after we return */
+static void
+get_3ddata(struct surface_points *this_plot)
+{
+ int xdatum = 0;
+ int ydatum = 0;
+ int j;
+ double v[MAXDATACOLS];
+ int pt_in_iso_crv = 0;
+ struct iso_curve *this_iso;
+
+ if (mapping3d == MAP3D_CARTESIAN) {
+ /* do this check only, if we have PM3D / PM3D-COLUMN not compiled in */
+ if (df_no_use_specs == 2)
+ int_error(this_plot->token, "Need 1 or 3 columns for cartesian data");
+ /* HBB NEW 20060427: if there's only one, explicit using
+ * column, it's z data. df_axis[] has to reflect that, so
+ * df_readline() will expect time/date input. */
+ if (df_no_use_specs == 1)
+ df_axis[0] = FIRST_Z_AXIS;
+ } 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;
+ this_plot->pm3d_color_from_column = FALSE;
+
+ /* 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 */
+
+#ifndef BINARY_DATA_FILE /* NO LONGER REQUIRED FOR GENERAL BINARY OR MATRIX BINARY DATA */
+ if (df_matrix)
+ xdatum = df_3dmatrix(this_plot, NEED_PALETTE(this_plot));
+ else {
+#else
+ if (df_matrix) {
+ this_plot->has_grid_topology = TRUE;
+ }
+
+ {
+#endif
+ /*{{{ read surface from text file */
+ struct iso_curve *local_this_iso = iso_alloc(samples_1);
+ struct coordinate GPHUGE *cp;
+ struct coordinate GPHUGE *cptail = NULL; /* Only for VECTOR plots */
+ double x, y, z;
+ double xtail, ytail, ztail;
+ double color = VERYLARGE;
+ int pm3d_color_from_column = FALSE;
+#define color_from_column(x) pm3d_color_from_column = x
+
+#ifdef EAM_DATASTRINGS
+ if (this_plot->plot_style == LABELPOINTS)
+ expect_string( 4 );
+#endif
+
+ if (this_plot->plot_style == VECTOR) {
+ local_this_iso->next = iso_alloc(samples_1);
+ local_this_iso->next->p_count = 0;
+ }
+
+ while ((j = df_readline(v,MAXDATACOLS)) != DF_EOF) {
+ if (j == DF_SECOND_BLANK)
+ break; /* two blank lines */
+ if (j == DF_FIRST_BLANK) {
+
+#if defined(WITH_IMAGE) && defined(BINARY_DATA_FILE)
+ /* Images are in a sense similar to isocurves.
+ * However, the routine for images is written to
+ * compute the two dimensions of coordinates by
+ * examining the data alone. That way it can be used
+ * in the 2D plots, for which there is no isoline
+ * record. So, toss out isoline information for
+ * images.
+ */
+ if ((this_plot->plot_style == IMAGE)
+ || (this_plot->plot_style == RGBIMAGE))
+ continue;
+#endif
+ if (this_plot->plot_style == VECTOR)
+ continue;
+
+ /* 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. Double the size, and add 1000 points, to
+ * avoid needlessly small steps. */
+ iso_extend(local_this_iso, xdatum + xdatum + 1000);
+ if (this_plot->plot_style == VECTOR) {
+ iso_extend(local_this_iso->next, xdatum + xdatum + 1000);
+ local_this_iso->next->p_count = 0;
+ }
+ }
+ cp = local_this_iso->points + xdatum;
+
+ if (this_plot->plot_style == VECTOR) {
+ if (j < 6) {
+ cp->type = UNDEFINED;
+ continue;
+ }
+ cptail = local_this_iso->next->points + xdatum;
+ }
+
+ if (j == DF_UNDEFINED || j == DF_MISSING) {
+ cp->type = UNDEFINED;
+ goto come_here_if_undefined;
+ }
+ cp->type = INRANGE; /* unless we find out different */
+
+ /* EAM Oct 2004 - Substantially rework this section */
+ /* now that there are many more plot types. */
+
+ x = y = z = 0.0;
+ xtail = ytail = ztail = 0.0;
+
+ /* The x, y, z coordinates depend on the mapping type */
+ switch (mapping3d) {
+
+ case MAP3D_CARTESIAN:
+ if (j == 1) {
+ x = xdatum;
+ y = ydatum;
+ z = v[0];
+ j = 3;
+ break;
+ }
+
+ if (j == 2) {
+ if (PM3DSURFACE != this_plot->plot_style)
+ int_error(this_plot->token,
+ "2 columns only possible with explicit pm3d style (line %d)",
+ df_line_number);
+ x = xdatum;
+ y = ydatum;
+ z = v[0];
+ color_from_column(TRUE);
+ color = v[1];
+ j = 3;
+ break;
+ }
+
+ /* Assume everybody agrees that x,y,z are the first three specs */
+ if (j >= 3) {
+ x = v[0];
+ y = v[1];
+ z = v[2];
+ break;
+ }
+
+ 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 */
+ j = 3;
+ }
+
+ /* 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 */
+ j = 3;
+ }
+
+ /* 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;
+ }
+
+ if (j < df_no_use_specs)
+ int_error(this_plot->token,
+ "Wrong number of columns in input data - line %d",
+ df_line_number);
+
+ /* After the first three columns it gets messy because */
+ /* different plot styles assume different contents in the columns */
+ if (j >= 4) {
+ if (( this_plot->plot_style == POINTSTYLE
+ || this_plot->plot_style == LINESPOINTS)
+ && this_plot->lp_properties.p_size == PTSZ_VARIABLE) {
+ cp->CRD_PTSIZE = v[3];
+ color = z;
+ color_from_column(FALSE);
+ }
+#ifdef EAM_DATASTRINGS
+ else if (this_plot->plot_style == LABELPOINTS) {
+ /* 4th column holds label text rather than color */
+ /* text = df_tokens[3]; */
+ color = z;
+ color_from_column(FALSE);
+ }
+#endif
+ else {
+ color = v[3];
+ color_from_column(TRUE);
+ }
+ }
+
+ if (j >= 5) {
+ if ((this_plot->plot_style == POINTSTYLE
+ || this_plot->plot_style == LINESPOINTS)
+ && this_plot->lp_properties.p_size == PTSZ_VARIABLE) {
+ color = v[4];
+ color_from_column(TRUE);
+ }
+#ifdef EAM_DATASTRINGS
+ if (this_plot->plot_style == LABELPOINTS) {
+ /* take color from an explicitly given 5th column */
+ color = v[4];
+ color_from_column(TRUE);
+ }
+#endif
+ }
+
+ if (j >= 6) {
+ if (this_plot->plot_style == VECTOR) {
+ xtail = x + v[3];
+ ytail = y + v[4];
+ ztail = z + v[5];
+ color_from_column(FALSE);
+ }
+ }
+#undef color_from_column
+
+
+ /* Adjust for logscales. Set min/max and point types. Store in cp.
+ * The macro cannot use continue, as it is wrapped in a loop.
+ * I regard this as correct goto use
+ */
+ cp->type = INRANGE;
+ 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);
+ if (this_plot->plot_style == VECTOR) {
+ cptail->type = INRANGE;
+ STORE_WITH_LOG_AND_UPDATE_RANGE(cptail->x, xtail, cp->type, x_axis, NOOP, goto come_here_if_undefined);
+ STORE_WITH_LOG_AND_UPDATE_RANGE(cptail->y, ytail, cp->type, y_axis, NOOP, goto come_here_if_undefined);
+ }
+
+ if (dgrid3d) {
+ /* HBB 20010424: in dgrid3d mode, delay log() taking
+ * and scaling until after the dgrid process. Only for
+ * z, not for x and y, so we can layout the newly
+ * created created grid more easily. */
+ cp->z = z;
+ if (this_plot->plot_style == VECTOR)
+ cptail->z = ztail;
+ } else {
+ STORE_WITH_LOG_AND_UPDATE_RANGE(cp->z, z, cp->type, z_axis, NOOP, goto come_here_if_undefined);
+ if (this_plot->plot_style == VECTOR)
+ STORE_WITH_LOG_AND_UPDATE_RANGE(cptail->z, ztail, cp->type, z_axis, NOOP, goto come_here_if_undefined);
+
+ if (NEED_PALETTE(this_plot)) {
+ if (pm3d_color_from_column) {
+ COLOR_STORE_WITH_LOG_AND_UPDATE_RANGE(cp->CRD_COLOR, color, cp->type, COLOR_AXIS, NOOP, goto come_here_if_undefined);
+ } else {
+ COLOR_STORE_WITH_LOG_AND_UPDATE_RANGE(cp->CRD_COLOR, z, cp->type, COLOR_AXIS, NOOP, goto come_here_if_undefined);
+ }
+ }
+ }
+
+#ifdef EAM_DATASTRINGS
+ /* At this point we have stored the point coordinates. Now we need to copy */
+ /* x,y,z into the text_label structure and add the actual text string. */
+ if (this_plot->plot_style == LABELPOINTS)
+ store_label(this_plot->labels, cp, xdatum, df_tokens[3], color);
+#endif
+
+ come_here_if_undefined:
+ /* some may complain, but I regard this as the correct use of goto */
+ ++xdatum;
+ } /* end of whileloop - end of surface */
+
+ if (pm3d_color_from_column) {
+ this_plot->pm3d_color_from_column = pm3d_color_from_column;
+ }
+
+ if (xdatum > 0) {
+ this_plot->num_iso_read++; /* Update last iso. */
+ local_this_iso->p_count = xdatum;
+
+ /* If this is a VECTOR plot then iso->next is already */
+ /* occupied by the vector tail coordinates. */
+ if (this_plot->plot_style != VECTOR)
+ 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 { /* Free last allocation */
+ if (this_plot->plot_style == VECTOR)
+ iso_free(local_this_iso->next);
+ iso_free(local_this_iso);
+ }
+
+ /*}}} */
+ }
+
+ if (dgrid3d && this_plot->num_iso_read > 0)
+ grid_nongrid_data(this_plot);
+
+ /* This check used to be done in graph3d */
+ if (X_AXIS.min == VERYLARGE || X_AXIS.max == -VERYLARGE ||
+ Y_AXIS.min == VERYLARGE || Y_AXIS.max == -VERYLARGE ||
+ Z_AXIS.min == VERYLARGE || Z_AXIS.max == -VERYLARGE)
+ int_warn(NO_CARET,
+ "Axis range undefined due to improper data values. NaN? Inf?");
+
+ if (this_plot->num_iso_read <= 1)
+ this_plot->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;
+ int i;
+
+ /* 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(int pcount)
+{
+ struct surface_points *this_plot;
+ int i, surface;
+ struct coordinate GPHUGE *point;
+ char *buffer = gp_alloc(150, "print_3dtable output buffer");
+ FILE *outfile = (table_outfile) ? table_outfile : gpoutfile;
+
+ for (surface = 0, this_plot = first_3dplot;
+ surface < pcount;
+ this_plot = this_plot->next_sp, surface++) {
+ fprintf(outfile, "\n#Surface %d of %d surfaces\n", surface, pcount);
+
+ if (draw_surface) {
+ struct iso_curve *icrvs;
+ int curve;
+
+ /* only the curves in one direction */
+ for (curve = 0, icrvs = this_plot->iso_crvs;
+ icrvs && curve < this_plot->num_iso_read;
+ icrvs = icrvs->next, curve++) {
+ fprintf(outfile, "\n#IsoCurve %d, %d points\n#x y z type\n",
+ curve, icrvs->p_count);
+ for (i = 0, point = icrvs->points;
+ i < icrvs->p_count;
+ i++, point++) {
+ /* HBB 20020405: new macro to use the 'set format'
+ * in their full flexibility */
+#define OUTPUT_NUMBER(field, axis) \
+ gprintf(buffer, 150, axis_array[axis].formatstring, \
+ 1.0, point->field); \
+ fputs(buffer, outfile); \
+ fputc(' ', outfile);
+ OUTPUT_NUMBER(x, FIRST_X_AXIS);
+ OUTPUT_NUMBER(y, FIRST_Y_AXIS);
+ OUTPUT_NUMBER(z, FIRST_Z_AXIS);
+ fprintf(outfile, "%c\n",
+ point->type == INRANGE
+ ? 'i' : point->type == OUTRANGE
+ ? 'o' : 'u');
+ } /* for(point) */
+ } /* for(icrvs) */
+ putc('\n', outfile);
+ } /* if(draw_surface) */
+
+ if (draw_contour) {
+ int number = 0;
+ struct gnuplot_contours *c = this_plot->contours;
+
+ while (c) {
+ int count = c->num_pts;
+ struct coordinate GPHUGE *point = c->coords;
+
+ if (c->isNewLevel)
+ /* don't 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(outfile, "\n# Contour %d, label: %s\n",
+ number++, c->label);
+
+ for (; --count >= 0; ++point) {
+ OUTPUT_NUMBER(x, FIRST_X_AXIS);
+ OUTPUT_NUMBER(y, FIRST_Y_AXIS);
+ OUTPUT_NUMBER(z, FIRST_Z_AXIS);
+ putc('\n', outfile);
+ }
+
+ /* blank line between segments of same contour */
+ putc('\n', outfile);
+ c = c->next;
+#undef OUTPUT_NUMBER
+ } /* while (contour) */
+ } /* if (draw_contour) */
+ } /* for(surface) */
+ fflush(outfile);
+
+ free(buffer);
+}
+
+/* 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(
+ 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,
+ TBOOLEAN need_palette)
+{
+ int i, j;
+ struct coordinate GPHUGE *points = (*this_iso)->points;
+ int do_update_color = need_palette && (!parametric || (parametric && value_axis == FIRST_Z_AXIS));
+
+ 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);
+ if (do_update_color) {
+ COLOR_STORE_WITH_LOG_AND_UPDATE_RANGE(points[i].CRD_COLOR, temp, points[i].type, COLOR_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;
+ legend_key *key = &keyT;
+
+ /* 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 plot expected");
+
+ start_token = c_token;
+
+ if (is_definition(c_token)) {
+ define();
+ } else {
+ int specs = -1;
+ struct surface_points *this_plot;
+
+ char *name_str;
+ TBOOLEAN duplication = FALSE;
+ TBOOLEAN set_title = FALSE, set_with = FALSE;
+ TBOOLEAN set_lpstyle = FALSE;
+ TBOOLEAN checked_once = FALSE;
+#ifdef EAM_DATASTRINGS
+ TBOOLEAN set_labelstyle = FALSE;
+#endif
+
+ dummy_func = &plot_func;
+ /* WARNING: do NOT free name_str */
+ /* FIXME: could also be saved in this_plot */
+ name_str = string_or_express(NULL);
+ dummy_func = NULL;
+ if (name_str) {
+ /*{{{ data file to plot */
+ 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 & AUTOSCALE_MIN) {
+ axis_array[FIRST_X_AXIS].min = VERYLARGE;
+ }
+ if (axis_array[FIRST_X_AXIS].autoscale & AUTOSCALE_MAX) {
+ axis_array[FIRST_X_AXIS].max = -VERYLARGE;
+ }
+ if (axis_array[FIRST_Y_AXIS].autoscale & AUTOSCALE_MIN) {
+ axis_array[FIRST_Y_AXIS].min = VERYLARGE;
+ }
+ if (axis_array[FIRST_Y_AXIS].autoscale & AUTOSCALE_MAX) {
+ 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;
+
+ df_set_plot_mode(MODE_SPLOT);
+ specs = df_open(name_str, MAXDATACOLS);
+#ifdef BINARY_DATA_FILE
+ if (df_matrix)
+ this_plot->has_grid_topology = TRUE;
+#endif
+
+ /* 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_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;
+ /* ignore it for now */
+ some_functions = TRUE;
+ end_token = c_token - 1;
+ /*}}} */
+
+ } /* end of IS THIS A FILE OR A FUNC block */
+
+ /* clear current title, if exist */
+ if (this_plot->title) {
+ free(this_plot->title);
+ this_plot->title = NULL;
+ }
+
+ /* default line and point types, no palette */
+ this_plot->lp_properties.use_palette = 0;
+ this_plot->lp_properties.l_type = line_num;
+ this_plot->lp_properties.p_type = point_num;
+
+ /* user may prefer explicit line styles */
+ if (prefer_line_styles)
+ lp_use_properties(&this_plot->lp_properties, line_num+1, TRUE);
+
+ /* pm 25.11.2001 allow any order of options */
+ while (!END_OF_COMMAND || !checked_once) {
+
+ /* deal with title */
+ if (almost_equals(c_token, "t$itle")) {
+ if (set_title) {
+ duplication=TRUE;
+ break;
+ }
+ this_plot->title_no_enhanced = !key->enhanced;
+ /* title can be enhanced if not explicitly disabled */
+ 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 . */
+ }
+ }
+ c_token++;
+ if (!(this_plot->title = try_to_get_string()))
+ int_error(c_token, "expecting \"title\" for plot");
+ set_title = TRUE;
+ continue;
+ }
+
+ if (almost_equals(c_token, "not$itle")) {
+ if (set_title) {
+ duplication=TRUE;
+ break;
+ }
+ c_token++;
+ if (isstringvalue(c_token))
+ try_to_get_string(); /* ignore optionally given title string */
+ if (xtitle != NULL)
+ xtitle[0] = '\0';
+ if (ytitle != NULL)
+ ytitle[0] = '\0';
+ /* this_plot->title = NULL; */
+ set_title = TRUE;
+ continue;
+ }
+
+ /* deal with style */
+ if (almost_equals(c_token, "w$ith")) {
+ if (set_with) {
+ duplication=TRUE;
+ break;
+ }
+ this_plot->plot_style = get_style();
+ if ((this_plot->plot_type == FUNC3D) &&
+ ((this_plot->plot_style & PLOT_STYLE_HAS_ERRORBAR)
+#ifdef EAM_DATASTRINGS
+ || (this_plot->plot_style == LABELPOINTS)
+#endif
+ )) {
+ int_warn(c_token, "This plot style is only for datafiles , reverting to \"points\"");
+ this_plot->plot_style = POINTSTYLE;
+ }
+
+ if ((this_plot->plot_style | data_style) & PM3DSURFACE) {
+ if (equals(c_token, "at")) {
+ /* option 'with pm3d [at ...]' is explicitly specified */
+ c_token++;
+ if (get_pm3d_at_option(&this_plot->pm3d_where[0]))
+ return; /* error */
+ }
+ }
+
+ set_with = TRUE;
+ continue;
+ }
+
+ /* EAM Dec 2005 - Hidden3D code by default includes points, */
+ /* labels and vectors in the hidden3d processing. Check here */
+ /* if this particular plot wants to be excluded. */
+ if (almost_equals(c_token, "nohidden$3d")) {
+ c_token++;
+ this_plot->opt_out_of_hidden3d = TRUE;
+ continue;
+ }
+
+ /* "set contour" is global. Allow individual plots to opt out */
+ if (almost_equals(c_token, "nocon$tours")) {
+ c_token++;
+ this_plot->opt_out_of_contours = TRUE;
+ continue;
+ }
+
+#ifdef EAM_DATASTRINGS
+ /* Labels can have font and text property info as plot options */
+ /* In any case we must allocate one instance of the text style */
+ /* that all labels in the plot will share. */
+ if (this_plot->plot_style == LABELPOINTS) {
+ int stored_token = c_token;
+ if (this_plot->labels == NULL) {
+ this_plot->labels = new_text_label(-1);
+ this_plot->labels->pos = JUST_CENTRE;
+ this_plot->labels->layer = LAYER_PLOTLABELS;
+ }
+ parse_label_options(this_plot->labels);
+ checked_once = TRUE;
+ if (stored_token != c_token) {
+ if (set_labelstyle) {
+ duplication = TRUE;
+ break;
+ } else {
+ set_labelstyle = TRUE;
+ continue;
+ }
+ }
+ }
+#endif
+
+ /* pick up line/point specs
+ * - point spec allowed if style uses points, ie style&2 != 0
+ * - keywords are optional
+ */
+ if (this_plot->plot_style == VECTOR) {
+ int stored_token = c_token;
+
+ if (!checked_once) {
+ default_arrow_style(&this_plot->arrow_properties);
+ this_plot->arrow_properties.lp_properties.l_type = line_num;
+ checked_once = TRUE;
+ }
+ arrow_parse(&this_plot->arrow_properties, TRUE);
+ if (stored_token != c_token) {
+ if (set_lpstyle) {
+ duplication = TRUE;
+ break;
+ } else {
+ set_lpstyle = TRUE;
+ this_plot->lp_properties = this_plot->arrow_properties.lp_properties;
+ continue;
+ }
+ }
+ } else {
+ int stored_token = c_token;
+ struct lp_style_type lp = DEFAULT_LP_STYLE_TYPE;
+
+ lp.l_type = line_num;
+ lp.p_type = point_num;
+
+ /* user may prefer explicit line styles */
+ if (prefer_line_styles)
+ lp_use_properties(&lp, line_num+1, TRUE);
+
+ lp_parse(&lp, TRUE,
+ this_plot->plot_style & PLOT_STYLE_HAS_POINT);
+ checked_once = TRUE;
+ if (stored_token != c_token) {
+ if (set_lpstyle) {
+ duplication=TRUE;
+ break;
+ } else {
+ this_plot->lp_properties = lp;
+ set_lpstyle = TRUE;
+ continue;
+ }
+ }
+ }
+
+ break; /* unknown option */
+
+ } /* while (!END_OF_COMMAND)*/
+
+ if (duplication)
+ int_error(c_token, "duplicated or contradicting arguments in plot options");
+
+ /* set default values for title if this has not been specified */
+ if (!set_title) {
+ this_plot->title_no_enhanced = TRUE; /* filename or function cannot be enhanced */
+ if (key->auto_titles == FILENAME_KEYTITLES) {
+#ifdef BINARY_DATA_FILE
+ /* DJS (20 Aug 2004) I'd prefer that the df_binary flag be discarded. There
+ * is nothing special about the file being binary that its title should be
+ * different. Can't the decision to do this be based on some other criteria,
+ * like the presence of a nonconventional `using`?
+ */
+#endif
+ if (this_plot->plot_type == DATA3D && df_binary==TRUE && end_token==start_token+1)
+ /* let default title for splot 'a.dat' binary is 'a.dat'
+ * while for 'a.dat' binary using 2:1:3 will be all 4 words */
+ m_capture(&(this_plot->title), start_token, start_token);
+ else
+ 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;
+ } else {
+ if (xtitle != NULL)
+ xtitle[0] = '\0';
+ if (ytitle != NULL)
+ ytitle[0] = '\0';
+ /* this_plot->title = NULL; */
+ }
+ }
+
+ /* No line/point style given. As lp_parse also supplies
+ * the defaults for linewidth and pointsize, call it now
+ * to define them. */
+ if (! set_lpstyle) {
+ if (this_plot->plot_style == VECTOR) {
+ this_plot->arrow_properties.lp_properties.l_type = line_num;
+ arrow_parse(&this_plot->arrow_properties, TRUE);
+ this_plot->lp_properties = this_plot->arrow_properties.lp_properties;
+ } else {
+ this_plot->lp_properties.l_type = line_num;
+ this_plot->lp_properties.l_width = 1.0;
+ this_plot->lp_properties.p_type = point_num;
+ this_plot->lp_properties.p_size = pointsize;
+ this_plot->lp_properties.use_palette = 0;
+
+ /* user may prefer explicit line styles */
+ if (prefer_line_styles)
+ lp_use_properties(&this_plot->lp_properties, line_num+1, TRUE);
+
+ lp_parse(&this_plot->lp_properties, TRUE,
+ this_plot->plot_style & PLOT_STYLE_HAS_POINT);
+ }
+
+#ifdef BACKWARDS_COMPATIBLE
+ /* allow old-style syntax - ignore case lt 3 4 for example */
+ if (!END_OF_COMMAND && isanumber(c_token)) {
+ struct value t;
+
+ this_plot->lp_properties.l_type =
+ this_plot->lp_properties.p_type =
+ (int) real(const_express(&t)) - 1;
+
+ if (isanumber(c_token))
+ this_plot->lp_properties.p_type =
+ (int) real(const_express(&t)) - 1;
+ }
+#endif
+
+ }
+
+ /* Rule out incompatible line/point/style options */
+ if (this_plot->plot_type == FUNC3D) {
+ if ((this_plot->plot_style & PLOT_STYLE_HAS_POINT)
+ && (this_plot->lp_properties.p_size == PTSZ_VARIABLE))
+ this_plot->lp_properties.p_size = 1;
+ }
+
+ if (crnt_param == 0
+ && this_plot->plot_style != PM3DSURFACE
+ /* don't increment the default line/point properties if
+ * this_plot is an EXPLICIT pm3d surface plot */
+#ifdef WITH_IMAGE
+ && this_plot->plot_style != IMAGE
+ && this_plot->plot_style != RGBIMAGE
+ /* same as above, for an (rgb)image plot */
+#endif
+ ) {
+ if (this_plot->plot_style & PLOT_STYLE_HAS_POINT)
+ point_num += 1 + (draw_contour != 0) + (hidden3d != 0);
+ line_num += 1 + (draw_contour != 0) + (hidden3d != 0);
+ }
+
+#ifdef WITH_IMAGE
+ /* Styles that utilize palettes. */
+ if (this_plot->plot_style == IMAGE)
+ this_plot->lp_properties.use_palette = 1;
+#endif
+
+ /* 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;
+ struct surface_points *first_dataset = this_plot;
+ /* pointer to the plot of the first dataset (surface) in the file */
+ int this_token = this_plot->token;
+
+ do {
+ 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)
+ this_plot->plot_type = NODATA;
+
+ if (this_plot != first_dataset)
+ /* copy (explicit) "with pm3d at ..." option from the first dataset in the file */
+ strcpy(this_plot->pm3d_where, first_dataset->pm3d_where);
+
+ /* 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;
+ } while (!df_eof);
+
+ 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.
+ */
+ 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 & AUTOSCALE_MIN)
+ axis_array[FIRST_X_AXIS].min = VERYLARGE;
+ if (axis_array[FIRST_X_AXIS].autoscale & AUTOSCALE_MAX)
+ axis_array[FIRST_X_AXIS].max = -VERYLARGE;
+ if (axis_array[FIRST_Y_AXIS].autoscale & AUTOSCALE_MIN)
+ axis_array[FIRST_Y_AXIS].min = VERYLARGE;
+ if (axis_array[FIRST_Y_AXIS].autoscale & AUTOSCALE_MAX)
+ 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 {
+ struct at_type *at_ptr;
+ char *name_str;
+
+ dummy_func = &plot_func;
+ name_str = string_or_express(&at_ptr);
+
+ if (!name_str) { /* 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;
+
+ plot_func.at = at_ptr;
+
+ 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,
+ NEED_PALETTE(this_plot));
+
+ 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,
+ NEED_PALETTE(this_plot));
+ }
+ /*}}} */
+ } /* 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");
+ if (splot_map)
+ axis_checked_extend_empty_range(FIRST_Z_AXIS, NULL); /* Suppress warning message */
+ else
+ 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);
+
+ set_plot_with_palette(plot_num, MODE_SPLOT);
+ if (is_plot_with_palette()) {
+ set_cbminmax();
+ axis_checked_extend_empty_range(COLOR_AXIS, "All points of colorbox value undefined");
+ setup_tics(COLOR_AXIS, 20);
+ /* axis_revert_and_unlog_range(COLOR_AXIS); */
+ /* fprintf(stderr,"plot3d.c: CB_AXIS.min=%g\tCB_AXIS.max=%g\n",CB_AXIS.min,CB_AXIS.max); */
+ }
+
+ 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);
+ }
+ this_plot->contours = NULL;
+ }
+
+ /* Allow individual surfaces to opt out of contouring */
+ if (this_plot->opt_out_of_contours)
+ continue;
+
+ /* Make sure this one can be contoured. */
+ if (!this_plot->has_grid_topology) {
+ int_warn(NO_CARET, "Notice: Cannot contour non grid data. Please use \"set dgrid3d\".");
+ } 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 */
+
+ /* the following ~9 lines were moved from the end of the
+ * function to here, as do_3dplot calles term->text, which
+ * itself might process input events in mouse enhanced
+ * terminals. For redrawing to work, line capturing and
+ * setting the plot3d_num must already be done before
+ * entering do_3dplot(). Thu Jan 27 23:54:49 2000 (joze) */
+
+ /* 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;
+ }
+
+ /* record that all went well */
+ plot3d_num=plot_num;
+
+ /* perform the plot */
+ if (table_mode)
+ 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;
+ /* update GPVAL_ variables available to user */
+ update_gpval_variables(1);
+ }
+}
+
+
+
+/*
+ * 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
+ *
+ * x and y ranges now fixed in eval_3dplots
+ */
+static void
+parametric_3dfixup(struct surface_points *start_plot, int *plot_num)
+{
+ 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;
+ struct coordinate GPHUGE *ypoints = yicrvs->points;
+ struct coordinate 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;
+ }
+
+#if 0 /* FIXME HBB 20001101: seems to cause a crash */
+ if (first_3dplot)
+ sp_free(first_3dplot);
+ first_3dplot = NULL;
+#endif
+
+ /* 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;
+}
+