PLplot: plcont.c Source File

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