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