PLplot  5.15.0
plot3d.c
Go to the documentation of this file.
1 // Copyright (C) 2004-2014 Alan W. Irwin
6 // Copyright (C) 2004 Joao Cardoso
7 // Copyright (C) 2004 Andrew Ross
8 //
9 // This file is part of PLplot.
10 //
11 // PLplot is free software; you can redistribute it and/or modify
12 // it under the terms of the GNU Library General Public License as published
13 // by the Free Software Foundation; either version 2 of the License, or
14 // (at your option) any later version.
15 //
16 // PLplot is distributed in the hope that it will be useful,
17 // but WITHOUT ANY WARRANTY; without even the implied warranty of
18 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
19 // GNU Library General Public License for more details.
20 //
21 // You should have received a copy of the GNU Library General Public License
22 // along with PLplot; if not, write to the Free Software
23 // Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
24 //
25 
26 #include "plplotP.h"
27 
28 // Internal constants
29 
30 #define BINC 50 // Block size for memory allocation
31 
32 static PLINT pl3mode = 0; // 0 3d solid; 1 mesh plot
33 static PLINT pl3upv = 1; // 1 update view; 0 no update
34 
35 static PLINT zbflg = 0, zbcol;
36 static PLFLT zbtck, zbwidth;
37 
38 static PLINT *oldhiview = NULL;
39 static PLINT *oldloview = NULL;
40 static PLINT *newhiview = NULL;
41 static PLINT *newloview = NULL;
42 static PLINT *utmp = NULL;
43 static PLINT *vtmp = NULL;
44 static PLFLT *ctmp = NULL;
45 
46 static PLINT mhi, xxhi, newhisize;
47 static PLINT mlo, xxlo, newlosize;
48 
49 // Light source for shading
50 static PLFLT xlight, ylight, zlight;
51 static PLINT falsecolor = 0;
52 static PLFLT fc_minz, fc_maxz;
53 
54 // Prototypes for static functions
55 
56 static void plgrid3( PLFLT );
57 static void plnxtv( PLINT *, PLINT *, PLFLT*, PLINT, PLINT );
58 static void
59 plside3( PLFLT_VECTOR x, PLFLT_VECTOR y, PLF2OPS zops, PLPointer zp, PLINT nx, PLINT ny, PLINT opt );
60 static void
61 plt3zz( PLINT x0, PLINT y0, PLINT dx, PLINT dy, PLINT flag, PLINT *init,
62  PLFLT_VECTOR x, PLFLT_VECTOR y, PLF2OPS zops, PLPointer zp, PLINT nx, PLINT ny,
63  PLINT *u, PLINT *v, PLFLT* c );
64 static void plnxtvhi( PLINT *, PLINT *, PLFLT*, PLINT, PLINT );
65 static void plnxtvlo( PLINT *, PLINT *, PLFLT*, PLINT, PLINT );
66 static void plnxtvhi_draw( PLINT *u, PLINT *v, PLFLT* c, PLINT n );
67 
68 static void savehipoint( PLINT, PLINT );
69 static void savelopoint( PLINT, PLINT );
70 static void swaphiview( void );
71 static void swaploview( void );
72 static void myexit( PLCHAR_VECTOR );
73 static void myabort( PLCHAR_VECTOR );
74 static void freework( void );
75 static int plabv( PLINT, PLINT, PLINT, PLINT, PLINT, PLINT );
76 static void pl3cut( PLINT, PLINT, PLINT, PLINT, PLINT,
77  PLINT, PLINT, PLINT, PLINT *, PLINT * );
78 static PLFLT plGetAngleToLight( PLFLT* x, PLFLT* y, PLFLT* z );
79 static void plP_draw3d( PLINT x, PLINT y, PLFLT *c, PLINT j, PLINT move );
80 //static void plxyindexlimits( PLINT instart, PLINT inn,
81 // PLINT *inarray_min, PLINT *inarray_max,
82 // PLINT *outstart, PLINT *outn, PLINT outnmax,
83 // PLINT *outarray_min, PLINT *outarray_max );
84 
85 
86 // #define MJL_HACK 1
87 #if MJL_HACK
88 static void plP_fill3( PLINT x0, PLINT y0, PLINT x1, PLINT y1,
89  PLINT x2, PLINT y2, PLINT j );
90 static void plP_fill4( PLINT x0, PLINT y0, PLINT x1, PLINT y1,
91  PLINT x2, PLINT y2, PLINT x3, PLINT y3, PLINT j );
92 #endif
93 
94 //--------------------------------------------------------------------------
95 // void plsetlightsource(x, y, z)
96 //
97 // Sets the position of the light source.
98 //--------------------------------------------------------------------------
99 
100 void
101 c_pllightsource( PLFLT x, PLFLT y, PLFLT z )
102 {
103  xlight = x;
104  ylight = y;
105  zlight = z;
106 }
107 
108 //--------------------------------------------------------------------------
109 // void plmesh(x, y, z, nx, ny, opt)
110 //
111 // Plots a mesh representation of the function z[x][y]. The x values
112 // are stored as x[0..nx-1], the y values as y[0..ny-1], and the
113 // z values are in the 2-d array z[][]. The integer "opt" specifies:
114 // see plmeshc() below.
115 //--------------------------------------------------------------------------
116 
117 void
118 c_plmesh( PLFLT_VECTOR x, PLFLT_VECTOR y, PLFLT_MATRIX z, PLINT nx, PLINT ny, PLINT opt )
119 {
120  plfplot3dc( x, y, plf2ops_c(), (PLPointer) z, nx, ny, opt | MESH, NULL, 0 );
121 }
122 
123 void
124 plfmesh( PLFLT_VECTOR x, PLFLT_VECTOR y, PLF2OPS zops, PLPointer zp,
125  PLINT nx, PLINT ny, PLINT opt )
126 {
127  plfplot3dc( x, y, zops, zp, nx, ny, opt | MESH, NULL, 0 );
128 }
129 
130 //--------------------------------------------------------------------------
131 // void plmeshc(x, y, z, nx, ny, opt, clevel, nlevel)
132 //
133 // Plots a mesh representation of the function z[x][y]. The x values
134 // are stored as x[0..nx-1], the y values as y[0..ny-1], and the
135 // z values are in the 2-d array z[][]. The integer "opt" specifies:
136 //
137 // DRAW_LINEX draw lines parallel to the X axis
138 // DRAW_LINEY draw lines parallel to the Y axis
139 // DRAW_LINEXY draw lines parallel to both the X and Y axis
140 // MAG_COLOR draw the mesh with a color dependent of the magnitude
141 // BASE_CONT draw contour plot at bottom xy plane
142 // TOP_CONT draw contour plot at top xy plane (not yet)
143 // DRAW_SIDES draw sides
144 //
145 // or any bitwise combination, e.g. "MAG_COLOR | DRAW_LINEX"
146 //
147 //--------------------------------------------------------------------------
148 
149 void
150 c_plmeshc( PLFLT_VECTOR x, PLFLT_VECTOR y, PLFLT_MATRIX z, PLINT nx, PLINT ny, PLINT opt,
151  PLFLT_VECTOR clevel, PLINT nlevel )
152 {
153  plfplot3dc( x, y, plf2ops_c(), (PLPointer) z, nx, ny, opt | MESH, clevel, nlevel );
154 }
155 
156 void
157 plfmeshc( PLFLT_VECTOR x, PLFLT_VECTOR y, PLF2OPS zops, PLPointer zp,
158  PLINT nx, PLINT ny, PLINT opt, PLFLT_VECTOR clevel, PLINT nlevel )
159 {
160  plfplot3dc( x, y, zops, zp, nx, ny, opt | MESH, clevel, nlevel );
161 }
162 
163 // clipping helper for 3d polygons
164 
165 int
166 plP_clip_poly( int Ni, PLFLT *Vi[3], int axis, PLFLT dir, PLFLT offset )
167 {
168  int anyout = 0;
169  PLFLT _in[PL_MAXPOLY], _T[3][PL_MAXPOLY];
170  PLFLT *in, *T[3], *TT = NULL;
171  int No = 0;
172  int i, j, k;
173 
174  if ( Ni > PL_MAXPOLY )
175  {
176  in = (PLFLT *) malloc( sizeof ( PLFLT ) * (size_t) Ni );
177  TT = (PLFLT *) malloc( 3 * sizeof ( PLFLT ) * (size_t) Ni );
178 
179  if ( in == NULL || TT == NULL )
180  {
181  plexit( "plP_clip_poly: insufficient memory for large polygon" );
182  }
183 
184  T[0] = &TT[0];
185  T[1] = &TT[Ni];
186  T[2] = &TT[2 * Ni];
187  }
188  else
189  {
190  in = _in;
191  T[0] = &_T[0][0];
192  T[1] = &_T[1][0];
193  T[2] = &_T[2][0];
194  }
195 
196  for ( i = 0; i < Ni; i++ )
197  {
198  in[i] = Vi[axis][i] * dir + offset;
199  anyout += in[i] < 0;
200  }
201 
202  // none out
203  if ( anyout == 0 )
204  return Ni;
205 
206  // all out
207  if ( anyout == Ni )
208  {
209  return 0;
210  }
211 
212  // some out
213  // copy over to a temp vector
214  for ( i = 0; i < 3; i++ )
215  {
216  for ( j = 0; j < Ni; j++ )
217  {
218  T[i][j] = Vi[i][j];
219  }
220  }
221  // copy back selectively
222  for ( i = 0; i < Ni; i++ )
223  {
224  j = ( i + 1 ) % Ni;
225 
226  if ( in[i] >= 0 && in[j] >= 0 )
227  {
228  for ( k = 0; k < 3; k++ )
229  Vi[k][No] = T[k][j];
230  No++;
231  }
232  else if ( in[i] >= 0 && in[j] < 0 )
233  {
234  PLFLT u = in[i] / ( in[i] - in[j] );
235  for ( k = 0; k < 3; k++ )
236  Vi[k][No] = T[k][i] * ( 1 - u ) + T[k][j] * u;
237  No++;
238  }
239  else if ( in[i] < 0 && in[j] >= 0 )
240  {
241  PLFLT u = in[i] / ( in[i] - in[j] );
242  for ( k = 0; k < 3; k++ )
243  Vi[k][No] = T[k][i] * ( 1 - u ) + T[k][j] * u;
244  No++;
245  for ( k = 0; k < 3; k++ )
246  Vi[k][No] = T[k][j];
247  No++;
248  }
249  }
250 
251  if ( Ni > PL_MAXPOLY )
252  {
253  free( in );
254  free( TT );
255  }
256 
257  return No;
258 }
259 
260 // helper for plsurf3d, similar to c_plfill3()
261 static void
262 shade_triangle( PLFLT x0, PLFLT y0, PLFLT z0,
263  PLFLT x1, PLFLT y1, PLFLT z1,
264  PLFLT x2, PLFLT y2, PLFLT z2 )
265 {
266  int i;
267  // arrays for interface to core functions
268  short u[6], v[6];
269  PLFLT x[6], y[6], z[6];
270  int n;
271  PLFLT xmin, xmax, ymin, ymax, zmin, zmax, zscale;
272  PLFLT *V[3];
273 
274  plP_gdom( &xmin, &xmax, &ymin, &ymax );
275  plP_grange( &zscale, &zmin, &zmax );
276 
277  x[0] = x0; x[1] = x1; x[2] = x2;
278  y[0] = y0; y[1] = y1; y[2] = y2;
279  z[0] = z0; z[1] = z1; z[2] = z2;
280  n = 3;
281 
282  V[0] = x; V[1] = y; V[2] = z;
283 
284  n = plP_clip_poly( n, V, 0, 1, -xmin );
285  n = plP_clip_poly( n, V, 0, -1, xmax );
286  n = plP_clip_poly( n, V, 1, 1, -ymin );
287  n = plP_clip_poly( n, V, 1, -1, ymax );
288  n = plP_clip_poly( n, V, 2, 1, -zmin );
289  n = plP_clip_poly( n, V, 2, -1, zmax );
290 
291  if ( n > 0 )
292  {
293  if ( falsecolor )
294  plcol1( ( ( z[0] + z[1] + z[2] ) / 3. - fc_minz ) / ( fc_maxz - fc_minz ) );
295  else
296  plcol1( plGetAngleToLight( x, y, z ) );
297 
298  for ( i = 0; i < n; i++ )
299  {
300  u[i] = (short) plP_wcpcx( plP_w3wcx( x[i], y[i], z[i] ) );
301  v[i] = (short) plP_wcpcy( plP_w3wcy( x[i], y[i], z[i] ) );
302  }
303  u[n] = u[0];
304  v[n] = v[0];
305 
306 #ifdef SHADE_DEBUG // show triangles
307  plcol0( 1 );
308  x[3] = x[0]; y[3] = y[0]; z[3] = z[0];
309  plline3( 4, x, y, z );
310 #else // fill triangles
311  plP_fill( u, v, n + 1 );
312 #endif
313  }
314 }
315 
316 //--------------------------------------------------------------------------
317 // void plsurf3d(x, y, z, nx, ny, opt, clevel, nlevel)
318 //
319 // Plots the 3-d surface representation of the function z[x][y].
320 // The x values are stored as x[0..nx-1], the y values as y[0..ny-1],
321 // and the z values are in the 2-d array z[][]. The integer "opt" specifies:
322 // see plsurf3dl() below.
323 //--------------------------------------------------------------------------
324 
325 void
326 c_plsurf3d( PLFLT_VECTOR x, PLFLT_VECTOR y, PLFLT_MATRIX z, PLINT nx, PLINT ny,
327  PLINT opt, PLFLT_VECTOR clevel, PLINT nlevel )
328 {
329  plfsurf3d( x, y, plf2ops_c(), (PLPointer) z, nx, ny,
330  opt, clevel, nlevel );
331 }
332 
333 void
334 plfsurf3d( PLFLT_VECTOR x, PLFLT_VECTOR y, PLF2OPS zops, PLPointer zp,
335  PLINT nx, PLINT ny, PLINT opt, PLFLT_VECTOR clevel, PLINT nlevel )
336 {
337  PLINT i;
338  PLINT *indexymin = (PLINT *) malloc( (size_t) nx * sizeof ( PLINT ) );
339  PLINT *indexymax = (PLINT *) malloc( (size_t) nx * sizeof ( PLINT ) );
340 
341  if ( !indexymin || !indexymax )
342  plexit( "plsurf3d: Out of memory." );
343  for ( i = 0; i < nx; i++ )
344  {
345  indexymin[i] = 0;
346  indexymax[i] = ny;
347  }
348  plfsurf3dl( x, y, zops, zp, nx, ny, opt, clevel, nlevel,
349  0, nx, indexymin, indexymax );
350  free_mem( indexymin );
351  free_mem( indexymax );
352 }
353 
354 //--------------------------------------------------------------------------
355 // void plsurf3dl(x, y, z, nx, ny, opt, clevel, nlevel,
356 // indexxmin, indexxmax, indexymin, indexymax)
357 //
358 // Plots the 3-d surface representation of the function z[x][y].
359 // The x values are stored as x[0..nx-1], the y values as y[0..ny-1],
360 // and the z values are in the 2-d array z[][].
361 //
362 //
363 // BASE_CONT draw contour plot at bottom xy plane
364 // TOP_CONT draw contour plot at top xy plane (not implemented)
365 // SURF_CONT draw contour at surface
366 // FACETED each square that makes up the surface is faceted
367 // DRAW_SIDES draw sides
368 // MAG_COLOR the surface is colored according to the value of z;
369 // if not defined, the surface is colored according to the
370 // intensity of the reflected light in the surface from a light
371 // source whose position is set using c_pllightsource()
372 //
373 // or any bitwise combination, e.g. "MAG_COLOR | SURF_CONT | BASE_CONT"
374 //
375 // indexymin and indexymax are arrays which specify the y index range
376 // (following the convention that the upper range limit is one more than
377 // actual index limit) for an x index range of indexxmin, indexxmax.
378 // This code is a complete departure from the approach taken in the old version
379 // of this routine. Formerly to code attempted to use the logic for the hidden
380 // line algorithm to draw the hidden surface. This was really hard. This code
381 // below uses a simple back to front (painters) algorithm. All the
382 // triangles are drawn.
383 //
384 // There are multitude of ways this code could be optimized. Given the
385 // problems with the old code, I tried to focus on clarity here.
386 //--------------------------------------------------------------------------
387 
388 void
389 c_plsurf3dl( PLFLT_VECTOR x, PLFLT_VECTOR y, PLFLT_MATRIX z, PLINT nx, PLINT ny,
390  PLINT opt, PLFLT_VECTOR clevel, PLINT nlevel,
391  PLINT indexxmin, PLINT indexxmax, PLINT_VECTOR indexymin, PLINT_VECTOR indexymax )
392 {
393  plfsurf3dl( x, y, plf2ops_c(), (PLPointer) z, nx, ny,
394  opt, clevel, nlevel, indexxmin, indexxmax, indexymin, indexymax );
395 }
396 
397 void
398 plfsurf3dl( PLFLT_VECTOR x, PLFLT_VECTOR y, PLF2OPS zops, PLPointer zp, PLINT nx, PLINT ny,
399  PLINT opt, PLFLT_VECTOR clevel, PLINT nlevel,
400  PLINT indexxmin, PLINT indexxmax, PLINT_VECTOR indexymin, PLINT_VECTOR indexymax )
401 {
402  PLFLT cxx, cxy, cyx, cyy, cyz;
403  PLINT i, j, k;
404  PLINT ixDir, ixOrigin, iyDir, iyOrigin, nFast, nSlow;
405  PLINT ixFast, ixSlow, iyFast, iySlow;
406  PLINT iFast, iSlow;
407  PLFLT xmin, xmax, ymin, ymax, zmin, zmax, zscale;
408  PLFLT xm, ym, zm;
409  PLINT ixmin = 0, ixmax = nx, iymin = 0, iymax = ny;
410  PLFLT xx[3], yy[3], zz[3];
411  PLFLT px[4], py[4], pz[4];
412  CONT_LEVEL *cont, *clev;
413  CONT_LINE *cline;
414  int ct, ix, iy, iftriangle;
415  PLINT color = plsc->icol0;
416  PLFLT width = plsc->width;
417  PLFLT ( *getz )( PLPointer, PLINT, PLINT ) = zops->get;
418 
419  if ( plsc->level < 3 )
420  {
421  myabort( "plsurf3dl: Please set up window first" );
422  return;
423  }
424 
425  if ( nx <= 0 || ny <= 0 )
426  {
427  myabort( "plsurf3dl: Bad array dimensions." );
428  return;
429  }
430 
431  //
432  // Don't use the data z value to scale the color, use the z axis
433  // values set by plw3d()
434  //
435  // plMinMax2dGrid(z, nx, ny, &fc_maxz, &fc_minz);
436  //
437 
438  fc_minz = plsc->ranmi;
439  fc_maxz = plsc->ranma;
440  if ( fc_maxz == fc_minz )
441  {
442  plwarn( "plsurf3dl: Maximum and minimum Z values are equal! \"fixing\"..." );
443  fc_maxz = fc_minz + 1e-6;
444  }
445 
446  if ( opt & MAG_COLOR )
447  falsecolor = 1;
448  else
449  falsecolor = 0;
450 
451  plP_gdom( &xmin, &xmax, &ymin, &ymax );
452  plP_grange( &zscale, &zmin, &zmax );
453  if ( zmin > zmax )
454  {
455  PLFLT t = zmin;
456  zmin = zmax;
457  zmax = t;
458  }
459 
460  // Check that points in x and in y are strictly increasing and in range
461 
462  for ( i = 0; i < nx - 1; i++ )
463  {
464  if ( x[i] >= x[i + 1] )
465  {
466  myabort( "plsurf3dl: X array must be strictly increasing" );
467  return;
468  }
469  if ( x[i] < xmin && x[i + 1] >= xmin )
470  ixmin = i;
471  if ( x[i + 1] > xmax && x[i] <= xmax )
472  ixmax = i + 2;
473  }
474  for ( i = 0; i < ny - 1; i++ )
475  {
476  if ( y[i] >= y[i + 1] )
477  {
478  myabort( "plsurf3dl: Y array must be strictly increasing" );
479  return;
480  }
481  if ( y[i] < ymin && y[i + 1] >= ymin )
482  iymin = i;
483  if ( y[i + 1] > ymax && y[i] <= ymax )
484  iymax = i + 2;
485  }
486 
487  // get the viewing parameters
488  plP_gw3wc( &cxx, &cxy, &cyx, &cyy, &cyz );
489 
490  // we're going to draw from back to front
491 
492  // iFast will index the dominant (fastest changing) dimension
493  // iSlow will index the slower changing dimension
494  //
495  // iX indexes the X dimension
496  // iY indexes the Y dimension
497 
498  // get direction for X
499  if ( cxy >= 0 )
500  {
501  ixDir = 1; // direction in X
502  ixOrigin = ixmin; // starting point
503  }
504  else
505  {
506  ixDir = -1;
507  ixOrigin = ixmax - 1;
508  }
509  // get direction for Y
510  if ( cxx >= 0 )
511  {
512  iyDir = -1;
513  iyOrigin = iymax - 1;
514  }
515  else
516  {
517  iyDir = 1;
518  iyOrigin = iymin;
519  }
520  // figure out which dimension is dominant
521  if ( fabs( cxx ) > fabs( cxy ) )
522  {
523  // X is dominant
524  nFast = ixmax - ixmin; // samples in the Fast direction
525  nSlow = iymax - iymin; // samples in the Slow direction
526 
527  ixFast = ixDir; ixSlow = 0;
528  iyFast = 0; iySlow = iyDir;
529  }
530  else
531  {
532  nFast = iymax - iymin;
533  nSlow = ixmax - ixmin;
534 
535  ixFast = 0; ixSlow = ixDir;
536  iyFast = iyDir; iySlow = 0;
537  }
538 
539  // we've got to draw the background grid first, hidden line code has to draw it last
540  if ( zbflg )
541  {
542  PLFLT bx[3], by[3], bz[3];
543  PLFLT tick = zbtck, tp;
544  PLINT nsub = 0;
545 
546  // get the tick spacing
547  pldtik( zmin, zmax, &tick, &nsub, FALSE );
548 
549  // determine the vertices for the background grid line
550  bx[0] = ( ixOrigin != ixmin && ixFast == 0 ) || ixFast > 0 ? xmax : xmin;
551  by[0] = ( iyOrigin != iymin && iyFast == 0 ) || iyFast > 0 ? ymax : ymin;
552  bx[1] = ixOrigin != ixmin ? xmax : xmin;
553  by[1] = iyOrigin != iymin ? ymax : ymin;
554  bx[2] = ( ixOrigin != ixmin && ixSlow == 0 ) || ixSlow > 0 ? xmax : xmin;
555  by[2] = ( iyOrigin != iymin && iySlow == 0 ) || iySlow > 0 ? ymax : ymin;
556 
557  plwidth( zbwidth );
558  plcol0( zbcol );
559  for ( tp = tick * floor( zmin / tick ) + tick; tp <= zmax; tp += tick )
560  {
561  bz[0] = bz[1] = bz[2] = tp;
562  plline3( 3, bx, by, bz );
563  }
564  // draw the vertical line at the back corner
565  bx[0] = bx[1];
566  by[0] = by[1];
567  bz[0] = zmin;
568  plline3( 2, bx, by, bz );
569  plwidth( width );
570  plcol0( color );
571  }
572 
573  // If enabled, draw the contour at the base
574 
575  // The contour plotted at the base will be identical to the one obtained
576  // with c_plcont(). The contour plotted at the surface is simple minded, but
577  // can be improved by using the contour data available.
578  //
579 
580  if ( clevel != NULL && opt & BASE_CONT )
581  {
582 #define NPTS 100
583  int np = NPTS;
584  PLFLT **zstore;
585  PLcGrid2 cgrid2;
586  PLFLT *zzloc = (PLFLT *) malloc( (size_t) NPTS * sizeof ( PLFLT ) );
587  if ( zzloc == NULL )
588  plexit( "plsurf3dl: Insufficient memory" );
589 
590  // get the contour lines
591 
592  // prepare cont_store input
593  cgrid2.nx = nx;
594  cgrid2.ny = ny;
595  plAlloc2dGrid( &cgrid2.xg, nx, ny );
596  plAlloc2dGrid( &cgrid2.yg, nx, ny );
597  plAlloc2dGrid( &zstore, nx, ny );
598 
599  for ( i = indexxmin; i < indexxmax; i++ )
600  {
601  for ( j = 0; j < indexymin[i]; j++ )
602  {
603  cgrid2.xg[i][j] = x[i];
604  cgrid2.yg[i][j] = y[indexymin[i]];
605  zstore[i][j] = getz( zp, i, indexymin[i] );
606  }
607  for ( j = indexymin[i]; j < indexymax[i]; j++ )
608  {
609  cgrid2.xg[i][j] = x[i];
610  cgrid2.yg[i][j] = y[j];
611  zstore[i][j] = getz( zp, i, j );
612  }
613  for ( j = indexymax[i]; j < ny; j++ )
614  {
615  cgrid2.xg[i][j] = x[i];
616  cgrid2.yg[i][j] = y[indexymax[i] - 1];
617  zstore[i][j] = getz( zp, i, indexymax[i] - 1 );
618  }
619  }
620  // Fill cont structure with contours.
621  cont_store( (PLFLT_MATRIX) zstore, nx, ny, indexxmin + 1, indexxmax, 1, ny,
622  clevel, nlevel, pltr2, (void *) &cgrid2, &cont );
623 
624  // Free the 2D input arrays to cont_store since not needed any more.
625  plFree2dGrid( zstore, nx, ny );
626  plFree2dGrid( cgrid2.xg, nx, ny );
627  plFree2dGrid( cgrid2.yg, nx, ny );
628 
629  // follow the contour levels and lines
630  clev = cont;
631  do // for each contour level
632  {
633  cline = clev->line;
634  do // there are several lines that make up the contour
635  {
636  if ( cline->npts > np )
637  {
638  np = cline->npts;
639  if ( ( zzloc = (PLFLT *) realloc( zzloc, (size_t) np * sizeof ( PLFLT ) ) ) == NULL )
640  {
641  plexit( "plsurf3dl: Insufficient memory" );
642  }
643  }
644  for ( j = 0; j < cline->npts; j++ )
645  zzloc[j] = plsc->ranmi;
646  if ( cline->npts > 0 )
647  {
648  plcol1( ( clev->level - fc_minz ) / ( fc_maxz - fc_minz ) );
649  plline3( cline->npts, cline->x, cline->y, zzloc );
650  }
651  cline = cline->next;
652  }
653  while ( cline != NULL );
654  clev = clev->next;
655  }
656  while ( clev != NULL );
657 
658  cont_clean_store( cont ); // now release the memory
659  free( zzloc );
660  }
661 
662  // Now we can iterate over the grid drawing the quads
663  for ( iSlow = 0; iSlow < nSlow - 1; iSlow++ )
664  {
665  for ( iFast = 0; iFast < nFast - 1; iFast++ )
666  {
667  // get the 4 corners of the Quad, which are
668  //
669  // 0--2
670  // | |
671  // 1--3
672  //
673 
674  xm = ym = zm = 0.;
675 
676  iftriangle = 1;
677  for ( i = 0; i < 2; i++ )
678  {
679  for ( j = 0; j < 2; j++ )
680  {
681  // we're transforming from Fast/Slow coordinates to x/y coordinates
682  // note, these are the indices, not the values
683  ix = ixFast * ( iFast + i ) + ixSlow * ( iSlow + j ) + ixOrigin;
684  iy = iyFast * ( iFast + i ) + iySlow * ( iSlow + j ) + iyOrigin;
685 
686  if ( indexxmin <= ix && ix < indexxmax &&
687  indexymin[ix] <= iy && iy < indexymax[ix] )
688  {
689  xm += px[2 * i + j] = x[ix];
690  ym += py[2 * i + j] = y[iy];
691  zm += pz[2 * i + j] = getz( zp, ix, iy );
692  }
693  else
694  {
695  iftriangle = 0;
696  break;
697  }
698  }
699  if ( iftriangle == 0 )
700  break;
701  }
702 
703  if ( iftriangle == 0 )
704  continue;
705  // the "mean point" of the quad, common to all four triangles
706  // -- perhaps not a good thing to do for the light shading
707 
708  xm /= 4.; ym /= 4.; zm /= 4.;
709 
710  // now draw the quad as four triangles
711 
712  for ( i = 1; i < 3; i++ )
713  {
714  for ( j = 0; j < 4; j += 3 )
715  {
716  shade_triangle( px[j], py[j], pz[j], xm, ym, zm, px[i], py[i], pz[i] );
717  }
718  }
719 
720  // After shading completed for a quad, render surface contours.
721  if ( clevel != NULL && ( opt & SURF_CONT ) )
722  {
723  for ( i = 1; i < 3; i++ )
724  {
725  for ( j = 0; j < 4; j += 3 )
726 #define min3( a, b, c ) ( MIN( ( MIN( a, b ) ), c ) )
727 #define max3( a, b, c ) ( MAX( ( MAX( a, b ) ), c ) )
728 
729  {
730  for ( k = 0; k < nlevel; k++ )
731  {
732  if ( clevel[k] >= min3( pz[i], zm, pz[j] ) && clevel[k] < max3( pz[i], zm, pz[j] ) )
733  {
734  ct = 0;
735  if ( clevel[k] >= MIN( pz[i], zm ) && clevel[k] < MAX( pz[i], zm ) ) // p0-pm
736  {
737  xx[ct] = ( ( clevel[k] - pz[i] ) * ( xm - px[i] ) ) / ( zm - pz[i] ) + px[i];
738  yy[ct] = ( ( clevel[k] - pz[i] ) * ( ym - py[i] ) ) / ( zm - pz[i] ) + py[i];
739  ct++;
740  }
741 
742  if ( clevel[k] >= MIN( pz[i], pz[j] ) && clevel[k] < MAX( pz[i], pz[j] ) ) // p0-p1
743  {
744  xx[ct] = ( ( clevel[k] - pz[i] ) * ( px[j] - px[i] ) ) / ( pz[j] - pz[i] ) + px[i];
745  yy[ct] = ( ( clevel[k] - pz[i] ) * ( py[j] - py[i] ) ) / ( pz[j] - pz[i] ) + py[i];
746  ct++;
747  }
748 
749  if ( clevel[k] >= MIN( pz[j], zm ) && clevel[k] < MAX( pz[j], zm ) ) // p1-pm
750  {
751  xx[ct] = ( ( clevel[k] - pz[j] ) * ( xm - px[j] ) ) / ( zm - pz[j] ) + px[j];
752  yy[ct] = ( ( clevel[k] - pz[j] ) * ( ym - py[j] ) ) / ( zm - pz[j] ) + py[j];
753  ct++;
754  }
755 
756  if ( ct == 2 )
757  {
758  // yes, xx and yy are the intersection points of the triangle with
759  // the contour line -- draw a straight line betweeen the points
760  // -- at the end this will make up the contour line
761 
762  // surface contour with color set by user
763  plcol0( color );
764  zz[0] = zz[1] = clevel[k];
765  plline3( 2, xx, yy, zz );
766 
767  // don't break; one triangle can span various contour levels
768  }
769  else
770  plwarn( "plsurf3dl: ***ERROR***\n" );
771  }
772  }
773  }
774  }
775  }
776  }
777  }
778 
779  if ( opt & FACETED )
780  {
781  plcol0( 0 );
782  plfplot3dcl( x, y, zops, zp, nx, ny, MESH | DRAW_LINEXY, NULL, 0,
783  indexxmin, indexxmax, indexymin, indexymax );
784  }
785 
786  if ( opt & DRAW_SIDES ) // the sides look ugly !!!
787  { // draw one more row with all the Z's set to zmin
788  plP_grange( &zscale, &zmin, &zmax );
789 
790  iSlow = nSlow - 1;
791  iftriangle = 1;
792  for ( iFast = 0; iFast < nFast - 1; iFast++ )
793  {
794  for ( i = 0; i < 2; i++ )
795  {
796  ix = ixFast * ( iFast + i ) + ixSlow * iSlow + ixOrigin;
797  iy = iyFast * ( iFast + i ) + iySlow * iSlow + iyOrigin;
798  if ( indexxmin <= ix && ix < indexxmax &&
799  indexymin[ix] <= iy && iy < indexymax[ix] )
800  {
801  px[2 * i] = x[ix];
802  py[2 * i] = y[iy];
803  pz[2 * i] = getz( zp, ix, iy );
804  }
805  else
806  {
807  iftriangle = 0;
808  break;
809  }
810  }
811  if ( iftriangle == 0 )
812  break;
813  // now draw the quad as two triangles (4 might be better)
814 
815  shade_triangle( px[0], py[0], pz[0], px[2], py[2], pz[2], px[0], py[0], zmin );
816  shade_triangle( px[2], py[2], pz[2], px[2], py[2], zmin, px[0], py[0], zmin );
817  }
818 
819  iFast = nFast - 1;
820  iftriangle = 1;
821  for ( iSlow = 0; iSlow < nSlow - 1; iSlow++ )
822  {
823  for ( i = 0; i < 2; i++ )
824  {
825  ix = ixFast * iFast + ixSlow * ( iSlow + i ) + ixOrigin;
826  iy = iyFast * iFast + iySlow * ( iSlow + i ) + iyOrigin;
827  if ( indexxmin <= ix && ix < indexxmax &&
828  indexymin[ix] <= iy && iy < indexymax[ix] )
829  {
830  px[2 * i] = x[ix];
831  py[2 * i] = y[iy];
832  pz[2 * i] = getz( zp, ix, iy );
833  }
834  else
835  {
836  iftriangle = 0;
837  break;
838  }
839  }
840  if ( iftriangle == 0 )
841  break;
842 
843  // now draw the quad as two triangles (4 might be better)
844  shade_triangle( px[0], py[0], pz[0], px[2], py[2], pz[2], px[0], py[0], zmin );
845  shade_triangle( px[2], py[2], pz[2], px[2], py[2], zmin, px[0], py[0], zmin );
846  }
847  }
848 }
849 
850 //--------------------------------------------------------------------------
851 // void plot3d(x, y, z, nx, ny, opt, side)
852 //
853 // Plots a 3-d representation of the function z[x][y]. The x values
854 // are stored as x[0..nx-1], the y values as y[0..ny-1], and the z
855 // values are in the 2-d array z[][]. The integer "opt" specifies:
856 // see plot3dcl() below
857 //--------------------------------------------------------------------------
858 
859 void
860 c_plot3d( PLFLT_VECTOR x, PLFLT_VECTOR y, PLFLT_MATRIX z,
861  PLINT nx, PLINT ny, PLINT opt, PLBOOL side )
862 {
863  plfplot3dc( x, y, plf2ops_c(), (PLPointer) z, nx, ny, opt | ( side != 0 ? DRAW_SIDES : 0 ), NULL, 0 );
864 }
865 
866 void
867 plfplot3d( PLFLT_VECTOR x, PLFLT_VECTOR y, PLF2OPS zops, PLPointer zp,
868  PLINT nx, PLINT ny, PLINT opt, PLBOOL side )
869 {
870  plfplot3dc( x, y, zops, zp, nx, ny, opt | ( side != 0 ? DRAW_SIDES : 0 ), NULL, 0 );
871 }
872 
873 //--------------------------------------------------------------------------
874 // void plot3dc(x, y, z, nx, ny, opt, clevel, nlevel)
875 //
876 // Plots a 3-d representation of the function z[x][y]. The x values
877 // are stored as x[0..nx-1], the y values as y[0..ny-1], and the z
878 // values are in the 2-d array z[][]. The integer "opt" specifies:
879 // see plot3dcl() below
880 //--------------------------------------------------------------------------
881 
882 void
883 c_plot3dc( PLFLT_VECTOR x, PLFLT_VECTOR y, PLFLT_MATRIX z,
884  PLINT nx, PLINT ny, PLINT opt,
885  PLFLT_VECTOR clevel, PLINT nlevel )
886 {
887  plfplot3dcl( x, y, plf2ops_c(), (PLPointer) z, nx, ny, opt, clevel, nlevel, 0, 0, NULL, NULL );
888 }
889 
890 void
891 plfplot3dc( PLFLT_VECTOR x, PLFLT_VECTOR y, PLF2OPS zops, PLPointer zp,
892  PLINT nx, PLINT ny, PLINT opt, PLFLT_VECTOR clevel, PLINT nlevel )
893 {
894  plfplot3dcl( x, y, zops, zp, nx, ny, opt, clevel, nlevel, 0, 0, NULL, NULL );
895 }
896 
897 //--------------------------------------------------------------------------
898 // void plot3dcl(x, y, z, nx, ny, opt, clevel, nlevel,
899 // indexxmin, indexxmax, indexymin, indexymax)
900 //
901 // Plots a 3-d representation of the function z[x][y]. The x values
902 // are stored as x[0..nx-1], the y values as y[0..ny-1], and the z
903 // values are in the 2-d array z[][]. The integer "opt" specifies:
904 //
905 // DRAW_LINEX : Draw lines parallel to x-axis
906 // DRAW_LINEY : Draw lines parallel to y-axis
907 // DRAW_LINEXY: Draw lines parallel to both axes
908 // MAG_COLOR: Magnitude coloring of wire frame
909 // BASE_CONT: Draw contour at bottom xy plane
910 // TOP_CONT: Draw contour at top xy plane (not yet)
911 // DRAW_SIDES: Draw sides around the plot
912 // MESH: Draw the "under" side of the plot
913 //
914 // or any bitwise combination, e.g. "MAG_COLOR | DRAW_LINEX"
915 // indexymin and indexymax are arrays which specify the y index limits
916 // (following the convention that the upper range limit is one more than
917 // actual index limit) for an x index range of indexxmin, indexxmax.
918 //--------------------------------------------------------------------------
919 
920 void
921 c_plot3dcl( PLFLT_VECTOR x, PLFLT_VECTOR y, PLFLT_MATRIX z,
922  PLINT nx, PLINT ny, PLINT opt,
923  PLFLT_VECTOR clevel, PLINT nlevel,
924  PLINT indexxmin, PLINT indexxmax, PLINT_VECTOR indexymin, PLINT_VECTOR indexymax )
925 {
926  plfplot3dcl( x, y, plf2ops_c(), (PLPointer) z, nx, ny,
927  opt, clevel, nlevel, indexxmin, indexxmax, indexymin, indexymax );
928 }
929 
930 //--------------------------------------------------------------------------
967 //--------------------------------------------------------------------------
968 
969 void
970 plfplot3dcl( PLFLT_VECTOR x, PLFLT_VECTOR y, PLF2OPS zops, PLPointer zp,
971  PLINT nx, PLINT ny, PLINT opt,
972  PLFLT_VECTOR clevel, PLINT nlevel,
973  PLINT PL_UNUSED( indexxmin ), PLINT PL_UNUSED( indexxmax ), PLINT_VECTOR PL_UNUSED( indexymin ), PLINT_VECTOR PL_UNUSED( indexymax ) )
974 {
975  PLFLT cxx, cxy, cyx, cyy, cyz;
976  PLINT init, ix, iy, color;
977  PLFLT width;
978  PLFLT xmin, xmax, ymin, ymax, zmin, zmax, zscale;
979  PLINT ixmin = 0, ixmax = nx - 1, iymin = 0, iymax = ny - 1;
980  PLINT clipped = 0, base_cont = 0, side = 0;
981  PLFLT ( *getz )( PLPointer, PLINT, PLINT ) = zops->get;
982  PLFLT *_x = NULL, *_y = NULL, **_z = NULL;
983  PLFLT_VECTOR x_modified, y_modified;
984  int i;
985 
986  pl3mode = 0;
987 
988  if ( plsc->level < 3 )
989  {
990  myabort( "plot3dcl: Please set up window first" );
991  return;
992  }
993 
994  if ( ( opt & 3 ) == 0 )
995  {
996  myabort( "plot3dcl: Bad option" );
997  return;
998  }
999 
1000  if ( nx <= 0 || ny <= 0 )
1001  {
1002  myabort( "plot3dcl: Bad array dimensions." );
1003  return;
1004  }
1005 
1006  plP_gdom( &xmin, &xmax, &ymin, &ymax );
1007  plP_grange( &zscale, &zmin, &zmax );
1008  if ( zmin > zmax )
1009  {
1010  PLFLT t = zmin;
1011  zmin = zmax;
1012  zmax = t;
1013  }
1014 
1015 // Check that points in x and in y are strictly increasing
1016 
1017  for ( i = 0; i < nx - 1; i++ )
1018  {
1019  if ( x[i] >= x[i + 1] )
1020  {
1021  myabort( "plot3dcl: X array must be strictly increasing" );
1022  return;
1023  }
1024  }
1025  for ( i = 0; i < ny - 1; i++ )
1026  {
1027  if ( y[i] >= y[i + 1] )
1028  {
1029  myabort( "plot3dcl: Y array must be strictly increasing" );
1030  return;
1031  }
1032  }
1033 
1034  if ( opt & MESH )
1035  pl3mode = 1;
1036 
1037  if ( opt & DRAW_SIDES )
1038  side = 1;
1039 
1040  // figure out the part of the data to use
1041  if ( xmin < x[0] )
1042  xmin = x[0];
1043  if ( xmax > x[nx - 1] )
1044  xmax = x[nx - 1];
1045  if ( ymin < y[0] )
1046  ymin = y[0];
1047  if ( ymax > y[ny - 1] )
1048  ymax = y[ny - 1];
1049  for ( ixmin = 0; ixmin < nx - 1 && x[ixmin + 1] < xmin; ixmin++ )
1050  {
1051  }
1052  for ( ixmax = nx - 1; ixmax > 0 && x[ixmax - 1] > xmax; ixmax-- )
1053  {
1054  }
1055  for ( iymin = 0; iymin < ny - 1 && y[iymin + 1] < ymin; iymin++ )
1056  {
1057  }
1058  for ( iymax = ny - 1; iymax > 0 && y[iymax - 1] > ymax; iymax-- )
1059  {
1060  }
1061  //fprintf(stderr, "(%d,%d) %d %d %d %d\n", nx, ny, ixmin, ixmax, iymin, iymax);
1062  // do we need to clip?
1063  if ( ixmin > 0 || ixmax < nx - 1 || iymin > 0 || iymax < ny - 1 )
1064  {
1065  // adjust the input so it stays within bounds
1066  int _nx = ixmax - ixmin + 1;
1067  int _ny = iymax - iymin + 1;
1068  PLFLT ty0, ty1, tx0, tx1;
1069  int j;
1070 
1071  if ( _nx <= 1 || _ny <= 1 )
1072  {
1073  myabort( "plot3dcl: selected x or y range has no data" );
1074  return;
1075  }
1076 
1077  // allocate storage for new versions of the input vectors
1078  if ( ( ( _x = (PLFLT *) malloc( (size_t) _nx * sizeof ( PLFLT ) ) ) == NULL ) ||
1079  ( ( _y = (PLFLT *) malloc( (size_t) _ny * sizeof ( PLFLT ) ) ) == NULL ) ||
1080  ( ( _z = (PLFLT **) malloc( (size_t) _nx * sizeof ( PLFLT* ) ) ) == NULL ) )
1081  {
1082  plexit( "c_plot3dcl: Insufficient memory" );
1083  }
1084 
1085  clipped = 1;
1086 
1087  // copy over the independent variables
1088  _x[0] = xmin;
1089  _x[_nx - 1] = xmax;
1090  for ( i = 1; i < _nx - 1; i++ )
1091  _x[i] = x[ixmin + i];
1092  _y[0] = ymin;
1093  _y[_ny - 1] = ymax;
1094  for ( i = 1; i < _ny - 1; i++ )
1095  _y[i] = y[iymin + i];
1096 
1097  // copy the data array so we can interpolate around the edges
1098  for ( i = 0; i < _nx; i++ )
1099  {
1100  if ( ( _z[i] = (PLFLT *) malloc( (size_t) _ny * sizeof ( PLFLT ) ) ) == NULL )
1101  {
1102  plexit( "c_plot3dcl: Insufficient memory" );
1103  }
1104  }
1105 
1106  // interpolation factors for the 4 edges
1107  ty0 = ( _y[0] - y[iymin] ) / ( y[iymin + 1] - y[iymin] );
1108  ty1 = ( _y[_ny - 1] - y[iymax - 1] ) / ( y[iymax] - y[iymax - 1] );
1109  tx0 = ( _x[0] - x[ixmin] ) / ( x[ixmin + 1] - x[ixmin] );
1110  tx1 = ( _x[_nx - 1] - x[ixmax - 1] ) / ( x[ixmax] - x[ixmax - 1] );
1111  for ( i = 0; i < _nx; i++ )
1112  {
1113  if ( i == 0 )
1114  {
1115  _z[i][0] = getz( zp, ixmin, iymin ) * ( 1 - ty0 ) * ( 1 - tx0 ) + getz( zp, ixmin, iymin + 1 ) * ( 1 - tx0 ) * ty0
1116  + getz( zp, ixmin + 1, iymin ) * tx0 * ( 1 - ty0 ) + getz( zp, ixmin + 1, iymin + 1 ) * tx0 * ty0;
1117  for ( j = 1; j < _ny - 1; j++ )
1118  _z[i][j] = getz( zp, ixmin, iymin + j ) * ( 1 - tx0 ) + getz( zp, ixmin + 1, iymin + j ) * tx0;
1119  _z[i][_ny - 1] = getz( zp, ixmin, iymax - 1 ) * ( 1 - tx0 ) * ( 1 - ty1 ) + getz( zp, ixmin, iymax ) * ( 1 - tx0 ) * ty1
1120  + getz( zp, ixmin + 1, iymax - 1 ) * tx0 * ( 1 - ty1 ) + getz( zp, ixmin + 1, iymax ) * tx0 * ty1;
1121  }
1122  else if ( i == _nx - 1 )
1123  {
1124  _z[i][0] = getz( zp, ixmax - 1, iymin ) * ( 1 - tx1 ) * ( 1 - ty0 ) + getz( zp, ixmax - 1, iymin + 1 ) * ( 1 - tx1 ) * ty0
1125  + getz( zp, ixmax, iymin ) * tx1 * ( 1 - ty0 ) + getz( zp, ixmax, iymin + 1 ) * tx1 * ty0;
1126  for ( j = 1; j < _ny - 1; j++ )
1127  _z[i][j] = getz( zp, ixmax - 1, iymin + j ) * ( 1 - tx1 ) + getz( zp, ixmax, iymin + j ) * tx1;
1128  _z[i][_ny - 1] = getz( zp, ixmax - 1, iymax - 1 ) * ( 1 - tx1 ) * ( 1 - ty1 ) + getz( zp, ixmax, iymax ) * ( 1 - tx1 ) * ty1
1129  + getz( zp, ixmax, iymax - 1 ) * tx1 * ( 1 - ty1 ) + getz( zp, ixmax, iymax ) * tx1 * ty1;
1130  }
1131  else
1132  {
1133  _z[i][0] = getz( zp, ixmin + i, iymin ) * ( 1 - ty0 ) + getz( zp, ixmin + i, iymin + 1 ) * ty0;
1134  for ( j = 1; j < _ny - 1; j++ )
1135  _z[i][j] = getz( zp, ixmin + i, iymin + j );
1136  _z[i][_ny - 1] = getz( zp, ixmin + i, iymax - 1 ) * ( 1 - ty1 ) + getz( zp, ixmin + i, iymax ) * ty1;
1137  }
1138  for ( j = 0; j < _ny; j++ )
1139  {
1140  if ( _z[i][j] < zmin )
1141  _z[i][j] = zmin;
1142  else if ( _z[i][j] > zmax )
1143  _z[i][j] = zmax;
1144  }
1145  }
1146  // replace the input with our clipped versions
1147  zp = (PLPointer) _z;
1148  getz = plf2ops_c()->get;
1149  nx = _nx;
1150  ny = _ny;
1151  // Do not want to modify input x and y (const modifier)
1152  x_modified = (PLFLT_VECTOR) _x;
1153  y_modified = (PLFLT_VECTOR) _y;
1154  }
1155  else
1156  {
1157  x_modified = x;
1158  y_modified = y;
1159  }
1160 
1161  // From here on must use x_modified and y_modified rather than
1162  // x and y.
1163  if ( ( opt & BASE_CONT ) || ( opt & TOP_CONT ) || ( opt & MAG_COLOR ) )
1164  {
1165  //
1166  // Don't use the data z value to scale the color, use the z axis
1167  // values set by plw3d()
1168  //
1169  // plMinMax2dGrid(z, nx, ny, &fc_maxz, &fc_minz);
1170  //
1171 
1172  fc_minz = plsc->ranmi;
1173  fc_maxz = plsc->ranma;
1174 
1175  if ( fc_maxz == fc_minz )
1176  {
1177  plwarn( "plot3dcl: Maximum and minimum Z values are equal! \"fixing\"..." );
1178  fc_maxz = fc_minz + 1e-6;
1179  }
1180  }
1181 
1182  if ( opt & BASE_CONT ) // If enabled, draw the contour at the base.
1183  {
1184  if ( clevel != NULL && nlevel != 0 )
1185  {
1186  base_cont = 1;
1187  // even if MESH is not set, "set it",
1188  // as the base contour can only be done in this case
1189  pl3mode = 1;
1190  }
1191  }
1192 
1193  if ( opt & MAG_COLOR ) // If enabled, use magnitude colored wireframe
1194  {
1195  if ( ( ctmp = (PLFLT *) malloc( (size_t) ( 2 * MAX( nx, ny ) ) * sizeof ( PLFLT ) ) ) == NULL )
1196  {
1197  plexit( "c_plot3dcl: Insufficient memory" );
1198  }
1199  }
1200  else
1201  ctmp = NULL;
1202 
1203  // next logic only knows opt = 1 | 2 | 3, make sure that it only gets that
1204  opt &= DRAW_LINEXY;
1205 
1206  // Allocate work arrays
1207 
1208  utmp = (PLINT *) malloc( (size_t) ( 2 * MAX( nx, ny ) ) * sizeof ( PLINT ) );
1209  vtmp = (PLINT *) malloc( (size_t) ( 2 * MAX( nx, ny ) ) * sizeof ( PLINT ) );
1210 
1211  if ( !utmp || !vtmp )
1212  myexit( "plot3dcl: Out of memory." );
1213 
1214  plP_gw3wc( &cxx, &cxy, &cyx, &cyy, &cyz );
1215  init = 1;
1216 // Call 3d line plotter. Each viewing quadrant
1217 // (perpendicular to x-y plane) must be handled separately.
1218  if ( cxx >= 0.0 && cxy <= 0.0 )
1219  {
1220  if ( opt == DRAW_LINEY )
1221  plt3zz( 1, ny, 1, -1, -opt, &init, x_modified, y_modified, zops, zp, nx, ny, utmp, vtmp, ctmp );
1222  else
1223  {
1224  for ( iy = 2; iy <= ny; iy++ )
1225  plt3zz( 1, iy, 1, -1, -opt, &init, x_modified, y_modified, zops, zp, nx, ny, utmp, vtmp, ctmp );
1226  }
1227  if ( opt == DRAW_LINEX )
1228  plt3zz( 1, ny, 1, -1, opt, &init, x_modified, y_modified, zops, zp, nx, ny, utmp, vtmp, ctmp );
1229  else
1230  {
1231  for ( ix = 1; ix <= nx - 1; ix++ )
1232  plt3zz( ix, ny, 1, -1, opt, &init, x_modified, y_modified, zops, zp, nx, ny, utmp, vtmp, ctmp );
1233  }
1234  }
1235 
1236  else if ( cxx <= 0.0 && cxy <= 0.0 )
1237  {
1238  if ( opt == DRAW_LINEX )
1239  plt3zz( nx, ny, -1, -1, opt, &init, x_modified, y_modified, zops, zp, nx, ny, utmp, vtmp, ctmp );
1240  else
1241  {
1242  for ( ix = 2; ix <= nx; ix++ )
1243  plt3zz( ix, ny, -1, -1, opt, &init, x_modified, y_modified, zops, zp, nx, ny, utmp, vtmp, ctmp );
1244  }
1245  if ( opt == DRAW_LINEY )
1246  plt3zz( nx, ny, -1, -1, -opt, &init, x_modified, y_modified, zops, zp, nx, ny, utmp, vtmp, ctmp );
1247  else
1248  {
1249  for ( iy = ny; iy >= 2; iy-- )
1250  plt3zz( nx, iy, -1, -1, -opt, &init, x_modified, y_modified, zops, zp, nx, ny, utmp, vtmp, ctmp );
1251  }
1252  }
1253 
1254  else if ( cxx <= 0.0 && cxy >= 0.0 )
1255  {
1256  if ( opt == DRAW_LINEY )
1257  plt3zz( nx, 1, -1, 1, -opt, &init, x_modified, y_modified, zops, zp, nx, ny, utmp, vtmp, ctmp );
1258  else
1259  {
1260  for ( iy = ny - 1; iy >= 1; iy-- )
1261  plt3zz( nx, iy, -1, 1, -opt, &init, x_modified, y_modified, zops, zp, nx, ny, utmp, vtmp, ctmp );
1262  }
1263  if ( opt == DRAW_LINEX )
1264  plt3zz( nx, 1, -1, 1, opt, &init, x_modified, y_modified, zops, zp, nx, ny, utmp, vtmp, ctmp );
1265  else
1266  {
1267  for ( ix = nx; ix >= 2; ix-- )
1268  plt3zz( ix, 1, -1, 1, opt, &init, x_modified, y_modified, zops, zp, nx, ny, utmp, vtmp, ctmp );
1269  }
1270  }
1271 
1272  else if ( cxx >= 0.0 && cxy >= 0.0 )
1273  {
1274  if ( opt == DRAW_LINEX )
1275  plt3zz( 1, 1, 1, 1, opt, &init, x_modified, y_modified, zops, zp, nx, ny, utmp, vtmp, ctmp );
1276  else
1277  {
1278  for ( ix = nx - 1; ix >= 1; ix-- )
1279  plt3zz( ix, 1, 1, 1, opt, &init, x_modified, y_modified, zops, zp, nx, ny, utmp, vtmp, ctmp );
1280  }
1281  if ( opt == DRAW_LINEY )
1282  plt3zz( 1, 1, 1, 1, -opt, &init, x_modified, y_modified, zops, zp, nx, ny, utmp, vtmp, ctmp );
1283  else
1284  {
1285  for ( iy = 1; iy <= ny - 1; iy++ )
1286  plt3zz( 1, iy, 1, 1, -opt, &init, x_modified, y_modified, zops, zp, nx, ny, utmp, vtmp, ctmp );
1287  }
1288  }
1289 
1290  // draw contour at the base. Not 100%! Why?
1291 
1292  if ( base_cont )
1293  {
1294  int np = NPTS, j;
1295  CONT_LEVEL *cont, *clev;
1296  CONT_LINE *cline;
1297 
1298  PLINT *uu = (PLINT *) malloc( (size_t) NPTS * sizeof ( PLINT ) );
1299  PLINT *vv = (PLINT *) malloc( (size_t) NPTS * sizeof ( PLINT ) );
1300  // prepare cont_store input
1301  PLFLT **zstore;
1302  PLcGrid2 cgrid2;
1303 
1304  if ( ( uu == NULL ) || ( vv == NULL ) )
1305  {
1306  plexit( "c_plot3dcl: Insufficient memory" );
1307  }
1308 
1309  cgrid2.nx = nx;
1310  cgrid2.ny = ny;
1311  plAlloc2dGrid( &cgrid2.xg, nx, ny );
1312  plAlloc2dGrid( &cgrid2.yg, nx, ny );
1313  plAlloc2dGrid( &zstore, nx, ny );
1314 
1315  for ( i = 0; i < nx; i++ )
1316  {
1317  for ( j = 0; j < ny; j++ )
1318  {
1319  cgrid2.xg[i][j] = x_modified[i];
1320  cgrid2.yg[i][j] = y_modified[j];
1321  zstore[i][j] = getz( zp, i, j );
1322  }
1323  }
1324 
1325  pl3upv = 0;
1326 
1327  // Fill cont structure with contours.
1328  cont_store( (PLFLT_MATRIX) zstore, nx, ny, 1, nx, 1, ny,
1329  clevel, nlevel, pltr2, (void *) &cgrid2, &cont );
1330 
1331  // Free the 2D input arrays to cont_store since not needed any more.
1332  plFree2dGrid( zstore, nx, ny );
1333  plFree2dGrid( cgrid2.xg, nx, ny );
1334  plFree2dGrid( cgrid2.yg, nx, ny );
1335 
1336  // follow the contour levels and lines
1337  clev = cont;
1338  do // for each contour level
1339  {
1340  cline = clev->line;
1341  do // there are several lines that make up each contour
1342  {
1343  int cx, k, l, m, start, end;
1344  PLFLT tx, ty;
1345  if ( cline->npts > np )
1346  {
1347  np = cline->npts;
1348  if ( ( ( uu = (PLINT *) realloc( uu, (size_t) np * sizeof ( PLINT ) ) ) == NULL ) ||
1349  ( ( vv = (PLINT *) realloc( vv, (size_t) np * sizeof ( PLINT ) ) ) == NULL ) )
1350  {
1351  plexit( "c_plot3dcl: Insufficient memory" );
1352  }
1353  }
1354 
1355  // the hidden line plotter plnxtv() only works OK if the x points are in increasing order.
1356  // As this does not always happens, the situation must be detected and the line segment
1357  // must be reversed before being plotted
1358  i = 0;
1359  if ( cline->npts > 0 )
1360  {
1361  do
1362  {
1363  plcol1( ( clev->level - fc_minz ) / ( fc_maxz - fc_minz ) );
1364  cx = plP_wcpcx( plP_w3wcx( cline->x[i], cline->y[i], plsc->ranmi ) );
1365  for ( j = i; j < cline->npts; j++ ) // convert to 2D coordinates
1366  {
1367  uu[j] = plP_wcpcx( plP_w3wcx( cline->x[j], cline->y[j], plsc->ranmi ) );
1368  vv[j] = plP_wcpcy( plP_w3wcy( cline->x[j], cline->y[j], plsc->ranmi ) );
1369  if ( uu[j] < cx ) // find turn back point
1370  break;
1371  else
1372  cx = uu[j];
1373  }
1374  plnxtv( &uu[i], &vv[i], NULL, j - i, 0 ); // plot line with increasing x
1375 
1376  if ( j < cline->npts ) // line not yet finished,
1377  {
1378  start = j - 1;
1379  for ( i = start; i < cline->npts; i++ ) // search turn forward point
1380  {
1381  uu[i] = plP_wcpcx( plP_w3wcx( cline->x[i], cline->y[i], plsc->ranmi ) );
1382  if ( uu[i] > cx )
1383  break;
1384  else
1385  cx = uu[i];
1386  }
1387  end = i - 1;
1388 
1389  for ( k = 0; k < ( end - start + 1 ) / 2; k++ ) // reverse line segment
1390  {
1391  l = start + k;
1392  m = end - k;
1393  tx = cline->x[l];
1394  ty = cline->y[l];
1395  cline->x[l] = cline->x[m];
1396  cline->y[l] = cline->y[m];
1397  cline->x[m] = tx;
1398  cline->y[m] = ty;
1399  }
1400 
1401  // convert to 2D coordinates
1402  for ( j = start; j <= end; j++ )
1403  {
1404  uu[j] = plP_wcpcx( plP_w3wcx( cline->x[j], cline->y[j], plsc->ranmi ) );
1405  vv[j] = plP_wcpcy( plP_w3wcy( cline->x[j], cline->y[j], plsc->ranmi ) );
1406  }
1407  plnxtv( &uu[start], &vv[start], NULL, end - start + 1, 0 ); // and plot it
1408 
1409  cline->x[end] = cline->x[start];
1410  cline->y[end] = cline->y[start];
1411  i = end; // restart where it was left
1412  }
1413  } while ( j < cline->npts && i < cline->npts );
1414  }
1415  cline = cline->next;
1416  }
1417  while ( cline != NULL );
1418  clev = clev->next;
1419  }
1420  while ( clev != NULL );
1421 
1422  cont_clean_store( cont ); // now release the contour memory
1423  pl3upv = 1;
1424  free( uu );
1425  free( vv );
1426  }
1427 
1428 // Finish up by drawing sides, background grid (both are optional)
1429 
1430  if ( side )
1431  plside3( x_modified, y_modified, zops, zp, nx, ny, opt );
1432 
1433  if ( zbflg )
1434  {
1435  color = plsc->icol0;
1436  width = plsc->width;
1437  plwidth( zbwidth );
1438  plcol0( zbcol );
1439  plgrid3( zbtck );
1440  plwidth( width );
1441  plcol0( color );
1442  }
1443 
1444  freework();
1445 
1446  if ( clipped )
1447  {
1448  free( _x );
1449  free( _y );
1450  for ( i = 0; i < nx; i++ )
1451  free( _z[i] );
1452  free( _z );
1453  }
1454 }
1455 
1456 //--------------------------------------------------------------------------
1457 // void plxyindexlimits()
1458 //
1459 // Transform from y array limits to corresponding x array limits (or vice
1460 // versa).
1461 //
1462 // N.B. we follow the convention here that all upper range limits are one
1463 // more than the actual last index.
1464 // instart (>= 0) through inn is the index range where
1465 // the input inarray_min and inarray_max arrays are defined.
1466 //
1467 // outstart (>= 0), through outn (with outn <= outnmax) is the index range
1468 // where the output outarray_min and outarray_max arrays are defined.
1469 //
1470 // In order to assure the transformation from y array limits to x array limits
1471 // (or vice versa) is single-valued, this programme plaborts if the
1472 // inarray_min array has a maximum or inarray_max array has a minimum.
1473 //--------------------------------------------------------------------------
1474 
1475 //static void
1476 //plxyindexlimits( PLINT instart, PLINT inn,
1477 // PLINT *inarray_min, PLINT *inarray_max,
1478 // PLINT *outstart, PLINT *outn, PLINT outnmax,
1479 // PLINT *outarray_min, PLINT *outarray_max )
1480 //{
1481 // PLINT i, j;
1482 // if ( inn < 0 )
1483 // {
1484 // myabort( "plxyindexlimits: Must have instart >= 0" );
1485 // return;
1486 // }
1487 // if ( inn - instart <= 0 )
1488 // {
1489 // myabort( "plxyindexlimits: Must have at least 1 defined point" );
1490 // return;
1491 // }
1492 // *outstart = inarray_min[instart];
1493 // *outn = inarray_max[instart];
1494 // for ( i = instart; i < inn; i++ )
1495 // {
1496 // *outstart = MIN( *outstart, inarray_min[i] );
1497 // *outn = MAX( *outn, inarray_max[i] );
1498 // if ( i + 2 < inn )
1499 // {
1500 // if ( inarray_min[i] < inarray_min[i + 1] &&
1501 // inarray_min[i + 1] > inarray_min[i + 2] )
1502 // {
1503 // myabort( "plxyindexlimits: inarray_min must not have a maximum" );
1504 // return;
1505 // }
1506 // if ( inarray_max[i] > inarray_max[i + 1] &&
1507 // inarray_max[i + 1] < inarray_max[i + 2] )
1508 // {
1509 // myabort( "plxyindexlimits: inarray_max must not have a minimum" );
1510 // return;
1511 // }
1512 // }
1513 // }
1514 // if ( *outstart < 0 )
1515 // {
1516 // myabort( "plxyindexlimits: Must have all elements of inarray_min >= 0" );
1517 // return;
1518 // }
1519 // if ( *outn > outnmax )
1520 // {
1521 // myabort( "plxyindexlimits: Must have all elements of inarray_max <= outnmax" );
1522 // return;
1523 // }
1524 // for ( j = *outstart; j < *outn; j++ )
1525 // {
1526 // i = instart;
1527 // // Find first valid i for this j.
1528 // while ( i < inn && !( inarray_min[i] <= j && j < inarray_max[i] ) )
1529 // i++;
1530 // if ( i < inn )
1531 // outarray_min[j] = i;
1532 // else
1533 // {
1534 // myabort( "plxyindexlimits: bad logic; invalid i should never happen" );
1535 // return;
1536 // }
1537 // // Find next invalid i for this j.
1538 // while ( i < inn && ( inarray_min[i] <= j && j < inarray_max[i] ) )
1539 // i++;
1540 // outarray_max[j] = i;
1541 // }
1542 //}
1543 
1544 //--------------------------------------------------------------------------
1545 // void plP_gzback()
1546 //
1547 // Get background parameters for 3d plot.
1548 //--------------------------------------------------------------------------
1549 
1550 void
1551 plP_gzback( PLINT **zbf, PLINT **zbc, PLFLT **zbt, PLFLT **zbw )
1552 {
1553  *zbf = &zbflg;
1554  *zbc = &zbcol;
1555  *zbt = &zbtck;
1556  *zbw = &zbwidth;
1557 }
1558 
1559 //--------------------------------------------------------------------------
1560 // PLFLT plGetAngleToLight()
1561 //
1562 // Gets cos of angle between normal to a polygon and a light source.
1563 // Requires at least 3 elements, forming non-parallel lines
1564 // in the arrays.
1565 //--------------------------------------------------------------------------
1566 
1567 static PLFLT
1568 plGetAngleToLight( PLFLT* x, PLFLT* y, PLFLT* z )
1569 {
1570  PLFLT vx1, vx2, vy1, vy2, vz1, vz2;
1571  PLFLT px, py, pz;
1572  PLFLT vlx, vly, vlz;
1573  PLFLT mag1, mag2;
1574  PLFLT cosangle;
1575 
1576  vx1 = x[1] - x[0];
1577  vx2 = x[2] - x[1];
1578  vy1 = y[1] - y[0];
1579  vy2 = y[2] - y[1];
1580  vz1 = z[1] - z[0];
1581  vz2 = z[2] - z[1];
1582 
1583 // Find vector perpendicular to the face
1584  px = vy1 * vz2 - vz1 * vy2;
1585  py = vz1 * vx2 - vx1 * vz2;
1586  pz = vx1 * vy2 - vy1 * vx2;
1587  mag1 = px * px + py * py + pz * pz;
1588 
1589 // Vectors were parallel!
1590  if ( mag1 == 0 )
1591  return 1;
1592 
1593  vlx = xlight - x[0];
1594  vly = ylight - y[0];
1595  vlz = zlight - z[0];
1596  mag2 = vlx * vlx + vly * vly + vlz * vlz;
1597  if ( mag2 == 0 )
1598  return 1;
1599 
1600 // Now have 3 vectors going through the first point on the given surface
1601  cosangle = fabs( ( vlx * px + vly * py + vlz * pz ) / sqrt( mag1 * mag2 ) );
1602 
1603 // In case of numerical rounding
1604  if ( cosangle > 1 )
1605  cosangle = 1;
1606  return cosangle;
1607 }
1608 
1609 //--------------------------------------------------------------------------
1610 // void plt3zz()
1611 //
1612 // Draws the next zig-zag line for a 3-d plot. The data is stored in array
1613 // z[][] as a function of x[] and y[], and is plotted out starting at index
1614 // (x0,y0).
1615 //
1616 // Depending on the state of "flag", the sequence of data points sent to
1617 // plnxtv is altered so as to allow cross-hatch plotting, or plotting
1618 // parallel to either the x-axis or the y-axis.
1619 //--------------------------------------------------------------------------
1620 
1621 static void
1622 plt3zz( PLINT x0, PLINT y0, PLINT dx, PLINT dy, PLINT flag, PLINT *init,
1623  PLFLT_VECTOR x, PLFLT_VECTOR y, PLF2OPS zops, PLPointer zp, PLINT nx, PLINT ny,
1624  PLINT *u, PLINT *v, PLFLT* c )
1625 {
1626  PLINT n = 0;
1627  PLFLT x2d, y2d;
1628  PLFLT ( *getz )( PLPointer, PLINT, PLINT ) = zops->get;
1629 
1630  while ( 1 <= x0 && x0 <= nx && 1 <= y0 && y0 <= ny )
1631  {
1632  x2d = plP_w3wcx( x[x0 - 1], y[y0 - 1], getz( zp, x0 - 1, y0 - 1 ) );
1633  y2d = plP_w3wcy( x[x0 - 1], y[y0 - 1], getz( zp, x0 - 1, y0 - 1 ) );
1634  u[n] = plP_wcpcx( x2d );
1635  v[n] = plP_wcpcy( y2d );
1636  if ( c != NULL )
1637  c[n] = ( getz( zp, x0 - 1, y0 - 1 ) - fc_minz ) / ( fc_maxz - fc_minz );
1638 
1639  switch ( flag )
1640  {
1641  case -3:
1642  y0 += dy;
1643  flag = -flag;
1644  break;
1645  case -2:
1646  y0 += dy;
1647  break;
1648  case -1:
1649  y0 += dy;
1650  flag = -flag;
1651  break;
1652  case 1:
1653  x0 += dx;
1654  break;
1655  case 2:
1656  x0 += dx;
1657  flag = -flag;
1658  break;
1659  case 3:
1660  x0 += dx;
1661  flag = -flag;
1662  break;
1663  }
1664  n++;
1665  }
1666 
1667  if ( flag == 1 || flag == -2 )
1668  {
1669  if ( flag == 1 )
1670  {
1671  x0 -= dx;
1672  y0 += dy;
1673  }
1674  else if ( flag == -2 )
1675  {
1676  y0 -= dy;
1677  x0 += dx;
1678  }
1679  if ( 1 <= x0 && x0 <= nx && 1 <= y0 && y0 <= ny )
1680  {
1681  x2d = plP_w3wcx( x[x0 - 1], y[y0 - 1], getz( zp, x0 - 1, y0 - 1 ) );
1682  y2d = plP_w3wcy( x[x0 - 1], y[y0 - 1], getz( zp, x0 - 1, y0 - 1 ) );
1683  u[n] = plP_wcpcx( x2d );
1684  v[n] = plP_wcpcy( y2d );
1685  if ( c != NULL )
1686  c[n] = ( getz( zp, x0 - 1, y0 - 1 ) - fc_minz ) / ( fc_maxz - fc_minz );
1687  n++;
1688  }
1689  }
1690 
1691 // All the setup is done. Time to do the work.
1692 
1693  plnxtv( u, v, c, n, *init );
1694  *init = 0;
1695 }
1696 
1697 //--------------------------------------------------------------------------
1698 // void plside3()
1699 //
1700 // This routine draws sides around the front of the 3d plot so that
1701 // it does not appear to float.
1702 //--------------------------------------------------------------------------
1703 
1704 static void
1705 plside3( PLFLT_VECTOR x, PLFLT_VECTOR y, PLF2OPS zops, PLPointer zp, PLINT nx, PLINT ny, PLINT opt )
1706 {
1707  PLINT i;
1708  PLFLT cxx, cxy, cyx, cyy, cyz;
1709  PLFLT xmin, ymin, zmin, xmax, ymax, zmax, zscale;
1710  PLFLT tx, ty, ux, uy;
1711  PLFLT ( *getz )( PLPointer, PLINT, PLINT ) = zops->get;
1712 
1713  plP_gw3wc( &cxx, &cxy, &cyx, &cyy, &cyz );
1714  plP_gdom( &xmin, &xmax, &ymin, &ymax );
1715  plP_grange( &zscale, &zmin, &zmax );
1716 
1717 // Get x, y coordinates of legs and plot
1718 
1719  if ( cxx >= 0.0 && cxy <= 0.0 )
1720  {
1721  if ( opt != 1 )
1722  {
1723  for ( i = 0; i < nx; i++ )
1724  {
1725  tx = plP_w3wcx( x[i], y[0], zmin );
1726  ty = plP_w3wcy( x[i], y[0], zmin );
1727  ux = plP_w3wcx( x[i], y[0], getz( zp, i, 0 ) );
1728  uy = plP_w3wcy( x[i], y[0], getz( zp, i, 0 ) );
1729  pljoin( tx, ty, ux, uy );
1730  }
1731  }
1732  if ( opt != 2 )
1733  {
1734  for ( i = 0; i < ny; i++ )
1735  {
1736  tx = plP_w3wcx( x[0], y[i], zmin );
1737  ty = plP_w3wcy( x[0], y[i], zmin );
1738  ux = plP_w3wcx( x[0], y[i], getz( zp, 0, i ) );
1739  uy = plP_w3wcy( x[0], y[i], getz( zp, 0, i ) );
1740  pljoin( tx, ty, ux, uy );
1741  }
1742  }
1743  }
1744  else if ( cxx <= 0.0 && cxy <= 0.0 )
1745  {
1746  if ( opt != 1 )
1747  {
1748  for ( i = 0; i < nx; i++ )
1749  {
1750  tx = plP_w3wcx( x[i], y[ny - 1], zmin );
1751  ty = plP_w3wcy( x[i], y[ny - 1], zmin );
1752  ux = plP_w3wcx( x[i], y[ny - 1], getz( zp, i, ny - 1 ) );
1753  uy = plP_w3wcy( x[i], y[ny - 1], getz( zp, i, ny - 1 ) );
1754  pljoin( tx, ty, ux, uy );
1755  }
1756  }
1757  if ( opt != 2 )
1758  {
1759  for ( i = 0; i < ny; i++ )
1760  {
1761  tx = plP_w3wcx( x[0], y[i], zmin );
1762  ty = plP_w3wcy( x[0], y[i], zmin );
1763  ux = plP_w3wcx( x[0], y[i], getz( zp, 0, i ) );
1764  uy = plP_w3wcy( x[0], y[i], getz( zp, 0, i ) );
1765  pljoin( tx, ty, ux, uy );
1766  }
1767  }
1768  }
1769  else if ( cxx <= 0.0 && cxy >= 0.0 )
1770  {
1771  if ( opt != 1 )
1772  {
1773  for ( i = 0; i < nx; i++ )
1774  {
1775  tx = plP_w3wcx( x[i], y[ny - 1], zmin );
1776  ty = plP_w3wcy( x[i], y[ny - 1], zmin );
1777  ux = plP_w3wcx( x[i], y[ny - 1], getz( zp, i, ny - 1 ) );
1778  uy = plP_w3wcy( x[i], y[ny - 1], getz( zp, i, ny - 1 ) );
1779  pljoin( tx, ty, ux, uy );
1780  }
1781  }
1782  if ( opt != 2 )
1783  {
1784  for ( i = 0; i < ny; i++ )
1785  {
1786  tx = plP_w3wcx( x[nx - 1], y[i], zmin );
1787  ty = plP_w3wcy( x[nx - 1], y[i], zmin );
1788  ux = plP_w3wcx( x[nx - 1], y[i], getz( zp, nx - 1, i ) );
1789  uy = plP_w3wcy( x[nx - 1], y[i], getz( zp, nx - 1, i ) );
1790  pljoin( tx, ty, ux, uy );
1791  }
1792  }
1793  }
1794  else if ( cxx >= 0.0 && cxy >= 0.0 )
1795  {
1796  if ( opt != 1 )
1797  {
1798  for ( i = 0; i < nx; i++ )
1799  {
1800  tx = plP_w3wcx( x[i], y[0], zmin );
1801  ty = plP_w3wcy( x[i], y[0], zmin );
1802  ux = plP_w3wcx( x[i], y[0], getz( zp, i, 0 ) );
1803  uy = plP_w3wcy( x[i], y[0], getz( zp, i, 0 ) );
1804  pljoin( tx, ty, ux, uy );
1805  }
1806  }
1807  if ( opt != 2 )
1808  {
1809  for ( i = 0; i < ny; i++ )
1810  {
1811  tx = plP_w3wcx( x[nx - 1], y[i], zmin );
1812  ty = plP_w3wcy( x[nx - 1], y[i], zmin );
1813  ux = plP_w3wcx( x[nx - 1], y[i], getz( zp, nx - 1, i ) );
1814  uy = plP_w3wcy( x[nx - 1], y[i], getz( zp, nx - 1, i ) );
1815  pljoin( tx, ty, ux, uy );
1816  }
1817  }
1818  }
1819 }
1820 
1821 //--------------------------------------------------------------------------
1822 // void plgrid3()
1823 //
1824 // This routine draws a grid around the back side of the 3d plot with
1825 // hidden line removal.
1826 //--------------------------------------------------------------------------
1827 
1828 static void
1829 plgrid3( PLFLT tick )
1830 {
1831  PLFLT xmin, ymin, zmin, xmax, ymax, zmax, zscale;
1832  PLFLT cxx, cxy, cyx, cyy, cyz, zmin_in, zmax_in;
1833  PLINT u[3], v[3];
1834  PLINT nsub = 0;
1835  PLFLT tp;
1836 
1837  plP_gw3wc( &cxx, &cxy, &cyx, &cyy, &cyz );
1838  plP_gdom( &xmin, &xmax, &ymin, &ymax );
1839  plP_grange( &zscale, &zmin_in, &zmax_in );
1840  zmin = ( zmax_in > zmin_in ) ? zmin_in : zmax_in;
1841  zmax = ( zmax_in > zmin_in ) ? zmax_in : zmin_in;
1842 
1843  pldtik( zmin, zmax, &tick, &nsub, FALSE );
1844  tp = tick * floor( zmin / tick ) + tick;
1845  pl3upv = 0;
1846 
1847  if ( cxx >= 0.0 && cxy <= 0.0 )
1848  {
1849  while ( tp <= zmax )
1850  {
1851  u[0] = plP_wcpcx( plP_w3wcx( xmin, ymax, tp ) );
1852  v[0] = plP_wcpcy( plP_w3wcy( xmin, ymax, tp ) );
1853  u[1] = plP_wcpcx( plP_w3wcx( xmax, ymax, tp ) );
1854  v[1] = plP_wcpcy( plP_w3wcy( xmax, ymax, tp ) );
1855  u[2] = plP_wcpcx( plP_w3wcx( xmax, ymin, tp ) );
1856  v[2] = plP_wcpcy( plP_w3wcy( xmax, ymin, tp ) );
1857  plnxtv( u, v, 0, 3, 0 );
1858 
1859  tp += tick;
1860  }
1861  u[0] = plP_wcpcx( plP_w3wcx( xmax, ymax, zmin ) );
1862  v[0] = plP_wcpcy( plP_w3wcy( xmax, ymax, zmin ) );
1863  u[1] = plP_wcpcx( plP_w3wcx( xmax, ymax, zmax ) );
1864  v[1] = plP_wcpcy( plP_w3wcy( xmax, ymax, zmax ) );
1865  plnxtv( u, v, 0, 2, 0 );
1866  }
1867  else if ( cxx <= 0.0 && cxy <= 0.0 )
1868  {
1869  while ( tp <= zmax )
1870  {
1871  u[0] = plP_wcpcx( plP_w3wcx( xmax, ymax, tp ) );
1872  v[0] = plP_wcpcy( plP_w3wcy( xmax, ymax, tp ) );
1873  u[1] = plP_wcpcx( plP_w3wcx( xmax, ymin, tp ) );
1874  v[1] = plP_wcpcy( plP_w3wcy( xmax, ymin, tp ) );
1875  u[2] = plP_wcpcx( plP_w3wcx( xmin, ymin, tp ) );
1876  v[2] = plP_wcpcy( plP_w3wcy( xmin, ymin, tp ) );
1877  plnxtv( u, v, 0, 3, 0 );
1878 
1879  tp += tick;
1880  }
1881  u[0] = plP_wcpcx( plP_w3wcx( xmax, ymin, zmin ) );
1882  v[0] = plP_wcpcy( plP_w3wcy( xmax, ymin, zmin ) );
1883  u[1] = plP_wcpcx( plP_w3wcx( xmax, ymin, zmax ) );
1884  v[1] = plP_wcpcy( plP_w3wcy( xmax, ymin, zmax ) );
1885  plnxtv( u, v, 0, 2, 0 );
1886  }
1887  else if ( cxx <= 0.0 && cxy >= 0.0 )
1888  {
1889  while ( tp <= zmax )
1890  {
1891  u[0] = plP_wcpcx( plP_w3wcx( xmax, ymin, tp ) );
1892  v[0] = plP_wcpcy( plP_w3wcy( xmax, ymin, tp ) );
1893  u[1] = plP_wcpcx( plP_w3wcx( xmin, ymin, tp ) );
1894  v[1] = plP_wcpcy( plP_w3wcy( xmin, ymin, tp ) );
1895  u[2] = plP_wcpcx( plP_w3wcx( xmin, ymax, tp ) );
1896  v[2] = plP_wcpcy( plP_w3wcy( xmin, ymax, tp ) );
1897  plnxtv( u, v, 0, 3, 0 );
1898 
1899  tp += tick;
1900  }
1901  u[0] = plP_wcpcx( plP_w3wcx( xmin, ymin, zmin ) );
1902  v[0] = plP_wcpcy( plP_w3wcy( xmin, ymin, zmin ) );
1903  u[1] = plP_wcpcx( plP_w3wcx( xmin, ymin, zmax ) );
1904  v[1] = plP_wcpcy( plP_w3wcy( xmin, ymin, zmax ) );
1905  plnxtv( u, v, 0, 2, 0 );
1906  }
1907  else if ( cxx >= 0.0 && cxy >= 0.0 )
1908  {
1909  while ( tp <= zmax )
1910  {
1911  u[0] = plP_wcpcx( plP_w3wcx( xmin, ymin, tp ) );
1912  v[0] = plP_wcpcy( plP_w3wcy( xmin, ymin, tp ) );
1913  u[1] = plP_wcpcx( plP_w3wcx( xmin, ymax, tp ) );
1914  v[1] = plP_wcpcy( plP_w3wcy( xmin, ymax, tp ) );
1915  u[2] = plP_wcpcx( plP_w3wcx( xmax, ymax, tp ) );
1916  v[2] = plP_wcpcy( plP_w3wcy( xmax, ymax, tp ) );
1917  plnxtv( u, v, 0, 3, 0 );
1918 
1919  tp += tick;
1920  }
1921  u[0] = plP_wcpcx( plP_w3wcx( xmin, ymax, zmin ) );
1922  v[0] = plP_wcpcy( plP_w3wcy( xmin, ymax, zmin ) );
1923  u[1] = plP_wcpcx( plP_w3wcx( xmin, ymax, zmax ) );
1924  v[1] = plP_wcpcy( plP_w3wcy( xmin, ymax, zmax ) );
1925  plnxtv( u, v, 0, 2, 0 );
1926  }
1927  pl3upv = 1;
1928 }
1929 
1930 //--------------------------------------------------------------------------
1931 // void plnxtv()
1932 //
1933 // Draw the next view of a 3-d plot. The physical coordinates of the
1934 // points for the next view are placed in the n points of arrays u and
1935 // v. The silhouette found so far is stored in the heap as a set of m peak
1936 // points.
1937 //
1938 // These routines dynamically allocate memory for hidden line removal.
1939 // Memory is allocated in blocks of 2*BINC*sizeof(PLINT) bytes. Large
1940 // values of BINC give better performance but also allocate more memory
1941 // than is needed. If your 3D plots are very "spiky" or you are working
1942 // with very large matrices then you will probably want to increase BINC.
1943 //--------------------------------------------------------------------------
1944 
1945 static void
1946 plnxtv( PLINT *u, PLINT *v, PLFLT* c, PLINT n, PLINT init )
1947 {
1948  plnxtvhi( u, v, c, n, init );
1949 
1950  if ( pl3mode )
1951  plnxtvlo( u, v, c, n, init );
1952 }
1953 
1954 //--------------------------------------------------------------------------
1955 // void plnxtvhi()
1956 //
1957 // Draw the top side of the 3-d plot.
1958 //--------------------------------------------------------------------------
1959 
1960 static void
1961 plnxtvhi( PLINT *u, PLINT *v, PLFLT* c, PLINT n, PLINT init )
1962 {
1963  //
1964  // For the initial set of points, just display them and store them as the
1965  // peak points.
1966  //
1967  if ( init == 1 )
1968  {
1969  int i;
1970  oldhiview = (PLINT *) malloc( (size_t) ( 2 * n ) * sizeof ( PLINT ) );
1971  if ( !oldhiview )
1972  myexit( "plnxtvhi: Out of memory." );
1973 
1974  oldhiview[0] = u[0];
1975  oldhiview[1] = v[0];
1976  plP_draw3d( u[0], v[0], c, 0, 1 );
1977  for ( i = 1; i < n; i++ )
1978  {
1979  oldhiview[2 * i] = u[i];
1980  oldhiview[2 * i + 1] = v[i];
1981  plP_draw3d( u[i], v[i], c, i, 0 );
1982  }
1983  mhi = n;
1984  return;
1985  }
1986 
1987  //
1988  // Otherwise, we need to consider hidden-line removal problem. We scan
1989  // over the points in both the old (i.e. oldhiview[]) and new (i.e. u[]
1990  // and v[]) arrays in order of increasing x coordinate. At each stage, we
1991  // find the line segment in the other array (if one exists) that straddles
1992  // the x coordinate of the point. We have to determine if the point lies
1993  // above or below the line segment, and to check if the below/above status
1994  // has changed since the last point.
1995  //
1996  // If pl3upv = 0 we do not update the view, this is useful for drawing
1997  // lines on the graph after we are done plotting points. Hidden line
1998  // removal is still done, but the view is not updated.
1999  //
2000  xxhi = 0;
2001  if ( pl3upv != 0 )
2002  {
2003  newhisize = 2 * ( mhi + BINC );
2004  if ( newhiview != NULL )
2005  {
2006  newhiview =
2007  (PLINT *) realloc( (void *) newhiview,
2008  (size_t) newhisize * sizeof ( PLINT ) );
2009  }
2010  else
2011  {
2012  newhiview =
2013  (PLINT *) malloc( (size_t) newhisize * sizeof ( PLINT ) );
2014  }
2015  if ( !newhiview )
2016  myexit( "plnxtvhi: Out of memory." );
2017  }
2018 
2019  // Do the draw or shading with hidden line removal
2020 
2021  plnxtvhi_draw( u, v, c, n );
2022 
2023  // Set oldhiview
2024 
2025  swaphiview();
2026 }
2027 
2028 //--------------------------------------------------------------------------
2029 // void plnxtvhi_draw()
2030 //
2031 // Draw the top side of the 3-d plot.
2032 //--------------------------------------------------------------------------
2033 
2034 static void
2035 plnxtvhi_draw( PLINT *u, PLINT *v, PLFLT* c, PLINT n )
2036 {
2037  PLINT i = 0, j = 0, first = 1;
2038  PLINT sx1 = 0, sx2 = 0, sy1 = 0, sy2 = 0;
2039  PLINT su1, su2, sv1, sv2;
2040  PLINT cx, cy, px, py;
2041  PLINT seg, ptold, lstold = 0, pthi, pnewhi = 0, newhi, change, ochange = 0;
2042 
2043 //
2044 // (oldhiview[2*i], oldhiview[2*i]) is the i'th point in the old array
2045 // (u[j], v[j]) is the j'th point in the new array
2046 //
2047 
2048 //
2049 // First attempt at 3d shading. It works ok for simple plots, but
2050 // will just not draw faces, or draw them overlapping for very
2051 // jagged plots
2052 //
2053 
2054  while ( i < mhi || j < n )
2055  {
2056  //
2057  // The coordinates of the point under consideration are (px,py). The
2058  // line segment joins (sx1,sy1) to (sx2,sy2). "ptold" is true if the
2059  // point lies in the old array. We set it by comparing the x coordinates
2060  // of the i'th old point and the j'th new point, being careful if we
2061  // have fallen past the edges. Having found the point, load up the point
2062  // and segment coordinates appropriately.
2063  //
2064 
2065  ptold = ( j >= n || ( i < mhi && oldhiview[2 * i] < u[j] ) );
2066  if ( ptold )
2067  {
2068  px = oldhiview[2 * i];
2069  py = oldhiview[2 * i + 1];
2070  seg = j > 0 && j < n;
2071  if ( seg )
2072  {
2073  sx1 = u[j - 1];
2074  sy1 = v[j - 1];
2075  sx2 = u[j];
2076  sy2 = v[j];
2077  }
2078  }
2079  else
2080  {
2081  px = u[j];
2082  py = v[j];
2083  seg = i > 0 && i < mhi;
2084  if ( seg )
2085  {
2086  sx1 = oldhiview[2 * ( i - 1 )];
2087  sy1 = oldhiview[2 * ( i - 1 ) + 1];
2088  sx2 = oldhiview[2 * i];
2089  sy2 = oldhiview[2 * i + 1];
2090  }
2091  }
2092 
2093  //
2094  // Now determine if the point is higher than the segment, using the
2095  // logical function "above". We also need to know if it is the old view
2096  // or the new view that is higher. "newhi" is set true if the new view
2097  // is higher than the old.
2098  //
2099  if ( seg )
2100  pthi = plabv( px, py, sx1, sy1, sx2, sy2 );
2101  else
2102  pthi = 1;
2103 
2104  newhi = ( ptold && !pthi ) || ( !ptold && pthi );
2105  //
2106  // The last point and this point lie on different sides of
2107  // the current silouette
2108  //
2109  change = ( newhi && !pnewhi ) || ( !newhi && pnewhi );
2110 
2111  //
2112  // There is a new intersection point to put in the peak array if the
2113  // state of "newhi" changes.
2114  //
2115  if ( first )
2116  {
2117  plP_draw3d( px, py, c, j, 1 );
2118  first = 0;
2119  lstold = ptold;
2120  savehipoint( px, py );
2121  pthi = 0;
2122  ochange = 0;
2123  }
2124  else if ( change )
2125  {
2126  //
2127  // Take care of special cases at end of arrays. If pl3upv is 0 the
2128  // endpoints are not connected to the old view.
2129  //
2130  if ( pl3upv == 0 && ( ( !ptold && j == 0 ) || ( ptold && i == 0 ) ) )
2131  {
2132  plP_draw3d( px, py, c, j, 1 );
2133  lstold = ptold;
2134  pthi = 0;
2135  ochange = 0;
2136  }
2137  else if ( pl3upv == 0 &&
2138  ( ( !ptold && i >= mhi ) || ( ptold && j >= n ) ) )
2139  {
2140  plP_draw3d( px, py, c, j, 1 );
2141  lstold = ptold;
2142  pthi = 0;
2143  ochange = 0;
2144  }
2145  else
2146  {
2147  //
2148  // If pl3upv is not 0 then we do want to connect the current line
2149  // with the previous view at the endpoints. Also find intersection
2150  // point with old view.
2151  //
2152  if ( i == 0 )
2153  {
2154  sx1 = oldhiview[0];
2155  sy1 = -1;
2156  sx2 = oldhiview[0];
2157  sy2 = oldhiview[1];
2158  }
2159  else if ( i >= mhi )
2160  {
2161  sx1 = oldhiview[2 * ( mhi - 1 )];
2162  sy1 = oldhiview[2 * ( mhi - 1 ) + 1];
2163  sx2 = oldhiview[2 * ( mhi - 1 )];
2164  sy2 = -1;
2165  }
2166  else
2167  {
2168  sx1 = oldhiview[2 * ( i - 1 )];
2169  sy1 = oldhiview[2 * ( i - 1 ) + 1];
2170  sx2 = oldhiview[2 * i];
2171  sy2 = oldhiview[2 * i + 1];
2172  }
2173 
2174  if ( j == 0 )
2175  {
2176  su1 = u[0];
2177  sv1 = -1;
2178  su2 = u[0];
2179  sv2 = v[0];
2180  }
2181  else if ( j >= n )
2182  {
2183  su1 = u[n - 1];
2184  sv1 = v[n - 1];
2185  su2 = u[n - 1];
2186  sv2 = -1;
2187  }
2188  else
2189  {
2190  su1 = u[j - 1];
2191  sv1 = v[j - 1];
2192  su2 = u[j];
2193  sv2 = v[j];
2194  }
2195 
2196  // Determine the intersection
2197 
2198  pl3cut( sx1, sy1, sx2, sy2, su1, sv1, su2, sv2, &cx, &cy );
2199  if ( cx == px && cy == py )
2200  {
2201  if ( lstold && !ochange )
2202  plP_draw3d( px, py, c, j, 1 );
2203  else
2204  plP_draw3d( px, py, c, j, 0 );
2205 
2206  savehipoint( px, py );
2207  lstold = 1;
2208  pthi = 0;
2209  }
2210  else
2211  {
2212  if ( lstold && !ochange )
2213  plP_draw3d( cx, cy, c, j, 1 );
2214  else
2215  plP_draw3d( cx, cy, c, j, 0 );
2216 
2217  lstold = 1;
2218  savehipoint( cx, cy );
2219  }
2220  ochange = 1;
2221  }
2222  }
2223 
2224  // If point is high then draw plot to point and update view.
2225 
2226  if ( pthi )
2227  {
2228  if ( lstold && ptold )
2229  plP_draw3d( px, py, c, j, 1 );
2230  else
2231  plP_draw3d( px, py, c, j, 0 );
2232 
2233  savehipoint( px, py );
2234  lstold = ptold;
2235  ochange = 0;
2236  }
2237  pnewhi = newhi;
2238 
2239  if ( ptold )
2240  i++;
2241  else
2242  j++;
2243  }
2244 }
2245 
2246 //--------------------------------------------------------------------------
2247 // void plP_draw3d()
2248 //
2249 // Does a simple move or line draw.
2250 //--------------------------------------------------------------------------
2251 
2252 static void
2253 plP_draw3d( PLINT x, PLINT y, PLFLT *c, PLINT j, PLINT move )
2254 {
2255  if ( move )
2256  plP_movphy( x, y );
2257  else
2258  {
2259  if ( c != NULL )
2260  plcol1( c[j - 1] );
2261  plP_draphy( x, y );
2262  }
2263 }
2264 
2265 //--------------------------------------------------------------------------
2266 // void plnxtvlo()
2267 //
2268 // Draw the bottom side of the 3-d plot.
2269 //--------------------------------------------------------------------------
2270 
2271 static void
2272 plnxtvlo( PLINT *u, PLINT *v, PLFLT*c, PLINT n, PLINT init )
2273 {
2274  PLINT i, j, first;
2275  PLINT sx1 = 0, sx2 = 0, sy1 = 0, sy2 = 0;
2276  PLINT su1, su2, sv1, sv2;
2277  PLINT cx, cy, px, py;
2278  PLINT seg, ptold, lstold = 0, ptlo, pnewlo, newlo, change, ochange = 0;
2279 
2280  first = 1;
2281  pnewlo = 0;
2282 
2283  //
2284  // For the initial set of points, just display them and store them as the
2285  // peak points.
2286  //
2287  if ( init == 1 )
2288  {
2289  oldloview = (PLINT *) malloc( (size_t) ( 2 * n ) * sizeof ( PLINT ) );
2290  if ( !oldloview )
2291  myexit( "\nplnxtvlo: Out of memory." );
2292 
2293  plP_draw3d( u[0], v[0], c, 0, 1 );
2294  oldloview[0] = u[0];
2295  oldloview[1] = v[0];
2296  for ( i = 1; i < n; i++ )
2297  {
2298  plP_draw3d( u[i], v[i], c, i, 0 );
2299  oldloview[2 * i] = u[i];
2300  oldloview[2 * i + 1] = v[i];
2301  }
2302  mlo = n;
2303  return;
2304  }
2305 
2306  //
2307  // Otherwise, we need to consider hidden-line removal problem. We scan
2308  // over the points in both the old (i.e. oldloview[]) and new (i.e. u[]
2309  // and v[]) arrays in order of increasing x coordinate. At each stage, we
2310  // find the line segment in the other array (if one exists) that straddles
2311  // the x coordinate of the point. We have to determine if the point lies
2312  // above or below the line segment, and to check if the below/above status
2313  // has changed since the last point.
2314  //
2315  // If pl3upv = 0 we do not update the view, this is useful for drawing
2316  // lines on the graph after we are done plotting points. Hidden line
2317  // removal is still done, but the view is not updated.
2318  //
2319  xxlo = 0;
2320  i = 0;
2321  j = 0;
2322  if ( pl3upv != 0 )
2323  {
2324  newlosize = 2 * ( mlo + BINC );
2325  if ( newloview != NULL )
2326  {
2327  newloview =
2328  (PLINT *) realloc( (void *) newloview,
2329  (size_t) newlosize * sizeof ( PLINT ) );
2330  }
2331  else
2332  {
2333  newloview =
2334  (PLINT *) malloc( (size_t) newlosize * sizeof ( PLINT ) );
2335  }
2336  if ( !newloview )
2337  myexit( "plnxtvlo: Out of memory." );
2338  }
2339 
2340  //
2341  // (oldloview[2*i], oldloview[2*i]) is the i'th point in the old array
2342  // (u[j], v[j]) is the j'th point in the new array.
2343  //
2344  while ( i < mlo || j < n )
2345  {
2346  //
2347  // The coordinates of the point under consideration are (px,py). The
2348  // line segment joins (sx1,sy1) to (sx2,sy2). "ptold" is true if the
2349  // point lies in the old array. We set it by comparing the x coordinates
2350  // of the i'th old point and the j'th new point, being careful if we
2351  // have fallen past the edges. Having found the point, load up the point
2352  // and segment coordinates appropriately.
2353  //
2354  ptold = ( j >= n || ( i < mlo && oldloview[2 * i] < u[j] ) );
2355  if ( ptold )
2356  {
2357  px = oldloview[2 * i];
2358  py = oldloview[2 * i + 1];
2359  seg = j > 0 && j < n;
2360  if ( seg )
2361  {
2362  sx1 = u[j - 1];
2363  sy1 = v[j - 1];
2364  sx2 = u[j];
2365  sy2 = v[j];
2366  }
2367  }
2368  else
2369  {
2370  px = u[j];
2371  py = v[j];
2372  seg = i > 0 && i < mlo;
2373  if ( seg )
2374  {
2375  sx1 = oldloview[2 * ( i - 1 )];
2376  sy1 = oldloview[2 * ( i - 1 ) + 1];
2377  sx2 = oldloview[2 * i];
2378  sy2 = oldloview[2 * i + 1];
2379  }
2380  }
2381 
2382  //
2383  // Now determine if the point is lower than the segment, using the
2384  // logical function "above". We also need to know if it is the old view
2385  // or the new view that is lower. "newlo" is set true if the new view is
2386  // lower than the old.
2387  //
2388  if ( seg )
2389  ptlo = !plabv( px, py, sx1, sy1, sx2, sy2 );
2390  else
2391  ptlo = 1;
2392 
2393  newlo = ( ptold && !ptlo ) || ( !ptold && ptlo );
2394  change = ( newlo && !pnewlo ) || ( !newlo && pnewlo );
2395 
2396  //
2397  // There is a new intersection point to put in the peak array if the
2398  // state of "newlo" changes.
2399  //
2400  if ( first )
2401  {
2402  plP_draw3d( px, py, c, j, 1 );
2403  first = 0;
2404  lstold = ptold;
2405  savelopoint( px, py );
2406  ptlo = 0;
2407  ochange = 0;
2408  }
2409  else if ( change )
2410  {
2411  //
2412  // Take care of special cases at end of arrays. If pl3upv is 0 the
2413  // endpoints are not connected to the old view.
2414  //
2415  if ( pl3upv == 0 && ( ( !ptold && j == 0 ) || ( ptold && i == 0 ) ) )
2416  {
2417  plP_draw3d( px, py, c, j, 1 );
2418  lstold = ptold;
2419  ptlo = 0;
2420  ochange = 0;
2421  }
2422  else if ( pl3upv == 0 &&
2423  ( ( !ptold && i >= mlo ) || ( ptold && j >= n ) ) )
2424  {
2425  plP_draw3d( px, py, c, j, 1 );
2426  lstold = ptold;
2427  ptlo = 0;
2428  ochange = 0;
2429  }
2430 
2431  //
2432  // If pl3upv is not 0 then we do want to connect the current line
2433  // with the previous view at the endpoints. Also find intersection
2434  // point with old view.
2435  //
2436  else
2437  {
2438  if ( i == 0 )
2439  {
2440  sx1 = oldloview[0];
2441  sy1 = 100000;
2442  sx2 = oldloview[0];
2443  sy2 = oldloview[1];
2444  }
2445  else if ( i >= mlo )
2446  {
2447  sx1 = oldloview[2 * ( mlo - 1 )];
2448  sy1 = oldloview[2 * ( mlo - 1 ) + 1];
2449  sx2 = oldloview[2 * ( mlo - 1 )];
2450  sy2 = 100000;
2451  }
2452  else
2453  {
2454  sx1 = oldloview[2 * ( i - 1 )];
2455  sy1 = oldloview[2 * ( i - 1 ) + 1];
2456  sx2 = oldloview[2 * i];
2457  sy2 = oldloview[2 * i + 1];
2458  }
2459 
2460  if ( j == 0 )
2461  {
2462  su1 = u[0];
2463  sv1 = 100000;
2464  su2 = u[0];
2465  sv2 = v[0];
2466  }
2467  else if ( j >= n )
2468  {
2469  su1 = u[n - 1];
2470  sv1 = v[n - 1];
2471  su2 = u[n];
2472  sv2 = 100000;
2473  }
2474  else
2475  {
2476  su1 = u[j - 1];
2477  sv1 = v[j - 1];
2478  su2 = u[j];
2479  sv2 = v[j];
2480  }
2481 
2482  // Determine the intersection
2483 
2484  pl3cut( sx1, sy1, sx2, sy2, su1, sv1, su2, sv2, &cx, &cy );
2485  if ( cx == px && cy == py )
2486  {
2487  if ( lstold && !ochange )
2488  plP_draw3d( px, py, c, j, 1 );
2489  else
2490  plP_draw3d( px, py, c, j, 0 );
2491 
2492  savelopoint( px, py );
2493  lstold = 1;
2494  ptlo = 0;
2495  }
2496  else
2497  {
2498  if ( lstold && !ochange )
2499  plP_draw3d( cx, cy, c, j, 1 );
2500  else
2501  plP_draw3d( cx, cy, c, j, 0 );
2502 
2503  lstold = 1;
2504  savelopoint( cx, cy );
2505  }
2506  ochange = 1;
2507  }
2508  }
2509 
2510  // If point is low then draw plot to point and update view.
2511 
2512  if ( ptlo )
2513  {
2514  if ( lstold && ptold )
2515  plP_draw3d( px, py, c, j, 1 );
2516  else
2517  plP_draw3d( px, py, c, j, 0 );
2518 
2519  savelopoint( px, py );
2520  lstold = ptold;
2521  ochange = 0;
2522  }
2523 
2524  pnewlo = newlo;
2525 
2526  if ( ptold )
2527  i = i + 1;
2528  else
2529  j = j + 1;
2530  }
2531 
2532  // Set oldloview
2533 
2534  swaploview();
2535 }
2536 
2537 //--------------------------------------------------------------------------
2538 // savehipoint
2539 // savelopoint
2540 //
2541 // Add a point to the list of currently visible peaks/valleys, when
2542 // drawing the top/bottom surface, respectively.
2543 //--------------------------------------------------------------------------
2544 
2545 static void
2546 savehipoint( PLINT px, PLINT py )
2547 {
2548  if ( pl3upv == 0 )
2549  return;
2550 
2551  if ( xxhi >= newhisize ) // allocate additional space
2552  {
2553  newhisize += 2 * BINC;
2554  newhiview = (PLINT *) realloc( (void *) newhiview,
2555  (size_t) newhisize * sizeof ( PLINT ) );
2556  if ( !newhiview )
2557  myexit( "savehipoint: Out of memory." );
2558  }
2559 
2560  newhiview[xxhi] = px;
2561  xxhi++;
2562  newhiview[xxhi] = py;
2563  xxhi++;
2564 }
2565 
2566 static void
2567 savelopoint( PLINT px, PLINT py )
2568 {
2569  if ( pl3upv == 0 )
2570  return;
2571 
2572  if ( xxlo >= newlosize ) // allocate additional space
2573  {
2574  newlosize += 2 * BINC;
2575  newloview = (PLINT *) realloc( (void *) newloview,
2576  (size_t) newlosize * sizeof ( PLINT ) );
2577  if ( !newloview )
2578  myexit( "savelopoint: Out of memory." );
2579  }
2580 
2581  newloview[xxlo] = px;
2582  xxlo++;
2583  newloview[xxlo] = py;
2584  xxlo++;
2585 }
2586 
2587 //--------------------------------------------------------------------------
2588 // swaphiview
2589 // swaploview
2590 //
2591 // Swaps the top/bottom views. Need to do a real swap so that the
2592 // memory cleanup routine really frees everything (and only once).
2593 //--------------------------------------------------------------------------
2594 
2595 static void
2596 swaphiview( void )
2597 {
2598  PLINT *tmp;
2599 
2600  if ( pl3upv != 0 )
2601  {
2602  mhi = xxhi / 2;
2603  tmp = oldhiview;
2604  oldhiview = newhiview;
2605  newhiview = tmp;
2606  }
2607 }
2608 
2609 static void
2610 swaploview( void )
2611 {
2612  PLINT *tmp;
2613 
2614  if ( pl3upv != 0 )
2615  {
2616  mlo = xxlo / 2;
2617  tmp = oldloview;
2618  oldloview = newloview;
2619  newloview = tmp;
2620  }
2621 }
2622 
2623 //--------------------------------------------------------------------------
2624 // freework
2625 //
2626 // Frees memory associated with work arrays
2627 //--------------------------------------------------------------------------
2628 
2629 static void
2630 freework( void )
2631 {
2632  free_mem( oldhiview );
2633  free_mem( oldloview );
2634  free_mem( newhiview );
2635  free_mem( newloview );
2636  free_mem( vtmp );
2637  free_mem( utmp );
2638  free_mem( ctmp );
2639 }
2640 
2641 //--------------------------------------------------------------------------
2642 // myexit
2643 //
2644 // Calls plexit, cleaning up first.
2645 //--------------------------------------------------------------------------
2646 
2647 static void
2648 myexit( PLCHAR_VECTOR msg )
2649 {
2650  freework();
2651  plexit( msg );
2652 }
2653 
2654 //--------------------------------------------------------------------------
2655 // myabort
2656 //
2657 // Calls plabort, cleaning up first.
2658 // Caller should return to the user program.
2659 //--------------------------------------------------------------------------
2660 
2661 static void
2662 myabort( PLCHAR_VECTOR msg )
2663 {
2664  freework();
2665  plabort( msg );
2666 }
2667 
2668 //--------------------------------------------------------------------------
2669 // int plabv()
2670 //
2671 // Determines if point (px,py) lies above the line joining (sx1,sy1) to
2672 // (sx2,sy2). It only works correctly if sx1 <= px <= sx2.
2673 //--------------------------------------------------------------------------
2674 
2675 static int
2676 plabv( PLINT px, PLINT py, PLINT sx1, PLINT sy1, PLINT sx2, PLINT sy2 )
2677 {
2678  int above;
2679 
2680  if ( py >= sy1 && py >= sy2 )
2681  above = 1;
2682  else if ( py < sy1 && py < sy2 )
2683  above = 0;
2684  else if ( (double) ( sx2 - sx1 ) * ( py - sy1 ) >=
2685  (double) ( px - sx1 ) * ( sy2 - sy1 ) )
2686  above = 1;
2687  else
2688  above = 0;
2689 
2690  return above;
2691 }
2692 
2693 //--------------------------------------------------------------------------
2694 // void pl3cut()
2695 //
2696 // Determines the point of intersection (cx,cy) between the line joining
2697 // (sx1,sy1) to (sx2,sy2) and the line joining (su1,sv1) to (su2,sv2).
2698 //--------------------------------------------------------------------------
2699 
2700 static void
2701 pl3cut( PLINT sx1, PLINT sy1, PLINT sx2, PLINT sy2,
2702  PLINT su1, PLINT sv1, PLINT su2, PLINT sv2, PLINT *cx, PLINT *cy )
2703 {
2704  PLINT x21, y21, u21, v21, yv1, xu1, a, b;
2705  double fa, fb;
2706 
2707  x21 = sx2 - sx1;
2708  y21 = sy2 - sy1;
2709  u21 = su2 - su1;
2710  v21 = sv2 - sv1;
2711  yv1 = sy1 - sv1;
2712  xu1 = sx1 - su1;
2713 
2714  a = x21 * v21 - y21 * u21;
2715  fa = (double) a;
2716  if ( a == 0 )
2717  {
2718  if ( sx2 < su2 )
2719  {
2720  *cx = sx2;
2721  *cy = sy2;
2722  }
2723  else
2724  {
2725  *cx = su2;
2726  *cy = sv2;
2727  }
2728  return;
2729  }
2730  else
2731  {
2732  b = yv1 * u21 - xu1 * v21;
2733  fb = (double) b;
2734  *cx = (PLINT) ( sx1 + ( fb * x21 ) / fa + .5 );
2735  *cy = (PLINT) ( sy1 + ( fb * y21 ) / fa + .5 );
2736  }
2737 }
2738 
2739 //--------------------------------------------------------------------------
2740 // void plRotationShear
2741 //
2742 // Calculates the rotation and shear angles from a plplot transformation matrix
2743 //
2744 // N.B. the plot transformation matrix
2745 //
2746 // [xFormMatrix[0] xFormMatrix[2]]
2747 // [xFormMatrix[1] xFormMatrix[3]]
2748 //
2749 // is calculated as
2750 //
2751 // [stride cos(t) stride sin(t)]
2752 // [sin(p-t) cos(p-t)]
2753 //
2754 // where t is the rotation angle and p is the shear angle.
2755 // The purpose of this routine is to determine stride, rotation angle,
2756 // and shear angle from xFormMatrix.
2757 //
2758 // For information only, xFormMatrix is the product of the following
2759 // rotation, skew(shear), and scale matrices:
2760 //
2761 // [stride 0] [1 0] [cos(t) sin(t)]
2762 // [0 cos(p)] [tan(p) 1] [-sin(t) cos(t)]
2763 //
2764 //--------------------------------------------------------------------------
2765 
2766 void
2767 plRotationShear( PLFLT *xFormMatrix, PLFLT *rotation, PLFLT *shear, PLFLT *stride )
2768 {
2769  PLFLT smr;
2770  *stride = sqrt( xFormMatrix[0] * xFormMatrix[0] + xFormMatrix[2] * xFormMatrix[2] );
2771 
2772  // Calculate rotation in range from -pi to pi.
2773  *rotation = atan2( xFormMatrix[2], xFormMatrix[0] );
2774 
2775  // Calculate shear - rotation in range from -pi to pi.
2776  smr = atan2( xFormMatrix[1], xFormMatrix[3] );
2777 
2778  // Calculate shear in range from -2 pi to 2 pi.
2779  *shear = smr + *rotation;
2780 
2781  // Calculate shear in range from -pi to pi.
2782  if ( *shear < -PI )
2783  *shear += 2. * PI;
2784  else if ( *shear > PI )
2785  *shear -= 2. * PI;
2786 
2787  // Actually must honour some convention to calculate the negative
2788  // of the shear angle instead of the shear angle. Why??
2789  *shear = -*shear;
2790  // Comment out the modified old logic which determines the negative
2791  // of the shear angle in a more complicated way. Note, the modification
2792  // to divide the asin argument by *stride which solved a long-standing
2793  // bug (as does the above logic in a simpler way).
2794  //
2795  //shear = -asin( (xFormMatrix[0] * xFormMatrix[1] +
2796  // xFormMatrix[2] * xFormMatrix[3] )/ *stride);
2797  //
2798 
2799  // Compute the cross product of the vectors [1,0] and [0,1] to
2800  // determine if we need to make a "quadrant 3,4" adjustment
2801  // to the shear angle.
2802 
2803  //
2804  // if ( xFormMatrix[0] * xFormMatrix[3] - xFormMatrix[1] * xFormMatrix[2] < 0.0 )
2805  // {
2806  //shear = -( M_PI + *shear );
2807  // }
2808  //
2809 }
2810 
plf2ops_t::get
PLFLT(* get)(PLPointer p, PLINT ix, PLINT iy)
Definition: plplot.h:603
plGetAngleToLight
static PLFLT plGetAngleToLight(PLFLT *x, PLFLT *y, PLFLT *z)
Definition: plot3d.c:1568
plnxtv
static void plnxtv(PLINT *, PLINT *, PLFLT *, PLINT, PLINT)
Definition: plot3d.c:1946
savelopoint
static void savelopoint(PLINT, PLINT)
Definition: plot3d.c:2567
plabv
static int plabv(PLINT, PLINT, PLINT, PLINT, PLINT, PLINT)
Definition: plot3d.c:2676
plexit
void plexit(PLCHAR_VECTOR errormsg)
Definition: plctrl.c:1958
c_plmeshc
void c_plmeshc(PLFLT_VECTOR x, PLFLT_VECTOR y, PLFLT_MATRIX z, PLINT nx, PLINT ny, PLINT opt, PLFLT_VECTOR clevel, PLINT nlevel)
Definition: plot3d.c:150
TOP_CONT
#define TOP_CONT
Definition: plplot.h:1508
shade_triangle
static void shade_triangle(PLFLT x0, PLFLT y0, PLFLT z0, PLFLT x1, PLFLT y1, PLFLT z1, PLFLT x2, PLFLT y2, PLFLT z2)
Definition: plot3d.c:262
plP_wcpcx
PLINT plP_wcpcx(PLFLT x)
Definition: plcvt.c:63
vtmp
static PLINT * vtmp
Definition: plot3d.c:43
plfmesh
void plfmesh(PLFLT_VECTOR x, PLFLT_VECTOR y, PLF2OPS zops, PLPointer zp, PLINT nx, PLINT ny, PLINT opt)
Definition: plot3d.c:124
c_plot3dcl
void c_plot3dcl(PLFLT_VECTOR x, PLFLT_VECTOR y, PLFLT_MATRIX z, PLINT nx, PLINT ny, PLINT opt, PLFLT_VECTOR clevel, PLINT nlevel, PLINT indexxmin, PLINT indexxmax, PLINT_VECTOR indexymin, PLINT_VECTOR indexymax)
Definition: plot3d.c:921
SURF_CONT
#define SURF_CONT
Definition: plplot.h:1509
min3
#define min3(a, b, c)
MAG_COLOR
#define MAG_COLOR
Definition: plplot.h:1506
swaploview
static void swaploview(void)
Definition: plot3d.c:2610
NPTS
#define NPTS
PLFLT_MATRIX
const PLFLT *const * PLFLT_MATRIX
Definition: plplot.h:253
newloview
static PLINT * newloview
Definition: plot3d.c:41
c_plsurf3d
void c_plsurf3d(PLFLT_VECTOR x, PLFLT_VECTOR y, PLFLT_MATRIX z, PLINT nx, PLINT ny, PLINT opt, PLFLT_VECTOR clevel, PLINT nlevel)
Definition: plot3d.c:326
plP_grange
void plP_grange(PLFLT *p_zscl, PLFLT *p_zmin, PLFLT *p_zmax)
Definition: plcore.c:4152
PLCHAR_VECTOR
const char * PLCHAR_VECTOR
Definition: plplot.h:243
plf2ops_t
Definition: plplot.h:601
savehipoint
static void savehipoint(PLINT, PLINT)
Definition: plot3d.c:2546
plt3zz
static void plt3zz(PLINT x0, PLINT y0, PLINT dx, PLINT dy, PLINT flag, PLINT *init, PLFLT_VECTOR x, PLFLT_VECTOR y, PLF2OPS zops, PLPointer zp, PLINT nx, PLINT ny, PLINT *u, PLINT *v, PLFLT *c)
Definition: plot3d.c:1622
cont_line::next
struct cont_line * next
Definition: plplotP.h:781
cont_level::line
struct cont_line * line
Definition: plplotP.h:788
MAX
#define MAX(a, b)
Definition: dsplint.c:28
cont_level
Definition: plplotP.h:784
_plm_format::offset
const size_t offset
Definition: plmetafile.c:89
c_plot3dc
void c_plot3dc(PLFLT_VECTOR x, PLFLT_VECTOR y, PLFLT_MATRIX z, PLINT nx, PLINT ny, PLINT opt, PLFLT_VECTOR clevel, PLINT nlevel)
Definition: plot3d.c:883
swaphiview
static void swaphiview(void)
Definition: plot3d.c:2596
c_plot3d
void c_plot3d(PLFLT_VECTOR x, PLFLT_VECTOR y, PLFLT_MATRIX z, PLINT nx, PLINT ny, PLINT opt, PLBOOL side)
Definition: plot3d.c:860
c_plsurf3dl
void c_plsurf3dl(PLFLT_VECTOR x, PLFLT_VECTOR y, PLFLT_MATRIX z, PLINT nx, PLINT ny, PLINT opt, PLFLT_VECTOR clevel, PLINT nlevel, PLINT indexxmin, PLINT indexxmax, PLINT_VECTOR indexymin, PLINT_VECTOR indexymax)
Definition: plot3d.c:389
plabort
void plabort(PLCHAR_VECTOR errormsg)
Definition: plctrl.c:1894
myexit
static void myexit(PLCHAR_VECTOR)
Definition: plot3d.c:2648
c_plmesh
void c_plmesh(PLFLT_VECTOR x, PLFLT_VECTOR y, PLFLT_MATRIX z, PLINT nx, PLINT ny, PLINT opt)
Definition: plot3d.c:118
PLPointer
void * PLPointer
Definition: plplot.h:209
plnxtvlo
static void plnxtvlo(PLINT *, PLINT *, PLFLT *, PLINT, PLINT)
Definition: plot3d.c:2272
pldtik
void pldtik(PLFLT vmin, PLFLT vmax, PLFLT *tick, PLINT *nsubt, PLBOOL ld)
Definition: pldtik.c:36
PLcGrid2::ny
PLINT ny
Definition: plplot.h:521
MESH
#define MESH
Definition: plplot.h:1512
newhiview
static PLINT * newhiview
Definition: plot3d.c:40
plfplot3dc
void plfplot3dc(PLFLT_VECTOR x, PLFLT_VECTOR y, PLF2OPS zops, PLPointer zp, PLINT nx, PLINT ny, PLINT opt, PLFLT_VECTOR clevel, PLINT nlevel)
Definition: plot3d.c:891
color
static int color
Definition: ps.c:78
PLINT
int PLINT
Definition: plplot.h:181
BASE_CONT
#define BASE_CONT
Definition: plplot.h:1507
MIN
#define MIN(a, b)
Definition: dsplint.c:29
PLBOOL
PLINT PLBOOL
Definition: plplot.h:204
cont_line::npts
PLINT npts
Definition: plplotP.h:780
plP_clip_poly
int plP_clip_poly(int Ni, PLFLT *Vi[3], int axis, PLFLT dir, PLFLT offset)
Definition: plot3d.c:166
PLcGrid2::yg
PLFLT_NC_MATRIX yg
Definition: plplot.h:520
pljoin
#define pljoin
Definition: plplot.h:756
FACETED
#define FACETED
Definition: plplot.h:1511
oldloview
static PLINT * oldloview
Definition: plot3d.c:39
max3
#define max3(a, b, c)
plP_draw3d
static void plP_draw3d(PLINT x, PLINT y, PLFLT *c, PLINT j, PLINT move)
Definition: plot3d.c:2253
plfplot3d
void plfplot3d(PLFLT_VECTOR x, PLFLT_VECTOR y, PLF2OPS zops, PLPointer zp, PLINT nx, PLINT ny, PLINT opt, PLBOOL side)
Definition: plot3d.c:867
plP_gzback
void plP_gzback(PLINT **zbf, PLINT **zbc, PLFLT **zbt, PLFLT **zbw)
Definition: plot3d.c:1551
plfsurf3d
void plfsurf3d(PLFLT_VECTOR x, PLFLT_VECTOR y, PLF2OPS zops, PLPointer zp, PLINT nx, PLINT ny, PLINT opt, PLFLT_VECTOR clevel, PLINT nlevel)
Definition: plot3d.c:334
pl3upv
static PLINT pl3upv
Definition: plot3d.c:33
c_pllightsource
void c_pllightsource(PLFLT x, PLFLT y, PLFLT z)
Definition: plot3d.c:101
plfplot3dcl
void plfplot3dcl(PLFLT_VECTOR x, PLFLT_VECTOR y, PLF2OPS zops, PLPointer zp, PLINT nx, PLINT ny, PLINT opt, PLFLT_VECTOR clevel, PLINT nlevel, PLINT PL_UNUSED(indexxmin), PLINT PL_UNUSED(indexxmax), PLINT_VECTOR PL_UNUSED(indexymin), PLINT_VECTOR PL_UNUSED(indexymax))
Definition: plot3d.c:970
cont_line::x
PLFLT * x
Definition: plplotP.h:778
plFree2dGrid
void plFree2dGrid(PLFLT **f, PLINT nx, PLINT PL_UNUSED(ny))
Definition: plmem.c:116
falsecolor
static PLINT falsecolor
Definition: plot3d.c:51
zbflg
static PLINT zbflg
Definition: plot3d.c:35
xxlo
static PLINT xxlo
Definition: plot3d.c:47
newlosize
static PLINT newlosize
Definition: plot3d.c:47
utmp
static PLINT * utmp
Definition: plot3d.c:42
FALSE
#define FALSE
Definition: plplotP.h:177
cont_line
Definition: plplotP.h:776
newhisize
static PLINT newhisize
Definition: plot3d.c:46
PLcGrid2::nx
PLINT nx
Definition: plplot.h:521
xlight
static PLFLT xlight
Definition: plot3d.c:50
plcol1
#define plcol1
Definition: plplot.h:703
mhi
static PLINT mhi
Definition: plot3d.c:46
plnxtvhi
static void plnxtvhi(PLINT *, PLINT *, PLFLT *, PLINT, PLINT)
Definition: plot3d.c:1961
ylight
static PLFLT ylight
Definition: plot3d.c:50
DRAW_SIDES
#define DRAW_SIDES
Definition: plplot.h:1510
plwidth
#define plwidth
Definition: plplot.h:863
PLINT_VECTOR
const PLINT * PLINT_VECTOR
Definition: plplot.h:241
cont_line::y
PLFLT * y
Definition: plplotP.h:779
plnxtvhi_draw
static void plnxtvhi_draw(PLINT *u, PLINT *v, PLFLT *c, PLINT n)
Definition: plot3d.c:2035
DRAW_LINEXY
#define DRAW_LINEXY
Definition: plplot.h:1505
plP_w3wcy
PLFLT plP_w3wcy(PLFLT x, PLFLT y, PLFLT z)
Definition: plcvt.c:222
plfsurf3dl
void plfsurf3dl(PLFLT_VECTOR x, PLFLT_VECTOR y, PLF2OPS zops, PLPointer zp, PLINT nx, PLINT ny, PLINT opt, PLFLT_VECTOR clevel, PLINT nlevel, PLINT indexxmin, PLINT indexxmax, PLINT_VECTOR indexymin, PLINT_VECTOR indexymax)
Definition: plot3d.c:398
plline3
#define plline3
Definition: plplot.h:762
pl3cut
static void pl3cut(PLINT, PLINT, PLINT, PLINT, PLINT, PLINT, PLINT, PLINT, PLINT *, PLINT *)
Definition: plot3d.c:2701
myabort
static void myabort(PLCHAR_VECTOR)
Definition: plot3d.c:2662
plRotationShear
void plRotationShear(PLFLT *xFormMatrix, PLFLT *rotation, PLFLT *shear, PLFLT *stride)
Definition: plot3d.c:2767
plP_draphy
void plP_draphy(PLINT x, PLINT y)
Definition: plline.c:472
DRAW_LINEX
#define DRAW_LINEX
Definition: plplot.h:1503
plP_fill
void plP_fill(short *x, short *y, PLINT npts)
Definition: plcore.c:451
mlo
static PLINT mlo
Definition: plot3d.c:47
PL_UNUSED
#define PL_UNUSED(x)
Definition: plplot.h:138
PLFLT
float PLFLT
Definition: plplot.h:163
plP_w3wcx
PLFLT plP_w3wcx(PLFLT x, PLFLT y, PLFLT PL_UNUSED(z))
Definition: plcvt.c:212
zbwidth
static PLFLT zbwidth
Definition: plot3d.c:36
plplotc.pltr2
def pltr2
Definition: plplotc.py:109
cont_clean_store
void cont_clean_store(CONT_LEVEL *ct)
Definition: plcont.c:173
fc_maxz
static PLFLT fc_maxz
Definition: plot3d.c:52
plP_movphy
void plP_movphy(PLINT x, PLINT y)
Definition: plline.c:459
free_mem
#define free_mem(a)
Definition: plplotP.h:182
plAlloc2dGrid
void plAlloc2dGrid(PLFLT ***f, PLINT nx, PLINT ny)
Definition: plmem.c:91
zlight
static PLFLT zlight
Definition: plot3d.c:50
PI
#define PI
Definition: plplotP.h:290
zbcol
static PLINT zbcol
Definition: plot3d.c:35
plcol0
#define plcol0
Definition: plplot.h:702
plP_wcpcy
PLINT plP_wcpcy(PLFLT y)
Definition: plcvt.c:73
oldhiview
static PLINT * oldhiview
Definition: plot3d.c:38
plP_gw3wc
void plP_gw3wc(PLFLT *p_dxx, PLFLT *p_dxy, PLFLT *p_dyx, PLFLT *p_dyy, PLFLT *p_dyz)
Definition: plcore.c:4162
ctmp
static PLFLT * ctmp
Definition: plot3d.c:44
freework
static void freework(void)
Definition: plot3d.c:2630
plf2ops_c
PLF2OPS plf2ops_c()
Definition: plf2ops.c:126
plwarn
void plwarn(PLCHAR_VECTOR errormsg)
Definition: plctrl.c:1863
plgrid3
static void plgrid3(PLFLT)
Definition: plot3d.c:1829
cont_store
void cont_store(PLFLT_MATRIX f, PLINT nx, PLINT ny, PLINT kx, PLINT lx, PLINT ky, PLINT ly, PLFLT_VECTOR clevel, PLINT nlevel, PLTRANSFORM_callback pltr, PLPointer pltr_data, CONT_LEVEL **contour)
Definition: plcont.c:486
plside3
static void plside3(PLFLT_VECTOR x, PLFLT_VECTOR y, PLF2OPS zops, PLPointer zp, PLINT nx, PLINT ny, PLINT opt)
Definition: plot3d.c:1705
plfmeshc
void plfmeshc(PLFLT_VECTOR x, PLFLT_VECTOR y, PLF2OPS zops, PLPointer zp, PLINT nx, PLINT ny, PLINT opt, PLFLT_VECTOR clevel, PLINT nlevel)
Definition: plot3d.c:157
fc_minz
static PLFLT fc_minz
Definition: plot3d.c:52
PLcGrid2::xg
PLFLT_NC_MATRIX xg
Definition: plplot.h:520
xxhi
static PLINT xxhi
Definition: plot3d.c:46
DRAW_LINEY
#define DRAW_LINEY
Definition: plplot.h:1504
zbtck
static PLFLT zbtck
Definition: plot3d.c:36
cont_level::level
PLFLT level
Definition: plplotP.h:787
PLFLT_VECTOR
const PLFLT * PLFLT_VECTOR
Definition: plplot.h:244
BINC
#define BINC
Definition: plot3d.c:30
pl3mode
static PLINT pl3mode
Definition: plot3d.c:32
cont_level::next
struct cont_level * next
Definition: plplotP.h:789
plP_gdom
void plP_gdom(PLFLT *p_xmin, PLFLT *p_xmax, PLFLT *p_ymin, PLFLT *p_ymax)
Definition: plcore.c:4141
PL_MAXPOLY
#define PL_MAXPOLY
Definition: plplotP.h:283