Back to index

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