PLplot  5.10.0
plot3d.c
Go to the documentation of this file.
00001 
00002 
00003 
00004 
00005 // Copyright (C) 2004-2014 Alan W. Irwin
00006 // Copyright (C) 2004  Joao Cardoso
00007 // Copyright (C) 2004  Andrew Ross
00008 //
00009 // This file is part of PLplot.
00010 //
00011 // PLplot is free software; you can redistribute it and/or modify
00012 // it under the terms of the GNU Library General Public License as published
00013 // by the Free Software Foundation; either version 2 of the License, or
00014 // (at your option) any later version.
00015 //
00016 // PLplot is distributed in the hope that it will be useful,
00017 // but WITHOUT ANY WARRANTY; without even the implied warranty of
00018 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00019 // GNU Library General Public License for more details.
00020 //
00021 // You should have received a copy of the GNU Library General Public License
00022 // along with PLplot; if not, write to the Free Software
00023 // Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
00024 //
00025 
00026 #include "plplotP.h"
00027 
00028 // Internal constants
00029 
00030 #define  BINC    50             // Block size for memory allocation
00031 
00032 static PLINT pl3mode = 0;       // 0 3d solid; 1 mesh plot
00033 static PLINT pl3upv  = 1;       // 1 update view; 0 no update
00034 
00035 static PLINT zbflg = 0, zbcol;
00036 static PLFLT zbtck, zbwidth;
00037 
00038 static PLINT *oldhiview = NULL;
00039 static PLINT *oldloview = NULL;
00040 static PLINT *newhiview = NULL;
00041 static PLINT *newloview = NULL;
00042 static PLINT *utmp      = NULL;
00043 static PLINT *vtmp      = NULL;
00044 static PLFLT *ctmp      = NULL;
00045 
00046 static PLINT mhi, xxhi, newhisize;
00047 static PLINT mlo, xxlo, newlosize;
00048 
00049 // Light source for shading
00050 static PLFLT xlight, ylight, zlight;
00051 static PLINT falsecolor = 0;
00052 static PLFLT fc_minz, fc_maxz;
00053 
00054 // Prototypes for static functions
00055 
00056 static void plgrid3( PLFLT );
00057 static void plnxtv( PLINT *, PLINT *, PLFLT*, PLINT, PLINT );
00058 static void plside3( const PLFLT *, const PLFLT *, PLF2OPS, PLPointer, PLINT, PLINT, PLINT );
00059 static void plt3zz( PLINT, PLINT, PLINT, PLINT,
00060                     PLINT, PLINT *, const PLFLT *, const PLFLT *, PLF2OPS, PLPointer,
00061                     PLINT, PLINT, PLINT *, PLINT *, PLFLT* );
00062 static void plnxtvhi( PLINT *, PLINT *, PLFLT*, PLINT, PLINT );
00063 static void plnxtvlo( PLINT *, PLINT *, PLFLT*, PLINT, PLINT );
00064 static void plnxtvhi_draw( PLINT *u, PLINT *v, PLFLT* c, PLINT n );
00065 
00066 static void savehipoint( PLINT, PLINT );
00067 static void savelopoint( PLINT, PLINT );
00068 static void swaphiview( void );
00069 static void swaploview( void );
00070 static void myexit( const char * );
00071 static void myabort( const char * );
00072 static void freework( void );
00073 static int  plabv( PLINT, PLINT, PLINT, PLINT, PLINT, PLINT );
00074 static void pl3cut( PLINT, PLINT, PLINT, PLINT, PLINT,
00075                     PLINT, PLINT, PLINT, PLINT *, PLINT * );
00076 static PLFLT plGetAngleToLight( PLFLT* x, PLFLT* y, PLFLT* z );
00077 static void plP_draw3d( PLINT x, PLINT y, PLFLT *c, PLINT j, PLINT move );
00078 //static void plxyindexlimits( PLINT instart, PLINT inn,
00079 //                             PLINT *inarray_min, PLINT *inarray_max,
00080 //                             PLINT *outstart, PLINT *outn, PLINT outnmax,
00081 //                             PLINT *outarray_min, PLINT *outarray_max );
00082 
00083 
00084 // #define MJL_HACK 1
00085 #if MJL_HACK
00086 static void plP_fill3( PLINT x0, PLINT y0, PLINT x1, PLINT y1,
00087                        PLINT x2, PLINT y2, PLINT j );
00088 static void plP_fill4( PLINT x0, PLINT y0, PLINT x1, PLINT y1,
00089                        PLINT x2, PLINT y2, PLINT x3, PLINT y3, PLINT j );
00090 #endif
00091 
00092 //--------------------------------------------------------------------------
00093 // void plsetlightsource(x, y, z)
00094 //
00095 // Sets the position of the light source.
00096 //--------------------------------------------------------------------------
00097 
00098 void
00099 c_pllightsource( PLFLT x, PLFLT y, PLFLT z )
00100 {
00101     xlight = x;
00102     ylight = y;
00103     zlight = z;
00104 }
00105 
00106 //--------------------------------------------------------------------------
00107 // void plmesh(x, y, z, nx, ny, opt)
00108 //
00109 // Plots a mesh representation of the function z[x][y]. The x values
00110 // are stored as x[0..nx-1], the y values as y[0..ny-1], and the
00111 // z values are in the 2-d array z[][]. The integer "opt" specifies:
00112 // see plmeshc() below.
00113 //--------------------------------------------------------------------------
00114 
00115 void
00116 c_plmesh( const PLFLT *x, const PLFLT *y, const PLFLT * const *z, PLINT nx, PLINT ny, PLINT opt )
00117 {
00118     plfplot3dc( x, y, plf2ops_c(), (PLPointer) z, nx, ny, opt | MESH, NULL, 0 );
00119 }
00120 
00121 void
00122 plfmesh( const PLFLT *x, const PLFLT *y, PLF2OPS zops, PLPointer zp,
00123          PLINT nx, PLINT ny, PLINT opt )
00124 {
00125     plfplot3dc( x, y, zops, zp, nx, ny, opt | MESH, NULL, 0 );
00126 }
00127 
00128 //--------------------------------------------------------------------------
00129 // void plmeshc(x, y, z, nx, ny, opt, clevel, nlevel)
00130 //
00131 // Plots a mesh representation of the function z[x][y]. The x values
00132 // are stored as x[0..nx-1], the y values as y[0..ny-1], and the
00133 // z values are in the 2-d array z[][]. The integer "opt" specifies:
00134 //
00135 // DRAW_LINEX   draw lines parallel to the X axis
00136 // DRAW_LINEY   draw lines parallel to the Y axis
00137 // DRAW_LINEXY  draw lines parallel to both the X and Y axis
00138 // MAG_COLOR    draw the mesh with a color dependent of the magnitude
00139 // BASE_CONT    draw contour plot at bottom xy plane
00140 // TOP_CONT     draw contour plot at top xy plane (not yet)
00141 // DRAW_SIDES   draw sides
00142 //
00143 // or any bitwise combination, e.g. "MAG_COLOR | DRAW_LINEX"
00144 //
00145 //--------------------------------------------------------------------------
00146 
00147 void
00148 c_plmeshc( const PLFLT *x, const PLFLT *y, const PLFLT * const *z, PLINT nx, PLINT ny, PLINT opt,
00149            const PLFLT *clevel, PLINT nlevel )
00150 {
00151     plfplot3dc( x, y, plf2ops_c(), (PLPointer) z, nx, ny, opt | MESH, clevel, nlevel );
00152 }
00153 
00154 void
00155 plfmeshc( const PLFLT *x, const PLFLT *y, PLF2OPS zops, PLPointer zp,
00156           PLINT nx, PLINT ny, PLINT opt, const PLFLT *clevel, PLINT nlevel )
00157 {
00158     plfplot3dc( x, y, zops, zp, nx, ny, opt | MESH, clevel, nlevel );
00159 }
00160 
00161 // clipping helper for 3d polygons
00162 
00163 int
00164 plP_clip_poly( int Ni, PLFLT *Vi[3], int axis, PLFLT dir, PLFLT offset )
00165 {
00166     int   anyout = 0;
00167     PLFLT _in[PL_MAXPOLY], _T[3][PL_MAXPOLY];
00168     PLFLT *in, *T[3], *TT = NULL;
00169     int   No = 0;
00170     int   i, j, k;
00171 
00172     if ( Ni > PL_MAXPOLY )
00173     {
00174         in = (PLFLT *) malloc( sizeof ( PLFLT ) * (size_t) Ni );
00175         TT = (PLFLT *) malloc( 3 * sizeof ( PLFLT ) * (size_t) Ni );
00176 
00177         if ( in == NULL || TT == NULL )
00178         {
00179             plexit( "plP_clip_poly: insufficient memory for large polygon" );
00180         }
00181 
00182         T[0] = &TT[0];
00183         T[1] = &TT[Ni];
00184         T[2] = &TT[2 * Ni];
00185     }
00186     else
00187     {
00188         in   = _in;
00189         T[0] = &_T[0][0];
00190         T[1] = &_T[1][0];
00191         T[2] = &_T[2][0];
00192     }
00193 
00194     for ( i = 0; i < Ni; i++ )
00195     {
00196         in[i]   = Vi[axis][i] * dir + offset;
00197         anyout += in[i] < 0;
00198     }
00199 
00200     // none out
00201     if ( anyout == 0 )
00202         return Ni;
00203 
00204     // all out
00205     if ( anyout == Ni )
00206     {
00207         return 0;
00208     }
00209 
00210     // some out
00211     // copy over to a temp vector
00212     for ( i = 0; i < 3; i++ )
00213     {
00214         for ( j = 0; j < Ni; j++ )
00215         {
00216             T[i][j] = Vi[i][j];
00217         }
00218     }
00219     // copy back selectively
00220     for ( i = 0; i < Ni; i++ )
00221     {
00222         j = ( i + 1 ) % Ni;
00223 
00224         if ( in[i] >= 0 && in[j] >= 0 )
00225         {
00226             for ( k = 0; k < 3; k++ )
00227                 Vi[k][No] = T[k][j];
00228             No++;
00229         }
00230         else if ( in[i] >= 0 && in[j] < 0 )
00231         {
00232             PLFLT u = in[i] / ( in[i] - in[j] );
00233             for ( k = 0; k < 3; k++ )
00234                 Vi[k][No] = T[k][i] * ( 1 - u ) + T[k][j] * u;
00235             No++;
00236         }
00237         else if ( in[i] < 0 && in[j] >= 0 )
00238         {
00239             PLFLT u = in[i] / ( in[i] - in[j] );
00240             for ( k = 0; k < 3; k++ )
00241                 Vi[k][No] = T[k][i] * ( 1 - u ) + T[k][j] * u;
00242             No++;
00243             for ( k = 0; k < 3; k++ )
00244                 Vi[k][No] = T[k][j];
00245             No++;
00246         }
00247     }
00248 
00249     if ( Ni > PL_MAXPOLY )
00250     {
00251         free( in );
00252         free( TT );
00253     }
00254 
00255     return No;
00256 }
00257 
00258 // helper for plsurf3d, similar to c_plfill3()
00259 static void
00260 shade_triangle( PLFLT x0, PLFLT y0, PLFLT z0,
00261                 PLFLT x1, PLFLT y1, PLFLT z1,
00262                 PLFLT x2, PLFLT y2, PLFLT z2 )
00263 {
00264     int   i;
00265     // arrays for interface to core functions
00266     short u[6], v[6];
00267     PLFLT x[6], y[6], z[6];
00268     int   n;
00269     PLFLT xmin, xmax, ymin, ymax, zmin, zmax, zscale;
00270     PLFLT *V[3];
00271 
00272     plP_gdom( &xmin, &xmax, &ymin, &ymax );
00273     plP_grange( &zscale, &zmin, &zmax );
00274 
00275     x[0] = x0; x[1] = x1; x[2] = x2;
00276     y[0] = y0; y[1] = y1; y[2] = y2;
00277     z[0] = z0; z[1] = z1; z[2] = z2;
00278     n    = 3;
00279 
00280     V[0] = x; V[1] = y; V[2] = z;
00281 
00282     n = plP_clip_poly( n, V, 0, 1, -xmin );
00283     n = plP_clip_poly( n, V, 0, -1, xmax );
00284     n = plP_clip_poly( n, V, 1, 1, -ymin );
00285     n = plP_clip_poly( n, V, 1, -1, ymax );
00286     n = plP_clip_poly( n, V, 2, 1, -zmin );
00287     n = plP_clip_poly( n, V, 2, -1, zmax );
00288 
00289     if ( n > 0 )
00290     {
00291         if ( falsecolor )
00292             plcol1( ( ( z[0] + z[1] + z[2] ) / 3. - fc_minz ) / ( fc_maxz - fc_minz ) );
00293         else
00294             plcol1( plGetAngleToLight( x, y, z ) );
00295 
00296         for ( i = 0; i < n; i++ )
00297         {
00298             u[i] = (short) plP_wcpcx( plP_w3wcx( x[i], y[i], z[i] ) );
00299             v[i] = (short) plP_wcpcy( plP_w3wcy( x[i], y[i], z[i] ) );
00300         }
00301         u[n] = u[0];
00302         v[n] = v[0];
00303 
00304 #ifdef SHADE_DEBUG // show triangles
00305         plcol0( 1 );
00306         x[3] = x[0]; y[3] = y[0]; z[3] = z[0];
00307         plline3( 4, x, y, z );
00308 #else   // fill triangles
00309         plP_fill( u, v, n + 1 );
00310 #endif
00311     }
00312 }
00313 
00314 //--------------------------------------------------------------------------
00315 // void plsurf3d(x, y, z, nx, ny, opt, clevel, nlevel)
00316 //
00317 // Plots the 3-d surface representation of the function z[x][y].
00318 // The x values are stored as x[0..nx-1], the y values as y[0..ny-1],
00319 //  and the z values are in the 2-d array z[][].  The integer "opt" specifies:
00320 // see plsurf3dl() below.
00321 //--------------------------------------------------------------------------
00322 
00323 void
00324 c_plsurf3d( const PLFLT *x, const PLFLT *y, const PLFLT * const *z, PLINT nx, PLINT ny,
00325             PLINT opt, const PLFLT *clevel, PLINT nlevel )
00326 {
00327     plfsurf3d( x, y, plf2ops_c(), (PLPointer) z, nx, ny,
00328         opt, clevel, nlevel );
00329 }
00330 
00331 void
00332 plfsurf3d( const PLFLT *x, const PLFLT *y, PLF2OPS zops, PLPointer zp,
00333            PLINT nx, PLINT ny, PLINT opt, const PLFLT *clevel, PLINT nlevel )
00334 {
00335     PLINT i;
00336     PLINT *indexymin = (PLINT *) malloc( (size_t) nx * sizeof ( PLINT ) );
00337     PLINT *indexymax = (PLINT *) malloc( (size_t) nx * sizeof ( PLINT ) );
00338 
00339     if ( !indexymin || !indexymax )
00340         plexit( "plsurf3d: Out of memory." );
00341     for ( i = 0; i < nx; i++ )
00342     {
00343         indexymin[i] = 0;
00344         indexymax[i] = ny;
00345     }
00346     plfsurf3dl( x, y, zops, zp, nx, ny, opt, clevel, nlevel,
00347         0, nx, indexymin, indexymax );
00348     free_mem( indexymin );
00349     free_mem( indexymax );
00350 }
00351 
00352 //--------------------------------------------------------------------------
00353 // void plsurf3dl(x, y, z, nx, ny, opt, clevel, nlevel,
00354 // indexxmin, indexxmax, indexymin, indexymax)
00355 //
00356 // Plots the 3-d surface representation of the function z[x][y].
00357 // The x values are stored as x[0..nx-1], the y values as y[0..ny-1],
00358 //  and the z values are in the 2-d array z[][].
00359 //
00360 //
00361 // BASE_CONT    draw contour plot at bottom xy plane
00362 // TOP_CONT     draw contour plot at top xy plane (not implemented)
00363 // SURF_CONT    draw contour at surface
00364 // FACETED      each square that makes up the surface is faceted
00365 // DRAW_SIDES   draw sides
00366 // MAG_COLOR    the surface is colored according to the value of z;
00367 //               if not defined, the surface is colored according to the
00368 //               intensity of the reflected light in the surface from a light
00369 //               source whose position is set using c_pllightsource()
00370 //
00371 // or any bitwise combination, e.g. "MAG_COLOR | SURF_CONT | SURF_BASE"
00372 //
00373 // indexymin and indexymax are arrays which specify the y index range
00374 // (following the convention that the upper range limit is one more than
00375 // actual index limit) for an x index range of indexxmin, indexxmax.
00376 // This code is a complete departure from the approach taken in the old version
00377 // of this routine. Formerly to code attempted to use the logic for the hidden
00378 // line algorithm to draw the hidden surface. This was really hard. This code
00379 // below uses a simple back to front (painters) algorithm. All the
00380 // triangles are drawn.
00381 //
00382 // There are multitude of ways this code could be optimized. Given the
00383 // problems with the old code, I tried to focus on clarity here.
00384 //--------------------------------------------------------------------------
00385 
00386 void
00387 c_plsurf3dl( const PLFLT *x, const PLFLT *y, const PLFLT * const *z, PLINT nx, PLINT ny,
00388              PLINT opt, const PLFLT *clevel, PLINT nlevel,
00389              PLINT indexxmin, PLINT indexxmax, const PLINT *indexymin, const PLINT *indexymax )
00390 {
00391     plfsurf3dl( x, y, plf2ops_c(), (PLPointer) z, nx, ny,
00392         opt, clevel, nlevel, indexxmin, indexxmax, indexymin, indexymax );
00393 }
00394 
00395 void
00396 plfsurf3dl( const PLFLT *x, const PLFLT *y, PLF2OPS zops, PLPointer zp, PLINT nx, PLINT ny,
00397             PLINT opt, const PLFLT *clevel, PLINT nlevel,
00398             PLINT indexxmin, PLINT indexxmax, const PLINT *indexymin, const PLINT *indexymax )
00399 {
00400     PLFLT      cxx, cxy, cyx, cyy, cyz;
00401     PLINT      i, j, k;
00402     PLINT      ixDir, ixOrigin, iyDir, iyOrigin, nFast, nSlow;
00403     PLINT      ixFast, ixSlow, iyFast, iySlow;
00404     PLINT      iFast, iSlow;
00405     PLFLT      xmin, xmax, ymin, ymax, zmin, zmax, zscale;
00406     PLFLT      xm, ym, zm;
00407     PLINT      ixmin = 0, ixmax = nx, iymin = 0, iymax = ny;
00408     PLFLT      xx[3], yy[3], zz[3];
00409     PLFLT      px[4], py[4], pz[4];
00410     CONT_LEVEL *cont, *clev;
00411     CONT_LINE  *cline;
00412     int        ct, ix, iy, iftriangle;
00413     PLINT      color = plsc->icol0;
00414     PLFLT      width = plsc->width;
00415     PLFLT      ( *getz )( PLPointer, PLINT, PLINT ) = zops->get;
00416 
00417     if ( plsc->level < 3 )
00418     {
00419         myabort( "plsurf3dl: Please set up window first" );
00420         return;
00421     }
00422 
00423     if ( nx <= 0 || ny <= 0 )
00424     {
00425         myabort( "plsurf3dl: Bad array dimensions." );
00426         return;
00427     }
00428 
00429     //
00430     // Don't use the data z value to scale the color, use the z axis
00431     // values set by plw3d()
00432     //
00433     // plMinMax2dGrid(z, nx, ny, &fc_maxz, &fc_minz);
00434     //
00435 
00436     fc_minz = plsc->ranmi;
00437     fc_maxz = plsc->ranma;
00438     if ( fc_maxz == fc_minz )
00439     {
00440         plwarn( "plsurf3dl: Maximum and minimum Z values are equal! \"fixing\"..." );
00441         fc_maxz = fc_minz + 1e-6;
00442     }
00443 
00444     if ( opt & MAG_COLOR )
00445         falsecolor = 1;
00446     else
00447         falsecolor = 0;
00448 
00449     plP_gdom( &xmin, &xmax, &ymin, &ymax );
00450     plP_grange( &zscale, &zmin, &zmax );
00451     if ( zmin > zmax )
00452     {
00453         PLFLT t = zmin;
00454         zmin = zmax;
00455         zmax = t;
00456     }
00457 
00458     // Check that points in x and in y are strictly increasing  and in range
00459 
00460     for ( i = 0; i < nx - 1; i++ )
00461     {
00462         if ( x[i] >= x[i + 1] )
00463         {
00464             myabort( "plsurf3dl: X array must be strictly increasing" );
00465             return;
00466         }
00467         if ( x[i] < xmin && x[i + 1] >= xmin )
00468             ixmin = i;
00469         if ( x[i + 1] > xmax && x[i] <= xmax )
00470             ixmax = i + 2;
00471     }
00472     for ( i = 0; i < ny - 1; i++ )
00473     {
00474         if ( y[i] >= y[i + 1] )
00475         {
00476             myabort( "plsurf3dl: Y array must be strictly increasing" );
00477             return;
00478         }
00479         if ( y[i] < ymin && y[i + 1] >= ymin )
00480             iymin = i;
00481         if ( y[i + 1] > ymax && y[i] <= ymax )
00482             iymax = i + 2;
00483     }
00484 
00485     // get the viewing parameters
00486     plP_gw3wc( &cxx, &cxy, &cyx, &cyy, &cyz );
00487 
00488     // we're going to draw from back to front
00489 
00490     // iFast will index the dominant (fastest changing) dimension
00491     // iSlow will index the slower changing dimension
00492     //
00493     // iX indexes the X dimension
00494     // iY indexes the Y dimension
00495 
00496     // get direction for X
00497     if ( cxy >= 0 )
00498     {
00499         ixDir    = 1;     // direction in X
00500         ixOrigin = ixmin; // starting point
00501     }
00502     else
00503     {
00504         ixDir    = -1;
00505         ixOrigin = ixmax - 1;
00506     }
00507     // get direction for Y
00508     if ( cxx >= 0 )
00509     {
00510         iyDir    = -1;
00511         iyOrigin = iymax - 1;
00512     }
00513     else
00514     {
00515         iyDir    = 1;
00516         iyOrigin = iymin;
00517     }
00518     // figure out which dimension is dominant
00519     if ( fabs( cxx ) > fabs( cxy ) )
00520     {
00521         // X is dominant
00522         nFast = ixmax - ixmin;  // samples in the Fast direction
00523         nSlow = iymax - iymin;  // samples in the Slow direction
00524 
00525         ixFast = ixDir; ixSlow = 0;
00526         iyFast = 0;     iySlow = iyDir;
00527     }
00528     else
00529     {
00530         nFast = iymax - iymin;
00531         nSlow = ixmax - ixmin;
00532 
00533         ixFast = 0;     ixSlow = ixDir;
00534         iyFast = iyDir; iySlow = 0;
00535     }
00536 
00537     // we've got to draw the background grid first, hidden line code has to draw it last
00538     if ( zbflg )
00539     {
00540         PLFLT bx[3], by[3], bz[3];
00541         PLFLT tick = zbtck, tp;
00542         PLINT nsub = 0;
00543 
00544         // get the tick spacing
00545         pldtik( zmin, zmax, &tick, &nsub, FALSE );
00546 
00547         // determine the vertices for the background grid line
00548         bx[0] = ( ixOrigin != ixmin && ixFast == 0 ) || ixFast > 0 ? xmax : xmin;
00549         by[0] = ( iyOrigin != iymin && iyFast == 0 ) || iyFast > 0 ? ymax : ymin;
00550         bx[1] = ixOrigin != ixmin ? xmax : xmin;
00551         by[1] = iyOrigin != iymin ? ymax : ymin;
00552         bx[2] = ( ixOrigin != ixmin && ixSlow == 0 ) || ixSlow > 0 ? xmax : xmin;
00553         by[2] = ( iyOrigin != iymin && iySlow == 0 ) || iySlow > 0 ? ymax : ymin;
00554 
00555         plwidth( zbwidth );
00556         plcol0( zbcol );
00557         for ( tp = tick * floor( zmin / tick ) + tick; tp <= zmax; tp += tick )
00558         {
00559             bz[0] = bz[1] = bz[2] = tp;
00560             plline3( 3, bx, by, bz );
00561         }
00562         // draw the vertical line at the back corner
00563         bx[0] = bx[1];
00564         by[0] = by[1];
00565         bz[0] = zmin;
00566         plline3( 2, bx, by, bz );
00567         plwidth( width );
00568         plcol0( color );
00569     }
00570 
00571     // If enabled, draw the contour at the base
00572 
00573     // The contour plotted at the base will be identical to the one obtained
00574     // with c_plcont(). The contour plotted at the surface is simple minded, but
00575     // can be improved by using the contour data available.
00576     //
00577 
00578     if ( clevel != NULL && opt & BASE_CONT )
00579     {
00580 #define NPTS    100
00581         int      np = NPTS;
00582         PLFLT    **zstore;
00583         PLcGrid2 cgrid2;
00584         PLFLT    *zzloc = (PLFLT *) malloc( (size_t) NPTS * sizeof ( PLFLT ) );
00585         if ( zzloc == NULL )
00586             plexit( "plsurf3dl: Insufficient memory" );
00587 
00588         // get the contour lines
00589 
00590         // prepare cont_store input
00591         cgrid2.nx = nx;
00592         cgrid2.ny = ny;
00593         plAlloc2dGrid( &cgrid2.xg, nx, ny );
00594         plAlloc2dGrid( &cgrid2.yg, nx, ny );
00595         plAlloc2dGrid( &zstore, nx, ny );
00596 
00597         for ( i = indexxmin; i < indexxmax; i++ )
00598         {
00599             for ( j = 0; j < indexymin[i]; j++ )
00600             {
00601                 cgrid2.xg[i][j] = x[i];
00602                 cgrid2.yg[i][j] = y[indexymin[i]];
00603                 zstore[i][j]    = getz( zp, i, indexymin[i] );
00604             }
00605             for ( j = indexymin[i]; j < indexymax[i]; j++ )
00606             {
00607                 cgrid2.xg[i][j] = x[i];
00608                 cgrid2.yg[i][j] = y[j];
00609                 zstore[i][j]    = getz( zp, i, j );
00610             }
00611             for ( j = indexymax[i]; j < ny; j++ )
00612             {
00613                 cgrid2.xg[i][j] = x[i];
00614                 cgrid2.yg[i][j] = y[indexymax[i] - 1];
00615                 zstore[i][j]    = getz( zp, i, indexymax[i] - 1 );
00616             }
00617         }
00618         // Fill cont structure with contours.
00619         cont_store( (const PLFLT * const *) zstore, nx, ny, indexxmin + 1, indexxmax, 1, ny,
00620             clevel, nlevel, pltr2, (void *) &cgrid2, &cont );
00621 
00622         // Free the 2D input arrays to cont_store since not needed any more.
00623         plFree2dGrid( zstore, nx, ny );
00624         plFree2dGrid( cgrid2.xg, nx, ny );
00625         plFree2dGrid( cgrid2.yg, nx, ny );
00626 
00627         // follow the contour levels and lines
00628         clev = cont;
00629         do  // for each contour level
00630         {
00631             cline = clev->line;
00632             do  // there are several lines that make up the contour
00633             {
00634                 if ( cline->npts > np )
00635                 {
00636                     np = cline->npts;
00637                     if ( ( zzloc = (PLFLT *) realloc( zzloc, (size_t) np * sizeof ( PLFLT ) ) ) == NULL )
00638                     {
00639                         plexit( "plsurf3dl: Insufficient memory" );
00640                     }
00641                 }
00642                 for ( j = 0; j < cline->npts; j++ )
00643                     zzloc[j] = plsc->ranmi;
00644                 if ( cline->npts > 0 )
00645                 {
00646                     plcol1( ( clev->level - fc_minz ) / ( fc_maxz - fc_minz ) );
00647                     plline3( cline->npts, cline->x, cline->y, zzloc );
00648                 }
00649                 cline = cline->next;
00650             }
00651             while ( cline != NULL );
00652             clev = clev->next;
00653         }
00654         while ( clev != NULL );
00655 
00656         cont_clean_store( cont ); // now release the memory
00657         free( zzloc );
00658     }
00659 
00660     // Now we can iterate over the grid drawing the quads
00661     for ( iSlow = 0; iSlow < nSlow - 1; iSlow++ )
00662     {
00663         for ( iFast = 0; iFast < nFast - 1; iFast++ )
00664         {
00665             // get the 4 corners of the Quad, which are
00666             //
00667             //       0--2
00668             //       |  |
00669             //       1--3
00670             //
00671 
00672             xm = ym = zm = 0.;
00673 
00674             iftriangle = 1;
00675             for ( i = 0; i < 2; i++ )
00676             {
00677                 for ( j = 0; j < 2; j++ )
00678                 {
00679                     // we're transforming from Fast/Slow coordinates to x/y coordinates
00680                     // note, these are the indices, not the values
00681                     ix = ixFast * ( iFast + i ) + ixSlow * ( iSlow + j ) + ixOrigin;
00682                     iy = iyFast * ( iFast + i ) + iySlow * ( iSlow + j ) + iyOrigin;
00683 
00684                     if ( indexxmin <= ix && ix < indexxmax &&
00685                          indexymin[ix] <= iy && iy < indexymax[ix] )
00686                     {
00687                         xm += px[2 * i + j] = x[ix];
00688                         ym += py[2 * i + j] = y[iy];
00689                         zm += pz[2 * i + j] = getz( zp, ix, iy );
00690                     }
00691                     else
00692                     {
00693                         iftriangle = 0;
00694                         break;
00695                     }
00696                 }
00697                 if ( iftriangle == 0 )
00698                     break;
00699             }
00700 
00701             if ( iftriangle == 0 )
00702                 continue;
00703             // the "mean point" of the quad, common to all four triangles
00704             // -- perhaps not a good thing to do for the light shading
00705 
00706             xm /= 4.; ym /= 4.; zm /= 4.;
00707 
00708             // now draw the quad as four triangles
00709 
00710             for ( i = 1; i < 3; i++ )
00711             {
00712                 for ( j = 0; j < 4; j += 3 )
00713                 {
00714                     shade_triangle( px[j], py[j], pz[j], xm, ym, zm, px[i], py[i], pz[i] );
00715 
00716                     // after shading, see if the triangle crosses       one contour plane
00717 
00718 #define min3( a, b, c )    ( MIN( ( MIN( a, b ) ), c ) )
00719 #define max3( a, b, c )    ( MAX( ( MAX( a, b ) ), c ) )
00720 
00721                     if ( clevel != NULL && ( opt & SURF_CONT ) )
00722                     {
00723                         for ( k = 0; k < nlevel; k++ )
00724                         {
00725                             if ( clevel[k] >= min3( pz[i], zm, pz[j] ) && clevel[k] < max3( pz[i], zm, pz[j] ) )
00726                             {
00727                                 ct = 0;
00728                                 if ( clevel[k] >= MIN( pz[i], zm ) && clevel[k] < MAX( pz[i], zm ) ) // p0-pm
00729                                 {
00730                                     xx[ct] = ( ( clevel[k] - pz[i] ) * ( xm - px[i] ) ) / ( zm - pz[i] ) + px[i];
00731                                     yy[ct] = ( ( clevel[k] - pz[i] ) * ( ym - py[i] ) ) / ( zm - pz[i] ) + py[i];
00732                                     ct++;
00733                                 }
00734 
00735                                 if ( clevel[k] >= MIN( pz[i], pz[j] ) && clevel[k] < MAX( pz[i], pz[j] ) ) // p0-p1
00736                                 {
00737                                     xx[ct] = ( ( clevel[k] - pz[i] ) * ( px[j] - px[i] ) ) / ( pz[j] - pz[i] ) + px[i];
00738                                     yy[ct] = ( ( clevel[k] - pz[i] ) * ( py[j] - py[i] ) ) / ( pz[j] - pz[i] ) + py[i];
00739                                     ct++;
00740                                 }
00741 
00742                                 if ( clevel[k] >= MIN( pz[j], zm ) && clevel[k] < MAX( pz[j], zm ) ) // p1-pm
00743                                 {
00744                                     xx[ct] = ( ( clevel[k] - pz[j] ) * ( xm - px[j] ) ) / ( zm - pz[j] ) + px[j];
00745                                     yy[ct] = ( ( clevel[k] - pz[j] ) * ( ym - py[j] ) ) / ( zm - pz[j] ) + py[j];
00746                                     ct++;
00747                                 }
00748 
00749                                 if ( ct == 2 )
00750                                 {
00751                                     // yes, xx and yy are the intersection points of the triangle with
00752                                     // the contour line -- draw a straight line betweeen the points
00753                                     // -- at the end this will make up the contour line
00754 
00755                                     if ( opt & SURF_CONT )
00756                                     {
00757                                         // surface contour with color set by user
00758                                         plcol0( color );
00759                                         zz[0] = zz[1] = clevel[k];
00760                                         plline3( 2, xx, yy, zz );
00761                                     }
00762 
00763                                     // don't break; one triangle can span various contour levels
00764                                 }
00765                                 else
00766                                     plwarn( "plsurf3dl: ***ERROR***\n" );
00767                             }
00768                         }
00769                     }
00770                 }
00771             }
00772         }
00773     }
00774 
00775     if ( opt & FACETED )
00776     {
00777         plcol0( 0 );
00778         plfplot3dcl( x, y, zops, zp, nx, ny, MESH | DRAW_LINEXY, NULL, 0,
00779             indexxmin, indexxmax, indexymin, indexymax );
00780     }
00781 
00782     if ( opt & DRAW_SIDES ) // the sides look ugly !!!
00783     {                       // draw one more row with all the Z's set to zmin
00784         plP_grange( &zscale, &zmin, &zmax );
00785 
00786         iSlow      = nSlow - 1;
00787         iftriangle = 1;
00788         for ( iFast = 0; iFast < nFast - 1; iFast++ )
00789         {
00790             for ( i = 0; i < 2; i++ )
00791             {
00792                 ix = ixFast * ( iFast + i ) + ixSlow * iSlow + ixOrigin;
00793                 iy = iyFast * ( iFast + i ) + iySlow * iSlow + iyOrigin;
00794                 if ( indexxmin <= ix && ix < indexxmax &&
00795                      indexymin[ix] <= iy && iy < indexymax[ix] )
00796                 {
00797                     px[2 * i] = x[ix];
00798                     py[2 * i] = y[iy];
00799                     pz[2 * i] = getz( zp, ix, iy );
00800                 }
00801                 else
00802                 {
00803                     iftriangle = 0;
00804                     break;
00805                 }
00806             }
00807             if ( iftriangle == 0 )
00808                 break;
00809             // now draw the quad as two triangles (4 might be better)
00810 
00811             shade_triangle( px[0], py[0], pz[0], px[2], py[2], pz[2], px[0], py[0], zmin );
00812             shade_triangle( px[2], py[2], pz[2], px[2], py[2], zmin, px[0], py[0], zmin );
00813         }
00814 
00815         iFast      = nFast - 1;
00816         iftriangle = 1;
00817         for ( iSlow = 0; iSlow < nSlow - 1; iSlow++ )
00818         {
00819             for ( i = 0; i < 2; i++ )
00820             {
00821                 ix = ixFast * iFast + ixSlow * ( iSlow + i ) + ixOrigin;
00822                 iy = iyFast * iFast + iySlow * ( iSlow + i ) + iyOrigin;
00823                 if ( indexxmin <= ix && ix < indexxmax &&
00824                      indexymin[ix] <= iy && iy < indexymax[ix] )
00825                 {
00826                     px[2 * i] = x[ix];
00827                     py[2 * i] = y[iy];
00828                     pz[2 * i] = getz( zp, ix, iy );
00829                 }
00830                 else
00831                 {
00832                     iftriangle = 0;
00833                     break;
00834                 }
00835             }
00836             if ( iftriangle == 0 )
00837                 break;
00838 
00839             // now draw the quad as two triangles (4 might be better)
00840             shade_triangle( px[0], py[0], pz[0], px[2], py[2], pz[2], px[0], py[0], zmin );
00841             shade_triangle( px[2], py[2], pz[2], px[2], py[2], zmin, px[0], py[0], zmin );
00842         }
00843     }
00844 }
00845 
00846 //--------------------------------------------------------------------------
00847 // void plot3d(x, y, z, nx, ny, opt, side)
00848 //
00849 // Plots a 3-d representation of the function z[x][y]. The x values
00850 // are stored as x[0..nx-1], the y values as y[0..ny-1], and the z
00851 // values are in the 2-d array z[][]. The integer "opt" specifies:
00852 // see plot3dcl() below
00853 //--------------------------------------------------------------------------
00854 
00855 void
00856 c_plot3d( const PLFLT *x, const PLFLT *y, const PLFLT * const *z,
00857           PLINT nx, PLINT ny, PLINT opt, PLBOOL side )
00858 {
00859     plfplot3dc( x, y, plf2ops_c(), (PLPointer) z, nx, ny, opt | ( side != 0 ? DRAW_SIDES : 0 ), NULL, 0 );
00860 }
00861 
00862 void
00863 plfplot3d( const PLFLT *x, const PLFLT *y, PLF2OPS zops, PLPointer zp,
00864            PLINT nx, PLINT ny, PLINT opt, PLBOOL side )
00865 {
00866     plfplot3dc( x, y, zops, zp, nx, ny, opt | ( side != 0 ? DRAW_SIDES : 0 ), NULL, 0 );
00867 }
00868 
00869 //--------------------------------------------------------------------------
00870 // void plot3dc(x, y, z, nx, ny, opt, clevel, nlevel)
00871 //
00872 // Plots a 3-d representation of the function z[x][y]. The x values
00873 // are stored as x[0..nx-1], the y values as y[0..ny-1], and the z
00874 // values are in the 2-d array z[][]. The integer "opt" specifies:
00875 // see plot3dcl() below
00876 //--------------------------------------------------------------------------
00877 
00878 void
00879 c_plot3dc( const PLFLT *x, const PLFLT *y, const PLFLT * const*z,
00880            PLINT nx, PLINT ny, PLINT opt,
00881            const PLFLT *clevel, PLINT nlevel )
00882 {
00883     plfplot3dcl( x, y, plf2ops_c(), (PLPointer) z, nx, ny, opt, clevel, nlevel, 0, 0, NULL, NULL );
00884 }
00885 
00886 void
00887 plfplot3dc( const PLFLT *x, const PLFLT *y, PLF2OPS zops, PLPointer zp,
00888             PLINT nx, PLINT ny, PLINT opt, const PLFLT *clevel, PLINT nlevel )
00889 {
00890     plfplot3dcl( x, y, zops, zp, nx, ny, opt, clevel, nlevel, 0, 0, NULL, NULL );
00891 }
00892 
00893 //--------------------------------------------------------------------------
00894 // void plot3dcl(x, y, z, nx, ny, opt, clevel, nlevel,
00895 //       indexxmin, indexxmax, indexymin, indexymax)
00896 //
00897 // Plots a 3-d representation of the function z[x][y]. The x values
00898 // are stored as x[0..nx-1], the y values as y[0..ny-1], and the z
00899 // values are in the 2-d array z[][]. The integer "opt" specifies:
00900 //
00901 //  DRAW_LINEX :  Draw lines parallel to x-axis
00902 //  DRAW_LINEY :  Draw lines parallel to y-axis
00903 //  DRAW_LINEXY:  Draw lines parallel to both axes
00904 //  MAG_COLOR:    Magnitude coloring of wire frame
00905 //  BASE_CONT:    Draw contour at bottom xy plane
00906 //  TOP_CONT:     Draw contour at top xy plane (not yet)
00907 //  DRAW_SIDES:   Draw sides around the plot
00908 //  MESH:       Draw the "under" side of the plot
00909 //
00910 // or any bitwise combination, e.g. "MAG_COLOR | DRAW_LINEX"
00911 // indexymin and indexymax are arrays which specify the y index limits
00912 // (following the convention that the upper range limit is one more than
00913 // actual index limit) for an x index range of indexxmin, indexxmax.
00914 //--------------------------------------------------------------------------
00915 
00916 void
00917 c_plot3dcl( const PLFLT *x, const PLFLT *y, const PLFLT * const *z,
00918             PLINT nx, PLINT ny, PLINT opt,
00919             const PLFLT *clevel, PLINT nlevel,
00920             PLINT indexxmin, PLINT indexxmax, const PLINT *indexymin, const PLINT *indexymax )
00921 {
00922     plfplot3dcl( x, y, plf2ops_c(), (PLPointer) z, nx, ny,
00923         opt, clevel, nlevel, indexxmin, indexxmax, indexymin, indexymax );
00924 }
00925 
00926 //--------------------------------------------------------------------------
00963 //--------------------------------------------------------------------------
00964 
00965 void
00966 plfplot3dcl( const PLFLT *x, const PLFLT *y, PLF2OPS zops, PLPointer zp,
00967              PLINT nx, PLINT ny, PLINT opt,
00968              const PLFLT *clevel, PLINT nlevel,
00969              PLINT PL_UNUSED( indexxmin ), PLINT PL_UNUSED( indexxmax ), const PLINT * PL_UNUSED( indexymin ), const PLINT * PL_UNUSED( indexymax ) )
00970 {
00971     PLFLT cxx, cxy, cyx, cyy, cyz;
00972     PLINT init, ix, iy, color;
00973     PLFLT width;
00974     PLFLT xmin, xmax, ymin, ymax, zmin, zmax, zscale;
00975     PLINT ixmin   = 0, ixmax = nx - 1, iymin = 0, iymax = ny - 1;
00976     PLINT clipped = 0, base_cont = 0, side = 0;
00977     PLFLT ( *getz )( PLPointer, PLINT, PLINT ) = zops->get;
00978     PLFLT *_x = NULL, *_y = NULL, **_z = NULL;
00979     const PLFLT *x_modified, *y_modified;
00980     int i;
00981 
00982     pl3mode = 0;
00983 
00984     if ( plsc->level < 3 )
00985     {
00986         myabort( "plot3dcl: Please set up window first" );
00987         return;
00988     }
00989 
00990     if ( ( opt & 3 ) == 0 )
00991     {
00992         myabort( "plot3dcl: Bad option" );
00993         return;
00994     }
00995 
00996     if ( nx <= 0 || ny <= 0 )
00997     {
00998         myabort( "plot3dcl: Bad array dimensions." );
00999         return;
01000     }
01001 
01002     plP_gdom( &xmin, &xmax, &ymin, &ymax );
01003     plP_grange( &zscale, &zmin, &zmax );
01004     if ( zmin > zmax )
01005     {
01006         PLFLT t = zmin;
01007         zmin = zmax;
01008         zmax = t;
01009     }
01010 
01011 // Check that points in x and in y are strictly increasing
01012 
01013     for ( i = 0; i < nx - 1; i++ )
01014     {
01015         if ( x[i] >= x[i + 1] )
01016         {
01017             myabort( "plot3dcl: X array must be strictly increasing" );
01018             return;
01019         }
01020     }
01021     for ( i = 0; i < ny - 1; i++ )
01022     {
01023         if ( y[i] >= y[i + 1] )
01024         {
01025             myabort( "plot3dcl: Y array must be strictly increasing" );
01026             return;
01027         }
01028     }
01029 
01030     if ( opt & MESH )
01031         pl3mode = 1;
01032 
01033     if ( opt & DRAW_SIDES )
01034         side = 1;
01035 
01036     // figure out the part of the data to use
01037     if ( xmin < x[0] )
01038         xmin = x[0];
01039     if ( xmax > x[nx - 1] )
01040         xmax = x[nx - 1];
01041     if ( ymin < y[0] )
01042         ymin = y[0];
01043     if ( ymax > y[ny - 1] )
01044         ymax = y[ny - 1];
01045     for ( ixmin = 0; ixmin < nx - 1 && x[ixmin + 1] < xmin; ixmin++ )
01046     {
01047     }
01048     for ( ixmax = nx - 1; ixmax > 0 && x[ixmax - 1] > xmax; ixmax-- )
01049     {
01050     }
01051     for ( iymin = 0; iymin < ny - 1 && y[iymin + 1] < ymin; iymin++ )
01052     {
01053     }
01054     for ( iymax = ny - 1; iymax > 0 && y[iymax - 1] > ymax; iymax-- )
01055     {
01056     }
01057     //fprintf(stderr, "(%d,%d) %d %d %d %d\n", nx, ny, ixmin, ixmax, iymin, iymax);
01058     // do we need to clip?
01059     if ( ixmin > 0 || ixmax < nx - 1 || iymin > 0 || iymax < ny - 1 )
01060     {
01061         // adjust the input so it stays within bounds
01062         int _nx = ixmax - ixmin + 1;
01063         int _ny = iymax - iymin + 1;
01064         PLFLT ty0, ty1, tx0, tx1;
01065         int j;
01066 
01067         if ( _nx <= 1 || _ny <= 1 )
01068         {
01069             myabort( "plot3dcl: selected x or y range has no data" );
01070             return;
01071         }
01072 
01073         // allocate storage for new versions of the input vectors
01074         if ( ( ( _x = (PLFLT *) malloc( (size_t) _nx * sizeof ( PLFLT ) ) ) == NULL ) ||
01075              ( ( _y = (PLFLT *) malloc( (size_t) _ny * sizeof ( PLFLT ) ) ) == NULL ) ||
01076              ( ( _z = (PLFLT **) malloc( (size_t) _nx * sizeof ( PLFLT* ) ) ) == NULL ) )
01077         {
01078             plexit( "c_plot3dcl: Insufficient memory" );
01079         }
01080 
01081         clipped = 1;
01082 
01083         // copy over the independent variables
01084         _x[0]       = xmin;
01085         _x[_nx - 1] = xmax;
01086         for ( i = 1; i < _nx - 1; i++ )
01087             _x[i] = x[ixmin + i];
01088         _y[0]       = ymin;
01089         _y[_ny - 1] = ymax;
01090         for ( i = 1; i < _ny - 1; i++ )
01091             _y[i] = y[iymin + i];
01092 
01093         // copy the data array so we can interpolate around the edges
01094         for ( i = 0; i < _nx; i++ )
01095         {
01096             if ( ( _z[i] = (PLFLT *) malloc( (size_t) _ny * sizeof ( PLFLT ) ) ) == NULL )
01097             {
01098                 plexit( "c_plot3dcl: Insufficient memory" );
01099             }
01100         }
01101 
01102         // interpolation factors for the 4 edges
01103         ty0 = ( _y[0] - y[iymin] ) / ( y[iymin + 1] - y[iymin] );
01104         ty1 = ( _y[_ny - 1] - y[iymax - 1] ) / ( y[iymax] - y[iymax - 1] );
01105         tx0 = ( _x[0] - x[ixmin] ) / ( x[ixmin + 1] - x[ixmin] );
01106         tx1 = ( _x[_nx - 1] - x[ixmax - 1] ) / ( x[ixmax] - x[ixmax - 1] );
01107         for ( i = 0; i < _nx; i++ )
01108         {
01109             if ( i == 0 )
01110             {
01111                 _z[i][0] = getz( zp, ixmin, iymin ) * ( 1 - ty0 ) * ( 1 - tx0 ) + getz( zp, ixmin, iymin + 1 ) * ( 1 - tx0 ) * ty0
01112                            + getz( zp, ixmin + 1, iymin ) * tx0 * ( 1 - ty0 ) + getz( zp, ixmin + 1, iymin + 1 ) * tx0 * ty0;
01113                 for ( j = 1; j < _ny - 1; j++ )
01114                     _z[i][j] = getz( zp, ixmin, iymin + j ) * ( 1 - tx0 ) + getz( zp, ixmin + 1, iymin + j ) * tx0;
01115                 _z[i][_ny - 1] = getz( zp, ixmin, iymax - 1 ) * ( 1 - tx0 ) * ( 1 - ty1 ) + getz( zp, ixmin, iymax ) * ( 1 - tx0 ) * ty1
01116                                  + getz( zp, ixmin + 1, iymax - 1 ) * tx0 * ( 1 - ty1 ) + getz( zp, ixmin + 1, iymax ) * tx0 * ty1;
01117             }
01118             else if ( i == _nx - 1 )
01119             {
01120                 _z[i][0] = getz( zp, ixmax - 1, iymin ) * ( 1 - tx1 ) * ( 1 - ty0 ) + getz( zp, ixmax - 1, iymin + 1 ) * ( 1 - tx1 ) * ty0
01121                            + getz( zp, ixmax, iymin ) * tx1 * ( 1 - ty0 ) + getz( zp, ixmax, iymin + 1 ) * tx1 * ty0;
01122                 for ( j = 1; j < _ny - 1; j++ )
01123                     _z[i][j] = getz( zp, ixmax - 1, iymin + j ) * ( 1 - tx1 ) + getz( zp, ixmax, iymin + j ) * tx1;
01124                 _z[i][_ny - 1] = getz( zp, ixmax - 1, iymax - 1 ) * ( 1 - tx1 ) * ( 1 - ty1 ) + getz( zp, ixmax, iymax ) * ( 1 - tx1 ) * ty1
01125                                  + getz( zp, ixmax, iymax - 1 ) * tx1 * ( 1 - ty1 ) + getz( zp, ixmax, iymax ) * tx1 * ty1;
01126             }
01127             else
01128             {
01129                 _z[i][0] = getz( zp, ixmin + i, iymin ) * ( 1 - ty0 ) + getz( zp, ixmin + i, iymin + 1 ) * ty0;
01130                 for ( j = 1; j < _ny - 1; j++ )
01131                     _z[i][j] = getz( zp, ixmin + i, iymin + j );
01132                 _z[i][_ny - 1] = getz( zp, ixmin + i, iymax - 1 ) * ( 1 - ty1 ) + getz( zp, ixmin + i, iymax ) * ty1;
01133             }
01134             for ( j = 0; j < _ny; j++ )
01135             {
01136                 if ( _z[i][j] < zmin )
01137                     _z[i][j] = zmin;
01138                 else if ( _z[i][j] > zmax )
01139                     _z[i][j] = zmax;
01140             }
01141         }
01142         // replace the input with our clipped versions
01143         zp   = (PLPointer) _z;
01144         getz = plf2ops_c()->get;
01145         nx   = _nx;
01146         ny   = _ny;
01147         // Do not want to modify input x and y (const modifier)
01148         x_modified = (const PLFLT *) _x;
01149         y_modified = (const PLFLT *) _y;
01150     }
01151     else
01152     {
01153         x_modified = x;
01154         y_modified = y;
01155     }
01156 
01157     // From here on must use x_modified and y_modified rather than
01158     // x and y.
01159     if ( ( opt & BASE_CONT ) || ( opt & TOP_CONT ) || ( opt && MAG_COLOR ) )
01160     {
01161         //
01162         // Don't use the data z value to scale the color, use the z axis
01163         // values set by plw3d()
01164         //
01165         // plMinMax2dGrid(z, nx, ny, &fc_maxz, &fc_minz);
01166         //
01167 
01168         fc_minz = plsc->ranmi;
01169         fc_maxz = plsc->ranma;
01170 
01171         if ( fc_maxz == fc_minz )
01172         {
01173             plwarn( "plot3dcl: Maximum and minimum Z values are equal! \"fixing\"..." );
01174             fc_maxz = fc_minz + 1e-6;
01175         }
01176     }
01177 
01178     if ( opt & BASE_CONT )    // If enabled, draw the contour at the base.
01179     {
01180         if ( clevel != NULL && nlevel != 0 )
01181         {
01182             base_cont = 1;
01183             // even if MESH is not set, "set it",
01184             // as the base contour can only be done in this case
01185             pl3mode = 1;
01186         }
01187     }
01188 
01189     if ( opt & MAG_COLOR )    // If enabled, use magnitude colored wireframe
01190     {
01191         if ( ( ctmp = (PLFLT *) malloc( (size_t) ( 2 * MAX( nx, ny ) ) * sizeof ( PLFLT ) ) ) == NULL )
01192         {
01193             plexit( "c_plot3dcl: Insufficient memory" );
01194         }
01195     }
01196     else
01197         ctmp = NULL;
01198 
01199     // next logic only knows opt = 1 | 2 | 3, make sure that it only gets that
01200     opt &= DRAW_LINEXY;
01201 
01202     // Allocate work arrays
01203 
01204     utmp = (PLINT *) malloc( (size_t) ( 2 * MAX( nx, ny ) ) * sizeof ( PLINT ) );
01205     vtmp = (PLINT *) malloc( (size_t) ( 2 * MAX( nx, ny ) ) * sizeof ( PLINT ) );
01206 
01207     if ( !utmp || !vtmp )
01208         myexit( "plot3dcl: Out of memory." );
01209 
01210     plP_gw3wc( &cxx, &cxy, &cyx, &cyy, &cyz );
01211     init = 1;
01212 // Call 3d line plotter.  Each viewing quadrant
01213 // (perpendicular to x-y plane) must be handled separately.
01214     if ( cxx >= 0.0 && cxy <= 0.0 )
01215     {
01216         if ( opt == DRAW_LINEY )
01217             plt3zz( 1, ny, 1, -1, -opt, &init, x_modified, y_modified, zops, zp, nx, ny, utmp, vtmp, ctmp );
01218         else
01219         {
01220             for ( iy = 2; iy <= ny; iy++ )
01221                 plt3zz( 1, iy, 1, -1, -opt, &init, x_modified, y_modified, zops, zp, nx, ny, utmp, vtmp, ctmp );
01222         }
01223         if ( opt == DRAW_LINEX )
01224             plt3zz( 1, ny, 1, -1, opt, &init, x_modified, y_modified, zops, zp, nx, ny, utmp, vtmp, ctmp );
01225         else
01226         {
01227             for ( ix = 1; ix <= nx - 1; ix++ )
01228                 plt3zz( ix, ny, 1, -1, opt, &init, x_modified, y_modified, zops, zp, nx, ny, utmp, vtmp, ctmp );
01229         }
01230     }
01231 
01232     else if ( cxx <= 0.0 && cxy <= 0.0 )
01233     {
01234         if ( opt == DRAW_LINEX )
01235             plt3zz( nx, ny, -1, -1, opt, &init, x_modified, y_modified, zops, zp, nx, ny, utmp, vtmp, ctmp );
01236         else
01237         {
01238             for ( ix = 2; ix <= nx; ix++ )
01239                 plt3zz( ix, ny, -1, -1, opt, &init, x_modified, y_modified, zops, zp, nx, ny, utmp, vtmp, ctmp );
01240         }
01241         if ( opt == DRAW_LINEY )
01242             plt3zz( nx, ny, -1, -1, -opt, &init, x_modified, y_modified, zops, zp, nx, ny, utmp, vtmp, ctmp );
01243         else
01244         {
01245             for ( iy = ny; iy >= 2; iy-- )
01246                 plt3zz( nx, iy, -1, -1, -opt, &init, x_modified, y_modified, zops, zp, nx, ny, utmp, vtmp, ctmp );
01247         }
01248     }
01249 
01250     else if ( cxx <= 0.0 && cxy >= 0.0 )
01251     {
01252         if ( opt == DRAW_LINEY )
01253             plt3zz( nx, 1, -1, 1, -opt, &init, x_modified, y_modified, zops, zp, nx, ny, utmp, vtmp, ctmp );
01254         else
01255         {
01256             for ( iy = ny - 1; iy >= 1; iy-- )
01257                 plt3zz( nx, iy, -1, 1, -opt, &init, x_modified, y_modified, zops, zp, nx, ny, utmp, vtmp, ctmp );
01258         }
01259         if ( opt == DRAW_LINEX )
01260             plt3zz( nx, 1, -1, 1, opt, &init, x_modified, y_modified, zops, zp, nx, ny, utmp, vtmp, ctmp );
01261         else
01262         {
01263             for ( ix = nx; ix >= 2; ix-- )
01264                 plt3zz( ix, 1, -1, 1, opt, &init, x_modified, y_modified, zops, zp, nx, ny, utmp, vtmp, ctmp );
01265         }
01266     }
01267 
01268     else if ( cxx >= 0.0 && cxy >= 0.0 )
01269     {
01270         if ( opt == DRAW_LINEX )
01271             plt3zz( 1, 1, 1, 1, opt, &init, x_modified, y_modified, zops, zp, nx, ny, utmp, vtmp, ctmp );
01272         else
01273         {
01274             for ( ix = nx - 1; ix >= 1; ix-- )
01275                 plt3zz( ix, 1, 1, 1, opt, &init, x_modified, y_modified, zops, zp, nx, ny, utmp, vtmp, ctmp );
01276         }
01277         if ( opt == DRAW_LINEY )
01278             plt3zz( 1, 1, 1, 1, -opt, &init, x_modified, y_modified, zops, zp, nx, ny, utmp, vtmp, ctmp );
01279         else
01280         {
01281             for ( iy = 1; iy <= ny - 1; iy++ )
01282                 plt3zz( 1, iy, 1, 1, -opt, &init, x_modified, y_modified, zops, zp, nx, ny, utmp, vtmp, ctmp );
01283         }
01284     }
01285 
01286     // draw contour at the base. Not 100%! Why?
01287 
01288     if ( base_cont )
01289     {
01290         int np = NPTS, j;
01291         CONT_LEVEL *cont, *clev;
01292         CONT_LINE *cline;
01293 
01294         PLINT *uu = (PLINT *) malloc( (size_t) NPTS * sizeof ( PLINT ) );
01295         PLINT *vv = (PLINT *) malloc( (size_t) NPTS * sizeof ( PLINT ) );
01296         // prepare cont_store input
01297         PLFLT **zstore;
01298         PLcGrid2 cgrid2;
01299 
01300         if ( ( uu == NULL ) || ( vv == NULL ) )
01301         {
01302             plexit( "c_plot3dcl: Insufficient memory" );
01303         }
01304 
01305         cgrid2.nx = nx;
01306         cgrid2.ny = ny;
01307         plAlloc2dGrid( &cgrid2.xg, nx, ny );
01308         plAlloc2dGrid( &cgrid2.yg, nx, ny );
01309         plAlloc2dGrid( &zstore, nx, ny );
01310 
01311         for ( i = 0; i < nx; i++ )
01312         {
01313             for ( j = 0; j < ny; j++ )
01314             {
01315                 cgrid2.xg[i][j] = x_modified[i];
01316                 cgrid2.yg[i][j] = y_modified[j];
01317                 zstore[i][j]    = getz( zp, i, j );
01318             }
01319         }
01320 
01321         pl3upv = 0;
01322 
01323         // Fill cont structure with contours.
01324         cont_store( (const PLFLT * const *) zstore, nx, ny, 1, nx, 1, ny,
01325             clevel, nlevel, pltr2, (void *) &cgrid2, &cont );
01326 
01327         // Free the 2D input arrays to cont_store since not needed any more.
01328         plFree2dGrid( zstore, nx, ny );
01329         plFree2dGrid( cgrid2.xg, nx, ny );
01330         plFree2dGrid( cgrid2.yg, nx, ny );
01331 
01332         // follow the contour levels and lines
01333         clev = cont;
01334         do  // for each contour level
01335         {
01336             cline = clev->line;
01337             do  // there are several lines that make up each contour
01338             {
01339                 int cx, k, l, m, start, end;
01340                 PLFLT tx, ty;
01341                 if ( cline->npts > np )
01342                 {
01343                     np = cline->npts;
01344                     if ( ( ( uu = (PLINT *) realloc( uu, (size_t) np * sizeof ( PLINT ) ) ) == NULL ) ||
01345                          ( ( vv = (PLINT *) realloc( vv, (size_t) np * sizeof ( PLINT ) ) ) == NULL ) )
01346                     {
01347                         plexit( "c_plot3dcl: Insufficient memory" );
01348                     }
01349                 }
01350 
01351                 // the hidden line plotter plnxtv() only works OK if the x points are in increasing order.
01352                 // As this does not always happens, the situation must be detected and the line segment
01353                 // must be reversed before being plotted
01354                 i = 0;
01355                 if ( cline->npts > 0 )
01356                 {
01357                     do
01358                     {
01359                         plcol1( ( clev->level - fc_minz ) / ( fc_maxz - fc_minz ) );
01360                         cx = plP_wcpcx( plP_w3wcx( cline->x[i], cline->y[i], plsc->ranmi ) );
01361                         for ( j = i; j < cline->npts; j++ ) // convert to 2D coordinates
01362                         {
01363                             uu[j] = plP_wcpcx( plP_w3wcx( cline->x[j], cline->y[j], plsc->ranmi ) );
01364                             vv[j] = plP_wcpcy( plP_w3wcy( cline->x[j], cline->y[j], plsc->ranmi ) );
01365                             if ( uu[j] < cx ) // find turn back point
01366                                 break;
01367                             else
01368                                 cx = uu[j];
01369                         }
01370                         plnxtv( &uu[i], &vv[i], NULL, j - i, 0 ); // plot line with increasing x
01371 
01372                         if ( j < cline->npts )                    // line not yet finished,
01373                         {
01374                             start = j - 1;
01375                             for ( i = start; i < cline->npts; i++ ) // search turn forward point
01376                             {
01377                                 uu[i] = plP_wcpcx( plP_w3wcx( cline->x[i], cline->y[i], plsc->ranmi ) );
01378                                 if ( uu[i] > cx )
01379                                     break;
01380                                 else
01381                                     cx = uu[i];
01382                             }
01383                             end = i - 1;
01384 
01385                             for ( k = 0; k < ( end - start + 1 ) / 2; k++ ) // reverse line segment
01386                             {
01387                                 l           = start + k;
01388                                 m           = end - k;
01389                                 tx          = cline->x[l];
01390                                 ty          = cline->y[l];
01391                                 cline->x[l] = cline->x[m];
01392                                 cline->y[l] = cline->y[m];
01393                                 cline->x[m] = tx;
01394                                 cline->y[m] = ty;
01395                             }
01396 
01397                             // convert to 2D coordinates
01398                             for ( j = start; j <= end; j++ )
01399                             {
01400                                 uu[j] = plP_wcpcx( plP_w3wcx( cline->x[j], cline->y[j], plsc->ranmi ) );
01401                                 vv[j] = plP_wcpcy( plP_w3wcy( cline->x[j], cline->y[j], plsc->ranmi ) );
01402                             }
01403                             plnxtv( &uu[start], &vv[start], NULL, end - start + 1, 0 ); // and plot it
01404 
01405                             cline->x[end] = cline->x[start];
01406                             cline->y[end] = cline->y[start];
01407                             i             = end; // restart where it was left
01408                         }
01409                     } while ( j < cline->npts && i < cline->npts );
01410                 }
01411                 cline = cline->next;
01412             }
01413             while ( cline != NULL );
01414             clev = clev->next;
01415         }
01416         while ( clev != NULL );
01417 
01418         cont_clean_store( cont ); // now release the contour memory
01419         pl3upv = 1;
01420         free( uu );
01421         free( vv );
01422     }
01423 
01424 // Finish up by drawing sides, background grid (both are optional)
01425 
01426     if ( side )
01427         plside3( x_modified, y_modified, zops, zp, nx, ny, opt );
01428 
01429     if ( zbflg )
01430     {
01431         color = plsc->icol0;
01432         width = plsc->width;
01433         plwidth( zbwidth );
01434         plcol0( zbcol );
01435         plgrid3( zbtck );
01436         plwidth( width );
01437         plcol0( color );
01438     }
01439 
01440     freework();
01441 
01442     if ( clipped )
01443     {
01444         free( _x );
01445         free( _y );
01446         for ( i = 0; i < nx; i++ )
01447             free( _z[i] );
01448         free( _z );
01449     }
01450 }
01451 
01452 //--------------------------------------------------------------------------
01453 // void plxyindexlimits()
01454 //
01455 // Transform from y array limits to corresponding x array limits (or vice
01456 // versa).
01457 //
01458 // N.B. we follow the convention here that all upper range limits are one
01459 // more than the actual last index.
01460 // instart (>= 0) through inn is the index range where
01461 // the input inarray_min and inarray_max arrays are defined.
01462 //
01463 // outstart (>= 0), through outn (with outn <= outnmax) is the index range
01464 // where the output outarray_min and outarray_max arrays are defined.
01465 //
01466 // In order to assure the transformation from y array limits to x array limits
01467 // (or vice versa) is single-valued, this programme plaborts if the
01468 // inarray_min array has a maximum or inarray_max array has a minimum.
01469 //--------------------------------------------------------------------------
01470 
01471 //static void
01472 //plxyindexlimits( PLINT instart, PLINT inn,
01473 //                 PLINT *inarray_min, PLINT *inarray_max,
01474 //                 PLINT *outstart, PLINT *outn, PLINT outnmax,
01475 //                 PLINT *outarray_min, PLINT *outarray_max )
01476 //{
01477 //    PLINT i, j;
01478 //    if ( inn < 0 )
01479 //    {
01480 //        myabort( "plxyindexlimits: Must have instart >= 0" );
01481 //        return;
01482 //    }
01483 //    if ( inn - instart <= 0 )
01484 //    {
01485 //        myabort( "plxyindexlimits: Must have at least 1 defined point" );
01486 //        return;
01487 //    }
01488 //    *outstart = inarray_min[instart];
01489 //    *outn     = inarray_max[instart];
01490 //    for ( i = instart; i < inn; i++ )
01491 //    {
01492 //        *outstart = MIN( *outstart, inarray_min[i] );
01493 //        *outn     = MAX( *outn, inarray_max[i] );
01494 //        if ( i + 2 < inn )
01495 //        {
01496 //            if ( inarray_min[i] < inarray_min[i + 1] &&
01497 //                 inarray_min[i + 1] > inarray_min[i + 2] )
01498 //            {
01499 //                myabort( "plxyindexlimits: inarray_min must not have a maximum" );
01500 //                return;
01501 //            }
01502 //            if ( inarray_max[i] > inarray_max[i + 1] &&
01503 //                 inarray_max[i + 1] < inarray_max[i + 2] )
01504 //            {
01505 //                myabort( "plxyindexlimits: inarray_max must not have a minimum" );
01506 //                return;
01507 //            }
01508 //        }
01509 //    }
01510 //    if ( *outstart < 0 )
01511 //    {
01512 //        myabort( "plxyindexlimits: Must have all elements of inarray_min >= 0" );
01513 //        return;
01514 //    }
01515 //    if ( *outn > outnmax )
01516 //    {
01517 //        myabort( "plxyindexlimits: Must have all elements of inarray_max <= outnmax" );
01518 //        return;
01519 //    }
01520 //    for ( j = *outstart; j < *outn; j++ )
01521 //    {
01522 //        i = instart;
01523 //        // Find first valid i for this j.
01524 //        while ( i < inn && !( inarray_min[i] <= j && j < inarray_max[i] ) )
01525 //            i++;
01526 //        if ( i < inn )
01527 //            outarray_min[j] = i;
01528 //        else
01529 //        {
01530 //            myabort( "plxyindexlimits: bad logic; invalid i should never happen" );
01531 //            return;
01532 //        }
01533 //        // Find next invalid i for this j.
01534 //        while ( i < inn && ( inarray_min[i] <= j && j < inarray_max[i] ) )
01535 //            i++;
01536 //        outarray_max[j] = i;
01537 //    }
01538 //}
01539 
01540 //--------------------------------------------------------------------------
01541 // void plP_gzback()
01542 //
01543 // Get background parameters for 3d plot.
01544 //--------------------------------------------------------------------------
01545 
01546 void
01547 plP_gzback( PLINT **zbf, PLINT **zbc, PLFLT **zbt, PLFLT **zbw )
01548 {
01549     *zbf = &zbflg;
01550     *zbc = &zbcol;
01551     *zbt = &zbtck;
01552     *zbw = &zbwidth;
01553 }
01554 
01555 //--------------------------------------------------------------------------
01556 // PLFLT plGetAngleToLight()
01557 //
01558 // Gets cos of angle between normal to a polygon and a light source.
01559 // Requires at least 3 elements, forming non-parallel lines
01560 // in the arrays.
01561 //--------------------------------------------------------------------------
01562 
01563 static PLFLT
01564 plGetAngleToLight( PLFLT* x, PLFLT* y, PLFLT* z )
01565 {
01566     PLFLT vx1, vx2, vy1, vy2, vz1, vz2;
01567     PLFLT px, py, pz;
01568     PLFLT vlx, vly, vlz;
01569     PLFLT mag1, mag2;
01570     PLFLT cosangle;
01571 
01572     vx1 = x[1] - x[0];
01573     vx2 = x[2] - x[1];
01574     vy1 = y[1] - y[0];
01575     vy2 = y[2] - y[1];
01576     vz1 = z[1] - z[0];
01577     vz2 = z[2] - z[1];
01578 
01579 // Find vector perpendicular to the face
01580     px   = vy1 * vz2 - vz1 * vy2;
01581     py   = vz1 * vx2 - vx1 * vz2;
01582     pz   = vx1 * vy2 - vy1 * vx2;
01583     mag1 = px * px + py * py + pz * pz;
01584 
01585 // Vectors were parallel!
01586     if ( mag1 == 0 )
01587         return 1;
01588 
01589     vlx  = xlight - x[0];
01590     vly  = ylight - y[0];
01591     vlz  = zlight - z[0];
01592     mag2 = vlx * vlx + vly * vly + vlz * vlz;
01593     if ( mag2 == 0 )
01594         return 1;
01595 
01596 // Now have 3 vectors going through the first point on the given surface
01597     cosangle = fabs( ( vlx * px + vly * py + vlz * pz ) / sqrt( mag1 * mag2 ) );
01598 
01599 // In case of numerical rounding
01600     if ( cosangle > 1 )
01601         cosangle = 1;
01602     return cosangle;
01603 }
01604 
01605 //--------------------------------------------------------------------------
01606 // void plt3zz()
01607 //
01608 // Draws the next zig-zag line for a 3-d plot.  The data is stored in array
01609 // z[][] as a function of x[] and y[], and is plotted out starting at index
01610 // (x0,y0).
01611 //
01612 // Depending on the state of "flag", the sequence of data points sent to
01613 // plnxtv is altered so as to allow cross-hatch plotting, or plotting
01614 // parallel to either the x-axis or the y-axis.
01615 //--------------------------------------------------------------------------
01616 
01617 static void
01618 plt3zz( PLINT x0, PLINT y0, PLINT dx, PLINT dy, PLINT flag, PLINT *init,
01619         const PLFLT *x, const PLFLT *y, PLF2OPS zops, PLPointer zp, PLINT nx, PLINT ny,
01620         PLINT *u, PLINT *v, PLFLT* c )
01621 {
01622     PLINT n = 0;
01623     PLFLT x2d, y2d;
01624     PLFLT ( *getz )( PLPointer, PLINT, PLINT ) = zops->get;
01625 
01626     while ( 1 <= x0 && x0 <= nx && 1 <= y0 && y0 <= ny )
01627     {
01628         x2d  = plP_w3wcx( x[x0 - 1], y[y0 - 1], getz( zp, x0 - 1, y0 - 1 ) );
01629         y2d  = plP_w3wcy( x[x0 - 1], y[y0 - 1], getz( zp, x0 - 1, y0 - 1 ) );
01630         u[n] = plP_wcpcx( x2d );
01631         v[n] = plP_wcpcy( y2d );
01632         if ( c != NULL )
01633             c[n] = ( getz( zp, x0 - 1, y0 - 1 ) - fc_minz ) / ( fc_maxz - fc_minz );
01634 
01635         switch ( flag )
01636         {
01637         case -3:
01638             y0  += dy;
01639             flag = -flag;
01640             break;
01641         case -2:
01642             y0 += dy;
01643             break;
01644         case -1:
01645             y0  += dy;
01646             flag = -flag;
01647             break;
01648         case 1:
01649             x0 += dx;
01650             break;
01651         case 2:
01652             x0  += dx;
01653             flag = -flag;
01654             break;
01655         case 3:
01656             x0  += dx;
01657             flag = -flag;
01658             break;
01659         }
01660         n++;
01661     }
01662 
01663     if ( flag == 1 || flag == -2 )
01664     {
01665         if ( flag == 1 )
01666         {
01667             x0 -= dx;
01668             y0 += dy;
01669         }
01670         else if ( flag == -2 )
01671         {
01672             y0 -= dy;
01673             x0 += dx;
01674         }
01675         if ( 1 <= x0 && x0 <= nx && 1 <= y0 && y0 <= ny )
01676         {
01677             x2d  = plP_w3wcx( x[x0 - 1], y[y0 - 1], getz( zp, x0 - 1, y0 - 1 ) );
01678             y2d  = plP_w3wcy( x[x0 - 1], y[y0 - 1], getz( zp, x0 - 1, y0 - 1 ) );
01679             u[n] = plP_wcpcx( x2d );
01680             v[n] = plP_wcpcy( y2d );
01681             if ( c != NULL )
01682                 c[n] = ( getz( zp, x0 - 1, y0 - 1 ) - fc_minz ) / ( fc_maxz - fc_minz );
01683             n++;
01684         }
01685     }
01686 
01687 // All the setup is done.  Time to do the work.
01688 
01689     plnxtv( u, v, c, n, *init );
01690     *init = 0;
01691 }
01692 
01693 //--------------------------------------------------------------------------
01694 // void plside3()
01695 //
01696 // This routine draws sides around the front of the 3d plot so that
01697 // it does not appear to float.
01698 //--------------------------------------------------------------------------
01699 
01700 static void
01701 plside3( const PLFLT *x, const PLFLT *y, PLF2OPS zops, PLPointer zp, PLINT nx, PLINT ny, PLINT opt )
01702 {
01703     PLINT i;
01704     PLFLT cxx, cxy, cyx, cyy, cyz;
01705     PLFLT xmin, ymin, zmin, xmax, ymax, zmax, zscale;
01706     PLFLT tx, ty, ux, uy;
01707     PLFLT ( *getz )( PLPointer, PLINT, PLINT ) = zops->get;
01708 
01709     plP_gw3wc( &cxx, &cxy, &cyx, &cyy, &cyz );
01710     plP_gdom( &xmin, &xmax, &ymin, &ymax );
01711     plP_grange( &zscale, &zmin, &zmax );
01712 
01713 // Get x, y coordinates of legs and plot
01714 
01715     if ( cxx >= 0.0 && cxy <= 0.0 )
01716     {
01717         if ( opt != 1 )
01718         {
01719             for ( i = 0; i < nx; i++ )
01720             {
01721                 tx = plP_w3wcx( x[i], y[0], zmin );
01722                 ty = plP_w3wcy( x[i], y[0], zmin );
01723                 ux = plP_w3wcx( x[i], y[0], getz( zp, i, 0 ) );
01724                 uy = plP_w3wcy( x[i], y[0], getz( zp, i, 0 ) );
01725                 pljoin( tx, ty, ux, uy );
01726             }
01727         }
01728         if ( opt != 2 )
01729         {
01730             for ( i = 0; i < ny; i++ )
01731             {
01732                 tx = plP_w3wcx( x[0], y[i], zmin );
01733                 ty = plP_w3wcy( x[0], y[i], zmin );
01734                 ux = plP_w3wcx( x[0], y[i], getz( zp, 0, i ) );
01735                 uy = plP_w3wcy( x[0], y[i], getz( zp, 0, i ) );
01736                 pljoin( tx, ty, ux, uy );
01737             }
01738         }
01739     }
01740     else if ( cxx <= 0.0 && cxy <= 0.0 )
01741     {
01742         if ( opt != 1 )
01743         {
01744             for ( i = 0; i < nx; i++ )
01745             {
01746                 tx = plP_w3wcx( x[i], y[ny - 1], zmin );
01747                 ty = plP_w3wcy( x[i], y[ny - 1], zmin );
01748                 ux = plP_w3wcx( x[i], y[ny - 1], getz( zp, i, ny - 1 ) );
01749                 uy = plP_w3wcy( x[i], y[ny - 1], getz( zp, i, ny - 1 ) );
01750                 pljoin( tx, ty, ux, uy );
01751             }
01752         }
01753         if ( opt != 2 )
01754         {
01755             for ( i = 0; i < ny; i++ )
01756             {
01757                 tx = plP_w3wcx( x[0], y[i], zmin );
01758                 ty = plP_w3wcy( x[0], y[i], zmin );
01759                 ux = plP_w3wcx( x[0], y[i], getz( zp, 0, i ) );
01760                 uy = plP_w3wcy( x[0], y[i], getz( zp, 0, i ) );
01761                 pljoin( tx, ty, ux, uy );
01762             }
01763         }
01764     }
01765     else if ( cxx <= 0.0 && cxy >= 0.0 )
01766     {
01767         if ( opt != 1 )
01768         {
01769             for ( i = 0; i < nx; i++ )
01770             {
01771                 tx = plP_w3wcx( x[i], y[ny - 1], zmin );
01772                 ty = plP_w3wcy( x[i], y[ny - 1], zmin );
01773                 ux = plP_w3wcx( x[i], y[ny - 1], getz( zp, i, ny - 1 ) );
01774                 uy = plP_w3wcy( x[i], y[ny - 1], getz( zp, i, ny - 1 ) );
01775                 pljoin( tx, ty, ux, uy );
01776             }
01777         }
01778         if ( opt != 2 )
01779         {
01780             for ( i = 0; i < ny; i++ )
01781             {
01782                 tx = plP_w3wcx( x[nx - 1], y[i], zmin );
01783                 ty = plP_w3wcy( x[nx - 1], y[i], zmin );
01784                 ux = plP_w3wcx( x[nx - 1], y[i], getz( zp, nx - 1, i ) );
01785                 uy = plP_w3wcy( x[nx - 1], y[i], getz( zp, nx - 1, i ) );
01786                 pljoin( tx, ty, ux, uy );
01787             }
01788         }
01789     }
01790     else if ( cxx >= 0.0 && cxy >= 0.0 )
01791     {
01792         if ( opt != 1 )
01793         {
01794             for ( i = 0; i < nx; i++ )
01795             {
01796                 tx = plP_w3wcx( x[i], y[0], zmin );
01797                 ty = plP_w3wcy( x[i], y[0], zmin );
01798                 ux = plP_w3wcx( x[i], y[0], getz( zp, i, 0 ) );
01799                 uy = plP_w3wcy( x[i], y[0], getz( zp, i, 0 ) );
01800                 pljoin( tx, ty, ux, uy );
01801             }
01802         }
01803         if ( opt != 2 )
01804         {
01805             for ( i = 0; i < ny; i++ )
01806             {
01807                 tx = plP_w3wcx( x[nx - 1], y[i], zmin );
01808                 ty = plP_w3wcy( x[nx - 1], y[i], zmin );
01809                 ux = plP_w3wcx( x[nx - 1], y[i], getz( zp, nx - 1, i ) );
01810                 uy = plP_w3wcy( x[nx - 1], y[i], getz( zp, nx - 1, i ) );
01811                 pljoin( tx, ty, ux, uy );
01812             }
01813         }
01814     }
01815 }
01816 
01817 //--------------------------------------------------------------------------
01818 // void plgrid3()
01819 //
01820 // This routine draws a grid around the back side of the 3d plot with
01821 // hidden line removal.
01822 //--------------------------------------------------------------------------
01823 
01824 static void
01825 plgrid3( PLFLT tick )
01826 {
01827     PLFLT xmin, ymin, zmin, xmax, ymax, zmax, zscale;
01828     PLFLT cxx, cxy, cyx, cyy, cyz, zmin_in, zmax_in;
01829     PLINT u[3], v[3];
01830     PLINT nsub = 0;
01831     PLFLT tp;
01832 
01833     plP_gw3wc( &cxx, &cxy, &cyx, &cyy, &cyz );
01834     plP_gdom( &xmin, &xmax, &ymin, &ymax );
01835     plP_grange( &zscale, &zmin_in, &zmax_in );
01836     zmin = ( zmax_in > zmin_in ) ? zmin_in : zmax_in;
01837     zmax = ( zmax_in > zmin_in ) ? zmax_in : zmin_in;
01838 
01839     pldtik( zmin, zmax, &tick, &nsub, FALSE );
01840     tp     = tick * floor( zmin / tick ) + tick;
01841     pl3upv = 0;
01842 
01843     if ( cxx >= 0.0 && cxy <= 0.0 )
01844     {
01845         while ( tp <= zmax )
01846         {
01847             u[0] = plP_wcpcx( plP_w3wcx( xmin, ymax, tp ) );
01848             v[0] = plP_wcpcy( plP_w3wcy( xmin, ymax, tp ) );
01849             u[1] = plP_wcpcx( plP_w3wcx( xmax, ymax, tp ) );
01850             v[1] = plP_wcpcy( plP_w3wcy( xmax, ymax, tp ) );
01851             u[2] = plP_wcpcx( plP_w3wcx( xmax, ymin, tp ) );
01852             v[2] = plP_wcpcy( plP_w3wcy( xmax, ymin, tp ) );
01853             plnxtv( u, v, 0, 3, 0 );
01854 
01855             tp += tick;
01856         }
01857         u[0] = plP_wcpcx( plP_w3wcx( xmax, ymax, zmin ) );
01858         v[0] = plP_wcpcy( plP_w3wcy( xmax, ymax, zmin ) );
01859         u[1] = plP_wcpcx( plP_w3wcx( xmax, ymax, zmax ) );
01860         v[1] = plP_wcpcy( plP_w3wcy( xmax, ymax, zmax ) );
01861         plnxtv( u, v, 0, 2, 0 );
01862     }
01863     else if ( cxx <= 0.0 && cxy <= 0.0 )
01864     {
01865         while ( tp <= zmax )
01866         {
01867             u[0] = plP_wcpcx( plP_w3wcx( xmax, ymax, tp ) );
01868             v[0] = plP_wcpcy( plP_w3wcy( xmax, ymax, tp ) );
01869             u[1] = plP_wcpcx( plP_w3wcx( xmax, ymin, tp ) );
01870             v[1] = plP_wcpcy( plP_w3wcy( xmax, ymin, tp ) );
01871             u[2] = plP_wcpcx( plP_w3wcx( xmin, ymin, tp ) );
01872             v[2] = plP_wcpcy( plP_w3wcy( xmin, ymin, tp ) );
01873             plnxtv( u, v, 0, 3, 0 );
01874 
01875             tp += tick;
01876         }
01877         u[0] = plP_wcpcx( plP_w3wcx( xmax, ymin, zmin ) );
01878         v[0] = plP_wcpcy( plP_w3wcy( xmax, ymin, zmin ) );
01879         u[1] = plP_wcpcx( plP_w3wcx( xmax, ymin, zmax ) );
01880         v[1] = plP_wcpcy( plP_w3wcy( xmax, ymin, zmax ) );
01881         plnxtv( u, v, 0, 2, 0 );
01882     }
01883     else if ( cxx <= 0.0 && cxy >= 0.0 )
01884     {
01885         while ( tp <= zmax )
01886         {
01887             u[0] = plP_wcpcx( plP_w3wcx( xmax, ymin, tp ) );
01888             v[0] = plP_wcpcy( plP_w3wcy( xmax, ymin, tp ) );
01889             u[1] = plP_wcpcx( plP_w3wcx( xmin, ymin, tp ) );
01890             v[1] = plP_wcpcy( plP_w3wcy( xmin, ymin, tp ) );
01891             u[2] = plP_wcpcx( plP_w3wcx( xmin, ymax, tp ) );
01892             v[2] = plP_wcpcy( plP_w3wcy( xmin, ymax, tp ) );
01893             plnxtv( u, v, 0, 3, 0 );
01894 
01895             tp += tick;
01896         }
01897         u[0] = plP_wcpcx( plP_w3wcx( xmin, ymin, zmin ) );
01898         v[0] = plP_wcpcy( plP_w3wcy( xmin, ymin, zmin ) );
01899         u[1] = plP_wcpcx( plP_w3wcx( xmin, ymin, zmax ) );
01900         v[1] = plP_wcpcy( plP_w3wcy( xmin, ymin, zmax ) );
01901         plnxtv( u, v, 0, 2, 0 );
01902     }
01903     else if ( cxx >= 0.0 && cxy >= 0.0 )
01904     {
01905         while ( tp <= zmax )
01906         {
01907             u[0] = plP_wcpcx( plP_w3wcx( xmin, ymin, tp ) );
01908             v[0] = plP_wcpcy( plP_w3wcy( xmin, ymin, tp ) );
01909             u[1] = plP_wcpcx( plP_w3wcx( xmin, ymax, tp ) );
01910             v[1] = plP_wcpcy( plP_w3wcy( xmin, ymax, tp ) );
01911             u[2] = plP_wcpcx( plP_w3wcx( xmax, ymax, tp ) );
01912             v[2] = plP_wcpcy( plP_w3wcy( xmax, ymax, tp ) );
01913             plnxtv( u, v, 0, 3, 0 );
01914 
01915             tp += tick;
01916         }
01917         u[0] = plP_wcpcx( plP_w3wcx( xmin, ymax, zmin ) );
01918         v[0] = plP_wcpcy( plP_w3wcy( xmin, ymax, zmin ) );
01919         u[1] = plP_wcpcx( plP_w3wcx( xmin, ymax, zmax ) );
01920         v[1] = plP_wcpcy( plP_w3wcy( xmin, ymax, zmax ) );
01921         plnxtv( u, v, 0, 2, 0 );
01922     }
01923     pl3upv = 1;
01924 }
01925 
01926 //--------------------------------------------------------------------------
01927 // void plnxtv()
01928 //
01929 // Draw the next view of a 3-d plot. The physical coordinates of the
01930 // points for the next view are placed in the n points of arrays u and
01931 // v. The silhouette found so far is stored in the heap as a set of m peak
01932 // points.
01933 //
01934 // These routines dynamically allocate memory for hidden line removal.
01935 // Memory is allocated in blocks of 2*BINC*sizeof(PLINT) bytes.  Large
01936 // values of BINC give better performance but also allocate more memory
01937 // than is needed. If your 3D plots are very "spiky" or you are working
01938 // with very large matrices then you will probably want to increase BINC.
01939 //--------------------------------------------------------------------------
01940 
01941 static void
01942 plnxtv( PLINT *u, PLINT *v, PLFLT* c, PLINT n, PLINT init )
01943 {
01944     plnxtvhi( u, v, c, n, init );
01945 
01946     if ( pl3mode )
01947         plnxtvlo( u, v, c, n, init );
01948 }
01949 
01950 //--------------------------------------------------------------------------
01951 // void plnxtvhi()
01952 //
01953 // Draw the top side of the 3-d plot.
01954 //--------------------------------------------------------------------------
01955 
01956 static void
01957 plnxtvhi( PLINT *u, PLINT *v, PLFLT* c, PLINT n, PLINT init )
01958 {
01959     //
01960     // For the initial set of points, just display them and store them as the
01961     // peak points.
01962     //
01963     if ( init == 1 )
01964     {
01965         int i;
01966         oldhiview = (PLINT *) malloc( (size_t) ( 2 * n ) * sizeof ( PLINT ) );
01967         if ( !oldhiview )
01968             myexit( "plnxtvhi: Out of memory." );
01969 
01970         oldhiview[0] = u[0];
01971         oldhiview[1] = v[0];
01972         plP_draw3d( u[0], v[0], c, 0, 1 );
01973         for ( i = 1; i < n; i++ )
01974         {
01975             oldhiview[2 * i]     = u[i];
01976             oldhiview[2 * i + 1] = v[i];
01977             plP_draw3d( u[i], v[i], c, i, 0 );
01978         }
01979         mhi = n;
01980         return;
01981     }
01982 
01983     //
01984     // Otherwise, we need to consider hidden-line removal problem. We scan
01985     // over the points in both the old (i.e. oldhiview[]) and new (i.e. u[]
01986     // and v[]) arrays in order of increasing x coordinate.  At each stage, we
01987     // find the line segment in the other array (if one exists) that straddles
01988     // the x coordinate of the point. We have to determine if the point lies
01989     // above or below the line segment, and to check if the below/above status
01990     // has changed since the last point.
01991     //
01992     // If pl3upv = 0 we do not update the view, this is useful for drawing
01993     // lines on the graph after we are done plotting points.  Hidden line
01994     // removal is still done, but the view is not updated.
01995     //
01996     xxhi = 0;
01997     if ( pl3upv != 0 )
01998     {
01999         newhisize = 2 * ( mhi + BINC );
02000         if ( newhiview != NULL )
02001         {
02002             newhiview =
02003                 (PLINT *) realloc( (void *) newhiview,
02004                     (size_t) newhisize * sizeof ( PLINT ) );
02005         }
02006         else
02007         {
02008             newhiview =
02009                 (PLINT *) malloc( (size_t) newhisize * sizeof ( PLINT ) );
02010         }
02011         if ( !newhiview )
02012             myexit( "plnxtvhi: Out of memory." );
02013     }
02014 
02015     // Do the draw or shading with hidden line removal
02016 
02017     plnxtvhi_draw( u, v, c, n );
02018 
02019     // Set oldhiview
02020 
02021     swaphiview();
02022 }
02023 
02024 //--------------------------------------------------------------------------
02025 // void plnxtvhi_draw()
02026 //
02027 // Draw the top side of the 3-d plot.
02028 //--------------------------------------------------------------------------
02029 
02030 static void
02031 plnxtvhi_draw( PLINT *u, PLINT *v, PLFLT* c, PLINT n )
02032 {
02033     PLINT i   = 0, j = 0, first = 1;
02034     PLINT sx1 = 0, sx2 = 0, sy1 = 0, sy2 = 0;
02035     PLINT su1, su2, sv1, sv2;
02036     PLINT cx, cy, px, py;
02037     PLINT seg, ptold, lstold = 0, pthi, pnewhi = 0, newhi, change, ochange = 0;
02038 
02039 //
02040 // (oldhiview[2*i], oldhiview[2*i]) is the i'th point in the old array
02041 // (u[j], v[j]) is the j'th point in the new array
02042 //
02043 
02044 //
02045 // First attempt at 3d shading.  It works ok for simple plots, but
02046 // will just not draw faces, or draw them overlapping for very
02047 // jagged plots
02048 //
02049 
02050     while ( i < mhi || j < n )
02051     {
02052         //
02053         // The coordinates of the point under consideration are (px,py).  The
02054         // line segment joins (sx1,sy1) to (sx2,sy2).  "ptold" is true if the
02055         // point lies in the old array. We set it by comparing the x coordinates
02056         // of the i'th old point and the j'th new point, being careful if we
02057         // have fallen past the edges. Having found the point, load up the point
02058         // and segment coordinates appropriately.
02059         //
02060 
02061         ptold = ( j >= n || ( i < mhi && oldhiview[2 * i] < u[j] ) );
02062         if ( ptold )
02063         {
02064             px  = oldhiview[2 * i];
02065             py  = oldhiview[2 * i + 1];
02066             seg = j > 0 && j < n;
02067             if ( seg )
02068             {
02069                 sx1 = u[j - 1];
02070                 sy1 = v[j - 1];
02071                 sx2 = u[j];
02072                 sy2 = v[j];
02073             }
02074         }
02075         else
02076         {
02077             px  = u[j];
02078             py  = v[j];
02079             seg = i > 0 && i < mhi;
02080             if ( seg )
02081             {
02082                 sx1 = oldhiview[2 * ( i - 1 )];
02083                 sy1 = oldhiview[2 * ( i - 1 ) + 1];
02084                 sx2 = oldhiview[2 * i];
02085                 sy2 = oldhiview[2 * i + 1];
02086             }
02087         }
02088 
02089         //
02090         // Now determine if the point is higher than the segment, using the
02091         // logical function "above". We also need to know if it is the old view
02092         // or the new view that is higher. "newhi" is set true if the new view
02093         // is higher than the old.
02094         //
02095         if ( seg )
02096             pthi = plabv( px, py, sx1, sy1, sx2, sy2 );
02097         else
02098             pthi = 1;
02099 
02100         newhi = ( ptold && !pthi ) || ( !ptold && pthi );
02101         //
02102         // The last point and this point lie on different sides of
02103         // the current silouette
02104         //
02105         change = ( newhi && !pnewhi ) || ( !newhi && pnewhi );
02106 
02107         //
02108         // There is a new intersection point to put in the peak array if the
02109         // state of "newhi" changes.
02110         //
02111         if ( first )
02112         {
02113             plP_draw3d( px, py, c, j, 1 );
02114             first  = 0;
02115             lstold = ptold;
02116             savehipoint( px, py );
02117             pthi    = 0;
02118             ochange = 0;
02119         }
02120         else if ( change )
02121         {
02122             //
02123             // Take care of special cases at end of arrays.  If pl3upv is 0 the
02124             // endpoints are not connected to the old view.
02125             //
02126             if ( pl3upv == 0 && ( ( !ptold && j == 0 ) || ( ptold && i == 0 ) ) )
02127             {
02128                 plP_draw3d( px, py, c, j, 1 );
02129                 lstold  = ptold;
02130                 pthi    = 0;
02131                 ochange = 0;
02132             }
02133             else if ( pl3upv == 0 &&
02134                       ( ( !ptold && i >= mhi ) || ( ptold && j >= n ) ) )
02135             {
02136                 plP_draw3d( px, py, c, j, 1 );
02137                 lstold  = ptold;
02138                 pthi    = 0;
02139                 ochange = 0;
02140             }
02141             else
02142             {
02143                 //
02144                 // If pl3upv is not 0 then we do want to connect the current line
02145                 // with the previous view at the endpoints.  Also find intersection
02146                 // point with old view.
02147                 //
02148                 if ( i == 0 )
02149                 {
02150                     sx1 = oldhiview[0];
02151                     sy1 = -1;
02152                     sx2 = oldhiview[0];
02153                     sy2 = oldhiview[1];
02154                 }
02155                 else if ( i >= mhi )
02156                 {
02157                     sx1 = oldhiview[2 * ( mhi - 1 )];
02158                     sy1 = oldhiview[2 * ( mhi - 1 ) + 1];
02159                     sx2 = oldhiview[2 * ( mhi - 1 )];
02160                     sy2 = -1;
02161                 }
02162                 else
02163                 {
02164                     sx1 = oldhiview[2 * ( i - 1 )];
02165                     sy1 = oldhiview[2 * ( i - 1 ) + 1];
02166                     sx2 = oldhiview[2 * i];
02167                     sy2 = oldhiview[2 * i + 1];
02168                 }
02169 
02170                 if ( j == 0 )
02171                 {
02172                     su1 = u[0];
02173                     sv1 = -1;
02174                     su2 = u[0];
02175                     sv2 = v[0];
02176                 }
02177                 else if ( j >= n )
02178                 {
02179                     su1 = u[n - 1];
02180                     sv1 = v[n - 1];
02181                     su2 = u[n - 1];
02182                     sv2 = -1;
02183                 }
02184                 else
02185                 {
02186                     su1 = u[j - 1];
02187                     sv1 = v[j - 1];
02188                     su2 = u[j];
02189                     sv2 = v[j];
02190                 }
02191 
02192                 // Determine the intersection
02193 
02194                 pl3cut( sx1, sy1, sx2, sy2, su1, sv1, su2, sv2, &cx, &cy );
02195                 if ( cx == px && cy == py )
02196                 {
02197                     if ( lstold && !ochange )
02198                         plP_draw3d( px, py, c, j, 1 );
02199                     else
02200                         plP_draw3d( px, py, c, j, 0 );
02201 
02202                     savehipoint( px, py );
02203                     lstold = 1;
02204                     pthi   = 0;
02205                 }
02206                 else
02207                 {
02208                     if ( lstold && !ochange )
02209                         plP_draw3d( cx, cy, c, j, 1 );
02210                     else
02211                         plP_draw3d( cx, cy, c, j, 0 );
02212 
02213                     lstold = 1;
02214                     savehipoint( cx, cy );
02215                 }
02216                 ochange = 1;
02217             }
02218         }
02219 
02220         // If point is high then draw plot to point and update view.
02221 
02222         if ( pthi )
02223         {
02224             if ( lstold && ptold )
02225                 plP_draw3d( px, py, c, j, 1 );
02226             else
02227                 plP_draw3d( px, py, c, j, 0 );
02228 
02229             savehipoint( px, py );
02230             lstold  = ptold;
02231             ochange = 0;
02232         }
02233         pnewhi = newhi;
02234 
02235         if ( ptold )
02236             i++;
02237         else
02238             j++;
02239     }
02240 }
02241 
02242 //--------------------------------------------------------------------------
02243 // void  plP_draw3d()
02244 //
02245 // Does a simple move or line draw.
02246 //--------------------------------------------------------------------------
02247 
02248 static void
02249 plP_draw3d( PLINT x, PLINT y, PLFLT *c, PLINT j, PLINT move )
02250 {
02251     if ( move )
02252         plP_movphy( x, y );
02253     else
02254     {
02255         if ( c != NULL )
02256             plcol1( c[j - 1] );
02257         plP_draphy( x, y );
02258     }
02259 }
02260 
02261 //--------------------------------------------------------------------------
02262 // void plnxtvlo()
02263 //
02264 // Draw the bottom side of the 3-d plot.
02265 //--------------------------------------------------------------------------
02266 
02267 static void
02268 plnxtvlo( PLINT *u, PLINT *v, PLFLT*c, PLINT n, PLINT init )
02269 {
02270     PLINT i, j, first;
02271     PLINT sx1 = 0, sx2 = 0, sy1 = 0, sy2 = 0;
02272     PLINT su1, su2, sv1, sv2;
02273     PLINT cx, cy, px, py;
02274     PLINT seg, ptold, lstold = 0, ptlo, pnewlo, newlo, change, ochange = 0;
02275 
02276     first  = 1;
02277     pnewlo = 0;
02278 
02279     //
02280     // For the initial set of points, just display them and store them as the
02281     // peak points.
02282     //
02283     if ( init == 1 )
02284     {
02285         oldloview = (PLINT *) malloc( (size_t) ( 2 * n ) * sizeof ( PLINT ) );
02286         if ( !oldloview )
02287             myexit( "\nplnxtvlo: Out of memory." );
02288 
02289         plP_draw3d( u[0], v[0], c, 0, 1 );
02290         oldloview[0] = u[0];
02291         oldloview[1] = v[0];
02292         for ( i = 1; i < n; i++ )
02293         {
02294             plP_draw3d( u[i], v[i], c, i, 0 );
02295             oldloview[2 * i]     = u[i];
02296             oldloview[2 * i + 1] = v[i];
02297         }
02298         mlo = n;
02299         return;
02300     }
02301 
02302     //
02303     // Otherwise, we need to consider hidden-line removal problem. We scan
02304     // over the points in both the old (i.e. oldloview[]) and new (i.e. u[]
02305     // and v[]) arrays in order of increasing x coordinate.  At each stage, we
02306     // find the line segment in the other array (if one exists) that straddles
02307     // the x coordinate of the point. We have to determine if the point lies
02308     // above or below the line segment, and to check if the below/above status
02309     // has changed since the last point.
02310     //
02311     // If pl3upv = 0 we do not update the view, this is useful for drawing
02312     // lines on the graph after we are done plotting points.  Hidden line
02313     // removal is still done, but the view is not updated.
02314     //
02315     xxlo = 0;
02316     i    = 0;
02317     j    = 0;
02318     if ( pl3upv != 0 )
02319     {
02320         newlosize = 2 * ( mlo + BINC );
02321         if ( newloview != NULL )
02322         {
02323             newloview =
02324                 (PLINT *) realloc( (void *) newloview,
02325                     (size_t) newlosize * sizeof ( PLINT ) );
02326         }
02327         else
02328         {
02329             newloview =
02330                 (PLINT *) malloc( (size_t) newlosize * sizeof ( PLINT ) );
02331         }
02332         if ( !newloview )
02333             myexit( "plnxtvlo: Out of memory." );
02334     }
02335 
02336     //
02337     // (oldloview[2*i], oldloview[2*i]) is the i'th point in the old array
02338     // (u[j], v[j]) is the j'th point in the new array.
02339     //
02340     while ( i < mlo || j < n )
02341     {
02342         //
02343         // The coordinates of the point under consideration are (px,py).  The
02344         // line segment joins (sx1,sy1) to (sx2,sy2).  "ptold" is true if the
02345         // point lies in the old array. We set it by comparing the x coordinates
02346         // of the i'th old point and the j'th new point, being careful if we
02347         // have fallen past the edges. Having found the point, load up the point
02348         // and segment coordinates appropriately.
02349         //
02350         ptold = ( j >= n || ( i < mlo && oldloview[2 * i] < u[j] ) );
02351         if ( ptold )
02352         {
02353             px  = oldloview[2 * i];
02354             py  = oldloview[2 * i + 1];
02355             seg = j > 0 && j < n;
02356             if ( seg )
02357             {
02358                 sx1 = u[j - 1];
02359                 sy1 = v[j - 1];
02360                 sx2 = u[j];
02361                 sy2 = v[j];
02362             }
02363         }
02364         else
02365         {
02366             px  = u[j];
02367             py  = v[j];
02368             seg = i > 0 && i < mlo;
02369             if ( seg )
02370             {
02371                 sx1 = oldloview[2 * ( i - 1 )];
02372                 sy1 = oldloview[2 * ( i - 1 ) + 1];
02373                 sx2 = oldloview[2 * i];
02374                 sy2 = oldloview[2 * i + 1];
02375             }
02376         }
02377 
02378         //
02379         // Now determine if the point is lower than the segment, using the
02380         // logical function "above". We also need to know if it is the old view
02381         // or the new view that is lower. "newlo" is set true if the new view is
02382         // lower than the old.
02383         //
02384         if ( seg )
02385             ptlo = !plabv( px, py, sx1, sy1, sx2, sy2 );
02386         else
02387             ptlo = 1;
02388 
02389         newlo  = ( ptold && !ptlo ) || ( !ptold && ptlo );
02390         change = ( newlo && !pnewlo ) || ( !newlo && pnewlo );
02391 
02392         //
02393         // There is a new intersection point to put in the peak array if the
02394         // state of "newlo" changes.
02395         //
02396         if ( first )
02397         {
02398             plP_draw3d( px, py, c, j, 1 );
02399             first  = 0;
02400             lstold = ptold;
02401             savelopoint( px, py );
02402             ptlo    = 0;
02403             ochange = 0;
02404         }
02405         else if ( change )
02406         {
02407             //
02408             // Take care of special cases at end of arrays.  If pl3upv is 0 the
02409             // endpoints are not connected to the old view.
02410             //
02411             if ( pl3upv == 0 && ( ( !ptold && j == 0 ) || ( ptold && i == 0 ) ) )
02412             {
02413                 plP_draw3d( px, py, c, j, 1 );
02414                 lstold  = ptold;
02415                 ptlo    = 0;
02416                 ochange = 0;
02417             }
02418             else if ( pl3upv == 0 &&
02419                       ( ( !ptold && i >= mlo ) || ( ptold && j >= n ) ) )
02420             {
02421                 plP_draw3d( px, py, c, j, 1 );
02422                 lstold  = ptold;
02423                 ptlo    = 0;
02424                 ochange = 0;
02425             }
02426 
02427             //
02428             // If pl3upv is not 0 then we do want to connect the current line
02429             // with the previous view at the endpoints.  Also find intersection
02430             // point with old view.
02431             //
02432             else
02433             {
02434                 if ( i == 0 )
02435                 {
02436                     sx1 = oldloview[0];
02437                     sy1 = 100000;
02438                     sx2 = oldloview[0];
02439                     sy2 = oldloview[1];
02440                 }
02441                 else if ( i >= mlo )
02442                 {
02443                     sx1 = oldloview[2 * ( mlo - 1 )];
02444                     sy1 = oldloview[2 * ( mlo - 1 ) + 1];
02445                     sx2 = oldloview[2 * ( mlo - 1 )];
02446                     sy2 = 100000;
02447                 }
02448                 else
02449                 {
02450                     sx1 = oldloview[2 * ( i - 1 )];
02451                     sy1 = oldloview[2 * ( i - 1 ) + 1];
02452                     sx2 = oldloview[2 * i];
02453                     sy2 = oldloview[2 * i + 1];
02454                 }
02455 
02456                 if ( j == 0 )
02457                 {
02458                     su1 = u[0];
02459                     sv1 = 100000;
02460                     su2 = u[0];
02461                     sv2 = v[0];
02462                 }
02463                 else if ( j >= n )
02464                 {
02465                     su1 = u[n - 1];
02466                     sv1 = v[n - 1];
02467                     su2 = u[n];
02468                     sv2 = 100000;
02469                 }
02470                 else
02471                 {
02472                     su1 = u[j - 1];
02473                     sv1 = v[j - 1];
02474                     su2 = u[j];
02475                     sv2 = v[j];
02476                 }
02477 
02478                 // Determine the intersection
02479 
02480                 pl3cut( sx1, sy1, sx2, sy2, su1, sv1, su2, sv2, &cx, &cy );
02481                 if ( cx == px && cy == py )
02482                 {
02483                     if ( lstold && !ochange )
02484                         plP_draw3d( px, py, c, j, 1 );
02485                     else
02486                         plP_draw3d( px, py, c, j, 0 );
02487 
02488                     savelopoint( px, py );
02489                     lstold = 1;
02490                     ptlo   = 0;
02491                 }
02492                 else
02493                 {
02494                     if ( lstold && !ochange )
02495                         plP_draw3d( cx, cy, c, j, 1 );
02496                     else
02497                         plP_draw3d( cx, cy, c, j, 0 );
02498 
02499                     lstold = 1;
02500                     savelopoint( cx, cy );
02501                 }
02502                 ochange = 1;
02503             }
02504         }
02505 
02506         // If point is low then draw plot to point and update view.
02507 
02508         if ( ptlo )
02509         {
02510             if ( lstold && ptold )
02511                 plP_draw3d( px, py, c, j, 1 );
02512             else
02513                 plP_draw3d( px, py, c, j, 0 );
02514 
02515             savelopoint( px, py );
02516             lstold  = ptold;
02517             ochange = 0;
02518         }
02519 
02520         pnewlo = newlo;
02521 
02522         if ( ptold )
02523             i = i + 1;
02524         else
02525             j = j + 1;
02526     }
02527 
02528     // Set oldloview
02529 
02530     swaploview();
02531 }
02532 
02533 //--------------------------------------------------------------------------
02534 // savehipoint
02535 // savelopoint
02536 //
02537 // Add a point to the list of currently visible peaks/valleys, when
02538 // drawing the top/bottom surface, respectively.
02539 //--------------------------------------------------------------------------
02540 
02541 static void
02542 savehipoint( PLINT px, PLINT py )
02543 {
02544     if ( pl3upv == 0 )
02545         return;
02546 
02547     if ( xxhi >= newhisize )      // allocate additional space
02548     {
02549         newhisize += 2 * BINC;
02550         newhiview  = (PLINT *) realloc( (void *) newhiview,
02551             (size_t) newhisize * sizeof ( PLINT ) );
02552         if ( !newhiview )
02553             myexit( "savehipoint: Out of memory." );
02554     }
02555 
02556     newhiview[xxhi] = px;
02557     xxhi++;
02558     newhiview[xxhi] = py;
02559     xxhi++;
02560 }
02561 
02562 static void
02563 savelopoint( PLINT px, PLINT py )
02564 {
02565     if ( pl3upv == 0 )
02566         return;
02567 
02568     if ( xxlo >= newlosize )      // allocate additional space
02569     {
02570         newlosize += 2 * BINC;
02571         newloview  = (PLINT *) realloc( (void *) newloview,
02572             (size_t) newlosize * sizeof ( PLINT ) );
02573         if ( !newloview )
02574             myexit( "savelopoint: Out of memory." );
02575     }
02576 
02577     newloview[xxlo] = px;
02578     xxlo++;
02579     newloview[xxlo] = py;
02580     xxlo++;
02581 }
02582 
02583 //--------------------------------------------------------------------------
02584 // swaphiview
02585 // swaploview
02586 //
02587 // Swaps the top/bottom views.  Need to do a real swap so that the
02588 // memory cleanup routine really frees everything (and only once).
02589 //--------------------------------------------------------------------------
02590 
02591 static void
02592 swaphiview( void )
02593 {
02594     PLINT *tmp;
02595 
02596     if ( pl3upv != 0 )
02597     {
02598         mhi       = xxhi / 2;
02599         tmp       = oldhiview;
02600         oldhiview = newhiview;
02601         newhiview = tmp;
02602     }
02603 }
02604 
02605 static void
02606 swaploview( void )
02607 {
02608     PLINT *tmp;
02609 
02610     if ( pl3upv != 0 )
02611     {
02612         mlo       = xxlo / 2;
02613         tmp       = oldloview;
02614         oldloview = newloview;
02615         newloview = tmp;
02616     }
02617 }
02618 
02619 //--------------------------------------------------------------------------
02620 // freework
02621 //
02622 // Frees memory associated with work arrays
02623 //--------------------------------------------------------------------------
02624 
02625 static void
02626 freework( void )
02627 {
02628     free_mem( oldhiview );
02629     free_mem( oldloview );
02630     free_mem( newhiview );
02631     free_mem( newloview );
02632     free_mem( vtmp );
02633     free_mem( utmp );
02634     free_mem( ctmp );
02635 }
02636 
02637 //--------------------------------------------------------------------------
02638 // myexit
02639 //
02640 // Calls plexit, cleaning up first.
02641 //--------------------------------------------------------------------------
02642 
02643 static void
02644 myexit( const char *msg )
02645 {
02646     freework();
02647     plexit( msg );
02648 }
02649 
02650 //--------------------------------------------------------------------------
02651 // myabort
02652 //
02653 // Calls plabort, cleaning up first.
02654 // Caller should return to the user program.
02655 //--------------------------------------------------------------------------
02656 
02657 static void
02658 myabort( const char *msg )
02659 {
02660     freework();
02661     plabort( msg );
02662 }
02663 
02664 //--------------------------------------------------------------------------
02665 // int plabv()
02666 //
02667 // Determines if point (px,py) lies above the line joining (sx1,sy1) to
02668 // (sx2,sy2). It only works correctly if sx1 <= px <= sx2.
02669 //--------------------------------------------------------------------------
02670 
02671 static int
02672 plabv( PLINT px, PLINT py, PLINT sx1, PLINT sy1, PLINT sx2, PLINT sy2 )
02673 {
02674     int above;
02675 
02676     if ( py >= sy1 && py >= sy2 )
02677         above = 1;
02678     else if ( py < sy1 && py < sy2 )
02679         above = 0;
02680     else if ( (double) ( sx2 - sx1 ) * ( py - sy1 ) >=
02681               (double) ( px - sx1 ) * ( sy2 - sy1 ) )
02682         above = 1;
02683     else
02684         above = 0;
02685 
02686     return above;
02687 }
02688 
02689 //--------------------------------------------------------------------------
02690 // void pl3cut()
02691 //
02692 // Determines the point of intersection (cx,cy) between the line joining
02693 // (sx1,sy1) to (sx2,sy2) and the line joining (su1,sv1) to (su2,sv2).
02694 //--------------------------------------------------------------------------
02695 
02696 static void
02697 pl3cut( PLINT sx1, PLINT sy1, PLINT sx2, PLINT sy2,
02698         PLINT su1, PLINT sv1, PLINT su2, PLINT sv2, PLINT *cx, PLINT *cy )
02699 {
02700     PLINT x21, y21, u21, v21, yv1, xu1, a, b;
02701     double fa, fb;
02702 
02703     x21 = sx2 - sx1;
02704     y21 = sy2 - sy1;
02705     u21 = su2 - su1;
02706     v21 = sv2 - sv1;
02707     yv1 = sy1 - sv1;
02708     xu1 = sx1 - su1;
02709 
02710     a  = x21 * v21 - y21 * u21;
02711     fa = (double) a;
02712     if ( a == 0 )
02713     {
02714         if ( sx2 < su2 )
02715         {
02716             *cx = sx2;
02717             *cy = sy2;
02718         }
02719         else
02720         {
02721             *cx = su2;
02722             *cy = sv2;
02723         }
02724         return;
02725     }
02726     else
02727     {
02728         b   = yv1 * u21 - xu1 * v21;
02729         fb  = (double) b;
02730         *cx = (PLINT) ( sx1 + ( fb * x21 ) / fa + .5 );
02731         *cy = (PLINT) ( sy1 + ( fb * y21 ) / fa + .5 );
02732     }
02733 }
02734 
02735 //--------------------------------------------------------------------------
02736 // void plRotationShear
02737 //
02738 // Calculates the rotation and shear angles from a plplot transformation matrix
02739 //
02740 // N.B. the plot transformation matrix
02741 //
02742 // [xFormMatrix[0] xFormMatrix[2]]
02743 // [xFormMatrix[1] xFormMatrix[3]]
02744 //
02745 // is calculated as
02746 //
02747 // [stride cos(t)    stride sin(t)]
02748 // [sin(p-t)              cos(p-t)]
02749 //
02750 // where t is the rotation angle and p is the shear angle.
02751 // The purpose of this routine is to determine stride, rotation angle,
02752 // and shear angle from xFormMatrix.
02753 //
02754 // For information only, xFormMatrix is the product of the following
02755 // rotation, skew(shear), and scale matrices:
02756 //
02757 //  [stride    0] [1      0] [cos(t)  sin(t)]
02758 //  [0    cos(p)] [tan(p) 1] [-sin(t) cos(t)]
02759 //
02760 //--------------------------------------------------------------------------
02761 
02762 void
02763 plRotationShear( PLFLT *xFormMatrix, PLFLT *rotation, PLFLT *shear, PLFLT *stride )
02764 {
02765     PLFLT smr;
02766     *stride = sqrt( xFormMatrix[0] * xFormMatrix[0] + xFormMatrix[2] * xFormMatrix[2] );
02767 
02768     // Calculate rotation in range from -pi to pi.
02769     *rotation = atan2( xFormMatrix[2], xFormMatrix[0] );
02770 
02771     // Calculate shear - rotation in range from -pi to pi.
02772     smr = atan2( xFormMatrix[1], xFormMatrix[3] );
02773 
02774     // Calculate shear in range from -2 pi to 2 pi.
02775     *shear = smr + *rotation;
02776 
02777     // Calculate shear in range from -pi to pi.
02778     if ( *shear < -PI )
02779         *shear += 2. * PI;
02780     else if ( *shear > PI )
02781         *shear -= 2. * PI;
02782 
02783     // Actually must honour some convention to calculate the negative
02784     // of the shear angle instead of the shear angle. Why??
02785     *shear = -*shear;
02786     // Comment out the modified old logic which determines the negative
02787     // of the shear angle in a more complicated way.  Note, the modification
02788     // to divide the asin argument by *stride which solved a long-standing
02789     // bug (as does the above logic in a simpler way).
02790     //
02791     //shear = -asin( (xFormMatrix[0] * xFormMatrix[1] +
02792     //               xFormMatrix[2] * xFormMatrix[3] )/ *stride);
02793     //
02794 
02795     // Compute the cross product of the vectors [1,0] and [0,1] to
02796     // determine if we need to make a "quadrant 3,4" adjustment
02797     // to the shear angle.
02798 
02799     //
02800     // if ( xFormMatrix[0] * xFormMatrix[3] - xFormMatrix[1] * xFormMatrix[2] < 0.0 )
02801     // {
02802     //shear = -( M_PI + *shear );
02803     // }
02804     //
02805 }
02806 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Defines