PLplot  5.10.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros
plot3d.c
Go to the documentation of this file.
1 // $Id: plot3d.c 12962 2014-01-27 08:07:41Z airwin $
2 //
7 // Copyright (C) 2004 Alan W. Irwin
8 // Copyright (C) 2004 Joao Cardoso
9 // Copyright (C) 2004 Andrew Ross
10 //
11 // This file is part of PLplot.
12 //
13 // PLplot is free software; you can redistribute it and/or modify
14 // it under the terms of the GNU Library General Public License as published
15 // by the Free Software Foundation; either version 2 of the License, or
16 // (at your option) any later version.
17 //
18 // PLplot is distributed in the hope that it will be useful,
19 // but WITHOUT ANY WARRANTY; without even the implied warranty of
20 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
21 // GNU Library General Public License for more details.
22 //
23 // You should have received a copy of the GNU Library General Public License
24 // along with PLplot; if not, write to the Free Software
25 // Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
26 //
27 
28 #include "plplotP.h"
29 
30 // Internal constants
31 
32 #define BINC 50 // Block size for memory allocation
33 
34 static PLINT pl3mode = 0; // 0 3d solid; 1 mesh plot
35 static PLINT pl3upv = 1; // 1 update view; 0 no update
36 
37 static PLINT zbflg = 0, zbcol;
39 
40 static PLINT *oldhiview = NULL;
41 static PLINT *oldloview = NULL;
42 static PLINT *newhiview = NULL;
43 static PLINT *newloview = NULL;
44 static PLINT *utmp = NULL;
45 static PLINT *vtmp = NULL;
46 static PLFLT *ctmp = NULL;
47 
50 
51 // Light source for shading
53 static PLINT falsecolor = 0;
55 
56 // Prototypes for static functions
57 
58 static void plgrid3( PLFLT );
59 static void plnxtv( PLINT *, PLINT *, PLFLT*, PLINT, PLINT );
60 static void plside3( const PLFLT *, const PLFLT *, PLF2OPS, PLPointer, PLINT, PLINT, PLINT );
61 static void plt3zz( PLINT, PLINT, PLINT, PLINT,
62  PLINT, PLINT *, const PLFLT *, const PLFLT *, PLF2OPS, PLPointer,
63  PLINT, PLINT, PLINT *, PLINT *, PLFLT* );
64 static void plnxtvhi( PLINT *, PLINT *, PLFLT*, PLINT, PLINT );
65 static void plnxtvlo( PLINT *, PLINT *, PLFLT*, PLINT, PLINT );
66 static void plnxtvhi_draw( PLINT *u, PLINT *v, PLFLT* c, PLINT n );
67 
68 static void savehipoint( PLINT, PLINT );
69 static void savelopoint( PLINT, PLINT );
70 static void swaphiview( void );
71 static void swaploview( void );
72 static void myexit( const char * );
73 static void myabort( const char * );
74 static void freework( void );
75 static int plabv( PLINT, PLINT, PLINT, PLINT, PLINT, PLINT );
76 static void pl3cut( PLINT, PLINT, PLINT, PLINT, PLINT,
77  PLINT, PLINT, PLINT, PLINT *, PLINT * );
78 static PLFLT plGetAngleToLight( PLFLT* x, PLFLT* y, PLFLT* z );
79 static void plP_draw3d( PLINT x, PLINT y, PLFLT *c, PLINT j, PLINT move );
80 //static void plxyindexlimits( PLINT instart, PLINT inn,
81 // PLINT *inarray_min, PLINT *inarray_max,
82 // PLINT *outstart, PLINT *outn, PLINT outnmax,
83 // PLINT *outarray_min, PLINT *outarray_max );
84 
85 
86 // #define MJL_HACK 1
87 #if MJL_HACK
88 static void plP_fill3( PLINT x0, PLINT y0, PLINT x1, PLINT y1,
89  PLINT x2, PLINT y2, PLINT j );
90 static void plP_fill4( PLINT x0, PLINT y0, PLINT x1, PLINT y1,
91  PLINT x2, PLINT y2, PLINT x3, PLINT y3, PLINT j );
92 #endif
93 
94 //--------------------------------------------------------------------------
95 // void plsetlightsource(x, y, z)
96 //
97 // Sets the position of the light source.
98 //--------------------------------------------------------------------------
99 
100 void
102 {
103  xlight = x;
104  ylight = y;
105  zlight = z;
106 }
107 
108 //--------------------------------------------------------------------------
109 // void plmesh(x, y, z, nx, ny, opt)
110 //
111 // Plots a mesh representation of the function z[x][y]. The x values
112 // are stored as x[0..nx-1], the y values as y[0..ny-1], and the
113 // z values are in the 2-d array z[][]. The integer "opt" specifies:
114 // see plmeshc() below.
115 //--------------------------------------------------------------------------
116 
117 void
118 c_plmesh( const PLFLT *x, const PLFLT *y, const PLFLT * const *z, PLINT nx, PLINT ny, PLINT opt )
119 {
120  plfplot3dc( x, y, plf2ops_c(), (PLPointer) z, nx, ny, opt | MESH, NULL, 0 );
121 }
122 
123 void
124 plfmesh( const PLFLT *x, const PLFLT *y, PLF2OPS zops, PLPointer zp,
125  PLINT nx, PLINT ny, PLINT opt )
126 {
127  plfplot3dc( x, y, zops, zp, nx, ny, opt | MESH, NULL, 0 );
128 }
129 
130 //--------------------------------------------------------------------------
131 // void plmeshc(x, y, z, nx, ny, opt, clevel, nlevel)
132 //
133 // Plots a mesh representation of the function z[x][y]. The x values
134 // are stored as x[0..nx-1], the y values as y[0..ny-1], and the
135 // z values are in the 2-d array z[][]. The integer "opt" specifies:
136 //
137 // DRAW_LINEX draw lines parallel to the X axis
138 // DRAW_LINEY draw lines parallel to the Y axis
139 // DRAW_LINEXY draw lines parallel to both the X and Y axis
140 // MAG_COLOR draw the mesh with a color dependent of the magnitude
141 // BASE_CONT draw contour plot at bottom xy plane
142 // TOP_CONT draw contour plot at top xy plane (not yet)
143 // DRAW_SIDES draw sides
144 //
145 // or any bitwise combination, e.g. "MAG_COLOR | DRAW_LINEX"
146 //
147 //--------------------------------------------------------------------------
148 
149 void
150 c_plmeshc( const PLFLT *x, const PLFLT *y, const PLFLT * const *z, PLINT nx, PLINT ny, PLINT opt,
151  const PLFLT *clevel, PLINT nlevel )
152 {
153  plfplot3dc( x, y, plf2ops_c(), (PLPointer) z, nx, ny, opt | MESH, clevel, nlevel );
154 }
155 
156 void
157 plfmeshc( const PLFLT *x, const PLFLT *y, PLF2OPS zops, PLPointer zp,
158  PLINT nx, PLINT ny, PLINT opt, const PLFLT *clevel, PLINT nlevel )
159 {
160  plfplot3dc( x, y, zops, zp, nx, ny, opt | MESH, clevel, nlevel );
161 }
162 
163 // clipping helper for 3d polygons
164 
165 int
166 plP_clip_poly( int Ni, PLFLT *Vi[3], int axis, PLFLT dir, PLFLT offset )
167 {
168  int anyout = 0;
169  PLFLT _in[PL_MAXPOLY], _T[3][PL_MAXPOLY];
170  PLFLT *in, *T[3], *TT = NULL;
171  int No = 0;
172  int i, j, k;
173 
174  if ( Ni > PL_MAXPOLY )
175  {
176  in = (PLFLT *) malloc( sizeof ( PLFLT ) * (size_t) Ni );
177  TT = (PLFLT *) malloc( 3 * sizeof ( PLFLT ) * (size_t) Ni );
178 
179  if ( in == NULL || TT == NULL )
180  {
181  plexit( "plP_clip_poly: insufficient memory for large polygon" );
182  }
183 
184  T[0] = &TT[0];
185  T[1] = &TT[Ni];
186  T[2] = &TT[2 * Ni];
187  }
188  else
189  {
190  in = _in;
191  T[0] = &_T[0][0];
192  T[1] = &_T[1][0];
193  T[2] = &_T[2][0];
194  }
195 
196  for ( i = 0; i < Ni; i++ )
197  {
198  in[i] = Vi[axis][i] * dir + offset;
199  anyout += in[i] < 0;
200  }
201 
202  // none out
203  if ( anyout == 0 )
204  return Ni;
205 
206  // all out
207  if ( anyout == Ni )
208  {
209  return 0;
210  }
211 
212  // some out
213  // copy over to a temp vector
214  for ( i = 0; i < 3; i++ )
215  {
216  for ( j = 0; j < Ni; j++ )
217  {
218  T[i][j] = Vi[i][j];
219  }
220  }
221  // copy back selectively
222  for ( i = 0; i < Ni; i++ )
223  {
224  j = ( i + 1 ) % Ni;
225 
226  if ( in[i] >= 0 && in[j] >= 0 )
227  {
228  for ( k = 0; k < 3; k++ )
229  Vi[k][No] = T[k][j];
230  No++;
231  }
232  else if ( in[i] >= 0 && in[j] < 0 )
233  {
234  PLFLT u = in[i] / ( in[i] - in[j] );
235  for ( k = 0; k < 3; k++ )
236  Vi[k][No] = T[k][i] * ( 1 - u ) + T[k][j] * u;
237  No++;
238  }
239  else if ( in[i] < 0 && in[j] >= 0 )
240  {
241  PLFLT u = in[i] / ( in[i] - in[j] );
242  for ( k = 0; k < 3; k++ )
243  Vi[k][No] = T[k][i] * ( 1 - u ) + T[k][j] * u;
244  No++;
245  for ( k = 0; k < 3; k++ )
246  Vi[k][No] = T[k][j];
247  No++;
248  }
249  }
250 
251  if ( Ni > PL_MAXPOLY )
252  {
253  free( in );
254  free( TT );
255  }
256 
257  return No;
258 }
259 
260 // helper for plsurf3d, similar to c_plfill3()
261 static void
263  PLFLT x1, PLFLT y1, PLFLT z1,
264  PLFLT x2, PLFLT y2, PLFLT z2 )
265 {
266  int i;
267  // arrays for interface to core functions
268  short u[6], v[6];
269  PLFLT x[6], y[6], z[6];
270  int n;
271  PLFLT xmin, xmax, ymin, ymax, zmin, zmax, zscale;
272  PLFLT *V[3];
273 
274  plP_gdom( &xmin, &xmax, &ymin, &ymax );
275  plP_grange( &zscale, &zmin, &zmax );
276 
277  x[0] = x0; x[1] = x1; x[2] = x2;
278  y[0] = y0; y[1] = y1; y[2] = y2;
279  z[0] = z0; z[1] = z1; z[2] = z2;
280  n = 3;
281 
282  V[0] = x; V[1] = y; V[2] = z;
283 
284  n = plP_clip_poly( n, V, 0, 1, -xmin );
285  n = plP_clip_poly( n, V, 0, -1, xmax );
286  n = plP_clip_poly( n, V, 1, 1, -ymin );
287  n = plP_clip_poly( n, V, 1, -1, ymax );
288  n = plP_clip_poly( n, V, 2, 1, -zmin );
289  n = plP_clip_poly( n, V, 2, -1, zmax );
290 
291  if ( n > 0 )
292  {
293  if ( falsecolor )
294  plcol1( ( ( z[0] + z[1] + z[2] ) / 3. - fc_minz ) / ( fc_maxz - fc_minz ) );
295  else
296  plcol1( plGetAngleToLight( x, y, z ) );
297 
298  for ( i = 0; i < n; i++ )
299  {
300  u[i] = (short) plP_wcpcx( plP_w3wcx( x[i], y[i], z[i] ) );
301  v[i] = (short) plP_wcpcy( plP_w3wcy( x[i], y[i], z[i] ) );
302  }
303  u[n] = u[0];
304  v[n] = v[0];
305 
306 #ifdef SHADE_DEBUG // show triangles
307  plcol0( 1 );
308  x[3] = x[0]; y[3] = y[0]; z[3] = z[0];
309  plline3( 4, x, y, z );
310 #else // fill triangles
311  plP_fill( u, v, n + 1 );
312 #endif
313  }
314 }
315 
316 //--------------------------------------------------------------------------
317 // void plsurf3d(x, y, z, nx, ny, opt, clevel, nlevel)
318 //
319 // Plots the 3-d surface representation of the function z[x][y].
320 // The x values are stored as x[0..nx-1], the y values as y[0..ny-1],
321 // and the z values are in the 2-d array z[][]. The integer "opt" specifies:
322 // see plsurf3dl() below.
323 //--------------------------------------------------------------------------
324 
325 void
326 c_plsurf3d( const PLFLT *x, const PLFLT *y, const PLFLT * const *z, PLINT nx, PLINT ny,
327  PLINT opt, const PLFLT *clevel, PLINT nlevel )
328 {
329  plfsurf3d( x, y, plf2ops_c(), (PLPointer) z, nx, ny,
330  opt, clevel, nlevel );
331 }
332 
333 void
334 plfsurf3d( const PLFLT *x, const PLFLT *y, PLF2OPS zops, PLPointer zp,
335  PLINT nx, PLINT ny, PLINT opt, const PLFLT *clevel, PLINT nlevel )
336 {
337  PLINT i;
338  PLINT *indexymin = (PLINT *) malloc( (size_t) nx * sizeof ( PLINT ) );
339  PLINT *indexymax = (PLINT *) malloc( (size_t) nx * sizeof ( PLINT ) );
340 
341  if ( !indexymin || !indexymax )
342  plexit( "plsurf3d: Out of memory." );
343  for ( i = 0; i < nx; i++ )
344  {
345  indexymin[i] = 0;
346  indexymax[i] = ny;
347  }
348  plfsurf3dl( x, y, zops, zp, nx, ny, opt, clevel, nlevel,
349  0, nx, indexymin, indexymax );
350  free_mem( indexymin );
351  free_mem( indexymax );
352 }
353 
354 //--------------------------------------------------------------------------
355 // void plsurf3dl(x, y, z, nx, ny, opt, clevel, nlevel,
356 // indexxmin, indexxmax, indexymin, indexymax)
357 //
358 // Plots the 3-d surface representation of the function z[x][y].
359 // The x values are stored as x[0..nx-1], the y values as y[0..ny-1],
360 // and the z values are in the 2-d array z[][].
361 //
362 //
363 // BASE_CONT draw contour plot at bottom xy plane
364 // TOP_CONT draw contour plot at top xy plane (not implemented)
365 // SURF_CONT draw contour at surface
366 // FACETED each square that makes up the surface is faceted
367 // DRAW_SIDES draw sides
368 // MAG_COLOR the surface is colored according to the value of z;
369 // if not defined, the surface is colored according to the
370 // intensity of the reflected light in the surface from a light
371 // source whose position is set using c_pllightsource()
372 //
373 // or any bitwise combination, e.g. "MAG_COLOR | SURF_CONT | SURF_BASE"
374 //
375 // indexymin and indexymax are arrays which specify the y index range
376 // (following the convention that the upper range limit is one more than
377 // actual index limit) for an x index range of indexxmin, indexxmax.
378 // This code is a complete departure from the approach taken in the old version
379 // of this routine. Formerly to code attempted to use the logic for the hidden
380 // line algorithm to draw the hidden surface. This was really hard. This code
381 // below uses a simple back to front (painters) algorithm. All the
382 // triangles are drawn.
383 //
384 // There are multitude of ways this code could be optimized. Given the
385 // problems with the old code, I tried to focus on clarity here.
386 //--------------------------------------------------------------------------
387 
388 void
389 c_plsurf3dl( const PLFLT *x, const PLFLT *y, const PLFLT * const *z, PLINT nx, PLINT ny,
390  PLINT opt, const PLFLT *clevel, PLINT nlevel,
391  PLINT indexxmin, PLINT indexxmax, const PLINT *indexymin, const PLINT *indexymax )
392 {
393  plfsurf3dl( x, y, plf2ops_c(), (PLPointer) z, nx, ny,
394  opt, clevel, nlevel, indexxmin, indexxmax, indexymin, indexymax );
395 }
396 
397 void
398 plfsurf3dl( const PLFLT *x, const PLFLT *y, PLF2OPS zops, PLPointer zp, PLINT nx, PLINT ny,
399  PLINT opt, const PLFLT *clevel, PLINT nlevel,
400  PLINT indexxmin, PLINT indexxmax, const PLINT *indexymin, const PLINT *indexymax )
401 {
402  PLFLT cxx, cxy, cyx, cyy, cyz;
403  PLINT i, j, k;
404  PLINT ixDir, ixOrigin, iyDir, iyOrigin, nFast, nSlow;
405  PLINT ixFast, ixSlow, iyFast, iySlow;
406  PLINT iFast, iSlow;
407  PLFLT xmin, xmax, ymin, ymax, zmin, zmax, zscale;
408  PLFLT xm, ym, zm;
409  PLINT ixmin = 0, ixmax = nx, iymin = 0, iymax = ny;
410  PLFLT xx[3], yy[3], zz[3];
411  PLFLT px[4], py[4], pz[4];
412  CONT_LEVEL *cont, *clev;
413  CONT_LINE *cline;
414  int ct, ix, iy, iftriangle;
415  PLINT color = plsc->icol0;
416  PLFLT width = plsc->width;
417  PLFLT ( *getz )( PLPointer, PLINT, PLINT ) = zops->get;
418 
419  if ( plsc->level < 3 )
420  {
421  myabort( "plsurf3dl: Please set up window first" );
422  return;
423  }
424 
425  if ( nx <= 0 || ny <= 0 )
426  {
427  myabort( "plsurf3dl: Bad array dimensions." );
428  return;
429  }
430 
431  //
432  // Don't use the data z value to scale the color, use the z axis
433  // values set by plw3d()
434  //
435  // plMinMax2dGrid(z, nx, ny, &fc_maxz, &fc_minz);
436  //
437 
438  fc_minz = plsc->ranmi;
439  fc_maxz = plsc->ranma;
440  if ( fc_maxz == fc_minz )
441  {
442  plwarn( "plsurf3dl: Maximum and minimum Z values are equal! \"fixing\"..." );
443  fc_maxz = fc_minz + 1e-6;
444  }
445 
446  if ( opt & MAG_COLOR )
447  falsecolor = 1;
448  else
449  falsecolor = 0;
450 
451  plP_gdom( &xmin, &xmax, &ymin, &ymax );
452  plP_grange( &zscale, &zmin, &zmax );
453  if ( zmin > zmax )
454  {
455  PLFLT t = zmin;
456  zmin = zmax;
457  zmax = t;
458  }
459 
460  // Check that points in x and in y are strictly increasing and in range
461 
462  for ( i = 0; i < nx - 1; i++ )
463  {
464  if ( x[i] >= x[i + 1] )
465  {
466  myabort( "plsurf3dl: X array must be strictly increasing" );
467  return;
468  }
469  if ( x[i] < xmin && x[i + 1] >= xmin )
470  ixmin = i;
471  if ( x[i + 1] > xmax && x[i] <= xmax )
472  ixmax = i + 2;
473  }
474  for ( i = 0; i < ny - 1; i++ )
475  {
476  if ( y[i] >= y[i + 1] )
477  {
478  myabort( "plsurf3dl: Y array must be strictly increasing" );
479  return;
480  }
481  if ( y[i] < ymin && y[i + 1] >= ymin )
482  iymin = i;
483  if ( y[i + 1] > ymax && y[i] <= ymax )
484  iymax = i + 2;
485  }
486 
487  // get the viewing parameters
488  plP_gw3wc( &cxx, &cxy, &cyx, &cyy, &cyz );
489 
490  // we're going to draw from back to front
491 
492  // iFast will index the dominant (fastest changing) dimension
493  // iSlow will index the slower changing dimension
494  //
495  // iX indexes the X dimension
496  // iY indexes the Y dimension
497 
498  // get direction for X
499  if ( cxy >= 0 )
500  {
501  ixDir = 1; // direction in X
502  ixOrigin = ixmin; // starting point
503  }
504  else
505  {
506  ixDir = -1;
507  ixOrigin = ixmax - 1;
508  }
509  // get direction for Y
510  if ( cxx >= 0 )
511  {
512  iyDir = -1;
513  iyOrigin = iymax - 1;
514  }
515  else
516  {
517  iyDir = 1;
518  iyOrigin = iymin;
519  }
520  // figure out which dimension is dominant
521  if ( fabs( cxx ) > fabs( cxy ) )
522  {
523  // X is dominant
524  nFast = ixmax - ixmin; // samples in the Fast direction
525  nSlow = iymax - iymin; // samples in the Slow direction
526 
527  ixFast = ixDir; ixSlow = 0;
528  iyFast = 0; iySlow = iyDir;
529  }
530  else
531  {
532  nFast = iymax - iymin;
533  nSlow = ixmax - ixmin;
534 
535  ixFast = 0; ixSlow = ixDir;
536  iyFast = iyDir; iySlow = 0;
537  }
538 
539  // we've got to draw the background grid first, hidden line code has to draw it last
540  if ( zbflg )
541  {
542  PLFLT bx[3], by[3], bz[3];
543  PLFLT tick = zbtck, tp;
544  PLINT nsub = 0;
545 
546  // get the tick spacing
547  pldtik( zmin, zmax, &tick, &nsub, FALSE );
548 
549  // determine the vertices for the background grid line
550  bx[0] = ( ixOrigin != ixmin && ixFast == 0 ) || ixFast > 0 ? xmax : xmin;
551  by[0] = ( iyOrigin != iymin && iyFast == 0 ) || iyFast > 0 ? ymax : ymin;
552  bx[1] = ixOrigin != ixmin ? xmax : xmin;
553  by[1] = iyOrigin != iymin ? ymax : ymin;
554  bx[2] = ( ixOrigin != ixmin && ixSlow == 0 ) || ixSlow > 0 ? xmax : xmin;
555  by[2] = ( iyOrigin != iymin && iySlow == 0 ) || iySlow > 0 ? ymax : ymin;
556 
557  plwidth( zbwidth );
558  plcol0( zbcol );
559  for ( tp = tick * floor( zmin / tick ) + tick; tp <= zmax; tp += tick )
560  {
561  bz[0] = bz[1] = bz[2] = tp;
562  plline3( 3, bx, by, bz );
563  }
564  // draw the vertical line at the back corner
565  bx[0] = bx[1];
566  by[0] = by[1];
567  bz[0] = zmin;
568  plline3( 2, bx, by, bz );
569  plwidth( width );
570  plcol0( color );
571  }
572 
573  // If enabled, draw the contour at the base
574 
575  // The contour plotted at the base will be identical to the one obtained
576  // with c_plcont(). The contour plotted at the surface is simple minded, but
577  // can be improved by using the contour data available.
578  //
579 
580  if ( clevel != NULL && opt & BASE_CONT )
581  {
582 #define NPTS 100
583  int np = NPTS;
584  PLFLT **zstore;
585  PLcGrid2 cgrid2;
586  PLFLT *zzloc = (PLFLT *) malloc( (size_t) NPTS * sizeof ( PLFLT ) );
587  if ( zzloc == NULL )
588  plexit( "plsurf3dl: Insufficient memory" );
589 
590  // get the contour lines
591 
592  // prepare cont_store input
593  cgrid2.nx = nx;
594  cgrid2.ny = ny;
595  plAlloc2dGrid( &cgrid2.xg, nx, ny );
596  plAlloc2dGrid( &cgrid2.yg, nx, ny );
597  plAlloc2dGrid( &zstore, nx, ny );
598 
599  for ( i = indexxmin; i < indexxmax; i++ )
600  {
601  for ( j = 0; j < indexymin[i]; j++ )
602  {
603  cgrid2.xg[i][j] = x[i];
604  cgrid2.yg[i][j] = y[indexymin[i]];
605  zstore[i][j] = getz( zp, i, indexymin[i] );
606  }
607  for ( j = indexymin[i]; j < indexymax[i]; j++ )
608  {
609  cgrid2.xg[i][j] = x[i];
610  cgrid2.yg[i][j] = y[j];
611  zstore[i][j] = getz( zp, i, j );
612  }
613  for ( j = indexymax[i]; j < ny; j++ )
614  {
615  cgrid2.xg[i][j] = x[i];
616  cgrid2.yg[i][j] = y[indexymax[i] - 1];
617  zstore[i][j] = getz( zp, i, indexymax[i] - 1 );
618  }
619  }
620  // Fill cont structure with contours.
621  cont_store( (const PLFLT * const *) zstore, nx, ny, indexxmin + 1, indexxmax, 1, ny,
622  clevel, nlevel, pltr2, (void *) &cgrid2, &cont );
623 
624  // Free the 2D input arrays to cont_store since not needed any more.
625  plFree2dGrid( zstore, nx, ny );
626  plFree2dGrid( cgrid2.xg, nx, ny );
627  plFree2dGrid( cgrid2.yg, nx, ny );
628 
629  // follow the contour levels and lines
630  clev = cont;
631  do // for each contour level
632  {
633  cline = clev->line;
634  do // there are several lines that make up the contour
635  {
636  if ( cline->npts > np )
637  {
638  np = cline->npts;
639  if ( ( zzloc = (PLFLT *) realloc( zzloc, (size_t) np * sizeof ( PLFLT ) ) ) == NULL )
640  {
641  plexit( "plsurf3dl: Insufficient memory" );
642  }
643  }
644  for ( j = 0; j < cline->npts; j++ )
645  zzloc[j] = plsc->ranmi;
646  if ( cline->npts > 0 )
647  {
648  plcol1( ( clev->level - fc_minz ) / ( fc_maxz - fc_minz ) );
649  plline3( cline->npts, cline->x, cline->y, zzloc );
650  }
651  cline = cline->next;
652  }
653  while ( cline != NULL );
654  clev = clev->next;
655  }
656  while ( clev != NULL );
657 
658  cont_clean_store( cont ); // now release the memory
659  free( zzloc );
660  }
661 
662  // Now we can iterate over the grid drawing the quads
663  for ( iSlow = 0; iSlow < nSlow - 1; iSlow++ )
664  {
665  for ( iFast = 0; iFast < nFast - 1; iFast++ )
666  {
667  // get the 4 corners of the Quad, which are
668  //
669  // 0--2
670  // | |
671  // 1--3
672  //
673 
674  xm = ym = zm = 0.;
675 
676  iftriangle = 1;
677  for ( i = 0; i < 2; i++ )
678  {
679  for ( j = 0; j < 2; j++ )
680  {
681  // we're transforming from Fast/Slow coordinates to x/y coordinates
682  // note, these are the indices, not the values
683  ix = ixFast * ( iFast + i ) + ixSlow * ( iSlow + j ) + ixOrigin;
684  iy = iyFast * ( iFast + i ) + iySlow * ( iSlow + j ) + iyOrigin;
685 
686  if ( indexxmin <= ix && ix < indexxmax &&
687  indexymin[ix] <= iy && iy < indexymax[ix] )
688  {
689  xm += px[2 * i + j] = x[ix];
690  ym += py[2 * i + j] = y[iy];
691  zm += pz[2 * i + j] = getz( zp, ix, iy );
692  }
693  else
694  {
695  iftriangle = 0;
696  break;
697  }
698  }
699  if ( iftriangle == 0 )
700  break;
701  }
702 
703  if ( iftriangle == 0 )
704  continue;
705  // the "mean point" of the quad, common to all four triangles
706  // -- perhaps not a good thing to do for the light shading
707 
708  xm /= 4.; ym /= 4.; zm /= 4.;
709 
710  // now draw the quad as four triangles
711 
712  for ( i = 1; i < 3; i++ )
713  {
714  for ( j = 0; j < 4; j += 3 )
715  {
716  shade_triangle( px[j], py[j], pz[j], xm, ym, zm, px[i], py[i], pz[i] );
717 
718  // after shading, see if the triangle crosses one contour plane
719 
720 #define min3( a, b, c ) ( MIN( ( MIN( a, b ) ), c ) )
721 #define max3( a, b, c ) ( MAX( ( MAX( a, b ) ), c ) )
722 
723  if ( clevel != NULL && ( opt & SURF_CONT ) )
724  {
725  for ( k = 0; k < nlevel; k++ )
726  {
727  if ( clevel[k] >= min3( pz[i], zm, pz[j] ) && clevel[k] < max3( pz[i], zm, pz[j] ) )
728  {
729  ct = 0;
730  if ( clevel[k] >= MIN( pz[i], zm ) && clevel[k] < MAX( pz[i], zm ) ) // p0-pm
731  {
732  xx[ct] = ( ( clevel[k] - pz[i] ) * ( xm - px[i] ) ) / ( zm - pz[i] ) + px[i];
733  yy[ct] = ( ( clevel[k] - pz[i] ) * ( ym - py[i] ) ) / ( zm - pz[i] ) + py[i];
734  ct++;
735  }
736 
737  if ( clevel[k] >= MIN( pz[i], pz[j] ) && clevel[k] < MAX( pz[i], pz[j] ) ) // p0-p1
738  {
739  xx[ct] = ( ( clevel[k] - pz[i] ) * ( px[j] - px[i] ) ) / ( pz[j] - pz[i] ) + px[i];
740  yy[ct] = ( ( clevel[k] - pz[i] ) * ( py[j] - py[i] ) ) / ( pz[j] - pz[i] ) + py[i];
741  ct++;
742  }
743 
744  if ( clevel[k] >= MIN( pz[j], zm ) && clevel[k] < MAX( pz[j], zm ) ) // p1-pm
745  {
746  xx[ct] = ( ( clevel[k] - pz[j] ) * ( xm - px[j] ) ) / ( zm - pz[j] ) + px[j];
747  yy[ct] = ( ( clevel[k] - pz[j] ) * ( ym - py[j] ) ) / ( zm - pz[j] ) + py[j];
748  ct++;
749  }
750 
751  if ( ct == 2 )
752  {
753  // yes, xx and yy are the intersection points of the triangle with
754  // the contour line -- draw a straight line betweeen the points
755  // -- at the end this will make up the contour line
756 
757  if ( opt & SURF_CONT )
758  {
759  // surface contour with color set by user
760  plcol0( color );
761  zz[0] = zz[1] = clevel[k];
762  plline3( 2, xx, yy, zz );
763  }
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