PLplot  5.10.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros
plcont.c
Go to the documentation of this file.
1 // Contour plotter.
2 //
3 // Copyright (C) 1995, 2000, 2001 Maurice LeBrun
4 // Copyright (C) 2000, 2002 Joao Cardoso
5 // Copyright (C) 2000-2014 Alan W. Irwin
6 // Copyright (C) 2004 Andrew Ross
7 //
8 // This file is part of PLplot.
9 //
10 // PLplot is free software; you can redistribute it and/or modify
11 // it under the terms of the GNU Library General Public License as published
12 // by the Free Software Foundation; either version 2 of the License, or
13 // (at your option) any later version.
14 //
15 // PLplot is distributed in the hope that it will be useful,
16 // but WITHOUT ANY WARRANTY; without even the implied warranty of
17 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18 // GNU Library General Public License for more details.
19 //
20 // You should have received a copy of the GNU Library General Public License
21 // along with PLplot; if not, write to the Free Software
22 // Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
23 //
24 
25 #include "plplotP.h"
26 
27 #ifdef MSDOS
28 #pragma optimize("",off)
29 #endif
30 
31 // Static function prototypes.
32 
33 static void
35  PLPointer plf2eval_data,
36  PLINT nx, PLINT ny, PLINT kx, PLINT lx,
37  PLINT ky, PLINT ly, PLFLT flev, PLINT **ipts,
38  void ( *pltr )( PLFLT, PLFLT, PLFLT *, PLFLT *, PLPointer ),
39  PLPointer pltr_data );
40 
41 static void
43  PLPointer plf2eval_data,
44  PLINT nx, PLINT ny, PLINT kx, PLINT lx,
45  PLINT ky, PLINT ly, PLFLT flev, char *flabel, PLINT kcol, PLINT krow,
46  PLFLT lastx, PLFLT lasty, PLINT startedge,
47  PLINT **ipts, PLFLT *distance, PLINT *lastindex,
48  void ( *pltr )( PLFLT, PLFLT, PLFLT *, PLFLT *, PLPointer ),
49  PLPointer pltr_data );
50 
51 static void
52 plfloatlabel( PLFLT value, char *string, PLINT len );
53 
54 static PLFLT
55 plP_pcwcx( PLINT x );
56 
57 static PLFLT
58 plP_pcwcy( PLINT y );
59 
60 static void
61 pl_drawcontlabel( PLFLT tpx, PLFLT tpy, char *flabel, PLFLT *distance, PLINT *lastindex );
62 
63 // Error flag for aborts
64 
65 static int error;
66 
67 //**************************************
68 //
69 // Defaults for contour label printing.
70 //
71 //**************************************
72 
73 // Font height for contour labels (normalized)
74 static PLFLT
76 
77 // Offset of label from contour line (if set to 0.0, labels are printed on the lines).
78 static PLFLT
80 
81 // Spacing parameter for contour labels
82 static PLFLT
84 
85 // Activate labels, default off
86 static PLINT
88 
89 // If the contour label exceed 10^(limexp) or 10^(-limexp), the exponential format is used
90 static PLINT
91  limexp = 4;
92 
93 // Number of significant digits
94 static PLINT
95  sigprec = 2;
96 
97 //******* contour lines storage ***************************
98 
99 static CONT_LEVEL *startlev = NULL;
102 
103 static int cont3d = 0;
104 
105 static CONT_LINE *
106 alloc_line( void )
107 {
108  CONT_LINE *line;
109 
110  if ( ( line = (CONT_LINE *) malloc( sizeof ( CONT_LINE ) ) ) == NULL )
111  {
112  plexit( "alloc_line: Insufficient memory" );
113  }
114 
115  line->x = (PLFLT *) malloc( LINE_ITEMS * sizeof ( PLFLT ) );
116  line->y = (PLFLT *) malloc( LINE_ITEMS * sizeof ( PLFLT ) );
117 
118  if ( ( line->x == NULL ) || ( line->y == NULL ) )
119  {
120  plexit( "alloc_line: Insufficient memory" );
121  }
122 
123  line->npts = 0;
124  line->next = NULL;
125 
126  return line;
127 }
128 
129 static CONT_LEVEL *
131 {
132  CONT_LEVEL *node;
133 
134  if ( ( node = (CONT_LEVEL *) malloc( sizeof ( CONT_LEVEL ) ) ) == NULL )
135  {
136  plexit( "alloc_level: Insufficient memory" );
137  }
138  node->level = level;
139  node->next = NULL;
140  node->line = alloc_line( );
141 
142  return node;
143 }
144 
145 static void
147 {
148  if ( ( ( line->x = (PLFLT *) realloc( line->x,
149  (size_t) ( line->npts + LINE_ITEMS ) * sizeof ( PLFLT ) ) ) == NULL ) ||
150  ( ( line->y = (PLFLT *) realloc( line->y,
151  (size_t) ( line->npts + LINE_ITEMS ) * sizeof ( PLFLT ) ) ) == NULL ) )
152  plexit( "realloc_line: Insufficient memory" );
153 }
154 
155 
156 // new contour level
157 static void
159 {
160  if ( cont3d )
161  {
162  if ( startlev == NULL )
163  {
164  startlev = alloc_level( level );
165  currlev = startlev;
166  }
167  else
168  {
169  currlev->next = alloc_level( level );
170  currlev = currlev->next;
171  }
172  currline = currlev->line;
173  }
174 }
175 
176 void
178 {
179  CONT_LINE *tline, *cline;
180  CONT_LEVEL *tlev, *clevel;
181 
182  if ( ct != NULL )
183  {
184  clevel = ct;
185 
186  do
187  {
188  cline = clevel->line;
189  do
190  {
191 #ifdef CONT_PLOT_DEBUG // for 2D plots. For 3D plots look at plot3.c:plotsh3di()
192  plP_movwor( cline->x[0], cline->y[0] );
193  for ( j = 1; j < cline->npts; j++ )
194  plP_drawor( cline->x[j], cline->y[j] );
195 #endif
196  tline = cline->next;
197  free( cline->x );
198  free( cline->y );
199  free( cline );
200  cline = tline;
201  }
202  while ( cline != NULL );
203  tlev = clevel->next;
204  free( clevel );
205  clevel = tlev;
206  }
207  while ( clevel != NULL );
208  startlev = NULL;
209  }
210 }
211 
212 static void
214 {
215  if ( cont3d )
216  {
217  PLINT pts = currline->npts;
218 
219  if ( pts % LINE_ITEMS == 0 )
220  realloc_line( currline );
221 
222  currline->x[pts] = xx;
223  currline->y[pts] = yy;
224  currline->npts++;
225  }
226  else
227  plP_drawor( xx, yy );
228 }
229 
230 static void
232 {
233  if ( cont3d )
234  {
235  if ( currline->npts != 0 ) // not an empty list, allocate new
236  {
237  currline->next = alloc_line( );
238  currline = currline->next;
239  }
240 
241  // and fill first element
242  currline->x[0] = xx;
243  currline->y[0] = yy;
244  currline->npts = 1;
245  }
246  else
247  plP_movwor( xx, yy );
248 }
249 
250 // small routine to set offset and spacing of contour labels, see desciption above
251 void c_pl_setcontlabelparam( PLFLT offset, PLFLT size, PLFLT spacing, PLINT active )
252 {
253  contlabel_offset = offset;
254  contlabel_size = size;
255  contlabel_space = spacing;
256  contlabel_active = active;
257 }
258 
259 // small routine to set the format of the contour labels, description of limexp and prec see above
260 void c_pl_setcontlabelformat( PLINT lexp, PLINT sigdig )
261 {
262  limexp = lexp;
263  sigprec = sigdig;
264 }
265 
266 static void pl_drawcontlabel( PLFLT tpx, PLFLT tpy, char *flabel, PLFLT *distance, PLINT *lastindex )
267 {
268  PLFLT delta_x, delta_y;
269  PLINT currx_old, curry_old;
270 
271  delta_x = plP_pcdcx( plsc->currx ) - plP_pcdcx( plP_wcpcx( tpx ) );
272  delta_y = plP_pcdcy( plsc->curry ) - plP_pcdcy( plP_wcpcy( tpy ) );
273 
274  currx_old = plsc->currx;
275  curry_old = plsc->curry;
276 
277  *distance += sqrt( delta_x * delta_x + delta_y * delta_y );
278 
279  plP_drawor( tpx, tpy );
280 
281  if ( (int) ( fabs( *distance / contlabel_space ) ) > *lastindex )
282  {
283  PLFLT scale, vec_x, vec_y, mx, my, dev_x, dev_y, off_x, off_y;
284 
285  vec_x = tpx - plP_pcwcx( currx_old );
286  vec_y = tpy - plP_pcwcy( curry_old );
287 
288  // Ensure labels appear the right way up
289  if ( vec_x < 0 )
290  {
291  vec_x = -vec_x;
292  vec_y = -vec_y;
293  }
294 
295  mx = (double) plsc->wpxscl / (double) plsc->phyxlen;
296  my = (double) plsc->wpyscl / (double) plsc->phyylen;
297 
298  dev_x = -my * vec_y / mx;
299  dev_y = mx * vec_x / my;
300 
301  scale = sqrt( ( mx * mx * dev_x * dev_x + my * my * dev_y * dev_y ) /
303 
304  off_x = dev_x / scale;
305  off_y = dev_y / scale;
306 
307  plptex( tpx + off_x, tpy + off_y, vec_x, vec_y, 0.5, flabel );
308  plP_movwor( tpx, tpy );
309  ( *lastindex )++;
310  }
311  else
312  plP_movwor( tpx, tpy );
313 }
314 
315 
316 // Format contour labels. Arguments:
317 // value: floating point number to be formatted
318 // string: the formatted label, plptex must be called with it to actually
319 // print the label
320 //
321 
322 static void plfloatlabel( PLFLT value, char *string, PLINT len )
323 {
324  PLINT setpre, precis;
325  // form[10] gives enough space for all non-malicious formats.
326  // tmpstring[15] gives enough room for 3 digits in a negative exponent
327  // or 4 digits in a positive exponent + null termination. That
328  // should be enough for all non-malicious use.
329  // Obviously there are security issues here that
330  // should be addressed as well.
331  //
332 #define FORM_LEN 10
333 #define TMPSTRING_LEN 15
334  char form[FORM_LEN], tmpstring[TMPSTRING_LEN];
335  PLINT exponent = 0;
336  PLFLT mant, tmp;
337 
338  PLINT prec = sigprec;
339 
340  plP_gprec( &setpre, &precis );
341 
342  if ( setpre )
343  prec = precis;
344 
345  if ( value > 0.0 )
346  tmp = log10( value );
347  else if ( value < 0.0 )
348  tmp = log10( -value );
349  else
350  tmp = 0;
351 
352  if ( tmp >= 0.0 )
353  exponent = (int) tmp;
354  else if ( tmp < 0.0 )
355  {
356  tmp = -tmp;
357  if ( floor( tmp ) < tmp )
358  exponent = -(int) ( floor( tmp ) + 1.0 );
359  else
360  exponent = -(int) ( floor( tmp ) );
361  }
362 
363  mant = value / pow( 10.0, exponent );
364 
365  if ( mant != 0.0 )
366  mant = (int) ( mant * pow( 10.0, prec - 1 ) + 0.5 * mant / fabs( mant ) ) / pow( 10.0, prec - 1 );
367 
368  snprintf( form, FORM_LEN, "%%.%df", prec - 1 );
369  snprintf( string, (size_t) len, form, mant );
370  snprintf( tmpstring, TMPSTRING_LEN, "#(229)10#u%d", exponent );
371  strncat( string, tmpstring, (size_t) len - strlen( string ) - 1 );
372 
373  if ( abs( exponent ) < limexp || value == 0.0 )
374  {
375  value = pow( 10.0, exponent ) * mant;
376 
377  if ( exponent >= 0 )
378  prec = prec - 1 - exponent;
379  else
380  prec = prec - 1 + abs( exponent );
381 
382  if ( prec < 0 )
383  prec = 0;
384 
385  snprintf( form, FORM_LEN, "%%.%df", (int) prec );
386  snprintf( string, (size_t) len, form, value );
387  }
388 }
389 
390 // physical coords (x) to world coords
391 
392 static PLFLT
394 {
395  return ( ( x - plsc->wpxoff ) / plsc->wpxscl );
396 }
397 
398 // physical coords (y) to world coords
399 
400 static PLFLT
402 {
403  return ( ( y - plsc->wpyoff ) / plsc->wpyscl );
404 }
405 
406 //--------------------------------------------------------------------------
407 // plf2eval1()
408 //
409 // Does a lookup from a 2d function array. Array is of type (PLFLT **),
410 // and is column dominant (normal C ordering).
411 //--------------------------------------------------------------------------
412 
413 PLFLT
414 plf2eval1( PLINT ix, PLINT iy, PLPointer plf2eval_data )
415 {
416  PLFLT value;
417  PLFLT **z = (PLFLT **) plf2eval_data;
418 
419  value = z[ix][iy];
420 
421  return value;
422 }
423 
424 //--------------------------------------------------------------------------
425 // plf2eval2()
426 //
427 // Does a lookup from a 2d function array. plf2eval_data is treated as type
428 // (PLfGrid2 *).
429 //--------------------------------------------------------------------------
430 
431 PLFLT
432 plf2eval2( PLINT ix, PLINT iy, PLPointer plf2eval_data )
433 {
434  PLFLT value;
435  PLfGrid2 *grid = (PLfGrid2 *) plf2eval_data;
436 
437  value = grid->f[ix][iy];
438 
439  return value;
440 }
441 
442 //--------------------------------------------------------------------------
443 // plf2eval()
444 //
445 // Does a lookup from a 2d function array. Array is of type (PLFLT *), and
446 // is column dominant (normal C ordering). You MUST fill the ny maximum
447 // array index entry in the PLfGrid struct.
448 //--------------------------------------------------------------------------
449 
450 PLFLT
451 plf2eval( PLINT ix, PLINT iy, PLPointer plf2eval_data )
452 {
453  PLFLT value;
454  PLfGrid *grid = (PLfGrid *) plf2eval_data;
455 
456  value = grid->f[ix * grid->ny + iy];
457 
458  return value;
459 }
460 
461 //--------------------------------------------------------------------------
462 // plf2evalr()
463 //
464 // Does a lookup from a 2d function array. Array is of type (PLFLT *), and
465 // is row dominant (Fortran ordering). You MUST fill the nx maximum array
466 // index entry in the PLfGrid struct.
467 //--------------------------------------------------------------------------
468 
469 PLFLT
470 plf2evalr( PLINT ix, PLINT iy, PLPointer plf2eval_data )
471 {
472  PLFLT value;
473  PLfGrid *grid = (PLfGrid *) plf2eval_data;
474 
475  value = grid->f[ix + iy * grid->nx];
476 
477  return value;
478 }
479 
480 //--------------------------------------------------------------------------
481 //
482 // cont_store:
483 //
484 // Draw contour lines in memory.
485 // cont_clean_store() must be called after use to release allocated memory.
486 //
487 //--------------------------------------------------------------------------
488 
489 void
490 cont_store( const PLFLT * const *f, PLINT nx, PLINT ny, PLINT kx, PLINT lx,
491  PLINT ky, PLINT ly, const PLFLT *clevel, PLINT nlevel,
492  void ( *pltr )( PLFLT, PLFLT, PLFLT *, PLFLT *, PLPointer ),
493  PLPointer pltr_data,
494  CONT_LEVEL **contour )
495 {
496  cont3d = 1;
497 
498  plcont( f, nx, ny, kx, lx, ky, ly, clevel, nlevel,
499  pltr, pltr_data );
500 
501  *contour = startlev;
502  cont3d = 0;
503 }
504 
505 //--------------------------------------------------------------------------
506 // void plcont()
507 //
508 // Draws a contour plot from data in f(nx,ny). Is just a front-end to
509 // plfcont, with a particular choice for f2eval and f2eval_data.
510 //--------------------------------------------------------------------------
511 
512 void
513 c_plcont( const PLFLT * const *f, PLINT nx, PLINT ny, PLINT kx, PLINT lx,
514  PLINT ky, PLINT ly, const PLFLT *clevel, PLINT nlevel,
515  void ( *pltr )( PLFLT, PLFLT, PLFLT *, PLFLT *, PLPointer ),
516  PLPointer pltr_data )
517 {
518  if ( pltr == NULL )
519  {
520  // If pltr is undefined, abort with an error.
521  plabort( "plcont: The pltr callback must be defined" );
522  return;
523  }
524 
526  nx, ny, kx, lx, ky, ly, clevel, nlevel,
527  pltr, pltr_data );
528 }
529 
530 //--------------------------------------------------------------------------
531 // void plfcont()
532 //
533 // Draws a contour plot using the function evaluator f2eval and data stored
534 // by way of the f2eval_data pointer. This allows arbitrary organizations
535 // of 2d array data to be used.
536 //
537 // The subrange of indices used for contouring is kx to lx in the x
538 // direction and from ky to ly in the y direction. The array of contour
539 // levels is clevel(nlevel), and "pltr" is the name of a function which
540 // transforms array indicies into world coordinates.
541 //
542 // Note that the fortran-like minimum and maximum indices (kx, lx, ky, ly)
543 // are translated into more C-like ones. I've only kept them as they are
544 // for the plcontf() argument list because of backward compatibility.
545 //--------------------------------------------------------------------------
546 
547 void
548 plfcont( PLFLT ( *f2eval )( PLINT, PLINT, PLPointer ),
549  PLPointer f2eval_data,
550  PLINT nx, PLINT ny, PLINT kx, PLINT lx,
551  PLINT ky, PLINT ly, const PLFLT *clevel, PLINT nlevel,
552  void ( *pltr )( PLFLT, PLFLT, PLFLT *, PLFLT *, PLPointer ),
553  PLPointer pltr_data )
554 {
555  PLINT i, **ipts;
556 
557  if ( kx < 1 || kx >= lx )
558  {
559  plabort( "plfcont: indices must satisfy 1 <= kx <= lx <= nx" );
560  return;
561  }
562  if ( ky < 1 || ky >= ly )
563  {
564  plabort( "plfcont: indices must satisfy 1 <= ky <= ly <= ny" );
565  return;
566  }
567 
568  if ( ( ipts = (PLINT **) malloc( (size_t) nx * sizeof ( PLINT * ) ) ) == NULL )
569  {
570  plexit( "plfcont: Insufficient memory" );
571  }
572 
573  for ( i = 0; i < nx; i++ )
574  {
575  if ( ( ipts[i] = (PLINT *) malloc( (size_t) ny * sizeof ( PLINT * ) ) ) == NULL )
576  {
577  plexit( "plfcont: Insufficient memory" );
578  }
579  }
580 
581  for ( i = 0; i < nlevel; i++ )
582  {
583  plcntr( f2eval, f2eval_data,
584  nx, ny, kx - 1, lx - 1, ky - 1, ly - 1, clevel[i], ipts,
585  pltr, pltr_data );
586 
587  if ( error )
588  {
589  error = 0;
590  goto done;
591  }
592  }
593 
594 done:
595  for ( i = 0; i < nx; i++ )
596  {
597  free( (void *) ipts[i] );
598  }
599  free( (void *) ipts );
600 }
601 
602 //--------------------------------------------------------------------------
603 // void plcntr()
604 //
605 // The contour for a given level is drawn here. Note iscan has nx
606 // elements. ixstor and iystor each have nstor elements.
607 //--------------------------------------------------------------------------
608 
609 static void
610 plcntr( PLFLT ( *f2eval )( PLINT, PLINT, PLPointer ),
611  PLPointer f2eval_data,
612  PLINT nx, PLINT ny, PLINT kx, PLINT lx,
613  PLINT ky, PLINT ly, PLFLT flev, PLINT **ipts,
614  void ( *pltr )( PLFLT, PLFLT, PLFLT *, PLFLT *, PLPointer ),
615  PLPointer pltr_data )
616 {
617  PLINT kcol, krow, lastindex;
618  PLFLT distance;
619  PLFLT save_def, save_scale;
620 
621  char flabel[30];
622  plgchr( &save_def, &save_scale );
623  save_scale = save_scale / save_def;
624 
625  cont_new_store( flev );
626 
627  // format contour label for plptex and define the font height of the labels
628  plfloatlabel( flev, flabel, 30 );
629  plschr( 0.0, contlabel_size );
630 
631  // Clear array for traversed squares
632  for ( kcol = kx; kcol < lx; kcol++ )
633  {
634  for ( krow = ky; krow < ly; krow++ )
635  {
636  ipts[kcol][krow] = 0;
637  }
638  }
639 
640 
641  for ( krow = ky; krow < ly; krow++ )
642  {
643  for ( kcol = kx; kcol < lx; kcol++ )
644  {
645  if ( ipts[kcol][krow] == 0 )
646  {
647  // Follow and draw a contour
648  pldrawcn( f2eval, f2eval_data,
649  nx, ny, kx, lx, ky, ly, flev, flabel, kcol, krow,
650  0.0, 0.0, -2, ipts, &distance, &lastindex,
651  pltr, pltr_data );
652 
653  if ( error )
654  return;
655  }
656  }
657  }
658  plschr( save_def, save_scale );
659 }
660 
661 //--------------------------------------------------------------------------
662 // void pldrawcn()
663 //
664 // Follow and draw a contour.
665 //--------------------------------------------------------------------------
666 
667 static void
668 pldrawcn( PLFLT ( *f2eval )( PLINT, PLINT, PLPointer ),
669  PLPointer f2eval_data,
670  PLINT nx, PLINT ny, PLINT kx, PLINT lx,
671  PLINT ky, PLINT ly, PLFLT flev, char *flabel, PLINT kcol, PLINT krow,
672  PLFLT lastx, PLFLT lasty, PLINT startedge, PLINT **ipts,
673  PLFLT *distance, PLINT *lastindex,
674  void ( *pltr )( PLFLT, PLFLT, PLFLT *, PLFLT *, PLPointer ),
675  PLPointer pltr_data )
676 {
677  PLFLT f[4];
678  PLFLT px[4], py[4], locx[4], locy[4];
679  PLINT iedge[4];
680  PLINT i, j, k, num, first, inext, kcolnext, krownext;
681 
682 
683  ( *pltr )( kcol, krow + 1, &px[0], &py[0], pltr_data );
684  ( *pltr )( kcol, krow, &px[1], &py[1], pltr_data );
685  ( *pltr )( kcol + 1, krow, &px[2], &py[2], pltr_data );
686  ( *pltr )( kcol + 1, krow + 1, &px[3], &py[3], pltr_data );
687 
688  f[0] = f2eval( kcol, krow + 1, f2eval_data ) - flev;
689  f[1] = f2eval( kcol, krow, f2eval_data ) - flev;
690  f[2] = f2eval( kcol + 1, krow, f2eval_data ) - flev;
691  f[3] = f2eval( kcol + 1, krow + 1, f2eval_data ) - flev;
692 
693  for ( i = 0, j = 1; i < 4; i++, j = ( j + 1 ) % 4 )
694  {
695  iedge[i] = ( f[i] * f[j] > 0.0 ) ? -1 : ( ( f[i] * f[j] < 0.0 ) ? 1 : 0 );
696  }
697 
698  // Mark this square as done
699  ipts[kcol][krow] = 1;
700 
701  // Check if no contour has been crossed i.e. iedge[i] = -1
702  if ( ( iedge[0] == -1 ) && ( iedge[1] == -1 ) && ( iedge[2] == -1 )
703  && ( iedge[3] == -1 ) )
704  return;
705 
706  // Check if this is a completely flat square - in which case
707  // ignore it
708  if ( ( f[0] == 0.0 ) && ( f[1] == 0.0 ) && ( f[2] == 0.0 ) &&
709  ( f[3] == 0.0 ) )
710  return;
711 
712  // Calculate intersection points
713  num = 0;
714  if ( startedge < 0 )
715  {
716  first = 1;
717  }
718  else
719  {
720  locx[num] = lastx;
721  locy[num] = lasty;
722  num++;
723  first = 0;
724  }
725  for ( k = 0, i = ( startedge < 0 ? 0 : startedge ); k < 4; k++, i = ( i + 1 ) % 4 )
726  {
727  if ( i == startedge )
728  continue;
729 
730  // If the contour is an edge check it hasn't already been done
731  if ( f[i] == 0.0 && f[( i + 1 ) % 4] == 0.0 )
732  {
733  kcolnext = kcol;
734  krownext = krow;
735  if ( i == 0 )
736  kcolnext--;
737  if ( i == 1 )
738  krownext--;
739  if ( i == 2 )
740  kcolnext++;
741  if ( i == 3 )
742  krownext++;
743  if ( ( kcolnext < kx ) || ( kcolnext >= lx ) ||
744  ( krownext < ky ) || ( krownext >= ly ) ||
745  ( ipts[kcolnext][krownext] == 1 ) )
746  continue;
747  }
748  if ( ( iedge[i] == 1 ) || ( f[i] == 0.0 ) )
749  {
750  j = ( i + 1 ) % 4;
751  if ( f[i] != 0.0 )
752  {
753  locx[num] = ( px[i] * fabs( f[j] ) + px[j] * fabs( f[i] ) ) / fabs( f[j] - f[i] );
754  locy[num] = ( py[i] * fabs( f[j] ) + py[j] * fabs( f[i] ) ) / fabs( f[j] - f[i] );
755  }
756  else
757  {
758  locx[num] = px[i];
759  locy[num] = py[i];
760  }
761  // If this is the start of the contour then move to the point
762  if ( first == 1 )
763  {
764  cont_mv_store( locx[num], locy[num] );
765  first = 0;
766  *distance = 0;
767  *lastindex = 0;
768  }
769  else
770  {
771  // Link to the next point on the contour
772  if ( contlabel_active )
773  pl_drawcontlabel( locx[num], locy[num], flabel, distance, lastindex );
774  else
775  cont_xy_store( locx[num], locy[num] );
776  // Need to follow contour into next grid box
777  // Easy case where contour does not pass through corner
778  if ( f[i] != 0.0 )
779  {
780  kcolnext = kcol;
781  krownext = krow;
782  inext = ( i + 2 ) % 4;
783  if ( i == 0 )
784  kcolnext--;
785  if ( i == 1 )
786  krownext--;
787  if ( i == 2 )
788  kcolnext++;
789  if ( i == 3 )
790  krownext++;
791  if ( ( kcolnext >= kx ) && ( kcolnext < lx ) &&
792  ( krownext >= ky ) && ( krownext < ly ) &&
793  ( ipts[kcolnext][krownext] == 0 ) )
794  {
795  pldrawcn( f2eval, f2eval_data,
796  nx, ny, kx, lx, ky, ly, flev, flabel,
797  kcolnext, krownext,
798  locx[num], locy[num], inext, ipts,
799  distance, lastindex,
800  pltr, pltr_data );
801  }
802  }
803  // Hard case where contour passes through corner
804  // This is still not perfect - it may lose the contour
805  // which won't upset the contour itself (we can find it
806  // again later) but might upset the labelling
807  else
808  {
809  kcolnext = kcol;
810  krownext = krow;
811  inext = ( i + 2 ) % 4;
812  if ( i == 0 )
813  {
814  kcolnext--; krownext++;
815  }
816  if ( i == 1 )
817  {
818  krownext--; kcolnext--;
819  }
820  if ( i == 2 )
821  {
822  kcolnext++; krownext--;
823  }
824  if ( i == 3 )
825  {
826  krownext++; kcolnext++;
827  }
828  if ( ( kcolnext >= kx ) && ( kcolnext < lx ) &&
829  ( krownext >= ky ) && ( krownext < ly ) &&
830  ( ipts[kcolnext][krownext] == 0 ) )
831  {
832  pldrawcn( f2eval, f2eval_data,
833  nx, ny, kx, lx, ky, ly, flev, flabel,
834  kcolnext, krownext,
835  locx[num], locy[num], inext, ipts,
836  distance, lastindex,
837  pltr, pltr_data );
838  }
839  }
840  if ( first == 1 )
841  {
842  // Move back to first point
843  cont_mv_store( locx[num], locy[num] );
844  first = 0;
845  *distance = 0;
846  *lastindex = 0;
847  first = 0;
848  }
849  else
850  {
851  first = 1;
852  }
853  num++;
854  }
855  }
856  }
857 }
858 
859 //--------------------------------------------------------------------------
860 // pltr0()
861 //
862 // Identity transformation.
863 //--------------------------------------------------------------------------
864 
865 void
866 pltr0( PLFLT x, PLFLT y, PLFLT *tx, PLFLT *ty, PLPointer PL_UNUSED( pltr_data ) )
867 {
868  *tx = x;
869  *ty = y;
870 }
871 
872 //--------------------------------------------------------------------------
873 // pltr1()
874 //
875 // Does linear interpolation from singly dimensioned coord arrays.
876 //
877 // Just abort for now if coordinates are out of bounds (don't think it's
878 // possible, but if so we could use linear extrapolation).
879 //--------------------------------------------------------------------------
880 
881 void
882 pltr1( PLFLT x, PLFLT y, PLFLT *tx, PLFLT *ty, PLPointer pltr_data )
883 {
884  PLINT ul, ur, vl, vr;
885  PLFLT du, dv;
886  PLFLT xl, xr, yl, yr;
887 
888  PLcGrid *grid = (PLcGrid *) pltr_data;
889  PLFLT *xg = grid->xg;
890  PLFLT *yg = grid->yg;
891  PLINT nx = grid->nx;
892  PLINT ny = grid->ny;
893 
894  ul = (PLINT) x;
895  ur = ul + 1;
896  du = x - ul;
897 
898  vl = (PLINT) y;
899  vr = vl + 1;
900  dv = y - vl;
901 
902  if ( x < 0 || x > nx - 1 || y < 0 || y > ny - 1 )
903  {
904  plexit( "pltr1: Invalid coordinates" );
905  }
906 
907 // Look up coordinates in row-dominant array.
908 // Have to handle right boundary specially -- if at the edge, we'd better
909 // not reference the out of bounds point.
910 //
911 
912  xl = xg[ul];
913  yl = yg[vl];
914 
915  if ( ur == nx )
916  {
917  *tx = xl;
918  }
919  else
920  {
921  xr = xg[ur];
922  *tx = xl * ( 1 - du ) + xr * du;
923  }
924  if ( vr == ny )
925  {
926  *ty = yl;
927  }
928  else
929  {
930  yr = yg[vr];
931  *ty = yl * ( 1 - dv ) + yr * dv;
932  }
933 }
934 
935 //--------------------------------------------------------------------------
936 // pltr2()
937 //
938 // Does linear interpolation from doubly dimensioned coord arrays (column
939 // dominant, as per normal C 2d arrays).
940 //
941 // This routine includes lots of checks for out of bounds. This would occur
942 // occasionally due to some bugs in the contour plotter (now fixed). If an
943 // out of bounds coordinate is obtained, the boundary value is provided
944 // along with a warning. These checks should stay since no harm is done if
945 // if everything works correctly.
946 //--------------------------------------------------------------------------
947 
948 void
949 pltr2( PLFLT x, PLFLT y, PLFLT *tx, PLFLT *ty, PLPointer pltr_data )
950 {
951  PLINT ul, ur, vl, vr;
952  PLFLT du, dv;
953  PLFLT xll, xlr, xrl, xrr;
954  PLFLT yll, ylr, yrl, yrr;
955  PLFLT xmin, xmax, ymin, ymax;
956 
957  PLcGrid2 *grid = (PLcGrid2 *) pltr_data;
958  PLFLT **xg = grid->xg;
959  PLFLT **yg = grid->yg;
960  PLINT nx = grid->nx;
961  PLINT ny = grid->ny;
962 
963  ul = (PLINT) x;
964  ur = ul + 1;
965  du = x - ul;
966 
967  vl = (PLINT) y;
968  vr = vl + 1;
969  dv = y - vl;
970 
971  xmin = 0;
972  xmax = nx - 1;
973  ymin = 0;
974  ymax = ny - 1;
975 
976  if ( x < xmin || x > xmax || y < ymin || y > ymax )
977  {
978  plwarn( "pltr2: Invalid coordinates" );
979  if ( x < xmin )
980  {
981  if ( y < ymin )
982  {
983  *tx = xg[0][0];
984  *ty = yg[0][0];
985  }
986  else if ( y > ymax )
987  {
988  *tx = xg[0][ny - 1];
989  *ty = yg[0][ny - 1];
990  }
991  else
992  {
993  xll = xg[0][vl];
994  yll = yg[0][vl];
995  xlr = xg[0][vr];
996  ylr = yg[0][vr];
997 
998  *tx = xll * ( 1 - dv ) + xlr * ( dv );
999  *ty = yll * ( 1 - dv ) + ylr * ( dv );
1000  }
1001  }
1002  else if ( x > xmax )
1003  {
1004  if ( y < ymin )
1005  {
1006  *tx = xg[nx - 1][0];
1007  *ty = yg[nx - 1][0];
1008  }
1009  else if ( y > ymax )
1010  {
1011  *tx = xg[nx - 1][ny - 1];
1012  *ty = yg[nx - 1][ny - 1];
1013  }
1014  else
1015  {
1016  xll = xg[nx - 1][vl];
1017  yll = yg[nx - 1][vl];
1018  xlr = xg[nx - 1][vr];
1019  ylr = yg[nx - 1][vr];
1020 
1021  *tx = xll * ( 1 - dv ) + xlr * ( dv );
1022  *ty = yll * ( 1 - dv ) + ylr * ( dv );
1023  }
1024  }
1025  else
1026  {
1027  if ( y < ymin )
1028  {
1029  xll = xg[ul][0];
1030  xrl = xg[ur][0];
1031  yll = yg[ul][0];
1032  yrl = yg[ur][0];
1033 
1034  *tx = xll * ( 1 - du ) + xrl * ( du );
1035  *ty = yll * ( 1 - du ) + yrl * ( du );
1036  }
1037  else if ( y > ymax )
1038  {
1039  xlr = xg[ul][ny - 1];
1040  xrr = xg[ur][ny - 1];
1041  ylr = yg[ul][ny - 1];
1042  yrr = yg[ur][ny - 1];
1043 
1044  *tx = xlr * ( 1 - du ) + xrr * ( du );
1045  *ty = ylr * ( 1 - du ) + yrr * ( du );
1046  }
1047  }
1048  }
1049 
1050 // Normal case.
1051 // Look up coordinates in row-dominant array.
1052 // Have to handle right boundary specially -- if at the edge, we'd
1053 // better not reference the out of bounds point.
1054 //
1055 
1056  else
1057  {
1058  xll = xg[ul][vl];
1059  yll = yg[ul][vl];
1060 
1061  // ur is out of bounds
1062 
1063  if ( ur == nx && vr < ny )
1064  {
1065  xlr = xg[ul][vr];
1066  ylr = yg[ul][vr];
1067 
1068  *tx = xll * ( 1 - dv ) + xlr * ( dv );
1069  *ty = yll * ( 1 - dv ) + ylr * ( dv );
1070  }
1071 
1072  // vr is out of bounds
1073 
1074  else if ( ur < nx && vr == ny )
1075  {
1076  xrl = xg[ur][vl];
1077  yrl = yg[ur][vl];
1078 
1079  *tx = xll * ( 1 - du ) + xrl * ( du );
1080  *ty = yll * ( 1 - du ) + yrl * ( du );
1081  }
1082 
1083  // both ur and vr are out of bounds
1084 
1085  else if ( ur == nx && vr == ny )
1086  {
1087  *tx = xll;
1088  *ty = yll;
1089  }
1090 
1091  // everything in bounds
1092 
1093  else
1094  {
1095  xrl = xg[ur][vl];
1096  xlr = xg[ul][vr];
1097  xrr = xg[ur][vr];
1098 
1099  yrl = yg[ur][vl];
1100  ylr = yg[ul][vr];
1101  yrr = yg[ur][vr];
1102 
1103  *tx = xll * ( 1 - du ) * ( 1 - dv ) + xlr * ( 1 - du ) * ( dv ) +
1104  xrl * ( du ) * ( 1 - dv ) + xrr * ( du ) * ( dv );
1105 
1106  *ty = yll * ( 1 - du ) * ( 1 - dv ) + ylr * ( 1 - du ) * ( dv ) +
1107  yrl * ( du ) * ( 1 - dv ) + yrr * ( du ) * ( dv );
1108  }
1109  }
1110 }
1111 
1112 //--------------------------------------------------------------------------
1113 // pltr2p()
1114 //
1115 // Just like pltr2() but uses pointer arithmetic to get coordinates from 2d
1116 // grid tables. This form of grid tables is compatible with those from
1117 // PLplot 4.0. The grid data must be pointed to by a PLcGrid structure.
1118 //--------------------------------------------------------------------------
1119 
1120 void
1121 pltr2p( PLFLT x, PLFLT y, PLFLT *tx, PLFLT *ty, PLPointer pltr_data )
1122 {
1123  PLINT ul, ur, vl, vr;
1124  PLFLT du, dv;
1125  PLFLT xll, xlr, xrl, xrr;
1126  PLFLT yll, ylr, yrl, yrr;
1127  PLFLT xmin, xmax, ymin, ymax;
1128 
1129  PLcGrid *grid = (PLcGrid *) pltr_data;
1130  PLFLT *xg = grid->xg;
1131  PLFLT *yg = grid->yg;
1132  PLINT nx = grid->nx;
1133  PLINT ny = grid->ny;
1134 
1135  ul = (PLINT) x;
1136  ur = ul + 1;
1137  du = x - ul;
1138 
1139  vl = (PLINT) y;
1140  vr = vl + 1;
1141  dv = y - vl;
1142 
1143  xmin = 0;
1144  xmax = nx - 1;
1145  ymin = 0;
1146  ymax = ny - 1;
1147 
1148  if ( x < xmin || x > xmax || y < ymin || y > ymax )
1149  {
1150  plwarn( "pltr2p: Invalid coordinates" );
1151  if ( x < xmin )
1152  {
1153  if ( y < ymin )
1154  {
1155  *tx = *xg;
1156  *ty = *yg;
1157  }
1158  else if ( y > ymax )
1159  {
1160  *tx = *( xg + ( ny - 1 ) );
1161  *ty = *( yg + ( ny - 1 ) );
1162  }
1163  else
1164  {
1165  ul = 0;
1166  xll = *( xg + ul * ny + vl );
1167  yll = *( yg + ul * ny + vl );
1168  xlr = *( xg + ul * ny + vr );
1169  ylr = *( yg + ul * ny + vr );
1170 
1171  *tx = xll * ( 1 - dv ) + xlr * ( dv );
1172  *ty = yll * ( 1 - dv ) + ylr * ( dv );
1173  }
1174  }
1175  else if ( x > xmax )
1176  {
1177  if ( y < ymin )
1178  {
1179  *tx = *( xg + ( ny - 1 ) * nx );
1180  *ty = *( yg + ( ny - 1 ) * nx );
1181  }
1182  else if ( y > ymax )
1183  {
1184  *tx = *( xg + ( ny - 1 ) + ( nx - 1 ) * ny );
1185  *ty = *( yg + ( ny - 1 ) + ( nx - 1 ) * ny );
1186  }
1187  else
1188  {
1189  ul = nx - 1;
1190  xll = *( xg + ul * ny + vl );
1191  yll = *( yg + ul * ny + vl );
1192  xlr = *( xg + ul * ny + vr );
1193  ylr = *( yg + ul * ny + vr );
1194 
1195  *tx = xll * ( 1 - dv ) + xlr * ( dv );
1196  *ty = yll * ( 1 - dv ) + ylr * ( dv );
1197  }
1198  }
1199  else
1200  {
1201  if ( y < ymin )
1202  {
1203  vl = 0;
1204  xll = *( xg + ul * ny + vl );
1205  xrl = *( xg + ur * ny + vl );
1206  yll = *( yg + ul * ny + vl );
1207  yrl = *( yg + ur * ny + vl );
1208 
1209  *tx = xll * ( 1 - du ) + xrl * ( du );
1210  *ty = yll * ( 1 - du ) + yrl * ( du );
1211  }
1212  else if ( y > ymax )
1213  {
1214  vr = ny - 1;
1215  xlr = *( xg + ul * ny + vr );
1216  xrr = *( xg + ur * ny + vr );
1217  ylr = *( yg + ul * ny + vr );
1218  yrr = *( yg + ur * ny + vr );
1219 
1220  *tx = xlr * ( 1 - du ) + xrr * ( du );
1221  *ty = ylr * ( 1 - du ) + yrr * ( du );
1222  }
1223  }
1224  }
1225 
1226 // Normal case.
1227 // Look up coordinates in row-dominant array.
1228 // Have to handle right boundary specially -- if at the edge, we'd better
1229 // not reference the out of bounds point.
1230 //
1231 
1232  else
1233  {
1234  xll = *( xg + ul * ny + vl );
1235  yll = *( yg + ul * ny + vl );
1236 
1237  // ur is out of bounds
1238 
1239  if ( ur == nx && vr < ny )
1240  {
1241  xlr = *( xg + ul * ny + vr );
1242  ylr = *( yg + ul * ny + vr );
1243 
1244  *tx = xll * ( 1 - dv ) + xlr * ( dv );
1245  *ty = yll * ( 1 - dv ) + ylr * ( dv );
1246  }
1247 
1248  // vr is out of bounds
1249 
1250  else if ( ur < nx && vr == ny )
1251  {
1252  xrl = *( xg + ur * ny + vl );
1253  yrl = *( yg + ur * ny + vl );
1254 
1255  *tx = xll * ( 1 - du ) + xrl * ( du );
1256  *ty = yll * ( 1 - du ) + yrl * ( du );
1257  }
1258 
1259  // both ur and vr are out of bounds
1260 
1261  else if ( ur == nx && vr == ny )
1262  {
1263  *tx = xll;
1264  *ty = yll;
1265  }
1266 
1267  // everything in bounds
1268 
1269  else
1270  {
1271  xrl = *( xg + ur * ny + vl );
1272  xlr = *( xg + ul * ny + vr );
1273  xrr = *( xg + ur * ny + vr );
1274 
1275  yrl = *( yg + ur * ny + vl );
1276  ylr = *( yg + ul * ny + vr );
1277  yrr = *( yg + ur * ny + vr );
1278 
1279  *tx = xll * ( 1 - du ) * ( 1 - dv ) + xlr * ( 1 - du ) * ( dv ) +
1280  xrl * ( du ) * ( 1 - dv ) + xrr * ( du ) * ( dv );
1281 
1282  *ty = yll * ( 1 - du ) * ( 1 - dv ) + ylr * ( 1 - du ) * ( dv ) +
1283  yrl * ( du ) * ( 1 - dv ) + yrr * ( du ) * ( dv );
1284  }
1285  }
1286 }
1287 
1288 //--------------------------------------------------------------------------
1289 // pltr2f()
1290 //
1291 // Does linear interpolation from doubly dimensioned coord arrays
1292 // (row dominant, i.e. Fortran ordering).
1293 //
1294 // This routine includes lots of checks for out of bounds. This would
1295 // occur occasionally due to a bug in the contour plotter that is now fixed.
1296 // If an out of bounds coordinate is obtained, the boundary value is provided
1297 // along with a warning. These checks should stay since no harm is done if
1298 // if everything works correctly.
1299 //--------------------------------------------------------------------------
1300 
1301 void
1302 pltr2f( PLFLT x, PLFLT y, PLFLT *tx, PLFLT *ty, void *pltr_data )
1303 {
1304  PLINT ul, ur, vl, vr;
1305  PLFLT du, dv;
1306  PLFLT xll, xlr, xrl, xrr;
1307  PLFLT yll, ylr, yrl, yrr;
1308  PLFLT xmin, xmax, ymin, ymax;
1309 
1310  PLcGrid *cgrid = (PLcGrid *) pltr_data;
1311  PLFLT *xg = cgrid->xg;
1312  PLFLT *yg = cgrid->yg;
1313  PLINT nx = cgrid->nx;
1314  PLINT ny = cgrid->ny;
1315 
1316  ul = (PLINT) x;
1317  ur = ul + 1;
1318  du = x - ul;
1319 
1320  vl = (PLINT) y;
1321  vr = vl + 1;
1322  dv = y - vl;
1323 
1324  xmin = 0;
1325  xmax = nx - 1;
1326  ymin = 0;
1327  ymax = ny - 1;
1328 
1329  if ( x < xmin || x > xmax || y < ymin || y > ymax )
1330  {
1331  plwarn( "pltr2f: Invalid coordinates" );
1332 
1333  if ( x < xmin )
1334  {
1335  if ( y < ymin )
1336  {
1337  *tx = *xg;
1338  *ty = *yg;
1339  }
1340  else if ( y > ymax )
1341  {
1342  *tx = *( xg + ( ny - 1 ) * nx );
1343  *ty = *( yg + ( ny - 1 ) * nx );
1344  }
1345  else
1346  {
1347  ul = 0;
1348  xll = *( xg + ul + vl * nx );
1349  yll = *( yg + ul + vl * nx );
1350  xlr = *( xg + ul + vr * nx );
1351  ylr = *( yg + ul + vr * nx );
1352 
1353  *tx = xll * ( 1 - dv ) + xlr * ( dv );
1354  *ty = yll * ( 1 - dv ) + ylr * ( dv );
1355  }
1356  }
1357  else if ( x > xmax )
1358  {
1359  if ( y < ymin )
1360  {
1361  *tx = *( xg + ( nx - 1 ) );
1362  *ty = *( yg + ( nx - 1 ) );
1363  }
1364  else if ( y > ymax )
1365  {
1366  *tx = *( xg + ( nx - 1 ) + ( ny - 1 ) * nx );
1367  *ty = *( yg + ( nx - 1 ) + ( ny - 1 ) * nx );
1368  }
1369  else
1370  {
1371  ul = nx - 1;
1372  xll = *( xg + ul + vl * nx );
1373  yll = *( yg + ul + vl * nx );
1374  xlr = *( xg + ul + vr * nx );
1375  ylr = *( yg + ul + vr * nx );
1376 
1377  *tx = xll * ( 1 - dv ) + xlr * ( dv );
1378  *ty = yll * ( 1 - dv ) + ylr * ( dv );
1379  }
1380  }
1381  else
1382  {
1383  if ( y < ymin )
1384  {
1385  vl = 0;
1386  xll = *( xg + ul + vl * nx );
1387  xrl = *( xg + ur + vl * nx );
1388  yll = *( yg + ul + vl * nx );
1389  yrl = *( yg + ur + vl * nx );
1390 
1391  *tx = xll * ( 1 - du ) + xrl * ( du );
1392  *ty = yll * ( 1 - du ) + yrl * ( du );
1393  }
1394  else if ( y > ymax )
1395  {
1396  vr = ny - 1;
1397  xlr = *( xg + ul + vr * nx );
1398  xrr = *( xg + ur + vr * nx );
1399  ylr = *( yg + ul + vr * nx );
1400  yrr = *( yg + ur + vr * nx );
1401 
1402  *tx = xlr * ( 1 - du ) + xrr * ( du );
1403  *ty = ylr * ( 1 - du ) + yrr * ( du );
1404  }
1405  }
1406  }
1407 
1408 // Normal case.
1409 // Look up coordinates in row-dominant array.
1410 // Have to handle right boundary specially -- if at the edge, we'd
1411 // better not reference the out of bounds point.
1412 
1413  else
1414  {
1415  xll = *( xg + ul + vl * nx );
1416  yll = *( yg + ul + vl * nx );
1417 
1418 // ur is out of bounds
1419 
1420  if ( ur == nx && vr < ny )
1421  {
1422  xlr = *( xg + ul + vr * nx );
1423  ylr = *( yg + ul + vr * nx );
1424 
1425  *tx = xll * ( 1 - dv ) + xlr * ( dv );
1426  *ty = yll * ( 1 - dv ) + ylr * ( dv );
1427  }
1428 
1429 // vr is out of bounds
1430 
1431  else if ( ur < nx && vr == ny )
1432  {
1433  xrl = *( xg + ur + vl * nx );
1434  yrl = *( yg + ur + vl * nx );
1435 
1436  *tx = xll * ( 1 - du ) + xrl * ( du );
1437  *ty = yll * ( 1 - du ) + yrl * ( du );
1438  }
1439 
1440 // both ur and vr are out of bounds
1441 
1442  else if ( ur == nx && vr == ny )
1443  {
1444  *tx = xll;
1445  *ty = yll;
1446  }
1447 
1448 // everything in bounds
1449 
1450  else
1451  {
1452  xrl = *( xg + ur + vl * nx );
1453  xlr = *( xg + ul + vr * nx );
1454  xrr = *( xg + ur + vr * nx );
1455 
1456  yrl = *( yg + ur + vl * nx );
1457  ylr = *( yg + ul + vr * nx );
1458  yrr = *( yg + ur + vr * nx );
1459  *tx = xll * ( 1 - du ) * ( 1 - dv ) + xlr * ( 1 - du ) * ( dv ) +
1460  xrl * ( du ) * ( 1 - dv ) + xrr * ( du ) * ( dv );
1461 
1462  *ty = yll * ( 1 - du ) * ( 1 - dv ) + ylr * ( 1 - du ) * ( dv ) +
1463  yrl * ( du ) * ( 1 - dv ) + yrr * ( du ) * ( dv );
1464  }
1465  }
1466 }