PLplot 5.9.6
|
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