Contents of file named plot3d.f

      program plot3D
C     Grid extends symmetrically from zero
C        for nxgr cells in X and nygr in Y
      parameter (nxgr=50,nygr=50)
      parameter (nx=1+2*nxgr, ny=1+2*nygr)
C     Parameters: depth of valley (fscl), X and Y ranges (4x2)
      parameter (fscl=1, xscl=4, yscl=2)
C     Data for grid.  Grid layout is
C        (*, 1) - lowest streak of y values
C        (*,ny) - highest streak of y values
C        (1, *) - lowest x value at any y position
C        (nx,*) - highest x value at any y position
      real data(nx,ny)

C     Create data:  Curved valley defined by a Rosenbrock fcn
      f(x,y) = (1-x)**2 + fscl*(y-x**2)**2

C     Run over grid and insert elevation information into it.

      do j=1,ny
         y = yscl*float(j-1-nygr)/nygr
         do i=1,nx
            x = xscl*float(i-1-nxgr)/nxgr
            data(i,j) = f(x,y)
         enddo
      enddo

C     Create default SAC file header; modify to make it 3-D data
      call newhdr

C     File type to XYZ for 3D data
      call setihv('IFTYPE', 'IXYZ', nerr)

C     Grid dimensions: NXSIZE, NYSIZE 
      call setnhv('NXSIZE', nx, nerr)
      call setnhv('NYSIZE', ny, nerr)

C     Grid scale:  XMINIMUM, XMAXIMUM, YMINIMUM, YMAXIMUM
      call setfhv('XMINIMUM', -xscl, nerr)
      call setfhv('XMAXIMUM',  xscl, nerr)
      call setfhv('YMINIMUM', -yscl, nerr)
      call setfhv('YMAXIMUM',  yscl, nerr)

C     Number of data points, begin, delta (required but irrelevant)
      call setfhv('B', 0.0, nerr)
      call setfhv('DELTA', 1.0, nerr)
      call setnhv('NPTS', nx*ny, nerr)

C     Write data to file using present header values
      call wsac0('/tmp/data3D.xyz',data,data,nerr)
      end