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