// Contour plotter. // // Copyright (C) 1995, 2000, 2001 Maurice LeBrun // Copyright (C) 2000, 2002 Joao Cardoso // Copyright (C) 2000-2014 Alan W. Irwin // Copyright (C) 2004 Andrew Ross // // This file is part of PLplot. // // PLplot is free software; you can redistribute it and/or modify // it under the terms of the GNU Library General Public License as published // by the Free Software Foundation; either version 2 of the License, or // (at your option) any later version. // // PLplot is distributed in the hope that it will be useful, // but WITHOUT ANY WARRANTY; without even the implied warranty of // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the // GNU Library General Public License for more details. // // You should have received a copy of the GNU Library General Public License // along with PLplot; if not, write to the Free Software // Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA // #include "plplotP.h" #ifdef MSDOS #pragma optimize("",off) #endif // Static function prototypes. static void plcntr( PLF2EVAL_callback plf2eval, PLPointer plf2eval_data, PLINT nx, PLINT ny, PLINT kx, PLINT lx, PLINT ky, PLINT ly, PLFLT flev, PLINT **ipts, PLTRANSFORM_callback pltr, PLPointer pltr_data ); static void pldrawcn( PLF2EVAL_callback plf2eval, PLPointer plf2eval_data, PLINT nx, PLINT ny, PLINT kx, PLINT lx, PLINT ky, PLINT ly, PLFLT flev, char *flabel, PLINT kcol, PLINT krow, PLFLT lastx, PLFLT lasty, PLINT startedge, PLINT **ipts, PLFLT *distance, PLINT *lastindex, PLTRANSFORM_callback pltr, PLPointer pltr_data ); static void plfloatlabel( PLFLT value, char *string, PLINT len ); static PLFLT plP_pcwcx( PLINT x ); static PLFLT plP_pcwcy( PLINT y ); static void pl_drawcontlabel( PLFLT tpx, PLFLT tpy, char *flabel, PLFLT *distance, PLINT *lastindex ); // Error flag for aborts static int error; //************************************** // // Defaults for contour label printing. // //************************************** // Font height for contour labels (normalized) static PLFLT contlabel_size = 0.3; // Offset of label from contour line (if set to 0.0, labels are printed on the lines). static PLFLT contlabel_offset = 0.006; // Spacing parameter for contour labels static PLFLT contlabel_space = 0.1; // Activate labels, default off static PLINT contlabel_active = 0; // If the contour label exceed 10^(limexp) or 10^(-limexp), the exponential format is used static PLINT limexp = 4; // Number of significant digits static PLINT sigprec = 2; //******* contour lines storage *************************** static CONT_LEVEL *startlev = NULL; static CONT_LEVEL *currlev; static CONT_LINE *currline; static int cont3d = 0; static CONT_LINE * alloc_line( void ) { CONT_LINE *line; if ( ( line = (CONT_LINE *) malloc( sizeof ( CONT_LINE ) ) ) == NULL ) { plexit( "alloc_line: Insufficient memory" ); } line->x = (PLFLT *) malloc( LINE_ITEMS * sizeof ( PLFLT ) ); line->y = (PLFLT *) malloc( LINE_ITEMS * sizeof ( PLFLT ) ); if ( ( line->x == NULL ) || ( line->y == NULL ) ) { plexit( "alloc_line: Insufficient memory" ); } line->npts = 0; line->next = NULL; return line; } static CONT_LEVEL * alloc_level( PLFLT level ) { CONT_LEVEL *node; if ( ( node = (CONT_LEVEL *) malloc( sizeof ( CONT_LEVEL ) ) ) == NULL ) { plexit( "alloc_level: Insufficient memory" ); } node->level = level; node->next = NULL; node->line = alloc_line( ); return node; } static void realloc_line( CONT_LINE *line ) { if ( ( ( line->x = (PLFLT *) realloc( line->x, (size_t) ( line->npts + LINE_ITEMS ) * sizeof ( PLFLT ) ) ) == NULL ) || ( ( line->y = (PLFLT *) realloc( line->y, (size_t) ( line->npts + LINE_ITEMS ) * sizeof ( PLFLT ) ) ) == NULL ) ) plexit( "realloc_line: Insufficient memory" ); } // new contour level static void cont_new_store( PLFLT level ) { if ( cont3d ) { if ( startlev == NULL ) { startlev = alloc_level( level ); currlev = startlev; } else { currlev->next = alloc_level( level ); currlev = currlev->next; } currline = currlev->line; } } void cont_clean_store( CONT_LEVEL *ct ) { CONT_LINE *tline, *cline; CONT_LEVEL *tlev, *clevel; if ( ct != NULL ) { clevel = ct; do { cline = clevel->line; do { #ifdef CONT_PLOT_DEBUG plP_movwor( cline->x[0], cline->y[0] ); for ( j = 1; j < cline->npts; j++ ) plP_drawor( cline->x[j], cline->y[j] ); #endif tline = cline->next; free( cline->x ); free( cline->y ); free( cline ); cline = tline; } while ( cline != NULL ); tlev = clevel->next; free( clevel ); clevel = tlev; } while ( clevel != NULL ); startlev = NULL; } } static void cont_xy_store( PLFLT xx, PLFLT yy ) { if ( cont3d ) { PLINT pts = currline->npts; if ( pts % LINE_ITEMS == 0 ) realloc_line( currline ); currline->x[pts] = xx; currline->y[pts] = yy; currline->npts++; } else plP_drawor( xx, yy ); } static void cont_mv_store( PLFLT xx, PLFLT yy ) { if ( cont3d ) { if ( currline->npts != 0 ) // not an empty list, allocate new { currline->next = alloc_line( ); currline = currline->next; } // and fill first element currline->x[0] = xx; currline->y[0] = yy; currline->npts = 1; } else plP_movwor( xx, yy ); } // small routine to set offset and spacing of contour labels, see desciption above void c_pl_setcontlabelparam( PLFLT offset, PLFLT size, PLFLT spacing, PLINT active ) { contlabel_offset = offset; contlabel_size = size; contlabel_space = spacing; contlabel_active = active; } // small routine to set the format of the contour labels, description of limexp and prec see above void c_pl_setcontlabelformat( PLINT lexp, PLINT sigdig ) { limexp = lexp; sigprec = sigdig; } static void pl_drawcontlabel( PLFLT tpx, PLFLT tpy, char *flabel, PLFLT *distance, PLINT *lastindex ) { PLFLT delta_x, delta_y; PLINT currx_old, curry_old; delta_x = plP_pcdcx( plsc->currx ) - plP_pcdcx( plP_wcpcx( tpx ) ); delta_y = plP_pcdcy( plsc->curry ) - plP_pcdcy( plP_wcpcy( tpy ) ); currx_old = plsc->currx; curry_old = plsc->curry; *distance += sqrt( delta_x * delta_x + delta_y * delta_y ); plP_drawor( tpx, tpy ); if ( (int) ( fabs( *distance / contlabel_space ) ) > *lastindex ) { PLFLT scale, vec_x, vec_y, mx, my, dev_x, dev_y, off_x, off_y; vec_x = tpx - plP_pcwcx( currx_old ); vec_y = tpy - plP_pcwcy( curry_old ); // Ensure labels appear the right way up if ( vec_x < 0 ) { vec_x = -vec_x; vec_y = -vec_y; } mx = (double) plsc->wpxscl / (double) plsc->phyxlen; my = (double) plsc->wpyscl / (double) plsc->phyylen; dev_x = -my * vec_y / mx; dev_y = mx * vec_x / my; scale = sqrt( ( mx * mx * dev_x * dev_x + my * my * dev_y * dev_y ) / ( contlabel_offset * contlabel_offset ) ); off_x = dev_x / scale; off_y = dev_y / scale; plptex( tpx + off_x, tpy + off_y, vec_x, vec_y, 0.5, flabel ); plP_movwor( tpx, tpy ); ( *lastindex )++; } else plP_movwor( tpx, tpy ); } // Format contour labels. Arguments: // value: floating point number to be formatted // string: the formatted label, plptex must be called with it to actually // print the label // static void plfloatlabel( PLFLT value, char *string, PLINT len ) { PLINT setpre, precis; // form[10] gives enough space for all non-malicious formats. // tmpstring[15] gives enough room for 3 digits in a negative exponent // or 4 digits in a positive exponent + null termination. That // should be enough for all non-malicious use. // Obviously there are security issues here that // should be addressed as well. // #define FORM_LEN 10 #define TMPSTRING_LEN 15 char form[FORM_LEN], tmpstring[TMPSTRING_LEN]; PLINT exponent = 0; PLFLT mant, tmp; PLINT prec = sigprec; plP_gprec( &setpre, &precis ); if ( setpre ) prec = precis; if ( value > 0.0 ) tmp = log10( value ); else if ( value < 0.0 ) tmp = log10( -value ); else tmp = 0; if ( tmp >= 0.0 ) exponent = (int) tmp; else if ( tmp < 0.0 ) { tmp = -tmp; if ( floor( tmp ) < tmp ) exponent = -(int) ( floor( tmp ) + 1.0 ); else exponent = -(int) ( floor( tmp ) ); } mant = value / pow( 10.0, exponent ); if ( mant != 0.0 ) mant = (int) ( mant * pow( 10.0, prec - 1 ) + 0.5 * mant / fabs( mant ) ) / pow( 10.0, prec - 1 ); snprintf( form, FORM_LEN, "%%.%df", prec - 1 ); snprintf( string, (size_t) len, form, mant ); snprintf( tmpstring, TMPSTRING_LEN, "#(229)10#u%d", exponent ); strncat( string, tmpstring, (size_t) len - strlen( string ) - 1 ); if ( abs( exponent ) < limexp || value == 0.0 ) { value = pow( 10.0, exponent ) * mant; if ( exponent >= 0 ) prec = prec - 1 - exponent; else prec = prec - 1 + abs( exponent ); if ( prec < 0 ) prec = 0; snprintf( form, FORM_LEN, "%%.%df", (int) prec ); snprintf( string, (size_t) len, form, value ); } } // physical coords (x) to world coords static PLFLT plP_pcwcx( PLINT x ) { return ( ( x - plsc->wpxoff ) / plsc->wpxscl ); } // physical coords (y) to world coords static PLFLT plP_pcwcy( PLINT y ) { return ( ( y - plsc->wpyoff ) / plsc->wpyscl ); } //-------------------------------------------------------------------------- // plf2eval1() // // Does a lookup from a 2d function array. Array is of type (PLFLT **), // and is column dominant (normal C ordering). //-------------------------------------------------------------------------- PLFLT plf2eval1( PLINT ix, PLINT iy, PLPointer plf2eval_data ) { PLFLT value; PLFLT **z = (PLFLT **) plf2eval_data; value = z[ix][iy]; return value; } //-------------------------------------------------------------------------- // plf2eval2() // // Does a lookup from a 2d function array. plf2eval_data is treated as type // (PLfGrid2 *). //-------------------------------------------------------------------------- PLFLT plf2eval2( PLINT ix, PLINT iy, PLPointer plf2eval_data ) { PLFLT value; PLfGrid2 *grid = (PLfGrid2 *) plf2eval_data; value = grid->f[ix][iy]; return value; } //-------------------------------------------------------------------------- // plf2eval() // // Does a lookup from a 2d function array. Array is of type (PLFLT *), and // is column dominant (normal C ordering). You MUST fill the ny maximum // array index entry in the PLfGrid struct. //-------------------------------------------------------------------------- PLFLT plf2eval( PLINT ix, PLINT iy, PLPointer plf2eval_data ) { PLFLT value; PLfGrid *grid = (PLfGrid *) plf2eval_data; value = grid->f[ix * grid->ny + iy]; return value; } //-------------------------------------------------------------------------- // plf2evalr() // // Does a lookup from a 2d function array. Array is of type (PLFLT *), and // is row dominant (Fortran ordering). You MUST fill the nx maximum array // index entry in the PLfGrid struct. //-------------------------------------------------------------------------- PLFLT plf2evalr( PLINT ix, PLINT iy, PLPointer plf2eval_data ) { PLFLT value; PLfGrid *grid = (PLfGrid *) plf2eval_data; value = grid->f[ix + iy * grid->nx]; return value; } //-------------------------------------------------------------------------- // // cont_store: // // Draw contour lines in memory. // cont_clean_store() must be called after use to release allocated memory. // //-------------------------------------------------------------------------- void cont_store( PLFLT_MATRIX f, PLINT nx, PLINT ny, PLINT kx, PLINT lx, PLINT ky, PLINT ly, PLFLT_VECTOR clevel, PLINT nlevel, PLTRANSFORM_callback pltr, PLPointer pltr_data, CONT_LEVEL **contour ) { cont3d = 1; plcont( f, nx, ny, kx, lx, ky, ly, clevel, nlevel, pltr, pltr_data ); *contour = startlev; cont3d = 0; } //-------------------------------------------------------------------------- // void plcont() // // Draws a contour plot from data in f(nx,ny). Is just a front-end to // plfcont, with a particular choice for f2eval and f2eval_data. //-------------------------------------------------------------------------- void c_plcont( PLFLT_MATRIX f, PLINT nx, PLINT ny, PLINT kx, PLINT lx, PLINT ky, PLINT ly, PLFLT_VECTOR clevel, PLINT nlevel, PLTRANSFORM_callback pltr, PLPointer pltr_data ) { plfcont( plf2eval1, (PLPointer) f, nx, ny, kx, lx, ky, ly, clevel, nlevel, pltr, pltr_data ); } //-------------------------------------------------------------------------- // void plfcont() // // Draws a contour plot using the function evaluator f2eval and data stored // by way of the f2eval_data pointer. This allows arbitrary organizations // of 2d array data to be used. // // The subrange of indices used for contouring is kx to lx in the x // direction and from ky to ly in the y direction. The array of contour // levels is clevel(nlevel), and "pltr" is the name of a function which // transforms array indicies into world coordinates. // // Note that the fortran-like minimum and maximum indices (kx, lx, ky, ly) // are translated into more C-like ones. I've only kept them as they are // for the plfcont() argument list because of backward compatibility. //-------------------------------------------------------------------------- void plfcont( PLF2EVAL_callback f2eval, PLPointer f2eval_data, PLINT nx, PLINT ny, PLINT kx, PLINT lx, PLINT ky, PLINT ly, PLFLT_VECTOR clevel, PLINT nlevel, PLTRANSFORM_callback pltr, PLPointer pltr_data ) { PLINT i, **ipts; if ( pltr == NULL ) { // If pltr is undefined, abort with an error. plabort( "plfcont: The pltr callback must be defined" ); return; } if ( kx < 1 || kx >= lx ) { plabort( "plfcont: indices must satisfy 1 <= kx <= lx <= nx" ); return; } if ( ky < 1 || ky >= ly ) { plabort( "plfcont: indices must satisfy 1 <= ky <= ly <= ny" ); return; } if ( ( ipts = (PLINT **) malloc( (size_t) nx * sizeof ( PLINT * ) ) ) == NULL ) { plexit( "plfcont: Insufficient memory" ); } for ( i = 0; i < nx; i++ ) { if ( ( ipts[i] = (PLINT *) malloc( (size_t) ny * sizeof ( PLINT * ) ) ) == NULL ) { plexit( "plfcont: Insufficient memory" ); } } for ( i = 0; i < nlevel; i++ ) { plcntr( f2eval, f2eval_data, nx, ny, kx - 1, lx - 1, ky - 1, ly - 1, clevel[i], ipts, pltr, pltr_data ); if ( error ) { error = 0; goto done; } } done: for ( i = 0; i < nx; i++ ) { free( (void *) ipts[i] ); } free( (void *) ipts ); } //-------------------------------------------------------------------------- // void plcntr() // // The contour for a given level is drawn here. Note iscan has nx // elements. ixstor and iystor each have nstor elements. //-------------------------------------------------------------------------- static void plcntr( PLF2EVAL_callback f2eval, PLPointer f2eval_data, PLINT nx, PLINT ny, PLINT kx, PLINT lx, PLINT ky, PLINT ly, PLFLT flev, PLINT **ipts, PLTRANSFORM_callback pltr, PLPointer pltr_data ) { PLINT kcol, krow, lastindex; PLFLT distance; PLFLT save_def, save_scale; char flabel[30]; plgchr( &save_def, &save_scale ); save_scale = save_scale / save_def; cont_new_store( flev ); // format contour label for plptex and define the font height of the labels plfloatlabel( flev, flabel, 30 ); plschr( 0.0, contlabel_size ); // Clear array for traversed squares for ( kcol = kx; kcol < lx; kcol++ ) { for ( krow = ky; krow < ly; krow++ ) { ipts[kcol][krow] = 0; } } for ( krow = ky; krow < ly; krow++ ) { for ( kcol = kx; kcol < lx; kcol++ ) { if ( ipts[kcol][krow] == 0 ) { // Follow and draw a contour pldrawcn( f2eval, f2eval_data, nx, ny, kx, lx, ky, ly, flev, flabel, kcol, krow, 0.0, 0.0, -2, ipts, &distance, &lastindex, pltr, pltr_data ); if ( error ) return; } } } plschr( save_def, save_scale ); } //-------------------------------------------------------------------------- // void pldrawcn() // // Follow and draw a contour. //-------------------------------------------------------------------------- static void pldrawcn( PLF2EVAL_callback f2eval, PLPointer f2eval_data, PLINT nx, PLINT ny, PLINT kx, PLINT lx, PLINT ky, PLINT ly, PLFLT flev, char *flabel, PLINT kcol, PLINT krow, PLFLT lastx, PLFLT lasty, PLINT startedge, PLINT **ipts, PLFLT *distance, PLINT *lastindex, PLTRANSFORM_callback pltr, PLPointer pltr_data ) { PLFLT f[4]; PLFLT px[4], py[4], locx[4], locy[4]; PLINT iedge[4]; PLINT i, j, k, num, first, inext, kcolnext, krownext, sfi, sfj; ( *pltr )( kcol, krow + 1, &px[0], &py[0], pltr_data ); ( *pltr )( kcol, krow, &px[1], &py[1], pltr_data ); ( *pltr )( kcol + 1, krow, &px[2], &py[2], pltr_data ); ( *pltr )( kcol + 1, krow + 1, &px[3], &py[3], pltr_data ); f[0] = f2eval( kcol, krow + 1, f2eval_data ) - flev; f[1] = f2eval( kcol, krow, f2eval_data ) - flev; f[2] = f2eval( kcol + 1, krow, f2eval_data ) - flev; f[3] = f2eval( kcol + 1, krow + 1, f2eval_data ) - flev; for ( i = 0, j = 1; i < 4; i++, j = ( j + 1 ) % 4 ) { // Use intermediates to avoid possible floating point // under / over flow during multiplication. sfi = ( f[i] > 0.0 ) ? 1 : ( ( f[i] < 0.0 ) ? -1 : 0 ); sfj = ( f[j] > 0.0 ) ? 1 : ( ( f[j] < 0.0 ) ? -1 : 0 ); iedge[i] = ( sfi * sfj > 0 ) ? -1 : ( ( sfi * sfj < 0 ) ? 1 : 0 ); } // Mark this square as done ipts[kcol][krow] = 1; // Check if no contour has been crossed i.e. iedge[i] = -1 if ( ( iedge[0] == -1 ) && ( iedge[1] == -1 ) && ( iedge[2] == -1 ) && ( iedge[3] == -1 ) ) return; // Check if this is a completely flat square - in which case // ignore it if ( ( f[0] == 0.0 ) && ( f[1] == 0.0 ) && ( f[2] == 0.0 ) && ( f[3] == 0.0 ) ) return; // Calculate intersection points num = 0; if ( startedge < 0 ) { first = 1; } else { locx[num] = lastx; locy[num] = lasty; num++; first = 0; } for ( k = 0, i = ( startedge < 0 ? 0 : startedge ); k < 4; k++, i = ( i + 1 ) % 4 ) { if ( i == startedge ) continue; // If the contour is an edge check it hasn't already been done if ( f[i] == 0.0 && f[( i + 1 ) % 4] == 0.0 ) { kcolnext = kcol; krownext = krow; if ( i == 0 ) kcolnext--; if ( i == 1 ) krownext--; if ( i == 2 ) kcolnext++; if ( i == 3 ) krownext++; if ( ( kcolnext < kx ) || ( kcolnext >= lx ) || ( krownext < ky ) || ( krownext >= ly ) || ( ipts[kcolnext][krownext] == 1 ) ) continue; } if ( ( iedge[i] == 1 ) || ( f[i] == 0.0 ) ) { j = ( i + 1 ) % 4; if ( f[i] != 0.0 ) { locx[num] = ( px[i] * fabs( f[j] ) + px[j] * fabs( f[i] ) ) / fabs( f[j] - f[i] ); locy[num] = ( py[i] * fabs( f[j] ) + py[j] * fabs( f[i] ) ) / fabs( f[j] - f[i] ); } else { locx[num] = px[i]; locy[num] = py[i]; } // If this is the start of the contour then move to the point if ( first == 1 ) { cont_mv_store( locx[num], locy[num] ); first = 0; *distance = 0; *lastindex = 0; } else { // Link to the next point on the contour if ( contlabel_active ) pl_drawcontlabel( locx[num], locy[num], flabel, distance, lastindex ); else cont_xy_store( locx[num], locy[num] ); // Need to follow contour into next grid box // Easy case where contour does not pass through corner if ( f[i] != 0.0 ) { kcolnext = kcol; krownext = krow; inext = ( i + 2 ) % 4; if ( i == 0 ) kcolnext--; if ( i == 1 ) krownext--; if ( i == 2 ) kcolnext++; if ( i == 3 ) krownext++; if ( ( kcolnext >= kx ) && ( kcolnext < lx ) && ( krownext >= ky ) && ( krownext < ly ) && ( ipts[kcolnext][krownext] == 0 ) ) { pldrawcn( f2eval, f2eval_data, nx, ny, kx, lx, ky, ly, flev, flabel, kcolnext, krownext, locx[num], locy[num], inext, ipts, distance, lastindex, pltr, pltr_data ); } } // Hard case where contour passes through corner // This is still not perfect - it may lose the contour // which won't upset the contour itself (we can find it // again later) but might upset the labelling else { kcolnext = kcol; krownext = krow; inext = ( i + 2 ) % 4; if ( i == 0 ) { kcolnext--; krownext++; } if ( i == 1 ) { krownext--; kcolnext--; } if ( i == 2 ) { kcolnext++; krownext--; } if ( i == 3 ) { krownext++; kcolnext++; } if ( ( kcolnext >= kx ) && ( kcolnext < lx ) && ( krownext >= ky ) && ( krownext < ly ) && ( ipts[kcolnext][krownext] == 0 ) ) { pldrawcn( f2eval, f2eval_data, nx, ny, kx, lx, ky, ly, flev, flabel, kcolnext, krownext, locx[num], locy[num], inext, ipts, distance, lastindex, pltr, pltr_data ); } } if ( first == 1 ) { // Move back to first point cont_mv_store( locx[num], locy[num] ); first = 0; *distance = 0; *lastindex = 0; first = 0; } else { first = 1; } num++; } } } } //-------------------------------------------------------------------------- // pltr0() // // Identity transformation. //-------------------------------------------------------------------------- void pltr0( PLFLT x, PLFLT y, PLFLT *tx, PLFLT *ty, PLPointer PL_UNUSED( pltr_data ) ) { *tx = x; *ty = y; } //-------------------------------------------------------------------------- // pltr1() // // Does linear interpolation from singly dimensioned coord arrays. // // Just abort for now if coordinates are out of bounds (don't think it's // possible, but if so we could use linear extrapolation). //-------------------------------------------------------------------------- void pltr1( PLFLT x, PLFLT y, PLFLT *tx, PLFLT *ty, PLPointer pltr_data ) { PLINT ul, ur, vl, vr; PLFLT du, dv; PLFLT xl, xr, yl, yr; PLcGrid *grid = (PLcGrid *) pltr_data; PLFLT *xg = grid->xg; PLFLT *yg = grid->yg; PLINT nx = grid->nx; PLINT ny = grid->ny; ul = (PLINT) x; ur = ul + 1; du = x - ul; vl = (PLINT) y; vr = vl + 1; dv = y - vl; if ( x < 0 || x > nx - 1 || y < 0 || y > ny - 1 ) { plexit( "pltr1: Invalid coordinates" ); } // Look up coordinates in row-dominant array. // Have to handle right boundary specially -- if at the edge, we'd better // not reference the out of bounds point. // xl = xg[ul]; yl = yg[vl]; if ( ur == nx ) { *tx = xl; } else { xr = xg[ur]; *tx = xl * ( 1 - du ) + xr * du; } if ( vr == ny ) { *ty = yl; } else { yr = yg[vr]; *ty = yl * ( 1 - dv ) + yr * dv; } } //-------------------------------------------------------------------------- // pltr2() // // Does linear interpolation from doubly dimensioned coord arrays (column // dominant, as per normal C 2d arrays). // // This routine includes lots of checks for out of bounds. This would occur // occasionally due to some bugs in the contour plotter (now fixed). If an // out of bounds coordinate is obtained, the boundary value is provided // along with a warning. These checks should stay since no harm is done if // if everything works correctly. //-------------------------------------------------------------------------- void pltr2( PLFLT x, PLFLT y, PLFLT *tx, PLFLT *ty, PLPointer pltr_data ) { PLINT ul, ur, vl, vr; PLFLT du, dv; PLFLT xll, xlr, xrl, xrr; PLFLT yll, ylr, yrl, yrr; PLFLT xmin, xmax, ymin, ymax; PLcGrid2 *grid = (PLcGrid2 *) pltr_data; PLFLT **xg = grid->xg; PLFLT **yg = grid->yg; PLINT nx = grid->nx; PLINT ny = grid->ny; ul = (PLINT) x; ur = ul + 1; du = x - ul; vl = (PLINT) y; vr = vl + 1; dv = y - vl; xmin = 0; xmax = nx - 1; ymin = 0; ymax = ny - 1; if ( x < xmin || x > xmax || y < ymin || y > ymax ) { plwarn( "pltr2: Invalid coordinates" ); if ( x < xmin ) { if ( y < ymin ) { *tx = xg[0][0]; *ty = yg[0][0]; } else if ( y > ymax ) { *tx = xg[0][ny - 1]; *ty = yg[0][ny - 1]; } else { xll = xg[0][vl]; yll = yg[0][vl]; xlr = xg[0][vr]; ylr = yg[0][vr]; *tx = xll * ( 1 - dv ) + xlr * ( dv ); *ty = yll * ( 1 - dv ) + ylr * ( dv ); } } else if ( x > xmax ) { if ( y < ymin ) { *tx = xg[nx - 1][0]; *ty = yg[nx - 1][0]; } else if ( y > ymax ) { *tx = xg[nx - 1][ny - 1]; *ty = yg[nx - 1][ny - 1]; } else { xll = xg[nx - 1][vl]; yll = yg[nx - 1][vl]; xlr = xg[nx - 1][vr]; ylr = yg[nx - 1][vr]; *tx = xll * ( 1 - dv ) + xlr * ( dv ); *ty = yll * ( 1 - dv ) + ylr * ( dv ); } } else { if ( y < ymin ) { xll = xg[ul][0]; xrl = xg[ur][0]; yll = yg[ul][0]; yrl = yg[ur][0]; *tx = xll * ( 1 - du ) + xrl * ( du ); *ty = yll * ( 1 - du ) + yrl * ( du ); } else if ( y > ymax ) { xlr = xg[ul][ny - 1]; xrr = xg[ur][ny - 1]; ylr = yg[ul][ny - 1]; yrr = yg[ur][ny - 1]; *tx = xlr * ( 1 - du ) + xrr * ( du ); *ty = ylr * ( 1 - du ) + yrr * ( du ); } } } // Normal case. // Look up coordinates in row-dominant array. // Have to handle right boundary specially -- if at the edge, we'd // better not reference the out of bounds point. // else { xll = xg[ul][vl]; yll = yg[ul][vl]; // ur is out of bounds if ( ur == nx && vr < ny ) { xlr = xg[ul][vr]; ylr = yg[ul][vr]; *tx = xll * ( 1 - dv ) + xlr * ( dv ); *ty = yll * ( 1 - dv ) + ylr * ( dv ); } // vr is out of bounds else if ( ur < nx && vr == ny ) { xrl = xg[ur][vl]; yrl = yg[ur][vl]; *tx = xll * ( 1 - du ) + xrl * ( du ); *ty = yll * ( 1 - du ) + yrl * ( du ); } // both ur and vr are out of bounds else if ( ur == nx && vr == ny ) { *tx = xll; *ty = yll; } // everything in bounds else { xrl = xg[ur][vl]; xlr = xg[ul][vr]; xrr = xg[ur][vr]; yrl = yg[ur][vl]; ylr = yg[ul][vr]; yrr = yg[ur][vr]; *tx = xll * ( 1 - du ) * ( 1 - dv ) + xlr * ( 1 - du ) * ( dv ) + xrl * ( du ) * ( 1 - dv ) + xrr * ( du ) * ( dv ); *ty = yll * ( 1 - du ) * ( 1 - dv ) + ylr * ( 1 - du ) * ( dv ) + yrl * ( du ) * ( 1 - dv ) + yrr * ( du ) * ( dv ); } } } //-------------------------------------------------------------------------- // pltr2p() // // Just like pltr2() but uses pointer arithmetic to get coordinates from 2d // grid tables. This form of grid tables is compatible with those from // PLplot 4.0. The grid data must be pointed to by a PLcGrid structure. //-------------------------------------------------------------------------- void pltr2p( PLFLT x, PLFLT y, PLFLT *tx, PLFLT *ty, PLPointer pltr_data ) { PLINT ul, ur, vl, vr; PLFLT du, dv; PLFLT xll, xlr, xrl, xrr; PLFLT yll, ylr, yrl, yrr; PLFLT xmin, xmax, ymin, ymax; PLcGrid *grid = (PLcGrid *) pltr_data; PLFLT *xg = grid->xg; PLFLT *yg = grid->yg; PLINT nx = grid->nx; PLINT ny = grid->ny; ul = (PLINT) x; ur = ul + 1; du = x - ul; vl = (PLINT) y; vr = vl + 1; dv = y - vl; xmin = 0; xmax = nx - 1; ymin = 0; ymax = ny - 1; if ( x < xmin || x > xmax || y < ymin || y > ymax ) { plwarn( "pltr2p: Invalid coordinates" ); if ( x < xmin ) { if ( y < ymin ) { *tx = *xg; *ty = *yg; } else if ( y > ymax ) { *tx = *( xg + ( ny - 1 ) ); *ty = *( yg + ( ny - 1 ) ); } else { ul = 0; xll = *( xg + ul * ny + vl ); yll = *( yg + ul * ny + vl ); xlr = *( xg + ul * ny + vr ); ylr = *( yg + ul * ny + vr ); *tx = xll * ( 1 - dv ) + xlr * ( dv ); *ty = yll * ( 1 - dv ) + ylr * ( dv ); } } else if ( x > xmax ) { if ( y < ymin ) { *tx = *( xg + ( ny - 1 ) * nx ); *ty = *( yg + ( ny - 1 ) * nx ); } else if ( y > ymax ) { *tx = *( xg + ( ny - 1 ) + ( nx - 1 ) * ny ); *ty = *( yg + ( ny - 1 ) + ( nx - 1 ) * ny ); } else { ul = nx - 1; xll = *( xg + ul * ny + vl ); yll = *( yg + ul * ny + vl ); xlr = *( xg + ul * ny + vr ); ylr = *( yg + ul * ny + vr ); *tx = xll * ( 1 - dv ) + xlr * ( dv ); *ty = yll * ( 1 - dv ) + ylr * ( dv ); } } else { if ( y < ymin ) { vl = 0; xll = *( xg + ul * ny + vl ); xrl = *( xg + ur * ny + vl ); yll = *( yg + ul * ny + vl ); yrl = *( yg + ur * ny + vl ); *tx = xll * ( 1 - du ) + xrl * ( du ); *ty = yll * ( 1 - du ) + yrl * ( du ); } else if ( y > ymax ) { vr = ny - 1; xlr = *( xg + ul * ny + vr ); xrr = *( xg + ur * ny + vr ); ylr = *( yg + ul * ny + vr ); yrr = *( yg + ur * ny + vr ); *tx = xlr * ( 1 - du ) + xrr * ( du ); *ty = ylr * ( 1 - du ) + yrr * ( du ); } } } // Normal case. // Look up coordinates in row-dominant array. // Have to handle right boundary specially -- if at the edge, we'd better // not reference the out of bounds point. // else { xll = *( xg + ul * ny + vl ); yll = *( yg + ul * ny + vl ); // ur is out of bounds if ( ur == nx && vr < ny ) { xlr = *( xg + ul * ny + vr ); ylr = *( yg + ul * ny + vr ); *tx = xll * ( 1 - dv ) + xlr * ( dv ); *ty = yll * ( 1 - dv ) + ylr * ( dv ); } // vr is out of bounds else if ( ur < nx && vr == ny ) { xrl = *( xg + ur * ny + vl ); yrl = *( yg + ur * ny + vl ); *tx = xll * ( 1 - du ) + xrl * ( du ); *ty = yll * ( 1 - du ) + yrl * ( du ); } // both ur and vr are out of bounds else if ( ur == nx && vr == ny ) { *tx = xll; *ty = yll; } // everything in bounds else { xrl = *( xg + ur * ny + vl ); xlr = *( xg + ul * ny + vr ); xrr = *( xg + ur * ny + vr ); yrl = *( yg + ur * ny + vl ); ylr = *( yg + ul * ny + vr ); yrr = *( yg + ur * ny + vr ); *tx = xll * ( 1 - du ) * ( 1 - dv ) + xlr * ( 1 - du ) * ( dv ) + xrl * ( du ) * ( 1 - dv ) + xrr * ( du ) * ( dv ); *ty = yll * ( 1 - du ) * ( 1 - dv ) + ylr * ( 1 - du ) * ( dv ) + yrl * ( du ) * ( 1 - dv ) + yrr * ( du ) * ( dv ); } } } //-------------------------------------------------------------------------- // pltr2f() // // Does linear interpolation from doubly dimensioned coord arrays // (row dominant, i.e. Fortran ordering). // // This routine includes lots of checks for out of bounds. This would // occur occasionally due to a bug in the contour plotter that is now fixed. // If an out of bounds coordinate is obtained, the boundary value is provided // along with a warning. These checks should stay since no harm is done if // if everything works correctly. //-------------------------------------------------------------------------- void pltr2f( PLFLT x, PLFLT y, PLFLT *tx, PLFLT *ty, void *pltr_data ) { PLINT ul, ur, vl, vr; PLFLT du, dv; PLFLT xll, xlr, xrl, xrr; PLFLT yll, ylr, yrl, yrr; PLFLT xmin, xmax, ymin, ymax; PLcGrid *cgrid = (PLcGrid *) pltr_data; PLFLT *xg = cgrid->xg; PLFLT *yg = cgrid->yg; PLINT nx = cgrid->nx; PLINT ny = cgrid->ny; ul = (PLINT) x; ur = ul + 1; du = x - ul; vl = (PLINT) y; vr = vl + 1; dv = y - vl; xmin = 0; xmax = nx - 1; ymin = 0; ymax = ny - 1; if ( x < xmin || x > xmax || y < ymin || y > ymax ) { plwarn( "pltr2f: Invalid coordinates" ); if ( x < xmin ) { if ( y < ymin ) { *tx = *xg; *ty = *yg; } else if ( y > ymax ) { *tx = *( xg + ( ny - 1 ) * nx ); *ty = *( yg + ( ny - 1 ) * nx ); } else { ul = 0; xll = *( xg + ul + vl * nx ); yll = *( yg + ul + vl * nx ); xlr = *( xg + ul + vr * nx ); ylr = *( yg + ul + vr * nx ); *tx = xll * ( 1 - dv ) + xlr * ( dv ); *ty = yll * ( 1 - dv ) + ylr * ( dv ); } } else if ( x > xmax ) { if ( y < ymin ) { *tx = *( xg + ( nx - 1 ) ); *ty = *( yg + ( nx - 1 ) ); } else if ( y > ymax ) { *tx = *( xg + ( nx - 1 ) + ( ny - 1 ) * nx ); *ty = *( yg + ( nx - 1 ) + ( ny - 1 ) * nx ); } else { ul = nx - 1; xll = *( xg + ul + vl * nx ); yll = *( yg + ul + vl * nx ); xlr = *( xg + ul + vr * nx ); ylr = *( yg + ul + vr * nx ); *tx = xll * ( 1 - dv ) + xlr * ( dv ); *ty = yll * ( 1 - dv ) + ylr * ( dv ); } } else { if ( y < ymin ) { vl = 0; xll = *( xg + ul + vl * nx ); xrl = *( xg + ur + vl * nx ); yll = *( yg + ul + vl * nx ); yrl = *( yg + ur + vl * nx ); *tx = xll * ( 1 - du ) + xrl * ( du ); *ty = yll * ( 1 - du ) + yrl * ( du ); } else if ( y > ymax ) { vr = ny - 1; xlr = *( xg + ul + vr * nx ); xrr = *( xg + ur + vr * nx ); ylr = *( yg + ul + vr * nx ); yrr = *( yg + ur + vr * nx ); *tx = xlr * ( 1 - du ) + xrr * ( du ); *ty = ylr * ( 1 - du ) + yrr * ( du ); } } } // Normal case. // Look up coordinates in row-dominant array. // Have to handle right boundary specially -- if at the edge, we'd // better not reference the out of bounds point. else { xll = *( xg + ul + vl * nx ); yll = *( yg + ul + vl * nx ); // ur is out of bounds if ( ur == nx && vr < ny ) { xlr = *( xg + ul + vr * nx ); ylr = *( yg + ul + vr * nx ); *tx = xll * ( 1 - dv ) + xlr * ( dv ); *ty = yll * ( 1 - dv ) + ylr * ( dv ); } // vr is out of bounds else if ( ur < nx && vr == ny ) { xrl = *( xg + ur + vl * nx ); yrl = *( yg + ur + vl * nx ); *tx = xll * ( 1 - du ) + xrl * ( du ); *ty = yll * ( 1 - du ) + yrl * ( du ); } // both ur and vr are out of bounds else if ( ur == nx && vr == ny ) { *tx = xll; *ty = yll; } // everything in bounds else { xrl = *( xg + ur + vl * nx ); xlr = *( xg + ul + vr * nx ); xrr = *( xg + ur + vr * nx ); yrl = *( yg + ur + vl * nx ); ylr = *( yg + ul + vr * nx ); yrr = *( yg + ur + vr * nx ); *tx = xll * ( 1 - du ) * ( 1 - dv ) + xlr * ( 1 - du ) * ( dv ) + xrl * ( du ) * ( 1 - dv ) + xrr * ( du ) * ( dv ); *ty = yll * ( 1 - du ) * ( 1 - dv ) + ylr * ( 1 - du ) * ( dv ) + yrl * ( du ) * ( 1 - dv ) + yrr * ( du ) * ( dv ); } } }