/* $Id: plot3d.c,v 1.3 2007/05/08 09:09:37 rice 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	(const char *);
static void myabort	(const 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);

/* pmr: not used ... removed 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 */
#ifdef MJL_HACK
#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
#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 xx0, PLFLT yy0, PLFLT zz0,
	       PLFLT xx1, PLFLT yy1, PLFLT zz1,
	       PLFLT xx2, PLFLT yy2, PLFLT zz2)
{
  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] = xx0; x[1] = xx1; x[2] = xx2;
  y[0] = yy0; y[1] = yy1; y[2] = yy2;
  z[0] = zz0; z[1] = zz1; z[2] = zz2;
  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;

  (void) ixstart;			/* pmr: make these used */
  (void) ixn;
  (void) indexymin;
  (void) indexymax;

  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 *myzz = (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;
	  myzz = (PLFLT *) realloc(zz, np*sizeof(PLFLT));
	}
	for (j=0; j<cline->npts; j++)
	  myzz[j] = plsc->ranmi;
	if (cline->npts > 0) {
	   plcol1((clev->level-fc_minz)/(fc_maxz-fc_minz));
	   plline3(cline->npts, cline->x, cline->y, myzz);
	}
	cline = cline->next;
      }
      while(cline != NULL);
      clev = clev->next;
    }
    while(clev != NULL);

    cont_clean_store(cont); /* now release the memory */
    free(myzz);
  }

  /* 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 myzscale, myzmin, myzmax;

    plP_grange(&myzscale, &myzmin, &myzmax);

    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], myzmin);
      shade_triangle(px[2], py[2], pz[2], px[2], py[2], myzmin,  px[0], py[0], myzmin);
    }

    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], myzmin);
      shade_triangle(px[2], py[2], pz[2], px[2], py[2], myzmin,  px[0], py[0], myzmin);
    }
  }
}

/*--------------------------------------------------------------------------*\
 * 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;

    (void) ixstart;			/* pmr: make these used */
    (void) ixn;
    (void) indexymin;
    (void) indexymax;

    pl3mode = 0;

    if (plsc->level < 3) {
	myabort("plot3dcl: Please set up window first");
	return;
    }

    if ((opt & 3) == 0) {
	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 ii, jj;

      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(ii=1; ii<_nx-1; ii++)
	_x[ii] = x[ixmin+ii];
      _y[0] = ymin;
      _y[_ny-1] = ymax;
      for(ii=1; ii<_ny-1; ii++)
	_y[ii] = y[iymin+ii];

      /* copy the data array so we can interpolate around the edges */
      for(ii=0; ii<_nx; ii++)
	_z[ii] = (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(ii=0; ii<_nx; ii++) {
	if(ii==0) {
	  _z[ii][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(jj=1; jj<_ny-1; jj++)
	    _z[ii][jj] = z[ixmin][iymin+jj]*(1-tx0) + z[ixmin+1][iymin+jj]*tx0;
	  _z[ii][_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(ii==_nx-1) {
	  _z[ii][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(jj=1; jj<_ny-1; jj++)
	    _z[ii][jj] = z[ixmax-1][iymin+jj]*(1-tx1) + z[ixmax][iymin+jj]*tx1;
	  _z[ii][_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[ii][0] = z[ixmin+ii][iymin]*(1-ty0) + z[ixmin+ii][iymin+1]*ty0;
	  for(jj=1; jj<_ny-1; jj++)
	    _z[ii][jj] = z[ixmin+ii][iymin+jj];
	  _z[ii][_ny-1] = z[ixmin+ii][iymax-1]*(1-ty1) + z[ixmin+ii][iymax]*ty1;
	}
	for(jj=0; jj<_ny; jj++) {
	  if(_z[ii][jj] < zmin)
	    _z[ii][jj] = zmin;
	  else if(_z[ii][jj] > zmax)
	    _z[ii][jj] = 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, jj;
      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 (jj = 0; jj < ny; jj++) {
	    cgrid2.xg[i][jj] = x[i];
	    cgrid2.yg[i][jj] = y[jj];
	    zstore[i][jj] = z[i][jj];
	 }
      }

      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, ii, 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 */
 	  ii = 0;
	  if (cline->npts > 0) {
	    do {
	      plcol1((clev->level-fc_minz)/(fc_maxz-fc_minz));
	      cx =  plP_wcpcx(plP_w3wcx(cline->x[ii],cline->y[ii], plsc->ranmi));
	      for (jj=ii; jj < cline->npts; jj++) {  /* convert to 2D coordinates */
		uu[jj] = plP_wcpcx(plP_w3wcx(cline->x[jj],cline->y[jj], plsc->ranmi));
		vv[jj] = plP_wcpcy(plP_w3wcy(cline->x[jj],cline->y[jj], plsc->ranmi));
		if (uu[jj] < cx) /* find turn back point */
		  break;
		else
		  cx = uu[jj];
	      }
	      plnxtv(&uu[ii], &vv[ii], NULL, jj-ii, 0); /* plot line with increasing x */
	      
	      if (jj < cline->npts) { /* line not yet finished, */
		start = jj-1;
		for (ii = start; ii < cline->npts; ii++) {  /* search turn forward point */
		  uu[ii] = plP_wcpcx(plP_w3wcx(cline->x[ii],cline->y[ii], plsc->ranmi));
		  if (uu[ii] > cx)
		    break;
		  else
		    cx = uu[ii];
		}
		end = ii-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 (jj=start; jj <= end; jj++) {
		  uu[jj] = plP_wcpcx(plP_w3wcx(cline->x[jj],cline->y[jj], plsc->ranmi));
		  vv[jj] = plP_wcpcy(plP_w3wcy(cline->x[jj],cline->y[jj], 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];
		ii = end; /* restart where it was left */
	      }
	    } while (jj < 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.
\*--------------------------------------------------------------------------*/

void					/* pmr: not used - removed static */
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 xx0, PLINT yy0, 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 <= xx0 && xx0 <= nx && 1 <= yy0 && yy0 <= ny) {
	x2d = plP_w3wcx(x[xx0 - 1], y[yy0 - 1], z[xx0 - 1][yy0 - 1]);
	y2d = plP_w3wcy(x[xx0 - 1], y[yy0 - 1], z[xx0 - 1][yy0 - 1]);
	u[n] = plP_wcpcx(x2d);
	v[n] = plP_wcpcy(y2d);
	if (c != NULL)
	  c[n] = (z[xx0 - 1][yy0 - 1] - fc_minz)/(fc_maxz-fc_minz);

	switch (flag) {
	case -3:
	    yy0 += dy;
	    flag = -flag;
	    break;
	case -2:
	    yy0 += dy;
	    break;
	case -1:
	    yy0 += dy;
	    flag = -flag;
	    break;
	case 1:
	    xx0 += dx;
	    break;
	case 2:
	    xx0 += dx;
	    flag = -flag;
	    break;
	case 3:
	    xx0 += dx;
	    flag = -flag;
	    break;
	}
	n++;
    }

    if (flag == 1 || flag == -2) {
	if (flag == 1) {
	    xx0 -= dx;
	    yy0 += dy;
	}
	else if (flag == -2) {
	    yy0 -= dy;
	    xx0 += dx;
	}
	if (1 <= xx0 && xx0 <= nx && 1 <= yy0 && yy0 <= ny) {
	    x2d = plP_w3wcx( x[xx0 - 1], y[yy0 - 1], z[xx0 - 1][yy0 - 1]);
	    y2d = plP_w3wcy( x[xx0 - 1], y[yy0 - 1], z[xx0 - 1][yy0 - 1]);
	    u[n] = plP_wcpcx(x2d);
	    v[n] = plP_wcpcy(y2d);
	    if (c != NULL)
	      c[n] = (z[xx0 - 1][yy0 - 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(const char *msg)
{
    freework();
    plexit(msg);
}

/*--------------------------------------------------------------------------*\
 * myabort
 *
 * Calls plabort, cleaning up first.
 * Caller should return to the user program.
\*--------------------------------------------------------------------------*/

static void
myabort(const 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);
    }
}