PLplot 5.9.6
|
00001 /* $Id$ 00002 * 00003 * Contour plotter. 00004 * 00005 * Copyright (C) 1995, 2000, 2001 Maurice LeBrun 00006 * Copyright (C) 2000, 2002 Joao Cardoso 00007 * Copyright (C) 2000, 2001, 2002, 2004 Alan W. Irwin 00008 * Copyright (C) 2004 Andrew Ross 00009 * 00010 * This file is part of PLplot. 00011 * 00012 * PLplot is free software; you can redistribute it and/or modify 00013 * it under the terms of the GNU General Library Public License as published 00014 * by the Free Software Foundation; either version 2 of the License, or 00015 * (at your option) any later version. 00016 * 00017 * PLplot is distributed in the hope that it will be useful, 00018 * but WITHOUT ANY WARRANTY; without even the implied warranty of 00019 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 00020 * GNU Library General Public License for more details. 00021 * 00022 * You should have received a copy of the GNU Library General Public License 00023 * along with PLplot; if not, write to the Free Software 00024 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA 00025 */ 00026 00027 #include "plplotP.h" 00028 00029 #ifdef MSDOS 00030 #pragma optimize("",off) 00031 #endif 00032 00033 /* Static function prototypes. */ 00034 00035 static void 00036 plcntr( PLFLT ( *plf2eval )( PLINT, PLINT, PLPointer ), 00037 PLPointer plf2eval_data, 00038 PLINT nx, PLINT ny, PLINT kx, PLINT lx, 00039 PLINT ky, PLINT ly, PLFLT flev, PLINT **ipts, 00040 void ( *pltr )( PLFLT, PLFLT, PLFLT *, PLFLT *, PLPointer ), 00041 PLPointer pltr_data ); 00042 00043 static void 00044 pldrawcn( PLFLT ( *plf2eval )( PLINT, PLINT, PLPointer ), 00045 PLPointer plf2eval_data, 00046 PLINT nx, PLINT ny, PLINT kx, PLINT lx, 00047 PLINT ky, PLINT ly, PLFLT flev, char *flabel, PLINT kcol, PLINT krow, 00048 PLFLT lastx, PLFLT lasty, PLINT startedge, 00049 PLINT **ipts, PLFLT *distance, PLINT *lastindex, 00050 void ( *pltr )( PLFLT, PLFLT, PLFLT *, PLFLT *, PLPointer ), 00051 PLPointer pltr_data ); 00052 00053 static void 00054 plfloatlabel( PLFLT value, char *string, PLINT len ); 00055 00056 static PLFLT 00057 plP_pcwcx( PLINT x ); 00058 00059 static PLFLT 00060 plP_pcwcy( PLINT y ); 00061 00062 static void 00063 pl_drawcontlabel( PLFLT tpx, PLFLT tpy, char *flabel, PLFLT *distance, PLINT *lastindex ); 00064 00065 /* Error flag for aborts */ 00066 00067 static int error; 00068 00069 /****************************************/ 00070 /* */ 00071 /* Defaults for contour label printing. */ 00072 /* */ 00073 /****************************************/ 00074 00075 /* Font height for contour labels (normalized) */ 00076 static PLFLT 00077 contlabel_size = 0.3; 00078 00079 /* Offset of label from contour line (if set to 0.0, labels are printed on the lines). */ 00080 static PLFLT 00081 contlabel_offset = 0.006; 00082 00083 /* Spacing parameter for contour labels */ 00084 static PLFLT 00085 contlabel_space = 0.1; 00086 00087 /* Activate labels, default off */ 00088 static PLINT 00089 contlabel_active = 0; 00090 00091 /* If the contour label exceed 10^(limexp) or 10^(-limexp), the exponential format is used */ 00092 static PLINT 00093 limexp = 4; 00094 00095 /* Number of significant digits */ 00096 static PLINT 00097 sigprec = 2; 00098 00099 /******** contour lines storage ****************************/ 00100 00101 static CONT_LEVEL *startlev = NULL; 00102 static CONT_LEVEL *currlev; 00103 static CONT_LINE *currline; 00104 00105 static int cont3d = 0; 00106 00107 static CONT_LINE * 00108 alloc_line( CONT_LEVEL *node ) 00109 { 00110 CONT_LINE *line; 00111 00112 if ( ( line = (CONT_LINE *) malloc( sizeof ( CONT_LINE ) ) ) == NULL ) 00113 { 00114 plexit( "alloc_line: Insufficient memory" ); 00115 } 00116 00117 line->x = (PLFLT *) malloc( LINE_ITEMS * sizeof ( PLFLT ) ); 00118 line->y = (PLFLT *) malloc( LINE_ITEMS * sizeof ( PLFLT ) ); 00119 00120 if ( ( line->x == NULL ) || ( line->y == NULL ) ) 00121 { 00122 plexit( "alloc_line: Insufficient memory" ); 00123 } 00124 00125 line->npts = 0; 00126 line->next = NULL; 00127 00128 return line; 00129 } 00130 00131 static CONT_LEVEL * 00132 alloc_level( PLFLT level ) 00133 { 00134 CONT_LEVEL *node; 00135 00136 if ( ( node = (CONT_LEVEL *) malloc( sizeof ( CONT_LEVEL ) ) ) == NULL ) 00137 { 00138 plexit( "alloc_level: Insufficient memory" ); 00139 } 00140 node->level = level; 00141 node->next = NULL; 00142 node->line = alloc_line( node ); 00143 00144 return node; 00145 } 00146 00147 static void 00148 realloc_line( CONT_LINE *line ) 00149 { 00150 if ( ( ( line->x = (PLFLT *) realloc( line->x, 00151 ( line->npts + LINE_ITEMS ) * sizeof ( PLFLT ) ) ) == NULL ) || 00152 ( ( line->y = (PLFLT *) realloc( line->y, 00153 ( line->npts + LINE_ITEMS ) * sizeof ( PLFLT ) ) ) == NULL ) ) 00154 plexit( "realloc_line: Insufficient memory" ); 00155 } 00156 00157 00158 /* new contour level */ 00159 static void 00160 cont_new_store( PLFLT level ) 00161 { 00162 if ( cont3d ) 00163 { 00164 if ( startlev == NULL ) 00165 { 00166 startlev = alloc_level( level ); 00167 currlev = startlev; 00168 } 00169 else 00170 { 00171 currlev->next = alloc_level( level ); 00172 currlev = currlev->next; 00173 } 00174 currline = currlev->line; 00175 } 00176 } 00177 00178 void 00179 cont_clean_store( CONT_LEVEL *ct ) 00180 { 00181 CONT_LINE *tline, *cline; 00182 CONT_LEVEL *tlev, *clevel; 00183 00184 if ( ct != NULL ) 00185 { 00186 clevel = ct; 00187 00188 do 00189 { 00190 cline = clevel->line; 00191 do 00192 { 00193 #ifdef CONT_PLOT_DEBUG /* for 2D plots. For 3D plots look at plot3.c:plotsh3di() */ 00194 plP_movwor( cline->x[0], cline->y[0] ); 00195 for ( j = 1; j < cline->npts; j++ ) 00196 plP_drawor( cline->x[j], cline->y[j] ); 00197 #endif 00198 tline = cline->next; 00199 free( cline->x ); 00200 free( cline->y ); 00201 free( cline ); 00202 cline = tline; 00203 } 00204 while ( cline != NULL ); 00205 tlev = clevel->next; 00206 free( clevel ); 00207 clevel = tlev; 00208 } 00209 while ( clevel != NULL ); 00210 startlev = NULL; 00211 } 00212 } 00213 00214 static void 00215 cont_xy_store( PLFLT xx, PLFLT yy ) 00216 { 00217 if ( cont3d ) 00218 { 00219 PLINT pts = currline->npts; 00220 00221 if ( pts % LINE_ITEMS == 0 ) 00222 realloc_line( currline ); 00223 00224 currline->x[pts] = xx; 00225 currline->y[pts] = yy; 00226 currline->npts++; 00227 } 00228 else 00229 plP_drawor( xx, yy ); 00230 } 00231 00232 static void 00233 cont_mv_store( PLFLT xx, PLFLT yy ) 00234 { 00235 if ( cont3d ) 00236 { 00237 if ( currline->npts != 0 ) /* not an empty list, allocate new */ 00238 { 00239 currline->next = alloc_line( currlev ); 00240 currline = currline->next; 00241 } 00242 00243 /* and fill first element */ 00244 currline->x[0] = xx; 00245 currline->y[0] = yy; 00246 currline->npts = 1; 00247 } 00248 else 00249 plP_movwor( xx, yy ); 00250 } 00251 00252 /* small routine to set offset and spacing of contour labels, see desciption above */ 00253 void c_pl_setcontlabelparam( PLFLT offset, PLFLT size, PLFLT spacing, PLINT active ) 00254 { 00255 contlabel_offset = offset; 00256 contlabel_size = size; 00257 contlabel_space = spacing; 00258 contlabel_active = active; 00259 } 00260 00261 /* small routine to set the format of the contour labels, description of limexp and prec see above */ 00262 void c_pl_setcontlabelformat( PLINT lexp, PLINT sigdig ) 00263 { 00264 limexp = lexp; 00265 sigprec = sigdig; 00266 } 00267 00268 static void pl_drawcontlabel( PLFLT tpx, PLFLT tpy, char *flabel, PLFLT *distance, PLINT *lastindex ) 00269 { 00270 PLFLT delta_x, delta_y; 00271 PLINT currx_old, curry_old; 00272 00273 delta_x = plP_pcdcx( plsc->currx ) - plP_pcdcx( plP_wcpcx( tpx ) ); 00274 delta_y = plP_pcdcy( plsc->curry ) - plP_pcdcy( plP_wcpcy( tpy ) ); 00275 00276 currx_old = plsc->currx; 00277 curry_old = plsc->curry; 00278 00279 *distance += sqrt( delta_x * delta_x + delta_y * delta_y ); 00280 00281 plP_drawor( tpx, tpy ); 00282 00283 if ( (int ) ( fabs( *distance / contlabel_space ) ) > *lastindex ) 00284 { 00285 PLFLT scale, vec_x, vec_y, mx, my, dev_x, dev_y, off_x, off_y; 00286 00287 vec_x = tpx - plP_pcwcx( currx_old ); 00288 vec_y = tpy - plP_pcwcy( curry_old ); 00289 00290 /* Ensure labels appear the right way up */ 00291 if ( vec_x < 0 ) 00292 { 00293 vec_x = -vec_x; 00294 vec_y = -vec_y; 00295 } 00296 00297 mx = (double ) plsc->wpxscl / (double ) plsc->phyxlen; 00298 my = (double ) plsc->wpyscl / (double ) plsc->phyylen; 00299 00300 dev_x = -my * vec_y / mx; 00301 dev_y = mx * vec_x / my; 00302 00303 scale = sqrt( ( mx * mx * dev_x * dev_x + my * my * dev_y * dev_y ) / 00304 ( contlabel_offset * contlabel_offset ) ); 00305 00306 off_x = dev_x / scale; 00307 off_y = dev_y / scale; 00308 00309 plptex( tpx + off_x, tpy + off_y, vec_x, vec_y, 0.5, flabel ); 00310 plP_movwor( tpx, tpy ); 00311 ( *lastindex )++; 00312 } 00313 else 00314 plP_movwor( tpx, tpy ); 00315 } 00316 00317 00318 /* Format contour labels. Arguments: 00319 * value: floating point number to be formatted 00320 * string: the formatted label, plptex must be called with it to actually 00321 * print the label 00322 */ 00323 00324 static void plfloatlabel( PLFLT value, char *string, PLINT len ) 00325 { 00326 PLINT setpre, precis; 00327 /* form[10] gives enough space for all non-malicious formats. 00328 * tmpstring[15] gives enough room for 3 digits in a negative exponent 00329 * or 4 digits in a positive exponent + null termination. That 00330 * should be enough for all non-malicious use. 00331 * Obviously there are security issues here that 00332 * should be addressed as well. 00333 */ 00334 #define FORM_LEN 10 00335 #define TMPSTRING_LEN 15 00336 char form[FORM_LEN], tmpstring[TMPSTRING_LEN]; 00337 PLINT exponent = 0; 00338 PLFLT mant, tmp; 00339 00340 PLINT prec = sigprec; 00341 00342 plP_gprec( &setpre, &precis ); 00343 00344 if ( setpre ) 00345 prec = precis; 00346 00347 if ( value > 0.0 ) 00348 tmp = log10( value ); 00349 else if ( value < 0.0 ) 00350 tmp = log10( -value ); 00351 else 00352 tmp = 0; 00353 00354 if ( tmp >= 0.0 ) 00355 exponent = (int ) tmp; 00356 else if ( tmp < 0.0 ) 00357 { 00358 tmp = -tmp; 00359 if ( floor( tmp ) < tmp ) 00360 exponent = -(int ) ( floor( tmp ) + 1.0 ); 00361 else 00362 exponent = -(int ) ( floor( tmp ) ); 00363 } 00364 00365 mant = value / pow( 10.0, exponent ); 00366 00367 if ( mant != 0.0 ) 00368 mant = (int ) ( mant * pow( 10.0, prec - 1 ) + 0.5 * mant / fabs( mant ) ) / pow( 10.0, prec - 1 ); 00369 00370 snprintf( form, FORM_LEN, "%%.%df", prec - 1 ); 00371 snprintf( string, len, form, mant ); 00372 snprintf( tmpstring, TMPSTRING_LEN, "#(229)10#u%d", exponent ); 00373 strncat( string, tmpstring, len - strlen( string ) - 1 ); 00374 00375 if ( abs( exponent ) < limexp || value == 0.0 ) 00376 { 00377 value = pow( 10.0, exponent ) * mant; 00378 00379 if ( exponent >= 0 ) 00380 prec = prec - 1 - exponent; 00381 else 00382 prec = prec - 1 + abs( exponent ); 00383 00384 if ( prec < 0 ) 00385 prec = 0; 00386 00387 snprintf( form, FORM_LEN, "%%.%df", (int) prec ); 00388 snprintf( string, len, form, value ); 00389 } 00390 } 00391 00392 /* physical coords (x) to world coords */ 00393 00394 static PLFLT 00395 plP_pcwcx( PLINT x ) 00396 { 00397 return ( ( x - plsc->wpxoff ) / plsc->wpxscl ); 00398 } 00399 00400 /* physical coords (y) to world coords */ 00401 00402 static PLFLT 00403 plP_pcwcy( PLINT y ) 00404 { 00405 return ( ( y - plsc->wpyoff ) / plsc->wpyscl ); 00406 } 00407 00408 /*--------------------------------------------------------------------------*\ 00409 * plf2eval1() 00410 * 00411 * Does a lookup from a 2d function array. Array is of type (PLFLT **), 00412 * and is column dominant (normal C ordering). 00413 \*--------------------------------------------------------------------------*/ 00414 00415 PLFLT 00416 plf2eval1( PLINT ix, PLINT iy, PLPointer plf2eval_data ) 00417 { 00418 PLFLT value; 00419 PLFLT **z = (PLFLT **) plf2eval_data; 00420 00421 value = z[ix][iy]; 00422 00423 return value; 00424 } 00425 00426 /*--------------------------------------------------------------------------*\ 00427 * plf2eval2() 00428 * 00429 * Does a lookup from a 2d function array. plf2eval_data is treated as type 00430 * (PLfGrid2 *). 00431 \*--------------------------------------------------------------------------*/ 00432 00433 PLFLT 00434 plf2eval2( PLINT ix, PLINT iy, PLPointer plf2eval_data ) 00435 { 00436 PLFLT value; 00437 PLfGrid2 *grid = (PLfGrid2 *) plf2eval_data; 00438 00439 value = grid->f[ix][iy]; 00440 00441 return value; 00442 } 00443 00444 /*--------------------------------------------------------------------------*\ 00445 * plf2eval() 00446 * 00447 * Does a lookup from a 2d function array. Array is of type (PLFLT *), and 00448 * is column dominant (normal C ordering). You MUST fill the ny maximum 00449 * array index entry in the PLfGrid struct. 00450 \*--------------------------------------------------------------------------*/ 00451 00452 PLFLT 00453 plf2eval( PLINT ix, PLINT iy, PLPointer plf2eval_data ) 00454 { 00455 PLFLT value; 00456 PLfGrid *grid = (PLfGrid *) plf2eval_data; 00457 00458 value = grid->f[ix * grid->ny + iy]; 00459 00460 return value; 00461 } 00462 00463 /*--------------------------------------------------------------------------*\ 00464 * plf2evalr() 00465 * 00466 * Does a lookup from a 2d function array. Array is of type (PLFLT *), and 00467 * is row dominant (Fortran ordering). You MUST fill the nx maximum array 00468 * index entry in the PLfGrid struct. 00469 \*--------------------------------------------------------------------------*/ 00470 00471 PLFLT 00472 plf2evalr( PLINT ix, PLINT iy, PLPointer plf2eval_data ) 00473 { 00474 PLFLT value; 00475 PLfGrid *grid = (PLfGrid *) plf2eval_data; 00476 00477 value = grid->f[ix + iy * grid->nx]; 00478 00479 return value; 00480 } 00481 00482 /*--------------------------------------------------------------------------*\ 00483 * 00484 * cont_store: 00485 * 00486 * Draw contour lines in memory. 00487 * cont_clean_store() must be called after use to release allocated memory. 00488 * 00489 \*--------------------------------------------------------------------------*/ 00490 00491 void 00492 cont_store( PLFLT **f, PLINT nx, PLINT ny, PLINT kx, PLINT lx, 00493 PLINT ky, PLINT ly, PLFLT *clevel, PLINT nlevel, 00494 void ( *pltr )( PLFLT, PLFLT, PLFLT *, PLFLT *, PLPointer ), 00495 PLPointer pltr_data, 00496 CONT_LEVEL **contour ) 00497 { 00498 cont3d = 1; 00499 00500 plcont( f, nx, ny, kx, lx, ky, ly, clevel, nlevel, 00501 pltr, pltr_data ); 00502 00503 *contour = startlev; 00504 cont3d = 0; 00505 } 00506 00507 /*--------------------------------------------------------------------------*\ 00508 * void plcont() 00509 * 00510 * Draws a contour plot from data in f(nx,ny). Is just a front-end to 00511 * plfcont, with a particular choice for f2eval and f2eval_data. 00512 \*--------------------------------------------------------------------------*/ 00513 00514 void 00515 c_plcont( PLFLT **f, PLINT nx, PLINT ny, PLINT kx, PLINT lx, 00516 PLINT ky, PLINT ly, PLFLT *clevel, PLINT nlevel, 00517 void ( *pltr )( PLFLT, PLFLT, PLFLT *, PLFLT *, PLPointer ), 00518 PLPointer pltr_data ) 00519 { 00520 if ( pltr == NULL ) 00521 { 00522 /* If pltr is undefined, abort with an error. */ 00523 plabort( "plcont: The pltr callback must be defined" ); 00524 return; 00525 } 00526 00527 plfcont( plf2eval1, ( PLPointer ) f, 00528 nx, ny, kx, lx, ky, ly, clevel, nlevel, 00529 pltr, pltr_data ); 00530 } 00531 00532 /*--------------------------------------------------------------------------*\ 00533 * void plfcont() 00534 * 00535 * Draws a contour plot using the function evaluator f2eval and data stored 00536 * by way of the f2eval_data pointer. This allows arbitrary organizations 00537 * of 2d array data to be used. 00538 * 00539 * The subrange of indices used for contouring is kx to lx in the x 00540 * direction and from ky to ly in the y direction. The array of contour 00541 * levels is clevel(nlevel), and "pltr" is the name of a function which 00542 * transforms array indicies into world coordinates. 00543 * 00544 * Note that the fortran-like minimum and maximum indices (kx, lx, ky, ly) 00545 * are translated into more C-like ones. I've only kept them as they are 00546 * for the plcontf() argument list because of backward compatibility. 00547 \*--------------------------------------------------------------------------*/ 00548 00549 void 00550 plfcont( PLFLT ( *f2eval )( PLINT, PLINT, PLPointer ), 00551 PLPointer f2eval_data, 00552 PLINT nx, PLINT ny, PLINT kx, PLINT lx, 00553 PLINT ky, PLINT ly, PLFLT *clevel, PLINT nlevel, 00554 void ( *pltr )( PLFLT, PLFLT, PLFLT *, PLFLT *, PLPointer ), 00555 PLPointer pltr_data ) 00556 { 00557 PLINT i, **ipts; 00558 00559 if ( kx < 1 || kx >= lx ) 00560 { 00561 plabort( "plfcont: indices must satisfy 1 <= kx <= lx <= nx" ); 00562 return; 00563 } 00564 if ( ky < 1 || ky >= ly ) 00565 { 00566 plabort( "plfcont: indices must satisfy 1 <= ky <= ly <= ny" ); 00567 return; 00568 } 00569 00570 if ( ( ipts = (PLINT **) malloc( nx * sizeof ( PLINT * ) ) ) == NULL ) 00571 { 00572 plexit( "plfcont: Insufficient memory" ); 00573 } 00574 00575 for ( i = 0; i < nx; i++ ) 00576 { 00577 if ( ( ipts[i] = (PLINT *) malloc( ny * sizeof ( PLINT * ) ) ) == NULL ) 00578 { 00579 plexit( "plfcont: Insufficient memory" ); 00580 } 00581 } 00582 00583 for ( i = 0; i < nlevel; i++ ) 00584 { 00585 plcntr( f2eval, f2eval_data, 00586 nx, ny, kx - 1, lx - 1, ky - 1, ly - 1, clevel[i], ipts, 00587 pltr, pltr_data ); 00588 00589 if ( error ) 00590 { 00591 error = 0; 00592 goto done; 00593 } 00594 } 00595 00596 done: 00597 for ( i = 0; i < nx; i++ ) 00598 { 00599 free( (void *) ipts[i] ); 00600 } 00601 free( (void *) ipts ); 00602 } 00603 00604 /*--------------------------------------------------------------------------*\ 00605 * void plcntr() 00606 * 00607 * The contour for a given level is drawn here. Note iscan has nx 00608 * elements. ixstor and iystor each have nstor elements. 00609 \*--------------------------------------------------------------------------*/ 00610 00611 static void 00612 plcntr( PLFLT ( *f2eval )( PLINT, PLINT, PLPointer ), 00613 PLPointer f2eval_data, 00614 PLINT nx, PLINT ny, PLINT kx, PLINT lx, 00615 PLINT ky, PLINT ly, PLFLT flev, PLINT **ipts, 00616 void ( *pltr )( PLFLT, PLFLT, PLFLT *, PLFLT *, PLPointer ), 00617 PLPointer pltr_data ) 00618 { 00619 PLINT kcol, krow, lastindex; 00620 PLFLT distance; 00621 PLFLT save_def, save_scale; 00622 00623 char flabel[30]; 00624 plgchr( &save_def, &save_scale ); 00625 save_scale = save_scale / save_def; 00626 00627 cont_new_store( flev ); 00628 00629 /* format contour label for plptex and define the font height of the labels */ 00630 plfloatlabel( flev, flabel, 30 ); 00631 plschr( 0.0, contlabel_size ); 00632 00633 /* Clear array for traversed squares */ 00634 for ( kcol = kx; kcol < lx; kcol++ ) 00635 { 00636 for ( krow = ky; krow < ly; krow++ ) 00637 { 00638 ipts[kcol][krow] = 0; 00639 } 00640 } 00641 00642 00643 for ( krow = ky; krow < ly; krow++ ) 00644 { 00645 for ( kcol = kx; kcol < lx; kcol++ ) 00646 { 00647 if ( ipts[kcol][krow] == 0 ) 00648 { 00649 /* Follow and draw a contour */ 00650 pldrawcn( f2eval, f2eval_data, 00651 nx, ny, kx, lx, ky, ly, flev, flabel, kcol, krow, 00652 0.0, 0.0, -2, ipts, &distance, &lastindex, 00653 pltr, pltr_data ); 00654 00655 if ( error ) 00656 return; 00657 } 00658 } 00659 } 00660 plschr( save_def, save_scale ); 00661 } 00662 00663 /*--------------------------------------------------------------------------*\ 00664 * void pldrawcn() 00665 * 00666 * Follow and draw a contour. 00667 \*--------------------------------------------------------------------------*/ 00668 00669 static void 00670 pldrawcn( PLFLT ( *f2eval )( PLINT, PLINT, PLPointer ), 00671 PLPointer f2eval_data, 00672 PLINT nx, PLINT ny, PLINT kx, PLINT lx, 00673 PLINT ky, PLINT ly, PLFLT flev, char *flabel, PLINT kcol, PLINT krow, 00674 PLFLT lastx, PLFLT lasty, PLINT startedge, PLINT **ipts, 00675 PLFLT *distance, PLINT *lastindex, 00676 void ( *pltr )( PLFLT, PLFLT, PLFLT *, PLFLT *, PLPointer ), 00677 PLPointer pltr_data ) 00678 { 00679 PLFLT f[4]; 00680 PLFLT px[4], py[4], locx[4], locy[4]; 00681 PLINT iedge[4]; 00682 PLINT i, j, k, num, first, inext, kcolnext, krownext; 00683 00684 00685 ( *pltr )( kcol, krow + 1, &px[0], &py[0], pltr_data ); 00686 ( *pltr )( kcol, krow, &px[1], &py[1], pltr_data ); 00687 ( *pltr )( kcol + 1, krow, &px[2], &py[2], pltr_data ); 00688 ( *pltr )( kcol + 1, krow + 1, &px[3], &py[3], pltr_data ); 00689 00690 f[0] = f2eval( kcol, krow + 1, f2eval_data ) - flev; 00691 f[1] = f2eval( kcol, krow, f2eval_data ) - flev; 00692 f[2] = f2eval( kcol + 1, krow, f2eval_data ) - flev; 00693 f[3] = f2eval( kcol + 1, krow + 1, f2eval_data ) - flev; 00694 00695 for ( i = 0, j = 1; i < 4; i++, j = ( j + 1 ) % 4 ) 00696 { 00697 iedge[i] = ( f[i] * f[j] > 0.0 ) ? -1 : ( ( f[i] * f[j] < 0.0 ) ? 1 : 0 ); 00698 } 00699 00700 /* Mark this square as done */ 00701 ipts[kcol][krow] = 1; 00702 00703 /* Check if no contour has been crossed i.e. iedge[i] = -1 */ 00704 if ( ( iedge[0] == -1 ) && ( iedge[1] == -1 ) && ( iedge[2] == -1 ) 00705 && ( iedge[3] == -1 ) ) 00706 return; 00707 00708 /* Check if this is a completely flat square - in which case 00709 * ignore it */ 00710 if ( ( f[0] == 0.0 ) && ( f[1] == 0.0 ) && ( f[2] == 0.0 ) && 00711 ( f[3] == 0.0 ) ) 00712 return; 00713 00714 /* Calculate intersection points */ 00715 num = 0; 00716 if ( startedge < 0 ) 00717 { 00718 first = 1; 00719 } 00720 else 00721 { 00722 locx[num] = lastx; 00723 locy[num] = lasty; 00724 num++; 00725 first = 0; 00726 } 00727 for ( k = 0, i = ( startedge < 0 ? 0 : startedge ); k < 4; k++, i = ( i + 1 ) % 4 ) 00728 { 00729 if ( i == startedge ) 00730 continue; 00731 00732 /* If the contour is an edge check it hasn't already been done */ 00733 if ( f[i] == 0.0 && f[( i + 1 ) % 4] == 0.0 ) 00734 { 00735 kcolnext = kcol; 00736 krownext = krow; 00737 if ( i == 0 ) 00738 kcolnext--; 00739 if ( i == 1 ) 00740 krownext--; 00741 if ( i == 2 ) 00742 kcolnext++; 00743 if ( i == 3 ) 00744 krownext++; 00745 if ( ( kcolnext < kx ) || ( kcolnext >= lx ) || 00746 ( krownext < ky ) || ( krownext >= ly ) || 00747 ( ipts[kcolnext][krownext] == 1 ) ) 00748 continue; 00749 } 00750 if ( ( iedge[i] == 1 ) || ( f[i] == 0.0 ) ) 00751 { 00752 j = ( i + 1 ) % 4; 00753 if ( f[i] != 0.0 ) 00754 { 00755 locx[num] = ( px[i] * fabs( f[j] ) + px[j] * fabs( f[i] ) ) / fabs( f[j] - f[i] ); 00756 locy[num] = ( py[i] * fabs( f[j] ) + py[j] * fabs( f[i] ) ) / fabs( f[j] - f[i] ); 00757 } 00758 else 00759 { 00760 locx[num] = px[i]; 00761 locy[num] = py[i]; 00762 } 00763 /* If this is the start of the contour then move to the point */ 00764 if ( first == 1 ) 00765 { 00766 cont_mv_store( locx[num], locy[num] ); 00767 first = 0; 00768 *distance = 0; 00769 *lastindex = 0; 00770 } 00771 else 00772 { 00773 /* Link to the next point on the contour */ 00774 if ( contlabel_active ) 00775 pl_drawcontlabel( locx[num], locy[num], flabel, distance, lastindex ); 00776 else 00777 cont_xy_store( locx[num], locy[num] ); 00778 /* Need to follow contour into next grid box */ 00779 /* Easy case where contour does not pass through corner */ 00780 if ( f[i] != 0.0 ) 00781 { 00782 kcolnext = kcol; 00783 krownext = krow; 00784 inext = ( i + 2 ) % 4; 00785 if ( i == 0 ) 00786 kcolnext--; 00787 if ( i == 1 ) 00788 krownext--; 00789 if ( i == 2 ) 00790 kcolnext++; 00791 if ( i == 3 ) 00792 krownext++; 00793 if ( ( kcolnext >= kx ) && ( kcolnext < lx ) && 00794 ( krownext >= ky ) && ( krownext < ly ) && 00795 ( ipts[kcolnext][krownext] == 0 ) ) 00796 { 00797 pldrawcn( f2eval, f2eval_data, 00798 nx, ny, kx, lx, ky, ly, flev, flabel, 00799 kcolnext, krownext, 00800 locx[num], locy[num], inext, ipts, 00801 distance, lastindex, 00802 pltr, pltr_data ); 00803 } 00804 } 00805 /* Hard case where contour passes through corner */ 00806 /* This is still not perfect - it may lose the contour 00807 * which won't upset the contour itself (we can find it 00808 * again later) but might upset the labelling */ 00809 else 00810 { 00811 kcolnext = kcol; 00812 krownext = krow; 00813 inext = ( i + 2 ) % 4; 00814 if ( i == 0 ) 00815 { 00816 kcolnext--; krownext++; 00817 } 00818 if ( i == 1 ) 00819 { 00820 krownext--; kcolnext--; 00821 } 00822 if ( i == 2 ) 00823 { 00824 kcolnext++; krownext--; 00825 } 00826 if ( i == 3 ) 00827 { 00828 krownext++; kcolnext++; 00829 } 00830 if ( ( kcolnext >= kx ) && ( kcolnext < lx ) && 00831 ( krownext >= ky ) && ( krownext < ly ) && 00832 ( ipts[kcolnext][krownext] == 0 ) ) 00833 { 00834 pldrawcn( f2eval, f2eval_data, 00835 nx, ny, kx, lx, ky, ly, flev, flabel, 00836 kcolnext, krownext, 00837 locx[num], locy[num], inext, ipts, 00838 distance, lastindex, 00839 pltr, pltr_data ); 00840 } 00841 } 00842 if ( first == 1 ) 00843 { 00844 /* Move back to first point */ 00845 cont_mv_store( locx[num], locy[num] ); 00846 first = 0; 00847 *distance = 0; 00848 *lastindex = 0; 00849 first = 0; 00850 } 00851 else 00852 { 00853 first = 1; 00854 } 00855 num++; 00856 } 00857 } 00858 } 00859 } 00860 00861 /*--------------------------------------------------------------------------*\ 00862 * pltr0() 00863 * 00864 * Identity transformation. 00865 \*--------------------------------------------------------------------------*/ 00866 00867 void 00868 pltr0( PLFLT x, PLFLT y, PLFLT *tx, PLFLT *ty, PLPointer pltr_data ) 00869 { 00870 *tx = x; 00871 *ty = y; 00872 } 00873 00874 /*--------------------------------------------------------------------------*\ 00875 * pltr1() 00876 * 00877 * Does linear interpolation from singly dimensioned coord arrays. 00878 * 00879 * Just abort for now if coordinates are out of bounds (don't think it's 00880 * possible, but if so we could use linear extrapolation). 00881 \*--------------------------------------------------------------------------*/ 00882 00883 void 00884 pltr1( PLFLT x, PLFLT y, PLFLT *tx, PLFLT *ty, PLPointer pltr_data ) 00885 { 00886 PLINT ul, ur, vl, vr; 00887 PLFLT du, dv; 00888 PLFLT xl, xr, yl, yr; 00889 00890 PLcGrid *grid = (PLcGrid *) pltr_data; 00891 PLFLT *xg = grid->xg; 00892 PLFLT *yg = grid->yg; 00893 PLINT nx = grid->nx; 00894 PLINT ny = grid->ny; 00895 00896 ul = (PLINT) x; 00897 ur = ul + 1; 00898 du = x - ul; 00899 00900 vl = (PLINT) y; 00901 vr = vl + 1; 00902 dv = y - vl; 00903 00904 if ( x < 0 || x > nx - 1 || y < 0 || y > ny - 1 ) 00905 { 00906 plexit( "pltr1: Invalid coordinates" ); 00907 } 00908 00909 /* Look up coordinates in row-dominant array. 00910 * Have to handle right boundary specially -- if at the edge, we'd better 00911 * not reference the out of bounds point. 00912 */ 00913 00914 xl = xg[ul]; 00915 yl = yg[vl]; 00916 00917 if ( ur == nx ) 00918 { 00919 *tx = xl; 00920 } 00921 else 00922 { 00923 xr = xg[ur]; 00924 *tx = xl * ( 1 - du ) + xr * du; 00925 } 00926 if ( vr == ny ) 00927 { 00928 *ty = yl; 00929 } 00930 else 00931 { 00932 yr = yg[vr]; 00933 *ty = yl * ( 1 - dv ) + yr * dv; 00934 } 00935 } 00936 00937 /*--------------------------------------------------------------------------*\ 00938 * pltr2() 00939 * 00940 * Does linear interpolation from doubly dimensioned coord arrays (column 00941 * dominant, as per normal C 2d arrays). 00942 * 00943 * This routine includes lots of checks for out of bounds. This would occur 00944 * occasionally due to some bugs in the contour plotter (now fixed). If an 00945 * out of bounds coordinate is obtained, the boundary value is provided 00946 * along with a warning. These checks should stay since no harm is done if 00947 * if everything works correctly. 00948 \*--------------------------------------------------------------------------*/ 00949 00950 void 00951 pltr2( PLFLT x, PLFLT y, PLFLT *tx, PLFLT *ty, PLPointer pltr_data ) 00952 { 00953 PLINT ul, ur, vl, vr; 00954 PLFLT du, dv; 00955 PLFLT xll, xlr, xrl, xrr; 00956 PLFLT yll, ylr, yrl, yrr; 00957 PLFLT xmin, xmax, ymin, ymax; 00958 00959 PLcGrid2 *grid = (PLcGrid2 *) pltr_data; 00960 PLFLT **xg = grid->xg; 00961 PLFLT **yg = grid->yg; 00962 PLINT nx = grid->nx; 00963 PLINT ny = grid->ny; 00964 00965 ul = (PLINT) x; 00966 ur = ul + 1; 00967 du = x - ul; 00968 00969 vl = (PLINT) y; 00970 vr = vl + 1; 00971 dv = y - vl; 00972 00973 xmin = 0; 00974 xmax = nx - 1; 00975 ymin = 0; 00976 ymax = ny - 1; 00977 00978 if ( x < xmin || x > xmax || y < ymin || y > ymax ) 00979 { 00980 plwarn( "pltr2: Invalid coordinates" ); 00981 if ( x < xmin ) 00982 { 00983 if ( y < ymin ) 00984 { 00985 *tx = xg[0][0]; 00986 *ty = yg[0][0]; 00987 } 00988 else if ( y > ymax ) 00989 { 00990 *tx = xg[0][ny - 1]; 00991 *ty = yg[0][ny - 1]; 00992 } 00993 else 00994 { 00995 xll = xg[0][vl]; 00996 yll = yg[0][vl]; 00997 xlr = xg[0][vr]; 00998 ylr = yg[0][vr]; 00999 01000 *tx = xll * ( 1 - dv ) + xlr * ( dv ); 01001 *ty = yll * ( 1 - dv ) + ylr * ( dv ); 01002 } 01003 } 01004 else if ( x > xmax ) 01005 { 01006 if ( y < ymin ) 01007 { 01008 *tx = xg[nx - 1][0]; 01009 *ty = yg[nx - 1][0]; 01010 } 01011 else if ( y > ymax ) 01012 { 01013 *tx = xg[nx - 1][ny - 1]; 01014 *ty = yg[nx - 1][ny - 1]; 01015 } 01016 else 01017 { 01018 xll = xg[nx - 1][vl]; 01019 yll = yg[nx - 1][vl]; 01020 xlr = xg[nx - 1][vr]; 01021 ylr = yg[nx - 1][vr]; 01022 01023 *tx = xll * ( 1 - dv ) + xlr * ( dv ); 01024 *ty = yll * ( 1 - dv ) + ylr * ( dv ); 01025 } 01026 } 01027 else 01028 { 01029 if ( y < ymin ) 01030 { 01031 xll = xg[ul][0]; 01032 xrl = xg[ur][0]; 01033 yll = yg[ul][0]; 01034 yrl = yg[ur][0]; 01035 01036 *tx = xll * ( 1 - du ) + xrl * ( du ); 01037 *ty = yll * ( 1 - du ) + yrl * ( du ); 01038 } 01039 else if ( y > ymax ) 01040 { 01041 xlr = xg[ul][ny - 1]; 01042 xrr = xg[ur][ny - 1]; 01043 ylr = yg[ul][ny - 1]; 01044 yrr = yg[ur][ny - 1]; 01045 01046 *tx = xlr * ( 1 - du ) + xrr * ( du ); 01047 *ty = ylr * ( 1 - du ) + yrr * ( du ); 01048 } 01049 } 01050 } 01051 01052 /* Normal case. 01053 * Look up coordinates in row-dominant array. 01054 * Have to handle right boundary specially -- if at the edge, we'd 01055 * better not reference the out of bounds point. 01056 */ 01057 01058 else 01059 { 01060 xll = xg[ul][vl]; 01061 yll = yg[ul][vl]; 01062 01063 /* ur is out of bounds */ 01064 01065 if ( ur == nx && vr < ny ) 01066 { 01067 xlr = xg[ul][vr]; 01068 ylr = yg[ul][vr]; 01069 01070 *tx = xll * ( 1 - dv ) + xlr * ( dv ); 01071 *ty = yll * ( 1 - dv ) + ylr * ( dv ); 01072 } 01073 01074 /* vr is out of bounds */ 01075 01076 else if ( ur < nx && vr == ny ) 01077 { 01078 xrl = xg[ur][vl]; 01079 yrl = yg[ur][vl]; 01080 01081 *tx = xll * ( 1 - du ) + xrl * ( du ); 01082 *ty = yll * ( 1 - du ) + yrl * ( du ); 01083 } 01084 01085 /* both ur and vr are out of bounds */ 01086 01087 else if ( ur == nx && vr == ny ) 01088 { 01089 *tx = xll; 01090 *ty = yll; 01091 } 01092 01093 /* everything in bounds */ 01094 01095 else 01096 { 01097 xrl = xg[ur][vl]; 01098 xlr = xg[ul][vr]; 01099 xrr = xg[ur][vr]; 01100 01101 yrl = yg[ur][vl]; 01102 ylr = yg[ul][vr]; 01103 yrr = yg[ur][vr]; 01104 01105 *tx = xll * ( 1 - du ) * ( 1 - dv ) + xlr * ( 1 - du ) * ( dv ) + 01106 xrl * ( du ) * ( 1 - dv ) + xrr * ( du ) * ( dv ); 01107 01108 *ty = yll * ( 1 - du ) * ( 1 - dv ) + ylr * ( 1 - du ) * ( dv ) + 01109 yrl * ( du ) * ( 1 - dv ) + yrr * ( du ) * ( dv ); 01110 } 01111 } 01112 } 01113 01114 /*--------------------------------------------------------------------------*\ 01115 * pltr2p() 01116 * 01117 * Just like pltr2() but uses pointer arithmetic to get coordinates from 2d 01118 * grid tables. This form of grid tables is compatible with those from 01119 * PLplot 4.0. The grid data must be pointed to by a PLcGrid structure. 01120 \*--------------------------------------------------------------------------*/ 01121 01122 void 01123 pltr2p( PLFLT x, PLFLT y, PLFLT *tx, PLFLT *ty, PLPointer pltr_data ) 01124 { 01125 PLINT ul, ur, vl, vr; 01126 PLFLT du, dv; 01127 PLFLT xll, xlr, xrl, xrr; 01128 PLFLT yll, ylr, yrl, yrr; 01129 PLFLT xmin, xmax, ymin, ymax; 01130 01131 PLcGrid *grid = (PLcGrid *) pltr_data; 01132 PLFLT *xg = grid->xg; 01133 PLFLT *yg = grid->yg; 01134 PLINT nx = grid->nx; 01135 PLINT ny = grid->ny; 01136 01137 ul = (PLINT) x; 01138 ur = ul + 1; 01139 du = x - ul; 01140 01141 vl = (PLINT) y; 01142 vr = vl + 1; 01143 dv = y - vl; 01144 01145 xmin = 0; 01146 xmax = nx - 1; 01147 ymin = 0; 01148 ymax = ny - 1; 01149 01150 if ( x < xmin || x > xmax || y < ymin || y > ymax ) 01151 { 01152 plwarn( "pltr2p: Invalid coordinates" ); 01153 if ( x < xmin ) 01154 { 01155 if ( y < ymin ) 01156 { 01157 *tx = *xg; 01158 *ty = *yg; 01159 } 01160 else if ( y > ymax ) 01161 { 01162 *tx = *( xg + ( ny - 1 ) ); 01163 *ty = *( yg + ( ny - 1 ) ); 01164 } 01165 else 01166 { 01167 ul = 0; 01168 xll = *( xg + ul * ny + vl ); 01169 yll = *( yg + ul * ny + vl ); 01170 xlr = *( xg + ul * ny + vr ); 01171 ylr = *( yg + ul * ny + vr ); 01172 01173 *tx = xll * ( 1 - dv ) + xlr * ( dv ); 01174 *ty = yll * ( 1 - dv ) + ylr * ( dv ); 01175 } 01176 } 01177 else if ( x > xmax ) 01178 { 01179 if ( y < ymin ) 01180 { 01181 *tx = *( xg + ( ny - 1 ) * nx ); 01182 *ty = *( yg + ( ny - 1 ) * nx ); 01183 } 01184 else if ( y > ymax ) 01185 { 01186 *tx = *( xg + ( ny - 1 ) + ( nx - 1 ) * ny ); 01187 *ty = *( yg + ( ny - 1 ) + ( nx - 1 ) * ny ); 01188 } 01189 else 01190 { 01191 ul = nx - 1; 01192 xll = *( xg + ul * ny + vl ); 01193 yll = *( yg + ul * ny + vl ); 01194 xlr = *( xg + ul * ny + vr ); 01195 ylr = *( yg + ul * ny + vr ); 01196 01197 *tx = xll * ( 1 - dv ) + xlr * ( dv ); 01198 *ty = yll * ( 1 - dv ) + ylr * ( dv ); 01199 } 01200 } 01201 else 01202 { 01203 if ( y < ymin ) 01204 { 01205 vl = 0; 01206 xll = *( xg + ul * ny + vl ); 01207 xrl = *( xg + ur * ny + vl ); 01208 yll = *( yg + ul * ny + vl ); 01209 yrl = *( yg + ur * ny + vl ); 01210 01211 *tx = xll * ( 1 - du ) + xrl * ( du ); 01212 *ty = yll * ( 1 - du ) + yrl * ( du ); 01213 } 01214 else if ( y > ymax ) 01215 { 01216 vr = ny - 1; 01217 xlr = *( xg + ul * ny + vr ); 01218 xrr = *( xg + ur * ny + vr ); 01219 ylr = *( yg + ul * ny + vr ); 01220 yrr = *( yg + ur * ny + vr ); 01221 01222 *tx = xlr * ( 1 - du ) + xrr * ( du ); 01223 *ty = ylr * ( 1 - du ) + yrr * ( du ); 01224 } 01225 } 01226 } 01227 01228 /* Normal case. 01229 * Look up coordinates in row-dominant array. 01230 * Have to handle right boundary specially -- if at the edge, we'd better 01231 * not reference the out of bounds point. 01232 */ 01233 01234 else 01235 { 01236 xll = *( xg + ul * ny + vl ); 01237 yll = *( yg + ul * ny + vl ); 01238 01239 /* ur is out of bounds */ 01240 01241 if ( ur == nx && vr < ny ) 01242 { 01243 xlr = *( xg + ul * ny + vr ); 01244 ylr = *( yg + ul * ny + vr ); 01245 01246 *tx = xll * ( 1 - dv ) + xlr * ( dv ); 01247 *ty = yll * ( 1 - dv ) + ylr * ( dv ); 01248 } 01249 01250 /* vr is out of bounds */ 01251 01252 else if ( ur < nx && vr == ny ) 01253 { 01254 xrl = *( xg + ur * ny + vl ); 01255 yrl = *( yg + ur * ny + vl ); 01256 01257 *tx = xll * ( 1 - du ) + xrl * ( du ); 01258 *ty = yll * ( 1 - du ) + yrl * ( du ); 01259 } 01260 01261 /* both ur and vr are out of bounds */ 01262 01263 else if ( ur == nx && vr == ny ) 01264 { 01265 *tx = xll; 01266 *ty = yll; 01267 } 01268 01269 /* everything in bounds */ 01270 01271 else 01272 { 01273 xrl = *( xg + ur * ny + vl ); 01274 xlr = *( xg + ul * ny + vr ); 01275 xrr = *( xg + ur * ny + vr ); 01276 01277 yrl = *( yg + ur * ny + vl ); 01278 ylr = *( yg + ul * ny + vr ); 01279 yrr = *( yg + ur * ny + vr ); 01280 01281 *tx = xll * ( 1 - du ) * ( 1 - dv ) + xlr * ( 1 - du ) * ( dv ) + 01282 xrl * ( du ) * ( 1 - dv ) + xrr * ( du ) * ( dv ); 01283 01284 *ty = yll * ( 1 - du ) * ( 1 - dv ) + ylr * ( 1 - du ) * ( dv ) + 01285 yrl * ( du ) * ( 1 - dv ) + yrr * ( du ) * ( dv ); 01286 } 01287 } 01288 } 01289 01290 /*----------------------------------------------------------------------*\ 01291 * pltr2f() 01292 * 01293 * Does linear interpolation from doubly dimensioned coord arrays 01294 * (row dominant, i.e. Fortran ordering). 01295 * 01296 * This routine includes lots of checks for out of bounds. This would 01297 * occur occasionally due to a bug in the contour plotter that is now fixed. 01298 * If an out of bounds coordinate is obtained, the boundary value is provided 01299 * along with a warning. These checks should stay since no harm is done if 01300 * if everything works correctly. 01301 \*----------------------------------------------------------------------*/ 01302 01303 void 01304 pltr2f( PLFLT x, PLFLT y, PLFLT *tx, PLFLT *ty, void *pltr_data ) 01305 { 01306 PLINT ul, ur, vl, vr; 01307 PLFLT du, dv; 01308 PLFLT xll, xlr, xrl, xrr; 01309 PLFLT yll, ylr, yrl, yrr; 01310 PLFLT xmin, xmax, ymin, ymax; 01311 01312 PLcGrid *cgrid = (PLcGrid *) pltr_data; 01313 PLFLT *xg = cgrid->xg; 01314 PLFLT *yg = cgrid->yg; 01315 PLINT nx = cgrid->nx; 01316 PLINT ny = cgrid->ny; 01317 01318 ul = (PLINT) x; 01319 ur = ul + 1; 01320 du = x - ul; 01321 01322 vl = (PLINT) y; 01323 vr = vl + 1; 01324 dv = y - vl; 01325 01326 xmin = 0; 01327 xmax = nx - 1; 01328 ymin = 0; 01329 ymax = ny - 1; 01330 01331 if ( x < xmin || x > xmax || y < ymin || y > ymax ) 01332 { 01333 plwarn( "pltr2f: Invalid coordinates" ); 01334 01335 if ( x < xmin ) 01336 { 01337 if ( y < ymin ) 01338 { 01339 *tx = *xg; 01340 *ty = *yg; 01341 } 01342 else if ( y > ymax ) 01343 { 01344 *tx = *( xg + ( ny - 1 ) * nx ); 01345 *ty = *( yg + ( ny - 1 ) * nx ); 01346 } 01347 else 01348 { 01349 ul = 0; 01350 xll = *( xg + ul + vl * nx ); 01351 yll = *( yg + ul + vl * nx ); 01352 xlr = *( xg + ul + vr * nx ); 01353 ylr = *( yg + ul + vr * nx ); 01354 01355 *tx = xll * ( 1 - dv ) + xlr * ( dv ); 01356 *ty = yll * ( 1 - dv ) + ylr * ( dv ); 01357 } 01358 } 01359 else if ( x > xmax ) 01360 { 01361 if ( y < ymin ) 01362 { 01363 *tx = *( xg + ( nx - 1 ) ); 01364 *ty = *( yg + ( nx - 1 ) ); 01365 } 01366 else if ( y > ymax ) 01367 { 01368 *tx = *( xg + ( nx - 1 ) + ( ny - 1 ) * nx ); 01369 *ty = *( yg + ( nx - 1 ) + ( ny - 1 ) * nx ); 01370 } 01371 else 01372 { 01373 ul = nx - 1; 01374 xll = *( xg + ul + vl * nx ); 01375 yll = *( yg + ul + vl * nx ); 01376 xlr = *( xg + ul + vr * nx ); 01377 ylr = *( yg + ul + vr * nx ); 01378 01379 *tx = xll * ( 1 - dv ) + xlr * ( dv ); 01380 *ty = yll * ( 1 - dv ) + ylr * ( dv ); 01381 } 01382 } 01383 else 01384 { 01385 if ( y < ymin ) 01386 { 01387 vl = 0; 01388 xll = *( xg + ul + vl * nx ); 01389 xrl = *( xg + ur + vl * nx ); 01390 yll = *( yg + ul + vl * nx ); 01391 yrl = *( yg + ur + vl * nx ); 01392 01393 *tx = xll * ( 1 - du ) + xrl * ( du ); 01394 *ty = yll * ( 1 - du ) + yrl * ( du ); 01395 } 01396 else if ( y > ymax ) 01397 { 01398 vr = ny - 1; 01399 xlr = *( xg + ul + vr * nx ); 01400 xrr = *( xg + ur + vr * nx ); 01401 ylr = *( yg + ul + vr * nx ); 01402 yrr = *( yg + ur + vr * nx ); 01403 01404 *tx = xlr * ( 1 - du ) + xrr * ( du ); 01405 *ty = ylr * ( 1 - du ) + yrr * ( du ); 01406 } 01407 } 01408 } 01409 01410 /* Normal case. 01411 * Look up coordinates in row-dominant array. 01412 * Have to handle right boundary specially -- if at the edge, we'd 01413 * better not reference the out of bounds point. */ 01414 01415 else 01416 { 01417 xll = *( xg + ul + vl * nx ); 01418 yll = *( yg + ul + vl * nx ); 01419 01420 /* ur is out of bounds */ 01421 01422 if ( ur == nx && vr < ny ) 01423 { 01424 xlr = *( xg + ul + vr * nx ); 01425 ylr = *( yg + ul + vr * nx ); 01426 01427 *tx = xll * ( 1 - dv ) + xlr * ( dv ); 01428 *ty = yll * ( 1 - dv ) + ylr * ( dv ); 01429 } 01430 01431 /* vr is out of bounds */ 01432 01433 else if ( ur < nx && vr == ny ) 01434 { 01435 xrl = *( xg + ur + vl * nx ); 01436 yrl = *( yg + ur + vl * nx ); 01437 01438 *tx = xll * ( 1 - du ) + xrl * ( du ); 01439 *ty = yll * ( 1 - du ) + yrl * ( du ); 01440 } 01441 01442 /* both ur and vr are out of bounds */ 01443 01444 else if ( ur == nx && vr == ny ) 01445 { 01446 *tx = xll; 01447 *ty = yll; 01448 } 01449 01450 /* everything in bounds */ 01451 01452 else 01453 { 01454 xrl = *( xg + ur + vl * nx ); 01455 xlr = *( xg + ul + vr * nx ); 01456 xrr = *( xg + ur + vr * nx ); 01457 01458 yrl = *( yg + ur + vl * nx ); 01459 ylr = *( yg + ul + vr * nx ); 01460 yrr = *( yg + ur + vr * nx ); 01461 /* INDENT OFF */ 01462 *tx = xll * ( 1 - du ) * ( 1 - dv ) + xlr * ( 1 - du ) * ( dv ) + 01463 xrl * ( du ) * ( 1 - dv ) + xrr * ( du ) * ( dv ); 01464 01465 *ty = yll * ( 1 - du ) * ( 1 - dv ) + ylr * ( 1 - du ) * ( dv ) + 01466 yrl * ( du ) * ( 1 - dv ) + yrr * ( du ) * ( dv ); 01467 /* INDENT ON */ 01468 } 01469 } 01470 }