# map3d6.rb

require "narray"
require "numru/dcl"
include NMath
include NumRu

nx=37; ny=37
xmin=0.0;  xmax=360.0;  ymin= -90.0;  ymax= 90.0
vxmin=-0.4; vxmax=0.4; vymin=-0.3; vymax=0.3
zmin=-1.9;  zmax=1.9
vzmin=-0.2; vzmax=0.2
drad=PI/180

alon = xmin + (xmax-xmin)/(nx-1) * NArray.sfloat(nx).indgen!
alat = ymin + (ymax-ymin)/(ny-1) * NArray.sfloat(ny).indgen!

slat = sin(alat*drad).reshape!(1,ny)
p = 3*sqrt(1.0-slat**2)*slat*cos(alon*drad) - 0.5*(3*slat**2-1)

iws = (ARGV[0] || (puts ' WORKSTATION ID (I)  ? ;'; DCL::sgpwsn; gets)).to_i
DCL::sgopn iws
DCL::sgfrm

#-- X-Y PLANE ----

DCL::sgswnd(  xmin,  xmax,  ymin,  ymax )
DCL::sgsvpt( vxmin, vxmax, vymin, vymax )
DCL::sgstrn( 1 )
DCL::sgstrf

DCL::scspln( 1, 2, vzmin)
DCL::scseye( -1.1, -1.1, 2.5 )
DCL::scsobj(  0.0,  0.0, 0.0 )
DCL::scsprj

DCL::uxaxdv( 'B', 10.0, 60.0 )
DCL::uxaxdv( 'T', 10.0, 60.0 )
DCL::uxsttl( 'B', 'LONGITUDE', 0.0 )
DCL::uyaxdv( 'L', 10.0, 30.0 )
DCL::uyaxdv( 'R', 10.0, 30.0 )
DCL::uysttl( 'L', 'LATITUDE', 0.0 )

#-- X-Z PLANE ----

DCL::sgswnd(  xmin,  xmax,  zmin,  zmax )
DCL::sgsvpt( vxmin, vxmax, vzmin, vzmax )
DCL::sgstrn( 1 )
DCL::sgstrf

DCL::scspln( 1, 3, vymax)
DCL::scsprj

DCL::uzinit
DCL::uxaxdv( 'T', 10.0, 60.0 )
DCL::uzlset( 'LABELXB', false )
DCL::uxaxdv( 'B', 10.0, 60.0 )
DCL::uyaxdv( 'L', 0.2, 1.0 )
DCL::uyaxdv( 'R', 0.2, 1.0 )
DCL::uysttl( 'L', 'AMPLITUDE', 0.0 )

#-- Y-Z PLANE ----

DCL::sgswnd(  ymin,  ymax,  zmin,  zmax )
DCL::sgsvpt( vymin, vymax, vzmin, vzmax )
DCL::sgstrn( 1 )
DCL::sgstrf

DCL::scspln( 2, 3, vxmax)
DCL::scsprj

DCL::uzinit
DCL::uxaxdv( 'T', 10.0, 30.0 )
DCL::uxaxdv( 'B', 10.0, 30.0 )
DCL::uzlset( 'LABELYL', false )
DCL::uyaxdv( 'L', 0.2, 1.0 )
DCL::uyaxdv( 'R', 0.2, 1.0 )

#---------------- 3-D ------------------

DCL::scsvpt(vxmin, vxmax, vymin, vymax, vzmin, vzmax)
DCL::scswnd( xmin,  xmax,  ymin,  ymax,  zmin,  zmax)
DCL::scslog(false, false, false)
DCL::scstrn(1)
DCL::scstrf

(ny-2).downto(0){|j|
  (nx-2).downto(0){|i|
    xp = [ alon[i], alon[i],   alon[i+1],  alon[i+1], alon[i] ]
    yp = [ alat[j], alat[j+1], alat[j+1],  alat[j],   alat[j] ]
    zp = [ p[i,j],  p[i,j+1],  p[i+1,j+1], p[i+1,j],  p[i,j] ]
    DCL::scplu(xp, yp, zp)
 }
}

DCL::sgcls