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

/* "work" is an integer work array of size at least 4*max(nx,ny) */

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

void plot3d(x,y,z,work,ly,nx,ny,opt)
int ly, nx, ny, opt;
float x[],y[],*z;
int work[];
{
      int b;
      float cxx, cxy, cyx, cyy, cyz;

      int init;
      int i, ix, iy;

      int level;
      glev(&level);
      if (level < 3) fatal("Please set up window before calling PLOT3D");
 
      if (opt<1 || opt>3) fatal("Bad option in PLOT3D");
      if (nx<=0 || ny<=0 || ly<ny) fatal("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]) fatal("X(*) must be strictly increasing in PLOT3D");

      for (i=0; i<ny-1; i++)
        if (y[i]>=y[i+1]) fatal("Y(*) must be strictly increasing 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]);
      }      
}