home *** CD-ROM | disk | FTP | other *** search
/ System Booster / System Booster.iso / Archives / GNU / GNUPLOTsrc.lha / plot3d.c < prev    next >
Encoding:
C/C++ Source or Header  |  1996-01-22  |  44.9 KB  |  1,507 lines
  1. #ifndef lint
  2. static char    *RCSid = "$Id: plot3d.c,v 1.15 1995/12/02 21:16:43 drd Exp $";
  3. #endif
  4.  
  5.  
  6. /* GNUPLOT - command.c */
  7. /*
  8.  * Copyright (C) 1986 - 1993   Thomas Williams, Colin Kelley
  9.  * 
  10.  * Permission to use, copy, and distribute this software and its documentation
  11.  * for any purpose with or without fee is hereby granted, provided that the
  12.  * above copyright notice appear in all copies and that both that copyright
  13.  * notice and this permission notice appear in supporting documentation.
  14.  * 
  15.  * Permission to modify the software is granted, but not the right to distribute
  16.  * the modified code.  Modifications are to be distributed as patches to
  17.  * released version.
  18.  * 
  19.  * This software is provided "as is" without express or implied warranty.
  20.  * 
  21.  * 
  22.  * AUTHORS
  23.  * 
  24.  * Original Software: Thomas Williams,  Colin Kelley.
  25.  * 
  26.  * Gnuplot 2.0 additions: Russell Lang, Dave Kotz, John Campbell.
  27.  * 
  28.  * Gnuplot 3.0 additions: Gershon Elber and many others.
  29.  * 
  30.  * formerly part of command.c
  31.  * 
  32.  * There is a mailing list for gnuplot users. Note, however, that the
  33.  * newsgroup 
  34.  *    comp.graphics.gnuplot 
  35.  * is identical to the mailing list (they
  36.  * both carry the same set of messages). We prefer that you read the
  37.  * messages through that newsgroup, to subscribing to the mailing list.
  38.  * (If you can read that newsgroup, and are already on the mailing list,
  39.  * please send a message info-gnuplot-request@dartmouth.edu, asking to be
  40.  * removed from the mailing list.)
  41.  *
  42.  * The address for mailing to list members is
  43.  *       info-gnuplot@dartmouth.edu
  44.  * and for mailing administrative requests is 
  45.  *       info-gnuplot-request@dartmouth.edu
  46.  * The mailing list for bug reports is 
  47.  *       bug-gnuplot@dartmouth.edu
  48.  * The list of those interested in beta-test versions is
  49.  *       info-gnuplot-beta@dartmouth.edu
  50.  */
  51.  
  52. #include <math.h>
  53. #include <ctype.h>
  54. #include <assert.h>
  55.  
  56. #include "plot.h"
  57. #include "setshow.h"
  58. #include "binary.h"
  59. #ifndef _Windows
  60. #include "help.h"
  61. #else
  62. #define MAXSTR 255
  63. #endif
  64.  
  65. #if defined(ATARI) || defined(MTOS)
  66. #ifdef __PUREC__
  67. #include <ext.h>
  68. #include <tos.h>
  69. #include <aes.h>
  70. #include <float.h>        /* get FLT_MAX */
  71. #else
  72. #include <osbind.h>
  73. #include <aesbind.h>
  74. #endif /* __PUREC__ */
  75. #endif /* ATARI || MTOS */
  76.  
  77. #ifndef STDOUT
  78. #define STDOUT 1
  79. #endif
  80.  
  81.  
  82. #define inrange(z,min,max) ((min<max) ? ((z>=min)&&(z<=max)) : ((z>=max)&&(z<=min)) )
  83.  
  84.  
  85. /* static prototypes */
  86.  
  87. static void get_3ddata __P((struct surface_points *this_plot));
  88. static void print_3dtable __P((int pcount));
  89. static void eval_3dplots __P((void));
  90. static void grid_nongrid_data __P((struct surface_points *this_plot));
  91. static void parametric_3dfixup __P((struct surface_points *start_plot, int *plot_num));
  92.  
  93. /* the curves/surfaces of the plot */
  94. struct surface_points *first_3dplot = NULL;
  95. static struct udft_entry plot_func;
  96.  
  97.  
  98. extern struct udft_entry *dummy_func;
  99. extern int datatype[];
  100. extern char timefmt[];
  101. extern TBOOLEAN         is_3d_plot;
  102. extern int plot_token;
  103.  
  104. /* in order to support multiple axes, and to
  105.  * simplify ranging in parametric plots, we use
  106.  * arrays to store some things. For 2d plots,
  107.  * elements are  y1=0 x1=1 y2=2 x2=3
  108.  * for 3d,  z=0, x=1, y=2
  109.  * these are given symbolic names in plot.h
  110.  */
  111.  
  112. extern double          min_array[AXIS_ARRAY_SIZE], max_array[AXIS_ARRAY_SIZE];
  113. extern int             auto_array[AXIS_ARRAY_SIZE];
  114. extern TBOOLEAN        log_array[AXIS_ARRAY_SIZE];
  115. extern double          base_array[AXIS_ARRAY_SIZE];
  116. extern double          log_base_array[AXIS_ARRAY_SIZE];
  117.  
  118. /* some file-wide variables to store which axis we are using */
  119. static int x_axis, y_axis, z_axis;
  120.  
  121.  
  122. /* info from datafile module */
  123. extern int df_datum;
  124. extern int df_line_number;
  125. extern int df_no_use_specs;
  126. extern int df_eof;
  127. extern int df_timecol[];
  128. extern TBOOLEAN df_binary;
  129.  
  130. #ifndef min
  131. #define min(a,b)    ((a)>(b) ? (b) : (a))
  132. #endif
  133. #define Inc_c_token if (++c_token >= num_tokens)    \
  134.                         int_error ("Syntax error", c_token);
  135.  
  136. extern int reverse_range[];
  137.  
  138. /*
  139.  * IMHO, code is getting too cluttered with repeated chunks of
  140.  * code. Some macros to simplify, I hope.
  141.  *
  142.  * do { } while(0) is comp.lang.c recommendation for complex macros
  143.  * also means that break can be specified as an action, and it will
  144.  * 
  145.  */
  146.  
  147. /*  copy scalar data to arrays
  148.  * optimiser should optimise infinite away
  149.  * dont know we have to support ranges [10:-10] - lets reverse
  150.  * it for now, then fix it at the end.
  151.  */
  152. #define INIT_ARRAYS(axis, min, max, auto, is_log, base, log_base, infinite) \
  153. do{if ((auto_array[axis]=auto)==0 && max<min) {\
  154.     min_array[axis]=max;\
  155.    max_array[axis]=min; /* we will fix later */ \
  156.  } else { \
  157.     min_array[axis]=(infinite && (auto&1)) ? VERYLARGE : min; \
  158.     max_array[axis]=(infinite && (auto&2)) ? -VERYLARGE : max; \
  159.  } \
  160.  log_array[axis]=is_log; base_array[axis]=base; log_base_array[axis]=log_base;\
  161. }while(0)
  162.  
  163. /* handle reversed ranges */
  164. #define CHECK_REVERSE(axis) \
  165. do{\
  166.  if (auto_array[axis]==0 && max_array[axis] < min_array[axis]) {\
  167.   double temp=min_array[axis]; min_array[axis]=max_array[axis]; max_array[axis]=temp;\
  168.   reverse_range[axis]=1; \
  169.  } else reverse_range[axis] = (range_flags[axis]&RANGE_REVERSE); \
  170. }while(0)
  171.  
  172.  
  173. /* get optional [min:max] */
  174. #define LOAD_RANGE(axis) \
  175. do {\
  176.  if (equals(c_token, "[")) { \
  177.   c_token++; \
  178.   auto_array[axis] = load_range(axis,&min_array[axis], &max_array[axis], auto_array[axis]);\
  179.   if (!equals(c_token, "]"))\
  180.    int_error("']' expected", c_token);\
  181.   c_token++;\
  182.  }\
  183. } while (0)
  184.  
  185.  
  186. /* store VALUE or log(VALUE) in STORE, set TYPE as appropriate
  187.  * Do OUT_ACTION or UNDEF_ACTION as appropriate
  188.  * adjust range provided type is INRANGE (ie dont adjust y if x is outrange
  189.  * VALUE must not be same as STORE
  190.  */
  191.  
  192. #define STORE_WITH_LOG_AND_FIXUP_RANGE(STORE, VALUE, TYPE, AXIS, OUT_ACTION, UNDEF_ACTION)\
  193. do { if (log_array[AXIS]) { if (VALUE<0.0) {TYPE=UNDEFINED; UNDEF_ACTION; break;} \
  194.               else if (VALUE==0.0){STORE=-VERYLARGE; TYPE=OUTRANGE; OUT_ACTION; break;} \
  195.               else { STORE=log(VALUE)/log_base_array[AXIS]; } \
  196.      } else STORE=VALUE; \
  197.      if (TYPE != INRANGE) break;  /* dont set y range if x is outrange, for example */ \
  198.      if ( VALUE<min_array[AXIS] ) \
  199.       if (auto_array[AXIS] & 1) min_array[AXIS]=VALUE; else { TYPE=OUTRANGE; OUT_ACTION; break; }  \
  200.      if ( VALUE>max_array[AXIS] ) \
  201.      if (auto_array[AXIS] & 2) max_array[AXIS]=VALUE; else { TYPE=OUTRANGE; OUT_ACTION; }   \
  202. } while(0)
  203.      
  204. /* use this instead empty macro arguments to work around NeXT cpp bug */
  205. /* if this fails on any system, we might use ((void)0) */
  206. #define NOOP /* */
  207.  
  208. /* check axis range is not too small -
  209.  * extend if you can (autoscale), else report error
  210.  */
  211. #ifdef ANSI_C
  212. # define STRINGIFY(x) #x
  213. # define RANGE_MSG(x) #x " range is less than `zero`"
  214. #else
  215. # define STRINGIFY(x) "x"
  216. # define RANGE_MSG(x) "x range is less than `zero`"
  217. #endif
  218.  
  219. #define FIXUP_RANGE(AXIS, WHICH) \
  220. do{if (fabs(max_array[AXIS] - min_array[AXIS]) < zero)    \
  221.     if (auto_array[AXIS]) { /* widen range */  \
  222.      fprintf(stderr, "Warning: empty %s range [%g:%g], ", STRINGIFY(WHICH), min_array[AXIS], max_array[AXIS]);      \
  223.      if (fabs(min_array[AXIS]) < zero) { \
  224.       if (auto_array[AXIS] & 1) min_array[AXIS] = -1.0; \
  225.       if (auto_array[AXIS] & 2) max_array[AXIS] = 1.0;   \
  226.      } else if (max_array[AXIS] < 0) { \
  227.       if (auto_array[AXIS] & 1) min_array[AXIS] *= 1.1; if (auto_array[AXIS] & 2) max_array[AXIS] *= 0.9;    \
  228.      } else { if (auto_array[AXIS] & 1) min_array[AXIS] *= 0.9; if (auto_array[AXIS] & 2) max_array[AXIS] *= 1.1;  }  \
  229.      fprintf(stderr, "adjusting to [%g:%g]\n", min_array[AXIS], max_array[AXIS]);          \
  230.     } else int_error(RANGE_MSG(WHICH), c_token);   \
  231. }while(0)
  232.  
  233. /* check range and take logs of min and max if logscale
  234.  * this also restores min and max for ranges like [10:-10]
  235.  */
  236.  
  237. #define FIXUP_RANGE_FOR_LOG(AXIS, WHICH) \
  238. do { if (reverse_range[AXIS]) { \
  239.       double temp = min_array[AXIS]; \
  240.       min_array[AXIS]=max_array[AXIS]; \
  241.       max_array[AXIS]=temp; \
  242.      }\
  243.      if (log_array[AXIS]) { \
  244.       if (min_array[AXIS]<=0.0 || max_array[AXIS]<=0.0) \
  245.        int_error(RANGE_MSG(WHICH), NO_CARET); \
  246.       min_array[AXIS] = log(min_array[AXIS])/log_base_array[AXIS]; \
  247.       max_array[AXIS] = log(max_array[AXIS])/log_base_array[AXIS];  \
  248. } } while(0)
  249.  
  250.  
  251.  
  252. /* support for dynamic size of input line */
  253.  
  254. void plot3drequest()
  255. /*
  256.  * in the parametric case we would say splot [u= -Pi:Pi] [v= 0:2*Pi] [-1:1]
  257.  * [-1:1] [-1:1] sin(v)*cos(u),sin(v)*cos(u),sin(u) in the non-parametric
  258.  * case we would say only splot [x= -2:2] [y= -5:5] sin(x)*cos(y)
  259.  * 
  260.  */
  261. {
  262.     TBOOLEAN         changed;
  263.     int             dummy_token0 = -1, dummy_token1 = -1;
  264.  
  265.     is_3d_plot = TRUE;
  266.  
  267.     if (parametric && strcmp(dummy_var[0], "t") == 0) {
  268.     strcpy(dummy_var[0], "u");
  269.     strcpy(dummy_var[1], "v");
  270.     }
  271.     autoscale_lx = autoscale_x;
  272.     autoscale_ly = autoscale_y;
  273.     autoscale_lz = autoscale_z;
  274.  
  275.     if (!term)            /* unknown */
  276.     int_error("use 'set term' to set terminal type first", c_token);
  277.  
  278.     if (equals(c_token, "[")) {
  279.     c_token++;
  280.     if (isletter(c_token)) {
  281.         if (equals(c_token + 1, "=")) {
  282.         dummy_token0 = c_token;
  283.         c_token += 2;
  284.         } else {
  285.         /* oops; probably an expression with a variable. */
  286.         /* Parse it as an xmin expression. */
  287.         /* used to be: int_error("'=' expected",c_token); */
  288.         }
  289.     }
  290.     changed = parametric ? load_range(U_AXIS,&umin, &umax, autoscale_lu) : load_range(FIRST_X_AXIS,&xmin, &xmax, autoscale_lx);
  291.     if (!equals(c_token, "]"))
  292.         int_error("']' expected", c_token);
  293.     c_token++;
  294.     /* if (changed) */
  295.         if(parametric) 
  296.             /* autoscale_lu = FALSE; */
  297.             autoscale_lu = changed;
  298.         else
  299.             /* autoscale_lx = FALSE; */
  300.             autoscale_lx = changed;
  301.     }
  302.     if (equals(c_token, "[")) {
  303.     c_token++;
  304.     if (isletter(c_token)) {
  305.         if (equals(c_token + 1, "=")) {
  306.         dummy_token1 = c_token;
  307.         c_token += 2;
  308.         } else {
  309.         /* oops; probably an expression with a variable. */
  310.         /* Parse it as an xmin expression. */
  311.         /* used to be: int_error("'=' expected",c_token); */
  312.         }
  313.     }
  314.     changed = parametric ? load_range(V_AXIS,&vmin, &vmax, autoscale_lv) : load_range(FIRST_Y_AXIS,&ymin, &ymax, autoscale_ly);
  315.     if (!equals(c_token, "]"))
  316.         int_error("']' expected", c_token);
  317.     c_token++;
  318.     /* if (changed) */
  319.         if(parametric) 
  320.             /* autoscale_lv = FALSE; */
  321.             autoscale_lv = changed;
  322.         else
  323.             /* autoscale_ly = FALSE; */
  324.             autoscale_ly = changed;
  325.     }
  326.  
  327.     if (parametric) {
  328.     if (equals(c_token, "[")) {    /* set optional x (parametric) or z ranges */
  329.     c_token++;
  330.     autoscale_lx = load_range(FIRST_X_AXIS,&xmin, &xmax, autoscale_lx);
  331.     if (!equals(c_token, "]"))
  332.         int_error("']' expected", c_token);
  333.     c_token++;
  334.     }
  335.     if (equals(c_token, "[")) {    /* set optional y ranges */
  336.     c_token++;
  337.     autoscale_ly = load_range(FIRST_Y_AXIS,&ymin, &ymax, autoscale_ly);
  338.     if (!equals(c_token, "]"))
  339.         int_error("']' expected", c_token);
  340.     c_token++;
  341.     }
  342.      } /* parametric */
  343.      
  344.     if (equals(c_token, "[")) {    /* set optional z ranges */
  345.     c_token++;
  346.     autoscale_lz = load_range(FIRST_Z_AXIS,&zmin, &zmax, autoscale_lz);
  347.     if (!equals(c_token, "]"))
  348.         int_error("']' expected", c_token);
  349.     c_token++;
  350.     }
  351.  
  352.     CHECK_REVERSE(FIRST_X_AXIS);
  353.     CHECK_REVERSE(FIRST_Y_AXIS);
  354.     CHECK_REVERSE(FIRST_Z_AXIS);
  355.  
  356.     /* use the default dummy variable unless changed */
  357.     if (dummy_token0 >= 0)
  358.     copy_str(c_dummy_var[0], dummy_token0, MAX_ID_LEN);
  359.     else
  360.     (void) strcpy(c_dummy_var[0], dummy_var[0]);
  361.  
  362.     if (dummy_token1 >= 0)
  363.     copy_str(c_dummy_var[1], dummy_token1, MAX_ID_LEN);
  364.     else
  365.     (void) strcpy(c_dummy_var[1], dummy_var[1]);
  366.  
  367.     eval_3dplots();
  368. }
  369.  
  370.  
  371. static void grid_nongrid_data(this_plot)
  372. struct surface_points *this_plot;
  373. {
  374.     int i, j, k;
  375.     double x, y, z, w, dx, dy, xmin, xmax, ymin, ymax;
  376.     struct iso_curve *old_iso_crvs = this_plot->iso_crvs;
  377.     struct iso_curve *icrv, *oicrv, *oicrvs;
  378.  
  379.     /* Compute XY bounding box on the original data. */
  380.     xmin = xmax = old_iso_crvs->points[0].x;
  381.     ymin = ymax = old_iso_crvs->points[0].y;
  382.     for (icrv = old_iso_crvs; icrv != NULL; icrv = icrv->next) {
  383.     struct coordinate GPHUGE *points = icrv->points;
  384.  
  385.     for (i = 0; i < icrv->p_count; i++, points++) {
  386.         if (xmin > points->x)
  387.         xmin = points->x;
  388.         if (xmax < points->x)
  389.         xmax = points->x;
  390.         if (ymin > points->y)
  391.         ymin = points->y;
  392.         if (ymax < points->y)
  393.         ymax = points->y;
  394.     }
  395.     }
  396.  
  397.     dx = (xmax - xmin) / (dgrid3d_col_fineness - 1);
  398.     dy = (ymax - ymin) / (dgrid3d_row_fineness - 1);
  399.  
  400.     /* Create the new grid structure, and compute the low pass filtering from
  401.      * non grid to grid structure.
  402.      */
  403.     this_plot->iso_crvs = NULL;
  404.     this_plot->num_iso_read = dgrid3d_col_fineness;
  405.     this_plot->has_grid_topology = TRUE;
  406.     for (i = 0, x = xmin; i < dgrid3d_col_fineness; i++, x += dx) {
  407.     struct coordinate GPHUGE *points;
  408.  
  409.     icrv = iso_alloc(dgrid3d_row_fineness + 1);
  410.     icrv->p_count = dgrid3d_row_fineness;
  411.     icrv->next = this_plot->iso_crvs;
  412.     this_plot->iso_crvs = icrv;
  413.     points = icrv->points;
  414.  
  415.     for (j = 0, y = ymin; j < dgrid3d_row_fineness; j++, y += dy, points++) {
  416.         z = w = 0.0;
  417.  
  418.         for (oicrv = old_iso_crvs; oicrv != NULL; oicrv = oicrv->next) {
  419.         struct coordinate GPHUGE *opoints = oicrv->points;
  420.         for (k = 0; k < oicrv->p_count; k++, opoints++) {
  421.             double dist,
  422.                dist_x = fabs( opoints->x - x ),
  423.                dist_y = fabs( opoints->y - y );
  424.  
  425.             switch (dgrid3d_norm_value) {
  426.             case 1:
  427.                 dist = dist_x + dist_y;
  428.                 break;
  429.             case 2:
  430.                 dist = dist_x * dist_x + dist_y * dist_y;
  431.                 break;
  432.             case 4:
  433.                 dist = dist_x * dist_x + dist_y * dist_y;
  434.                 dist *= dist;
  435.                 break;
  436.             case 8:
  437.                 dist = dist_x * dist_x + dist_y * dist_y;
  438.                 dist *= dist;
  439.                 dist *= dist;
  440.                 break;
  441.             case 16:
  442.                 dist = dist_x * dist_x + dist_y * dist_y;
  443.                 dist *= dist;
  444.                 dist *= dist;
  445.                 dist *= dist;
  446.                 break;
  447.             default:
  448.                             dist = pow( dist_x, (double)dgrid3d_norm_value ) +
  449.                                    pow( dist_y, (double)dgrid3d_norm_value );
  450.                 break;
  451.             }
  452.  
  453.             /* The weight of this point is inverse proportional
  454.              * to the distance.
  455.              */
  456.             if ( dist == 0.0 )
  457. #if !defined(AMIGA_SC_6_1) && !defined(__PUREC__)
  458.             dist = VERYLARGE;
  459. #else /* !AMIGA_SC_6_1 && !__PUREC__ */
  460.             /* Multiplying VERYLARGE by opoints->z below
  461.              * might yield Inf (i.e. a number that can't
  462.              * be represented on the machine). This will
  463.              * result in points->z being set to NaN. It's
  464.              * better to have a pretty large number that is
  465.              * also on the safe side... The numbers that are
  466.              * read by gnuplot are float values anyway, so
  467.              * they can't be bigger than FLT_MAX. So setting
  468.              * dist to FLT_MAX^2 will make dist pretty large
  469.              * with respect to any value that has been read. */
  470.             dist = ((double)FLT_MAX)*((double)FLT_MAX);
  471. #endif /* !AMIGA_SC_6_1 && !__PUREC__ */
  472.             else
  473.             dist = 1.0 / dist;
  474.  
  475.             z += opoints->z * dist;
  476.             w += dist;
  477.         }
  478.         }
  479.  
  480.         points->x = x;
  481.         points->y = y;
  482.         points->z = z / w;
  483.         points->type = INRANGE;
  484.     }
  485.     }
  486.     
  487.     /* Delete the old non grid data. */
  488.     for (oicrvs = old_iso_crvs; oicrvs != NULL;) {
  489.     oicrv = oicrvs;
  490.     oicrvs = oicrvs->next;
  491.     iso_free(oicrv);
  492.     }
  493. }
  494.  
  495. static void get_3ddata(this_plot)
  496.     struct surface_points *this_plot;
  497. /* this_plot->token is end of datafile spec, before title etc
  498.  * will be moved passed title etc after we return
  499.  */
  500. {
  501.     int xdatum=0;
  502.     int ydatum=0;
  503.     int i,j;
  504.     double v[3];
  505.     int pt_in_iso_crv=0;
  506.     struct iso_curve *this_iso;
  507.  
  508.     if (mapping3d == MAP3D_CARTESIAN)
  509.     {     if (df_no_use_specs == 2)
  510.             int_error("Need 1 or 3 columns for cartesian data", this_plot->token);
  511.     }
  512.     else
  513.     {
  514.         if (df_no_use_specs == 1)
  515.             int_error("Need 2 or 3 columns for polar data", this_plot->token);
  516.     }
  517.  
  518.     this_plot->num_iso_read = 0;
  519.     this_plot->has_grid_topology = TRUE;
  520.  
  521.     /* we ought to keep old memory - most likely case
  522.      * is a replot, so it will probably exactly fit into
  523.      * memory already allocated ?
  524.      */
  525.     if (this_plot->iso_crvs != NULL) {
  526.         struct iso_curve *icrv, *icrvs = this_plot->iso_crvs;
  527.  
  528.         while (icrvs) {
  529.             icrv = icrvs;
  530.             icrvs = icrvs->next;
  531.             iso_free(icrv);
  532.         }
  533.         this_plot->iso_crvs = NULL;
  534.     }
  535.  
  536.     /* data file is already open */
  537.  
  538.     if (df_binary)
  539.         xdatum=df_3dbinary(this_plot);
  540.     else    {
  541.         /*{{{  read surface from text file*/
  542.         struct iso_curve *this_iso=iso_alloc(samples);
  543.         struct coordinate GPHUGE *cp;
  544.         double x,y,z;
  545.         
  546.         while ((j=df_readline(v,3)) != DF_EOF)  {
  547.             if (j==DF_SECOND_BLANK)
  548.                 break;  /* two blank lines */
  549.             if (j==DF_FIRST_BLANK) {
  550.                 /* one blank line */
  551.                 if (pt_in_iso_crv == 0) {
  552.                     if (xdatum == 0)
  553.                     continue;
  554.                     pt_in_iso_crv = xdatum;
  555.                 }
  556.                 if (xdatum > 0) {
  557.                     this_iso->p_count = xdatum;
  558.                     this_iso->next = this_plot->iso_crvs;
  559.                     this_plot->iso_crvs = this_iso;
  560.                     this_plot->num_iso_read++;
  561.         
  562.                     if (xdatum != pt_in_iso_crv)
  563.                     this_plot->has_grid_topology = FALSE;
  564.         
  565.                     this_iso = iso_alloc(pt_in_iso_crv);
  566.                     xdatum = 0;
  567.                     ydatum++;
  568.                 }
  569.                 continue;
  570.             }
  571.         
  572.             /* its a data point or undefined */
  573.         
  574.             if (xdatum >= this_iso->p_max) {
  575.             /*
  576.              * overflow about to occur. Extend size of points[] array. We
  577.              * either double the size, or add 1000 points, whichever is a
  578.              * smaller increment. Note i=p_max.
  579.              */
  580.             iso_extend(this_iso,
  581.                    xdatum + (xdatum < 1000 ? xdatum : 1000));
  582.             }
  583.   
  584.         cp = this_iso->points + xdatum;
  585.   
  586.         if (j==DF_UNDEFINED) {
  587.             cp->type=UNDEFINED;
  588.             continue;
  589.         }
  590.      
  591.         cp->type=INRANGE; /* unless we find out different */
  592.   
  593.         switch(mapping3d) {
  594.             case MAP3D_CARTESIAN:
  595.                 switch(j) {
  596.                     case 1:
  597.                             x = xdatum;
  598.                             y = ydatum;
  599.                             z = v[0];
  600.                         break;
  601.                     case 3:
  602.                             x = v[0];
  603.                             y = v[1];
  604.                             z = v[2];
  605.                         break;
  606.                     default:
  607.                         {
  608.                             char msg[80];
  609.                             sprintf(msg, "Need 1 or 3 columns - line %d", df_line_number);
  610.                             int_error(msg, this_plot->token);
  611.                             return; /* avoid gcc -Wuninitialised for x,y,z */
  612.                         }
  613.                 }
  614.                 break;
  615.             case MAP3D_SPHERICAL:
  616.                 if (j<2)
  617.                     int_error("Need 2 or 3 columns", this_plot->token);
  618.                 if (j<3)
  619.                     v[2]=1; /* default radius */
  620.                 if (angles_format == ANGLES_DEGREES) {
  621.                     v[0] *= DEG2RAD;    /* Convert to radians. */
  622.                     v[1] *= DEG2RAD;
  623.                 }
  624.                 x = v[2]*cos(v[0]) * cos(v[1]);
  625.                 y = v[2]*sin(v[0]) * cos(v[1]);
  626.                 z = v[2]*sin(v[1]);
  627.                 break;
  628.             case MAP3D_CYLINDRICAL:
  629.                 if (j<2)
  630.                     int_error("Need 2 or 3 columns", this_plot->token);
  631.                 if (j<3)
  632.                     v[2]=1; /* default radius */
  633.                 if (angles_format == ANGLES_DEGREES) {
  634.                     v[0] *= DEG2RAD;    /* Convert to radians. */
  635.                 }
  636.                 x = v[2]*cos(v[0]);
  637.                 y = v[2]*sin(v[0]);
  638.                 z = v[1];
  639.                 break;
  640.              default:
  641.                 int_error("Internal error : Unknown mapping type", NO_CARET);
  642.                 return;
  643.         }
  644.   
  645.         /* adjust for logscales. Set min/max and point types. store in cp */
  646.         cp->type=INRANGE;
  647.         /* cannot use continue, as macro is wrapped in a loop. I regard this as correct goto use */
  648.         STORE_WITH_LOG_AND_FIXUP_RANGE(cp->x, x, cp->type, x_axis, NOOP, goto come_here_if_undefined );
  649.         STORE_WITH_LOG_AND_FIXUP_RANGE(cp->y, y, cp->type, y_axis, NOOP, goto come_here_if_undefined );
  650.         STORE_WITH_LOG_AND_FIXUP_RANGE(cp->z, z, cp->type, z_axis, NOOP, goto come_here_if_undefined );
  651.   
  652.             /* some may complain, but I regard this as the correct use of goto */
  653.             come_here_if_undefined:
  654.             ++xdatum;
  655.         }  /* end of whileloop - end of surface */
  656.         
  657.         if (xdatum > 0) {
  658.             this_plot->num_iso_read++;    /* Update last iso. */
  659.             this_iso->p_count = xdatum;
  660.         
  661.             this_iso->next = this_plot->iso_crvs;
  662.             this_plot->iso_crvs = this_iso;
  663.         
  664.             if (xdatum != pt_in_iso_crv)
  665.             this_plot->has_grid_topology = FALSE;
  666.         
  667.         } else {
  668.             iso_free(this_iso);    /* Free last allocation. */
  669.         }
  670.         
  671.         if (dgrid3d) grid_nongrid_data(this_plot);
  672.         /*}}}*/
  673.     }
  674.  
  675.     if (this_plot->num_iso_read <= 1)
  676.     this_plot->has_grid_topology = FALSE;
  677.     if (this_plot->has_grid_topology && !hidden3d) {
  678.     struct iso_curve *new_icrvs = NULL;
  679.     int             num_new_iso = this_plot->iso_crvs->p_count, len_new_iso = this_plot->num_iso_read;
  680.  
  681.     /* Now we need to set the other direction (pseudo) isolines. */
  682.     for (i = 0; i < num_new_iso; i++) {
  683.         struct iso_curve *new_icrv = iso_alloc(len_new_iso);
  684.  
  685.         new_icrv->p_count = len_new_iso;
  686.  
  687.         for (j = 0, this_iso = this_plot->iso_crvs;
  688.          this_iso != NULL;
  689.          j++, this_iso = this_iso->next) {
  690.         /* copy whole point struct to get type too.
  691.          * wasteful for windows, with padding */
  692.         /* more efficient would be extra pointer to same struct */
  693.         new_icrv->points[j]=this_iso->points[i];
  694.         }
  695.  
  696.         new_icrv->next = new_icrvs;
  697.         new_icrvs = new_icrv;
  698.     }
  699.  
  700.     /* Append the new iso curves after the read ones. */
  701.     for (this_iso = this_plot->iso_crvs;
  702.          this_iso->next != NULL;
  703.          this_iso = this_iso->next);
  704.         this_iso->next = new_icrvs;
  705.     }
  706.  
  707. }
  708.  
  709.  
  710.  
  711. static void print_3dtable(pcount)
  712. int pcount;
  713. {
  714.     register struct surface_points *this_plot;
  715.     int             i, curve,surface;
  716.     struct iso_curve *icrvs;
  717.     struct coordinate GPHUGE *points;
  718.  
  719.     for (surface=0, this_plot=first_3dplot ; surface < pcount; 
  720.         this_plot=this_plot->next_sp, surface++){
  721.         fprintf(outfile, "\n#Surface %d of %d surfaces\n", surface, pcount);
  722.         icrvs = this_plot->iso_crvs;
  723.         curve = 0;
  724.  
  725.         if (draw_surface) {
  726.             /* only the curves in one direction */
  727.             while(icrvs && curve < this_plot->num_iso_read){
  728.             fprintf(outfile, "\n#IsoCurve %d, %d points\n#x y z type\n", curve, icrvs->p_count);
  729.             for(i=0, points = icrvs->points; i < icrvs->p_count; i++){
  730.                 fprintf(outfile, "%g %g %g %c\n",
  731.                 points[i].x,
  732.                 points[i].y,
  733.                 points[i].z,
  734.                 points[i].type == INRANGE ? 'i'
  735.                 : points[i].type == OUTRANGE ? 'o'
  736.                 : 'u');
  737.             }
  738.             icrvs = icrvs->next;
  739.             curve++;
  740.             }
  741.             putc('\n', outfile);
  742.         }
  743.  
  744.         if (draw_contour) {
  745.             int number=0;
  746.             struct gnuplot_contours *c=this_plot->contours;
  747.             while (c) {
  748.                 int count=c->num_pts;
  749.                 struct coordinate GPHUGE *p=c->coords;
  750.                 if (c->isNewLevel)
  751.                     /* dont display count - contour split across chunks */
  752.                     /* put # in case user wants to use it for a plot */
  753.                     /* double blank line to allow plot ... index ... */
  754.                     fprintf(outfile, "\n# Contour %d, label:%s\n", number++, c->label);
  755.                 for ( ; --count >= 0 ; ++p)
  756.                     fprintf(outfile, "%f %f %f\n", p->x, p->y, p->z);
  757.                 putc('\n', outfile);  /* blank line between segments of same contour */
  758.                 c=c->next;
  759.             }
  760.         }            
  761.     }
  762.     fflush(outfile);
  763. }
  764.  
  765.  
  766.  
  767. #define SET_DUMMY_RANGE(AXIS) \
  768. do{\
  769.  if (parametric || polar) { \
  770.   t_min = tmin; t_max = tmax;\
  771.  } else if (log_array[AXIS]) {\
  772.   if (min_array[AXIS] <= 0.0 || max_array[AXIS] <= 0.0)\
  773.    int_error("x/x2 range must be greater than 0 for log scale!", NO_CARET);\
  774.   t_min = log(min_array[AXIS])/log_base_array[AXIS]; t_max = log(max_array[AXIS])/log_base_array[AXIS];\
  775.  } else {\
  776.   t_min = min_array[AXIS]; t_max = max_array[AXIS];\
  777.  }\
  778.  t_step = (t_max - t_min) / (samples - 1); \
  779. }while(0)
  780.  
  781.  
  782. /*
  783.  * This parses the splot command after any range specifications. To support
  784.  * autoscaling on the x/z axis, we want any data files to define the x/y
  785.  * range, then to plot any functions using that range. We thus parse the
  786.  * input twice, once to pick up the data files, and again to pick up the
  787.  * functions. Definitions are processed twice, but that won't hurt.
  788.  * div - okay, it doesn't hurt, but every time an option as added for
  789.  * datafiles, code to parse it has to be added here. Change so that
  790.  * we store starting-token in the plot structure.
  791.  */
  792. static void eval_3dplots()
  793. {
  794.     register int    i, j;
  795.     register struct surface_points *this_plot=NULL, **tp_3d_ptr;
  796.     register int    start_token, end_token;
  797.     register int    begin_token;
  798.     TBOOLEAN        some_data_files = FALSE, some_functions=FALSE;
  799.     int             plot_num, line_num, point_num, crnt_param = 0;    /* 0=z, 1=y, 2=x */
  800.     char           *xtitle;
  801.     char           *ytitle;
  802.  
  803.     /* Reset first_3dplot. This is usually done at the end of this function.
  804.        If there is an error within this function, the memory is left allocated,
  805.        since we cannot call sp_free if the list is incomplete */
  806.     first_3dplot=NULL;
  807.  
  808.     /* put stuff into arrays to simplify access */
  809.     INIT_ARRAYS(FIRST_X_AXIS, xmin, xmax, autoscale_lx, is_log_x, base_log_x, log_base_log_x, 0);
  810.     INIT_ARRAYS(FIRST_Y_AXIS, ymin, ymax, autoscale_ly, is_log_y, base_log_y, log_base_log_y, 0);
  811.     INIT_ARRAYS(FIRST_Z_AXIS, zmin, zmax, autoscale_lz, is_log_z, base_log_z, log_base_log_z, 1);
  812.  
  813.     x_axis=FIRST_X_AXIS;
  814.     y_axis=FIRST_Y_AXIS;
  815.     z_axis=FIRST_Z_AXIS;
  816.  
  817.     tp_3d_ptr = &(first_3dplot);
  818.     plot_num = 0;
  819.     line_num = 0;        /* default line type */
  820.     point_num = 0;        /* default point type */
  821.  
  822.     xtitle = NULL;
  823.     ytitle = NULL;
  824.  
  825.     begin_token = c_token;
  826.  
  827.     /*** First Pass: Read through data files ***/
  828.     /*
  829.      * This pass serves to set the x/yranges and to parse the command, as
  830.      * well as filling in every thing except the function data. That is done
  831.      * after the x/yrange is defined.
  832.      */
  833.     while (TRUE) {
  834.     if (END_OF_COMMAND)
  835.         int_error("function to plt3d expected", c_token);
  836.  
  837.     start_token = c_token;
  838.  
  839.     if (is_definition(c_token)) {
  840.         define();
  841.     } else {
  842.         int specs;
  843.         if (isstring(c_token)) {    /* data file to plot */
  844.  
  845.         if (parametric && crnt_param != 0)
  846.             int_error("previous parametric function not fully specified",
  847.                   c_token);
  848.  
  849.         if (!some_data_files) {
  850.             if (autoscale_lx & 1) {
  851.             min_array[FIRST_X_AXIS] = VERYLARGE;
  852.             }
  853.             if (autoscale_lx & 2) {
  854.             max_array[FIRST_X_AXIS] = -VERYLARGE;
  855.             }
  856.             if (autoscale_ly & 1) {
  857.             min_array[FIRST_Y_AXIS] = VERYLARGE;
  858.             }
  859.             if (autoscale_ly & 2) {
  860.             max_array[FIRST_Y_AXIS] = -VERYLARGE;
  861.             }
  862.         }
  863.         some_data_files = TRUE;
  864.  
  865.         if (*tp_3d_ptr)
  866.             this_plot = *tp_3d_ptr;
  867.         else {        /* no memory malloc()'d there yet */
  868.             /* Allocate enough isosamples and samples */
  869.             this_plot = sp_alloc(0, 0, 0, 0);
  870.             *tp_3d_ptr = this_plot;
  871.         }
  872.  
  873.         this_plot->plot_type = DATA3D;
  874.         this_plot->plot_style = data_style;
  875.  
  876.         specs = df_open(3);
  877.         /* parses all datafile-specific modifiers */
  878.         /* we will load the data after parsing title,with,... */
  879.         this_plot->token = end_token = c_token-1;  /* for capture to key */
  880.         /* this_plot->token is temporary, for errors in get_3ddata() */
  881.  
  882.         if (datatype[FIRST_X_AXIS]==TIME)
  883.         {
  884.             if (specs<3)
  885.                 int_error("Need full using spec for x time data", c_token);
  886.             df_timecol[0]=1;
  887.         }
  888.  
  889.         if (datatype[FIRST_Y_AXIS]==TIME)
  890.         {
  891.             if (specs<3)
  892.                 int_error("Need full using spec for y time data", c_token);
  893.             df_timecol[1]=1;
  894.         }
  895.         
  896.         if (datatype[FIRST_Z_AXIS]==TIME)
  897.         {
  898.             if (specs<3)
  899.                 df_timecol[0]=1;
  900.             else
  901.                 df_timecol[2]=1;
  902.         }
  903.  
  904.         
  905.         } else {        /* function to plot */
  906.         ++plot_num;
  907.         if (parametric)    /* Rotate between x/y/z axes */
  908.             crnt_param = (crnt_param + 2) % 3;  /* +2 same as -1, but beats -ve problem */
  909.         if (*tp_3d_ptr) {
  910.             this_plot = *tp_3d_ptr;
  911.             if (!hidden3d)
  912.             sp_replace(this_plot, samples_1, iso_samples_1,
  913.                                    samples_2, iso_samples_2);
  914.             else
  915.             sp_replace(this_plot, iso_samples_1, 0,
  916.                                    0, iso_samples_2);
  917.         } else {    /* no memory malloc()'d there yet */
  918.             /* Allocate enough isosamples and samples */
  919.             if (!hidden3d)
  920.             this_plot = sp_alloc(samples_1, iso_samples_1,
  921.                                              samples_2, iso_samples_2);
  922.             else
  923.             this_plot = sp_alloc(iso_samples_1, 0,
  924.                                              0, iso_samples_2);
  925.             *tp_3d_ptr = this_plot;
  926.         }
  927.  
  928.         this_plot->plot_type = FUNC3D;
  929.         this_plot->has_grid_topology = TRUE;
  930.         this_plot->plot_style = func_style;
  931.         this_plot->num_iso_read = iso_samples_2;
  932.         dummy_func = &plot_func;
  933.         plot_func.at = temp_at();
  934.         dummy_func=NULL;
  935.         /* ignore it for now */
  936.         some_functions=TRUE;
  937.         end_token = c_token - 1;
  938.     }  /* end of IS THIS A FILE OR A FUNC block */
  939.  
  940.         if (this_plot->title) {
  941.             free(this_plot->title);
  942.             this_plot->title=NULL;
  943.         }
  944.  
  945.     
  946.         if (almost_equals(c_token, "t$itle")) {
  947.             if (!isstring(++c_token))
  948.                 int_error("Expected title", c_token);
  949.             m_quote_capture(&(this_plot->title), c_token, c_token);
  950.             ++c_token;
  951.         } else if (almost_equals(c_token,"not$itle")) {
  952.             this_plot->title=NULL;
  953.             ++c_token;
  954.         } else {
  955.                 m_capture(&(this_plot->title), start_token, end_token);
  956.         }
  957.  
  958.         this_plot->line_type = line_num;
  959.         this_plot->point_type = point_num;
  960.  
  961.         if (almost_equals(c_token, "w$ith")) {
  962.             this_plot->plot_style = get_style();
  963.         }
  964.         if (!equals(c_token, ",") && !END_OF_COMMAND) {
  965.             struct value    t;
  966.             this_plot->line_type = (int) real(const_express(&t)) - 1;
  967.         }
  968.         if (!equals(c_token, ",") && !END_OF_COMMAND) {
  969.             struct value    t;
  970.             this_plot->point_type = (int) real(const_express(&t)) - 1;
  971.         }
  972.  
  973.  
  974.         if ((this_plot->plot_style == POINTSTYLE) ||
  975.         (this_plot->plot_style == LINESPOINTS) ||
  976.                 (this_plot->plot_style == YERRORBARS) ||
  977.                 (this_plot->plot_style == XERRORBARS) ||
  978.                 (this_plot->plot_style == XYERRORBARS) ||
  979.                 (this_plot->plot_style == BOXXYERROR))
  980.         if (crnt_param == 0)
  981.             point_num +=
  982.             1 + (draw_contour != 0)
  983.             + (hidden3d != 0);
  984.         if (crnt_param == 0)
  985.         line_num += 1 + (draw_contour != 0)
  986.             + (hidden3d != 0);
  987.     }
  988.  
  989.  
  990.     /* now get the data... having to think hard here...
  991.      * first time through, we fill in this_plot. For second
  992.      * surface in file, we have to allocate another surface
  993.      * struct. BUT we may allocate this store only to
  994.      * find that it is merely some blank lines at end of file
  995.      * tp_3d_ptr is still pointing at next field of prev. plot,
  996.          * before :    prev_or_first -> this_plot -> possible_preallocated_store
  997.      *                tp_3d_ptr--^
  998.      * after  :    prev_or_first -> first -> second -> last -> possibly_more_store
  999.      *                                        tp_3d_ptr ----^
  1000.      * if file is empty, tp_3d_ptr is not moved. this_plot continues
  1001.      * to point at allocated storage, but that will be reused later
  1002.      */
  1003.  
  1004.     assert(this_plot == *tp_3d_ptr);
  1005.  
  1006.     if (this_plot->plot_type == DATA3D) {
  1007.         /* remember settings for second surface in file */
  1008.         int this_line=this_plot->line_type;
  1009.         int this_point=this_plot->point_type;
  1010.         int this_style=this_plot->plot_style;
  1011.         int this_token=this_plot->token;
  1012.         while (!df_eof)
  1013.         {
  1014.             this_plot = *tp_3d_ptr;
  1015.             assert(this_plot != NULL);
  1016.  
  1017.             /* dont move tp_3d_ptr until we are sure we
  1018.              * have read a surface
  1019.              */
  1020.             this_plot->token=this_token; /* used by get_3ddata() */
  1021.             get_3ddata(this_plot);
  1022.             this_plot->token=c_token; /* for second pass */
  1023.  
  1024.             if (this_plot->num_iso_read==0)
  1025.                 /* probably df_eof, in which case we
  1026.                  * will leave loop. if not eof, then
  1027.                  * how come we got no surface ? - retry
  1028.                  * in neither case do we update tp_3d_ptr
  1029.                  */
  1030.                 continue;
  1031.  
  1032.             /* okay, we have read a surface */
  1033.             ++plot_num;
  1034.             tp_3d_ptr = &(this_plot->next_sp);
  1035.             if (df_eof)
  1036.                 break;
  1037.  
  1038.             /* there might be another surface */
  1039.             if ( (this_plot=*tp_3d_ptr) != NULL ) {
  1040.                 if (this_plot->title) {
  1041.                     free(this_plot->title);
  1042.                     this_plot->title = NULL;
  1043.                 }
  1044.             } else {
  1045.                 /* Allocate enough isosamples and samples */
  1046.                 this_plot = *tp_3d_ptr = sp_alloc(0, 0, 0, 0);
  1047.             }
  1048.             this_plot->plot_type=DATA3D;
  1049.             this_plot->line_type=this_line;
  1050.             this_plot->point_type=this_point;
  1051.             this_plot->plot_style=this_style;
  1052.         }
  1053.         df_close();
  1054.     } else { /* not a data file */
  1055.         tp_3d_ptr = &(this_plot->next_sp);
  1056.         this_plot->token = c_token;  /* store for second pass */
  1057.     }
  1058.  
  1059.     /* DO NOT USE this_plot AFTER HERE - it may point at an unused data plot */
  1060.  
  1061.     if (equals(c_token, ","))
  1062.         c_token++;
  1063.     else
  1064.         break;
  1065.     }
  1066.  
  1067.     if (parametric && crnt_param != 0)
  1068.     int_error("parametric function not fully specified", NO_CARET);
  1069.  
  1070.  
  1071.     /*** Second Pass: Evaluate the functions ***/
  1072.     /*
  1073.      * Everything is defined now, except the function data. We expect no
  1074.      * syntax errors, etc, since the above parsed it all. This makes the code
  1075.      * below simpler. If autoscale_ly, the yrange may still change.
  1076.      * - eh ?  - z can still change.  x/y/z can change if we are parametric ??
  1077.      */
  1078.  
  1079.   if (some_functions) {
  1080.  
  1081.     /* I've changed the controlled variable in fn plots to u_min etc since
  1082.      *it's easier for me to think parametric - 'normal' plot is after all
  1083.      * a special case. I was confused about x_min being both minimum of
  1084.      * x values found, and starting value for fn plots.
  1085.      */
  1086.     register double u_min, u_max, u_step, v_min, v_max, v_step;
  1087.     double uisodiff, visodiff;
  1088.  
  1089.     if (!parametric) {
  1090.  
  1091.     /* give error if xrange badly set from missing datafile error
  1092.      * parametric fn can still set ranges
  1093.      * if there are no fns, we'll report it later as 'nothing to plot'
  1094.      */
  1095.  
  1096.     if (min_array[FIRST_X_AXIS] == VERYLARGE || max_array[FIRST_X_AXIS] == -VERYLARGE) {
  1097.         int_error("x range is invalid", c_token);
  1098.     }
  1099.     if (min_array[FIRST_Y_AXIS] == VERYLARGE || max_array[FIRST_Y_AXIS] == -VERYLARGE) {
  1100.         int_error("y range is invalid", c_token);
  1101.     }
  1102.  
  1103.     /* check that xmin -> xmax is not too small */
  1104.     FIXUP_RANGE(FIRST_X_AXIS, x);
  1105.     FIXUP_RANGE(FIRST_Y_AXIS, y);
  1106.     }
  1107.  
  1108.    if (parametric && !some_data_files) {
  1109.     /* parametric fn can still change x/y range */
  1110.     if (autoscale_lx & 1)
  1111.         min_array[FIRST_X_AXIS] = VERYLARGE;
  1112.     if (autoscale_lx & 2)
  1113.         max_array[FIRST_X_AXIS] = -VERYLARGE;
  1114.     if (autoscale_ly & 1)
  1115.         min_array[FIRST_Y_AXIS] = VERYLARGE;
  1116.     if (autoscale_ly & 2)
  1117.         max_array[FIRST_Y_AXIS] = -VERYLARGE;
  1118.     }
  1119.  
  1120.  
  1121.     if (parametric) {
  1122.     u_min=umin;
  1123.     u_max=umax;
  1124.     v_min=vmin;
  1125.     v_max=vmax;
  1126.     } else {
  1127.     if (is_log_x) {
  1128.         if (min_array[FIRST_X_AXIS] <= 0.0 || max_array[FIRST_X_AXIS] <= 0.0)
  1129.         int_error("x range must be greater than 0 for log scale!", NO_CARET);
  1130.         u_min = log(min_array[FIRST_X_AXIS])/log_base_log_x;
  1131.         u_max = log(max_array[FIRST_X_AXIS])/log_base_log_x;
  1132.         } else {
  1133.         u_min = min_array[FIRST_X_AXIS];
  1134.         u_max = max_array[FIRST_X_AXIS];
  1135.         }
  1136.  
  1137.     if (is_log_y) {
  1138.         if (min_array[FIRST_Y_AXIS] <= 0.0 || max_array[FIRST_Y_AXIS] <= 0.0)
  1139.         int_error("y range must be greater than 0 for log scale!", NO_CARET);
  1140.         v_min = log(min_array[FIRST_Y_AXIS])/log_base_log_y;
  1141.         v_max = log(max_array[FIRST_Y_AXIS])/log_base_log_y;
  1142.     } else {
  1143.         v_min = min_array[FIRST_Y_AXIS];
  1144.         v_max = max_array[FIRST_Y_AXIS];
  1145.     }
  1146.     }
  1147.  
  1148.  
  1149.     if (samples_1 < 2 || samples_2 < 2 || iso_samples_1 < 2 || iso_samples_2 < 2)
  1150.     int_error("samples or iso_samples < 2. Must be at least 2.", NO_CARET);
  1151.  
  1152.     if (this_plot && this_plot->has_grid_topology && hidden3d) {
  1153.     u_step = (u_max - u_min) / (iso_samples_1 - 1);
  1154.     v_step = (v_max - v_min) / (iso_samples_2 - 1);
  1155.     } else {
  1156.     u_step = (u_max - u_min) / (samples_1 - 1);
  1157.     v_step = (v_max - v_min) / (samples_2 - 1);
  1158.     }
  1159.     uisodiff = (u_max - u_min) / (iso_samples_1 - 1);
  1160.     visodiff = (v_max - v_min) / (iso_samples_2 - 1);
  1161.  
  1162.     this_plot = first_3dplot;
  1163.     c_token = begin_token;    /* start over */
  1164.  
  1165.     /* Read through functions */
  1166.     while (TRUE) {
  1167.     if (is_definition(c_token)) {
  1168.         define();
  1169.     } else {
  1170.         if (!isstring(c_token)) { /* func to plot */
  1171.         struct iso_curve *this_iso = this_plot->iso_crvs;
  1172.         struct coordinate GPHUGE *points = this_iso->points;
  1173.         int num_sam_to_use, num_iso_to_use;
  1174.  
  1175.         if (parametric)
  1176.             crnt_param = (crnt_param + 2) % 3;
  1177.         dummy_func = &plot_func;
  1178.         plot_func.at = temp_at();    /* reparse function */
  1179.         dummy_func=NULL;
  1180.         num_iso_to_use = iso_samples_2;
  1181.         if (!(this_plot->has_grid_topology && hidden3d))
  1182.             num_sam_to_use = samples_1;
  1183.         else
  1184.             num_sam_to_use = iso_samples_1;
  1185.  
  1186.         for (j = 0; j < num_iso_to_use; j++) {
  1187.             double y = v_min + j * visodiff;
  1188.             /* if (is_log_y) PEM fix logscale y axis */
  1189.             /* y = pow(log_base_log_y,y); 26-Sep-89 */
  1190.             /* parametric => NOT a log quantity (?) */
  1191.             (void) Gcomplex(&plot_func.dummy_values[1],
  1192.                    !parametric && is_log_y ? pow(base_log_y, y) : y,
  1193.                    0.0);
  1194.  
  1195.             for (i = 0; i < num_sam_to_use; i++) {
  1196.             double x = u_min + i * u_step;
  1197.             struct value a;
  1198.             double temp;
  1199.             
  1200.             /* if (is_log_x) PEM fix logscale x axis */
  1201.             /* x = pow(base_log_x,x); 26-Sep-89 */
  1202.             /* parametric => NOT a log quantity (?) */
  1203.             (void) Gcomplex(&plot_func.dummy_values[0],
  1204.                        !parametric && is_log_x ? pow(base_log_x, x) : x,
  1205.                        0.0);
  1206.  
  1207.             points[i].x = x;
  1208.             points[i].y = y;
  1209.  
  1210.             evaluate_at(plot_func.at, &a);
  1211.  
  1212.             if (undefined || (fabs(imag(&a)) > zero)) {
  1213.                 points[i].type = UNDEFINED;
  1214.                 continue;
  1215.             }
  1216.             temp = real(&a);
  1217.  
  1218.             points[i].type=INRANGE;
  1219.             STORE_WITH_LOG_AND_FIXUP_RANGE(points[i].z, temp, points[i].type,
  1220.                 crnt_param, NOOP, NOOP);
  1221.  
  1222.             }
  1223.             this_iso->p_count = num_sam_to_use;
  1224.             this_iso = this_iso->next;
  1225.             points = this_iso? this_iso->points: NULL;
  1226.         }
  1227.  
  1228.         if (!(this_plot->has_grid_topology && hidden3d)) {
  1229.             num_iso_to_use = iso_samples_1;
  1230.             num_sam_to_use = samples_2;
  1231.             for (i = 0; i < num_iso_to_use; i++) {
  1232.             double x = u_min + i * uisodiff;
  1233.             /* if (is_log_x) PEM fix logscale x axis */
  1234.             /* x = pow(base_log_x,x); 26-Sep-89 */
  1235.             /* if parametric, no logs involved - 3.6 */
  1236.             (void) Gcomplex(&plot_func.dummy_values[0],
  1237.                        (!parametric && is_log_x) ? pow(base_log_x, x) : x,
  1238.                        0.0);
  1239.  
  1240.             for (j = 0; j < num_sam_to_use; j++) {
  1241.                 double y = v_min + j * v_step;
  1242.                 struct value a;
  1243.                 double temp;
  1244.                 /* if (is_log_y) PEM fix logscale y axis */
  1245.                 /* y = pow(base_log_y,y); 26-Sep-89 */
  1246.                 (void) Gcomplex(&plot_func.dummy_values[1],
  1247.                        (!parametric && is_log_y) ? pow(base_log_y, y) : y,
  1248.                        0.0);
  1249.  
  1250.                 points[j].x = x;
  1251.                 points[j].y = y;
  1252.  
  1253.                 evaluate_at(plot_func.at, &a);
  1254.  
  1255.                 if (undefined || (fabs(imag(&a)) > zero)) {
  1256.                 points[j].type = UNDEFINED;
  1257.                 continue;
  1258.                 }
  1259.                 temp = real(&a);
  1260.  
  1261.                 points[j].type=INRANGE;
  1262.                 STORE_WITH_LOG_AND_FIXUP_RANGE(points[j].z, temp, points[j].type,
  1263.                     crnt_param, NOOP, NOOP);
  1264.             }
  1265.             this_iso->p_count = num_sam_to_use;
  1266.             this_iso = this_iso->next;
  1267.             points = this_iso ? this_iso->points : NULL;
  1268.             }
  1269.         }
  1270.         }  /* end of ITS A FUNCTION TO PLOT */
  1271.         c_token = this_plot->token;  /* we saved it from first pass */
  1272.  
  1273.         do
  1274.         this_plot = this_plot->next_sp;
  1275.         while (this_plot && this_plot->token == c_token);  /* one data file can make several plots */
  1276.     }
  1277.  
  1278.     if (equals(c_token, ","))
  1279.         c_token++;
  1280.     else
  1281.         break;
  1282.     }
  1283.  
  1284.     if (parametric) {
  1285.     /* Now actually fix the plot triplets to be single plots. */
  1286.     parametric_3dfixup(first_3dplot, &plot_num);
  1287.  
  1288.     FIXUP_RANGE(FIRST_X_AXIS, x);
  1289.     FIXUP_RANGE(FIRST_Y_AXIS, y);
  1290.     }
  1291.  
  1292.   }  /* some functions */
  1293.  
  1294.  
  1295.     /* if first_3dplot is NULL, we have no functions or data at all. This can
  1296.        happen, if you type "splot x=5", since x=5 is a variable assignment */
  1297.  
  1298.     if(plot_num == 0 || first_3dplot==NULL) {
  1299.     int_error("no functions or data to plot", c_token);
  1300.     }
  1301.  
  1302.     if (min_array[FIRST_X_AXIS] ==  VERYLARGE ||
  1303.         max_array[FIRST_X_AXIS] == -VERYLARGE ||
  1304.         min_array[FIRST_Y_AXIS] ==  VERYLARGE ||
  1305.         max_array[FIRST_Y_AXIS] == -VERYLARGE ||
  1306.         min_array[FIRST_Z_AXIS] ==  VERYLARGE ||
  1307.         max_array[FIRST_Z_AXIS] == -VERYLARGE)
  1308.        int_error("All points undefined", NO_CARET);
  1309.  
  1310.     FIXUP_RANGE(FIRST_Z_AXIS, z);
  1311.  
  1312.     FIXUP_RANGE_FOR_LOG(FIRST_X_AXIS, x);
  1313.     FIXUP_RANGE_FOR_LOG(FIRST_Y_AXIS, y);
  1314.     FIXUP_RANGE_FOR_LOG(FIRST_Z_AXIS, z);
  1315.  
  1316.     /* last parameter should take plot size into effect...
  1317.      * probably needs to be moved to graph3d.c
  1318.      * in the meantime, a value of 20 gives same behaviour
  1319.      * as 3.5 which will do for the moment
  1320.      */
  1321.  
  1322.     if (xtics) setup_tics(FIRST_X_AXIS, &xticdef, xformat, 20);
  1323.     if (ytics) setup_tics(FIRST_Y_AXIS, &yticdef, yformat, 20);
  1324.     if (ztics) setup_tics(FIRST_Z_AXIS, &zticdef, zformat, 20);
  1325.  
  1326. #define WRITEBACK(axis,min,max) \
  1327. if(range_flags[axis]&RANGE_WRITEBACK) \
  1328.   {if (auto_array[axis]&1) min=min_array[axis]; \
  1329.    if (auto_array[axis]&2) max=max_array[axis]; \
  1330.   }
  1331.  
  1332.     WRITEBACK(FIRST_X_AXIS,xmin,xmax)
  1333.     WRITEBACK(FIRST_Y_AXIS,ymin,ymax)
  1334.     WRITEBACK(FIRST_Z_AXIS,zmin,zmax)
  1335.  
  1336.     if (plot_token != -1) {
  1337.     /* note that m_capture also frees the old replot_line */
  1338.     m_capture(&replot_line, plot_token, c_token);
  1339.     plot_token = -1;
  1340.     }
  1341.  
  1342.     if(plot_num==0 || first_3dplot==NULL) {
  1343.     int_error("no functions or data to plot", c_token);
  1344.     }
  1345.  
  1346.     /* Creates contours if contours are to be plotted as well. */
  1347.     if (draw_contour) {
  1348.     for (this_plot = first_3dplot, i = 0;
  1349.          i < plot_num;
  1350.          this_plot = this_plot->next_sp, i++) {
  1351.         if (this_plot->contours) {
  1352.         struct gnuplot_contours *cntr, *cntrs = this_plot->contours;
  1353.  
  1354.         while (cntrs) {
  1355.             cntr = cntrs;
  1356.             cntrs = cntrs->next;
  1357.             free(cntr->coords);
  1358.             free(cntr);
  1359.         }
  1360.         }
  1361.         /* Make sure this one can be contoured. */
  1362.         if (!this_plot->has_grid_topology) {
  1363.         this_plot->contours = NULL;
  1364.         fprintf(stderr,"Notice: cannot contour non grid data!\n");
  1365.         /* changed from int_error by recommendation of rkc@xn.ll.mit.edu */
  1366.         }
  1367.         else if (this_plot->plot_type == DATA3D)
  1368.         this_plot->contours = contour(
  1369.                          this_plot->num_iso_read,
  1370.                          this_plot->iso_crvs,
  1371.                          contour_levels, contour_pts,
  1372.                          contour_kind, contour_order,
  1373.                          levels_kind, levels_list);
  1374.         else
  1375.         this_plot->contours = contour(iso_samples_2,
  1376.                           this_plot->iso_crvs,
  1377.                           contour_levels, contour_pts,
  1378.                           contour_kind, contour_order,
  1379.                           levels_kind, levels_list);
  1380.     }
  1381.     }
  1382.     if (strcmp(term->name, "table") == 0)
  1383.     print_3dtable(plot_num);
  1384.     else
  1385.     do_3dplot(first_3dplot, plot_num);
  1386.     sp_free(first_3dplot);
  1387.     first_3dplot = NULL;
  1388. }
  1389.  
  1390. static void 
  1391. parametric_3dfixup(start_plot, plot_num)
  1392.     struct surface_points *start_plot;
  1393.     int            *plot_num;
  1394. /*
  1395.  * The hardest part of this routine is collapsing the FUNC plot types in the
  1396.  * list (which are gauranteed to occur in (x,y,z) triplets while preserving
  1397.  * the non-FUNC type plots intact.  This means we have to work our way
  1398.  * through various lists.  Examples (hand checked):
  1399.  * start_plot:F1->F2->F3->NULL ==> F3->NULL
  1400.  * start_plot:F1->F2->F3->F4->F5->F6->NULL ==> F3->F6->NULL
  1401.  * start_plot:F1->F2->F3->D1->D2->F4->F5->F6->D3->NULL ==>
  1402.  * F3->D1->D2->F6->D3->NULL
  1403.  */
  1404. {
  1405. /*
  1406.  * I initialized *free_list with NULL, because my compiler warns some lines
  1407.  * later that it might be uninited. The code however seems to not access that
  1408.  * line in that case, but if I'm right, my change is OK and if not, this is a
  1409.  * serious bug in the code.
  1410.  *
  1411.  * x and y ranges now fixed in eval_3dplots
  1412.  */
  1413.     struct surface_points *xp, *new_list, *free_list = NULL;
  1414.     struct surface_points **last_pointer=&new_list;
  1415.  
  1416.     int             i, tlen, surface;
  1417.     char           *new_title;
  1418.  
  1419.     /*
  1420.      * Ok, go through all the plots and move FUNC3D types together.  Note:
  1421.      * this originally was written to look for a NULL next pointer, but
  1422.      * gnuplot wants to be sticky in grabbing memory and the right number of
  1423.      * items in the plot list is controlled by the plot_num variable.
  1424.      * 
  1425.      * Since gnuplot wants to do this sticky business, a free_list of
  1426.      * surface_points is kept and then tagged onto the end of the plot list
  1427.      * as this seems more in the spirit of the original memory behavior than
  1428.      * simply freeing the memory.  I'm personally not convinced this sort of
  1429.      * concern is worth it since the time spent computing points seems to
  1430.      * dominate any garbage collecting that might be saved here...
  1431.      */
  1432.     new_list = xp = start_plot;
  1433.     for (surface = 0; surface < *plot_num; surface++) {
  1434.       if (xp->plot_type == FUNC3D) {
  1435.     struct surface_points *yp = xp->next_sp;
  1436.     struct surface_points *zp = yp->next_sp;
  1437.  
  1438.     /* Here's a FUNC3D parametric function defined as three parts.
  1439.      * Go through all the points and assign the x's and y's from xp and
  1440.      * yp to zp. min/max already done
  1441.      */
  1442.     struct iso_curve *xicrvs = xp->iso_crvs;
  1443.     struct iso_curve *yicrvs = yp->iso_crvs;
  1444.     struct iso_curve *zicrvs = zp->iso_crvs;
  1445.  
  1446.     (*plot_num) -= 2;
  1447.  
  1448.     assert(INRANGE < OUTRANGE && OUTRANGE < UNDEFINED);
  1449.  
  1450.     while (zicrvs) {
  1451.         struct coordinate GPHUGE *xpoints = xicrvs->points, GPHUGE *ypoints = yicrvs->points, GPHUGE *zpoints = zicrvs->points;
  1452.         for (i = 0; i < zicrvs->p_count; ++i) {
  1453.         zpoints[i].x = xpoints[i].z;
  1454.         zpoints[i].y = ypoints[i].z;
  1455.         if (zpoints[i].type < xpoints[i].type) zpoints[i].type = xpoints[i].type;
  1456.         if (zpoints[i].type < ypoints[i].type) zpoints[i].type = ypoints[i].type;
  1457.  
  1458.         }
  1459.         xicrvs = xicrvs->next;
  1460.         yicrvs = yicrvs->next;
  1461.         zicrvs = zicrvs->next;
  1462.     }
  1463.  
  1464.     /* Ok, fix up the title to include xp and yp plots. */
  1465.     if ((xp->title && xp->title[0] != '\0') ||
  1466.         (yp->title && yp->title[0] != '\0')) {
  1467.         tlen = (xp->title ? strlen(xp->title) : 0) +
  1468.         (yp->title ? strlen(yp->title) : 0) +
  1469.         (zp->title ? strlen(zp->title) : 0) + 5;
  1470.         new_title = alloc((unsigned long) tlen, "string");
  1471.         new_title[0] = 0;
  1472.         if (xp->title) {
  1473.         strcat(new_title, xp->title);
  1474.         strcat(new_title, ", ");    /* + 2 */
  1475.         }
  1476.         if (yp->title) {
  1477.         strcat(new_title, yp->title);
  1478.         strcat(new_title, ", ");    /* + 2 */
  1479.         }
  1480.         if (zp->title) {
  1481.         strcat(new_title, zp->title);
  1482.         }
  1483.         free(zp->title);
  1484.         zp->title = new_title;
  1485.     }
  1486.  
  1487.     /* add xp and yp to head of free list */
  1488.     assert(xp->next_sp == yp);
  1489.     yp->next_sp = free_list;
  1490.     free_list = xp;
  1491.  
  1492.     /* add zp to tail of new_list */
  1493.     *last_pointer = zp;
  1494.     last_pointer = &(zp->next_sp);
  1495.     xp=zp->next_sp;
  1496.       } else {  /* its a data plot */
  1497.     assert(*last_pointer == xp);  /* think this is true ! */
  1498.     last_pointer = &(xp->next_sp);
  1499.     xp=xp->next_sp;
  1500.       }
  1501.     }
  1502.  
  1503.     /* Ok, append free list and write first_plot */
  1504.     *last_pointer = free_list;
  1505.     first_3dplot=new_list;
  1506. }
  1507.