git.maemo.org Git - gnuplot/blob - src/plot3d.c

   1 #ifndef lint
   2 static char *RCSid() { return RCSid("$Id: plot3d.c,v 1.131.2.6 2008/03/20 09:39:46 sfeam Exp $"); }
   3 #endif
   4
   5 /* GNUPLOT - plot3d.c */
   6
   7 /*[
   8  * Copyright 1986 - 1993, 1998, 2004   Thomas Williams, Colin Kelley
   9  *
  10  * Permission to use, copy, and distribute this software and its
  11  * documentation for any purpose with or without fee is hereby granted,
  12  * provided that the above copyright notice appear in all copies and
  13  * that both that copyright notice and this permission notice appear
  14  * in supporting documentation.
  15  *
  16  * Permission to modify the software is granted, but not the right to
  17  * distribute the complete modified source code.  Modifications are to
  18  * be distributed as patches to the released version.  Permission to
  19  * distribute binaries produced by compiling modified sources is granted,
  20  * provided you
  21  *   1. distribute the corresponding source modifications from the
  22  *    released version in the form of a patch file along with the binaries,
  23  *   2. add special version identification to distinguish your version
  24  *    in addition to the base release version number,
  25  *   3. provide your name and address as the primary contact for the
  26  *    support of your modified version, and
  27  *   4. retain our contact information in regard to use of the base
  28  *    software.
  29  * Permission to distribute the released version of the source code along
  30  * with corresponding source modifications in the form of a patch file is
  31  * granted with same provisions 2 through 4 for binary distributions.
  32  *
  33  * This software is provided "as is" without express or implied warranty
  34  * to the extent permitted by applicable law.
  35 ]*/
  36
  37 #include "plot3d.h"
  38 #include "gp_types.h"
  39
  40 #include "alloc.h"
  41 #include "axis.h"
  42 #include "binary.h"
  43 #include "command.h"
  44 #include "contour.h"
  45 #include "datafile.h"
  46 #include "eval.h"
  47 #include "graph3d.h"
  48 #include "misc.h"
  49 #include "parse.h"
  50 #include "setshow.h"
  51 #include "term_api.h"
  52 #include "util.h"
  53 #include "pm3d.h"
  54 #ifdef BINARY_DATA_FILE
  55 #include "plot.h"
  56 #endif
  57
  58 #ifdef EAM_DATASTRINGS
  59 #include "plot2d.h" /* Only for store_label() */
  60 #endif
  61
  62 #ifdef THIN_PLATE_SPLINES_GRID
  63 # include "matrix.h"
  64 #endif
  65
  66 #ifndef _Windows
  67 # include "help.h"
  68 #endif
  69
  70 /* global variables exported by this module */
  71
  72 t_data_mapping mapping3d = MAP3D_CARTESIAN;
  73
  74 int dgrid3d_row_fineness = 10;
  75 int dgrid3d_col_fineness = 10;
  76 int dgrid3d_norm_value = 1;
  77 TBOOLEAN dgrid3d = FALSE;
  78
  79 /* static prototypes */
  80
  81 static void calculate_set_of_isolines __PROTO((AXIS_INDEX value_axis, TBOOLEAN cross, struct iso_curve **this_iso,
  82                                                AXIS_INDEX iso_axis, double iso_min, double iso_step, int num_iso_to_use,
  83                                                AXIS_INDEX sam_axis, double sam_min, double sam_step, int num_sam_to_use,
  84                                                TBOOLEAN need_palette));
  85 static void get_3ddata __PROTO((struct surface_points * this_plot));
  86 static void print_3dtable __PROTO((int pcount));
  87 static void eval_3dplots __PROTO((void));
  88 static void grid_nongrid_data __PROTO((struct surface_points * this_plot));
  89 static void parametric_3dfixup __PROTO((struct surface_points * start_plot, int *plot_num));
  90 static struct surface_points * sp_alloc __PROTO((int num_samp_1, int num_iso_1, int num_samp_2, int num_iso_2));
  91 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));
  92 #ifdef THIN_PLATE_SPLINES_GRID
  93 static double splines_kernel __PROTO((double h));
  94 #endif
  95
  96 /* the curves/surfaces of the plot */
  97 struct surface_points *first_3dplot = NULL;
  98 static struct udft_entry plot_func;
  99
 100 int plot3d_num=0;
 101
 102 /* HBB 20000508: moved these functions to the only module that uses them
 103  * so they can be turned 'static' */
 104 /*
 105  * sp_alloc() allocates a surface_points structure that can hold 'num_iso_1'
 106  * iso-curves with 'num_samp_2' samples and 'num_iso_2' iso-curves with
 107  * 'num_samp_1' samples.
 108  * If, however num_iso_2 or num_samp_1 is zero no iso curves are allocated.
 109  */
 110 static struct surface_points *
 111 sp_alloc(int num_samp_1, int num_iso_1, int num_samp_2, int num_iso_2)
 112 {
 113     struct lp_style_type default_lp_properties = DEFAULT_LP_STYLE_TYPE;
 114
 115     struct surface_points *sp = gp_alloc(sizeof(*sp), "surface");
 116     memset(sp,0,sizeof(struct surface_points));
 117
 118     /* Initialize various fields */
 119     sp->lp_properties = default_lp_properties;
 120     default_arrow_style(&(sp->arrow_properties));
 121
 122     if (num_iso_2 > 0 && num_samp_1 > 0) {
 123         int i;
 124         struct iso_curve *icrv;
 125
 126         for (i = 0; i < num_iso_1; i++) {
 127             icrv = iso_alloc(num_samp_2);
 128             icrv->next = sp->iso_crvs;
 129             sp->iso_crvs = icrv;
 130         }
 131         for (i = 0; i < num_iso_2; i++) {
 132             icrv = iso_alloc(num_samp_1);
 133             icrv->next = sp->iso_crvs;
 134             sp->iso_crvs = icrv;
 135         }
 136     }
 137
 138     return (sp);
 139 }
 140
 141 /*
 142  * sp_replace() updates a surface_points structure so it can hold 'num_iso_1'
 143  * iso-curves with 'num_samp_2' samples and 'num_iso_2' iso-curves with
 144  * 'num_samp_1' samples.
 145  * If, however num_iso_2 or num_samp_1 is zero no iso curves are allocated.
 146  */
 147 static void
 148 sp_replace(
 149     struct surface_points *sp,
 150     int num_samp_1, int num_iso_1, int num_samp_2, int num_iso_2)
 151 {
 152     int i;
 153     struct iso_curve *icrv, *icrvs = sp->iso_crvs;
 154
 155     while (icrvs) {
 156         icrv = icrvs;
 157         icrvs = icrvs->next;
 158         iso_free(icrv);
 159     }
 160     sp->iso_crvs = NULL;
 161
 162     if (num_iso_2 > 0 && num_samp_1 > 0) {
 163         for (i = 0; i < num_iso_1; i++) {
 164             icrv = iso_alloc(num_samp_2);
 165             icrv->next = sp->iso_crvs;
 166             sp->iso_crvs = icrv;
 167         }
 168         for (i = 0; i < num_iso_2; i++) {
 169             icrv = iso_alloc(num_samp_1);
 170             icrv->next = sp->iso_crvs;
 171             sp->iso_crvs = icrv;
 172         }
 173     } else
 174         sp->iso_crvs = NULL;
 175 }
 176
 177 /*
 178  * sp_free() releases any memory which was previously malloc()'d to hold
 179  *   surface points.
 180  */
 181 /* HBB 20000506: don't risk stack havoc by recursion, use iterative list
 182  * cleanup unstead */
 183 void
 184 sp_free(struct surface_points *sp)
 185 {
 186     while (sp) {
 187         struct surface_points *next = sp->next_sp;
 188         if (sp->title)
 189             free(sp->title);
 190
 191         while (sp->contours) {
 192             struct gnuplot_contours *next_cntrs = sp->contours->next;
 193
 194             free(sp->contours->coords);
 195             free(sp->contours);
 196             sp->contours = next_cntrs;
 197         }
 198
 199         while (sp->iso_crvs) {
 200             struct iso_curve *next_icrvs = sp->iso_crvs->next;
 201
 202             iso_free(sp->iso_crvs);
 203             sp->iso_crvs = next_icrvs;
 204         }
 205
 206 #ifdef EAM_DATASTRINGS
 207         if (sp->labels) {
 208             free_labels(sp->labels);
 209             sp->labels = (struct text_label *)NULL;
 210         }
 211 #endif
 212
 213         free(sp);
 214         sp = next;
 215     }
 216 }
 217
 218
 219
 220 /* support for dynamic size of input line */
 221
 222 void
 223 plot3drequest()
 224 /*
 225  * in the parametric case we would say splot [u= -Pi:Pi] [v= 0:2*Pi] [-1:1]
 226  * [-1:1] [-1:1] sin(v)*cos(u),sin(v)*cos(u),sin(u) in the non-parametric
 227  * case we would say only splot [x= -2:2] [y= -5:5] sin(x)*cos(y)
 228  *
 229  */
 230 {
 231     int dummy_token0 = -1, dummy_token1 = -1;
 232     AXIS_INDEX u_axis, v_axis;
 233
 234     is_3d_plot = TRUE;
 235
 236     /* change view to become map if requested by 'set view map' */
 237     if (splot_map == TRUE)
 238         splot_map_activate();
 239
 240     if (parametric && strcmp(set_dummy_var[0], "t") == 0) {
 241         strcpy(set_dummy_var[0], "u");
 242         strcpy(set_dummy_var[1], "v");
 243     }
 244
 245     /* put stuff into arrays to simplify access */
 246     AXIS_INIT3D(FIRST_X_AXIS, 0, 0);
 247     AXIS_INIT3D(FIRST_Y_AXIS, 0, 0);
 248     AXIS_INIT3D(FIRST_Z_AXIS, 0, 1);
 249     AXIS_INIT3D(U_AXIS, 1, 0);
 250     AXIS_INIT3D(V_AXIS, 1, 0);
 251     AXIS_INIT3D(COLOR_AXIS, 0, 1);
 252
 253     if (!term)                  /* unknown */
 254         int_error(c_token, "use 'set term' to set terminal type first");
 255
 256     u_axis = (parametric ? U_AXIS : FIRST_X_AXIS);
 257     v_axis = (parametric ? V_AXIS : FIRST_Y_AXIS);
 258
 259     PARSE_NAMED_RANGE(u_axis, dummy_token0);
 260     if (splot_map == TRUE && !parametric) /* v_axis==FIRST_Y_AXIS */
 261         splot_map_deactivate();
 262     PARSE_NAMED_RANGE(v_axis, dummy_token1);
 263     if (splot_map == TRUE && !parametric) /* v_axis==FIRST_Y_AXIS */
 264         splot_map_activate();
 265
 266     if (parametric) {
 267         PARSE_RANGE(FIRST_X_AXIS);
 268         if (splot_map == TRUE)
 269             splot_map_deactivate();
 270         PARSE_RANGE(FIRST_Y_AXIS);
 271         if (splot_map == TRUE)
 272             splot_map_activate();
 273     }                           /* parametric */
 274     PARSE_RANGE(FIRST_Z_AXIS);
 275     CHECK_REVERSE(FIRST_X_AXIS);
 276     CHECK_REVERSE(FIRST_Y_AXIS);
 277     CHECK_REVERSE(FIRST_Z_AXIS);
 278
 279     /* use the default dummy variable unless changed */
 280     if (dummy_token0 >= 0)
 281         copy_str(c_dummy_var[0], dummy_token0, MAX_ID_LEN);
 282     else
 283         (void) strcpy(c_dummy_var[0], set_dummy_var[0]);
 284
 285     if (dummy_token1 >= 0)
 286         copy_str(c_dummy_var[1], dummy_token1, MAX_ID_LEN);
 287     else
 288         (void) strcpy(c_dummy_var[1], set_dummy_var[1]);
 289
 290     eval_3dplots();
 291 }
 292
 293 #ifdef THIN_PLATE_SPLINES_GRID
 294
 295 static double
 296 splines_kernel(double h)
 297 {
 298 #if 0
 299     /* this is normaly not useful ... */
 300     h = fabs(h);
 301
 302     if (h != 0.0) {
 303 #else
 304     if (h > 0) {
 305 #endif
 306         return h * h * log(h);
 307     } else {
 308         return 0;
 309     }
 310 }
 311
 312 #endif
 313
 314 static void
 315 grid_nongrid_data(struct surface_points *this_plot)
 316 {
 317     int i, j, k;
 318     double x, y, z, w, dx, dy, xmin, xmax, ymin, ymax;
 319     struct iso_curve *old_iso_crvs = this_plot->iso_crvs;
 320     struct iso_curve *icrv, *oicrv, *oicrvs;
 321
 322 #ifdef THIN_PLATE_SPLINES_GRID
 323     double *b, **K, *xx, *yy, *zz, d;
 324     int *indx, numpoints;
 325 #endif
 326
 327     /* Compute XY bounding box on the original data. */
 328     /* FIXME HBB 20010424: Does this make any sense? Shouldn't we just
 329      * use whatever the x and y ranges have been found to be, and
 330      * that's that? The largest difference this is going to make is if
 331      * we plot a datafile that doesn't span the whole x/y range
 332      * used. Do we want a dgrid3d over the actual data rectangle, or
 333      * over the xrange/yrange area? */
 334     xmin = xmax = old_iso_crvs->points[0].x;
 335     ymin = ymax = old_iso_crvs->points[0].y;
 336     for (icrv = old_iso_crvs; icrv != NULL; icrv = icrv->next) {
 337         struct coordinate GPHUGE *points = icrv->points;
 338
 339         for (i = 0; i < icrv->p_count; i++, points++) {
 340             /* HBB 20010424: avoid crashing for undefined input */
 341             if (points->type == UNDEFINED)
 342                 continue;
 343             if (xmin > points->x)
 344                 xmin = points->x;
 345             if (xmax < points->x)
 346                 xmax = points->x;
 347             if (ymin > points->y)
 348                 ymin = points->y;
 349             if (ymax < points->y)
 350                 ymax = points->y;
 351         }
 352     }
 353
 354     dx = (xmax - xmin) / (dgrid3d_col_fineness - 1);
 355     dy = (ymax - ymin) / (dgrid3d_row_fineness - 1);
 356
 357     /* Create the new grid structure, and compute the low pass filtering from
 358      * non grid to grid structure.
 359      */
 360     this_plot->iso_crvs = NULL;
 361     this_plot->num_iso_read = dgrid3d_col_fineness;
 362     this_plot->has_grid_topology = TRUE;
 363
 364 #ifdef THIN_PLATE_SPLINES_GRID
 365     numpoints = 0;
 366     for (oicrv = old_iso_crvs; oicrv != NULL; oicrv = oicrv->next) {
 367         numpoints += oicrv->p_count;
 368     }
 369     xx = gp_alloc(sizeof(xx[0]) * (numpoints + 3) * (numpoints + 8),
 370                   "thin plate splines in dgrid3d");
 371     /* the memory needed is not really (n+3)*(n+8) for now,
 372        but might be if I take into account errors ... */
 373     K = gp_alloc(sizeof(K[0]) * (numpoints + 3),
 374                  "matrix : thin plate splines 2d");
 375     yy = xx + numpoints;
 376     zz = yy + numpoints;
 377     b = zz + numpoints;
 378
 379     /* HBB 20010424: Count actual input points without the UNDEFINED
 380      * ones, as we copy them */
 381     numpoints = 0;
 382     for (oicrv = old_iso_crvs; oicrv != NULL; oicrv = oicrv->next) {
 383         struct coordinate GPHUGE *opoints = oicrv->points;
 384
 385         for (k = 0; k < oicrv->p_count; k++, opoints++) {
 386             /* HBB 20010424: avoid crashing for undefined input */
 387             if (opoints->type == UNDEFINED)
 388                 continue;
 389             xx[numpoints] = opoints->x;
 390             yy[numpoints] = opoints->y;
 391             zz[numpoints] = opoints->z;
 392             numpoints++;
 393         }
 394     }
 395
 396     for (i = 0; i < numpoints + 3; i++) {
 397         K[i] = b + (numpoints + 3) * (i + 1);
 398     }
 399
 400     for (i = 0; i < numpoints; i++) {
 401         for (j = i + 1; j < numpoints; j++) {
 402             double dx = xx[i] - xx[j], dy = yy[i] - yy[j];
 403             K[i][j] = K[j][i] = -splines_kernel(sqrt(dx * dx + dy * dy));
 404         }
 405         K[i][i] = 0.0;          /* here will come the weights for errors */
 406         b[i] = zz[i];
 407     }
 408     for (i = 0; i < numpoints; i++) {
 409         K[i][numpoints] = K[numpoints][i] = 1.0;
 410         K[i][numpoints + 1] = K[numpoints + 1][i] = xx[i];
 411         K[i][numpoints + 2] = K[numpoints + 2][i] = yy[i];
 412     }
 413     b[numpoints] = 0.0;
 414     b[numpoints + 1] = 0.0;
 415     b[numpoints + 2] = 0.0;
 416     K[numpoints][numpoints] = 0.0;
 417     K[numpoints][numpoints + 1] = 0.0;
 418     K[numpoints][numpoints + 2] = 0.0;
 419     K[numpoints + 1][numpoints] = 0.0;
 420     K[numpoints + 1][numpoints + 1] = 0.0;
 421     K[numpoints + 1][numpoints + 2] = 0.0;
 422     K[numpoints + 2][numpoints] = 0.0;
 423     K[numpoints + 2][numpoints + 1] = 0.0;
 424     K[numpoints + 2][numpoints + 2] = 0.0;
 425     indx = gp_alloc(sizeof(indx[0]) * (numpoints + 3), "indexes lu");
 426     /* actually, K is *not* positive definite, but
 427        has only non zero real eigenvalues ->
 428        we can use an lu_decomp safely */
 429     lu_decomp(K, numpoints + 3, indx, &d);
 430     lu_backsubst(K, numpoints + 3, indx, b);
 431 #endif /* THIN_PLATE_SPLINES_GRID */
 432
 433     for (i = 0, x = xmin; i < dgrid3d_col_fineness; i++, x += dx) {
 434         struct coordinate GPHUGE *points;
 435
 436         icrv = iso_alloc(dgrid3d_row_fineness + 1);
 437         icrv->p_count = dgrid3d_row_fineness;
 438         icrv->next = this_plot->iso_crvs;
 439         this_plot->iso_crvs = icrv;
 440         points = icrv->points;
 441
 442         for (j = 0, y = ymin; j < dgrid3d_row_fineness; j++, y += dy, points++) {
 443             z = w = 0.0;
 444
 445             /* as soon as ->type is changed to UNDEFINED, break out of
 446              * two inner loops! */
 447             points->type = INRANGE;
 448
 449 #ifdef THIN_PLATE_SPLINES_GRID
 450             z = b[numpoints];
 451             for (k = 0; k < numpoints; k++) {
 452                 double dx = xx[k] - x, dy = yy[k] - y;
 453                 z = z - b[k] * splines_kernel(sqrt(dx * dx + dy * dy));
 454             }
 455             z = z + b[numpoints + 1] * x + b[numpoints + 2] * y;
 456 #else
 457             for (oicrv = old_iso_crvs; oicrv != NULL; oicrv = oicrv->next) {
 458                 struct coordinate GPHUGE *opoints = oicrv->points;
 459                 for (k = 0; k < oicrv->p_count; k++, opoints++) {
 460                     double dist, dist_x = fabs(opoints->x - x), dist_y = fabs(opoints->y - y);
 461                     switch (dgrid3d_norm_value) {
 462                     case 1:
 463                         dist = dist_x + dist_y;
 464                         break;
 465                     case 2:
 466                         dist = dist_x * dist_x + dist_y * dist_y;
 467                         break;
 468                     case 4:
 469                         dist = dist_x * dist_x + dist_y * dist_y;
 470                         dist *= dist;
 471                         break;
 472                     case 8:
 473                         dist = dist_x * dist_x + dist_y * dist_y;
 474                         dist *= dist;
 475                         dist *= dist;
 476                         break;
 477                     case 16:
 478                         dist = dist_x * dist_x + dist_y * dist_y;
 479                         dist *= dist;
 480                         dist *= dist;
 481                         dist *= dist;
 482                         break;
 483                     default:
 484                         dist = pow(dist_x, (double) dgrid3d_norm_value) +
 485                             pow(dist_y, (double) dgrid3d_norm_value);
 486                         break;
 487                     }
 488
 489                     /* The weight of this point is inverse proportional
 490                      * to the distance.
 491                      */
 492                     if (dist == 0.0) {
 493                         /* HBB 981209: revised flagging as undefined */
 494                         /* Supporting all those infinities on various
 495                          * platforms becomes tiresome, to say the least :-(
 496                          * Let's just return the first z where this happens,
 497                          * unchanged, and be done with this, period. */
 498                         points->type = UNDEFINED;
 499                         z = opoints->z;
 500                         w = 1.0;
 501                         break;  /* out of inner loop */
 502                     } else {
 503                         dist = 1.0 / dist;
 504                         z += opoints->z * dist;
 505                         w += dist;
 506                     }
 507                 }
 508                 if (points->type != INRANGE)
 509                     break;      /* out of the second-inner loop as well ... */
 510             }
 511 #endif /* THIN_PLATE_SPLINES_GRID */
 512
 513             /* Now that we've escaped the loops safely, we know that we
 514              * do have a good value in z and w, so we can proceed just as
 515              * if nothing had happened at all. Nice, isn't it? */
 516             points->type = INRANGE;
 517
 518             /* HBB 20010424: if log x or log y axis, we don't want to
 519              * log() the value again --> just store it, and trust that
 520              * it's always inrange */
 521             points->x = x;
 522             points->y = y;
 523
 524 #ifndef THIN_PLATE_SPLINES_GRID
 525             z = z / w;
 526 #endif
 527
 528             STORE_WITH_LOG_AND_UPDATE_RANGE(points->z, z, points->type, z_axis, NOOP, continue);
 529
 530             if (this_plot->pm3d_color_from_column)
 531                 int_error(NO_CARET, "Gridding of the color column is not implemented");
 532             else {
 533                 COLOR_STORE_WITH_LOG_AND_UPDATE_RANGE(points->CRD_COLOR, z, points->type, COLOR_AXIS, NOOP, continue);
 534             }
 535         }
 536     }
 537
 538 #ifdef THIN_PLATE_SPLINES_GRID
 539     free(K);
 540     free(xx);
 541     free(indx);
 542 #endif
 543
 544     /* Delete the old non grid data. */
 545     for (oicrvs = old_iso_crvs; oicrvs != NULL;) {
 546         oicrv = oicrvs;
 547         oicrvs = oicrvs->next;
 548         iso_free(oicrv);
 549     }
 550 }
 551
 552 /* Get 3D data from file, and store into this_plot data
 553  * structure. Takes care of 'set mapping' and 'set dgrid3d'.
 554  *
 555  * Notice: this_plot->token is end of datafile spec, before title etc
 556  * will be moved past title etc after we return */
 557 static void
 558 get_3ddata(struct surface_points *this_plot)
 559 {
 560     int xdatum = 0;
 561     int ydatum = 0;
 562     int j;
 563     double v[MAXDATACOLS];
 564     int pt_in_iso_crv = 0;
 565     struct iso_curve *this_iso;
 566
 567     if (mapping3d == MAP3D_CARTESIAN) {
 568         /* do this check only, if we have PM3D / PM3D-COLUMN not compiled in */
 569         if (df_no_use_specs == 2)
 570             int_error(this_plot->token, "Need 1 or 3 columns for cartesian data");
 571         /* HBB NEW 20060427: if there's only one, explicit using
 572          * column, it's z data.  df_axis[] has to reflect that, so
 573          * df_readline() will expect time/date input. */
 574         if (df_no_use_specs == 1)
 575             df_axis[0] = FIRST_Z_AXIS;
 576     } else {
 577         if (df_no_use_specs == 1)
 578             int_error(this_plot->token, "Need 2 or 3 columns for polar data");
 579     }
 580
 581     this_plot->num_iso_read = 0;
 582     this_plot->has_grid_topology = TRUE;
 583     this_plot->pm3d_color_from_column = FALSE;
 584
 585     /* we ought to keep old memory - most likely case
 586      * is a replot, so it will probably exactly fit into
 587      * memory already allocated ?
 588      */
 589     if (this_plot->iso_crvs != NULL) {
 590         struct iso_curve *icrv, *icrvs = this_plot->iso_crvs;
 591
 592         while (icrvs) {
 593             icrv = icrvs;
 594             icrvs = icrvs->next;
 595             iso_free(icrv);
 596         }
 597         this_plot->iso_crvs = NULL;
 598     }
 599     /* data file is already open */
 600
 601 #ifndef BINARY_DATA_FILE /* NO LONGER REQUIRED FOR GENERAL BINARY OR MATRIX BINARY DATA */
 602     if (df_matrix)
 603         xdatum = df_3dmatrix(this_plot, NEED_PALETTE(this_plot));
 604     else {
 605 #else
 606     if (df_matrix) {
 607         this_plot->has_grid_topology = TRUE;
 608     }
 609
 610     {
 611 #endif
 612         /*{{{  read surface from text file */
 613         struct iso_curve *local_this_iso = iso_alloc(samples_1);
 614         struct coordinate GPHUGE *cp;
 615         struct coordinate GPHUGE *cptail = NULL; /* Only for VECTOR plots */
 616         double x, y, z;
 617         double xtail, ytail, ztail;
 618         double color = VERYLARGE;
 619         int pm3d_color_from_column = FALSE;
 620 #define color_from_column(x) pm3d_color_from_column = x
 621
 622 #ifdef EAM_DATASTRINGS
 623         if (this_plot->plot_style == LABELPOINTS)
 624             expect_string( 4 );
 625 #endif
 626
 627         if (this_plot->plot_style == VECTOR) {
 628             local_this_iso->next = iso_alloc(samples_1);
 629             local_this_iso->next->p_count = 0;
 630         }
 631
 632         while ((j = df_readline(v,MAXDATACOLS)) != DF_EOF) {
 633             if (j == DF_SECOND_BLANK)
 634                 break;          /* two blank lines */
 635             if (j == DF_FIRST_BLANK) {
 636
 637 #if defined(WITH_IMAGE) && defined(BINARY_DATA_FILE)
 638                 /* Images are in a sense similar to isocurves.
 639                  * However, the routine for images is written to
 640                  * compute the two dimensions of coordinates by
 641                  * examining the data alone.  That way it can be used
 642                  * in the 2D plots, for which there is no isoline
 643                  * record.  So, toss out isoline information for
 644                  * images.
 645                  */
 646                 if ((this_plot->plot_style == IMAGE)
 647                     || (this_plot->plot_style == RGBIMAGE))
 648                     continue;
 649 #endif
 650                 if (this_plot->plot_style == VECTOR)
 651                     continue;
 652
 653                 /* one blank line */
 654                 if (pt_in_iso_crv == 0) {
 655                     if (xdatum == 0)
 656                         continue;
 657                     pt_in_iso_crv = xdatum;
 658                 }
 659                 if (xdatum > 0) {
 660                     local_this_iso->p_count = xdatum;
 661                     local_this_iso->next = this_plot->iso_crvs;
 662                     this_plot->iso_crvs = local_this_iso;
 663                     this_plot->num_iso_read++;
 664
 665                     if (xdatum != pt_in_iso_crv)
 666                         this_plot->has_grid_topology = FALSE;
 667
 668                     local_this_iso = iso_alloc(pt_in_iso_crv);
 669                     xdatum = 0;
 670                     ydatum++;
 671                 }
 672                 continue;
 673             }
 674
 675             /* its a data point or undefined */
 676             if (xdatum >= local_this_iso->p_max) {
 677                 /* overflow about to occur. Extend size of points[]
 678                  * array. Double the size, and add 1000 points, to
 679                  * avoid needlessly small steps. */
 680                 iso_extend(local_this_iso, xdatum + xdatum + 1000);
 681                 if (this_plot->plot_style == VECTOR) {
 682                     iso_extend(local_this_iso->next, xdatum + xdatum + 1000);
 683                     local_this_iso->next->p_count = 0;
 684                 }
 685             }
 686             cp = local_this_iso->points + xdatum;
 687
 688             if (this_plot->plot_style == VECTOR) {
 689                 if (j < 6) {
 690                     cp->type = UNDEFINED;
 691                     continue;
 692                 }
 693                 cptail = local_this_iso->next->points + xdatum;
 694             }
 695
 696             if (j == DF_UNDEFINED || j == DF_MISSING) {
 697                 cp->type = UNDEFINED;
 698                 goto come_here_if_undefined;
 699             }
 700             cp->type = INRANGE; /* unless we find out different */
 701
 702             /* EAM Oct 2004 - Substantially rework this section */
 703             /* now that there are many more plot types.         */
 704
 705             x = y = z = 0.0;
 706             xtail = ytail = ztail = 0.0;
 707
 708             /* The x, y, z coordinates depend on the mapping type */
 709             switch (mapping3d) {
 710
 711             case MAP3D_CARTESIAN:
 712                 if (j == 1) {
 713                     x = xdatum;
 714                     y = ydatum;
 715                     z = v[0];
 716                     j = 3;
 717                     break;
 718                 }
 719
 720                 if (j == 2) {
 721                     if (PM3DSURFACE != this_plot->plot_style)
 722                         int_error(this_plot->token,
 723                                   "2 columns only possible with explicit pm3d style (line %d)",
 724                                   df_line_number);
 725                     x = xdatum;
 726                     y = ydatum;
 727                     z = v[0];
 728                     color_from_column(TRUE);
 729                     color = v[1];
 730                     j = 3;
 731                     break;
 732                 }
 733
 734                 /* Assume everybody agrees that x,y,z are the first three specs */
 735                 if (j >= 3) {
 736                     x = v[0];
 737                     y = v[1];
 738                     z = v[2];
 739                     break;
 740                 }
 741
 742                 break;
 743
 744             case MAP3D_SPHERICAL:
 745                 if (j < 2)
 746                     int_error(this_plot->token, "Need 2 or 3 columns");
 747                 if (j < 3) {
 748                     v[2] = 1;   /* default radius */
 749                     j = 3;
 750                 }
 751
 752                 /* Convert to radians. */
 753                 v[0] *= ang2rad;
 754                 v[1] *= ang2rad;
 755
 756                 x = v[2] * cos(v[0]) * cos(v[1]);
 757                 y = v[2] * sin(v[0]) * cos(v[1]);
 758                 z = v[2] * sin(v[1]);
 759
 760                 break;
 761
 762             case MAP3D_CYLINDRICAL:
 763                 if (j < 2)
 764                     int_error(this_plot->token, "Need 2 or 3 columns");
 765                 if (j < 3) {
 766                     v[2] = 1;   /* default radius */
 767                     j = 3;
 768                 }
 769
 770                 /* Convert to radians. */
 771                 v[0] *= ang2rad;
 772
 773                 x = v[2] * cos(v[0]);
 774                 y = v[2] * sin(v[0]);
 775                 z = v[1];
 776
 777                 break;
 778
 779             default:
 780                 int_error(NO_CARET, "Internal error: Unknown mapping type");
 781                 return;
 782             }
 783
 784             if (j < df_no_use_specs)
 785                 int_error(this_plot->token,
 786                         "Wrong number of columns in input data - line %d",
 787                         df_line_number);
 788
 789             /* After the first three columns it gets messy because */
 790             /* different plot styles assume different contents in the columns */
 791             if (j >= 4) {
 792                 if (( this_plot->plot_style == POINTSTYLE
 793                    || this_plot->plot_style == LINESPOINTS)
 794                 &&  this_plot->lp_properties.p_size == PTSZ_VARIABLE) {
 795                     cp->CRD_PTSIZE = v[3];
 796                     color = z;
 797                     color_from_column(FALSE);
 798                 }
 799 #ifdef EAM_DATASTRINGS
 800                 else if (this_plot->plot_style == LABELPOINTS) {
 801                 /* 4th column holds label text rather than color */
 802                 /* text = df_tokens[3]; */
 803                     color = z;
 804                     color_from_column(FALSE);
 805                 }
 806 #endif
 807                 else {
 808                     color = v[3];
 809                     color_from_column(TRUE);
 810                 }
 811             }
 812
 813             if (j >= 5) {
 814                 if ((this_plot->plot_style == POINTSTYLE
 815                    || this_plot->plot_style == LINESPOINTS)
 816                 &&  this_plot->lp_properties.p_size == PTSZ_VARIABLE) {
 817                     color = v[4];
 818                     color_from_column(TRUE);
 819                 }
 820 #ifdef EAM_DATASTRINGS
 821                 if (this_plot->plot_style == LABELPOINTS) {
 822                     /* take color from an explicitly given 5th column */
 823                     color = v[4];
 824                     color_from_column(TRUE);
 825                 }
 826 #endif
 827             }
 828
 829             if (j >= 6) {
 830                 if (this_plot->plot_style == VECTOR) {
 831                     xtail = x + v[3];
 832                     ytail = y + v[4];
 833                     ztail = z + v[5];
 834                     color_from_column(FALSE);
 835                 }
 836             }
 837 #undef color_from_column
 838
 839
 840             /* Adjust for logscales. Set min/max and point types. Store in cp.
 841              * The macro cannot use continue, as it is wrapped in a loop.
 842              * I regard this as correct goto use
 843              */
 844             cp->type = INRANGE;
 845             STORE_WITH_LOG_AND_UPDATE_RANGE(cp->x, x, cp->type, x_axis, NOOP, goto come_here_if_undefined);
 846             STORE_WITH_LOG_AND_UPDATE_RANGE(cp->y, y, cp->type, y_axis, NOOP, goto come_here_if_undefined);
 847             if (this_plot->plot_style == VECTOR) {
 848                 cptail->type = INRANGE;
 849                 STORE_WITH_LOG_AND_UPDATE_RANGE(cptail->x, xtail, cp->type, x_axis, NOOP, goto come_here_if_undefined);
 850                 STORE_WITH_LOG_AND_UPDATE_RANGE(cptail->y, ytail, cp->type, y_axis, NOOP, goto come_here_if_undefined);
 851             }
 852
 853             if (dgrid3d) {
 854                 /* HBB 20010424: in dgrid3d mode, delay log() taking
 855                  * and scaling until after the dgrid process. Only for
 856                  * z, not for x and y, so we can layout the newly
 857                  * created created grid more easily. */
 858                 cp->z = z;
 859                 if (this_plot->plot_style == VECTOR)
 860                     cptail->z = ztail;
 861             } else {
 862                 STORE_WITH_LOG_AND_UPDATE_RANGE(cp->z, z, cp->type, z_axis, NOOP, goto come_here_if_undefined);
 863                 if (this_plot->plot_style == VECTOR)
 864                     STORE_WITH_LOG_AND_UPDATE_RANGE(cptail->z, ztail, cp->type, z_axis, NOOP, goto come_here_if_undefined);
 865
 866                 if (NEED_PALETTE(this_plot)) {
 867                     if (pm3d_color_from_column) {
 868                         COLOR_STORE_WITH_LOG_AND_UPDATE_RANGE(cp->CRD_COLOR, color, cp->type, COLOR_AXIS, NOOP, goto come_here_if_undefined);
 869                     } else {
 870                         COLOR_STORE_WITH_LOG_AND_UPDATE_RANGE(cp->CRD_COLOR, z, cp->type, COLOR_AXIS, NOOP, goto come_here_if_undefined);
 871                     }
 872                 }
 873             }
 874
 875 #ifdef EAM_DATASTRINGS
 876             /* At this point we have stored the point coordinates. Now we need to copy */
 877             /* x,y,z into the text_label structure and add the actual text string.     */
 878             if (this_plot->plot_style == LABELPOINTS)
 879                 store_label(this_plot->labels, cp, xdatum, df_tokens[3], color);
 880 #endif
 881
 882         come_here_if_undefined:
 883             /* some may complain, but I regard this as the correct use of goto */
 884             ++xdatum;
 885         }                       /* end of whileloop - end of surface */
 886
 887         if (pm3d_color_from_column) {
 888             this_plot->pm3d_color_from_column = pm3d_color_from_column;
 889         }
 890
 891         if (xdatum > 0) {
 892             this_plot->num_iso_read++;  /* Update last iso. */
 893             local_this_iso->p_count = xdatum;
 894
 895             /* If this is a VECTOR plot then iso->next is already */
 896             /* occupied by the vector tail coordinates.           */
 897             if (this_plot->plot_style != VECTOR)
 898                 local_this_iso->next = this_plot->iso_crvs;
 899             this_plot->iso_crvs = local_this_iso;
 900
 901             if (xdatum != pt_in_iso_crv)
 902                 this_plot->has_grid_topology = FALSE;
 903
 904         } else { /* Free last allocation */
 905             if (this_plot->plot_style == VECTOR)
 906                 iso_free(local_this_iso->next);
 907             iso_free(local_this_iso);
 908         }
 909
 910         /*}}} */
 911     }
 912
 913     if (dgrid3d && this_plot->num_iso_read > 0)
 914         grid_nongrid_data(this_plot);
 915
 916     /* This check used to be done in graph3d */
 917         if (X_AXIS.min == VERYLARGE || X_AXIS.max == -VERYLARGE ||
 918             Y_AXIS.min == VERYLARGE || Y_AXIS.max == -VERYLARGE ||
 919             Z_AXIS.min == VERYLARGE || Z_AXIS.max == -VERYLARGE)
 920                 int_warn(NO_CARET,
 921                     "Axis range undefined due to improper data values. NaN? Inf?");
 922
 923     if (this_plot->num_iso_read <= 1)
 924         this_plot->has_grid_topology = FALSE;
 925     if (this_plot->has_grid_topology && !hidden3d) {
 926         struct iso_curve *new_icrvs = NULL;
 927         int num_new_iso = this_plot->iso_crvs->p_count, len_new_iso = this_plot->num_iso_read;
 928         int i;
 929
 930         /* Now we need to set the other direction (pseudo) isolines. */
 931         for (i = 0; i < num_new_iso; i++) {
 932             struct iso_curve *new_icrv = iso_alloc(len_new_iso);
 933
 934             new_icrv->p_count = len_new_iso;
 935
 936             for (j = 0, this_iso = this_plot->iso_crvs;
 937                  this_iso != NULL;
 938                  j++, this_iso = this_iso->next) {
 939                 /* copy whole point struct to get type too.
 940                  * wasteful for windows, with padding */
 941                 /* more efficient would be extra pointer to same struct */
 942                 new_icrv->points[j] = this_iso->points[i];
 943             }
 944
 945             new_icrv->next = new_icrvs;
 946             new_icrvs = new_icrv;
 947         }
 948
 949         /* Append the new iso curves after the read ones. */
 950         for (this_iso = this_plot->iso_crvs;
 951              this_iso->next != NULL;
 952              this_iso = this_iso->next);
 953         this_iso->next = new_icrvs;
 954     }
 955 }
 956
 957 static void
 958 print_3dtable(int pcount)
 959 {
 960     struct surface_points *this_plot;
 961     int i, surface;
 962     struct coordinate GPHUGE *point;
 963     char *buffer = gp_alloc(150, "print_3dtable output buffer");
 964     FILE *outfile = (table_outfile) ? table_outfile : gpoutfile;
 965
 966     for (surface = 0, this_plot = first_3dplot;
 967          surface < pcount;
 968          this_plot = this_plot->next_sp, surface++) {
 969         fprintf(outfile, "\n#Surface %d of %d surfaces\n", surface, pcount);
 970
 971         if (draw_surface) {
 972             struct iso_curve *icrvs;
 973             int curve;
 974
 975             /* only the curves in one direction */
 976             for (curve = 0, icrvs = this_plot->iso_crvs;
 977                  icrvs && curve < this_plot->num_iso_read;
 978                  icrvs = icrvs->next, curve++) {
 979                 fprintf(outfile, "\n#IsoCurve %d, %d points\n#x y z type\n",
 980                         curve, icrvs->p_count);
 981                 for (i = 0, point = icrvs->points;
 982                      i < icrvs->p_count;
 983                      i++, point++) {
 984                     /* HBB 20020405: new macro to use the 'set format'
 985                      * in their full flexibility */
 986 #define OUTPUT_NUMBER(field, axis)                                      \
 987                     gprintf(buffer, 150, axis_array[axis].formatstring, \
 988                             1.0, point->field);                         \
 989                     fputs(buffer, outfile);                             \
 990                     fputc(' ', outfile);
 991                     OUTPUT_NUMBER(x, FIRST_X_AXIS);
 992                     OUTPUT_NUMBER(y, FIRST_Y_AXIS);
 993                     OUTPUT_NUMBER(z, FIRST_Z_AXIS);
 994                     fprintf(outfile, "%c\n",
 995                             point->type == INRANGE
 996                             ? 'i' : point->type == OUTRANGE
 997                             ? 'o' : 'u');
 998                 } /* for(point) */
 999             } /* for(icrvs) */
1000             putc('\n', outfile);
1001         } /* if(draw_surface) */
1002
1003         if (draw_contour) {
1004             int number = 0;
1005             struct gnuplot_contours *c = this_plot->contours;
1006
1007             while (c) {
1008                 int count = c->num_pts;
1009                 struct coordinate GPHUGE *point = c->coords;
1010
1011                 if (c->isNewLevel)
1012                     /* don't display count - contour split across chunks */
1013                     /* put # in case user wants to use it for a plot */
1014                     /* double blank line to allow plot ... index ... */
1015                     fprintf(outfile, "\n# Contour %d, label: %s\n",
1016                             number++, c->label);
1017
1018                 for (; --count >= 0; ++point) {
1019                     OUTPUT_NUMBER(x, FIRST_X_AXIS);
1020                     OUTPUT_NUMBER(y, FIRST_Y_AXIS);
1021                     OUTPUT_NUMBER(z, FIRST_Z_AXIS);
1022                     putc('\n', outfile);
1023                 }
1024
1025                 /* blank line between segments of same contour */
1026                 putc('\n', outfile);
1027                 c = c->next;
1028 #undef OUTPUT_NUMBER
1029             } /* while (contour) */
1030         } /* if (draw_contour) */
1031     } /* for(surface) */
1032     fflush(outfile);
1033
1034     free(buffer);
1035 }
1036
1037 /* HBB 20000501: code isolated from eval_3dplots(), where practically
1038  * identical code occured twice, for direct and crossing isolines,
1039  * respectively.  The latter only are done for in non-hidden3d
1040  * mode. */
1041 static void
1042 calculate_set_of_isolines(
1043     AXIS_INDEX value_axis,
1044     TBOOLEAN cross,
1045     struct iso_curve **this_iso,
1046     AXIS_INDEX iso_axis,
1047     double iso_min, double iso_step,
1048     int num_iso_to_use,
1049     AXIS_INDEX sam_axis,
1050     double sam_min, double sam_step,
1051     int num_sam_to_use,
1052     TBOOLEAN need_palette)
1053 {
1054     int i, j;
1055     struct coordinate GPHUGE *points = (*this_iso)->points;
1056     int do_update_color = need_palette && (!parametric || (parametric && value_axis == FIRST_Z_AXIS));
1057
1058     for (j = 0; j < num_iso_to_use; j++) {
1059         double iso = iso_min + j * iso_step;
1060         /* HBB 20000501: with the new code, it should
1061          * be safe to rely on the actual 'v' axis not
1062          * to be improperly logscaled... */
1063         (void) Gcomplex(&plot_func.dummy_values[cross ? 0 : 1],
1064                         AXIS_DE_LOG_VALUE(iso_axis, iso), 0.0);
1065
1066         for (i = 0; i < num_sam_to_use; i++) {
1067             double sam = sam_min + i * sam_step;
1068             struct value a;
1069             double temp;
1070
1071             (void) Gcomplex(&plot_func.dummy_values[cross ? 1 : 0],
1072                             AXIS_DE_LOG_VALUE(sam_axis, sam), 0.0);
1073
1074             if (cross) {
1075                 points[i].x = iso;
1076                 points[i].y = sam;
1077             } else {
1078                 points[i].x = sam;
1079                 points[i].y = iso;
1080             }
1081
1082             evaluate_at(plot_func.at, &a);
1083
1084             if (undefined || (fabs(imag(&a)) > zero)) {
1085                 points[i].type = UNDEFINED;
1086                 continue;
1087             }
1088
1089             temp = real(&a);
1090             points[i].type = INRANGE;
1091             STORE_WITH_LOG_AND_UPDATE_RANGE(points[i].z, temp, points[i].type,
1092                                             value_axis, NOOP, NOOP);
1093             if (do_update_color) {
1094                 COLOR_STORE_WITH_LOG_AND_UPDATE_RANGE(points[i].CRD_COLOR, temp, points[i].type, COLOR_AXIS, NOOP, NOOP);
1095             }
1096         }
1097         (*this_iso)->p_count = num_sam_to_use;
1098         *this_iso = (*this_iso)->next;
1099         points = (*this_iso) ? (*this_iso)->points : NULL;
1100     }
1101 }
1102
1103
1104 /*
1105  * This parses the splot command after any range specifications. To support
1106  * autoscaling on the x/z axis, we want any data files to define the x/y
1107  * range, then to plot any functions using that range. We thus parse the
1108  * input twice, once to pick up the data files, and again to pick up the
1109  * functions. Definitions are processed twice, but that won't hurt.
1110  * div - okay, it doesn't hurt, but every time an option as added for
1111  * datafiles, code to parse it has to be added here. Change so that
1112  * we store starting-token in the plot structure.
1113  */
1114 static void
1115 eval_3dplots()
1116 {
1117     int i;
1118     struct surface_points **tp_3d_ptr;
1119     int start_token, end_token;
1120     int begin_token;
1121     TBOOLEAN some_data_files = FALSE, some_functions = FALSE;
1122     int plot_num, line_num, point_num;
1123     /* part number of parametric function triplet: 0 = z, 1 = y, 2 = x */
1124     int crnt_param = 0;
1125     char *xtitle;
1126     char *ytitle;
1127     legend_key *key = &keyT;
1128
1129     /* Reset first_3dplot. This is usually done at the end of this function.
1130      * If there is an error within this function, the memory is left allocated,
1131      * since we cannot call sp_free if the list is incomplete
1132      */
1133     if (first_3dplot && plot3d_num>0)
1134       sp_free(first_3dplot);
1135     plot3d_num=0;
1136     first_3dplot = NULL;
1137
1138     x_axis = FIRST_X_AXIS;
1139     y_axis = FIRST_Y_AXIS;
1140     z_axis = FIRST_Z_AXIS;
1141
1142     tp_3d_ptr = &(first_3dplot);
1143     plot_num = 0;
1144     line_num = 0;               /* default line type */
1145     point_num = 0;              /* default point type */
1146
1147     xtitle = NULL;
1148     ytitle = NULL;
1149
1150     begin_token = c_token;
1151
1152 /*** First Pass: Read through data files ***/
1153     /*
1154      * This pass serves to set the x/yranges and to parse the command, as
1155      * well as filling in every thing except the function data. That is done
1156      * after the x/yrange is defined.
1157      */
1158     while (TRUE) {
1159         if (END_OF_COMMAND)
1160             int_error(c_token, "function to plot expected");
1161
1162         start_token = c_token;
1163
1164         if (is_definition(c_token)) {
1165             define();
1166         } else {
1167             int specs = -1;
1168             struct surface_points *this_plot;
1169
1170             char *name_str;
1171             TBOOLEAN duplication = FALSE;
1172             TBOOLEAN set_title = FALSE, set_with = FALSE;
1173             TBOOLEAN set_lpstyle = FALSE;
1174             TBOOLEAN checked_once = FALSE;
1175 #ifdef EAM_DATASTRINGS
1176             TBOOLEAN set_labelstyle = FALSE;
1177 #endif
1178
1179             dummy_func = &plot_func;
1180             /* WARNING: do NOT free name_str */
1181             /* FIXME: could also be saved in this_plot */
1182             name_str = string_or_express(NULL);
1183             dummy_func = NULL;
1184             if (name_str) {
1185                 /*{{{  data file to plot */
1186                 if (parametric && crnt_param != 0)
1187                     int_error(c_token, "previous parametric function not fully specified");
1188
1189                 if (!some_data_files) {
1190                     if (axis_array[FIRST_X_AXIS].autoscale & AUTOSCALE_MIN) {
1191                         axis_array[FIRST_X_AXIS].min = VERYLARGE;
1192                     }
1193                     if (axis_array[FIRST_X_AXIS].autoscale & AUTOSCALE_MAX) {
1194                         axis_array[FIRST_X_AXIS].max = -VERYLARGE;
1195                     }
1196                     if (axis_array[FIRST_Y_AXIS].autoscale & AUTOSCALE_MIN) {
1197                         axis_array[FIRST_Y_AXIS].min = VERYLARGE;
1198                     }
1199                     if (axis_array[FIRST_Y_AXIS].autoscale & AUTOSCALE_MAX) {
1200                         axis_array[FIRST_Y_AXIS].max = -VERYLARGE;
1201                     }
1202                     some_data_files = TRUE;
1203                 }
1204                 if (*tp_3d_ptr)
1205                     this_plot = *tp_3d_ptr;
1206                 else {          /* no memory malloc()'d there yet */
1207                     /* Allocate enough isosamples and samples */
1208                     this_plot = sp_alloc(0, 0, 0, 0);
1209                     *tp_3d_ptr = this_plot;
1210                 }
1211
1212                 this_plot->plot_type = DATA3D;
1213                 this_plot->plot_style = data_style;
1214
1215                 df_set_plot_mode(MODE_SPLOT);
1216                 specs = df_open(name_str, MAXDATACOLS);
1217 #ifdef BINARY_DATA_FILE
1218                 if (df_matrix)
1219                     this_plot->has_grid_topology = TRUE;
1220 #endif
1221
1222                 /* parses all datafile-specific modifiers */
1223                 /* we will load the data after parsing title,with,... */
1224
1225                 /* for capture to key */
1226                 this_plot->token = end_token = c_token - 1;
1227
1228                 /* this_plot->token is temporary, for errors in get_3ddata() */
1229
1230                 if (specs < 3) {
1231                     if (axis_array[FIRST_X_AXIS].is_timedata) {
1232                         int_error(c_token, "Need full using spec for x time data");
1233                     }
1234                     if (axis_array[FIRST_Y_AXIS].is_timedata) {
1235                         int_error(c_token, "Need full using spec for y time data");
1236                     }
1237                 }
1238                 df_axis[0] = FIRST_X_AXIS;
1239                 df_axis[1] = FIRST_Y_AXIS;
1240                 df_axis[2] = FIRST_Z_AXIS;
1241
1242                 /*}}} */
1243
1244             } else {            /* function to plot */
1245
1246                 /*{{{  function */
1247                 ++plot_num;
1248                 if (parametric) {
1249                     /* Rotate between x/y/z axes */
1250                     /* +2 same as -1, but beats -ve problem */
1251                     crnt_param = (crnt_param + 2) % 3;
1252                 }
1253                 if (*tp_3d_ptr) {
1254                     this_plot = *tp_3d_ptr;
1255                     if (!hidden3d)
1256                         sp_replace(this_plot, samples_1, iso_samples_1,
1257                                    samples_2, iso_samples_2);
1258                     else
1259                         sp_replace(this_plot, iso_samples_1, 0,
1260                                    0, iso_samples_2);
1261
1262                 } else {        /* no memory malloc()'d there yet */
1263                     /* Allocate enough isosamples and samples */
1264                     if (!hidden3d)
1265                         this_plot = sp_alloc(samples_1, iso_samples_1,
1266                                              samples_2, iso_samples_2);
1267                     else
1268                         this_plot = sp_alloc(iso_samples_1, 0,
1269                                              0, iso_samples_2);
1270                     *tp_3d_ptr = this_plot;
1271                 }
1272
1273                 this_plot->plot_type = FUNC3D;
1274                 this_plot->has_grid_topology = TRUE;
1275                 this_plot->plot_style = func_style;
1276                 this_plot->num_iso_read = iso_samples_2;
1277                 /* ignore it for now */
1278                 some_functions = TRUE;
1279                 end_token = c_token - 1;
1280                 /*}}} */
1281
1282             }                   /* end of IS THIS A FILE OR A FUNC block */
1283
1284             /* clear current title, if exist */
1285             if (this_plot->title) {
1286                 free(this_plot->title);
1287                 this_plot->title = NULL;
1288             }
1289
1290             /* default line and point types, no palette */
1291             this_plot->lp_properties.use_palette = 0;
1292             this_plot->lp_properties.l_type = line_num;
1293             this_plot->lp_properties.p_type = point_num;
1294
1295             /* user may prefer explicit line styles */
1296             if (prefer_line_styles)
1297                 lp_use_properties(&this_plot->lp_properties, line_num+1, TRUE);
1298
1299             /* pm 25.11.2001 allow any order of options */
1300             while (!END_OF_COMMAND || !checked_once) {
1301
1302                 /* deal with title */
1303                 if (almost_equals(c_token, "t$itle")) {
1304                     if (set_title) {
1305                         duplication=TRUE;
1306                         break;
1307                     }
1308                     this_plot->title_no_enhanced = !key->enhanced;
1309                         /* title can be enhanced if not explicitly disabled */
1310                     if (parametric) {
1311                         if (crnt_param != 0)
1312                             int_error(c_token, "\"title\" allowed only after parametric function fully specified");
1313                         else {
1314                             if (xtitle != NULL)
1315                                 xtitle[0] = NUL;        /* Remove default title . */
1316                             if (ytitle != NULL)
1317                                 ytitle[0] = NUL;        /* Remove default title . */
1318                         }
1319                     }
1320                     c_token++;
1321                     if (!(this_plot->title = try_to_get_string()))
1322                         int_error(c_token, "expecting \"title\" for plot");
1323                     set_title = TRUE;
1324                     continue;
1325                 }
1326
1327                 if (almost_equals(c_token, "not$itle")) {
1328                     if (set_title) {
1329                         duplication=TRUE;
1330                         break;
1331                     }
1332                     c_token++;
1333                     if (isstringvalue(c_token))
1334                         try_to_get_string(); /* ignore optionally given title string */
1335                     if (xtitle != NULL)
1336                         xtitle[0] = '\0';
1337                     if (ytitle != NULL)
1338                         ytitle[0] = '\0';
1339                     /*   this_plot->title = NULL;   */
1340                     set_title = TRUE;
1341                     continue;
1342                 }
1343
1344                 /* deal with style */
1345                 if (almost_equals(c_token, "w$ith")) {
1346                     if (set_with) {
1347                         duplication=TRUE;
1348                         break;
1349                     }
1350                     this_plot->plot_style = get_style();
1351                     if ((this_plot->plot_type == FUNC3D) &&
1352                         ((this_plot->plot_style & PLOT_STYLE_HAS_ERRORBAR)
1353 #ifdef EAM_DATASTRINGS
1354                         || (this_plot->plot_style == LABELPOINTS)
1355 #endif
1356                         )) {
1357                         int_warn(c_token, "This plot style is only for datafiles , reverting to \"points\"");
1358                         this_plot->plot_style = POINTSTYLE;
1359                     }
1360
1361                     if ((this_plot->plot_style | data_style) & PM3DSURFACE) {
1362                         if (equals(c_token, "at")) {
1363                         /* option 'with pm3d [at ...]' is explicitly specified */
1364                         c_token++;
1365                         if (get_pm3d_at_option(&this_plot->pm3d_where[0]))
1366                             return; /* error */
1367                         }
1368                     }
1369
1370                     set_with = TRUE;
1371                     continue;
1372                 }
1373
1374                 /* EAM Dec 2005 - Hidden3D code by default includes points,  */
1375                 /* labels and vectors in the hidden3d processing. Check here */
1376                 /* if this particular plot wants to be excluded.             */
1377                 if (almost_equals(c_token, "nohidden$3d")) {
1378                     c_token++;
1379                     this_plot->opt_out_of_hidden3d = TRUE;
1380                     continue;
1381                 }
1382
1383                 /* "set contour" is global. Allow individual plots to opt out */
1384                 if (almost_equals(c_token, "nocon$tours")) {
1385                     c_token++;
1386                     this_plot->opt_out_of_contours = TRUE;
1387                     continue;
1388                 }
1389
1390 #ifdef EAM_DATASTRINGS
1391                 /* Labels can have font and text property info as plot options */
1392                 /* In any case we must allocate one instance of the text style */
1393                 /* that all labels in the plot will share.                     */
1394                 if (this_plot->plot_style == LABELPOINTS)  {
1395                     int stored_token = c_token;
1396                     if (this_plot->labels == NULL) {
1397                         this_plot->labels = new_text_label(-1);
1398                         this_plot->labels->pos = JUST_CENTRE;
1399                         this_plot->labels->layer = LAYER_PLOTLABELS;
1400                     }
1401                     parse_label_options(this_plot->labels);
1402                     checked_once = TRUE;
1403                     if (stored_token != c_token) {
1404                         if (set_labelstyle) {
1405                             duplication = TRUE;
1406                             break;
1407                         } else {
1408                             set_labelstyle = TRUE;
1409                             continue;
1410                         }
1411                     }
1412                 }
1413 #endif
1414
1415                 /* pick up line/point specs
1416                  * - point spec allowed if style uses points, ie style&2 != 0
1417                  * - keywords are optional
1418                  */
1419                 if (this_plot->plot_style == VECTOR) {
1420                     int stored_token = c_token;
1421
1422                     if (!checked_once) {
1423                         default_arrow_style(&this_plot->arrow_properties);
1424                         this_plot->arrow_properties.lp_properties.l_type = line_num;
1425                         checked_once = TRUE;
1426                     }
1427                     arrow_parse(&this_plot->arrow_properties, TRUE);
1428                     if (stored_token != c_token) {
1429                          if (set_lpstyle) {
1430                             duplication = TRUE;
1431                             break;
1432                          } else {
1433                             set_lpstyle = TRUE;
1434                             this_plot->lp_properties = this_plot->arrow_properties.lp_properties;
1435                             continue;
1436                         }
1437                     }
1438                 } else {
1439                     int stored_token = c_token;
1440                     struct lp_style_type lp = DEFAULT_LP_STYLE_TYPE;
1441
1442                     lp.l_type = line_num;
1443                     lp.p_type = point_num;
1444
1445                     /* user may prefer explicit line styles */
1446                     if (prefer_line_styles)
1447                         lp_use_properties(&lp, line_num+1, TRUE);
1448
1449                     lp_parse(&lp, TRUE,
1450                              this_plot->plot_style & PLOT_STYLE_HAS_POINT);
1451                     checked_once = TRUE;
1452                     if (stored_token != c_token) {
1453                         if (set_lpstyle) {
1454                             duplication=TRUE;
1455                             break;
1456                         } else {
1457                             this_plot->lp_properties = lp;
1458                             set_lpstyle = TRUE;
1459                             continue;
1460                         }
1461                     }
1462                 }
1463
1464                 break; /* unknown option */
1465
1466             }  /* while (!END_OF_COMMAND)*/
1467
1468             if (duplication)
1469                 int_error(c_token, "duplicated or contradicting arguments in plot options");
1470
1471             /* set default values for title if this has not been specified */
1472             if (!set_title) {
1473                 this_plot->title_no_enhanced = TRUE; /* filename or function cannot be enhanced */
1474                 if (key->auto_titles == FILENAME_KEYTITLES) {
1475 #ifdef BINARY_DATA_FILE
1476                     /* DJS (20 Aug 2004) I'd prefer that the df_binary flag be discarded.  There
1477                      * is nothing special about the file being binary that its title should be
1478                      * different.  Can't the decision to do this be based on some other criteria,
1479                      * like the presence of a nonconventional `using`?
1480                      */
1481 #endif
1482                     if (this_plot->plot_type == DATA3D && df_binary==TRUE && end_token==start_token+1)
1483                         /* let default title for  splot 'a.dat' binary  is 'a.dat'
1484                          * while for  'a.dat' binary using 2:1:3  will be all 4 words */
1485                         m_capture(&(this_plot->title), start_token, start_token);
1486                     else
1487                         m_capture(&(this_plot->title), start_token, end_token);
1488                     if (crnt_param == 2)
1489                         xtitle = this_plot->title;
1490                     else if (crnt_param == 1)
1491                         ytitle = this_plot->title;
1492                 } else {
1493                     if (xtitle != NULL)
1494                         xtitle[0] = '\0';
1495                     if (ytitle != NULL)
1496                         ytitle[0] = '\0';
1497                     /*   this_plot->title = NULL;   */
1498                 }
1499             }
1500
1501             /* No line/point style given. As lp_parse also supplies
1502              * the defaults for linewidth and pointsize, call it now
1503              * to define them. */
1504             if (! set_lpstyle) {
1505                 if (this_plot->plot_style == VECTOR) {
1506                     this_plot->arrow_properties.lp_properties.l_type = line_num;
1507                     arrow_parse(&this_plot->arrow_properties, TRUE);
1508                     this_plot->lp_properties = this_plot->arrow_properties.lp_properties;
1509                 } else {
1510                     this_plot->lp_properties.l_type = line_num;
1511                     this_plot->lp_properties.l_width = 1.0;
1512                     this_plot->lp_properties.p_type = point_num;
1513                     this_plot->lp_properties.p_size = pointsize;
1514                     this_plot->lp_properties.use_palette = 0;
1515
1516                     /* user may prefer explicit line styles */
1517                     if (prefer_line_styles)
1518                         lp_use_properties(&this_plot->lp_properties, line_num+1, TRUE);
1519
1520                     lp_parse(&this_plot->lp_properties, TRUE,
1521                          this_plot->plot_style & PLOT_STYLE_HAS_POINT);
1522                 }
1523
1524 #ifdef BACKWARDS_COMPATIBLE
1525                 /* allow old-style syntax - ignore case lt 3 4 for example */
1526                 if (!END_OF_COMMAND && isanumber(c_token)) {
1527                     struct value t;
1528
1529                     this_plot->lp_properties.l_type =
1530                         this_plot->lp_properties.p_type =
1531                         (int) real(const_express(&t)) - 1;
1532
1533                     if (isanumber(c_token))
1534                         this_plot->lp_properties.p_type =
1535                             (int) real(const_express(&t)) - 1;
1536                 }
1537 #endif
1538
1539             }
1540
1541             /* Rule out incompatible line/point/style options */
1542             if (this_plot->plot_type == FUNC3D) {
1543                 if ((this_plot->plot_style & PLOT_STYLE_HAS_POINT)
1544                 &&  (this_plot->lp_properties.p_size == PTSZ_VARIABLE))
1545                     this_plot->lp_properties.p_size = 1;
1546             }
1547
1548             if (crnt_param == 0
1549                 && this_plot->plot_style != PM3DSURFACE
1550                 /* don't increment the default line/point properties if
1551                  * this_plot is an EXPLICIT pm3d surface plot */
1552 #ifdef WITH_IMAGE
1553                 && this_plot->plot_style != IMAGE
1554                 && this_plot->plot_style != RGBIMAGE
1555                 /* same as above, for an (rgb)image plot */
1556 #endif
1557                 ) {
1558                 if (this_plot->plot_style & PLOT_STYLE_HAS_POINT)
1559                     point_num += 1 + (draw_contour != 0) + (hidden3d != 0);
1560                 line_num += 1 + (draw_contour != 0) + (hidden3d != 0);
1561             }
1562
1563 #ifdef WITH_IMAGE
1564             /* Styles that utilize palettes. */
1565             if (this_plot->plot_style == IMAGE)
1566                 this_plot->lp_properties.use_palette = 1;
1567 #endif
1568
1569             /* now get the data... having to think hard here...
1570              * first time through, we fill in this_plot. For second
1571              * surface in file, we have to allocate another surface
1572              * struct. BUT we may allocate this store only to
1573              * find that it is merely some blank lines at end of file
1574              * tp_3d_ptr is still pointing at next field of prev. plot,
1575              * before :    prev_or_first -> this_plot -> possible_preallocated_store
1576              *                tp_3d_ptr--^
1577              * after  :    prev_or_first -> first -> second -> last -> possibly_more_store
1578              *                                        tp_3d_ptr ----^
1579              * if file is empty, tp_3d_ptr is not moved. this_plot continues
1580              * to point at allocated storage, but that will be reused later
1581              */
1582
1583             assert(this_plot == *tp_3d_ptr);
1584
1585             if (this_plot->plot_type == DATA3D) {
1586                 /*{{{  read data */
1587                 /* remember settings for second surface in file */
1588                 struct lp_style_type *these_props = &(this_plot->lp_properties);
1589                 enum PLOT_STYLE this_style = this_plot->plot_style;
1590                 struct surface_points *first_dataset = this_plot;
1591                     /* pointer to the plot of the first dataset (surface) in the file */
1592                 int this_token = this_plot->token;
1593
1594                 do {
1595                     this_plot = *tp_3d_ptr;
1596                     assert(this_plot != NULL);
1597
1598                     /* dont move tp_3d_ptr until we are sure we
1599                      * have read a surface
1600                      */
1601
1602                     /* used by get_3ddata() */
1603                     this_plot->token = this_token;
1604
1605                     get_3ddata(this_plot);
1606                     /* for second pass */
1607                     this_plot->token = c_token;
1608
1609                     if (this_plot->num_iso_read == 0)
1610                         this_plot->plot_type = NODATA;
1611
1612                     if (this_plot != first_dataset)
1613                         /* copy (explicit) "with pm3d at ..." option from the first dataset in the file */
1614                         strcpy(this_plot->pm3d_where, first_dataset->pm3d_where);
1615
1616                     /* okay, we have read a surface */
1617                     ++plot_num;
1618                     tp_3d_ptr = &(this_plot->next_sp);
1619                     if (df_eof)
1620                         break;
1621
1622                     /* there might be another surface so allocate
1623                      * and prepare another surface structure
1624                      * This does no harm if in fact there are
1625                      * no more surfaces to read
1626                      */
1627
1628                     if ((this_plot = *tp_3d_ptr) != NULL) {
1629                         if (this_plot->title) {
1630                             free(this_plot->title);
1631                             this_plot->title = NULL;
1632                         }
1633                     } else {
1634                         /* Allocate enough isosamples and samples */
1635                         this_plot = *tp_3d_ptr = sp_alloc(0, 0, 0, 0);
1636                     }
1637
1638                     this_plot->plot_type = DATA3D;
1639                     this_plot->plot_style = this_style;
1640                     /* Struct copy */
1641                     this_plot->lp_properties = *these_props;
1642                 } while (!df_eof);
1643
1644                 df_close();
1645                 /*}}} */
1646
1647             } else {            /* not a data file */
1648                 tp_3d_ptr = &(this_plot->next_sp);
1649                 this_plot->token = c_token;     /* store for second pass */
1650             }
1651
1652         }                       /* !is_definition() : end of scope of this_plot */
1653
1654
1655         if (equals(c_token, ","))
1656             c_token++;
1657         else
1658             break;
1659
1660     }                           /* while(TRUE), ie first pass */
1661
1662
1663     if (parametric && crnt_param != 0)
1664         int_error(NO_CARET, "parametric function not fully specified");
1665
1666
1667 /*** Second Pass: Evaluate the functions ***/
1668     /*
1669      * Everything is defined now, except the function data. We expect no
1670      * syntax errors, etc, since the above parsed it all. This makes the code
1671      * below simpler. If axis_array[FIRST_Y_AXIS].autoscale, the yrange may still change.
1672      * - eh ?  - z can still change.  x/y/z can change if we are parametric ??
1673      */
1674
1675     if (some_functions) {
1676         /* I've changed the controlled variable in fn plots to u_min etc since
1677          * it's easier for me to think parametric - 'normal' plot is after all
1678          * a special case. I was confused about x_min being both minimum of
1679          * x values found, and starting value for fn plots.
1680          */
1681         double u_min, u_max, u_step, v_min, v_max, v_step;
1682         double u_isostep, v_isostep;
1683         AXIS_INDEX u_axis, v_axis;
1684         struct surface_points *this_plot;
1685
1686         /* Make these point out the right 'u' and 'v' axis. In
1687          * non-parametric mode, x is used as u, and y as v */
1688         u_axis = parametric ? U_AXIS : FIRST_X_AXIS;
1689         v_axis = parametric ? V_AXIS : FIRST_Y_AXIS;
1690
1691         if (!parametric) {
1692             /*{{{  check ranges */
1693             /* give error if xrange badly set from missing datafile error
1694              * parametric fn can still set ranges
1695              * if there are no fns, we'll report it later as 'nothing to plot'
1696              */
1697
1698             /* check that xmin -> xmax is not too small */
1699             axis_checked_extend_empty_range(FIRST_X_AXIS, "x range is invalid");
1700             axis_checked_extend_empty_range(FIRST_Y_AXIS, "y range is invalid");
1701             /*}}} */
1702         }
1703         if (parametric && !some_data_files) {
1704             /*{{{  set ranges */
1705             /* parametric fn can still change x/y range */
1706             if (axis_array[FIRST_X_AXIS].autoscale & AUTOSCALE_MIN)
1707                 axis_array[FIRST_X_AXIS].min = VERYLARGE;
1708             if (axis_array[FIRST_X_AXIS].autoscale & AUTOSCALE_MAX)
1709                 axis_array[FIRST_X_AXIS].max = -VERYLARGE;
1710             if (axis_array[FIRST_Y_AXIS].autoscale & AUTOSCALE_MIN)
1711                 axis_array[FIRST_Y_AXIS].min = VERYLARGE;
1712             if (axis_array[FIRST_Y_AXIS].autoscale & AUTOSCALE_MAX)
1713                 axis_array[FIRST_Y_AXIS].max = -VERYLARGE;
1714             /*}}} */
1715         }
1716
1717         /*{{{  figure ranges, taking logs etc into account */
1718         u_min = axis_log_value_checked(u_axis, axis_array[u_axis].min, "x range");
1719         u_max = axis_log_value_checked(u_axis, axis_array[u_axis].max, "x range");
1720         v_min = axis_log_value_checked(v_axis, axis_array[v_axis].min, "y range");
1721         v_max = axis_log_value_checked(v_axis, axis_array[v_axis].max, "y range");
1722         /*}}} */
1723
1724
1725         if (samples_1 < 2 || samples_2 < 2 || iso_samples_1 < 2 ||
1726             iso_samples_2 < 2) {
1727             int_error(NO_CARET, "samples or iso_samples < 2. Must be at least 2.");
1728         }
1729
1730         /* start over */
1731         this_plot = first_3dplot;
1732         c_token = begin_token;
1733
1734         /* why do attributes of this_plot matter ? */
1735         /* FIXME HBB 20000501: I think they don't, actually. I'm
1736          * taking out references to has_grid_topology in this part of
1737          * the code.  It only deals with function, which always is
1738          * gridded */
1739
1740         if (hidden3d) {
1741             u_step = (u_max - u_min) / (iso_samples_1 - 1);
1742             v_step = (v_max - v_min) / (iso_samples_2 - 1);
1743         } else {
1744             u_step = (u_max - u_min) / (samples_1 - 1);
1745             v_step = (v_max - v_min) / (samples_2 - 1);
1746         }
1747
1748         u_isostep = (u_max - u_min) / (iso_samples_1 - 1);
1749         v_isostep = (v_max - v_min) / (iso_samples_2 - 1);
1750
1751
1752         /* Read through functions */
1753         while (TRUE) {
1754             if (is_definition(c_token)) {
1755                 define();
1756             } else {
1757                 struct at_type *at_ptr;
1758                 char *name_str;
1759
1760                 dummy_func = &plot_func;
1761                 name_str = string_or_express(&at_ptr);
1762
1763                 if (!name_str) {                /* func to plot */
1764                     /*{{{  evaluate function */
1765                     struct iso_curve *this_iso = this_plot->iso_crvs;
1766                     int num_sam_to_use, num_iso_to_use;
1767
1768                     /* crnt_param is used as the axis number.  As the
1769                      * axis array indices are ordered z, y, x, we have
1770                      * to count *backwards*, starting starting at 2,
1771                      * to properly store away contents to x, y and
1772                      * z. The following little gimmick does that. */
1773                     if (parametric)
1774                         crnt_param = (crnt_param + 2) % 3;
1775
1776                     plot_func.at = at_ptr;
1777
1778                     num_iso_to_use = iso_samples_2;
1779                     num_sam_to_use = hidden3d ? iso_samples_1 : samples_1;
1780
1781                     calculate_set_of_isolines(crnt_param, FALSE, &this_iso,
1782                                               v_axis, v_min, v_isostep,
1783                                               num_iso_to_use,
1784                                               u_axis, u_min, u_step,
1785                                               num_sam_to_use,
1786                                               NEED_PALETTE(this_plot));
1787
1788                     if (!hidden3d) {
1789                         num_iso_to_use = iso_samples_1;
1790                         num_sam_to_use = samples_2;
1791
1792                         calculate_set_of_isolines(crnt_param, TRUE, &this_iso,
1793                                                   u_axis, u_min, u_isostep,
1794                                                   num_iso_to_use,
1795                                                   v_axis, v_min, v_step,
1796                                                   num_sam_to_use,
1797                                                   NEED_PALETTE(this_plot));
1798                     }
1799                     /*}}} */
1800                 }               /* end of ITS A FUNCTION TO PLOT */
1801                 /* we saved it from first pass */
1802                 c_token = this_plot->token;
1803
1804                 /* one data file can make several plots */
1805                 do
1806                     this_plot = this_plot->next_sp;
1807                 while (this_plot && this_plot->token == c_token);
1808
1809             }                   /* !is_definition */
1810
1811             if (equals(c_token, ","))
1812                 c_token++;
1813             else
1814                 break;
1815
1816         }                       /* while(TRUE) */
1817
1818
1819         if (parametric) {
1820             /* Now actually fix the plot triplets to be single plots. */
1821             parametric_3dfixup(first_3dplot, &plot_num);
1822         }
1823     }                           /* some functions */
1824
1825     /* if first_3dplot is NULL, we have no functions or data at all.
1826        * This can happen, if you type "splot x=5", since x=5 is a
1827        * variable assignment
1828      */
1829     if (plot_num == 0 || first_3dplot == NULL) {
1830         int_error(c_token, "no functions or data to plot");
1831     }
1832
1833     axis_checked_extend_empty_range(FIRST_X_AXIS, "All points x value undefined");
1834     axis_checked_extend_empty_range(FIRST_Y_AXIS, "All points y value undefined");
1835     if (splot_map)
1836         axis_checked_extend_empty_range(FIRST_Z_AXIS, NULL); /* Suppress warning message */
1837     else
1838         axis_checked_extend_empty_range(FIRST_Z_AXIS, "All points z value undefined");
1839
1840     axis_revert_and_unlog_range(FIRST_X_AXIS);
1841     axis_revert_and_unlog_range(FIRST_Y_AXIS);
1842     axis_revert_and_unlog_range(FIRST_Z_AXIS);
1843
1844     setup_tics(FIRST_X_AXIS, 20);
1845     setup_tics(FIRST_Y_AXIS, 20);
1846     setup_tics(FIRST_Z_AXIS, 20);
1847
1848     set_plot_with_palette(plot_num, MODE_SPLOT);
1849     if (is_plot_with_palette()) {
1850         set_cbminmax();
1851         axis_checked_extend_empty_range(COLOR_AXIS, "All points of colorbox value undefined");
1852         setup_tics(COLOR_AXIS, 20);
1853         /* axis_revert_and_unlog_range(COLOR_AXIS); */
1854         /* fprintf(stderr,"plot3d.c: CB_AXIS.min=%g\tCB_AXIS.max=%g\n",CB_AXIS.min,CB_AXIS.max); */
1855     }
1856
1857     AXIS_WRITEBACK(FIRST_X_AXIS);
1858
1859     if (plot_num == 0 || first_3dplot == NULL) {
1860         int_error(c_token, "no functions or data to plot");
1861     }
1862     /* Creates contours if contours are to be plotted as well. */
1863
1864     if (draw_contour) {
1865         struct surface_points *this_plot;
1866         for (this_plot = first_3dplot, i = 0;
1867              i < plot_num;
1868              this_plot = this_plot->next_sp, i++) {
1869
1870             if (this_plot->contours) {
1871                 struct gnuplot_contours *cntrs = this_plot->contours;
1872
1873                 while (cntrs) {
1874                     struct gnuplot_contours *cntr = cntrs;
1875                     cntrs = cntrs->next;
1876                     free(cntr->coords);
1877                     free(cntr);
1878                 }
1879                 this_plot->contours = NULL;
1880             }
1881
1882             /* Allow individual surfaces to opt out of contouring */
1883             if (this_plot->opt_out_of_contours)
1884                 continue;
1885
1886             /* Make sure this one can be contoured. */
1887             if (!this_plot->has_grid_topology) {
1888                 int_warn(NO_CARET, "Notice: Cannot contour non grid data. Please use \"set dgrid3d\".");
1889             } else if (this_plot->plot_type == DATA3D) {
1890                 this_plot->contours = contour(this_plot->num_iso_read,
1891                                               this_plot->iso_crvs);
1892             } else {
1893                 this_plot->contours = contour(iso_samples_2,
1894                                               this_plot->iso_crvs);
1895             }
1896         }
1897     }                           /* draw_contour */
1898
1899     /* the following ~9 lines were moved from the end of the
1900      * function to here, as do_3dplot calles term->text, which
1901      * itself might process input events in mouse enhanced
1902      * terminals. For redrawing to work, line capturing and
1903      * setting the plot3d_num must already be done before
1904      * entering do_3dplot(). Thu Jan 27 23:54:49 2000 (joze) */
1905
1906     /* if we get here, all went well, so record the line for replot */
1907     if (plot_token != -1) {
1908         /* note that m_capture also frees the old replot_line */
1909         m_capture(&replot_line, plot_token, c_token - 1);
1910         plot_token = -1;
1911     }
1912
1913     /* record that all went well */
1914     plot3d_num=plot_num;
1915
1916     /* perform the plot */
1917     if (table_mode)
1918         print_3dtable(plot_num);
1919     else {
1920         START_LEAK_CHECK();     /* assert no memory leaks here ! */
1921         do_3dplot(first_3dplot, plot_num, 0);
1922         END_LEAK_CHECK();
1923
1924         /* after do_3dplot(), axis_array[] and max_array[].min
1925          * contain the plotting range actually used (rounded
1926          * to tic marks, not only the min/max data values)
1927          * --> save them now for writeback if requested
1928          */
1929         SAVE_WRITEBACK_ALL_AXES;
1930         /* update GPVAL_ variables available to user */
1931         update_gpval_variables(1);
1932     }
1933 }
1934
1935
1936
1937 /*
1938  * The hardest part of this routine is collapsing the FUNC plot types in the
1939  * list (which are gauranteed to occur in (x,y,z) triplets while preserving
1940  * the non-FUNC type plots intact.  This means we have to work our way
1941  * through various lists.  Examples (hand checked):
1942  * start_plot:F1->F2->F3->NULL ==> F3->NULL
1943  * start_plot:F1->F2->F3->F4->F5->F6->NULL ==> F3->F6->NULL
1944  * start_plot:F1->F2->F3->D1->D2->F4->F5->F6->D3->NULL ==>
1945  * F3->D1->D2->F6->D3->NULL
1946  *
1947  * x and y ranges now fixed in eval_3dplots
1948  */
1949 static void
1950 parametric_3dfixup(struct surface_points *start_plot, int *plot_num)
1951 {
1952     struct surface_points *xp, *new_list, *free_list = NULL;
1953     struct surface_points **last_pointer = &new_list;
1954     size_t tlen;
1955     int i, surface;
1956     char *new_title;
1957
1958     /*
1959      * Ok, go through all the plots and move FUNC3D types together.  Note:
1960      * this originally was written to look for a NULL next pointer, but
1961      * gnuplot wants to be sticky in grabbing memory and the right number of
1962      * items in the plot list is controlled by the plot_num variable.
1963      *
1964      * Since gnuplot wants to do this sticky business, a free_list of
1965      * surface_points is kept and then tagged onto the end of the plot list
1966      * as this seems more in the spirit of the original memory behavior than
1967      * simply freeing the memory.  I'm personally not convinced this sort of
1968      * concern is worth it since the time spent computing points seems to
1969      * dominate any garbage collecting that might be saved here...
1970      */
1971     new_list = xp = start_plot;
1972     for (surface = 0; surface < *plot_num; surface++) {
1973         if (xp->plot_type == FUNC3D) {
1974             struct surface_points *yp = xp->next_sp;
1975             struct surface_points *zp = yp->next_sp;
1976
1977             /* Here's a FUNC3D parametric function defined as three parts.
1978              * Go through all the points and assign the x's and y's from xp and
1979              * yp to zp. min/max already done
1980              */
1981             struct iso_curve *xicrvs = xp->iso_crvs;
1982             struct iso_curve *yicrvs = yp->iso_crvs;
1983             struct iso_curve *zicrvs = zp->iso_crvs;
1984
1985             (*plot_num) -= 2;
1986
1987             assert(INRANGE < OUTRANGE && OUTRANGE < UNDEFINED);
1988
1989             while (zicrvs) {
1990                 struct coordinate GPHUGE *xpoints = xicrvs->points;
1991                 struct coordinate GPHUGE *ypoints = yicrvs->points;
1992                 struct coordinate GPHUGE *zpoints = zicrvs->points;
1993
1994                 for (i = 0; i < zicrvs->p_count; ++i) {
1995                     zpoints[i].x = xpoints[i].z;
1996                     zpoints[i].y = ypoints[i].z;
1997                     if (zpoints[i].type < xpoints[i].type)
1998                         zpoints[i].type = xpoints[i].type;
1999                     if (zpoints[i].type < ypoints[i].type)
2000                         zpoints[i].type = ypoints[i].type;
2001
2002                 }
2003                 xicrvs = xicrvs->next;
2004                 yicrvs = yicrvs->next;
2005                 zicrvs = zicrvs->next;
2006             }
2007
2008 #if 0 /* FIXME HBB 20001101: seems to cause a crash */
2009             if (first_3dplot)
2010                 sp_free(first_3dplot);
2011             first_3dplot = NULL;
2012 #endif
2013
2014             /* Ok, fix up the title to include xp and yp plots. */
2015             if (((xp->title && xp->title[0] != '\0') ||
2016                  (yp->title && yp->title[0] != '\0')) && zp->title) {
2017                 tlen = (xp->title ? strlen(xp->title) : 0) +
2018                     (yp->title ? strlen(yp->title) : 0) +
2019                     (zp->title ? strlen(zp->title) : 0) + 5;
2020                 new_title = gp_alloc(tlen, "string");
2021                 new_title[0] = 0;
2022                 if (xp->title && xp->title[0] != '\0') {
2023                     strcat(new_title, xp->title);
2024                     strcat(new_title, ", ");    /* + 2 */
2025                 }
2026                 if (yp->title && yp->title[0] != '\0') {
2027                     strcat(new_title, yp->title);
2028                     strcat(new_title, ", ");    /* + 2 */
2029                 }
2030                 strcat(new_title, zp->title);
2031                 free(zp->title);
2032                 zp->title = new_title;
2033             }
2034             /* add xp and yp to head of free list */
2035             assert(xp->next_sp == yp);
2036             yp->next_sp = free_list;
2037             free_list = xp;
2038
2039             /* add zp to tail of new_list */
2040             *last_pointer = zp;
2041             last_pointer = &(zp->next_sp);
2042             xp = zp->next_sp;
2043         } else {                /* its a data plot */
2044             assert(*last_pointer == xp);        /* think this is true ! */
2045             last_pointer = &(xp->next_sp);
2046             xp = xp->next_sp;
2047         }
2048     }
2049
2050     /* Ok, append free list and write first_plot */
2051     *last_pointer = free_list;
2052     first_3dplot = new_list;
2053 }
2054