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