#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 <= i - 1; j++) sum -= *ac++ * *bp++; } else if (sum) memi = i; b[i] = sum; } ar--; for (i = n - 1; i >= 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 <= 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; /* 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; }