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