/* 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[][0..ly-1]. The integer "opt" */
/* specifies: */
/*  opt = 1:  Draw lines parallel to x-axis */
/*  opt = 2:  Draw lines parallel to y-axis */
/*  opt = 3:  Draw lines parallel to both axes */

#include "plplot.h"
#include <math.h>

#ifdef PLSTDC
#include <stdlib.h>
#else
extern char *malloc();
extern void free();
#endif

extern PLINT *oldhiview;
PLINT zbackflag=0, zbcolor;
PLFLT zticksp;

void plot3d(x,y,z,ly,nx,ny,opt,side)
PLINT ly, nx, ny, opt, side;
PLFLT *x, *y, *z;
{
      PLINT b, color, font;
      PLFLT cxx, cxy, cyx, cyy, cyz;
      PLINT init;
      PLINT *work;
      PLINT i, ix, iy;
      PLINT level;

      glev(&level);
      if (level < 3) plexit("Please set up window before calling plot3d");

      if (opt<1 || opt>3) plexit("Bad option in plot3d");
      if (nx<=0 || ny<=0 || ly<ny) plexit("Bad array dimensions in plot3d.");

      /* Check that points in x and in y are strictly increasing */

      for (i=0; i<nx-1; i++)
        if (x[i]>=x[i+1]) plexit("X array must be strictly increasing in plot3d");

      for (i=0; i<ny-1; i++)
        if (y[i]>=y[i+1]) plexit("Y array must be strictly increasing in plot3d");

      work = (PLINT *)malloc(4*max(nx,ny)*sizeof(PLINT));
      if(!work)
         plexit("Out of memory in plot3d.");
      b = 2*max(nx,ny)+1;
      gw3wc(&cxx,&cxy,&cyx,&cyy,&cyz);
      init = 1;

      if (cxx >= 0.0 && cxy <= 0.0) {
        if (opt != 2) {
          for (iy=2; iy<=ny; iy++)  {
            plt3zz(1,iy,1,-1,-opt,init,x,y,z,ly,nx,ny,&work[0],&work[b-1]);
            init = 0;
          }
        }
        else  {
          plt3zz(1,ny,1,-1,-opt,init,x,y,z,ly,nx,ny,&work[0],&work[b-1]);
          init = 0;
        }
        if (opt != 1)
          for (ix=1; ix<=nx-1; ix++)
            plt3zz(ix,ny,1,-1,opt,init,x,y,z,ly,nx,ny,&work[0],&work[b-1]);
        else
          plt3zz(1,ny,1,-1,opt,init,x,y,z,ly,nx,ny,&work[0],&work[b-1]);
      }
      else if (cxx <= 0.0 && cxy <= 0.0) {
        if (opt != 1) {
          for (ix=2; ix<=nx; ix++) {
            plt3zz(ix,ny,-1,-1,opt,init,x,y,z,ly,nx,ny,&work[0],&work[b-1]);
            init = 0;
          }
        }
        else  {
          plt3zz(nx,ny,-1,-1,opt,init,x,y,z,ly,nx,ny,&work[0],&work[b-1]);
          init = 0;
        }
        if (opt != 2)
          for (iy=ny; iy>=2; iy--)
            plt3zz(nx,iy,-1,-1,-opt,init,x,y,z,ly,nx,ny,&work[0],&work[b-1]);
        else
          plt3zz(nx,ny,-1,-1,-opt,init,x,y,z,ly,nx,ny,&work[0],&work[b-1]);
      }
      else if (cxx <= 0.0 && cxy >= 0.0) {
        if (opt != 2) {
          for (iy=ny-1; iy>=1; iy--) {
            plt3zz(nx,iy,-1,1,-opt,init,x,y,z,ly,nx,ny,&work[0],&work[b-1]);
            init = 0;
          }
        }
        else {
          plt3zz(nx,1,-1,1,-opt,init,x,y,z,ly,nx,ny,&work[0],&work[b-1]);
          init = 0;
        }
        if (opt != 1)
          for (ix=nx; ix>=2; ix--)
            plt3zz(ix,1,-1,1,opt,init,x,y,z,ly,nx,ny,&work[0],&work[b-1]);
        else
          plt3zz(nx,1,-1,1,opt,init,x,y,z,ly,nx,ny,&work[0],&work[b-1]);
      }
      else if (cxx >= 0.0 && cxy >= 0.0) {
        if (opt != 1) {
          for (ix=nx-1; ix>=1; ix--) {
            plt3zz(ix,1,1,1,opt,init,x,y,z,ly,nx,ny,&work[0],&work[b-1]);
            init = 0;
          }
        }
        else  {
          plt3zz(1,1,1,1,opt,init,x,y,z,ly,nx,ny,&work[0],&work[b-1]);
          init = 0;
        }
        if (opt != 2)
          for (iy=1; iy<=ny-1; iy++)
            plt3zz(1,iy,1,1,-opt,init,x,y,z,ly,nx,ny,&work[0],&work[b-1]);
        else
          plt3zz(1,1,1,1,-opt,init,x,y,z,ly,nx,ny,&work[0],&work[b-1]);
      }

      free((VOID *)work);

      if(side)
         plside3a(x,y,z,ly,nx,ny,opt);

      if(zbackflag) {
         gatt(&font,&color);
         plcol(zbcolor);
         plgrid3a(zticksp);
         plcol(color);
      }

      free((VOID *)oldhiview);
}