/* $Id: plot3d.c,v 1.66 2006/05/12 23:53:41 airwin Exp $
3d plot routines.
Copyright (C) 2004 Alan W. Irwin
Copyright (C) 2004 Joao Cardoso
Copyright (C) 2004 Andrew Ross
This file is part of PLplot.
PLplot is free software; you can redistribute it and/or modify
it under the terms of the GNU Library General Public License as published
by the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
PLplot is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU Library General Public License for more details.
You should have received a copy of the GNU Library General Public License
along with PLplot; if not, write to the Free Software
Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
*/
#include "plplotP.h"
/* Internal constants */
#define BINC 50 /* Block size for memory allocation */
static PLINT pl3mode = 0; /* 0 3d solid; 1 mesh plot */
static PLINT pl3upv = 1; /* 1 update view; 0 no update */
static PLINT zbflg = 0, zbcol, zbwidth;
static PLFLT zbtck;
static PLINT *oldhiview = NULL;
static PLINT *oldloview = NULL;
static PLINT *newhiview = NULL;
static PLINT *newloview = NULL;
static PLINT *utmp = NULL;
static PLINT *vtmp = NULL;
static PLFLT *ctmp = NULL;
static PLINT mhi, xxhi, newhisize;
static PLINT mlo, xxlo, newlosize;
/* Light source for shading */
static PLFLT xlight, ylight, zlight;
static PLINT falsecolor=0;
static PLFLT fc_minz, fc_maxz;
/* Prototypes for static functions */
static void plgrid3 (PLFLT);
static void plnxtv (PLINT *, PLINT *, PLFLT*, PLINT, PLINT);
static void plside3 (PLFLT *, PLFLT *, PLFLT **, PLINT, PLINT, PLINT);
static void plt3zz (PLINT, PLINT, PLINT, PLINT,
PLINT, PLINT *, PLFLT *, PLFLT *, PLFLT **,
PLINT, PLINT, PLINT *, PLINT *, PLFLT*);
static void plnxtvhi (PLINT *, PLINT *, PLFLT*, PLINT, PLINT);
static void plnxtvlo (PLINT *, PLINT *, PLFLT*, PLINT, PLINT);
static void plnxtvhi_draw(PLINT *u, PLINT *v, PLFLT* c, PLINT n);
static void savehipoint (PLINT, PLINT);
static void savelopoint (PLINT, PLINT);
static void swaphiview (void);
static void swaploview (void);
static void myexit (char *);
static void myabort (char *);
static void freework (void);
static int plabv (PLINT, PLINT, PLINT, PLINT, PLINT, PLINT);
static void pl3cut (PLINT, PLINT, PLINT, PLINT, PLINT,
PLINT, PLINT, PLINT, PLINT *, PLINT *);
static PLFLT plGetAngleToLight(PLFLT* x, PLFLT* y, PLFLT* z);
static void plP_draw3d(PLINT x, PLINT y, PLFLT *c, PLINT j, PLINT move);
static void plxyindexlimits(PLINT instart, PLINT inn,
PLINT *inarray_min, PLINT *inarray_max,
PLINT *outstart, PLINT *outn, PLINT outnmax,
PLINT *outarray_min, PLINT *outarray_max);
/* #define MJL_HACK 1 */
#if MJL_HACK
static void plP_fill3(PLINT x0, PLINT y0, PLINT x1, PLINT y1,
PLINT x2, PLINT y2, PLINT j);
static void plP_fill4(PLINT x0, PLINT y0, PLINT x1, PLINT y1,
PLINT x2, PLINT y2, PLINT x3, PLINT y3, PLINT j);
#endif
/*--------------------------------------------------------------------------*\
* void plsetlightsource(x, y, z)
*
* Sets the position of the light source.
\*--------------------------------------------------------------------------*/
void
c_pllightsource(PLFLT x, PLFLT y, PLFLT z)
{
xlight = x;
ylight = y;
zlight = z;
}
/*--------------------------------------------------------------------------*\
* void plmesh(x, y, z, nx, ny, opt)
*
* Plots a mesh representation of the function z[x][y]. The x values
* are stored as x[0..nx-1], the y values as y[0..ny-1], and the
* z values are in the 2-d array z[][]. The integer "opt" specifies:
* see plmeshc() below.
\*--------------------------------------------------------------------------*/
void
c_plmesh(PLFLT *x, PLFLT *y, PLFLT **z, PLINT nx, PLINT ny, PLINT opt)
{
c_plot3dc(x, y, z, nx, ny, opt | MESH, NULL, 0);
}
/*--------------------------------------------------------------------------*\
* void plmeshc(x, y, z, nx, ny, opt, clevel, nlevel)
*
* Plots a mesh representation of the function z[x][y]. The x values
* are stored as x[0..nx-1], the y values as y[0..ny-1], and the
* z values are in the 2-d array z[][]. The integer "opt" specifies:
*
* DRAW_LINEX draw lines parallel to the X axis
* DRAW_LINEY draw lines parallel to the Y axis
* DRAW_LINEXY draw lines parallel to both the X and Y axis
* MAG_COLOR draw the mesh with a color dependent of the magnitude
* BASE_CONT draw contour plot at bottom xy plane
* TOP_CONT draw contour plot at top xy plane (not yet)
* DRAW_SIDES draw sides
*
* or any bitwise combination, e.g. "MAG_COLOR | DRAW_LINEX"
*
\*--------------------------------------------------------------------------*/
void
c_plmeshc(PLFLT *x, PLFLT *y, PLFLT **z, PLINT nx, PLINT ny, PLINT opt,
PLFLT *clevel, PLINT nlevel)
{
c_plot3dc(x, y, z, nx, ny, opt | MESH, clevel, nlevel);
}
/* clipping helper for 3d polygons */
int
plP_clip_poly(int Ni, PLFLT *Vi[3], int axis, PLFLT dir, PLFLT offset)
{
int anyout = 0;
PLFLT in[PL_MAXPOLY], T[3][PL_MAXPOLY];
int No = 0;
int i, j, k;
for(i=0; i<Ni; i++) {
in[i] = Vi[axis][i] * dir + offset;
anyout += in[i] < 0;
}
/* none out */
if(anyout == 0)
return Ni;
/* all out */
if(anyout == Ni) {
return 0;
}
/* some out */
/* copy over to a temp vector */
for(i=0; i<3; i++) {
for(j=0; j<Ni; j++) {
T[i][j] = Vi[i][j];
}
}
/* copy back selectively */
for(i=0; i<Ni; i++) {
j = (i+1) % Ni;
if(in[i]>=0 && in[j]>=0) {
for(k=0; k<3; k++)
Vi[k][No] = T[k][j];
No++;
} else if(in[i]>=0 && in[j]<0) {
PLFLT u = in[i] / (in[i] - in[j]);
for(k = 0; k<3; k++)
Vi[k][No] = T[k][i]*(1-u) + T[k][j]*u;
No++;
} else if(in[i]<0 && in[j]>=0) {
PLFLT u = in[i] / (in[i] - in[j]);
for(k = 0; k<3; k++)
Vi[k][No] = T[k][i]*(1-u) + T[k][j]*u;
No++;
for(k=0; k<3; k++)
Vi[k][No] = T[k][j];
No++;
}
}
return No;
}
/* helper for plsurf3d, similar to c_plfill3() */
static void
shade_triangle(PLFLT x0, PLFLT y0, PLFLT z0,
PLFLT x1, PLFLT y1, PLFLT z1,
PLFLT x2, PLFLT y2, PLFLT z2)
{
int i;
/* arrays for interface to core functions */
short u[6], v[6];
PLFLT x[6], y[6], z[6];
int n;
PLFLT xmin, xmax, ymin, ymax, zmin, zmax, zscale;
PLFLT *V[3];
plP_gdom(&xmin, &xmax, &ymin, &ymax);
plP_grange(&zscale, &zmin, &zmax);
x[0] = x0; x[1] = x1; x[2] = x2;
y[0] = y0; y[1] = y1; y[2] = y2;
z[0] = z0; z[1] = z1; z[2] = z2;
n = 3;
V[0] = x; V[1] = y; V[2] = z;
n = plP_clip_poly(n, V, 0, 1, -xmin);
n = plP_clip_poly(n, V, 0, -1, xmax);
n = plP_clip_poly(n, V, 1, 1, -ymin);
n = plP_clip_poly(n, V, 1, -1, ymax);
n = plP_clip_poly(n, V, 2, 1, -zmin);
n = plP_clip_poly(n, V, 2, -1, zmax);
if(n > 0) {
if (falsecolor)
plcol1(((z[0] + z[1] + z[2]) /3. - fc_minz) / (fc_maxz - fc_minz));
else
plcol1(plGetAngleToLight(x, y, z));
for(i=0; i<n; i++) {
u[i] = plP_wcpcx(plP_w3wcx(x[i], y[i], z[i]));
v[i] = plP_wcpcy(plP_w3wcy(x[i], y[i], z[i]));
}
u[n] = u[0];
v[n] = v[0];
#ifdef SHADE_DEBUG /* show triangles */
plcol0(1);
x[3]=x[0]; y[3]=y[0]; z[3]=z[0];
plline3(4,x,y,z);
#else /* fill triangles */
plP_fill(u, v, n+1);
#endif
}
}
/*--------------------------------------------------------------------------*\
* void plsurf3d(x, y, z, nx, ny, opt, clevel, nlevel)
*
* Plots the 3-d surface representation of the function z[x][y].
* The x values are stored as x[0..nx-1], the y values as y[0..ny-1],
* and the z values are in the 2-d array z[][]. The integer "opt" specifies:
* see plsurf3dl() below.
\*--------------------------------------------------------------------------*/
void
c_plsurf3d(PLFLT *x, PLFLT *y, PLFLT **z, PLINT nx, PLINT ny,
PLINT opt, PLFLT *clevel, PLINT nlevel)
{
PLINT i;
PLINT *indexymin = (PLINT *) malloc((size_t) (nx*sizeof(PLINT)));
PLINT *indexymax = (PLINT *) malloc((size_t) (nx*sizeof(PLINT)));
if ( ! indexymin || ! indexymax)
plexit("plsurf3d: Out of memory.");
for (i = 0; i < nx; i++) {
indexymin[i] = 0;
indexymax[i] = ny;
}
c_plsurf3dl(x, y, z, nx, ny, opt, clevel, nlevel,
0, nx, indexymin, indexymax);
free_mem(indexymin);
free_mem(indexymax);
}
/*--------------------------------------------------------------------------*\
* void plsurf3dl(x, y, z, nx, ny, opt, clevel, nlevel,
* ixstart, ixn, indexymin, indexymax)
*
* Plots the 3-d surface representation of the function z[x][y].
* The x values are stored as x[0..nx-1], the y values as y[0..ny-1],
* and the z values are in the 2-d array z[][].
*
*
* BASE_CONT draw contour plot at bottom xy plane
* TOP_CONT draw contour plot at top xy plane (not implemented)
* SURF_CONT draw contour at surface
* FACETED each square that makes up the surface is faceted
* DRAW_SIDES draw sides
* MAG_COLOR the surface is colored according to the value of z;
* if not defined, the surface is colored according to the
* intensity of the reflected light in the surface from a light
* source whose position is set using c_pllightsource()
*
* or any bitwise combination, e.g. "MAG_COLOR | SURF_CONT | SURF_BASE"
*
* indexymin and indexymax are arrays which specify the y index range
* (following the convention that the upper range limit is one more than
* actual index limit) for an x index range of ixstart, ixn.
* This code is a complete departure from the approach taken in the old version
* of this routine. Formerly to code attempted to use the logic for the hidden
* line algorithm to draw the hidden surface. This was really hard. This code
* below uses a simple back to front (painters) algorithm. All the
* triangles are drawn.
*
* There are multitude of ways this code could be optimized. Given the
* problems with the old code, I tried to focus on clarity here.
\*--------------------------------------------------------------------------*/
void
c_plsurf3dl(PLFLT *x, PLFLT *y, PLFLT **z, PLINT nx, PLINT ny,
PLINT opt, PLFLT *clevel, PLINT nlevel,
PLINT ixstart, PLINT ixn, PLINT *indexymin, PLINT *indexymax)
{
PLFLT cxx, cxy, cyx, cyy, cyz;
PLINT i, j, k;
PLINT ixDir, ixOrigin, iyDir, iyOrigin, nFast, nSlow;
PLINT ixFast, ixSlow, iyFast, iySlow;
PLINT iFast, iSlow;
PLFLT xmin, xmax, ymin, ymax, zmin, zmax, zscale;
PLFLT xm, ym, zm;
PLINT ixmin=0, ixmax=nx, iymin=0, iymax=ny;
PLFLT xx[3],yy[3],zz[3];
PLFLT px[4], py[4], pz[4];
CONT_LEVEL *cont, *clev;
CONT_LINE *cline;
int ct, ix, iy, iftriangle;
PLINT color = plsc->icol0, width = plsc->width;
if (plsc->level < 3) {
myabort("plsurf3dl: Please set up window first");
return;
}
if (nx <= 0 || ny <= 0) {
myabort("plsurf3dl: Bad array dimensions.");
return;
}
/*
* Don't use the data z value to scale the color, use the z axis
* values set by plw3d()
*
* plMinMax2dGrid(z, nx, ny, &fc_maxz, &fc_minz);
*/
fc_minz = plsc->ranmi;
fc_maxz = plsc->ranma;
if (fc_maxz == fc_minz) {
plwarn("plsurf3dl: Maximum and minimum Z values are equal! \"fixing\"...");
fc_maxz = fc_minz + 1e-6;
}
if (opt & MAG_COLOR)
falsecolor = 1;
else
falsecolor = 0;
plP_gdom(&xmin, &xmax, &ymin, &ymax);
plP_grange(&zscale, &zmin, &zmax);
if(zmin > zmax) {
PLFLT t = zmin;
zmin = zmax;
zmax = t;
}
/* Check that points in x and in y are strictly increasing and in range */
for (i = 0; i < nx - 1; i++) {
if (x[i] >= x[i + 1]) {
myabort("plsurf3dl: X array must be strictly increasing");
return;
}
if (x[i] < xmin && x[i+1] >= xmin)
ixmin = i;
if (x[i+1] > xmax && x[i] <= xmax)
ixmax = i+2;
}
for (i = 0; i < ny - 1; i++) {
if (y[i] >= y[i + 1]) {
myabort("plsurf3dl: Y array must be strictly increasing");
return;
}
if (y[i] < ymin && y[i+1] >= ymin)
iymin = i;
if (y[i+1] > ymax && y[i] <= ymax)
iymax = i+2;
}
/* get the viewing parameters */
plP_gw3wc(&cxx, &cxy, &cyx, &cyy, &cyz);
/* we're going to draw from back to front */
/* iFast will index the dominant (fastest changing) dimension
iSlow will index the slower changing dimension
iX indexes the X dimension
iY indexes the Y dimension */
/* get direction for X */
if(cxy >= 0) {
ixDir = 1; /* direction in X */
ixOrigin = ixmin; /* starting point */
} else {
ixDir = -1;
ixOrigin = ixmax - 1;
}
/* get direction for Y */
if(cxx >= 0) {
iyDir = -1;
iyOrigin = iymax - 1;
} else {
iyDir = 1;
iyOrigin = iymin;
}
/* figure out which dimension is dominant */
if (fabs(cxx) > fabs(cxy)) {
/* X is dominant */
nFast = ixmax - ixmin; /* samples in the Fast direction */
nSlow = iymax - iymin; /* samples in the Slow direction */
ixFast = ixDir; ixSlow = 0;
iyFast = 0; iySlow = iyDir;
} else {
nFast = iymax - iymin;
nSlow = ixmax - ixmin;
ixFast = 0; ixSlow = ixDir;
iyFast = iyDir; iySlow = 0;
}
/* we've got to draw the background grid first, hidden line code has to draw it last */
if (zbflg) {
PLFLT bx[3], by[3], bz[3];
PLFLT tick=zbtck, tp;
PLINT nsub=0;
/* get the tick spacing */
pldtik(zmin, zmax, &tick, &nsub);
/* determine the vertices for the background grid line */
bx[0] = (ixOrigin!=ixmin && ixFast==0) || ixFast > 0 ? xmax : xmin;
by[0] = (iyOrigin!=iymin && iyFast==0) || iyFast > 0 ? ymax : ymin;
bx[1] = ixOrigin!=ixmin ? xmax : xmin;
by[1] = iyOrigin!=iymin ? ymax : ymin;
bx[2] = (ixOrigin!=ixmin && ixSlow==0) || ixSlow > 0 ? xmax : xmin;
by[2] = (iyOrigin!=iymin && iySlow==0) || iySlow > 0 ? ymax : ymin;
plwid(zbwidth);
plcol0(zbcol);
for(tp = tick * floor(zmin / tick) + tick; tp <= zmax; tp += tick) {
bz[0] = bz[1] = bz[2] = tp;
plline3(3, bx, by, bz);
}
/* draw the vertical line at the back corner */
bx[0] = bx[1];
by[0] = by[1];
bz[0] = zmin;
plline3(2, bx, by, bz);
plwid(width);
plcol0(color);
}
/* If enabled, draw the contour at the base */
/* The contour plotted at the base will be identical to the one obtained
* with c_plcont(). The contour plotted at the surface is simple minded, but
* can be improved by using the contour data available.
*/
if (clevel != NULL && opt & BASE_CONT) {
#define NPTS 100
int np = NPTS;
PLFLT *zz = (PLFLT *) malloc(NPTS*sizeof(PLFLT));
/* get the contour lines */
/* prepare cont_store input */
PLFLT **zstore;
PLcGrid2 cgrid2;
cgrid2.nx = nx;
cgrid2.ny = ny;
plAlloc2dGrid(&cgrid2.xg, nx, ny);
plAlloc2dGrid(&cgrid2.yg, nx, ny);
plAlloc2dGrid(&zstore, nx, ny);
for (i = ixstart; i < ixn; i++) {
for (j = 0; j < indexymin[i]; j++) {
cgrid2.xg[i][j] = x[i];
cgrid2.yg[i][j] = y[indexymin[i]];
zstore[i][j] = z[i][indexymin[i]];
}
for (j = indexymin[i]; j < indexymax[i]; j++) {
cgrid2.xg[i][j] = x[i];
cgrid2.yg[i][j] = y[j];
zstore[i][j] = z[i][j];
}
for (j = indexymax[i]; j < ny; j++) {
cgrid2.xg[i][j] = x[i];
cgrid2.yg[i][j] = y[indexymax[i]-1];
zstore[i][j] = z[i][indexymax[i]-1];
}
}
/* Fill cont structure with contours. */
cont_store(zstore, nx, ny, ixstart+1, ixn, 1, ny,
clevel, nlevel, pltr2, (void *) &cgrid2, &cont);
/* Free the 2D input arrays to cont_store since not needed any more. */
plFree2dGrid(zstore, nx, ny);
plFree2dGrid(cgrid2.xg, nx, ny);
plFree2dGrid(cgrid2.yg, nx, ny);
/* follow the contour levels and lines */
clev = cont;
do { /* for each contour level */
cline = clev->line;
do { /* there are several lines that make up the contour */
if (cline->npts > np) {
np = cline->npts;
zz = (PLFLT *) realloc(zz, np*sizeof(PLFLT));
}
for (j=0; j<cline->npts; j++)
zz[j] = plsc->ranmi;
if (cline->npts > 0) {
plcol1((clev->level-fc_minz)/(fc_maxz-fc_minz));
plline3(cline->npts, cline->x, cline->y, zz);
}
cline = cline->next;
}
while(cline != NULL);
clev = clev->next;
}
while(clev != NULL);
cont_clean_store(cont); /* now release the memory */
free(zz);
}
/* Now we can iterate over the grid drawing the quads */
for(iSlow=0; iSlow < nSlow-1; iSlow++) {
for(iFast=0; iFast < nFast-1; iFast++) {
/* get the 4 corners of the Quad, which are
*
* 0--2
* | |
* 1--3
*/
xm = ym = zm = 0.;
iftriangle = 1;
for(i=0; i<2; i++) {
for(j=0; j<2; j++) {
/* we're transforming from Fast/Slow coordinates to x/y coordinates */
/* note, these are the indices, not the values */
ix = ixFast * (iFast+i) + ixSlow * (iSlow+j) + ixOrigin;
iy = iyFast * (iFast+i) + iySlow * (iSlow+j) + iyOrigin;
if (ixstart <= ix && ix < ixn &&
indexymin[ix] <= iy && iy < indexymax[ix]) {
xm += px[2*i+j] = x[ix];
ym += py[2*i+j] = y[iy];
zm += pz[2*i+j] = z[ix][iy];
}
else {
iftriangle = 0;
break;
}
}
if (iftriangle == 0)
break;
}
if (iftriangle == 0)
continue;
/* the "mean point" of the quad, common to all four triangles
-- perhaps not a good thing to do for the light shading */
xm /= 4.; ym /= 4.; zm /= 4.;
/* now draw the quad as four triangles */
for (i=1; i<3; i++) {
for (j=0; j<4; j+=3) {
shade_triangle(px[j], py[j], pz[j], xm, ym, zm, px[i], py[i], pz[i]);
/* after shading, see if the triangle crosses one contour plane */
#define min3(a,b,c) (MIN((MIN(a,b)),c))
#define max3(a,b,c) (MAX((MAX(a,b)),c))
if (clevel != NULL && (opt & SURF_CONT)) {
for (k=0; k<nlevel; k++) {
if (clevel[k] >= min3(pz[i],zm,pz[j]) && clevel[k] < max3(pz[i],zm,pz[j])) {
ct = 0;
if (clevel[k] >= MIN(pz[i],zm) && clevel[k] < MAX(pz[i],zm)) { /* p0-pm */
xx[ct] = ((clevel[k] - pz[i])*(xm-px[i]))/(zm-pz[i]) + px[i];
yy[ct] = ((clevel[k] - pz[i])*(ym-py[i]))/(zm-pz[i]) + py[i];
ct++;
}
if (clevel[k] >= MIN(pz[i],pz[j]) && clevel[k] < MAX(pz[i],pz[j])) { /* p0-p1 */
xx[ct] = ((clevel[k] - pz[i])*(px[j]-px[i]))/(pz[j]-pz[i]) + px[i];
yy[ct] = ((clevel[k] - pz[i])*(py[j]-py[i]))/(pz[j]-pz[i]) + py[i];
ct++;
}
if (clevel[k] >= MIN(pz[j],zm) && clevel[k] < MAX(pz[j],zm)) { /* p1-pm */
xx[ct] = ((clevel[k] - pz[j])*(xm-px[j]))/(zm-pz[j]) + px[j];
yy[ct] = ((clevel[k] - pz[j])*(ym-py[j]))/(zm-pz[j]) + py[j];
ct++;
}
if (ct == 2) {
/* yes, xx and yy are the intersection points of the triangle with
* the contour line -- draw a straight line betweeen the points
* -- at the end this will make up the contour line */
if (opt & SURF_CONT) {
/* surface contour with color set by user */
plcol0(color);
zz[0] = zz[1] = clevel[k];
plline3(2, xx, yy, zz);
}
/* don't break; one triangle can span various contour levels */
} else
plwarn("plsurf3dl: ***ERROR***\n");
}
}
}
}
}
}
}
if (opt & FACETED) {
plcol0(0);
c_plot3dcl(x, y, z, nx, ny, MESH | DRAW_LINEXY, NULL, 0,
ixstart, ixn, indexymin, indexymax);
}
if (opt & DRAW_SIDES) { /* the sides look ugly !!! */
/* draw one more row with all the Z's set to zmin */
PLFLT zscale, zmin, zmax;
plP_grange(&zscale, &zmin, &zmax);
iSlow = nSlow-1;
iftriangle = 1;
for(iFast=0; iFast < nFast-1; iFast++) {
for(i=0; i<2; i++) {
ix = ixFast * (iFast+i) + ixSlow * iSlow + ixOrigin;
iy = iyFast * (iFast+i) + iySlow * iSlow + iyOrigin;
if (ixstart <= ix && ix < ixn &&
indexymin[ix] <= iy && iy < indexymax[ix]) {
px[2*i] = x[ix];
py[2*i] = y[iy];
pz[2*i] = z[ix][iy];
}
else {
iftriangle = 0;
break;
}
}
if (iftriangle == 0)
break;
/* now draw the quad as two triangles (4 might be better) */
shade_triangle(px[0], py[0], pz[0], px[2], py[2], pz[2], px[0], py[0], zmin);
shade_triangle(px[2], py[2], pz[2], px[2], py[2], zmin, px[0], py[0], zmin);
}
iFast = nFast-1;
iftriangle = 1;
for(iSlow=0; iSlow < nSlow-1; iSlow++) {
for(i=0; i<2; i++) {
ix = ixFast * iFast + ixSlow * (iSlow+i) + ixOrigin;
iy = iyFast * iFast + iySlow * (iSlow+i) + iyOrigin;
if (ixstart <= ix && ix < ixn &&
indexymin[ix] <= iy && iy < indexymax[ix]) {
px[2*i] = x[ix];
py[2*i] = y[iy];
pz[2*i] = z[ix][iy];
}
else {
iftriangle = 0;
break;
}
}
if (iftriangle == 0)
break;
/* now draw the quad as two triangles (4 might be better) */
shade_triangle(px[0], py[0], pz[0], px[2], py[2], pz[2], px[0], py[0], zmin);
shade_triangle(px[2], py[2], pz[2], px[2], py[2], zmin, px[0], py[0], zmin);
}
}
}
/*--------------------------------------------------------------------------*\
* void plot3d(x, y, z, nx, ny, opt, side)
*
* Plots a 3-d representation of the function z[x][y]. The x values
* are stored as x[0..nx-1], the y values as y[0..ny-1], and the z
* values are in the 2-d array z[][]. The integer "opt" specifies:
* see plot3dcl() below
\*--------------------------------------------------------------------------*/
void
c_plot3d(PLFLT *x, PLFLT *y, PLFLT **z,
PLINT nx, PLINT ny, PLINT opt, PLBOOL side)
{
c_plot3dc( x, y, z, nx, ny, opt | (side != 0 ? DRAW_SIDES : 0), NULL, 0);
}
/*--------------------------------------------------------------------------*\
* void plot3dc(x, y, z, nx, ny, opt, clevel, nlevel)
*
* Plots a 3-d representation of the function z[x][y]. The x values
* are stored as x[0..nx-1], the y values as y[0..ny-1], and the z
* values are in the 2-d array z[][]. The integer "opt" specifies:
* see plot3dcl() below
\*--------------------------------------------------------------------------*/
void
c_plot3dc(PLFLT *x, PLFLT *y, PLFLT **z,
PLINT nx, PLINT ny, PLINT opt,
PLFLT *clevel, PLINT nlevel)
{
c_plot3dcl( x, y, z, nx, ny, opt, clevel, nlevel, 0, 0, NULL, NULL);
}
/*--------------------------------------------------------------------------*\
* void plot3dcl(x, y, z, nx, ny, opt, clevel, nlevel,
* ixstart, ixn, indexymin, indexymax)
*
* Plots a 3-d representation of the function z[x][y]. The x values
* are stored as x[0..nx-1], the y values as y[0..ny-1], and the z
* values are in the 2-d array z[][]. The integer "opt" specifies:
*
* DRAW_LINEX : Draw lines parallel to x-axis
* DRAW_LINEY : Draw lines parallel to y-axis
* DRAW_LINEXY: Draw lines parallel to both axes
* MAG_COLOR: Magnitude coloring of wire frame
* BASE_CONT: Draw contour at bottom xy plane
* TOP_CONT: Draw contour at top xy plane (not yet)
* DRAW_SIDES: Draw sides around the plot
* MESH: Draw the "under" side of the plot
*
* or any bitwise combination, e.g. "MAG_COLOR | DRAW_LINEX"
* indexymin and indexymax are arrays which specify the y index limits
* (following the convention that the upper range limit is one more than
* actual index limit) for an x index range of ixstart, ixn.
\*--------------------------------------------------------------------------*/
void
c_plot3dcl(PLFLT *x, PLFLT *y, PLFLT **z,
PLINT nx, PLINT ny, PLINT opt,
PLFLT *clevel, PLINT nlevel,
PLINT ixstart, PLINT ixn, PLINT *indexymin, PLINT *indexymax)
{
PLFLT cxx, cxy, cyx, cyy, cyz;
PLINT init, i, ix, iy, color, width;
PLFLT xmin, xmax, ymin, ymax, zmin, zmax, zscale;
PLINT ixmin=0, ixmax=nx-1, iymin=0, iymax=ny-1;
PLINT clipped = 0, base_cont = 0, side = 0;
pl3mode = 0;
if (plsc->level < 3) {
myabort("plot3dcl: Please set up window first");
return;
}
if (opt < 1) {
myabort("plot3dcl: Bad option");
return;
}
if (nx <= 0 || ny <= 0) {
myabort("plot3dcl: Bad array dimensions.");
return;
}
plP_gdom(&xmin, &xmax, &ymin, &ymax);
plP_grange(&zscale, &zmin, &zmax);
if(zmin > zmax) {
PLFLT t = zmin;
zmin = zmax;
zmax = t;
}
/* Check that points in x and in y are strictly increasing */
for (i = 0; i < nx - 1; i++) {
if (x[i] >= x[i + 1]) {
myabort("plot3dcl: X array must be strictly increasing");
return;
}
}
for (i = 0; i < ny - 1; i++) {
if (y[i] >= y[i + 1]) {
myabort("plot3dcl: Y array must be strictly increasing");
return;
}
}
if (opt & MESH)
pl3mode = 1;
if (opt & DRAW_SIDES)
side = 1;
/* figure out the part of the data to use */
if (xmin < x[0])
xmin = x[0];
if (xmax > x[nx-1])
xmax = x[nx-1];
if (ymin < y[0])
ymin = y[0];
if (ymax > y[ny-1])
ymax = y[ny-1];
for (ixmin = 0; ixmin < nx-1 && x[ixmin+1] < xmin; ixmin++) {}
for (ixmax = nx-1; ixmax > 0 && x[ixmax-1] > xmax; ixmax--) {}
for (iymin = 0; iymin < ny-1 && y[iymin+1] < ymin; iymin++) {}
for (iymax = ny-1; iymax > 0 && y[iymax-1] > ymax; iymax--) {}
/*fprintf(stderr, "(%d,%d) %d %d %d %d\n", nx, ny, ixmin, ixmax, iymin, iymax);*/
/* do we need to clip? */
if(ixmin > 0 || ixmax < nx-1 || iymin > 0 || iymax < ny-1) {
/* adjust the input so it stays within bounds */
int _nx = ixmax - ixmin + 1;
int _ny = iymax - iymin + 1;
PLFLT *_x, *_y, **_z;
PLFLT ty0, ty1, tx0, tx1;
int i, j;
if(_nx <= 1 || _ny <= 1) {
myabort("plot3dcl: selected x or y range has no data");
return;
}
/* allocate storage for new versions of the input vectors */
_x = (PLFLT*)malloc(_nx * sizeof(PLFLT));
_y = (PLFLT*)malloc(_ny * sizeof(PLFLT));
_z = (PLFLT**)malloc(_nx * sizeof(PLFLT*));
clipped = 1;
/* copy over the independent variables */
_x[0] = xmin;
_x[_nx-1] = xmax;
for(i=1; i<_nx-1; i++)
_x[i] = x[ixmin+i];
_y[0] = ymin;
_y[_ny-1] = ymax;
for(i=1; i<_ny-1; i++)
_y[i] = y[iymin+i];
/* copy the data array so we can interpolate around the edges */
for(i=0; i<_nx; i++)
_z[i] = (PLFLT*)malloc(_ny * sizeof(PLFLT));
/* interpolation factors for the 4 edges */
ty0 = (_y[0] - y[iymin]) / (y[iymin+1] - y[iymin]);
ty1 = (_y[_ny-1] - y[iymax-1]) / (y[iymax] - y[iymax-1]);
tx0 = (_x[0] - x[ixmin]) / (x[ixmin+1] - x[ixmin]);
tx1 = (_x[_nx-1] - x[ixmax-1]) / (x[ixmax] - x[ixmax-1]);
for(i=0; i<_nx; i++) {
if(i==0) {
_z[i][0] = z[ixmin][iymin]*(1-ty0)*(1-tx0) + z[ixmin][iymin+1]*(1-tx0)*ty0
+ z[ixmin+1][iymin]*tx0*(1-ty0) + z[ixmin+1][iymin+1]*tx0*ty0;
for(j=1; j<_ny-1; j++)
_z[i][j] = z[ixmin][iymin+j]*(1-tx0) + z[ixmin+1][iymin+j]*tx0;
_z[i][_ny-1] = z[ixmin][iymax-1]*(1-tx0)*(1-ty1) + z[ixmin][iymax]*(1-tx0)*ty1
+ z[ixmin+1][iymax-1]*tx0*(1-ty1) + z[ixmin+1][iymax]*tx0*ty1;
} else if(i==_nx-1) {
_z[i][0] = z[ixmax-1][iymin]*(1-tx1)*(1-ty0) + z[ixmax-1][iymin+1]*(1-tx1)*ty0
+ z[ixmax][iymin]*tx1*(1-ty0) + z[ixmax][iymin+1]*tx1*ty0;
for(j=1; j<_ny-1; j++)
_z[i][j] = z[ixmax-1][iymin+j]*(1-tx1) + z[ixmax][iymin+j]*tx1;
_z[i][_ny-1] = z[ixmax-1][iymax-1]*(1-tx1)*(1-ty1) + z[ixmax][iymax]*(1-tx1)*ty1
+ z[ixmax][iymax-1]*tx1*(1-ty1) + z[ixmax][iymax]*tx1*ty1;
} else {
_z[i][0] = z[ixmin+i][iymin]*(1-ty0) + z[ixmin+i][iymin+1]*ty0;
for(j=1; j<_ny-1; j++)
_z[i][j] = z[ixmin+i][iymin+j];
_z[i][_ny-1] = z[ixmin+i][iymax-1]*(1-ty1) + z[ixmin+i][iymax]*ty1;
}
for(j=0; j<_ny; j++) {
if(_z[i][j] < zmin)
_z[i][j] = zmin;
else if(_z[i][j] > zmax)
_z[i][j] = zmax;
}
}
/* replace the input with our clipped versions */
x = _x;
y = _y;
z = _z;
nx = _nx;
ny = _ny;
}
if ((opt & BASE_CONT) || (opt & TOP_CONT) || (opt && MAG_COLOR)) {
/*
* Don't use the data z value to scale the color, use the z axis
* values set by plw3d()
*
* plMinMax2dGrid(z, nx, ny, &fc_maxz, &fc_minz);
*/
fc_minz = plsc->ranmi;
fc_maxz = plsc->ranma;
if (fc_maxz == fc_minz) {
plwarn("plot3dcl: Maximum and minimum Z values are equal! \"fixing\"...");
fc_maxz = fc_minz + 1e-6;
}
}
if (opt & BASE_CONT) { /* If enabled, draw the contour at the base. */
if (clevel != NULL && nlevel != 0) {
base_cont = 1;
/* even if MESH is not set, "set it",
as the base contour can only be done in this case */
pl3mode = 1;
}
}
if (opt & MAG_COLOR) { /* If enabled, use magnitude colored wireframe */
ctmp = (PLFLT *) malloc((size_t) (2 * MAX(nx, ny) * sizeof(PLFLT)));
} else
ctmp = NULL;
/* next logic only knows opt = 1 | 2 | 3, make sure that it only gets that */
opt &= DRAW_LINEXY;
/* Allocate work arrays */
utmp = (PLINT *) malloc((size_t) (2 * MAX(nx, ny) * sizeof(PLINT)));
vtmp = (PLINT *) malloc((size_t) (2 * MAX(nx, ny) * sizeof(PLINT)));
if ( ! utmp || ! vtmp)
myexit("plot3dcl: Out of memory.");
plP_gw3wc(&cxx, &cxy, &cyx, &cyy, &cyz);
init = 1;
/* Call 3d line plotter. Each viewing quadrant
(perpendicular to x-y plane) must be handled separately. */
if (cxx >= 0.0 && cxy <= 0.0) {
if (opt == DRAW_LINEY)
plt3zz(1, ny, 1, -1, -opt, &init, x, y, z, nx, ny, utmp, vtmp,ctmp);
else {
for (iy = 2; iy <= ny; iy++)
plt3zz(1, iy, 1, -1, -opt, &init, x, y, z, nx, ny, utmp, vtmp,ctmp);
}
if (opt == DRAW_LINEX)
plt3zz(1, ny, 1, -1, opt, &init, x, y, z, nx, ny, utmp, vtmp, ctmp);
else {
for (ix = 1; ix <= nx - 1; ix++)
plt3zz(ix, ny, 1, -1, opt, &init, x, y, z, nx, ny, utmp, vtmp, ctmp);
}
}
else if (cxx <= 0.0 && cxy <= 0.0) {
if (opt == DRAW_LINEX)
plt3zz(nx, ny, -1, -1, opt, &init, x, y, z, nx, ny, utmp, vtmp, ctmp);
else {
for (ix = 2; ix <= nx; ix++)
plt3zz(ix, ny, -1, -1, opt, &init, x, y, z, nx, ny, utmp, vtmp, ctmp);
}
if (opt == DRAW_LINEY)
plt3zz(nx, ny, -1, -1, -opt, &init, x, y, z, nx, ny, utmp, vtmp, ctmp);
else {
for (iy = ny; iy >= 2; iy--)
plt3zz(nx, iy, -1, -1, -opt, &init, x, y, z, nx, ny, utmp, vtmp, ctmp);
}
}
else if (cxx <= 0.0 && cxy >= 0.0) {
if (opt == DRAW_LINEY)
plt3zz(nx, 1, -1, 1, -opt, &init, x, y, z, nx, ny, utmp, vtmp, ctmp);
else {
for (iy = ny - 1; iy >= 1; iy--)
plt3zz(nx, iy, -1, 1, -opt, &init, x, y, z, nx, ny, utmp, vtmp, ctmp);
}
if (opt == DRAW_LINEX)
plt3zz(nx, 1, -1, 1, opt, &init, x, y, z, nx, ny, utmp, vtmp, ctmp);
else {
for (ix = nx; ix >= 2; ix--)
plt3zz(ix, 1, -1, 1, opt, &init, x, y, z, nx, ny, utmp, vtmp, ctmp);
}
}
else if (cxx >= 0.0 && cxy >= 0.0) {
if (opt == DRAW_LINEX)
plt3zz(1, 1, 1, 1, opt, &init, x, y, z, nx, ny, utmp, vtmp, ctmp);
else {
for (ix = nx - 1; ix >= 1; ix--)
plt3zz(ix, 1, 1, 1, opt, &init, x, y, z, nx, ny, utmp, vtmp, ctmp);
}
if (opt == DRAW_LINEY)
plt3zz(1, 1, 1, 1, -opt, &init, x, y, z, nx, ny, utmp, vtmp, ctmp);
else {
for (iy = 1; iy <= ny - 1; iy++)
plt3zz(1, iy, 1, 1, -opt, &init, x, y, z, nx, ny, utmp, vtmp, ctmp);
}
}
/* draw contour at the base. Not 100%! Why? */
if (base_cont){
int np = NPTS, j;
CONT_LEVEL *cont, *clev;
CONT_LINE *cline;
PLINT *uu = (PLINT *) malloc(NPTS*sizeof(PLINT));
PLINT *vv = (PLINT *) malloc(NPTS*sizeof(PLINT));
/* prepare cont_store input */
PLFLT **zstore;
PLcGrid2 cgrid2;
cgrid2.nx = nx;
cgrid2.ny = ny;
plAlloc2dGrid(&cgrid2.xg, nx, ny);
plAlloc2dGrid(&cgrid2.yg, nx, ny);
plAlloc2dGrid(&zstore, nx, ny);
for (i = 0; i < nx; i++) {
for (j = 0; j < ny; j++) {
cgrid2.xg[i][j] = x[i];
cgrid2.yg[i][j] = y[j];
zstore[i][j] = z[i][j];
}
}
pl3upv = 0;
/* Fill cont structure with contours. */
cont_store(zstore, nx, ny, 1, nx, 1, ny,
clevel, nlevel, pltr2, (void *) &cgrid2, &cont);
/* Free the 2D input arrays to cont_store since not needed any more. */
plFree2dGrid(zstore, nx, ny);
plFree2dGrid(cgrid2.xg, nx, ny);
plFree2dGrid(cgrid2.yg, nx, ny);
/* follow the contour levels and lines */
clev = cont;
do { /* for each contour level */
cline = clev->line;
do { /* there are several lines that make up each contour */
int cx, i, k, l, m, start, end;
PLFLT tx, ty;
if (cline->npts > np) {
np = cline->npts;
uu = (PLINT *) realloc(uu, np*sizeof(PLINT));
vv = (PLINT *) realloc(vv, np*sizeof(PLINT));
}
/* the hidden line plotter plnxtv() only works OK if the x points are in increasing order.
As this does not always happens, the situation must be detected and the line segment
must be reversed before being plotted */
i = 0;
if (cline->npts > 0) {
do {
plcol1((clev->level-fc_minz)/(fc_maxz-fc_minz));
cx = plP_wcpcx(plP_w3wcx(cline->x[i],cline->y[i], plsc->ranmi));
for (j=i; j < cline->npts; j++) { /* convert to 2D coordinates */
uu[j] = plP_wcpcx(plP_w3wcx(cline->x[j],cline->y[j], plsc->ranmi));
vv[j] = plP_wcpcy(plP_w3wcy(cline->x[j],cline->y[j], plsc->ranmi));
if (uu[j] < cx) /* find turn back point */
break;
else
cx = uu[j];
}
plnxtv(&uu[i], &vv[i], NULL, j-i, 0); /* plot line with increasing x */
if (j < cline->npts) { /* line not yet finished, */
start = j-1;
for (i = start; i < cline->npts; i++) { /* search turn forward point */
uu[i] = plP_wcpcx(plP_w3wcx(cline->x[i],cline->y[i], plsc->ranmi));
if (uu[i] > cx)
break;
else
cx = uu[i];
}
end = i-1;
for (k=0; k < (end-start+1)/2; k++) { /* reverse line segment */
l = start+k;
m = end-k;
tx = cline->x[l];
ty = cline->y[l];
cline->x[l] = cline->x[m];
cline->y[l] = cline->y[m];
cline->x[m] = tx;
cline->y[m] = ty;
}
/* convert to 2D coordinates */
for (j=start; j <= end; j++) {
uu[j] = plP_wcpcx(plP_w3wcx(cline->x[j],cline->y[j], plsc->ranmi));
vv[j] = plP_wcpcy(plP_w3wcy(cline->x[j],cline->y[j], plsc->ranmi));
}
plnxtv(&uu[start], &vv[start], NULL, end-start+1, 0); /* and plot it */
cline->x[end] = cline->x[start];
cline->y[end] = cline->y[start];
i = end; /* restart where it was left */
}
} while (j < cline->npts && i < cline->npts);
}
cline = cline->next;
}
while(cline != NULL);
clev = clev->next;
}
while(clev != NULL);
cont_clean_store(cont); /* now release the contour memory */
pl3upv = 1;
free(uu);
free(vv);
}
/* Finish up by drawing sides, background grid (both are optional) */
if (side)
plside3(x, y, z, nx, ny, opt);
if (zbflg) {
color = plsc->icol0;
width = plsc->width;
plwid(zbwidth);
plcol0(zbcol);
plgrid3(zbtck);
plwid(width);
plcol0(color);
}
freework();
if(clipped) {
free(x);
free(y);
for(i=0; i<nx; i++)
free(z[i]);
free(z);
}
}
/*--------------------------------------------------------------------------*\
* void plxyindexlimits()
*
* Transform from y array limits to corresponding x array limits (or vice
* versa).
*
* N.B. we follow the convention here that all upper range limits are one
* more than the actual last index.
* instart (>= 0) through inn is the index range where
* the input inarray_min and inarray_max arrays are defined.
*
* outstart (>= 0), through outn (with outn <= outnmax) is the index range
* where the output outarray_min and outarray_max arrays are defined.
*
* In order to assure the transformation from y array limits to x array limits
* (or vice versa) is single-valued, this programme plaborts if the
* inarray_min array has a maximum or inarray_max array has a minimum.
\*--------------------------------------------------------------------------*/
static void
plxyindexlimits(PLINT instart, PLINT inn,
PLINT *inarray_min, PLINT *inarray_max,
PLINT *outstart, PLINT *outn, PLINT outnmax,
PLINT *outarray_min, PLINT *outarray_max)
{
PLINT i, j;
if (inn < 0) {
myabort("plxyindexlimits: Must have instart >= 0");
return;
}
if (inn - instart <= 0) {
myabort("plxyindexlimits: Must have at least 1 defined point");
return;
}
*outstart = inarray_min[instart];
*outn = inarray_max[instart];
for (i=instart; i<inn; i++) {
*outstart = MIN(*outstart, inarray_min[i]);
*outn = MAX(*outn, inarray_max[i]);
if (i + 2 < inn) {
if (inarray_min[i] < inarray_min[i+1] &&
inarray_min[i+1] > inarray_min[i+2]) {
myabort("plxyindexlimits: inarray_min must not have a maximum");
return;
}
if (inarray_max[i] > inarray_max[i+1] &&
inarray_max[i+1] < inarray_max[i+2]) {
myabort("plxyindexlimits: inarray_max must not have a minimum");
return;
}
}
}
if (*outstart < 0) {
myabort("plxyindexlimits: Must have all elements of inarray_min >= 0");
return;
}
if (*outn > outnmax) {
myabort("plxyindexlimits: Must have all elements of inarray_max <= outnmax");
return;
}
for (j=*outstart; j < *outn; j++) {
i = instart;
/* Find first valid i for this j. */
while (i < inn && !(inarray_min[i] <= j && j < inarray_max[i]))
i++;
if (i < inn)
outarray_min[j] = i;
else {
myabort("plxyindexlimits: bad logic; invalid i should never happen");
return;
}
/* Find next invalid i for this j. */
while (i < inn && (inarray_min[i] <= j && j < inarray_max[i]))
i++;
outarray_max[j] = i;
}
}
/*--------------------------------------------------------------------------*\
* void plP_gzback()
*
* Get background parameters for 3d plot.
\*--------------------------------------------------------------------------*/
void
plP_gzback(PLINT **zbf, PLINT **zbc, PLFLT **zbt, PLINT **zbw)
{
*zbf = &zbflg;
*zbc = &zbcol;
*zbt = &zbtck;
*zbw = &zbwidth;
}
/*--------------------------------------------------------------------------*\
* PLFLT plGetAngleToLight()
*
* Gets cos of angle between normal to a polygon and a light source.
* Requires at least 3 elements, forming non-parallel lines
* in the arrays.
\*--------------------------------------------------------------------------*/
static PLFLT
plGetAngleToLight(PLFLT* x, PLFLT* y, PLFLT* z)
{
PLFLT vx1, vx2, vy1, vy2, vz1, vz2;
PLFLT px, py, pz;
PLFLT vlx, vly, vlz;
PLFLT mag1, mag2;
PLFLT cosangle;
vx1 = x[1] - x[0];
vx2 = x[2] - x[1];
vy1 = y[1] - y[0];
vy2 = y[2] - y[1];
vz1 = z[1] - z[0];
vz2 = z[2] - z[1];
/* Find vector perpendicular to the face */
px = vy1*vz2 - vz1*vy2;
py = vz1*vx2 - vx1*vz2;
pz = vx1*vy2 - vy1*vx2;
mag1 = px*px + py*py + pz*pz;
/* Vectors were parallel! */
if (mag1 == 0)
return 1;
vlx = xlight - x[0];
vly = ylight - y[0];
vlz = zlight - z[0];
mag2 = vlx*vlx + vly*vly + vlz*vlz;
if (mag2 ==0)
return 1;
/* Now have 3 vectors going through the first point on the given surface */
cosangle = fabs( (vlx*px + vly*py + vlz*pz) / sqrt(mag1*mag2) );
/* In case of numerical rounding */
if (cosangle > 1) cosangle = 1;
return cosangle;
}
/*--------------------------------------------------------------------------*\
* void plt3zz()
*
* Draws the next zig-zag line for a 3-d plot. The data is stored in array
* z[][] as a function of x[] and y[], and is plotted out starting at index
* (x0,y0).
*
* Depending on the state of "flag", the sequence of data points sent to
* plnxtv is altered so as to allow cross-hatch plotting, or plotting
* parallel to either the x-axis or the y-axis.
\*--------------------------------------------------------------------------*/
static void
plt3zz(PLINT x0, PLINT y0, PLINT dx, PLINT dy, PLINT flag, PLINT *init,
PLFLT *x, PLFLT *y, PLFLT **z, PLINT nx, PLINT ny,
PLINT *u, PLINT *v, PLFLT* c)
{
PLINT n = 0;
PLFLT x2d, y2d;
while (1 <= x0 && x0 <= nx && 1 <= y0 && y0 <= ny) {
x2d = plP_w3wcx(x[x0 - 1], y[y0 - 1], z[x0 - 1][y0 - 1]);
y2d = plP_w3wcy(x[x0 - 1], y[y0 - 1], z[x0 - 1][y0 - 1]);
u[n] = plP_wcpcx(x2d);
v[n] = plP_wcpcy(y2d);
if (c != NULL)
c[n] = (z[x0 - 1][y0 - 1] - fc_minz)/(fc_maxz-fc_minz);
switch (flag) {
case -3:
y0 += dy;
flag = -flag;
break;
case -2:
y0 += dy;
break;
case -1:
y0 += dy;
flag = -flag;
break;
case 1:
x0 += dx;
break;
case 2:
x0 += dx;
flag = -flag;
break;
case 3:
x0 += dx;
flag = -flag;
break;
}
n++;
}
if (flag == 1 || flag == -2) {
if (flag == 1) {
x0 -= dx;
y0 += dy;
}
else if (flag == -2) {
y0 -= dy;
x0 += dx;
}
if (1 <= x0 && x0 <= nx && 1 <= y0 && y0 <= ny) {
x2d = plP_w3wcx( x[x0 - 1], y[y0 - 1], z[x0 - 1][y0 - 1]);
y2d = plP_w3wcy( x[x0 - 1], y[y0 - 1], z[x0 - 1][y0 - 1]);
u[n] = plP_wcpcx(x2d);
v[n] = plP_wcpcy(y2d);
if (c != NULL)
c[n] = (z[x0 - 1][y0 - 1] - fc_minz)/(fc_maxz-fc_minz);
n++;
}
}
/* All the setup is done. Time to do the work. */
plnxtv(u, v, c, n, *init);
*init = 0;
}
/*--------------------------------------------------------------------------*\
* void plside3()
*
* This routine draws sides around the front of the 3d plot so that
* it does not appear to float.
\*--------------------------------------------------------------------------*/
static void
plside3(PLFLT *x, PLFLT *y, PLFLT **z, PLINT nx, PLINT ny, PLINT opt)
{
PLINT i;
PLFLT cxx, cxy, cyx, cyy, cyz;
PLFLT xmin, ymin, zmin, xmax, ymax, zmax, zscale;
PLFLT tx, ty, ux, uy;
plP_gw3wc(&cxx, &cxy, &cyx, &cyy, &cyz);
plP_gdom(&xmin, &xmax, &ymin, &ymax);
plP_grange(&zscale, &zmin, &zmax);
/* Get x, y coordinates of legs and plot */
if (cxx >= 0.0 && cxy <= 0.0) {
if (opt != 1) {
for (i = 0; i < nx; i++) {
tx = plP_w3wcx(x[i], y[0], zmin);
ty = plP_w3wcy(x[i], y[0], zmin);
ux = plP_w3wcx(x[i], y[0], z[i][0]);
uy = plP_w3wcy(x[i], y[0], z[i][0]);
pljoin(tx, ty, ux, uy);
}
}
if (opt != 2) {
for (i = 0; i < ny; i++) {
tx = plP_w3wcx(x[0], y[i], zmin);
ty = plP_w3wcy(x[0], y[i], zmin);
ux = plP_w3wcx(x[0], y[i], z[0][i]);
uy = plP_w3wcy(x[0], y[i], z[0][i]);
pljoin(tx, ty, ux, uy);
}
}
}
else if (cxx <= 0.0 && cxy <= 0.0) {
if (opt != 1) {
for (i = 0; i < nx; i++) {
tx = plP_w3wcx(x[i], y[ny - 1], zmin);
ty = plP_w3wcy(x[i], y[ny - 1], zmin);
ux = plP_w3wcx(x[i], y[ny - 1], z[i][ny - 1]);
uy = plP_w3wcy(x[i], y[ny - 1], z[i][ny - 1]);
pljoin(tx, ty, ux, uy);
}
}
if (opt != 2) {
for (i = 0; i < ny; i++) {
tx = plP_w3wcx(x[0], y[i], zmin);
ty = plP_w3wcy(x[0], y[i], zmin);
ux = plP_w3wcx(x[0], y[i], z[0][i]);
uy = plP_w3wcy(x[0], y[i], z[0][i]);
pljoin(tx, ty, ux, uy);
}
}
}
else if (cxx <= 0.0 && cxy >= 0.0) {
if (opt != 1) {
for (i = 0; i < nx; i++) {
tx = plP_w3wcx(x[i], y[ny - 1], zmin);
ty = plP_w3wcy(x[i], y[ny - 1], zmin);
ux = plP_w3wcx(x[i], y[ny - 1], z[i][ny - 1]);
uy = plP_w3wcy(x[i], y[ny - 1], z[i][ny - 1]);
pljoin(tx, ty, ux, uy);
}
}
if (opt != 2) {
for (i = 0; i < ny; i++) {
tx = plP_w3wcx(x[nx - 1], y[i], zmin);
ty = plP_w3wcy(x[nx - 1], y[i], zmin);
ux = plP_w3wcx(x[nx - 1], y[i], z[nx - 1][i]);
uy = plP_w3wcy(x[nx - 1], y[i], z[nx - 1][i]);
pljoin(tx, ty, ux, uy);
}
}
}
else if (cxx >= 0.0 && cxy >= 0.0) {
if (opt != 1) {
for (i = 0; i < nx; i++) {
tx = plP_w3wcx(x[i], y[0], zmin);
ty = plP_w3wcy(x[i], y[0], zmin);
ux = plP_w3wcx(x[i], y[0], z[i][0]);
uy = plP_w3wcy(x[i], y[0], z[i][0]);
pljoin(tx, ty, ux, uy);
}
}
if (opt != 2) {
for (i = 0; i < ny; i++) {
tx = plP_w3wcx(x[nx - 1], y[i], zmin);
ty = plP_w3wcy(x[nx - 1], y[i], zmin);
ux = plP_w3wcx(x[nx - 1], y[i], z[nx - 1][i]);
uy = plP_w3wcy(x[nx - 1], y[i], z[nx - 1][i]);
pljoin(tx, ty, ux, uy);
}
}
}
}
/*--------------------------------------------------------------------------*\
* void plgrid3()
*
* This routine draws a grid around the back side of the 3d plot with
* hidden line removal.
\*--------------------------------------------------------------------------*/
static void
plgrid3(PLFLT tick)
{
PLFLT xmin, ymin, zmin, xmax, ymax, zmax, zscale;
PLFLT cxx, cxy, cyx, cyy, cyz, zmin_in, zmax_in;
PLINT u[3], v[3];
PLINT nsub = 0;
PLFLT tp;
plP_gw3wc(&cxx, &cxy, &cyx, &cyy, &cyz);
plP_gdom(&xmin, &xmax, &ymin, &ymax);
plP_grange(&zscale, &zmin_in, &zmax_in);
zmin = (zmax_in > zmin_in) ? zmin_in: zmax_in;
zmax = (zmax_in > zmin_in) ? zmax_in: zmin_in;
pldtik(zmin, zmax, &tick, &nsub);
tp = tick * floor(zmin / tick) + tick;
pl3upv = 0;
if (cxx >= 0.0 && cxy <= 0.0) {
while (tp <= zmax) {
u[0] = plP_wcpcx(plP_w3wcx(xmin, ymax, tp));
v[0] = plP_wcpcy(plP_w3wcy(xmin, ymax, tp));
u[1] = plP_wcpcx(plP_w3wcx(xmax, ymax, tp));
v[1] = plP_wcpcy(plP_w3wcy(xmax, ymax, tp));
u[2] = plP_wcpcx(plP_w3wcx(xmax, ymin, tp));
v[2] = plP_wcpcy(plP_w3wcy(xmax, ymin, tp));
plnxtv(u, v, 0, 3, 0);
tp += tick;
}
u[0] = plP_wcpcx(plP_w3wcx(xmax, ymax, zmin));
v[0] = plP_wcpcy(plP_w3wcy(xmax, ymax, zmin));
u[1] = plP_wcpcx(plP_w3wcx(xmax, ymax, zmax));
v[1] = plP_wcpcy(plP_w3wcy(xmax, ymax, zmax));
plnxtv(u, v, 0, 2, 0);
}
else if (cxx <= 0.0 && cxy <= 0.0) {
while (tp <= zmax) {
u[0] = plP_wcpcx(plP_w3wcx(xmax, ymax, tp));
v[0] = plP_wcpcy(plP_w3wcy(xmax, ymax, tp));
u[1] = plP_wcpcx(plP_w3wcx(xmax, ymin, tp));
v[1] = plP_wcpcy(plP_w3wcy(xmax, ymin, tp));
u[2] = plP_wcpcx(plP_w3wcx(xmin, ymin, tp));
v[2] = plP_wcpcy(plP_w3wcy(xmin, ymin, tp));
plnxtv(u, v, 0,3, 0);
tp += tick;
}
u[0] = plP_wcpcx(plP_w3wcx(xmax, ymin, zmin));
v[0] = plP_wcpcy(plP_w3wcy(xmax, ymin, zmin));
u[1] = plP_wcpcx(plP_w3wcx(xmax, ymin, zmax));
v[1] = plP_wcpcy(plP_w3wcy(xmax, ymin, zmax));
plnxtv(u, v, 0,2, 0);
}
else if (cxx <= 0.0 && cxy >= 0.0) {
while (tp <= zmax) {
u[0] = plP_wcpcx(plP_w3wcx(xmax, ymin, tp));
v[0] = plP_wcpcy(plP_w3wcy(xmax, ymin, tp));
u[1] = plP_wcpcx(plP_w3wcx(xmin, ymin, tp));
v[1] = plP_wcpcy(plP_w3wcy(xmin, ymin, tp));
u[2] = plP_wcpcx(plP_w3wcx(xmin, ymax, tp));
v[2] = plP_wcpcy(plP_w3wcy(xmin, ymax, tp));
plnxtv(u, v, 0,3, 0);
tp += tick;
}
u[0] = plP_wcpcx(plP_w3wcx(xmin, ymin, zmin));
v[0] = plP_wcpcy(plP_w3wcy(xmin, ymin, zmin));
u[1] = plP_wcpcx(plP_w3wcx(xmin, ymin, zmax));
v[1] = plP_wcpcy(plP_w3wcy(xmin, ymin, zmax));
plnxtv(u, v, 0,2, 0);
}
else if (cxx >= 0.0 && cxy >= 0.0) {
while (tp <= zmax) {
u[0] = plP_wcpcx(plP_w3wcx(xmin, ymin, tp));
v[0] = plP_wcpcy(plP_w3wcy(xmin, ymin, tp));
u[1] = plP_wcpcx(plP_w3wcx(xmin, ymax, tp));
v[1] = plP_wcpcy(plP_w3wcy(xmin, ymax, tp));
u[2] = plP_wcpcx(plP_w3wcx(xmax, ymax, tp));
v[2] = plP_wcpcy(plP_w3wcy(xmax, ymax, tp));
plnxtv(u, v, 0,3, 0);
tp += tick;
}
u[0] = plP_wcpcx(plP_w3wcx(xmin, ymax, zmin));
v[0] = plP_wcpcy(plP_w3wcy(xmin, ymax, zmin));
u[1] = plP_wcpcx(plP_w3wcx(xmin, ymax, zmax));
v[1] = plP_wcpcy(plP_w3wcy(xmin, ymax, zmax));
plnxtv(u, v, 0,2, 0);
}
pl3upv = 1;
}
/*--------------------------------------------------------------------------*\
* void plnxtv()
*
* Draw the next view of a 3-d plot. The physical coordinates of the
* points for the next view are placed in the n points of arrays u and
* v. The silhouette found so far is stored in the heap as a set of m peak
* points.
*
* These routines dynamically allocate memory for hidden line removal.
* Memory is allocated in blocks of 2*BINC*sizeof(PLINT) bytes. Large
* values of BINC give better performance but also allocate more memory
* than is needed. If your 3D plots are very "spiky" or you are working
* with very large matrices then you will probably want to increase BINC.
\*--------------------------------------------------------------------------*/
static void
plnxtv(PLINT *u, PLINT *v, PLFLT* c, PLINT n, PLINT init)
{
plnxtvhi(u, v, c, n, init);
if (pl3mode)
plnxtvlo(u, v, c, n, init);
}
/*--------------------------------------------------------------------------*\
* void plnxtvhi()
*
* Draw the top side of the 3-d plot.
\*--------------------------------------------------------------------------*/
static void
plnxtvhi(PLINT *u, PLINT *v, PLFLT* c, PLINT n, PLINT init)
{
/*
* For the initial set of points, just display them and store them as the
* peak points.
*/
if (init == 1) {
int i;
oldhiview = (PLINT *) malloc((size_t) (2 * n * sizeof(PLINT)));
if ( ! oldhiview)
myexit("plnxtvhi: Out of memory.");
oldhiview[0] = u[0];
oldhiview[1] = v[0];
plP_draw3d(u[0], v[0], c, 0, 1);
for (i = 1; i < n; i++) {
oldhiview[2 * i] = u[i];
oldhiview[2 * i + 1] = v[i];
plP_draw3d(u[i], v[i], c, i, 0);
}
mhi = n;
return;
}
/*
* Otherwise, we need to consider hidden-line removal problem. We scan
* over the points in both the old (i.e. oldhiview[]) and new (i.e. u[]
* and v[]) arrays in order of increasing x coordinate. At each stage, we
* find the line segment in the other array (if one exists) that straddles
* the x coordinate of the point. We have to determine if the point lies
* above or below the line segment, and to check if the below/above status
* has changed since the last point.
*
* If pl3upv = 0 we do not update the view, this is useful for drawing
* lines on the graph after we are done plotting points. Hidden line
* removal is still done, but the view is not updated.
*/
xxhi = 0;
if (pl3upv != 0) {
newhisize = 2 * (mhi + BINC);
if (newhiview != NULL) {
newhiview =
(PLINT *) realloc((void *) newhiview,
(size_t) (newhisize * sizeof(PLINT)));
}
else {
newhiview =
(PLINT *) malloc((size_t) (newhisize * sizeof(PLINT)));
}
if ( ! newhiview)
myexit("plnxtvhi: Out of memory.");
}
/* Do the draw or shading with hidden line removal */
plnxtvhi_draw(u, v, c, n);
/* Set oldhiview */
swaphiview();
}
/*--------------------------------------------------------------------------*\
* void plnxtvhi_draw()
*
* Draw the top side of the 3-d plot.
\*--------------------------------------------------------------------------*/
static void
plnxtvhi_draw(PLINT *u, PLINT *v, PLFLT* c, PLINT n)
{
PLINT i = 0, j = 0, first = 1;
PLINT sx1 = 0, sx2 = 0, sy1 = 0, sy2 = 0;
PLINT su1, su2, sv1, sv2;
PLINT cx, cy, px, py;
PLINT seg, ptold, lstold = 0, pthi, pnewhi = 0, newhi, change, ochange = 0;
/*
* (oldhiview[2*i], oldhiview[2*i]) is the i'th point in the old array
* (u[j], v[j]) is the j'th point in the new array
*/
/*
* First attempt at 3d shading. It works ok for simple plots, but
* will just not draw faces, or draw them overlapping for very
* jagged plots
*/
while (i < mhi || j < n) {
/*
* The coordinates of the point under consideration are (px,py). The
* line segment joins (sx1,sy1) to (sx2,sy2). "ptold" is true if the
* point lies in the old array. We set it by comparing the x coordinates
* of the i'th old point and the j'th new point, being careful if we
* have fallen past the edges. Having found the point, load up the point
* and segment coordinates appropriately.
*/
ptold = (j >= n || (i < mhi && oldhiview[2 * i] < u[j]));
if (ptold) {
px = oldhiview[2 * i];
py = oldhiview[2 * i + 1];
seg = j > 0 && j < n;
if (seg) {
sx1 = u[j - 1];
sy1 = v[j - 1];
sx2 = u[j];
sy2 = v[j];
}
} else {
px = u[j];
py = v[j];
seg = i > 0 && i < mhi;
if (seg) {
sx1 = oldhiview[2 * (i - 1)];
sy1 = oldhiview[2 * (i - 1) + 1];
sx2 = oldhiview[2 * i];
sy2 = oldhiview[2 * i + 1];
}
}
/*
* Now determine if the point is higher than the segment, using the
* logical function "above". We also need to know if it is the old view
* or the new view that is higher. "newhi" is set true if the new view
* is higher than the old.
*/
if (seg)
pthi = plabv(px, py, sx1, sy1, sx2, sy2);
else
pthi = 1;
newhi = (ptold && !pthi) || (!ptold && pthi);
/*
* The last point and this point lie on different sides of
* the current silouette
*/
change = (newhi && !pnewhi) || (!newhi && pnewhi);
/*
* There is a new intersection point to put in the peak array if the
* state of "newhi" changes.
*/
if (first) {
plP_draw3d(px, py, c, j, 1);
first = 0;
lstold = ptold;
savehipoint(px, py);
pthi = 0;
ochange = 0;
} else if (change) {
/*
* Take care of special cases at end of arrays. If pl3upv is 0 the
* endpoints are not connected to the old view.
*/
if (pl3upv == 0 && ((!ptold && j == 0) || (ptold && i == 0))) {
plP_draw3d(px, py, c, j, 1);
lstold = ptold;
pthi = 0;
ochange = 0;
} else if (pl3upv == 0 &&
(( ! ptold && i >= mhi) || (ptold && j >= n))) {
plP_draw3d(px, py, c, j, 1);
lstold = ptold;
pthi = 0;
ochange = 0;
} else {
/*
* If pl3upv is not 0 then we do want to connect the current line
* with the previous view at the endpoints. Also find intersection
* point with old view.
*/
if (i == 0) {
sx1 = oldhiview[0];
sy1 = -1;
sx2 = oldhiview[0];
sy2 = oldhiview[1];
} else if (i >= mhi) {
sx1 = oldhiview[2 * (mhi - 1)];
sy1 = oldhiview[2 * (mhi - 1) + 1];
sx2 = oldhiview[2 * (mhi - 1)];
sy2 = -1;
} else {
sx1 = oldhiview[2 * (i - 1)];
sy1 = oldhiview[2 * (i - 1) + 1];
sx2 = oldhiview[2 * i];
sy2 = oldhiview[2 * i + 1];
}
if (j == 0) {
su1 = u[0];
sv1 = -1;
su2 = u[0];
sv2 = v[0];
} else if (j >= n) {
su1 = u[n - 1];
sv1 = v[n - 1];
su2 = u[n - 1];
sv2 = -1;
} else {
su1 = u[j - 1];
sv1 = v[j - 1];
su2 = u[j];
sv2 = v[j];
}
/* Determine the intersection */
pl3cut(sx1, sy1, sx2, sy2, su1, sv1, su2, sv2, &cx, &cy);
if (cx == px && cy == py) {
if (lstold && !ochange)
plP_draw3d(px, py, c, j, 1);
else
plP_draw3d(px, py, c, j, 0);
savehipoint(px, py);
lstold = 1;
pthi = 0;
} else {
if (lstold && !ochange)
plP_draw3d(cx, cy, c, j, 1);
else
plP_draw3d(cx, cy, c, j, 0);
lstold = 1;
savehipoint(cx, cy);
}
ochange = 1;
}
}
/* If point is high then draw plot to point and update view. */
if (pthi) {
if (lstold && ptold)
plP_draw3d(px, py, c, j, 1);
else
plP_draw3d(px, py, c, j, 0);
savehipoint(px, py);
lstold = ptold;
ochange = 0;
}
pnewhi = newhi;
if (ptold)
i++;
else
j++;
}
}
/*--------------------------------------------------------------------------*\
* void plP_draw3d()
*
* Does a simple move or line draw.
\*--------------------------------------------------------------------------*/
static void
plP_draw3d(PLINT x, PLINT y, PLFLT *c, PLINT j, PLINT move)
{
if (move)
plP_movphy(x, y);
else {
if (c != NULL)
plcol1(c[j-1]);
plP_draphy(x, y);
}
}
/*--------------------------------------------------------------------------*\
* void plnxtvlo()
*
* Draw the bottom side of the 3-d plot.
\*--------------------------------------------------------------------------*/
static void
plnxtvlo(PLINT *u, PLINT *v, PLFLT*c, PLINT n, PLINT init)
{
PLINT i, j, first;
PLINT sx1 = 0, sx2 = 0, sy1 = 0, sy2 = 0;
PLINT su1, su2, sv1, sv2;
PLINT cx, cy, px, py;
PLINT seg, ptold, lstold = 0, ptlo, pnewlo, newlo, change, ochange = 0;
first = 1;
pnewlo = 0;
/*
* For the initial set of points, just display them and store them as the
* peak points.
*/
if (init == 1) {
oldloview = (PLINT *) malloc((size_t) (2 * n * sizeof(PLINT)));
if ( ! oldloview)
myexit("\nplnxtvlo: Out of memory.");
plP_draw3d(u[0], v[0], c, 0, 1);
oldloview[0] = u[0];
oldloview[1] = v[0];
for (i = 1; i < n; i++) {
plP_draw3d(u[i], v[i], c, i, 0);
oldloview[2 * i] = u[i];
oldloview[2 * i + 1] = v[i];
}
mlo = n;
return;
}
/*
* Otherwise, we need to consider hidden-line removal problem. We scan
* over the points in both the old (i.e. oldloview[]) and new (i.e. u[]
* and v[]) arrays in order of increasing x coordinate. At each stage, we
* find the line segment in the other array (if one exists) that straddles
* the x coordinate of the point. We have to determine if the point lies
* above or below the line segment, and to check if the below/above status
* has changed since the last point.
*
* If pl3upv = 0 we do not update the view, this is useful for drawing
* lines on the graph after we are done plotting points. Hidden line
* removal is still done, but the view is not updated.
*/
xxlo = 0;
i = 0;
j = 0;
if (pl3upv != 0) {
newlosize = 2 * (mlo + BINC);
if (newloview != NULL) {
newloview =
(PLINT *) realloc((void *) newloview,
(size_t) (newlosize * sizeof(PLINT)));
}
else {
newloview =
(PLINT *) malloc((size_t) (newlosize * sizeof(PLINT)));
}
if ( ! newloview)
myexit("plnxtvlo: Out of memory.");
}
/*
* (oldloview[2*i], oldloview[2*i]) is the i'th point in the old array
* (u[j], v[j]) is the j'th point in the new array.
*/
while (i < mlo || j < n) {
/*
* The coordinates of the point under consideration are (px,py). The
* line segment joins (sx1,sy1) to (sx2,sy2). "ptold" is true if the
* point lies in the old array. We set it by comparing the x coordinates
* of the i'th old point and the j'th new point, being careful if we
* have fallen past the edges. Having found the point, load up the point
* and segment coordinates appropriately.
*/
ptold = (j >= n || (i < mlo && oldloview[2 * i] < u[j]));
if (ptold) {
px = oldloview[2 * i];
py = oldloview[2 * i + 1];
seg = j > 0 && j < n;
if (seg) {
sx1 = u[j - 1];
sy1 = v[j - 1];
sx2 = u[j];
sy2 = v[j];
}
}
else {
px = u[j];
py = v[j];
seg = i > 0 && i < mlo;
if (seg) {
sx1 = oldloview[2 * (i - 1)];
sy1 = oldloview[2 * (i - 1) + 1];
sx2 = oldloview[2 * i];
sy2 = oldloview[2 * i + 1];
}
}
/*
* Now determine if the point is lower than the segment, using the
* logical function "above". We also need to know if it is the old view
* or the new view that is lower. "newlo" is set true if the new view is
* lower than the old.
*/
if (seg)
ptlo = !plabv(px, py, sx1, sy1, sx2, sy2);
else
ptlo = 1;
newlo = (ptold && !ptlo) || (!ptold && ptlo);
change = (newlo && !pnewlo) || (!newlo && pnewlo);
/*
* There is a new intersection point to put in the peak array if the
* state of "newlo" changes.
*/
if (first) {
plP_draw3d(px, py, c, j, 1);
first = 0;
lstold = ptold;
savelopoint(px, py);
ptlo = 0;
ochange = 0;
}
else if (change) {
/*
* Take care of special cases at end of arrays. If pl3upv is 0 the
* endpoints are not connected to the old view.
*/
if (pl3upv == 0 && ((!ptold && j == 0) || (ptold && i == 0))) {
plP_draw3d(px, py, c, j, 1);
lstold = ptold;
ptlo = 0;
ochange = 0;
}
else if (pl3upv == 0 &&
(( ! ptold && i >= mlo) || (ptold && j >= n))) {
plP_draw3d(px, py, c, j, 1);
lstold = ptold;
ptlo = 0;
ochange = 0;
}
/*
* If pl3upv is not 0 then we do want to connect the current line
* with the previous view at the endpoints. Also find intersection
* point with old view.
*/
else {
if (i == 0) {
sx1 = oldloview[0];
sy1 = 100000;
sx2 = oldloview[0];
sy2 = oldloview[1];
}
else if (i >= mlo) {
sx1 = oldloview[2 * (mlo - 1)];
sy1 = oldloview[2 * (mlo - 1) + 1];
sx2 = oldloview[2 * (mlo - 1)];
sy2 = 100000;
}
else {
sx1 = oldloview[2 * (i - 1)];
sy1 = oldloview[2 * (i - 1) + 1];
sx2 = oldloview[2 * i];
sy2 = oldloview[2 * i + 1];
}
if (j == 0) {
su1 = u[0];
sv1 = 100000;
su2 = u[0];
sv2 = v[0];
}
else if (j >= n) {
su1 = u[n - 1];
sv1 = v[n - 1];
su2 = u[n];
sv2 = 100000;
}
else {
su1 = u[j - 1];
sv1 = v[j - 1];
su2 = u[j];
sv2 = v[j];
}
/* Determine the intersection */
pl3cut(sx1, sy1, sx2, sy2, su1, sv1, su2, sv2, &cx, &cy);
if (cx == px && cy == py) {
if (lstold && !ochange)
plP_draw3d(px, py, c, j, 1);
else
plP_draw3d(px, py, c, j, 0);
savelopoint(px, py);
lstold = 1;
ptlo = 0;
}
else {
if (lstold && !ochange)
plP_draw3d(cx, cy, c, j, 1);
else
plP_draw3d(cx, cy, c, j, 0);
lstold = 1;
savelopoint(cx, cy);
}
ochange = 1;
}
}
/* If point is low then draw plot to point and update view. */
if (ptlo) {
if (lstold && ptold)
plP_draw3d(px, py, c, j, 1);
else
plP_draw3d(px, py, c, j, 0);
savelopoint(px, py);
lstold = ptold;
ochange = 0;
}
pnewlo = newlo;
if (ptold)
i = i + 1;
else
j = j + 1;
}
/* Set oldloview */
swaploview();
}
/*--------------------------------------------------------------------------*\
* savehipoint
* savelopoint
*
* Add a point to the list of currently visible peaks/valleys, when
* drawing the top/bottom surface, respectively.
\*--------------------------------------------------------------------------*/
static void
savehipoint(PLINT px, PLINT py)
{
if (pl3upv == 0)
return;
if (xxhi >= newhisize) { /* allocate additional space */
newhisize += 2 * BINC;
newhiview = (PLINT *) realloc((void *) newhiview,
(size_t) (newhisize * sizeof(PLINT)));
if ( ! newhiview)
myexit("savehipoint: Out of memory.");
}
newhiview[xxhi] = px;
xxhi++;
newhiview[xxhi] = py;
xxhi++;
}
static void
savelopoint(PLINT px, PLINT py)
{
if (pl3upv == 0)
return;
if (xxlo >= newlosize) { /* allocate additional space */
newlosize += 2 * BINC;
newloview = (PLINT *) realloc((void *) newloview,
(size_t) (newlosize * sizeof(PLINT)));
if ( ! newloview)
myexit("savelopoint: Out of memory.");
}
newloview[xxlo] = px;
xxlo++;
newloview[xxlo] = py;
xxlo++;
}
/*--------------------------------------------------------------------------*\
* swaphiview
* swaploview
*
* Swaps the top/bottom views. Need to do a real swap so that the
* memory cleanup routine really frees everything (and only once).
\*--------------------------------------------------------------------------*/
static void
swaphiview(void)
{
PLINT *tmp;
if (pl3upv != 0) {
mhi = xxhi / 2;
tmp = oldhiview;
oldhiview = newhiview;
newhiview = tmp;
}
}
static void
swaploview(void)
{
PLINT *tmp;
if (pl3upv != 0) {
mlo = xxlo / 2;
tmp = oldloview;
oldloview = newloview;
newloview = tmp;
}
}
/*--------------------------------------------------------------------------*\
* freework
*
* Frees memory associated with work arrays
\*--------------------------------------------------------------------------*/
static void
freework(void)
{
free_mem(oldhiview);
free_mem(oldloview);
free_mem(newhiview);
free_mem(newloview);
free_mem(vtmp);
free_mem(utmp);
free_mem(ctmp);
}
/*--------------------------------------------------------------------------*\
* myexit
*
* Calls plexit, cleaning up first.
\*--------------------------------------------------------------------------*/
static void
myexit(char *msg)
{
freework();
plexit(msg);
}
/*--------------------------------------------------------------------------*\
* myabort
*
* Calls plabort, cleaning up first.
* Caller should return to the user program.
\*--------------------------------------------------------------------------*/
static void
myabort(char *msg)
{
freework();
plabort(msg);
}
/*--------------------------------------------------------------------------*\
* int plabv()
*
* Determines if point (px,py) lies above the line joining (sx1,sy1) to
* (sx2,sy2). It only works correctly if sx1 <= px <= sx2.
\*--------------------------------------------------------------------------*/
static int
plabv(PLINT px, PLINT py, PLINT sx1, PLINT sy1, PLINT sx2, PLINT sy2)
{
int above;
if (py >= sy1 && py >= sy2)
above = 1;
else if (py < sy1 && py < sy2)
above = 0;
else if ((double) (sx2 - sx1) * (py - sy1) >=
(double) (px - sx1) * (sy2 - sy1))
above = 1;
else
above = 0;
return above;
}
/*--------------------------------------------------------------------------*\
* void pl3cut()
*
* Determines the point of intersection (cx,cy) between the line joining
* (sx1,sy1) to (sx2,sy2) and the line joining (su1,sv1) to (su2,sv2).
\*--------------------------------------------------------------------------*/
static void
pl3cut(PLINT sx1, PLINT sy1, PLINT sx2, PLINT sy2,
PLINT su1, PLINT sv1, PLINT su2, PLINT sv2, PLINT *cx, PLINT *cy)
{
PLINT x21, y21, u21, v21, yv1, xu1, a, b;
double fa, fb;
x21 = sx2 - sx1;
y21 = sy2 - sy1;
u21 = su2 - su1;
v21 = sv2 - sv1;
yv1 = sy1 - sv1;
xu1 = sx1 - su1;
a = x21 * v21 - y21 * u21;
fa = (double) a;
if (a == 0) {
if (sx2 < su2) {
*cx = sx2;
*cy = sy2;
}
else {
*cx = su2;
*cy = sv2;
}
return;
}
else {
b = yv1 * u21 - xu1 * v21;
fb = (double) b;
*cx = (PLINT) (sx1 + (fb * x21) / fa + .5);
*cy = (PLINT) (sy1 + (fb * y21) / fa + .5);
}
}
syntax highlighted by Code2HTML, v. 0.9.1