Back to index

plt-scheme  4.2.1
plot3d.c
Go to the documentation of this file.
00001 /* $Id: plot3d.c,v 1.2 2005/03/17 21:39:21 eli Exp $
00002 
00003        3d plot routines.
00004 */
00005 
00006 #include "plplotP.h"
00007 
00008 /* Internal constants */
00009 
00010 #define  BINC      50              /* Block size for memory allocation */
00011 
00012 static PLINT pl3mode = 0;   /* 0 3d solid; 1 mesh plot */
00013 static PLINT pl3upv = 1;    /* 1 update view; 0 no update */
00014 
00015 static PLINT zbflg = 0, zbcol;
00016 static PLFLT zbtck;
00017 
00018 static PLINT *oldhiview = NULL;
00019 static PLINT *oldloview = NULL;
00020 static PLINT *newhiview = NULL;
00021 static PLINT *newloview = NULL;
00022 static PLINT *utmp = NULL;
00023 static PLINT *vtmp = NULL;
00024 static PLFLT *ctmp = NULL;
00025 
00026 static PLINT mhi, xxhi, newhisize;
00027 static PLINT mlo, xxlo, newlosize;
00028 
00029 /* Light source for shading */
00030 static PLFLT xlight, ylight, zlight;
00031 static PLINT falsecolor=0;
00032 static PLFLT fc_minz, fc_maxz;
00033 
00034 /* Prototypes for static functions */
00035 
00036 static void plgrid3  (PLFLT);
00037 static void plnxtv (PLINT *, PLINT *, PLFLT*, PLINT, PLINT);
00038 static void plside3  (PLFLT *, PLFLT *, PLFLT **, PLINT, PLINT, PLINT);
00039 static void plt3zz   (PLINT, PLINT, PLINT, PLINT, 
00040                       PLINT, PLINT *, PLFLT *, PLFLT *, PLFLT **, 
00041                         PLINT, PLINT, PLINT *, PLINT *, PLFLT*);
00042 static void plnxtvhi (PLINT *, PLINT *, PLFLT*, PLINT, PLINT);
00043 static void plnxtvlo (PLINT *, PLINT *, PLFLT*, PLINT, PLINT);
00044 static void plnxtvhi_draw(PLINT *u, PLINT *v, PLFLT* c, PLINT n);
00045 
00046 static void savehipoint     (PLINT, PLINT);
00047 static void savelopoint     (PLINT, PLINT);
00048 static void swaphiview      (void);
00049 static void swaploview      (void);
00050 static void myexit   (char *);
00051 static void myabort  (char *);
00052 static void freework (void);
00053 static int  plabv    (PLINT, PLINT, PLINT, PLINT, PLINT, PLINT);
00054 static void pl3cut   (PLINT, PLINT, PLINT, PLINT, PLINT, 
00055                             PLINT, PLINT, PLINT, PLINT *, PLINT *);
00056 static PLFLT plGetAngleToLight(PLFLT* x, PLFLT* y, PLFLT* z);
00057 static void plP_draw3d(PLINT x, PLINT y, PLFLT *c, PLINT j, PLINT move);
00058 
00059 /* #define MJL_HACK 1 */
00060 #if MJL_HACK
00061 static void plP_fill3(PLINT x0, PLINT y0, PLINT x1, PLINT y1,
00062                     PLINT x2, PLINT y2, PLINT j);
00063 static void plP_fill4(PLINT x0, PLINT y0, PLINT x1, PLINT y1,
00064                     PLINT x2, PLINT y2, PLINT x3, PLINT y3, PLINT j);
00065 #endif
00066 
00067 /*--------------------------------------------------------------------------*\
00068  * void plsetlightsource(x, y, z)
00069  *
00070  * Sets the position of the light source.
00071 \*--------------------------------------------------------------------------*/
00072 
00073 void
00074 c_pllightsource(PLFLT x, PLFLT y, PLFLT z)
00075 {
00076     xlight = x;
00077     ylight = y;
00078     zlight = z;
00079 }
00080 
00081 /*--------------------------------------------------------------------------*\
00082  * void plmesh(x, y, z, nx, ny, opt)
00083  *
00084  * Plots a mesh representation of the function z[x][y]. The x values
00085  * are stored as x[0..nx-1], the y values as y[0..ny-1], and the
00086  * z values are in the 2-d array z[][]. The integer "opt" specifies:
00087  * see plmeshc() bellow.
00088 \*--------------------------------------------------------------------------*/
00089 
00090 void
00091 c_plmesh(PLFLT *x, PLFLT *y, PLFLT **z, PLINT nx, PLINT ny, PLINT opt)
00092 {
00093   c_plot3dc(x, y, z, nx, ny, opt | MESH, NULL, 0);
00094 }
00095 
00096 /*--------------------------------------------------------------------------*\
00097  * void plmeshc(x, y, z, nx, ny, opt, clevel, nlevel)
00098  *
00099  * Plots a mesh representation of the function z[x][y]. The x values
00100  * are stored as x[0..nx-1], the y values as y[0..ny-1], and the
00101  * z values are in the 2-d array z[][]. The integer "opt" specifies:
00102  *
00103  * DRAW_LINEX   draw lines parallel to the X axis
00104  * DRAW_LINEY   draw lines parallel to the Y axis
00105  * DRAW_LINEXY  draw lines parallel to both the X and Y axis
00106  * MAG_COLOR    draw the mesh with a color dependent of the magnitude
00107  * BASE_CONT    draw contour plot at bottom xy plane
00108  * TOP_CONT     draw contour plot at top xy plane (not yet)
00109  * DRAW_SIDES   draw sides
00110  * 
00111  * or any bitwise combination, e.g. "MAG_COLOR | DRAW_LINEX"
00112  *
00113 \*--------------------------------------------------------------------------*/
00114 
00115 MZ_DLLEXPORT
00116 void
00117 c_plmeshc(PLFLT *x, PLFLT *y, PLFLT **z, PLINT nx, PLINT ny, PLINT opt,
00118         PLFLT *clevel, PLINT nlevel)
00119 {
00120     plot3dc(x, y, z, nx, ny, opt | MESH, clevel, nlevel);
00121 }
00122 
00123 /* clipping helper for 3d polygons */
00124 
00125 int
00126 plP_clip_poly(int Ni, PLFLT *Vi[3], int axis, PLFLT dir, PLFLT offset)
00127 {
00128   int anyout = 0;
00129   PLFLT in[PL_MAXPOLY], T[3][PL_MAXPOLY];
00130   int No = 0;
00131   int i, j, k;
00132 
00133   for(i=0; i<Ni; i++) {
00134     in[i] = Vi[axis][i] * dir + offset;
00135     anyout += in[i] < 0;
00136   }
00137 
00138   /* none out */
00139   if(anyout == 0)
00140     return Ni;
00141 
00142   /* all out */
00143   if(anyout == Ni) {
00144     return 0;
00145   }
00146 
00147   /* some out */
00148   /* copy over to a temp vector */
00149   for(i=0; i<3; i++) {
00150     for(j=0; j<Ni; j++) {
00151       T[i][j] = Vi[i][j];
00152     }
00153   }
00154   /* copy back selectively */
00155   for(i=0; i<Ni; i++) {
00156     j = (i+1) % Ni;
00157 
00158     if(in[i]>=0 && in[j]>=0) {
00159       for(k=0; k<3; k++)
00160        Vi[k][No] = T[k][j];
00161       No++;
00162     } else if(in[i]>=0 && in[j]<0) {
00163       PLFLT u = in[i] / (in[i] - in[j]);
00164       for(k = 0; k<3; k++)
00165        Vi[k][No] = T[k][i]*(1-u) + T[k][j]*u;
00166       No++;
00167     } else if(in[i]<0 && in[j]>=0) {
00168       PLFLT u = in[i] / (in[i] - in[j]);
00169       for(k = 0; k<3; k++)
00170        Vi[k][No] = T[k][i]*(1-u) + T[k][j]*u;
00171       No++;
00172       for(k=0; k<3; k++)
00173        Vi[k][No] = T[k][j];
00174       No++;
00175     }
00176   }
00177   return No;
00178 }
00179 
00180 /* helper for plsurf3d, similar to c_plfill3() */
00181 static void
00182 shade_triangle(PLFLT x0, PLFLT y0, PLFLT z0,
00183                         PLFLT x1, PLFLT y1, PLFLT z1,
00184                         PLFLT x2, PLFLT y2, PLFLT z2)
00185 {
00186   int i;
00187   /* arrays for interface to core functions */
00188   short u[6], v[6];
00189   PLFLT x[6], y[6], z[6], c;
00190   int n;
00191   PLFLT xmin, xmax, ymin, ymax, zmin, zmax, zscale;
00192   PLFLT *V[3];
00193 
00194   plP_gdom(&xmin, &xmax, &ymin, &ymax);
00195   plP_grange(&zscale, &zmin, &zmax);
00196 
00197   x[0] = x0; x[1] = x1; x[2] = x2;
00198   y[0] = y0; y[1] = y1; y[2] = y2;
00199   z[0] = z0; z[1] = z1; z[2] = z2;
00200   n = 3;
00201 
00202   if (falsecolor)
00203       plcol1(((z[0] + z[1] + z[2]) /3. - fc_minz) / (fc_maxz - fc_minz));
00204   else
00205     plcol1(plGetAngleToLight(x, y, z));
00206 
00207   V[0] = x; V[1] = y; V[2] = z;
00208 
00209   n = plP_clip_poly(n, V, 0,  1, -xmin);
00210   n = plP_clip_poly(n, V, 0, -1,  xmax);
00211   n = plP_clip_poly(n, V, 1,  1, -ymin);
00212   n = plP_clip_poly(n, V, 1, -1,  ymax);
00213   n = plP_clip_poly(n, V, 2,  1, -zmin);
00214   n = plP_clip_poly(n, V, 2, -1,  zmax);
00215 
00216   for(i=0; i<n; i++) {
00217     u[i] = plP_wcpcx(plP_w3wcx(x[i], y[i], z[i]));
00218     v[i] = plP_wcpcy(plP_w3wcy(x[i], y[i], z[i]));
00219     }
00220   u[n] = u[0];
00221   v[n] = v[0];
00222 
00223 #ifdef SHADE_DEBUG /* show triangles */
00224   plcol0(1);
00225   x[3]=x[0]; y[3]=y[0]; z[3]=z[0]; 
00226   plline3(4,x,y,z);
00227 #else /* fill triangles */
00228   plP_fill(u, v, n+1);
00229 #endif
00230 }
00231 
00232 /*--------------------------------------------------------------------------*\
00233  * void plotsurf3d(x, y, z, nx, ny, opt, clevel, nlevel)
00234  *
00235  * Plots the 3-d surface representation of the function z[x][y].
00236  * The x values are stored as x[0..nx-1], the y values as y[0..ny-1],
00237  *  and the z values are in the 2-d array z[][].
00238  *
00239  *
00240  * DRAW_LINEX   draw lines parallel to the X axis
00241  * DRAW_LINEY   draw lines parallel to the Y axis
00242  * DRAW_LINEXY  draw lines parallel to both the X and Y axis
00243  * BASE_CONT    draw contour plot at bottom xy plane
00244  * TOP_CONT     draw contour plot at top xy plane
00245  * SURF_CONT    draw contour at surface
00246  * FACETED      each square that makes up the surface is faceted
00247  * DRAW_SIDES   draw sides
00248  * MAG_COLOR    the surface is colored according to the value of z;
00249  *               if not defined, the surface is colored according to the
00250  *               intensity of the reflected light in the surface from a light
00251  *               source whose position is set using c_pllightsource()
00252  *
00253  * or any bitwise combination, e.g. "MAG_COLOR | SURF_CONT | SURF_BASE"
00254  *
00255  * This code is a complete departure from the approach taken in the old version
00256  * of this routine. Formerly to code attempted to use the logic for the hidden
00257  * line algorithm to draw the hidden surface. This was really hard. This code
00258  * below uses a simple back to front (painters) algorithm. All the
00259  * triangles are drawn.
00260  *
00261  * There are multitude of ways this code could be optimized. Given the
00262  * problems with the old code, I tried to focus on clarity here. 
00263 \*--------------------------------------------------------------------------*/
00264 
00265 void
00266 plsurf3d(PLFLT *x, PLFLT *y, PLFLT **z, PLINT nx, PLINT ny,
00267          PLINT opt, PLFLT *clevel, PLINT nlevel)
00268 {
00269   PLFLT cxx, cxy, cyx, cyy, cyz;
00270   PLINT i, j, k;
00271   PLINT ixDir, ixOrigin, iyDir, iyOrigin, nFast, nSlow;
00272   PLINT ixFast, ixSlow, iyFast, iySlow;
00273   PLINT iFast, iSlow;
00274   PLFLT xmin, xmax, ymin, ymax, zmin, zmax, zscale;
00275   PLFLT xm, ym, zm;
00276   PLINT ixmin=0, ixmax=nx, iymin=0, iymax=ny;
00277   PLFLT xx[3],yy[3],zz[3];
00278   PLFLT px[4], py[4], pz[4];
00279   CONT_LEVEL *cont, *clev;
00280   CONT_LINE *cline;
00281   int   ct, ix, iy;
00282 
00283   if (plsc->level < 3) {
00284     myabort("plsurf3d: Please set up window first");
00285     return;
00286   }
00287     
00288   if (nx <= 0 || ny <= 0) {
00289     myabort("plsurf3d: Bad array dimensions.");
00290     return;
00291   }
00292 
00293   /*  
00294    * Don't use the data z value to scale the color, use the z axis
00295    * values set by plw3d()
00296    *
00297    * plMinMax2dGrid(z, nx, ny, &fc_maxz, &fc_minz);
00298    */
00299 
00300   fc_minz = plsc->ranmi;
00301   fc_maxz = plsc->ranma;
00302   if (fc_maxz == fc_minz) {
00303     plwarn("plsurf3d.c: Maximum and minimum Z values are equal! \"fixing\"...");
00304     fc_maxz = fc_minz + 1e-6;
00305   }
00306 
00307   if (opt & MAG_COLOR)
00308     falsecolor = 1;
00309   else
00310     falsecolor = 0;
00311     
00312   plP_gdom(&xmin, &xmax, &ymin, &ymax);
00313   plP_grange(&zscale, &zmin, &zmax);
00314   if(zmin > zmax) {
00315     PLFLT t = zmin;
00316     zmin = zmax;
00317     zmax = t;
00318   }
00319 
00320   /* Check that points in x and in y are strictly increasing  and in range */
00321 
00322   for (i = 0; i < nx - 1; i++) {
00323     if (x[i] >= x[i + 1]) {
00324       myabort("plsurf3d: X array must be strictly increasing");
00325       return;
00326     }
00327     if (x[i] < xmin && x[i+1] >= xmin)
00328       ixmin = i;
00329     if (x[i+1] > xmax && x[i] <= xmax)
00330       ixmax = i+2;
00331   }
00332   for (i = 0; i < ny - 1; i++) {
00333     if (y[i] >= y[i + 1]) {
00334       myabort("plsurf3d: Y array must be strictly increasing");
00335       return;
00336     }
00337     if (y[i] < ymin && y[i+1] >= ymin)
00338       iymin = i;
00339     if (y[i+1] > ymax && y[i] <= ymax)
00340       iymax = i+2;
00341   }
00342 
00343   /* get the viewing parameters */
00344   plP_gw3wc(&cxx, &cxy, &cyx, &cyy, &cyz);
00345 
00346   /* we're going to draw from back to front */
00347 
00348   /* iFast will index the dominant (fastest changing) dimension
00349      iSlow will index the slower changing dimension
00350 
00351      iX indexes the X dimension
00352      iY indexes the Y dimension */
00353 
00354   /* get direction for X */
00355   if(cxy >= 0) {
00356     ixDir = 1;       /* direction in X */
00357     ixOrigin = ixmin;       /* starting point */
00358   } else {
00359     ixDir = -1;
00360     ixOrigin = ixmax - 1;
00361   }
00362   /* get direction for Y */
00363   if(cxx >= 0) {
00364     iyDir = -1;
00365     iyOrigin = iymax - 1;
00366   } else {
00367     iyDir = 1;
00368     iyOrigin = iymin;
00369   }
00370   /* figure out which dimension is dominant */
00371   if (fabs(cxx) > fabs(cxy)) {
00372     /* X is dominant */
00373     nFast = ixmax - ixmin;  /* samples in the Fast direction */
00374     nSlow = iymax - iymin;  /* samples in the Slow direction */
00375       
00376     ixFast = ixDir; ixSlow = 0;
00377     iyFast = 0;     iySlow = iyDir;
00378   } else {
00379     nFast = iymax - iymin;
00380     nSlow = ixmax - ixmin;
00381 
00382     ixFast = 0;     ixSlow = ixDir;
00383     iyFast = iyDir; iySlow = 0;
00384   }
00385 
00386   /* we've got to draw the background grid first, hidden line code has to draw it last */
00387   if (zbflg) {
00388     PLINT color = plsc->icol0;
00389     PLFLT bx[3], by[3], bz[3];
00390     PLFLT tick=0, tp;
00391     PLINT nsub=0;
00392 
00393     /* get the tick spacing */
00394     pldtik(zmin, zmax, &tick, &nsub);
00395 
00396     /* determine the vertices for the background grid line */
00397     bx[0] = (ixOrigin!=ixmin && ixFast==0) || ixFast > 0 ? xmax : xmin;
00398     by[0] = (iyOrigin!=iymin && iyFast==0) || iyFast > 0 ? ymax : ymin;
00399     bx[1] = ixOrigin!=ixmin ? xmax : xmin;
00400     by[1] = iyOrigin!=iymin ? ymax : ymin;
00401     bx[2] = (ixOrigin!=ixmin && ixSlow==0) || ixSlow > 0 ? xmax : xmin;
00402     by[2] = (iyOrigin!=iymin && iySlow==0) || iySlow > 0 ? ymax : ymin;
00403 
00404     plcol(zbcol);
00405     for(tp = tick * floor(zmin / tick) + tick; tp <= zmax; tp += tick) {
00406       bz[0] = bz[1] = bz[2] = tp;
00407       plline3(3, bx, by, bz);      
00408     }
00409     /* draw the vertical line at the back corner */
00410     bx[0] = bx[1];
00411     by[0] = by[1];
00412     bz[0] = zmin;
00413     plline3(2, bx, by, bz);
00414     plcol(color);
00415   }
00416 
00417   /* If enabled, draw the contour at the base */
00418 
00419   /* The contour ploted at the base will be identical to the one obtained
00420    * with c_plcont(). The contour ploted at the surface is simple minded, but
00421    * can be improved by using the contour data available.
00422    */
00423   
00424   if (clevel != NULL && opt & BASE_CONT) {
00425 #define NPTS 100
00426     int np = NPTS;
00427     PLFLT *zz = (PLFLT *) malloc(NPTS*sizeof(PLFLT));
00428 
00429     /* get the contour lines */
00430     cont_store(x, y, z,  nx, ny, 1, nx, 1, ny, clevel, nlevel, &cont);
00431 
00432     /* follow the contour levels and lines */
00433     clev = cont;
00434     do { /* for each contour level */
00435       cline = clev->line;
00436       do { /* there are several lines that make up the contour */
00437        if (cline->npts > np) {
00438          np = cline->npts;
00439          zz = (PLFLT *) realloc(zz, np*sizeof(PLFLT));
00440        }
00441        for (j=0; j<cline->npts; j++)
00442          zz[j] = plsc->ranmi;
00443        plcol1((clev->level-fc_minz)/(fc_maxz-fc_minz));
00444        plline3(cline->npts, cline->x, cline->y, zz);
00445        cline = cline->next;
00446       }
00447       while(cline != NULL);
00448       clev = clev->next;
00449     }
00450     while(clev != NULL);
00451   
00452     cont_clean_store(cont); /* now release the memory */
00453     free(zz);
00454   }
00455 
00456   /* Now we can iterate over the grid drawing the quads */
00457   for(iSlow=0; iSlow < nSlow-1; iSlow++) {
00458     for(iFast=0; iFast < nFast-1; iFast++) {
00459       /* get the 4 corners of the Quad, which are
00460        *      
00461        *       0--2
00462        *       |  |
00463        *       1--3
00464        */
00465 
00466       xm = ym = zm = 0.;
00467 
00468       for(i=0; i<2; i++) {
00469        for(j=0; j<2; j++) {
00470          /* we're transforming from Fast/Slow coordinates to x/y coordinates */
00471          /* note, these are the indices, not the values */
00472          ix = ixFast * (iFast+i) + ixSlow * (iSlow+j) + ixOrigin;
00473          iy = iyFast * (iFast+i) + iySlow * (iSlow+j) + iyOrigin;
00474 
00475          xm += px[2*i+j] = x[ix];
00476          ym += py[2*i+j] = y[iy];
00477          zm += pz[2*i+j] = z[ix][iy];
00478        }
00479       }
00480 
00481       /* the "mean point" of the quad, common to all four triangles
00482         -- perhaps not a good thing to do for the light shading */
00483 
00484       xm /= 4.; ym /= 4.; zm /= 4.;
00485 
00486       /* now draw the quad as four triangles */
00487 
00488       for (i=1; i<3; i++) {
00489        for (j=0; j<4; j+=3) {
00490          shade_triangle(px[j], py[j], pz[j], xm, ym, zm, px[i], py[i], pz[i]);
00491 
00492          /* after shading, see if the triangle crosses  one contour plane */
00493 
00494 #define min3(a,b,c) (MIN((MIN(a,b)),c))
00495 #define max3(a,b,c) (MAX((MAX(a,b)),c))
00496 
00497          if (clevel != NULL && (opt & SURF_CONT)) {
00498            for (k=0; k<nlevel; k++) {
00499              if (clevel[k] >= min3(pz[i],zm,pz[j]) && clevel[k] < max3(pz[i],zm,pz[j])) {
00500               ct = 0;
00501               if (clevel[k] >= MIN(pz[i],zm) && clevel[k] < MAX(pz[i],zm)) { /* p0-pm */
00502                 xx[ct] = ((clevel[k] - pz[i])*(xm-px[i]))/(zm-pz[i]) + px[i];
00503                 yy[ct] = ((clevel[k] - pz[i])*(ym-py[i]))/(zm-pz[i]) + py[i];
00504                 ct++;
00505               }
00506 
00507               if (clevel[k] >= MIN(pz[i],pz[j]) && clevel[k] < MAX(pz[i],pz[j])) { /* p0-p1 */
00508                 xx[ct] = ((clevel[k] - pz[i])*(px[j]-px[i]))/(pz[j]-pz[i]) + px[i];
00509                 yy[ct] = ((clevel[k] - pz[i])*(py[j]-py[i]))/(pz[j]-pz[i]) + py[i];
00510                 ct++;
00511               }
00512 
00513               if (clevel[k] >= MIN(pz[j],zm) && clevel[k] < MAX(pz[j],zm)) { /* p1-pm */
00514                 xx[ct] = ((clevel[k] - pz[j])*(xm-px[j]))/(zm-pz[j]) + px[j];
00515                 yy[ct] = ((clevel[k] - pz[j])*(ym-py[j]))/(zm-pz[j]) + py[j];
00516                 ct++;
00517               }
00518               
00519               if (ct == 2) {
00520 
00521                 /* yes, xx and yy are the intersection points of the triangle with
00522                  * the contour line -- draw a straight line betweeen the points
00523                  * -- at the end this will make up the contour line */
00524 
00525                 if (opt & SURF_CONT) {
00526                   /* surface contour with black color (suggestions?) */
00527                   plcol0(0);
00528                   zz[0] = zz[1] = clevel[k]; 
00529                   plline3(2, xx, yy, zz);
00530                 }
00531 
00532                 /* don't break; one triangle can span various contour levels */
00533 
00534               } else
00535                 plwarn("plot3d.c:plsurf3d() ***ERROR***\n");
00536              } 
00537            }
00538          }
00539        }
00540       }
00541     }
00542   }
00543 
00544   if (opt & FACETED) {
00545     plcol0(0);
00546     plot3dc(x, y, z, nx, ny, MESH | DRAW_LINEXY, NULL, 0);
00547   }
00548 
00549   if (opt & DRAW_SIDES) { /* the sides look ugly !!! */
00550     /* draw one more row with all the Z's set to zmin */
00551     PLFLT zscale, zmin, zmax;
00552 
00553     plP_grange(&zscale, &zmin, &zmax);
00554     
00555     iSlow = nSlow-1;
00556     for(iFast=0; iFast < nFast-1; iFast++) {
00557       for(i=0; i<2; i++) {
00558        ix = ixFast * (iFast+i) + ixSlow * iSlow + ixOrigin;
00559        iy = iyFast * (iFast+i) + iySlow * iSlow + iyOrigin;
00560        px[2*i] = x[ix]; 
00561        py[2*i] = y[iy];
00562        pz[2*i] = z[ix][iy];
00563       }
00564       /* now draw the quad as two triangles (4 might be better) */
00565         
00566       shade_triangle(px[0], py[0], pz[0], px[2], py[2], pz[2], px[0], py[0], zmin);
00567       shade_triangle(px[2], py[2], pz[2], px[2], py[2], zmin,  px[0], py[0], zmin);
00568     }
00569     
00570     iFast = nFast-1;
00571     for(iSlow=0; iSlow < nSlow-1; iSlow++) {
00572       for(i=0; i<2; i++) {
00573          ix = ixFast * iFast + ixSlow * (iSlow+i) + ixOrigin;
00574          iy = iyFast * iFast + iySlow * (iSlow+i) + iyOrigin;
00575          px[2*i] = x[ix];
00576          py[2*i] = y[iy];
00577          pz[2*i] = z[ix][iy];
00578        }
00579 
00580       /* now draw the quad as two triangles (4 might be better) */
00581       shade_triangle(px[0], py[0], pz[0], px[2], py[2], pz[2], px[0], py[0], zmin);
00582       shade_triangle(px[2], py[2], pz[2], px[2], py[2], zmin,  px[0], py[0], zmin);
00583     }
00584   }
00585 }
00586 
00587 /*--------------------------------------------------------------------------*\
00588  * void plot3d(x, y, z, nx, ny, opt, side)
00589  *
00590  * Plots a 3-d representation of the function z[x][y]. The x values
00591  * are stored as x[0..nx-1], the y values as y[0..ny-1], and the z
00592  * values are in the 2-d array z[][]. The integer "opt" specifies:
00593  * see plot3dc() bellow
00594 \*--------------------------------------------------------------------------*/
00595 
00596 MZ_DLLEXPORT
00597 void
00598 c_plot3d(PLFLT *x, PLFLT *y, PLFLT **z,
00599         PLINT nx, PLINT ny, PLINT opt, PLINT side)
00600 {
00601   c_plot3dc( x, y, z, nx, ny, opt | (side == 1 ? DRAW_SIDES : 0), NULL, 0);
00602 }
00603 
00604 /*--------------------------------------------------------------------------*\
00605  * void plot3dc(x, y, z, nx, ny, opt, clevel, nlevel)
00606  *
00607  * Plots a 3-d representation of the function z[x][y]. The x values
00608  * are stored as x[0..nx-1], the y values as y[0..ny-1], and the z
00609  * values are in the 2-d array z[][]. The integer "opt" specifies:
00610  *
00611  *  DRAW_LINEX :  Draw lines parallel to x-axis
00612  *  DRAW_LINEY :  Draw lines parallel to y-axis
00613  *  DRAW_LINEXY:  Draw lines parallel to both axes
00614  *  MAG_COLOR:    Magnitude coloring of wire frame
00615  *  BASE_CONT:    Draw contour at bottom xy plane 
00616  *  TOP_CONT:     Draw contour at top xy plane (not yet)
00617  *  DRAW_SIDES:   Draw sides around the plot
00618  *  MESH:       Draw the "under" side of the plot
00619  *
00620  * or any bitwise combination, e.g. "MAG_COLOR | DRAW_LINEX"
00621  *
00622 \*--------------------------------------------------------------------------*/
00623 
00624 void
00625 c_plot3dc(PLFLT *x, PLFLT *y, PLFLT **z,
00626         PLINT nx, PLINT ny, PLINT opt,
00627         PLFLT *clevel, PLINT nlevel)
00628 {
00629     PLFLT cxx, cxy, cyx, cyy, cyz;
00630     PLINT init, i, ix, iy, color;
00631     PLFLT xmin, xmax, ymin, ymax, zmin, zmax, zscale;
00632     PLINT ixmin=0, ixmax=nx-1, iymin=0, iymax=ny-1;
00633     PLINT clipped = 0, base_cont = 0, side = 0;
00634 
00635     pl3mode = 0;
00636 
00637     if (plsc->level < 3) {
00638        myabort("plot3dc: Please set up window first");
00639        return;
00640     }
00641 
00642     if (opt < 1) {
00643        myabort("plot3dc: Bad option");
00644        return;
00645     }
00646 
00647     if (nx <= 0 || ny <= 0) {
00648        myabort("plot3dc: Bad array dimensions.");
00649        return;
00650     }
00651 
00652     plP_gdom(&xmin, &xmax, &ymin, &ymax);
00653     plP_grange(&zscale, &zmin, &zmax);
00654     if(zmin > zmax) {
00655       PLFLT t = zmin;
00656       zmin = zmax;
00657       zmax = t;
00658     }
00659 
00660 /* Check that points in x and in y are strictly increasing */
00661 
00662     for (i = 0; i < nx - 1; i++) {
00663        if (x[i] >= x[i + 1]) {
00664            myabort("plot3dc: X array must be strictly increasing");
00665            return;
00666        }
00667     }
00668     for (i = 0; i < ny - 1; i++) {
00669        if (y[i] >= y[i + 1]) {
00670            myabort("plot3dc: Y array must be strictly increasing");
00671            return;
00672        }
00673     }
00674 
00675     if (opt & MESH)
00676       pl3mode = 1;
00677 
00678     if (opt & DRAW_SIDES)
00679       side = 1;
00680 
00681     /* figure out the part of the data to use */
00682     if (xmin < x[0])
00683       xmin = x[0];
00684     if (xmax > x[nx-1])
00685       xmax = x[nx-1];
00686     if (ymin < y[0])
00687       ymin = y[0];
00688     if (ymax > y[ny-1])
00689       ymax = y[ny-1];
00690     for (ixmin = 0; ixmin < nx-1 && x[ixmin+1] < xmin; ixmin++) {}
00691     for (ixmax = nx-1; ixmax > 0 && x[ixmax-1] > xmax; ixmax--) {}
00692     for (iymin = 0; iymin < ny-1 && y[iymin+1] < ymin; iymin++) {}
00693     for (iymax = ny-1; iymax > 0 && y[iymax-1] > ymax; iymax--) {}
00694     /*fprintf(stderr, "(%d,%d) %d %d %d %d\n", nx, ny, ixmin, ixmax, iymin, iymax);*/
00695     /* do we need to clip? */
00696     if(ixmin > 0 || ixmax < nx-1 || iymin > 0 || iymax < ny-1) {
00697       /* adjust the input so it stays within bounds */
00698       int _nx = ixmax - ixmin + 1;
00699       int _ny = iymax - iymin + 1;
00700       PLFLT *_x, *_y, **_z;
00701       PLFLT ty0, ty1, tx0, tx1;
00702       int i, j;
00703 
00704       if(_nx <= 1 || _ny <= 1) {
00705        fprintf(stderr, "bail\n");
00706        return;
00707       }
00708 
00709       /* allocate storage for new versions of the input vectors */
00710       _x = (PLFLT*)malloc(_nx * sizeof(PLFLT));
00711       _y = (PLFLT*)malloc(_ny * sizeof(PLFLT));
00712       _z = (PLFLT**)malloc(_nx * sizeof(PLFLT*));
00713 
00714       clipped = 1;
00715 
00716       /* copy over the independent variables */
00717       _x[0] = xmin;
00718       _x[_nx-1] = xmax;
00719       for(i=1; i<_nx-1; i++)
00720        _x[i] = x[ixmin+i];
00721       _y[0] = ymin;
00722       _y[_ny-1] = ymax;
00723       for(i=1; i<_ny-1; i++)
00724        _y[i] = y[iymin+i];
00725 
00726       /* copy the data array so we can interpolate around the edges */
00727       for(i=0; i<_nx; i++)
00728        _z[i] = (PLFLT*)malloc(_ny * sizeof(PLFLT));
00729 
00730       /* interpolation factors for the 4 edges */
00731       ty0 = (_y[0] - y[iymin]) / (y[iymin+1] - y[iymin]);
00732       ty1 = (_y[_ny-1] - y[iymax-1]) / (y[iymax] - y[iymax-1]);
00733       tx0 = (_x[0] - x[ixmin]) / (x[ixmin+1] - x[ixmin]);
00734       tx1 = (_x[_nx-1] - x[ixmax-1]) / (x[ixmax] - x[ixmax-1]);
00735       for(i=0; i<_nx; i++) {
00736        if(i==0) {
00737          _z[i][0] = z[ixmin][iymin]*(1-ty0)*(1-tx0) + z[ixmin][iymin+1]*(1-tx0)*ty0
00738            + z[ixmin+1][iymin]*tx0*(1-ty0) + z[ixmin+1][iymin+1]*tx0*ty0;
00739          for(j=1; j<_ny-1; j++)
00740            _z[i][j] = z[ixmin][iymin+j]*(1-tx0) + z[ixmin+1][iymin+j]*tx0;
00741          _z[i][_ny-1] = z[ixmin][iymax-1]*(1-tx0)*(1-ty1) + z[ixmin][iymax]*(1-tx0)*ty1
00742            + z[ixmin+1][iymax-1]*tx0*(1-ty1) + z[ixmin+1][iymax]*tx0*ty1;
00743        } else if(i==_nx-1) {
00744          _z[i][0] = z[ixmax-1][iymin]*(1-tx1)*(1-ty0) + z[ixmax-1][iymin+1]*(1-tx1)*ty0
00745            + z[ixmax][iymin]*tx1*(1-ty0) + z[ixmax][iymin+1]*tx1*ty0;
00746          for(j=1; j<_ny-1; j++)
00747            _z[i][j] = z[ixmax-1][iymin+j]*(1-tx1) + z[ixmax][iymin+j]*tx1;
00748          _z[i][_ny-1] = z[ixmax-1][iymax-1]*(1-tx1)*(1-ty1) + z[ixmax][iymax]*(1-tx1)*ty1
00749            + z[ixmax][iymax-1]*tx1*(1-ty1) + z[ixmax][iymax]*tx1*ty1;
00750        } else {
00751          _z[i][0] = z[ixmin+i][iymin]*(1-ty0) + z[ixmin+i][iymin+1]*ty0;
00752          for(j=1; j<_ny-1; j++)
00753            _z[i][j] = z[ixmin+i][iymin+j];
00754          _z[i][_ny-1] = z[ixmin+i][iymax-1]*(1-ty1) + z[ixmin+i][iymax]*ty1;
00755        }
00756        for(j=0; j<_ny; j++) {
00757          if(_z[i][j] < zmin)
00758            _z[i][j] = zmin;
00759          else if(_z[i][j] > zmax)
00760            _z[i][j] = zmax;
00761        }
00762       }
00763       /* replace the input with our clipped versions */
00764       x = _x;
00765       y = _y;
00766       z = _z;
00767       nx = _nx;
00768       ny = _ny;        
00769     }
00770 
00771    if (opt & BASE_CONT || opt & TOP_CONT || opt && MAG_COLOR ) { 
00772      /*  
00773       * Don't use the data z value to scale the color, use the z axis
00774       * values set by plw3d()
00775       *
00776       * plMinMax2dGrid(z, nx, ny, &fc_maxz, &fc_minz);
00777       */
00778      
00779      fc_minz = plsc->ranmi;
00780      fc_maxz = plsc->ranma;     
00781 
00782      if (fc_maxz == fc_minz) {
00783        plwarn("plot3dc: Maximum and minimum Z values are equal! \"fixing\"...");
00784        fc_maxz = fc_minz + 1e-6;
00785      }
00786    }
00787  
00788    if (opt & BASE_CONT) {     /* If enabled, draw the contour at the base.  */
00789      if (clevel != NULL && nlevel != 0) {
00790        base_cont = 1;
00791        /* even if MESH is not set, "set it",
00792          as the base contour can only be done in this case */
00793        pl3mode = 1; 
00794      }
00795    }
00796 
00797    if (opt & MAG_COLOR) {     /* If enabled, use magnitude colored wireframe  */
00798      ctmp = (PLFLT *) malloc((size_t) (2 * MAX(nx, ny) * sizeof(PLFLT)));
00799    } else
00800      ctmp = NULL;
00801 
00802    /* next logic only knows opt = 1 | 2 | 3, make sure that it only gets that */
00803    opt &= DRAW_LINEXY;
00804 
00805     /* Allocate work arrays */
00806 
00807     utmp = (PLINT *) malloc((size_t) (2 * MAX(nx, ny) * sizeof(PLINT)));
00808     vtmp = (PLINT *) malloc((size_t) (2 * MAX(nx, ny) * sizeof(PLINT)));
00809 
00810     if ( ! utmp || ! vtmp)
00811        myexit("plot3dc: Out of memory.");
00812 
00813     plP_gw3wc(&cxx, &cxy, &cyx, &cyy, &cyz);
00814     init = 1;
00815 /* Call 3d line plotter.  Each viewing quadrant 
00816    (perpendicular to x-y plane) must be handled separately. */ 
00817 
00818     if (cxx >= 0.0 && cxy <= 0.0) {
00819        if (opt == DRAW_LINEY) 
00820            plt3zz(1, ny, 1, -1, -opt, &init, x, y, z, nx, ny, utmp, vtmp,ctmp);
00821        else {
00822            for (iy = 2; iy <= ny; iy++)
00823               plt3zz(1, iy, 1, -1, -opt, &init, x, y, z, nx, ny, utmp, vtmp,ctmp);
00824        }
00825        if (opt == DRAW_LINEX)
00826            plt3zz(1, ny, 1, -1, opt, &init, x, y, z, nx, ny, utmp, vtmp, ctmp);
00827        else {
00828            for (ix = 1; ix <= nx - 1; ix++)
00829               plt3zz(ix, ny, 1, -1, opt, &init, x, y, z, nx, ny, utmp, vtmp, ctmp);
00830        }
00831     }
00832 
00833     else if (cxx <= 0.0 && cxy <= 0.0) {
00834         if (opt == DRAW_LINEX)
00835            plt3zz(nx, ny, -1, -1, opt, &init, x, y, z, nx, ny, utmp, vtmp, ctmp);
00836        else {
00837            for (ix = 2; ix <= nx; ix++) 
00838               plt3zz(ix, ny, -1, -1, opt, &init, x, y, z, nx, ny, utmp, vtmp, ctmp);
00839        }
00840        if (opt == DRAW_LINEY)
00841            plt3zz(nx, ny, -1, -1, -opt, &init, x, y, z, nx, ny, utmp, vtmp, ctmp);
00842        else {
00843            for (iy = ny; iy >= 2; iy--)
00844              plt3zz(nx, iy, -1, -1, -opt, &init, x, y, z, nx, ny, utmp, vtmp, ctmp);
00845        }
00846     }
00847 
00848     else if (cxx <= 0.0 && cxy >= 0.0) {
00849        if (opt == DRAW_LINEY) 
00850            plt3zz(nx, 1, -1, 1, -opt, &init, x, y, z, nx, ny, utmp, vtmp, ctmp);
00851        else {
00852            for (iy = ny - 1; iy >= 1; iy--) 
00853               plt3zz(nx, iy, -1, 1, -opt, &init, x, y, z, nx, ny, utmp, vtmp, ctmp);
00854        }
00855        if (opt == DRAW_LINEX)
00856            plt3zz(nx, 1, -1, 1, opt, &init, x, y, z, nx, ny, utmp, vtmp, ctmp);
00857        else {
00858            for (ix = nx; ix >= 2; ix--)
00859               plt3zz(ix, 1, -1, 1, opt, &init, x, y, z, nx, ny, utmp, vtmp, ctmp);
00860        }
00861     }
00862 
00863     else if (cxx >= 0.0 && cxy >= 0.0) {
00864        if (opt == DRAW_LINEX) 
00865            plt3zz(1, 1, 1, 1, opt, &init, x, y, z, nx, ny, utmp, vtmp, ctmp);
00866        else {
00867            for (ix = nx - 1; ix >= 1; ix--) 
00868               plt3zz(ix, 1, 1, 1, opt, &init, x, y, z, nx, ny, utmp, vtmp, ctmp);
00869        }
00870        if (opt == DRAW_LINEY)
00871            plt3zz(1, 1, 1, 1, -opt, &init, x, y, z, nx, ny, utmp, vtmp, ctmp);
00872        else {
00873            for (iy = 1; iy <= ny - 1; iy++)
00874               plt3zz(1, iy, 1, 1, -opt, &init, x, y, z, nx, ny, utmp, vtmp, ctmp);
00875        }
00876     }
00877 
00878     /* draw contour at the base. Not 100%! Why? */
00879  
00880     if (base_cont){
00881       int np = NPTS, j;
00882       CONT_LEVEL *cont, *clev;
00883       CONT_LINE *cline;
00884 
00885       PLINT *uu = (PLINT *) malloc(NPTS*sizeof(PLINT));
00886       PLINT *vv = (PLINT *) malloc(NPTS*sizeof(PLINT));
00887 
00888       pl3upv = 0;
00889 
00890       /* get the contour lines */
00891       cont_store(x, y, z,  nx, ny, 1, nx, 1, ny, clevel, nlevel, &cont);
00892 
00893       /* follow the contour levels and lines */
00894       clev = cont;
00895       do { /* for each contour level */
00896        cline = clev->line;
00897        do { /* there are several lines that make up each contour */
00898          int cx, i, k, l, m, start, end;
00899          PLFLT tx, ty;
00900          if (cline->npts > np) {
00901            np = cline->npts;
00902            uu = (PLINT *) realloc(uu, np*sizeof(PLINT));
00903            vv = (PLINT *) realloc(vv, np*sizeof(PLINT));
00904          }
00905 
00906          /* the hidden line plotter plnxtv() only works OK if the x points are in increasing order.
00907             As this does not always happens, the situation must be detected and the line segment
00908             must be reversed before being plotted */
00909 
00910          plcol1((clev->level-fc_minz)/(fc_maxz-fc_minz));
00911          i = 0;
00912          do { 
00913            cx =  plP_wcpcx(plP_w3wcx(cline->x[i],cline->y[i], plsc->ranmi));
00914            for (j=i; j < cline->npts; j++) {  /* convert to 2D coordinates */
00915              uu[j] = plP_wcpcx(plP_w3wcx(cline->x[j],cline->y[j], plsc->ranmi)); 
00916              vv[j] = plP_wcpcy(plP_w3wcy(cline->x[j],cline->y[j], plsc->ranmi)); 
00917              if (uu[j] < cx) /* find turn back point */
00918               break;
00919              else
00920               cx = uu[j];
00921            }
00922            plnxtv(&uu[i], &vv[i], NULL, j-i, 0); /* plot line with increasing x */
00923 
00924            if (j < cline->npts) { /* line not yet finished, */
00925              start = j-1;
00926              for (i = start; i < cline->npts; i++) {  /* search turn forward point */
00927               uu[i] = plP_wcpcx(plP_w3wcx(cline->x[i],cline->y[i], plsc->ranmi));
00928               if (uu[i] > cx)
00929                 break;
00930               else
00931                 cx = uu[i];
00932              }
00933              end = i-1;
00934 
00935              for (k=0; k < (end-start+1)/2; k++) { /* reverse line segment */
00936               l = start+k;
00937               m = end-k;
00938               tx = cline->x[l];
00939               ty = cline->y[l];
00940               cline->x[l] = cline->x[m];
00941               cline->y[l] = cline->y[m];
00942               cline->x[m] = tx;
00943               cline->y[m] = ty;
00944              }
00945 
00946              /* convert to 2D coordinates */
00947              for (j=start; j <= end; j++) {
00948               uu[j] = plP_wcpcx(plP_w3wcx(cline->x[j],cline->y[j], plsc->ranmi)); 
00949               vv[j] = plP_wcpcy(plP_w3wcy(cline->x[j],cline->y[j], plsc->ranmi)); 
00950              }
00951              plnxtv(&uu[start], &vv[start], NULL, end-start+1, 0); /* and plot it */
00952 
00953              i = end+1; /* restart where it was left */
00954            }
00955          } while (j < cline->npts && i < cline->npts);
00956          cline = cline->next;
00957        }
00958        while(cline != NULL);
00959        clev = clev->next;
00960       }
00961       while(clev != NULL);
00962 
00963       cont_clean_store(cont); /* now release the contour memory */
00964       pl3upv = 1;
00965       free(uu);
00966       free(vv);
00967     }
00968 
00969 /* Finish up by drawing sides, background grid (both are optional) */
00970 
00971     if (side)
00972        plside3(x, y, z, nx, ny, opt);
00973 
00974     if (zbflg) {
00975        color = plsc->icol0;
00976        plcol(zbcol);
00977        plgrid3(zbtck);
00978        plcol(color);
00979     }
00980 
00981     freework();
00982 
00983     if(clipped) {
00984       free(x);
00985       free(y);
00986       for(i=0; i<nx; i++)
00987        free(z[i]);
00988       free(z);
00989     }
00990 }
00991 
00992 /*--------------------------------------------------------------------------*\
00993  * void plP_gzback()
00994  *
00995  * Get background parameters for 3d plot.
00996 \*--------------------------------------------------------------------------*/
00997 
00998 void
00999 plP_gzback(PLINT **zbf, PLINT **zbc, PLFLT **zbt)
01000 {
01001     *zbf = &zbflg;
01002     *zbc = &zbcol;
01003     *zbt = &zbtck;
01004 }
01005 
01006 /*--------------------------------------------------------------------------*\
01007  * PLFLT plGetAngleToLight()
01008  *
01009  * Gets cos of angle between normal to a polygon and a light source.
01010  * Requires at least 3 elements, forming non-parallel lines 
01011  * in the arrays.
01012 \*--------------------------------------------------------------------------*/
01013 
01014 static PLFLT
01015 plGetAngleToLight(PLFLT* x, PLFLT* y, PLFLT* z)
01016 {
01017     PLFLT vx1, vx2, vy1, vy2, vz1, vz2;
01018     PLFLT px, py, pz;
01019     PLFLT vlx, vly, vlz;
01020     PLFLT mag1, mag2;
01021     PLFLT cosangle;
01022 
01023     vx1 = x[1] - x[0];
01024     vx2 = x[2] - x[1];
01025     vy1 = y[1] - y[0];
01026     vy2 = y[2] - y[1];
01027     vz1 = z[1] - z[0];
01028     vz2 = z[2] - z[1];
01029 
01030 /* Find vector perpendicular to the face */
01031     px = vy1*vz2 - vz1*vy2;
01032     py = vz1*vx2 - vx1*vz2;
01033     pz = vx1*vy2 - vy1*vx2;
01034     mag1 = px*px + py*py + pz*pz;
01035 
01036 /* Vectors were parallel! */
01037     if (mag1 == 0) 
01038        return 1;
01039 
01040     vlx = xlight - x[0];
01041     vly = ylight - y[0];
01042     vlz = zlight - z[0];
01043     mag2 = vlx*vlx + vly*vly + vlz*vlz;
01044     if (mag2 ==0) 
01045        return 1;
01046 
01047 /* Now have 3 vectors going through the first point on the given surface */
01048     cosangle = fabs( (vlx*px + vly*py + vlz*pz) / sqrt(mag1*mag2) );
01049 
01050 /* In case of numerical rounding */
01051     if (cosangle > 1) cosangle = 1;
01052     return cosangle;
01053 }
01054 
01055 /*--------------------------------------------------------------------------*\
01056  * void plt3zz()
01057  *
01058  * Draws the next zig-zag line for a 3-d plot.  The data is stored in array
01059  * z[][] as a function of x[] and y[], and is plotted out starting at index
01060  * (x0,y0).
01061  *
01062  * Depending on the state of "flag", the sequence of data points sent to
01063  * plnxtv is altered so as to allow cross-hatch plotting, or plotting
01064  * parallel to either the x-axis or the y-axis.
01065 \*--------------------------------------------------------------------------*/
01066 
01067 static void
01068 plt3zz(PLINT x0, PLINT y0, PLINT dx, PLINT dy, PLINT flag, PLINT *init,
01069        PLFLT *x, PLFLT *y, PLFLT **z, PLINT nx, PLINT ny, 
01070        PLINT *u, PLINT *v, PLFLT* c)
01071 {
01072     PLINT n = 0;
01073     PLFLT x2d, y2d;
01074 
01075     while (1 <= x0 && x0 <= nx && 1 <= y0 && y0 <= ny) {
01076        x2d = plP_w3wcx(x[x0 - 1], y[y0 - 1], z[x0 - 1][y0 - 1]);
01077        y2d = plP_w3wcy(x[x0 - 1], y[y0 - 1], z[x0 - 1][y0 - 1]);
01078        u[n] = plP_wcpcx(x2d);
01079        v[n] = plP_wcpcy(y2d);
01080        if (c != NULL)
01081          c[n] = (z[x0 - 1][y0 - 1] - fc_minz)/(fc_maxz-fc_minz);
01082 
01083        switch (flag) {
01084        case -3:
01085            y0 += dy;
01086            flag = -flag;
01087            break;
01088        case -2:
01089            y0 += dy;
01090            break;
01091        case -1:
01092            y0 += dy;
01093            flag = -flag;
01094            break;
01095        case 1:
01096            x0 += dx;
01097            break;
01098        case 2:
01099            x0 += dx;
01100            flag = -flag;
01101            break;
01102        case 3:
01103            x0 += dx;
01104            flag = -flag;
01105            break;
01106        }
01107        n++;
01108     }
01109 
01110     if (flag == 1 || flag == -2) {
01111        if (flag == 1) {
01112            x0 -= dx;
01113            y0 += dy;
01114        }
01115        else if (flag == -2) {
01116            y0 -= dy;
01117            x0 += dx;
01118        }
01119        if (1 <= x0 && x0 <= nx && 1 <= y0 && y0 <= ny) {
01120            x2d = plP_w3wcx( x[x0 - 1], y[y0 - 1], z[x0 - 1][y0 - 1]);
01121            y2d = plP_w3wcy( x[x0 - 1], y[y0 - 1], z[x0 - 1][y0 - 1]);
01122            u[n] = plP_wcpcx(x2d);
01123            v[n] = plP_wcpcy(y2d);
01124            if (c != NULL)
01125              c[n] = (z[x0 - 1][y0 - 1] - fc_minz)/(fc_maxz-fc_minz);
01126            n++;
01127        }
01128     }
01129 
01130 /* All the setup is done.  Time to do the work. */
01131 
01132     plnxtv(u, v, c, n, *init);
01133     *init = 0;
01134 }
01135 
01136 /*--------------------------------------------------------------------------*\
01137  * void plside3()
01138  *
01139  * This routine draws sides around the front of the 3d plot so that
01140  * it does not appear to float.
01141 \*--------------------------------------------------------------------------*/
01142 
01143 static void
01144 plside3(PLFLT *x, PLFLT *y, PLFLT **z, PLINT nx, PLINT ny, PLINT opt)
01145 {
01146     PLINT i;
01147     PLFLT cxx, cxy, cyx, cyy, cyz;
01148     PLFLT xmin, ymin, zmin, xmax, ymax, zmax, zscale;
01149     PLFLT tx, ty, ux, uy;
01150 
01151     plP_gw3wc(&cxx, &cxy, &cyx, &cyy, &cyz);
01152     plP_gdom(&xmin, &xmax, &ymin, &ymax);
01153     plP_grange(&zscale, &zmin, &zmax);
01154 
01155 /* Get x, y coordinates of legs and plot */
01156 
01157     if (cxx >= 0.0 && cxy <= 0.0) {
01158        if (opt != 1) {
01159            for (i = 0; i < nx; i++) {
01160               tx = plP_w3wcx(x[i], y[0], zmin);
01161               ty = plP_w3wcy(x[i], y[0], zmin);
01162               ux = plP_w3wcx(x[i], y[0], z[i][0]);
01163               uy = plP_w3wcy(x[i], y[0], z[i][0]);
01164               pljoin(tx, ty, ux, uy);
01165            }
01166        }
01167        if (opt != 2) {
01168            for (i = 0; i < ny; i++) {
01169               tx = plP_w3wcx(x[0], y[i], zmin);
01170               ty = plP_w3wcy(x[0], y[i], zmin);
01171               ux = plP_w3wcx(x[0], y[i], z[0][i]);
01172               uy = plP_w3wcy(x[0], y[i], z[0][i]);
01173               pljoin(tx, ty, ux, uy);
01174            }
01175        }
01176     }
01177     else if (cxx <= 0.0 && cxy <= 0.0) {
01178        if (opt != 1) {
01179            for (i = 0; i < nx; i++) {
01180               tx = plP_w3wcx(x[i], y[ny - 1], zmin);
01181               ty = plP_w3wcy(x[i], y[ny - 1], zmin);
01182               ux = plP_w3wcx(x[i], y[ny - 1], z[i][ny - 1]);
01183               uy = plP_w3wcy(x[i], y[ny - 1], z[i][ny - 1]);
01184               pljoin(tx, ty, ux, uy);
01185            }
01186        }
01187        if (opt != 2) {
01188            for (i = 0; i < ny; i++) {
01189               tx = plP_w3wcx(x[0], y[i], zmin);
01190               ty = plP_w3wcy(x[0], y[i], zmin);
01191               ux = plP_w3wcx(x[0], y[i], z[0][i]);
01192               uy = plP_w3wcy(x[0], y[i], z[0][i]);
01193               pljoin(tx, ty, ux, uy);
01194            }
01195        }
01196     }
01197     else if (cxx <= 0.0 && cxy >= 0.0) {
01198        if (opt != 1) {
01199            for (i = 0; i < nx; i++) {
01200               tx = plP_w3wcx(x[i], y[ny - 1], zmin);
01201               ty = plP_w3wcy(x[i], y[ny - 1], zmin);
01202               ux = plP_w3wcx(x[i], y[ny - 1], z[i][ny - 1]);
01203               uy = plP_w3wcy(x[i], y[ny - 1], z[i][ny - 1]);
01204               pljoin(tx, ty, ux, uy);
01205            }
01206        }
01207        if (opt != 2) {
01208            for (i = 0; i < ny; i++) {
01209               tx = plP_w3wcx(x[nx - 1], y[i], zmin);
01210               ty = plP_w3wcy(x[nx - 1], y[i], zmin);
01211               ux = plP_w3wcx(x[nx - 1], y[i], z[nx - 1][i]);
01212               uy = plP_w3wcy(x[nx - 1], y[i], z[nx - 1][i]);
01213               pljoin(tx, ty, ux, uy);
01214            }
01215        }
01216     }
01217     else if (cxx >= 0.0 && cxy >= 0.0) {
01218        if (opt != 1) {
01219            for (i = 0; i < nx; i++) {
01220               tx = plP_w3wcx(x[i], y[0], zmin);
01221               ty = plP_w3wcy(x[i], y[0], zmin);
01222               ux = plP_w3wcx(x[i], y[0], z[i][0]);
01223               uy = plP_w3wcy(x[i], y[0], z[i][0]);
01224               pljoin(tx, ty, ux, uy);
01225            }
01226        }
01227        if (opt != 2) {
01228            for (i = 0; i < ny; i++) {
01229               tx = plP_w3wcx(x[nx - 1], y[i], zmin);
01230               ty = plP_w3wcy(x[nx - 1], y[i], zmin);
01231               ux = plP_w3wcx(x[nx - 1], y[i], z[nx - 1][i]);
01232               uy = plP_w3wcy(x[nx - 1], y[i], z[nx - 1][i]);
01233               pljoin(tx, ty, ux, uy);
01234            }
01235        }
01236     }
01237 }
01238 
01239 /*--------------------------------------------------------------------------*\
01240  * void plgrid3()
01241  *
01242  * This routine draws a grid around the back side of the 3d plot with
01243  * hidden line removal.
01244 \*--------------------------------------------------------------------------*/
01245 
01246 static void
01247 plgrid3(PLFLT tick)
01248 {
01249     PLFLT xmin, ymin, zmin, xmax, ymax, zmax, zscale;
01250     PLFLT cxx, cxy, cyx, cyy, cyz, zmin_in, zmax_in;
01251     PLINT u[3], v[3];
01252     PLINT nsub = 0;
01253     PLFLT tp;
01254 
01255     plP_gw3wc(&cxx, &cxy, &cyx, &cyy, &cyz);
01256     plP_gdom(&xmin, &xmax, &ymin, &ymax);
01257     plP_grange(&zscale, &zmin_in, &zmax_in);
01258     zmin = (zmax_in > zmin_in) ? zmin_in: zmax_in;
01259     zmax = (zmax_in > zmin_in) ? zmax_in: zmin_in;
01260 
01261     pldtik(zmin, zmax, &tick, &nsub);
01262     tp = tick * floor(zmin / tick) + tick;
01263     pl3upv = 0;
01264 
01265     if (cxx >= 0.0 && cxy <= 0.0) {
01266        while (tp <= zmax) {
01267            u[0] = plP_wcpcx(plP_w3wcx(xmin, ymax, tp));
01268            v[0] = plP_wcpcy(plP_w3wcy(xmin, ymax, tp));
01269            u[1] = plP_wcpcx(plP_w3wcx(xmax, ymax, tp));
01270            v[1] = plP_wcpcy(plP_w3wcy(xmax, ymax, tp));
01271            u[2] = plP_wcpcx(plP_w3wcx(xmax, ymin, tp));
01272            v[2] = plP_wcpcy(plP_w3wcy(xmax, ymin, tp));
01273            plnxtv(u, v, 0, 3, 0);
01274 
01275            tp += tick;
01276        }
01277        u[0] = plP_wcpcx(plP_w3wcx(xmax, ymax, zmin));
01278        v[0] = plP_wcpcy(plP_w3wcy(xmax, ymax, zmin));
01279        u[1] = plP_wcpcx(plP_w3wcx(xmax, ymax, zmax));
01280        v[1] = plP_wcpcy(plP_w3wcy(xmax, ymax, zmax));
01281        plnxtv(u, v, 0, 2, 0);
01282     }
01283     else if (cxx <= 0.0 && cxy <= 0.0) {
01284        while (tp <= zmax) {
01285            u[0] = plP_wcpcx(plP_w3wcx(xmax, ymax, tp));
01286            v[0] = plP_wcpcy(plP_w3wcy(xmax, ymax, tp));
01287            u[1] = plP_wcpcx(plP_w3wcx(xmax, ymin, tp));
01288            v[1] = plP_wcpcy(plP_w3wcy(xmax, ymin, tp));
01289            u[2] = plP_wcpcx(plP_w3wcx(xmin, ymin, tp));
01290            v[2] = plP_wcpcy(plP_w3wcy(xmin, ymin, tp));
01291            plnxtv(u, v, 0,3, 0);
01292 
01293            tp += tick;
01294        }
01295        u[0] = plP_wcpcx(plP_w3wcx(xmax, ymin, zmin));
01296        v[0] = plP_wcpcy(plP_w3wcy(xmax, ymin, zmin));
01297        u[1] = plP_wcpcx(plP_w3wcx(xmax, ymin, zmax));
01298        v[1] = plP_wcpcy(plP_w3wcy(xmax, ymin, zmax));
01299        plnxtv(u, v, 0,2, 0);
01300     }
01301     else if (cxx <= 0.0 && cxy >= 0.0) {
01302        while (tp <= zmax) {
01303            u[0] = plP_wcpcx(plP_w3wcx(xmax, ymin, tp));
01304            v[0] = plP_wcpcy(plP_w3wcy(xmax, ymin, tp));
01305            u[1] = plP_wcpcx(plP_w3wcx(xmin, ymin, tp));
01306            v[1] = plP_wcpcy(plP_w3wcy(xmin, ymin, tp));
01307            u[2] = plP_wcpcx(plP_w3wcx(xmin, ymax, tp));
01308            v[2] = plP_wcpcy(plP_w3wcy(xmin, ymax, tp));
01309            plnxtv(u, v, 0,3, 0);
01310 
01311            tp += tick;
01312        }
01313        u[0] = plP_wcpcx(plP_w3wcx(xmin, ymin, zmin));
01314        v[0] = plP_wcpcy(plP_w3wcy(xmin, ymin, zmin));
01315        u[1] = plP_wcpcx(plP_w3wcx(xmin, ymin, zmax));
01316        v[1] = plP_wcpcy(plP_w3wcy(xmin, ymin, zmax));
01317        plnxtv(u, v, 0,2, 0);
01318     }
01319     else if (cxx >= 0.0 && cxy >= 0.0) {
01320        while (tp <= zmax) {
01321            u[0] = plP_wcpcx(plP_w3wcx(xmin, ymin, tp));
01322            v[0] = plP_wcpcy(plP_w3wcy(xmin, ymin, tp));
01323            u[1] = plP_wcpcx(plP_w3wcx(xmin, ymax, tp));
01324            v[1] = plP_wcpcy(plP_w3wcy(xmin, ymax, tp));
01325            u[2] = plP_wcpcx(plP_w3wcx(xmax, ymax, tp));
01326            v[2] = plP_wcpcy(plP_w3wcy(xmax, ymax, tp));
01327            plnxtv(u, v, 0,3, 0);
01328 
01329            tp += tick;
01330        }
01331        u[0] = plP_wcpcx(plP_w3wcx(xmin, ymax, zmin));
01332        v[0] = plP_wcpcy(plP_w3wcy(xmin, ymax, zmin));
01333        u[1] = plP_wcpcx(plP_w3wcx(xmin, ymax, zmax));
01334        v[1] = plP_wcpcy(plP_w3wcy(xmin, ymax, zmax));
01335        plnxtv(u, v, 0,2, 0);
01336     }
01337     pl3upv = 1;
01338 }
01339 
01340 /*--------------------------------------------------------------------------*\
01341  * void plnxtv()
01342  *
01343  * Draw the next view of a 3-d plot. The physical coordinates of the
01344  * points for the next view are placed in the n points of arrays u and
01345  * v. The silhouette found so far is stored in the heap as a set of m peak
01346  * points.
01347  *
01348  * These routines dynamically allocate memory for hidden line removal.
01349  * Memory is allocated in blocks of 2*BINC*sizeof(PLINT) bytes.  Large
01350  * values of BINC give better performance but also allocate more memory
01351  * than is needed. If your 3D plots are very "spiky" or you are working
01352  * with very large matrices then you will probably want to increase BINC.
01353 \*--------------------------------------------------------------------------*/
01354 
01355 static void
01356 plnxtv(PLINT *u, PLINT *v, PLFLT* c, PLINT n, PLINT init)
01357 {
01358     plnxtvhi(u, v, c, n, init);
01359 
01360     if (pl3mode)
01361        plnxtvlo(u, v, c, n, init);
01362 }
01363 
01364 /*--------------------------------------------------------------------------*\
01365  * void plnxtvhi()
01366  *
01367  * Draw the top side of the 3-d plot.
01368 \*--------------------------------------------------------------------------*/
01369 
01370 static void
01371 plnxtvhi(PLINT *u, PLINT *v, PLFLT* c, PLINT n, PLINT init)
01372 {
01373   /*
01374    * For the initial set of points, just display them and store them as the
01375    * peak points.
01376    */
01377   if (init == 1) {
01378     int i;
01379     oldhiview = (PLINT *) malloc((size_t) (2 * n * sizeof(PLINT)));
01380     if ( ! oldhiview)
01381       myexit("plnxtvhi: Out of memory.");
01382 
01383     oldhiview[0] = u[0];
01384     oldhiview[1] = v[0];
01385     plP_draw3d(u[0], v[0], c, 0, 1);
01386     for (i = 1; i < n; i++) {
01387       oldhiview[2 * i] = u[i];
01388       oldhiview[2 * i + 1] = v[i];
01389       plP_draw3d(u[i], v[i], c, i, 0);
01390     }
01391     mhi = n;
01392     return;
01393   }
01394 
01395   /*
01396    * Otherwise, we need to consider hidden-line removal problem. We scan
01397    * over the points in both the old (i.e. oldhiview[]) and new (i.e. u[]
01398    * and v[]) arrays in order of increasing x coordinate.  At each stage, we
01399    * find the line segment in the other array (if one exists) that straddles
01400    * the x coordinate of the point. We have to determine if the point lies
01401    * above or below the line segment, and to check if the below/above status
01402    * has changed since the last point.
01403    *
01404    * If pl3upv = 0 we do not update the view, this is useful for drawing
01405    * lines on the graph after we are done plotting points.  Hidden line
01406    * removal is still done, but the view is not updated.
01407    */
01408   xxhi = 0;
01409   if (pl3upv != 0) {
01410     newhisize = 2 * (mhi + BINC);
01411     if (newhiview != NULL) {
01412       newhiview = 
01413        (PLINT *) realloc((void *) newhiview,
01414                        (size_t) (newhisize * sizeof(PLINT)));
01415     }
01416     else {
01417       newhiview = 
01418        (PLINT *) malloc((size_t) (newhisize * sizeof(PLINT)));
01419     }
01420     if ( ! newhiview)
01421       myexit("plnxtvhi: Out of memory.");
01422   }
01423 
01424   /* Do the draw or shading with hidden line removal */
01425 
01426   plnxtvhi_draw(u, v, c, n);
01427 
01428   /* Set oldhiview */
01429 
01430   swaphiview();
01431 }
01432 
01433 /*--------------------------------------------------------------------------*\
01434  * void plnxtvhi_draw()
01435  *
01436  * Draw the top side of the 3-d plot.
01437 \*--------------------------------------------------------------------------*/
01438 
01439 static void
01440 plnxtvhi_draw(PLINT *u, PLINT *v, PLFLT* c, PLINT n)
01441 {
01442     PLINT i = 0, j = 0, first = 1;
01443     PLINT sx1 = 0, sx2 = 0, sy1 = 0, sy2 = 0;
01444     PLINT su1, su2, sv1, sv2;
01445     PLINT cx, cy, px, py;
01446     PLINT seg, ptold, lstold = 0, pthi, pnewhi = 0, newhi, change, ochange = 0;
01447 
01448 /*
01449  * (oldhiview[2*i], oldhiview[2*i]) is the i'th point in the old array
01450  * (u[j], v[j]) is the j'th point in the new array
01451  */
01452 
01453 /* 
01454  * First attempt at 3d shading.  It works ok for simple plots, but
01455  * will just not draw faces, or draw them overlapping for very
01456  * jagged plots
01457  */
01458 
01459     while (i < mhi || j < n) {
01460 
01461     /*
01462      * The coordinates of the point under consideration are (px,py).  The
01463      * line segment joins (sx1,sy1) to (sx2,sy2).  "ptold" is true if the
01464      * point lies in the old array. We set it by comparing the x coordinates
01465      * of the i'th old point and the j'th new point, being careful if we
01466      * have fallen past the edges. Having found the point, load up the point
01467      * and segment coordinates appropriately.
01468      */
01469 
01470     ptold = (j >= n || (i < mhi && oldhiview[2 * i] < u[j]));
01471     if (ptold) {
01472       px = oldhiview[2 * i];
01473       py = oldhiview[2 * i + 1];
01474       seg = j > 0 && j < n;
01475       if (seg) {
01476        sx1 = u[j - 1];
01477        sy1 = v[j - 1];
01478        sx2 = u[j];
01479        sy2 = v[j];
01480       }
01481     } else {
01482       px = u[j];
01483       py = v[j];
01484       seg = i > 0 && i < mhi;
01485       if (seg) {
01486        sx1 = oldhiview[2 * (i - 1)];
01487        sy1 = oldhiview[2 * (i - 1) + 1];
01488        sx2 = oldhiview[2 * i];
01489        sy2 = oldhiview[2 * i + 1];
01490       }
01491     }
01492 
01493     /*
01494      * Now determine if the point is higher than the segment, using the
01495      * logical function "above". We also need to know if it is the old view
01496      * or the new view that is higher. "newhi" is set true if the new view
01497      * is higher than the old.
01498      */
01499     if (seg)
01500       pthi = plabv(px, py, sx1, sy1, sx2, sy2);
01501     else
01502       pthi = 1;
01503 
01504     newhi = (ptold && !pthi) || (!ptold && pthi);
01505     /*
01506      * The last point and this point lie on different sides of
01507      * the current silouette
01508      */
01509     change = (newhi && !pnewhi) || (!newhi && pnewhi);
01510 
01511     /*
01512      * There is a new intersection point to put in the peak array if the
01513      * state of "newhi" changes.
01514      */
01515     if (first) {
01516       plP_draw3d(px, py, c, j, 1);
01517       first = 0;
01518       lstold = ptold;
01519       savehipoint(px, py);
01520       pthi = 0;
01521       ochange = 0;
01522     } else if (change) {
01523 
01524       /*
01525        * Take care of special cases at end of arrays.  If pl3upv is 0 the
01526        * endpoints are not connected to the old view.
01527        */
01528       if (pl3upv == 0 && ((!ptold && j == 0) || (ptold && i == 0))) {
01529        plP_draw3d(px, py, c, j, 1);
01530        lstold = ptold;
01531        pthi = 0;
01532        ochange = 0;
01533       } else if (pl3upv == 0 &&
01534                (( ! ptold && i >= mhi) || (ptold && j >= n))) {
01535        plP_draw3d(px, py, c, j, 1);
01536        lstold = ptold;
01537        pthi = 0;
01538        ochange = 0;
01539       } else {
01540 
01541        /*
01542         * If pl3upv is not 0 then we do want to connect the current line
01543         * with the previous view at the endpoints.  Also find intersection
01544         * point with old view.
01545         */
01546        if (i == 0) {
01547          sx1 = oldhiview[0];
01548          sy1 = -1;
01549          sx2 = oldhiview[0];
01550          sy2 = oldhiview[1];
01551        } else if (i >= mhi) {
01552          sx1 = oldhiview[2 * (mhi - 1)];
01553          sy1 = oldhiview[2 * (mhi - 1) + 1];
01554          sx2 = oldhiview[2 * (mhi - 1)];
01555          sy2 = -1;
01556        } else {
01557          sx1 = oldhiview[2 * (i - 1)];
01558          sy1 = oldhiview[2 * (i - 1) + 1];
01559          sx2 = oldhiview[2 * i];
01560          sy2 = oldhiview[2 * i + 1];
01561        }
01562 
01563        if (j == 0) {
01564          su1 = u[0];
01565          sv1 = -1;
01566          su2 = u[0];
01567          sv2 = v[0];
01568        } else if (j >= n) {
01569          su1 = u[n - 1];
01570          sv1 = v[n - 1];
01571          su2 = u[n - 1];
01572          sv2 = -1;
01573        } else {
01574          su1 = u[j - 1];
01575          sv1 = v[j - 1];
01576          su2 = u[j];
01577          sv2 = v[j];
01578        }
01579 
01580        /* Determine the intersection */
01581 
01582        pl3cut(sx1, sy1, sx2, sy2, su1, sv1, su2, sv2, &cx, &cy);
01583        if (cx == px && cy == py) {
01584          if (lstold && !ochange)
01585            plP_draw3d(px, py, c, j, 1);
01586          else
01587            plP_draw3d(px, py, c, j, 0);
01588 
01589          savehipoint(px, py);
01590          lstold = 1;
01591          pthi = 0;
01592        } else {
01593          if (lstold && !ochange)
01594            plP_draw3d(cx, cy, c, j, 1);
01595          else
01596            plP_draw3d(cx, cy, c, j, 0);
01597 
01598          lstold = 1;
01599          savehipoint(cx, cy);
01600        }
01601        ochange = 1;
01602       }
01603     }
01604 
01605     /* If point is high then draw plot to point and update view. */
01606 
01607     if (pthi) {
01608       if (lstold && ptold)
01609        plP_draw3d(px, py, c, j, 1);
01610       else
01611        plP_draw3d(px, py, c, j, 0);
01612 
01613       savehipoint(px, py);
01614       lstold = ptold;
01615       ochange = 0;
01616     }
01617     pnewhi = newhi;
01618 
01619     if (ptold)
01620       i++;
01621     else
01622       j++;
01623   }
01624 }
01625 
01626 /*--------------------------------------------------------------------------*\
01627  * void  plP_draw3d()
01628  *
01629  * Does a simple move or line draw.
01630 \*--------------------------------------------------------------------------*/
01631 
01632 static void
01633 plP_draw3d(PLINT x, PLINT y, PLFLT *c, PLINT j, PLINT move)
01634 {
01635     if (move)
01636       plP_movphy(x, y);
01637     else {
01638       if (c != NULL)
01639        plcol1(c[j-1]);
01640       plP_draphy(x, y);
01641     }
01642 }
01643 
01644 /*--------------------------------------------------------------------------*\
01645  * void plnxtvlo()
01646  *
01647  * Draw the bottom side of the 3-d plot.
01648 \*--------------------------------------------------------------------------*/
01649 
01650 static void
01651 plnxtvlo(PLINT *u, PLINT *v, PLFLT*c, PLINT n, PLINT init)
01652 {
01653   PLINT i, j, first;
01654   PLINT sx1 = 0, sx2 = 0, sy1 = 0, sy2 = 0;
01655   PLINT su1, su2, sv1, sv2;
01656   PLINT cx, cy, px, py;
01657   PLINT seg, ptold, lstold = 0, ptlo, pnewlo, newlo, change, ochange = 0;
01658 
01659   first = 1;
01660   pnewlo = 0;
01661 
01662   /*
01663    * For the initial set of points, just display them and store them as the
01664    * peak points.
01665    */
01666   if (init == 1) {
01667 
01668     oldloview = (PLINT *) malloc((size_t) (2 * n * sizeof(PLINT)));
01669     if ( ! oldloview)
01670       myexit("\nplnxtvlo: Out of memory.");
01671 
01672     plP_draw3d(u[0], v[0], c, 0, 1);
01673     oldloview[0] = u[0];
01674     oldloview[1] = v[0];
01675     for (i = 1; i < n; i++) {
01676       plP_draw3d(u[i], v[i], c, i, 0);
01677       oldloview[2 * i] = u[i];
01678       oldloview[2 * i + 1] = v[i];
01679     }
01680     mlo = n;
01681     return;
01682   }
01683 
01684   /*
01685    * Otherwise, we need to consider hidden-line removal problem. We scan
01686    * over the points in both the old (i.e. oldloview[]) and new (i.e. u[]
01687    * and v[]) arrays in order of increasing x coordinate.  At each stage, we
01688    * find the line segment in the other array (if one exists) that straddles
01689    * the x coordinate of the point. We have to determine if the point lies
01690    * above or below the line segment, and to check if the below/above status
01691    * has changed since the last point.
01692    *
01693    * If pl3upv = 0 we do not update the view, this is useful for drawing
01694    * lines on the graph after we are done plotting points.  Hidden line
01695    * removal is still done, but the view is not updated.
01696    */
01697   xxlo = 0;
01698   i = 0;
01699   j = 0;
01700   if (pl3upv != 0) {
01701     newlosize = 2 * (mlo + BINC);
01702     if (newloview != NULL) {
01703       newloview = 
01704        (PLINT *) realloc((void *) newloview,
01705                        (size_t) (newlosize * sizeof(PLINT)));
01706     }
01707     else {
01708       newloview = 
01709        (PLINT *) malloc((size_t) (newlosize * sizeof(PLINT)));
01710     }
01711     if ( ! newloview)
01712       myexit("plnxtvlo: Out of memory.");
01713   }
01714 
01715   /*
01716    * (oldloview[2*i], oldloview[2*i]) is the i'th point in the old array
01717    * (u[j], v[j]) is the j'th point in the new array.
01718    */
01719   while (i < mlo || j < n) {
01720 
01721     /*
01722      * The coordinates of the point under consideration are (px,py).  The
01723      * line segment joins (sx1,sy1) to (sx2,sy2).  "ptold" is true if the
01724      * point lies in the old array. We set it by comparing the x coordinates
01725      * of the i'th old point and the j'th new point, being careful if we
01726      * have fallen past the edges. Having found the point, load up the point
01727      * and segment coordinates appropriately.
01728      */
01729     ptold = (j >= n || (i < mlo && oldloview[2 * i] < u[j]));
01730     if (ptold) {
01731       px = oldloview[2 * i];
01732       py = oldloview[2 * i + 1];
01733       seg = j > 0 && j < n;
01734       if (seg) {
01735        sx1 = u[j - 1];
01736        sy1 = v[j - 1];
01737        sx2 = u[j];
01738        sy2 = v[j];
01739       }
01740     }
01741     else {
01742       px = u[j];
01743       py = v[j];
01744       seg = i > 0 && i < mlo;
01745       if (seg) {
01746        sx1 = oldloview[2 * (i - 1)];
01747        sy1 = oldloview[2 * (i - 1) + 1];
01748        sx2 = oldloview[2 * i];
01749        sy2 = oldloview[2 * i + 1];
01750       }
01751     }
01752 
01753     /*
01754      * Now determine if the point is lower than the segment, using the
01755      * logical function "above". We also need to know if it is the old view
01756      * or the new view that is lower. "newlo" is set true if the new view is
01757      * lower than the old.
01758      */
01759     if (seg)
01760       ptlo = !plabv(px, py, sx1, sy1, sx2, sy2);
01761     else
01762       ptlo = 1;
01763 
01764     newlo = (ptold && !ptlo) || (!ptold && ptlo);
01765     change = (newlo && !pnewlo) || (!newlo && pnewlo);
01766 
01767     /*
01768      * There is a new intersection point to put in the peak array if the
01769      * state of "newlo" changes.
01770      */
01771     if (first) {
01772       plP_draw3d(px, py, c, j, 1);
01773       first = 0;
01774       lstold = ptold;
01775       savelopoint(px, py);
01776       ptlo = 0;
01777       ochange = 0;
01778     }
01779     else if (change) {
01780 
01781       /*
01782        * Take care of special cases at end of arrays.  If pl3upv is 0 the
01783        * endpoints are not connected to the old view.
01784        */
01785       if (pl3upv == 0 && ((!ptold && j == 0) || (ptold && i == 0))) {
01786        plP_draw3d(px, py, c, j, 1);
01787        lstold = ptold;
01788        ptlo = 0;
01789        ochange = 0;
01790       }
01791       else if (pl3upv == 0 &&
01792               (( ! ptold && i >= mlo) || (ptold && j >= n))) {
01793        plP_draw3d(px, py, c, j, 1);
01794        lstold = ptold;
01795        ptlo = 0;
01796        ochange = 0;
01797       }
01798 
01799       /*
01800        * If pl3upv is not 0 then we do want to connect the current line
01801        * with the previous view at the endpoints.  Also find intersection
01802        * point with old view.
01803        */
01804       else {
01805        if (i == 0) {
01806          sx1 = oldloview[0];
01807          sy1 = 100000;
01808          sx2 = oldloview[0];
01809          sy2 = oldloview[1];
01810        }
01811        else if (i >= mlo) {
01812          sx1 = oldloview[2 * (mlo - 1)];
01813          sy1 = oldloview[2 * (mlo - 1) + 1];
01814          sx2 = oldloview[2 * (mlo - 1)];
01815          sy2 = 100000;
01816        }
01817        else {
01818          sx1 = oldloview[2 * (i - 1)];
01819          sy1 = oldloview[2 * (i - 1) + 1];
01820          sx2 = oldloview[2 * i];
01821          sy2 = oldloview[2 * i + 1];
01822        }
01823 
01824        if (j == 0) {
01825          su1 = u[0];
01826          sv1 = 100000;
01827          su2 = u[0];
01828          sv2 = v[0];
01829        }
01830        else if (j >= n) {
01831          su1 = u[n - 1];
01832          sv1 = v[n - 1];
01833          su2 = u[n];
01834          sv2 = 100000;
01835        }
01836        else {
01837          su1 = u[j - 1];
01838          sv1 = v[j - 1];
01839          su2 = u[j];
01840          sv2 = v[j];
01841        }
01842 
01843        /* Determine the intersection */
01844 
01845        pl3cut(sx1, sy1, sx2, sy2, su1, sv1, su2, sv2, &cx, &cy);
01846        if (cx == px && cy == py) {
01847          if (lstold && !ochange)
01848            plP_draw3d(px, py, c, j, 1);
01849          else
01850            plP_draw3d(px, py, c, j, 0);
01851 
01852          savelopoint(px, py);
01853          lstold = 1;
01854          ptlo = 0;
01855        }
01856        else {
01857          if (lstold && !ochange)
01858            plP_draw3d(cx, cy, c, j, 1);
01859          else
01860            plP_draw3d(cx, cy, c, j, 0);
01861 
01862          lstold = 1;
01863          savelopoint(cx, cy);
01864        }
01865        ochange = 1;
01866       }
01867     }
01868 
01869     /* If point is low then draw plot to point and update view. */
01870 
01871     if (ptlo) {
01872       if (lstold && ptold)
01873        plP_draw3d(px, py, c, j, 1);
01874       else
01875        plP_draw3d(px, py, c, j, 0);
01876 
01877       savelopoint(px, py);
01878       lstold = ptold;
01879       ochange = 0;
01880     }
01881 
01882     pnewlo = newlo;
01883 
01884     if (ptold)
01885       i = i + 1;
01886     else
01887       j = j + 1;
01888   }
01889 
01890   /* Set oldloview */
01891 
01892   swaploview();
01893 }
01894 
01895 /*--------------------------------------------------------------------------*\
01896  * savehipoint
01897  * savelopoint
01898  *
01899  * Add a point to the list of currently visible peaks/valleys, when
01900  * drawing the top/bottom surface, respectively.
01901 \*--------------------------------------------------------------------------*/
01902 
01903 static void
01904 savehipoint(PLINT px, PLINT py)
01905 {
01906     if (pl3upv == 0)
01907        return;
01908 
01909     if (xxhi >= newhisize) {       /* allocate additional space */
01910        newhisize += 2 * BINC;
01911        newhiview = (PLINT *) realloc((void *) newhiview,
01912                                   (size_t) (newhisize * sizeof(PLINT)));
01913        if ( ! newhiview) 
01914            myexit("savehipoint: Out of memory.");
01915     }
01916 
01917     newhiview[xxhi] = px;
01918     xxhi++;
01919     newhiview[xxhi] = py;
01920     xxhi++;
01921 }
01922 
01923 static void
01924 savelopoint(PLINT px, PLINT py)
01925 {
01926     if (pl3upv == 0)
01927        return;
01928 
01929     if (xxlo >= newlosize) {       /* allocate additional space */
01930        newlosize += 2 * BINC;
01931        newloview = (PLINT *) realloc((void *) newloview,
01932                                   (size_t) (newlosize * sizeof(PLINT)));
01933        if ( ! newloview)
01934            myexit("savelopoint: Out of memory.");
01935     }
01936 
01937     newloview[xxlo] = px;
01938     xxlo++;
01939     newloview[xxlo] = py;
01940     xxlo++;
01941 }
01942 
01943 /*--------------------------------------------------------------------------*\
01944  * swaphiview
01945  * swaploview
01946  *
01947  * Swaps the top/bottom views.  Need to do a real swap so that the 
01948  * memory cleanup routine really frees everything (and only once).
01949 \*--------------------------------------------------------------------------*/
01950 
01951 static void
01952 swaphiview(void)
01953 {
01954     PLINT *tmp;
01955 
01956     if (pl3upv != 0) {
01957        mhi = xxhi / 2;
01958        tmp = oldhiview;
01959        oldhiview = newhiview;
01960        newhiview = tmp;
01961     }
01962 }
01963 
01964 static void
01965 swaploview(void)
01966 {
01967     PLINT *tmp;
01968 
01969     if (pl3upv != 0) {
01970        mlo = xxlo / 2;
01971        tmp = oldloview;
01972        oldloview = newloview;
01973        newloview = tmp;
01974     }
01975 }
01976 
01977 /*--------------------------------------------------------------------------*\
01978  * freework
01979  *
01980  * Frees memory associated with work arrays
01981 \*--------------------------------------------------------------------------*/
01982 
01983 static void
01984 freework(void)
01985 {
01986     free_mem(oldhiview);
01987     free_mem(oldloview);
01988     free_mem(newhiview);
01989     free_mem(newloview);
01990     free_mem(vtmp);
01991     free_mem(utmp);
01992     free_mem(ctmp);
01993 }
01994 
01995 /*--------------------------------------------------------------------------*\
01996  * myexit
01997  *
01998  * Calls plexit, cleaning up first.
01999 \*--------------------------------------------------------------------------*/
02000 
02001 static void
02002 myexit(char *msg)
02003 {
02004     freework();
02005     plexit(msg);
02006 }
02007 
02008 /*--------------------------------------------------------------------------*\
02009  * myabort
02010  *
02011  * Calls plabort, cleaning up first.
02012  * Caller should return to the user program.
02013 \*--------------------------------------------------------------------------*/
02014 
02015 static void
02016 myabort(char *msg)
02017 {
02018     freework();
02019     plabort(msg);
02020 }
02021 
02022 /*--------------------------------------------------------------------------*\
02023  * int plabv()
02024  *
02025  * Determines if point (px,py) lies above the line joining (sx1,sy1) to
02026  * (sx2,sy2). It only works correctly if sx1 <= px <= sx2.
02027 \*--------------------------------------------------------------------------*/
02028 
02029 static int
02030 plabv(PLINT px, PLINT py, PLINT sx1, PLINT sy1, PLINT sx2, PLINT sy2)
02031 {
02032     int above;
02033 
02034     if (py >= sy1 && py >= sy2)
02035        above = 1;
02036     else if (py < sy1 && py < sy2)
02037        above = 0;
02038     else if ((double) (sx2 - sx1) * (py - sy1) >=
02039             (double) (px - sx1) * (sy2 - sy1))
02040        above = 1;
02041     else
02042        above = 0;
02043 
02044     return above;
02045 }
02046 
02047 /*--------------------------------------------------------------------------*\
02048  * void pl3cut()
02049  *
02050  * Determines the point of intersection (cx,cy) between the line joining
02051  * (sx1,sy1) to (sx2,sy2) and the line joining (su1,sv1) to (su2,sv2).
02052 \*--------------------------------------------------------------------------*/
02053 
02054 static void
02055 pl3cut(PLINT sx1, PLINT sy1, PLINT sx2, PLINT sy2,
02056        PLINT su1, PLINT sv1, PLINT su2, PLINT sv2, PLINT *cx, PLINT *cy)
02057 {
02058     PLINT x21, y21, u21, v21, yv1, xu1, a, b;
02059     double fa, fb;
02060 
02061     x21 = sx2 - sx1;
02062     y21 = sy2 - sy1;
02063     u21 = su2 - su1;
02064     v21 = sv2 - sv1;
02065     yv1 = sy1 - sv1;
02066     xu1 = sx1 - su1;
02067 
02068     a = x21 * v21 - y21 * u21;
02069     fa = (double) a;
02070     if (a == 0) {
02071        if (sx2 < su2) {
02072            *cx = sx2;
02073            *cy = sy2;
02074        }
02075        else {
02076            *cx = su2;
02077            *cy = sv2;
02078        }
02079        return;
02080     }
02081     else {
02082        b = yv1 * u21 - xu1 * v21;
02083        fb = (double) b;
02084        *cx = (PLINT) (sx1 + (fb * x21) / fa + .5);
02085        *cy = (PLINT) (sy1 + (fb * y21) / fa + .5);
02086     }
02087 }