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