PLplot  5.11.1
 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 | BASE_CONT"
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  }
717 
718  // After shading completed for a quad, render surface contours.
719  if ( clevel != NULL && ( opt & SURF_CONT ) )
720  {
721  for ( i = 1; i < 3; i++ )
722  {
723  for ( j = 0; j < 4; j += 3 )
724 #define min3( a, b, c ) ( MIN( ( MIN( a, b ) ), c ) )
725 #define max3( a, b, c ) ( MAX( ( MAX( a, b ) ), c ) )
726 
727  {
728  for ( k = 0; k < nlevel; k++ )
729  {
730  if ( clevel[k] >= min3( pz[i], zm, pz[j] ) && clevel[k] < max3( pz[i], zm, pz[j] ) )
731  {
732  ct = 0;
733  if ( clevel[k] >= MIN( pz[i], zm ) && clevel[k] < MAX( pz[i], zm ) ) // p0-pm
734  {
735  xx[ct] = ( ( clevel[k] - pz[i] ) * ( xm - px[i] ) ) / ( zm - pz[i] ) + px[i];
736  yy[ct] = ( ( clevel[k] - pz[i] ) * ( ym - py[i] ) ) / ( zm - pz[i] ) + py[i];
737  ct++;
738  }
739 
740  if ( clevel[k] >= MIN( pz[i], pz[j] ) && clevel[k] < MAX( pz[i], pz[j] ) ) // p0-p1
741  {
742  xx[ct] = ( ( clevel[k] - pz[i] ) * ( px[j] - px[i] ) ) / ( pz[j] - pz[i] ) + px[i];
743  yy[ct] = ( ( clevel[k] - pz[i] ) * ( py[j] - py[i] ) ) / ( pz[j] - pz[i] ) + py[i];
744  ct++;
745  }
746 
747  if ( clevel[k] >= MIN( pz[j], zm ) && clevel[k] < MAX( pz[j], zm ) ) // p1-pm
748  {
749  xx[ct] = ( ( clevel[k] - pz[j] ) * ( xm - px[j] ) ) / ( zm - pz[j] ) + px[j];
750  yy[ct] = ( ( clevel[k] - pz[j] ) * ( ym - py[j] ) ) / ( zm - pz[j] ) + py[j];
751  ct++;
752  }
753 
754  if ( ct == 2 )
755  {
756  // yes, xx and yy are the intersection points of the triangle with
757  // the contour line -- draw a straight line betweeen the points
758  // -- at the end this will make up the contour line
759 
760  // surface contour with color set by user
761  plcol0( color );
762  zz[0] = zz[1] = clevel[k];
763  plline3( 2, xx, yy, zz );
764 
765  // don't break; one triangle can span various contour levels
766  }
767  else
768  plwarn( "plsurf3dl: ***ERROR***\n" );
769  }
770  }
771  }
772  }
773  }
774  }
775  }
776 
777  if ( opt & FACETED )
778  {
779  plcol0( 0 );
780  plfplot3dcl( x, y, zops, zp, nx, ny, MESH | DRAW_LINEXY, NULL, 0,
781  indexxmin, indexxmax, indexymin, indexymax );
782  }
783 
784  if ( opt & DRAW_SIDES ) // the sides look ugly !!!
785  { // draw one more row with all the Z's set to zmin
786  plP_grange( &zscale, &zmin, &zmax );
787 
788  iSlow = nSlow - 1;
789  iftriangle = 1;
790  for ( iFast = 0; iFast < nFast - 1; iFast++ )
791  {
792  for ( i = 0; i < 2; i++ )
793  {
794  ix = ixFast * ( iFast + i ) + ixSlow * iSlow + ixOrigin;
795  iy = iyFast * ( iFast + i ) + iySlow * iSlow + iyOrigin;
796  if ( indexxmin <= ix && ix < indexxmax &&
797  indexymin[ix] <= iy && iy < indexymax[ix] )
798  {
799  px[2 * i] = x[ix];
800  py[2 * i] = y[iy];
801  pz[2 * i] = getz( zp, ix, iy );
802  }
803  else
804  {
805  iftriangle = 0;
806  break;
807  }
808  }
809  if ( iftriangle == 0 )
810  break;
811  // now draw the quad as two triangles (4 might be better)
812 
813  shade_triangle( px[0], py[0], pz[0], px[2], py[2], pz[2], px[0], py[0], zmin );
814  shade_triangle( px[2], py[2], pz[2], px[2], py[2], zmin, px[0], py[0], zmin );
815  }
816 
817  iFast = nFast - 1;
818  iftriangle = 1;
819  for ( iSlow = 0; iSlow < nSlow - 1; iSlow++ )
820  {
821  for ( i = 0; i < 2; i++ )
822  {
823  ix = ixFast * iFast + ixSlow * ( iSlow + i ) + ixOrigin;
824  iy = iyFast * iFast + iySlow * ( iSlow + i ) + iyOrigin;
825  if ( indexxmin <= ix && ix < indexxmax &&
826  indexymin[ix] <= iy && iy < indexymax[ix] )
827  {
828  px[2 * i] = x[ix];
829  py[2 * i] = y[iy];
830  pz[2 * i] = getz( zp, ix, iy );
831  }
832  else
833  {
834  iftriangle = 0;
835  break;
836  }
837  }
838  if ( iftriangle == 0 )
839  break;
840 
841  // now draw the quad as two triangles (4 might be better)
842  shade_triangle( px[0], py[0], pz[0], px[2], py[2], pz[2], px[0], py[0], zmin );
843  shade_triangle( px[2], py[2], pz[2], px[2], py[2], zmin, px[0], py[0], zmin );
844  }
845  }
846 }
847 
848 //--------------------------------------------------------------------------
849 // void plot3d(x, y, z, nx, ny, opt, side)
850 //
851 // Plots a 3-d representation of the function z[x][y]. The x values
852 // are stored as x[0..nx-1], the y values as y[0..ny-1], and the z
853 // values are in the 2-d array z[][]. The integer "opt" specifies:
854 // see plot3dcl() below
855 //--------------------------------------------------------------------------
856 
857 void
858 c_plot3d( const PLFLT *x, const PLFLT *y, const PLFLT * const *z,
859  PLINT nx, PLINT ny, PLINT opt, PLBOOL side )
860 {
861  plfplot3dc( x, y, plf2ops_c(), (PLPointer) z, nx, ny, opt | ( side != 0 ? DRAW_SIDES : 0 ), NULL, 0 );
862 }
863 
864 void
865 plfplot3d( const PLFLT *x, const PLFLT *y, PLF2OPS zops, PLPointer zp,
866  PLINT nx, PLINT ny, PLINT opt, PLBOOL side )
867 {
868  plfplot3dc( x, y, zops, zp, nx, ny, opt | ( side != 0 ? DRAW_SIDES : 0 ), NULL, 0 );
869 }
870 
871 //--------------------------------------------------------------------------
872 // void plot3dc(x, y, z, nx, ny, opt, clevel, nlevel)
873 //
874 // Plots a 3-d representation of the function z[x][y]. The x values
875 // are stored as x[0..nx-1], the y values as y[0..ny-1], and the z
876 // values are in the 2-d array z[][]. The integer "opt" specifies:
877 // see plot3dcl() below
878 //--------------------------------------------------------------------------
879 
880 void
881 c_plot3dc( const PLFLT *x, const PLFLT *y, const PLFLT * const*z,
882  PLINT nx, PLINT ny, PLINT opt,
883  const PLFLT *clevel, PLINT nlevel )
884 {
885  plfplot3dcl( x, y, plf2ops_c(), (PLPointer) z, nx, ny, opt, clevel, nlevel, 0, 0, NULL, NULL );
886 }
887 
888 void
889 plfplot3dc( const PLFLT *x, const PLFLT *y, PLF2OPS zops, PLPointer zp,
890  PLINT nx, PLINT ny, PLINT opt, const PLFLT *clevel, PLINT nlevel )
891 {
892  plfplot3dcl( x, y, zops, zp, nx, ny, opt, clevel, nlevel, 0, 0, NULL, NULL );
893 }
894 
895 //--------------------------------------------------------------------------
896 // void plot3dcl(x, y, z, nx, ny, opt, clevel, nlevel,
897 // indexxmin, indexxmax, indexymin, indexymax)
898 //
899 // Plots a 3-d representation of the function z[x][y]. The x values
900 // are stored as x[0..nx-1], the y values as y[0..ny-1], and the z
901 // values are in the 2-d array z[][]. The integer "opt" specifies:
902 //
903 // DRAW_LINEX : Draw lines parallel to x-axis
904 // DRAW_LINEY : Draw lines parallel to y-axis
905 // DRAW_LINEXY: Draw lines parallel to both axes
906 // MAG_COLOR: Magnitude coloring of wire frame
907 // BASE_CONT: Draw contour at bottom xy plane
908 // TOP_CONT: Draw contour at top xy plane (not yet)
909 // DRAW_SIDES: Draw sides around the plot
910 // MESH: Draw the "under" side of the plot
911 //
912 // or any bitwise combination, e.g. "MAG_COLOR | DRAW_LINEX"
913 // indexymin and indexymax are arrays which specify the y index limits
914 // (following the convention that the upper range limit is one more than
915 // actual index limit) for an x index range of indexxmin, indexxmax.
916 //--------------------------------------------------------------------------
917 
918 void
919 c_plot3dcl( const PLFLT *x, const PLFLT *y, const PLFLT * const *z,
920  PLINT nx, PLINT ny, PLINT opt,
921  const PLFLT *clevel, PLINT nlevel,
922  PLINT indexxmin, PLINT indexxmax, const PLINT *indexymin, const PLINT *indexymax )
923 {
924  plfplot3dcl( x, y, plf2ops_c(), (PLPointer) z, nx, ny,
925  opt, clevel, nlevel, indexxmin, indexxmax, indexymin, indexymax );
926 }
927 
928 //--------------------------------------------------------------------------
965 //--------------------------------------------------------------------------
966 
967 void
968 plfplot3dcl( const PLFLT *x, const PLFLT *y, PLF2OPS zops, PLPointer zp,
969  PLINT nx, PLINT ny, PLINT opt,
970  const PLFLT *clevel, PLINT nlevel,
971  PLINT PL_UNUSED( indexxmin ), PLINT PL_UNUSED( indexxmax ), const PLINT * PL_UNUSED( indexymin ), const PLINT * PL_UNUSED( indexymax ) )
972 {
973  PLFLT cxx, cxy, cyx, cyy, cyz;
974  PLINT init, ix, iy, color;
975  PLFLT width;
976  PLFLT xmin, xmax, ymin, ymax, zmin, zmax, zscale;
977  PLINT ixmin = 0, ixmax = nx - 1, iymin = 0, iymax = ny - 1;
978  PLINT clipped = 0, base_cont = 0, side = 0;
979  PLFLT ( *getz )( PLPointer, PLINT, PLINT ) = zops->get;
980  PLFLT *_x = NULL, *_y = NULL, **_z = NULL;
981  const PLFLT *x_modified, *y_modified;
982  int i;
983 
984  pl3mode = 0;
985 
986  if ( plsc->level < 3 )
987  {
988  myabort( "plot3dcl: Please set up window first" );
989  return;
990  }
991 
992  if ( ( opt & 3 ) == 0 )
993  {
994  myabort( "plot3dcl: Bad option" );
995  return;
996  }
997 
998  if ( nx <= 0 || ny <= 0 )
999  {
1000  myabort( "plot3dcl: Bad array dimensions." );
1001  return;
1002  }
1003 
1004  plP_gdom( &xmin, &xmax, &ymin, &ymax );
1005  plP_grange( &zscale, &zmin, &zmax );
1006  if ( zmin > zmax )
1007  {
1008  PLFLT t = zmin;
1009  zmin = zmax;
1010  zmax = t;
1011  }
1012 
1013 // Check that points in x and in y are strictly increasing
1014 
1015  for ( i = 0; i < nx - 1; i++ )
1016  {
1017  if ( x[i] >= x[i + 1] )
1018  {
1019  myabort( "plot3dcl: X array must be strictly increasing" );
1020  return;
1021  }
1022  }
1023  for ( i = 0; i < ny - 1; i++ )
1024  {
1025  if ( y[i] >= y[i + 1] )
1026  {
1027  myabort( "plot3dcl: Y array must be strictly increasing" );
1028  return;
1029  }
1030  }
1031 
1032  if ( opt & MESH )
1033  pl3mode = 1;
1034 
1035  if ( opt & DRAW_SIDES )
1036  side = 1;
1037 
1038  // figure out the part of the data to use
1039  if ( xmin < x[0] )
1040  xmin = x[0];
1041  if ( xmax > x[nx - 1] )
1042  xmax = x[nx - 1];
1043  if ( ymin < y[0] )
1044  ymin = y[0];
1045  if ( ymax > y[ny - 1] )
1046  ymax = y[ny - 1];
1047  for ( ixmin = 0; ixmin < nx - 1 && x[ixmin + 1] < xmin; ixmin++ )
1048  {
1049  }
1050  for ( ixmax = nx - 1; ixmax > 0 && x[ixmax - 1] > xmax; ixmax-- )
1051  {
1052  }
1053  for ( iymin = 0; iymin < ny - 1 && y[iymin + 1] < ymin; iymin++ )
1054  {
1055  }
1056  for ( iymax = ny - 1; iymax > 0 && y[iymax - 1] > ymax; iymax-- )
1057  {
1058  }
1059  //fprintf(stderr, "(%d,%d) %d %d %d %d\n", nx, ny, ixmin, ixmax, iymin, iymax);
1060  // do we need to clip?
1061  if ( ixmin > 0 || ixmax < nx - 1 || iymin > 0 || iymax < ny - 1 )
1062  {
1063  // adjust the input so it stays within bounds
1064  int _nx = ixmax - ixmin + 1;
1065  int _ny = iymax - iymin + 1;
1066  PLFLT ty0, ty1, tx0, tx1;
1067  int j;
1068 
1069  if ( _nx <= 1 || _ny <= 1 )
1070  {
1071  myabort( "plot3dcl: selected x or y range has no data" );
1072  return;
1073  }
1074 
1075  // allocate storage for new versions of the input vectors
1076  if ( ( ( _x = (PLFLT *) malloc( (size_t) _nx * sizeof ( PLFLT ) ) ) == NULL ) ||
1077  ( ( _y = (PLFLT *) malloc( (size_t) _ny * sizeof ( PLFLT ) ) ) == NULL ) ||
1078  ( ( _z = (PLFLT **) malloc( (size_t) _nx * sizeof ( PLFLT* ) ) ) == NULL ) )
1079  {
1080  plexit( "c_plot3dcl: Insufficient memory" );
1081  }
1082 
1083  clipped = 1;
1084 
1085  // copy over the independent variables
1086  _x[0] = xmin;
1087  _x[_nx - 1] = xmax;
1088  for ( i = 1; i < _nx - 1; i++ )
1089  _x[i] = x[ixmin + i];
1090  _y[0] = ymin;
1091  _y[_ny - 1] = ymax;
1092  for ( i = 1; i < _ny - 1; i++ )
1093  _y[i] = y[iymin + i];
1094 
1095  // copy the data array so we can interpolate around the edges
1096  for ( i = 0; i < _nx; i++ )
1097  {
1098  if ( ( _z[i] = (PLFLT *) malloc( (size_t) _ny * sizeof ( PLFLT ) ) ) == NULL )
1099  {
1100  plexit( "c_plot3dcl: Insufficient memory" );
1101  }
1102  }
1103 
1104  // interpolation factors for the 4 edges
1105  ty0 = ( _y[0] - y[iymin] ) / ( y[iymin + 1] - y[iymin] );
1106  ty1 = ( _y[_ny - 1] - y[iymax - 1] ) / ( y[iymax] - y[iymax - 1] );
1107  tx0 = ( _x[0] - x[ixmin] ) / ( x[ixmin + 1] - x[ixmin] );
1108  tx1 = ( _x[_nx - 1] - x[ixmax - 1] ) / ( x[ixmax] - x[ixmax - 1] );
1109  for ( i = 0; i < _nx; i++ )
1110  {
1111  if ( i == 0 )
1112  {
1113  _z[i][0] = getz( zp, ixmin, iymin ) * ( 1 - ty0 ) * ( 1 - tx0 ) + getz( zp, ixmin, iymin + 1 ) * ( 1 - tx0 ) * ty0
1114  + getz( zp, ixmin + 1, iymin ) * tx0 * ( 1 - ty0 ) + getz( zp, ixmin + 1, iymin + 1 ) * tx0 * ty0;
1115  for ( j = 1; j < _ny - 1; j++ )
1116  _z[i][j] = getz( zp, ixmin, iymin + j ) * ( 1 - tx0 ) + getz( zp, ixmin + 1, iymin + j ) * tx0;
1117  _z[i][_ny - 1] = getz( zp, ixmin, iymax - 1 ) * ( 1 - tx0 ) * ( 1 - ty1 ) + getz( zp, ixmin, iymax ) * ( 1 - tx0 ) * ty1
1118  + getz( zp, ixmin + 1, iymax - 1 ) * tx0 * ( 1 - ty1 ) + getz( zp, ixmin + 1, iymax ) * tx0 * ty1;
1119  }
1120  else if ( i == _nx - 1 )
1121  {
1122  _z[i][0] = getz( zp, ixmax - 1, iymin ) * ( 1 - tx1 ) * ( 1 - ty0 ) + getz( zp, ixmax - 1, iymin + 1 ) * ( 1 - tx1 ) * ty0
1123  + getz( zp, ixmax, iymin ) * tx1 * ( 1 - ty0 ) + getz( zp, ixmax, iymin + 1 ) * tx1 * ty0;
1124  for ( j = 1; j < _ny - 1; j++ )
1125  _z[i][j] = getz( zp, ixmax - 1, iymin + j ) * ( 1 - tx1 ) + getz( zp, ixmax, iymin + j ) * tx1;
1126  _z[i][_ny - 1] = getz( zp, ixmax - 1, iymax - 1 ) * ( 1 - tx1 ) * ( 1 - ty1 ) + getz( zp, ixmax, iymax ) * ( 1 - tx1 ) * ty1
1127  + getz( zp, ixmax, iymax - 1 ) * tx1 * ( 1 - ty1 ) + getz( zp, ixmax, iymax ) * tx1 * ty1;
1128  }
1129  else
1130  {
1131  _z[i][0] = getz( zp, ixmin + i, iymin ) * ( 1 - ty0 ) + getz( zp, ixmin + i, iymin + 1 ) * ty0;
1132  for ( j = 1; j < _ny - 1; j++ )
1133  _z[i][j] = getz( zp, ixmin + i, iymin + j );
1134  _z[i][_ny - 1] = getz( zp, ixmin + i, iymax - 1 ) * ( 1 - ty1 ) + getz( zp, ixmin + i, iymax ) * ty1;
1135  }
1136  for ( j = 0; j < _ny; j++ )
1137  {
1138  if ( _z[i][j] < zmin )
1139  _z[i][j] = zmin;
1140  else if ( _z[i][j] > zmax )
1141  _z[i][j] = zmax;
1142  }
1143  }
1144  // replace the input with our clipped versions
1145  zp = (PLPointer) _z;
1146  getz = plf2ops_c()->get;
1147  nx = _nx;
1148  ny = _ny;
1149  // Do not want to modify input x and y (const modifier)
1150  x_modified = (const PLFLT *) _x;
1151  y_modified = (const PLFLT *) _y;
1152  }
1153  else
1154  {
1155  x_modified = x;
1156  y_modified = y;
1157  }
1158 
1159  // From here on must use x_modified and y_modified rather than
1160  // x and y.
1161  if ( ( opt & BASE_CONT ) || ( opt & TOP_CONT ) || ( opt & MAG_COLOR ) )
1162  {
1163  //
1164  // Don't use the data z value to scale the color, use the z axis
1165  // values set by plw3d()
1166  //
1167  // plMinMax2dGrid(z, nx, ny, &fc_maxz, &fc_minz);
1168  //
1169 
1170  fc_minz = plsc->ranmi;
1171  fc_maxz = plsc->ranma;
1172 
1173  if ( fc_maxz == fc_minz )
1174  {
1175  plwarn( "plot3dcl: Maximum and minimum Z values are equal! \"fixing\"..." );
1176  fc_maxz = fc_minz + 1e-6;
1177  }
1178  }
1179 
1180  if ( opt & BASE_CONT ) // If enabled, draw the contour at the base.
1181  {
1182  if ( clevel != NULL && nlevel != 0 )
1183  {
1184  base_cont = 1;
1185  // even if MESH is not set, "set it",
1186  // as the base contour can only be done in this case
1187  pl3mode = 1;
1188  }
1189  }
1190 
1191  if ( opt & MAG_COLOR ) // If enabled, use magnitude colored wireframe
1192  {
1193  if ( ( ctmp = (PLFLT *) malloc( (size_t) ( 2 * MAX( nx, ny ) ) * sizeof ( PLFLT ) ) ) == NULL )
1194  {
1195  plexit( "c_plot3dcl: Insufficient memory" );
1196  }
1197  }
1198  else
1199  ctmp = NULL;
1200 
1201  // next logic only knows opt = 1 | 2 | 3, make sure that it only gets that
1202  opt &= DRAW_LINEXY;
1203 
1204  // Allocate work arrays
1205 
1206  utmp = (PLINT *) malloc( (size_t) ( 2 * MAX( nx, ny ) ) * sizeof ( PLINT ) );
1207  vtmp = (PLINT *) malloc( (size_t) ( 2 * MAX( nx, ny ) ) * sizeof ( PLINT ) );
1208 
1209  if ( !utmp || !vtmp )
1210  myexit( "plot3dcl: Out of memory." );
1211 
1212  plP_gw3wc( &cxx, &cxy, &cyx, &cyy, &cyz );
1213  init = 1;
1214 // Call 3d line plotter. Each viewing quadrant
1215 // (perpendicular to x-y plane) must be handled separately.
1216  if ( cxx >= 0.0 && cxy <= 0.0 )
1217  {
1218  if ( opt == DRAW_LINEY )
1219  plt3zz( 1, ny, 1, -1, -opt, &init, x_modified, y_modified, zops, zp, nx, ny, utmp, vtmp, ctmp );
1220  else
1221  {
1222  for ( iy = 2; iy <= ny; iy++ )
1223  plt3zz( 1, iy, 1, -1, -opt, &init, x_modified, y_modified, zops, zp, nx, ny, utmp, vtmp, ctmp );
1224  }
1225  if ( opt == DRAW_LINEX )
1226  plt3zz( 1, ny, 1, -1, opt, &init, x_modified, y_modified, zops, zp, nx, ny, utmp, vtmp, ctmp );
1227  else
1228  {
1229  for ( ix = 1; ix <= nx - 1; ix++ )
1230  plt3zz( ix, ny, 1, -1, opt, &init, x_modified, y_modified, zops, zp, nx, ny, utmp, vtmp, ctmp );
1231  }
1232  }
1233 
1234  else if ( cxx <= 0.0 && cxy <= 0.0 )
1235  {
1236  if ( opt == DRAW_LINEX )
1237  plt3zz( nx, ny, -1, -1, opt, &init, x_modified, y_modified, zops, zp, nx, ny, utmp, vtmp, ctmp );
1238  else
1239  {
1240  for ( ix = 2; ix <= nx; ix++ )
1241  plt3zz( ix, ny, -1, -1, opt, &init, x_modified, y_modified, zops, zp, nx, ny, utmp, vtmp, ctmp );
1242  }
1243  if ( opt == DRAW_LINEY )
1244  plt3zz( nx, ny, -1, -1, -opt, &init, x_modified, y_modified, zops, zp, nx, ny, utmp, vtmp, ctmp );
1245  else
1246  {
1247  for ( iy = ny; iy >= 2; iy-- )
1248  plt3zz( nx, iy, -1, -1, -opt, &init, x_modified, y_modified, zops, zp, nx, ny, utmp, vtmp, ctmp );
1249  }
1250  }
1251 
1252  else if ( cxx <= 0.0 && cxy >= 0.0 )
1253  {
1254  if ( opt == DRAW_LINEY )
1255  plt3zz( nx, 1, -1, 1, -opt, &init, x_modified, y_modified, zops, zp, nx, ny, utmp, vtmp, ctmp );
1256  else
1257  {
1258  for ( iy = ny - 1; iy >= 1; iy-- )
1259  plt3zz( nx, iy, -1, 1, -opt, &init, x_modified, y_modified, zops, zp, nx, ny, utmp, vtmp, ctmp );
1260  }
1261  if ( opt == DRAW_LINEX )
1262  plt3zz( nx, 1, -1, 1, opt, &init, x_modified, y_modified, zops, zp, nx, ny, utmp, vtmp, ctmp );
1263  else
1264  {
1265  for ( ix = nx; ix >= 2; ix-- )
1266  plt3zz( ix, 1, -1, 1, opt, &init, x_modified, y_modified, zops, zp, nx, ny, utmp, vtmp, ctmp );
1267  }
1268  }
1269 
1270  else if ( cxx >= 0.0 && cxy >= 0.0 )
1271  {
1272  if ( opt == DRAW_LINEX )
1273  plt3zz( 1, 1, 1, 1, opt, &init, x_modified, y_modified, zops, zp, nx, ny, utmp, vtmp, ctmp );
1274  else
1275  {
1276  for ( ix = nx - 1; ix >= 1; ix-- )
1277  plt3zz( ix, 1, 1, 1, opt, &init, x_modified, y_modified, zops, zp, nx, ny, utmp, vtmp, ctmp );
1278  }
1279  if ( opt == DRAW_LINEY )
1280  plt3zz( 1, 1, 1, 1, -opt, &init, x_modified, y_modified, zops, zp, nx, ny, utmp, vtmp, ctmp );
1281  else
1282  {
1283  for ( iy = 1; iy <= ny - 1; iy++ )
1284  plt3zz( 1, iy, 1, 1, -opt, &init, x_modified, y_modified, zops, zp, nx, ny, utmp, vtmp, ctmp );
1285  }
1286  }
1287 
1288  // draw contour at the base. Not 100%! Why?
1289 
1290  if ( base_cont )
1291  {
1292  int np = NPTS, j;
1293  CONT_LEVEL *cont, *clev;
1294  CONT_LINE *cline;
1295 
1296  PLINT *uu = (PLINT *) malloc( (size_t) NPTS * sizeof ( PLINT ) );
1297  PLINT *vv = (PLINT *) malloc( (size_t) NPTS * sizeof ( PLINT ) );
1298  // prepare cont_store input
1299  PLFLT **zstore;
1300  PLcGrid2 cgrid2;
1301 
1302  if ( ( uu == NULL ) || ( vv == NULL ) )
1303  {
1304  plexit( "c_plot3dcl: Insufficient memory" );
1305  }
1306 
1307  cgrid2.nx = nx;
1308  cgrid2.ny = ny;
1309  plAlloc2dGrid( &cgrid2.xg, nx, ny );
1310  plAlloc2dGrid( &cgrid2.yg, nx, ny );
1311  plAlloc2dGrid( &zstore, nx, ny );
1312 
1313  for ( i = 0; i < nx; i++ )
1314  {
1315  for ( j = 0; j < ny; j++ )
1316  {
1317  cgrid2.xg[i][j] = x_modified[i];
1318  cgrid2.yg[i][j] = y_modified[j];
1319  zstore[i][j] = getz( zp, i, j );
1320  }
1321  }
1322 
1323  pl3upv = 0;
1324 
1325  // Fill cont structure with contours.
1326  cont_store( (const PLFLT * const *) zstore, nx, ny, 1, nx, 1, ny,
1327  clevel, nlevel, pltr2, (void *) &cgrid2, &cont );
1328 
1329  // Free the 2D input arrays to cont_store since not needed any more.
1330  plFree2dGrid( zstore, nx, ny );
1331  plFree2dGrid( cgrid2.xg, nx, ny );
1332  plFree2dGrid( cgrid2.yg, nx, ny );
1333 
1334  // follow the contour levels and lines
1335  clev = cont;
1336  do // for each contour level
1337  {
1338  cline = clev->line;
1339  do // there are several lines that make up each contour
1340  {
1341  int cx, k, l, m, start, end;
1342  PLFLT tx, ty;
1343  if ( cline->npts > np )
1344  {
1345  np = cline->npts;
1346  if ( ( ( uu = (PLINT *) realloc( uu, (size_t) np * sizeof ( PLINT ) ) ) == NULL ) ||
1347  ( ( vv = (PLINT *) realloc( vv, (size_t) np * sizeof ( PLINT ) ) ) == NULL ) )
1348  {
1349  plexit( "c_plot3dcl: Insufficient memory" );
1350  }
1351  }
1352 
1353  // the hidden line plotter plnxtv() only works OK if the x points are in increasing order.
1354  // As this does not always happens, the situation must be detected and the line segment
1355  // must be reversed before being plotted
1356  i = 0;
1357  if ( cline->npts > 0 )
1358  {
1359  do
1360  {
1361  plcol1( ( clev->level - fc_minz ) / ( fc_maxz - fc_minz ) );
1362  cx = plP_wcpcx( plP_w3wcx( cline->x[i], cline->y[i], plsc->ranmi ) );
1363  for ( j = i; j < cline->npts; j++ ) // convert to 2D coordinates
1364  {
1365  uu[j] = plP_wcpcx( plP_w3wcx( cline->x[j], cline->y[j], plsc->ranmi ) );
1366  vv[j] = plP_wcpcy( plP_w3wcy( cline->x[j], cline->y[j], plsc->ranmi ) );
1367  if ( uu[j] < cx ) // find turn back point
1368  break;
1369  else
1370  cx = uu[j];
1371  }
1372  plnxtv( &uu[i], &vv[i], NULL, j - i, 0 ); // plot line with increasing x
1373 
1374  if ( j < cline->npts ) // line not yet finished,
1375  {
1376  start = j - 1;
1377  for ( i = start; i < cline->npts; i++ ) // search turn forward point
1378  {
1379  uu[i] = plP_wcpcx( plP_w3wcx( cline->x[i], cline->y[i], plsc->ranmi ) );
1380  if ( uu[i] > cx )
1381  break;
1382  else
1383  cx = uu[i];
1384  }
1385  end = i - 1;
1386 
1387  for ( k = 0; k < ( end - start + 1 ) / 2; k++ ) // reverse line segment
1388  {
1389  l = start + k;
1390  m = end - k;
1391  tx = cline->x[l];
1392  ty = cline->y[l];
1393  cline->x[l] = cline->x[m];
1394  cline->y[l] = cline->y[m];
1395  cline->x[m] = tx;
1396  cline->y[m] = ty;
1397  }
1398 
1399  // convert to 2D coordinates
1400  for ( j = start; j <= end; j++ )
1401  {
1402  uu[j] = plP_wcpcx( plP_w3wcx( cline->x[j], cline->y[j], plsc->ranmi ) );
1403  vv[j] = plP_wcpcy( plP_w3wcy( cline->x[j], cline->y[j], plsc->ranmi ) );
1404  }
1405  plnxtv( &uu[start], &vv[start], NULL, end - start + 1, 0 ); // and plot it
1406 
1407  cline->x[end] = cline->x[start];
1408  cline->y[end] = cline->y[start];
1409  i = end; // restart where it was left
1410  }
1411  } while ( j < cline->npts && i < cline->npts );
1412  }
1413  cline = cline->next;
1414  }
1415  while ( cline != NULL );
1416  clev = clev->next;
1417  }
1418  while ( clev != NULL );
1419 
1420  cont_clean_store( cont ); // now release the contour memory
1421  pl3upv = 1;
1422  free( uu );
1423  free( vv );
1424  }
1425 
1426 // Finish up by drawing sides, background grid (both are optional)
1427 
1428  if ( side )
1429  plside3( x_modified, y_modified, zops, zp, nx, ny, opt );
1430 
1431  if ( zbflg )
1432  {
1433  color = plsc->icol0;
1434  width = plsc->width;
1435  plwidth( zbwidth );
1436  plcol0( zbcol );
1437  plgrid3( zbtck );
1438  plwidth( width );
1439  plcol0( color );
1440  }
1441 
1442  freework();
1443 
1444  if ( clipped )
1445  {
1446  free( _x );
1447  free( _y );
1448  for ( i = 0; i < nx; i++ )
1449  free( _z[i] );
1450  free( _z );
1451  }
1452 }
1453 
1454 //--------------------------------------------------------------------------
1455 // void plxyindexlimits()
1456 //
1457 // Transform from y array limits to corresponding x array limits (or vice
1458 // versa).
1459 //
1460 // N.B. we follow the convention here that all upper range limits are one
1461 // more than the actual last index.
1462 // instart (>= 0) through inn is the index range where
1463 // the input inarray_min and inarray_max arrays are defined.
1464 //
1465 // outstart (>= 0), through outn (with outn <= outnmax) is the index range
1466 // where the output outarray_min and outarray_max arrays are defined.
1467 //
1468 // In order to assure the transformation from y array limits to x array limits
1469 // (or vice versa) is single-valued, this programme plaborts if the
1470 // inarray_min array has a maximum or inarray_max array has a minimum.
1471 //--------------------------------------------------------------------------
1472 
1473 //static void
1474 //plxyindexlimits( PLINT instart, PLINT inn,
1475 // PLINT *inarray_min, PLINT *inarray_max,
1476 // PLINT *outstart, PLINT *outn, PLINT outnmax,
1477 // PLINT *outarray_min, PLINT *outarray_max )
1478 //{
1479 // PLINT i, j;
1480 // if ( inn < 0 )
1481 // {
1482 // myabort( "plxyindexlimits: Must have instart >= 0" );
1483 // return;
1484 // }
1485 // if ( inn - instart <= 0 )
1486 // {
1487 // myabort( "plxyindexlimits: Must have at least 1 defined point" );
1488 // return;
1489 // }
1490 // *outstart = inarray_min[instart];
1491 // *outn = inarray_max[instart];
1492 // for ( i = instart; i < inn; i++ )
1493 // {
1494 // *outstart = MIN( *outstart, inarray_min[i] );
1495 // *outn = MAX( *outn, inarray_max[i] );
1496 // if ( i + 2 < inn )
1497 // {
1498 // if ( inarray_min[i] < inarray_min[i + 1] &&
1499 // inarray_min[i + 1] > inarray_min[i + 2] )
1500 // {
1501 // myabort( "plxyindexlimits: inarray_min must not have a maximum" );
1502 // return;
1503 // }
1504 // if ( inarray_max[i] > inarray_max[i + 1] &&
1505 // inarray_max[i + 1] < inarray_max[i + 2] )
1506 // {
1507 // myabort( "plxyindexlimits: inarray_max must not have a minimum" );
1508 // return;
1509 // }
1510 // }
1511 // }
1512 // if ( *outstart < 0 )
1513 // {
1514 // myabort( "plxyindexlimits: Must have all elements of inarray_min >= 0" );
1515 // return;
1516 // }
1517 // if ( *outn > outnmax )
1518 // {
1519 // myabort( "plxyindexlimits: Must have all elements of inarray_max <= outnmax" );
1520 // return;
1521 // }
1522 // for ( j = *outstart; j < *outn; j++ )
1523 // {
1524 // i = instart;
1525 // // Find first valid i for this j.
1526 // while ( i < inn && !( inarray_min[i] <= j && j < inarray_max[i] ) )
1527 // i++;
1528 // if ( i < inn )
1529 // outarray_min[j] = i;
1530 // else
1531 // {
1532 // myabort( "plxyindexlimits: bad logic; invalid i should never happen" );
1533 // return;
1534 // }
1535 // // Find next invalid i for this j.
1536 // while ( i < inn && ( inarray_min[i] <= j && j < inarray_max[i] ) )
1537 // i++;
1538 // outarray_max[j] = i;
1539 // }
1540 //}
1541 
1542 //--------------------------------------------------------------------------
1543 // void plP_gzback()
1544 //
1545 // Get background parameters for 3d plot.
1546 //--------------------------------------------------------------------------
1547 
1548 void
1549 plP_gzback( PLINT **zbf, PLINT **zbc, PLFLT **zbt, PLFLT **zbw )
1550 {
1551  *zbf = &zbflg;
1552  *zbc = &zbcol;
1553  *zbt = &zbtck;
1554  *zbw = &zbwidth;
1555 }
1556 
1557 //--------------------------------------------------------------------------
1558 // PLFLT plGetAngleToLight()
1559 //
1560 // Gets cos of angle between normal to a polygon and a light source.
1561 // Requires at least 3 elements, forming non-parallel lines
1562 // in the arrays.
1563 //--------------------------------------------------------------------------
1564 
1565 static PLFLT
1567 {
1568  PLFLT vx1, vx2, vy1, vy2, vz1, vz2;
1569  PLFLT px, py, pz;
1570  PLFLT vlx, vly, vlz;
1571  PLFLT mag1, mag2;
1572  PLFLT cosangle;
1573 
1574  vx1 = x[1] - x[0];
1575  vx2 = x[2] - x[1];
1576  vy1 = y[1] - y[0];
1577  vy2 = y[2] - y[1];
1578  vz1 = z[1] - z[0];
1579  vz2 = z[2] - z[1];
1580 
1581 // Find vector perpendicular to the face
1582  px = vy1 * vz2 - vz1 * vy2;
1583  py = vz1 * vx2 - vx1 * vz2;
1584  pz = vx1 * vy2 - vy1 * vx2;
1585  mag1 = px * px + py * py + pz * pz;
1586 
1587 // Vectors were parallel!
1588  if ( mag1 == 0 )
1589  return 1;
1590 
1591  vlx = xlight - x[0];
1592  vly = ylight - y[0];
1593  vlz = zlight - z[0];
1594  mag2 = vlx * vlx + vly * vly + vlz * vlz;
1595  if ( mag2 == 0 )
1596  return 1;
1597 
1598 // Now have 3 vectors going through the first point on the given surface
1599  cosangle = fabs( ( vlx * px + vly * py + vlz * pz ) / sqrt( mag1 * mag2 ) );
1600 
1601 // In case of numerical rounding
1602  if ( cosangle > 1 )
1603  cosangle = 1;
1604  return cosangle;
1605 }
1606 
1607 //--------------------------------------------------------------------------
1608 // void plt3zz()
1609 //
1610 // Draws the next zig-zag line for a 3-d plot. The data is stored in array
1611 // z[][] as a function of x[] and y[], and is plotted out starting at index
1612 // (x0,y0).
1613 //
1614 // Depending on the state of "flag", the sequence of data points sent to
1615 // plnxtv is altered so as to allow cross-hatch plotting, or plotting
1616 // parallel to either the x-axis or the y-axis.
1617 //--------------------------------------------------------------------------
1618 
1619 static void
1620 plt3zz( PLINT x0, PLINT y0, PLINT dx, PLINT dy, PLINT flag, PLINT *init,
1621  const PLFLT *x, const PLFLT *y, PLF2OPS zops, PLPointer zp, PLINT nx, PLINT ny,
1622  PLINT *u, PLINT *v, PLFLT* c )
1623 {
1624  PLINT n = 0;
1625  PLFLT x2d, y2d;
1626  PLFLT ( *getz )( PLPointer, PLINT, PLINT ) = zops->get;
1627 
1628  while ( 1 <= x0 && x0 <= nx && 1 <= y0 && y0 <= ny )
1629  {
1630  x2d = plP_w3wcx( x[x0 - 1], y[y0 - 1], getz( zp, x0 - 1, y0 - 1 ) );
1631  y2d = plP_w3wcy( x[x0 - 1], y[y0 - 1], getz( zp, x0 - 1, y0 - 1 ) );
1632  u[n] = plP_wcpcx( x2d );
1633  v[n] = plP_wcpcy( y2d );
1634  if ( c != NULL )
1635  c[n] = ( getz( zp, x0 - 1, y0 - 1 ) - fc_minz ) / ( fc_maxz - fc_minz );
1636 
1637  switch ( flag )
1638  {
1639  case -3:
1640  y0 += dy;
1641  flag = -flag;
1642  break;
1643  case -2:
1644  y0 += dy;
1645  break;
1646  case -1:
1647  y0 += dy;
1648  flag = -flag;
1649  break;
1650  case 1:
1651  x0 += dx;
1652  break;
1653  case 2:
1654  x0 += dx;
1655  flag = -flag;
1656  break;
1657  case 3:
1658  x0 += dx;
1659  flag = -flag;
1660  break;
1661  }
1662  n++;
1663  }
1664 
1665  if ( flag == 1 || flag == -2 )
1666  {
1667  if ( flag == 1 )
1668  {
1669  x0 -= dx;
1670  y0 += dy;
1671  }
1672  else if ( flag == -2 )
1673  {
1674  y0 -= dy;
1675  x0 += dx;
1676  }
1677  if ( 1 <= x0 && x0 <= nx && 1 <= y0 && y0 <= ny )
1678  {
1679  x2d = plP_w3wcx( x[x0 - 1], y[y0 - 1], getz( zp, x0 - 1, y0 - 1 ) );
1680  y2d = plP_w3wcy( x[x0 - 1], y[y0 - 1], getz( zp, x0 - 1, y0 - 1 ) );
1681  u[n] = plP_wcpcx( x2d );
1682  v[n] = plP_wcpcy( y2d );
1683  if ( c != NULL )
1684  c[n] = ( getz( zp, x0 - 1, y0 - 1 ) - fc_minz ) / ( fc_maxz - fc_minz );
1685  n++;
1686  }
1687  }
1688 
1689 // All the setup is done. Time to do the work.
1690 
1691  plnxtv( u, v, c, n, *init );
1692  *init = 0;
1693 }
1694 
1695 //--------------------------------------------------------------------------
1696 // void plside3()
1697 //
1698 // This routine draws sides around the front of the 3d plot so that
1699 // it does not appear to float.
1700 //--------------------------------------------------------------------------
1701 
1702 static void
1703 plside3( const PLFLT *x, const PLFLT *y, PLF2OPS zops, PLPointer zp, PLINT nx, PLINT ny, PLINT opt )
1704 {
1705  PLINT i;
1706  PLFLT cxx, cxy, cyx, cyy, cyz;
1707  PLFLT xmin, ymin, zmin, xmax, ymax, zmax, zscale;
1708  PLFLT tx, ty, ux, uy;
1709  PLFLT ( *getz )( PLPointer, PLINT, PLINT ) = zops->get;
1710 
1711  plP_gw3wc( &cxx, &cxy, &cyx, &cyy, &cyz );
1712  plP_gdom( &xmin, &xmax, &ymin, &ymax );
1713  plP_grange( &zscale, &zmin, &zmax );
1714 
1715 // Get x, y coordinates of legs and plot
1716 
1717  if ( cxx >= 0.0 && cxy <= 0.0 )
1718  {
1719  if ( opt != 1 )
1720  {
1721  for ( i = 0; i < nx; i++ )
1722  {
1723  tx = plP_w3wcx( x[i], y[0], zmin );
1724  ty = plP_w3wcy( x[i], y[0], zmin );
1725  ux = plP_w3wcx( x[i], y[0], getz( zp, i, 0 ) );
1726  uy = plP_w3wcy( x[i], y[0], getz( zp, i, 0 ) );
1727  pljoin( tx, ty, ux, uy );
1728  }
1729  }
1730  if ( opt != 2 )
1731  {
1732  for ( i = 0; i < ny; i++ )
1733  {
1734  tx = plP_w3wcx( x[0], y[i], zmin );
1735  ty = plP_w3wcy( x[0], y[i], zmin );
1736  ux = plP_w3wcx( x[0], y[i], getz( zp, 0, i ) );
1737  uy = plP_w3wcy( x[0], y[i], getz( zp, 0, i ) );
1738  pljoin( tx, ty, ux, uy );
1739  }
1740  }
1741  }
1742  else if ( cxx <= 0.0 && cxy <= 0.0 )
1743  {
1744  if ( opt != 1 )
1745  {
1746  for ( i = 0; i < nx; i++ )
1747  {
1748  tx = plP_w3wcx( x[i], y[ny - 1], zmin );
1749  ty = plP_w3wcy( x[i], y[ny - 1], zmin );
1750  ux = plP_w3wcx( x[i], y[ny - 1], getz( zp, i, ny - 1 ) );
1751  uy = plP_w3wcy( x[i], y[ny - 1], getz( zp, i, ny - 1 ) );
1752  pljoin( tx, ty, ux, uy );
1753  }
1754  }
1755  if ( opt != 2 )
1756  {
1757  for ( i = 0; i < ny; i++ )
1758  {
1759  tx = plP_w3wcx( x[0], y[i], zmin );
1760  ty = plP_w3wcy( x[0], y[i], zmin );
1761  ux = plP_w3wcx( x[0], y[i], getz( zp, 0, i ) );
1762  uy = plP_w3wcy( x[0], y[i], getz( zp, 0, i ) );
1763  pljoin( tx, ty, ux, uy );
1764  }
1765  }
1766  }
1767  else if ( cxx <= 0.0 && cxy >= 0.0 )
1768  {
1769  if ( opt != 1 )
1770  {
1771  for ( i = 0; i < nx; i++ )
1772  {
1773  tx = plP_w3wcx( x[i], y[ny - 1], zmin );
1774  ty = plP_w3wcy( x[i], y[ny - 1], zmin );
1775  ux = plP_w3wcx( x[i], y[ny - 1], getz( zp, i, ny - 1 ) );
1776  uy = plP_w3wcy( x[i], y[ny - 1], getz( zp, i, ny - 1 ) );
1777  pljoin( tx, ty, ux, uy );
1778  }
1779  }
1780  if ( opt != 2 )
1781  {
1782  for ( i = 0; i < ny; i++ )
1783  {
1784  tx = plP_w3wcx( x[nx - 1], y[i], zmin );
1785  ty = plP_w3wcy( x[nx - 1], y[i], zmin );
1786  ux = plP_w3wcx( x[nx - 1], y[i], getz( zp, nx - 1, i ) );
1787  uy = plP_w3wcy( x[nx - 1], y[i], getz( zp, nx - 1, i ) );
1788  pljoin( tx, ty, ux, uy );
1789  }
1790  }
1791  }
1792  else if ( cxx >= 0.0 && cxy >= 0.0 )
1793  {
1794  if ( opt != 1 )
1795  {
1796  for ( i = 0; i < nx; i++ )
1797  {
1798  tx = plP_w3wcx( x[i], y[0], zmin );
1799  ty = plP_w3wcy( x[i], y[0], zmin );
1800  ux = plP_w3wcx( x[i], y[0], getz( zp, i, 0 ) );
1801  uy = plP_w3wcy( x[i], y[0], getz( zp, i, 0 ) );
1802  pljoin( tx, ty, ux, uy );
1803  }
1804  }
1805  if ( opt != 2 )
1806  {
1807  for ( i = 0; i < ny; i++ )
1808  {
1809  tx = plP_w3wcx( x[nx - 1], y[i], zmin );
1810  ty = plP_w3wcy( x[nx - 1], y[i], zmin );
1811  ux = plP_w3wcx( x[nx - 1], y[i], getz( zp, nx - 1, i ) );
1812  uy = plP_w3wcy( x[nx - 1], y[i], getz( zp, nx - 1, i ) );
1813  pljoin( tx, ty, ux, uy );
1814  }
1815  }
1816  }
1817 }
1818 
1819 //--------------------------------------------------------------------------
1820 // void plgrid3()
1821 //
1822 // This routine draws a grid around the back side of the 3d plot with
1823 // hidden line removal.
1824 //--------------------------------------------------------------------------
1825 
1826 static void
1828 {
1829  PLFLT xmin, ymin, zmin, xmax, ymax, zmax, zscale;
1830  PLFLT cxx, cxy, cyx, cyy, cyz, zmin_in, zmax_in;
1831  PLINT u[3], v[3];
1832  PLINT nsub = 0;
1833  PLFLT tp;
1834 
1835  plP_gw3wc( &cxx, &cxy, &cyx, &cyy, &cyz );
1836  plP_gdom( &xmin, &xmax, &ymin, &ymax );
1837  plP_grange( &zscale, &zmin_in, &zmax_in );
1838  zmin = ( zmax_in > zmin_in ) ? zmin_in : zmax_in;
1839  zmax = ( zmax_in > zmin_in ) ? zmax_in : zmin_in;
1840 
1841  pldtik( zmin, zmax, &tick, &nsub, FALSE );
1842  tp = tick * floor( zmin / tick ) + tick;
1843  pl3upv = 0;
1844 
1845  if ( cxx >= 0.0 && cxy <= 0.0 )
1846  {
1847  while ( tp <= zmax )
1848  {
1849  u[0] = plP_wcpcx( plP_w3wcx( xmin, ymax, tp ) );
1850  v[0] = plP_wcpcy( plP_w3wcy( xmin, ymax, tp ) );
1851  u[1] = plP_wcpcx( plP_w3wcx( xmax, ymax, tp ) );
1852  v[1] = plP_wcpcy( plP_w3wcy( xmax, ymax, tp ) );
1853  u[2] = plP_wcpcx( plP_w3wcx( xmax, ymin, tp ) );
1854  v[2] = plP_wcpcy( plP_w3wcy( xmax, ymin, tp ) );
1855  plnxtv( u, v, 0, 3, 0 );
1856 
1857  tp += tick;
1858  }
1859  u[0] = plP_wcpcx( plP_w3wcx( xmax, ymax, zmin ) );
1860  v[0] = plP_wcpcy( plP_w3wcy( xmax, ymax, zmin ) );
1861  u[1] = plP_wcpcx( plP_w3wcx( xmax, ymax, zmax ) );
1862  v[1] = plP_wcpcy( plP_w3wcy( xmax, ymax, zmax ) );
1863  plnxtv( u, v, 0, 2, 0 );
1864  }
1865  else if ( cxx <= 0.0 && cxy <= 0.0 )
1866  {
1867  while ( tp <= zmax )
1868  {
1869  u[0] = plP_wcpcx( plP_w3wcx( xmax, ymax, tp ) );
1870  v[0] = plP_wcpcy( plP_w3wcy( xmax, ymax, tp ) );
1871  u[1] = plP_wcpcx( plP_w3wcx( xmax, ymin, tp ) );
1872  v[1] = plP_wcpcy( plP_w3wcy( xmax, ymin, tp ) );
1873  u[2] = plP_wcpcx( plP_w3wcx( xmin, ymin, tp ) );
1874  v[2] = plP_wcpcy( plP_w3wcy( xmin, ymin, tp ) );
1875  plnxtv( u, v, 0, 3, 0 );
1876 
1877  tp += tick;
1878  }
1879  u[0] = plP_wcpcx( plP_w3wcx( xmax, ymin, zmin ) );
1880  v[0] = plP_wcpcy( plP_w3wcy( xmax, ymin, zmin ) );
1881  u[1] = plP_wcpcx( plP_w3wcx( xmax, ymin, zmax ) );
1882  v[1] = plP_wcpcy( plP_w3wcy( xmax, ymin, zmax ) );
1883  plnxtv( u, v, 0, 2, 0 );
1884  }
1885  else if ( cxx <= 0.0 && cxy >= 0.0 )
1886  {
1887  while ( tp <= zmax )
1888  {
1889  u[0] = plP_wcpcx( plP_w3wcx( xmax, ymin, tp ) );
1890  v[0] = plP_wcpcy( plP_w3wcy( xmax, ymin, tp ) );
1891  u[1] = plP_wcpcx( plP_w3wcx( xmin, ymin, tp ) );
1892  v[1] = plP_wcpcy( plP_w3wcy( xmin, ymin, tp ) );
1893  u[2] = plP_wcpcx( plP_w3wcx( xmin, ymax, tp ) );
1894  v[2] = plP_wcpcy( plP_w3wcy( xmin, ymax, tp ) );
1895  plnxtv( u, v, 0, 3, 0 );
1896 
1897  tp += tick;
1898  }
1899  u[0] = plP_wcpcx( plP_w3wcx( xmin, ymin, zmin ) );
1900  v[0] = plP_wcpcy( plP_w3wcy( xmin, ymin, zmin ) );
1901  u[1] = plP_wcpcx( plP_w3wcx( xmin, ymin, zmax ) );
1902  v[1] = plP_wcpcy( plP_w3wcy( xmin, ymin, zmax ) );
1903  plnxtv( u, v, 0, 2, 0 );
1904  }
1905  else if ( cxx >= 0.0 && cxy >= 0.0 )
1906  {
1907  while ( tp <= zmax )
1908  {
1909  u[0] = plP_wcpcx( plP_w3wcx( xmin, ymin, tp ) );
1910  v[0] = plP_wcpcy( plP_w3wcy( xmin, ymin, tp ) );
1911  u[1] = plP_wcpcx( plP_w3wcx( xmin, ymax, tp ) );
1912  v[1] = plP_wcpcy( plP_w3wcy( xmin, ymax, tp ) );
1913  u[2] = plP_wcpcx( plP_w3wcx( xmax, ymax, tp ) );
1914  v[2] = plP_wcpcy( plP_w3wcy( xmax, ymax, tp ) );
1915  plnxtv( u, v, 0, 3, 0 );
1916 
1917  tp += tick;
1918  }
1919  u[0] = plP_wcpcx( plP_w3wcx( xmin, ymax, zmin ) );
1920  v[0] = plP_wcpcy( plP_w3wcy( xmin, ymax, zmin ) );
1921  u[1] = plP_wcpcx( plP_w3wcx( xmin, ymax, zmax ) );
1922  v[1] = plP_wcpcy( plP_w3wcy( xmin, ymax, zmax ) );
1923  plnxtv( u, v, 0, 2, 0 );
1924  }
1925  pl3upv = 1;
1926 }
1927 
1928 //--------------------------------------------------------------------------
1929 // void plnxtv()
1930 //
1931 // Draw the next view of a 3-d plot. The physical coordinates of the
1932 // points for the next view are placed in the n points of arrays u and
1933 // v. The silhouette found so far is stored in the heap as a set of m peak
1934 // points.
1935 //
1936 // These routines dynamically allocate memory for hidden line removal.
1937 // Memory is allocated in blocks of 2*BINC*sizeof(PLINT) bytes. Large
1938 // values of BINC give better performance but also allocate more memory
1939 // than is needed. If your 3D plots are very "spiky" or you are working
1940 // with very large matrices then you will probably want to increase BINC.
1941 //--------------------------------------------------------------------------
1942 
1943 static void
1944 plnxtv( PLINT *u, PLINT *v, PLFLT* c, PLINT n, PLINT init )
1945 {
1946  plnxtvhi( u, v, c, n, init );
1947 
1948  if ( pl3mode )
1949  plnxtvlo( u, v, c, n, init );
1950 }
1951 
1952 //--------------------------------------------------------------------------
1953 // void plnxtvhi()
1954 //
1955 // Draw the top side of the 3-d plot.
1956 //--------------------------------------------------------------------------
1957 
1958 static void
1959 plnxtvhi( PLINT *u, PLINT *v, PLFLT* c, PLINT n, PLINT init )
1960 {
1961  //
1962  // For the initial set of points, just display them and store them as the
1963  // peak points.
1964  //
1965  if ( init == 1 )
1966  {
1967  int i;
1968  oldhiview = (PLINT *) malloc( (size_t) ( 2 * n ) * sizeof ( PLINT ) );
1969  if ( !oldhiview )
1970  myexit( "plnxtvhi: Out of memory." );
1971 
1972  oldhiview[0] = u[0];
1973  oldhiview[1] = v[0];
1974  plP_draw3d( u[0], v[0], c, 0, 1 );
1975  for ( i = 1; i < n; i++ )
1976  {
1977  oldhiview[2 * i] = u[i];
1978  oldhiview[2 * i + 1] = v[i];
1979  plP_draw3d( u[i], v[i], c, i, 0 );
1980  }
1981  mhi = n;
1982  return;
1983  }
1984 
1985  //
1986  // Otherwise, we need to consider hidden-line removal problem. We scan
1987  // over the points in both the old (i.e. oldhiview[]) and new (i.e. u[]
1988  // and v[]) arrays in order of increasing x coordinate. At each stage, we
1989  // find the line segment in the other array (if one exists) that straddles
1990  // the x coordinate of the point. We have to determine if the point lies
1991  // above or below the line segment, and to check if the below/above status
1992  // has changed since the last point.
1993  //
1994  // If pl3upv = 0 we do not update the view, this is useful for drawing
1995  // lines on the graph after we are done plotting points. Hidden line
1996  // removal is still done, but the view is not updated.
1997  //
1998  xxhi = 0;
1999  if ( pl3upv != 0 )
2000  {
2001  newhisize = 2 * ( mhi + BINC );
2002  if ( newhiview != NULL )
2003  {
2004  newhiview =
2005  (PLINT *) realloc( (void *) newhiview,
2006  (size_t) newhisize * sizeof ( PLINT ) );
2007  }
2008  else
2009  {
2010  newhiview =
2011  (PLINT *) malloc( (size_t) newhisize * sizeof ( PLINT ) );
2012  }
2013  if ( !newhiview )
2014  myexit( "plnxtvhi: Out of memory." );
2015  }
2016 
2017  // Do the draw or shading with hidden line removal
2018 
2019  plnxtvhi_draw( u, v, c, n );
2020 
2021  // Set oldhiview
2022 
2023  swaphiview();
2024 }
2025 
2026 //--------------------------------------------------------------------------
2027 // void plnxtvhi_draw()
2028 //
2029 // Draw the top side of the 3-d plot.
2030 //--------------------------------------------------------------------------
2031 
2032 static void
2034 {
2035  PLINT i = 0, j = 0, first = 1;
2036  PLINT sx1 = 0, sx2 = 0, sy1 = 0, sy2 = 0;
2037  PLINT su1, su2, sv1, sv2;
2038  PLINT cx, cy, px, py;
2039  PLINT seg, ptold, lstold = 0, pthi, pnewhi = 0, newhi, change, ochange = 0;
2040 
2041 //
2042 // (oldhiview[2*i], oldhiview[2*i]) is the i'th point in the old array
2043 // (u[j], v[j]) is the j'th point in the new array
2044 //
2045 
2046 //
2047 // First attempt at 3d shading. It works ok for simple plots, but
2048 // will just not draw faces, or draw them overlapping for very
2049 // jagged plots
2050 //
2051 
2052  while ( i < mhi || j < n )
2053  {
2054  //
2055  // The coordinates of the point under consideration are (px,py). The
2056  // line segment joins (sx1,sy1) to (sx2,sy2). "ptold" is true if the
2057  // point lies in the old array. We set it by comparing the x coordinates
2058  // of the i'th old point and the j'th new point, being careful if we
2059  // have fallen past the edges. Having found the point, load up the point
2060  // and segment coordinates appropriately.
2061  //
2062 
2063  ptold = ( j >= n || ( i < mhi && oldhiview[2 * i] < u[j] ) );
2064  if ( ptold )
2065  {
2066  px = oldhiview[2 * i];
2067  py = oldhiview[2 * i + 1];
2068  seg = j > 0 && j < n;
2069  if ( seg )
2070  {
2071  sx1 = u[j - 1];
2072  sy1 = v[j - 1];
2073  sx2 = u[j];
2074  sy2 = v[j];
2075  }
2076  }
2077  else
2078  {
2079  px = u[j];
2080  py = v[j];
2081  seg = i > 0 && i < mhi;
2082  if ( seg )
2083  {
2084  sx1 = oldhiview[2 * ( i - 1 )];
2085  sy1 = oldhiview[2 * ( i - 1 ) + 1];
2086  sx2 = oldhiview[2 * i];
2087  sy2 = oldhiview[2 * i + 1];
2088  }
2089  }
2090 
2091  //
2092  // Now determine if the point is higher than the segment, using the
2093  // logical function "above". We also need to know if it is the old view
2094  // or the new view that is higher. "newhi" is set true if the new view
2095  // is higher than the old.
2096  //
2097  if ( seg )
2098  pthi = plabv( px, py, sx1, sy1, sx2, sy2 );
2099  else
2100  pthi = 1;
2101 
2102  newhi = ( ptold && !pthi ) || ( !ptold && pthi );
2103  //
2104  // The last point and this point lie on different sides of
2105  // the current silouette
2106  //
2107  change = ( newhi && !pnewhi ) || ( !newhi && pnewhi );
2108 
2109  //
2110  // There is a new intersection point to put in the peak array if the
2111  // state of "newhi" changes.
2112  //
2113  if ( first )
2114  {
2115  plP_draw3d( px, py, c, j, 1 );
2116  first = 0;
2117  lstold = ptold;
2118  savehipoint( px, py );
2119  pthi = 0;
2120  ochange = 0;
2121  }
2122  else if ( change )
2123  {
2124  //
2125  // Take care of special cases at end of arrays. If pl3upv is 0 the
2126  // endpoints are not connected to the old view.
2127  //
2128  if ( pl3upv == 0 && ( ( !ptold && j == 0 ) || ( ptold && i == 0 ) ) )
2129  {
2130  plP_draw3d( px, py, c, j, 1 );
2131  lstold = ptold;
2132  pthi = 0;
2133  ochange = 0;
2134  }
2135  else if ( pl3upv == 0 &&
2136  ( ( !ptold && i >= mhi ) || ( ptold && j >= n ) ) )
2137  {
2138  plP_draw3d( px, py, c, j, 1 );
2139  lstold = ptold;
2140  pthi = 0;
2141  ochange = 0;
2142  }
2143  else
2144  {
2145  //
2146  // If pl3upv is not 0 then we do want to connect the current line
2147  // with the previous view at the endpoints. Also find intersection
2148  // point with old view.
2149  //
2150  if ( i == 0 )
2151  {
2152  sx1 = oldhiview[0];
2153  sy1 = -1;
2154  sx2 = oldhiview[0];
2155  sy2 = oldhiview[1];
2156  }
2157  else if ( i >= mhi )
2158  {
2159  sx1 = oldhiview[2 * ( mhi - 1 )];
2160  sy1 = oldhiview[2 * ( mhi - 1 ) + 1];
2161  sx2 = oldhiview[2 * ( mhi - 1 )];
2162  sy2 = -1;
2163  }
2164  else
2165  {
2166  sx1 = oldhiview[2 * ( i - 1 )];
2167  sy1 = oldhiview[2 * ( i - 1 ) + 1];
2168  sx2 = oldhiview[2 * i];
2169  sy2 = oldhiview[2 * i + 1];
2170  }
2171 
2172  if ( j == 0 )
2173  {
2174  su1 = u[0];
2175  sv1 = -1;
2176  su2 = u[0];
2177  sv2 = v[0];
2178  }
2179  else if ( j >= n )
2180  {
2181  su1 = u[n - 1];
2182  sv1 = v[n - 1];
2183  su2 = u[n - 1];
2184  sv2 = -1;
2185  }
2186  else
2187  {
2188  su1 = u[j - 1];
2189  sv1 = v[j - 1];
2190  su2 = u[j];
2191  sv2 = v[j];
2192  }
2193 
2194  // Determine the intersection
2195 
2196  pl3cut( sx1, sy1, sx2, sy2, su1, sv1, su2, sv2, &cx, &cy );
2197  if ( cx == px && cy == py )
2198  {
2199  if ( lstold && !ochange )
2200  plP_draw3d( px, py, c, j, 1 );
2201  else
2202  plP_draw3d( px, py, c, j, 0 );
2203 
2204  savehipoint( px, py );
2205  lstold = 1;
2206  pthi = 0;
2207  }
2208  else
2209  {
2210  if ( lstold && !ochange )
2211  plP_draw3d( cx, cy, c, j, 1 );
2212  else
2213  plP_draw3d( cx, cy, c, j, 0 );
2214 
2215  lstold = 1;
2216  savehipoint( cx, cy );
2217  }
2218  ochange = 1;
2219  }
2220  }
2221 
2222  // If point is high then draw plot to point and update view.
2223 
2224  if ( pthi )
2225  {
2226  if ( lstold && ptold )
2227  plP_draw3d( px, py, c, j, 1 );
2228  else
2229  plP_draw3d( px, py, c, j, 0 );
2230 
2231  savehipoint( px, py );
2232  lstold = ptold;
2233  ochange = 0;
2234  }
2235  pnewhi = newhi;
2236 
2237  if ( ptold )
2238  i++;
2239  else
2240  j++;
2241  }
2242 }
2243 
2244 //--------------------------------------------------------------------------
2245 // void plP_draw3d()
2246 //
2247 // Does a simple move or line draw.
2248 //--------------------------------------------------------------------------
2249 
2250 static void
2251 plP_draw3d( PLINT x, PLINT y, PLFLT *c, PLINT j, PLINT move )
2252 {
2253  if ( move )
2254  plP_movphy( x, y );
2255  else
2256  {
2257  if ( c != NULL )
2258  plcol1( c[j - 1] );
2259  plP_draphy( x, y );
2260  }
2261 }
2262 
2263 //--------------------------------------------------------------------------
2264 // void plnxtvlo()
2265 //
2266 // Draw the bottom side of the 3-d plot.
2267 //--------------------------------------------------------------------------
2268 
2269 static void
2270 plnxtvlo( PLINT *u, PLINT *v, PLFLT*c, PLINT n, PLINT init )
2271 {
2272  PLINT i, j, first;
2273  PLINT sx1 = 0, sx2 = 0, sy1 = 0, sy2 = 0;
2274  PLINT su1, su2, sv1, sv2;
2275  PLINT cx, cy, px, py;
2276  PLINT seg, ptold, lstold = 0, ptlo, pnewlo, newlo, change, ochange = 0;
2277 
2278  first = 1;
2279  pnewlo = 0;
2280 
2281  //
2282  // For the initial set of points, just display them and store them as the
2283  // peak points.
2284  //
2285  if ( init == 1 )
2286  {
2287  oldloview = (PLINT *) malloc( (size_t) ( 2 * n ) * sizeof ( PLINT ) );
2288  if ( !oldloview )
2289  myexit( "\nplnxtvlo: Out of memory." );
2290 
2291  plP_draw3d( u[0], v[0], c, 0, 1 );
2292  oldloview[0] = u[0];
2293  oldloview[1] = v[0];
2294  for ( i = 1; i < n; i++ )
2295  {
2296  plP_draw3d( u[i], v[i], c, i, 0 );
2297  oldloview[2 * i] = u[i];
2298  oldloview[2 * i + 1] = v[i];
2299  }
2300  mlo = n;
2301  return;
2302  }
2303 
2304  //
2305  // Otherwise, we need to consider hidden-line removal problem. We scan
2306  // over the points in both the old (i.e. oldloview[]) and new (i.e. u[]
2307  // and v[]) arrays in order of increasing x coordinate. At each stage, we
2308  // find the line segment in the other array (if one exists) that straddles
2309  // the x coordinate of the point. We have to determine if the point lies
2310  // above or below the line segment, and to check if the below/above status
2311  // has changed since the last point.
2312  //
2313  // If pl3upv = 0 we do not update the view, this is useful for drawing
2314  // lines on the graph after we are done plotting points. Hidden line
2315  // removal is still done, but the view is not updated.
2316  //
2317  xxlo = 0;
2318  i = 0;
2319  j = 0;
2320  if ( pl3upv != 0 )
2321  {
2322  newlosize = 2 * ( mlo + BINC );
2323  if ( newloview != NULL )
2324  {
2325  newloview =
2326  (PLINT *) realloc( (void *) newloview,
2327  (size_t) newlosize * sizeof ( PLINT ) );
2328  }
2329  else
2330  {
2331  newloview =
2332  (PLINT *) malloc( (size_t) newlosize * sizeof ( PLINT ) );
2333  }
2334  if ( !newloview )
2335  myexit( "plnxtvlo: Out of memory." );
2336  }
2337 
2338  //
2339  // (oldloview[2*i], oldloview[2*i]) is the i'th point in the old array
2340  // (u[j], v[j]) is the j'th point in the new array.
2341  //
2342  while ( i < mlo || j < n )
2343  {
2344  //
2345  // The coordinates of the point under consideration are (px,py). The
2346  // line segment joins (sx1,sy1) to (sx2,sy2). "ptold" is true if the
2347  // point lies in the old array. We set it by comparing the x coordinates
2348  // of the i'th old point and the j'th new point, being careful if we
2349  // have fallen past the edges. Having found the point, load up the point
2350  // and segment coordinates appropriately.
2351  //
2352  ptold = ( j >= n || ( i < mlo && oldloview[2 * i] < u[j] ) );
2353  if ( ptold )
2354  {
2355  px = oldloview[2 * i];
2356  py = oldloview[2 * i + 1];
2357  seg = j > 0 && j < n;
2358  if ( seg )
2359  {
2360  sx1 = u[j - 1];
2361  sy1 = v[j - 1];
2362  sx2 = u[j];
2363  sy2 = v[j];
2364  }
2365  }
2366  else
2367  {
2368  px = u[j];
2369  py = v[j];
2370  seg = i > 0 && i < mlo;
2371  if ( seg )
2372  {
2373  sx1 = oldloview[2 * ( i - 1 )];
2374  sy1 = oldloview[2 * ( i - 1 ) + 1];
2375  sx2 = oldloview[2 * i];
2376  sy2 = oldloview[2 * i + 1];
2377  }
2378  }
2379 
2380  //
2381  // Now determine if the point is lower than the segment, using the
2382  // logical function "above". We also need to know if it is the old view
2383  // or the new view that is lower. "newlo" is set true if the new view is
2384  // lower than the old.
2385  //
2386  if ( seg )
2387  ptlo = !plabv( px, py, sx1, sy1, sx2, sy2 );
2388  else
2389  ptlo = 1;
2390 
2391  newlo = ( ptold && !ptlo ) || ( !ptold && ptlo );
2392  change = ( newlo && !pnewlo ) || ( !newlo && pnewlo );
2393 
2394  //
2395  // There is a new intersection point to put in the peak array if the
2396  // state of "newlo" changes.
2397  //
2398  if ( first )
2399  {
2400  plP_draw3d( px, py, c, j, 1 );
2401  first = 0;
2402  lstold = ptold;
2403  savelopoint( px, py );
2404  ptlo = 0;
2405  ochange = 0;
2406  }
2407  else if ( change )
2408  {
2409  //
2410  // Take care of special cases at end of arrays. If pl3upv is 0 the
2411  // endpoints are not connected to the old view.
2412  //
2413  if ( pl3upv == 0 && ( ( !ptold && j == 0 ) || ( ptold && i == 0 ) ) )
2414  {
2415  plP_draw3d( px, py, c, j, 1 );
2416  lstold = ptold;
2417  ptlo = 0;
2418  ochange = 0;
2419  }
2420  else if ( pl3upv == 0 &&
2421  ( ( !ptold && i >= mlo ) || ( ptold && j >= n ) ) )
2422  {
2423  plP_draw3d( px, py, c, j, 1 );
2424  lstold = ptold;
2425  ptlo = 0;
2426  ochange = 0;
2427  }
2428 
2429  //
2430  // If pl3upv is not 0 then we do want to connect the current line
2431  // with the previous view at the endpoints. Also find intersection
2432  // point with old view.
2433  //
2434  else
2435  {
2436  if ( i == 0 )
2437  {
2438  sx1 = oldloview[0];
2439  sy1 = 100000;
2440  sx2 = oldloview[0];
2441  sy2 = oldloview[1];
2442  }
2443  else if ( i >= mlo )
2444  {
2445  sx1 = oldloview[2 * ( mlo - 1 )];
2446  sy1 = oldloview[2 * ( mlo - 1 ) + 1];
2447  sx2 = oldloview[2 * ( mlo - 1 )];
2448  sy2 = 100000;
2449  }
2450  else
2451  {
2452  sx1 = oldloview[2 * ( i - 1 )];
2453  sy1 = oldloview[2 * ( i - 1 ) + 1];
2454  sx2 = oldloview[2 * i];
2455  sy2 = oldloview[2 * i + 1];
2456  }
2457 
2458  if ( j == 0 )
2459  {
2460  su1 = u[0];
2461  sv1 = 100000;
2462  su2 = u[0];
2463  sv2 = v[0];
2464  }
2465  else if ( j >= n )
2466  {
2467  su1 = u[n - 1];
2468  sv1 = v[n - 1];
2469  su2 = u[n];
2470  sv2 = 100000;
2471  }
2472  else
2473  {
2474  su1 = u[j - 1];
2475  sv1 = v[j - 1];
2476  su2 = u[j];
2477  sv2 = v[j];
2478  }
2479 
2480  // Determine the intersection
2481 
2482  pl3cut( sx1, sy1, sx2, sy2, su1, sv1, su2, sv2, &cx, &cy );
2483  if ( cx == px && cy == py )
2484  {
2485  if ( lstold && !ochange )
2486  plP_draw3d( px, py, c, j, 1 );
2487  else
2488  plP_draw3d( px, py, c, j, 0 );
2489 
2490  savelopoint( px, py );
2491  lstold = 1;
2492  ptlo = 0;
2493  }
2494  else
2495  {
2496  if ( lstold && !ochange )
2497  plP_draw3d( cx, cy, c, j, 1 );
2498  else
2499  plP_draw3d( cx, cy, c, j, 0 );
2500 
2501  lstold = 1;
2502  savelopoint( cx, cy );
2503  }
2504  ochange = 1;
2505  }
2506  }
2507 
2508  // If point is low then draw plot to point and update view.
2509 
2510  if ( ptlo )
2511  {
2512  if ( lstold && ptold )
2513  plP_draw3d( px, py, c, j, 1 );
2514  else
2515  plP_draw3d( px, py, c, j, 0 );
2516 
2517  savelopoint( px, py );
2518  lstold = ptold;
2519  ochange = 0;
2520  }
2521 
2522  pnewlo = newlo;
2523 
2524  if ( ptold )
2525  i = i + 1;
2526  else
2527  j = j + 1;
2528  }
2529 
2530  // Set oldloview
2531 
2532  swaploview();
2533 }
2534 
2535 //--------------------------------------------------------------------------
2536 // savehipoint
2537 // savelopoint
2538 //
2539 // Add a point to the list of currently visible peaks/valleys, when
2540 // drawing the top/bottom surface, respectively.
2541 //--------------------------------------------------------------------------
2542 
2543 static void
2545 {
2546  if ( pl3upv == 0 )
2547  return;
2548 
2549  if ( xxhi >= newhisize ) // allocate additional space
2550  {
2551  newhisize += 2 * BINC;
2552  newhiview = (PLINT *) realloc( (void *) newhiview,
2553  (size_t) newhisize * sizeof ( PLINT ) );
2554  if ( !newhiview )
2555  myexit( "savehipoint: Out of memory." );
2556  }
2557 
2558  newhiview[xxhi] = px;
2559  xxhi++;
2560  newhiview[xxhi] = py;
2561  xxhi++;
2562 }
2563 
2564 static void
2566 {
2567  if ( pl3upv == 0 )
2568  return;
2569 
2570  if ( xxlo >= newlosize ) // allocate additional space
2571  {
2572  newlosize += 2 * BINC;
2573  newloview = (PLINT *) realloc( (void *) newloview,
2574  (size_t) newlosize * sizeof ( PLINT ) );
2575  if ( !newloview )
2576  myexit( "savelopoint: Out of memory." );
2577  }
2578 
2579  newloview[xxlo] = px;
2580  xxlo++;
2581  newloview[xxlo] = py;
2582  xxlo++;
2583 }
2584 
2585 //--------------------------------------------------------------------------
2586 // swaphiview
2587 // swaploview
2588 //
2589 // Swaps the top/bottom views. Need to do a real swap so that the
2590 // memory cleanup routine really frees everything (and only once).
2591 //--------------------------------------------------------------------------
2592 
2593 static void
2594 swaphiview( void )
2595 {
2596  PLINT *tmp;
2597 
2598  if ( pl3upv != 0 )
2599  {
2600  mhi = xxhi / 2;
2601  tmp = oldhiview;
2602  oldhiview = newhiview;
2603  newhiview = tmp;
2604  }
2605 }
2606 
2607 static void
2608 swaploview( void )
2609 {
2610  PLINT *tmp;
2611 
2612  if ( pl3upv != 0 )
2613  {
2614  mlo = xxlo / 2;
2615  tmp = oldloview;
2616  oldloview = newloview;
2617  newloview = tmp;
2618  }
2619 }
2620 
2621 //--------------------------------------------------------------------------
2622 // freework
2623 //
2624 // Frees memory associated with work arrays
2625 //--------------------------------------------------------------------------
2626 
2627 static void
2628 freework( void )
2629 {
2630  free_mem( oldhiview );
2631  free_mem( oldloview );
2632  free_mem( newhiview );
2633  free_mem( newloview );
2634  free_mem( vtmp );
2635  free_mem( utmp );
2636  free_mem( ctmp );
2637 }
2638 
2639 //--------------------------------------------------------------------------
2640 // myexit
2641 //
2642 // Calls plexit, cleaning up first.
2643 //--------------------------------------------------------------------------
2644 
2645 static void
2646 myexit( const char *msg )
2647 {
2648  freework();
2649  plexit( msg );
2650 }
2651 
2652 //--------------------------------------------------------------------------
2653 // myabort
2654 //
2655 // Calls plabort, cleaning up first.
2656 // Caller should return to the user program.
2657 //--------------------------------------------------------------------------
2658 
2659 static void
2660 myabort( const char *msg )
2661 {
2662  freework();
2663  plabort( msg );
2664 }
2665 
2666 //--------------------------------------------------------------------------
2667 // int plabv()
2668 //
2669 // Determines if point (px,py) lies above the line joining (sx1,sy1) to
2670 // (sx2,sy2). It only works correctly if sx1 <= px <= sx2.
2671 //--------------------------------------------------------------------------
2672 
2673 static int
2674 plabv( PLINT px, PLINT py, PLINT sx1, PLINT sy1, PLINT sx2, PLINT sy2 )
2675 {
2676  int above;
2677 
2678  if ( py >= sy1 && py >= sy2 )
2679  above = 1;
2680  else if ( py < sy1 && py < sy2 )
2681  above = 0;
2682  else if ( (double) ( sx2 - sx1 ) * ( py - sy1 ) >=
2683  (double) ( px - sx1 ) * ( sy2 - sy1 ) )
2684  above = 1;
2685  else
2686  above = 0;
2687 
2688  return above;
2689 }
2690 
2691 //--------------------------------------------------------------------------
2692 // void pl3cut()
2693 //
2694 // Determines the point of intersection (cx,cy) between the line joining
2695 // (sx1,sy1) to (sx2,sy2) and the line joining (su1,sv1) to (su2,sv2).
2696 //--------------------------------------------------------------------------
2697 
2698 static void
2699 pl3cut( PLINT sx1, PLINT sy1, PLINT sx2, PLINT sy2,
2700  PLINT su1, PLINT sv1, PLINT su2, PLINT sv2, PLINT *cx, PLINT *cy )
2701 {
2702  PLINT x21, y21, u21, v21, yv1, xu1, a, b;
2703  double fa, fb;
2704 
2705  x21 = sx2 - sx1;
2706  y21 = sy2 - sy1;
2707  u21 = su2 - su1;
2708  v21 = sv2 - sv1;
2709  yv1 = sy1 - sv1;
2710  xu1 = sx1 - su1;
2711 
2712  a = x21 * v21 - y21 * u21;
2713  fa = (double) a;
2714  if ( a == 0 )
2715  {
2716  if ( sx2 < su2 )
2717  {
2718  *cx = sx2;
2719  *cy = sy2;
2720  }
2721  else
2722  {
2723  *cx = su2;
2724  *cy = sv2;
2725  }
2726  return;
2727  }
2728  else
2729  {
2730  b = yv1 * u21 - xu1 * v21;
2731  fb = (double) b;
2732  *cx = (PLINT) ( sx1 + ( fb * x21 ) / fa + .5 );
2733  *cy = (PLINT) ( sy1 + ( fb * y21 ) / fa + .5 );
2734  }
2735 }
2736 
2737 //--------------------------------------------------------------------------
2738 // void plRotationShear
2739 //
2740 // Calculates the rotation and shear angles from a plplot transformation matrix
2741 //
2742 // N.B. the plot transformation matrix
2743 //
2744 // [xFormMatrix[0] xFormMatrix[2]]
2745 // [xFormMatrix[1] xFormMatrix[3]]
2746 //
2747 // is calculated as
2748 //
2749 // [stride cos(t) stride sin(t)]
2750 // [sin(p-t) cos(p-t)]
2751 //
2752 // where t is the rotation angle and p is the shear angle.
2753 // The purpose of this routine is to determine stride, rotation angle,
2754 // and shear angle from xFormMatrix.
2755 //
2756 // For information only, xFormMatrix is the product of the following
2757 // rotation, skew(shear), and scale matrices:
2758 //
2759 // [stride 0] [1 0] [cos(t) sin(t)]
2760 // [0 cos(p)] [tan(p) 1] [-sin(t) cos(t)]
2761 //
2762 //--------------------------------------------------------------------------
2763 
2764 void
2765 plRotationShear( PLFLT *xFormMatrix, PLFLT *rotation, PLFLT *shear, PLFLT *stride )
2766 {
2767  PLFLT smr;
2768  *stride = sqrt( xFormMatrix[0] * xFormMatrix[0] + xFormMatrix[2] * xFormMatrix[2] );
2769 
2770  // Calculate rotation in range from -pi to pi.
2771  *rotation = atan2( xFormMatrix[2], xFormMatrix[0] );
2772 
2773  // Calculate shear - rotation in range from -pi to pi.
2774  smr = atan2( xFormMatrix[1], xFormMatrix[3] );
2775 
2776  // Calculate shear in range from -2 pi to 2 pi.
2777  *shear = smr + *rotation;
2778 
2779  // Calculate shear in range from -pi to pi.
2780  if ( *shear < -PI )
2781  *shear += 2. * PI;
2782  else if ( *shear > PI )
2783  *shear -= 2. * PI;
2784 
2785  // Actually must honour some convention to calculate the negative
2786  // of the shear angle instead of the shear angle. Why??
2787  *shear = -*shear;
2788  // Comment out the modified old logic which determines the negative
2789  // of the shear angle in a more complicated way. Note, the modification
2790  // to divide the asin argument by *stride which solved a long-standing
2791  // bug (as does the above logic in a simpler way).
2792  //
2793  //shear = -asin( (xFormMatrix[0] * xFormMatrix[1] +
2794  // xFormMatrix[2] * xFormMatrix[3] )/ *stride);
2795  //
2796 
2797  // Compute the cross product of the vectors [1,0] and [0,1] to
2798  // determine if we need to make a "quadrant 3,4" adjustment
2799  // to the shear angle.
2800 
2801  //
2802  // if ( xFormMatrix[0] * xFormMatrix[3] - xFormMatrix[1] * xFormMatrix[2] < 0.0 )
2803  // {
2804  //shear = -( M_PI + *shear );
2805  // }
2806  //
2807 }
2808