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