/* CMD: Surface Plot
 *
 * 2-D Function Surface Maker for Modeler
 * Originally by Arnie Cachelin © 1992 NewTek Inc., Sat Jul 11 1992
 *
 * Modified to use Modeler-based UI by Stuart Ferguson, 11/92
 * Further modified by Arnie, Sat May  8 16:04:25 1993
 */

arg bas
call addlib "LWModelerARexx.port", 0
address "LWModelerARexx.port"

SIGNAL ON ERROR
SIGNAL ON SYNTAX
MATHLIB="rexxmathlib.library"
IF POS(MATHLIB , SHOW('L')) = 0 THEN
  IF ~ADDLIB(MATHLIB , 0 , -30 , 0) THEN DO
	call notify(1,"!Can't find "MATHLIB)
	exit
	END
call addlib "rexxsupport.library", 0, -30, 0

sysnam = 'Plot 2D Function'
filnam = 'ENV:plot2d.state'
version = 'Plot 2D v1.1'

/* Setup state.  Read stored one, if any.
 */
NSX  = 20
NSY  = 20
xmin = -10
xmax = 10
ymin = -10
ymax = 10
func = 'x*y/10'
ufunc=func
fn = 6
flip = 1
tri = 1
if bas="" then bas = 0
bas = 0
typ=2
if (exists(filnam)) then do
    if (~open(state, filnam, 'R')) then break
    if (readln(state) ~= version) then break
    parse value readln(state) with nsx nsy xmin xmax ymin ymax fn tri typ .
    func = readln(state)
    call close state
end

FnList.1='RandomSheet'
FnList.2='Wave'
FnList.3='RadialWave'
FnList.4='Gaussian'
FnList.5='Interfere'
FnList.6='Custom'
FList= FnList.1 FnList.2 FnList.3 FnList.4 FnList.5 FnList.6

call req_begin sysnam

id_fnl = req_addcontrol("f(x,y)", 'CH',FList)
id_lox = req_addcontrol("Low X", 'n', 1)
id_hix = req_addcontrol("High X", 'n', 1)
id_nsx = req_addcontrol("X Segments", 'n')
id_loy = req_addcontrol("Low Y", 'n', 1)
id_hiy = req_addcontrol("High Y", 'n', 1)
id_nsy = req_addcontrol("Y Segments", 'n')
id_fun = req_addcontrol("Custom Function", 's', 35)
/* id_hgt = req_addcontrol("Height Map", 'b') */
id_tri = req_addcontrol("Triangles", 'b')
id_typ = req_addcontrol("Build: ","CH","Points Polys Curves")
/* id_bas = req_addcontrol("Base", 'b') */

call req_setval id_lox, xmin, -10
call req_setval id_loy, ymin, -10
call req_setval id_hix, xmax, 10
call req_setval id_hiy, ymax, 10
call req_setval id_nsx, nsx, 20
call req_setval id_nsy, nsy, 20
call req_setval id_fun, func
call req_setval id_fnl, fn,6
/* call req_setval id_hgt, flip */
call req_setval id_tri, tri,tri
call req_setval id_typ, typ,typ
/* call req_setval id_bas, bas,bas */

if (~req_post()) then do
    call req_end
    exit
end

NSX = req_getval(id_nsx) % 1
NSY = req_getval(id_nsx) % 1
xmin = req_getval(id_lox)
xmax = req_getval(id_hix)
ymin = req_getval(id_loy)
ymax = req_getval(id_hiy)
fn   = req_getval(id_fnl)
typ   = req_getval(id_typ)
if fn=6 then func = req_getval(id_fun)
else func = FnList.fn'(x,y)'
/* flip = req_getval(id_hgt) */
tri = req_getval(id_tri)

call req_end

/* Save state now, in case something fails.
 */
if (open(state, filnam, 'W')) then do
    call writeln state, version
    call writeln state, nsx nsy xmin xmax ymin ymax fn tri typ
    call writeln state, ufunc
    call close state
end

/* Prep
 */
xrange = xmax - xmin
yrange = ymax - ymin
xmesh = xrange / NSX
ymesh = yrange / NSY
rscale = sqrt(xrange*yrange)
tri_height = sqrt(3)/2
ifunc = "z =" func
say ifunc
z=0
zmax=z
zmin=z
call randu(time('s'))  /* Seed random number generator */

if typ=3 then call RectCurves
else
if tri then do
  if typ=1 then call TriPoints
  else call TriMesh
  end
else do
  if typ=1 then call RectPoints
  else call RectMesh
  end

/* if Bas totalpoints=totalpoints+MakeBase() */

l1 = "Points created:" totalpoints
l2 = "Polygons created:" poly
l3 = "Z ranges between" zmin "and" zmax
call notify 1, '!'sysnam, l1, l2, l3

exit

RectMesh:
    totalpoints = (NSX+1) * (NSY+1)
    totalpolys  = NSX * NSY
    call add_begin
    call meter_begin totalpoints+2, sysnam, "Computing "totalpoints" points for "totalpolys" squares"
    do y=ymin to ymax by ymesh
      if y=ymax then TopCorner.4=totalpoints
      do x=xmin to xmax by xmesh
        interpret ifunc
        if z<zmin then zmin=z
        if z>zmax then zmax=z     /* Just some silly stats for later */
        if (flip) then vec =x z y
        else vec = x y z
        call add_point(vec)
        call meter_step
      end
      if y=ymin then TopCorner.2=totalpoints
      if y=ymax then TopCorner.3=totalpoints
    end

    point=1
    poly=0
    call meter_begin totalpolys, sysnam, "Generating "totalpolys" Polygon Mesh"
    do y=ymin to ymax-ymesh by ymesh   /* Don't make wrap-around polygon */
      do x=xmin to xmax by xmesh
        if x<xmax then do   /* Again, Don't make wrap-around polygon! */
          if (flip) then
            call add_quad point point+NSX+1 point+NSX+2 point+1
          else
            call add_quad point point+1 point+NSX+2 point+NSX+1
          poly = poly + 1
          call meter_step
        end
        point = point + 1
      end
    end
    call meter_end
    call add_end
return totalpoints
/*  */

RectPoints:
    totalpoints = (NSX+1) * (NSY+1)
    poly=totalpoints
    call add_begin
    call meter_begin totalpoints+2, sysnam, "Computing "totalpoints" points for 1-point polygons"
    point=1
    do y=ymin to ymax by ymesh
      do x=xmin to xmax by xmesh
        interpret ifunc
        if z<zmin then zmin=z
        if z>zmax then zmax=z     /* Just some silly stats for later */
        if (flip) then vec =x z y
        else vec = x y z
        call add_point(vec)
        call add_polygon(point)
        point=point+1
        call meter_step
      end
    end
    call add_end
return totalpoints
/*  */

RectCurves:
    totalpoints = (NSX+1) * (NSY+1)
    poly=NSX+1 + NSY+1
    call add_begin
    call meter_begin totalpoints+NSX+2, sysnam, "Computing "totalpoints" points for "poly" Curves"
    crv=""
    point=1
    do y=ymin to ymax by ymesh
      do x=xmin to xmax by xmesh

        interpret ifunc

        if z<zmin then zmin=z
        if z>zmax then zmax=z

        if (flip) then vec =x z y
        else vec = x y z
        call add_point(vec)
        crv=crv point
        point=point+1
        call meter_step
      end
      call Add_Curve(crv)
      crv=""
    end
    do p=1 to NSX+1
      do o=0 to NSY
        crv=crv p+o*(NSX+1)
      end
      call meter_step
      call Add_Curve(crv)
      crv=""
    end
    call add_end
return totalpoints
/*  */

TriMesh:
    totalpoints = (NSX+1) * (NSY+1)
    totalpolys  = NSX * NSY * 2
    call add_begin
    call meter_begin totalpoints*2, sysnam, "Computing "totalpoints" points"
    offset=0
    rows=0
    totalpoints=0

    do y=ymin to ymax by ymesh* tri_height
    	rows=rows+1
	    columns=0
      if y=ymax then TopCorner.4=totalpoints
      do x=xmin+offset to xmax+offset by xmesh

     		columns = columns + 1
        interpret ifunc

        if z<zmin then zmin=z
        if z>zmax then zmax=z     /* Just some silly stats for later */

        if (flip) then vec =x z y
        else vec = x y z
        call add_point vec
        totalpoints = totalpoints + 1
        call meter_step
      end
      if y=ymin then TopCorner.2=totalpoints
      if y=ymax then TopCorner.3=totalpoints
	    if offset=0 then offset=.5 * xmesh  /* offset alternate lines */
	    else offset=0
    end
    call meter_end

    point=1
    poly=0
    off=0
    call meter_begin totalpolys, sysnam, "Generating "totalpolys" Polygon Mesh"
    do row=0 to rows-2
      if off=0 then off=1  /* Boy this feels kludgey!!! */
      else off=0
     	do col=1 to columns - 1
        if (flip) then do
          call add_quad col+row*columns col+(row*columns)+1 col+((row+1)*columns)+abs(off-1)
          call add_quad col+(row*columns)+off col+((row+1)*columns)+1 col+((row+1)*columns)
          poly=poly+2
          end
        else do
          call add_quad col+row*columns col+((row+1)*columns)+abs(off-1) col+(row*columns)+1
          call add_quad col+(row*columns)+off col+((row+1)*columns) col+((row+1)*columns)+1
          poly=poly+2
          end
        call meter_step
      end
    end
    call meter_end
    call add_end
return totalpoints
/*  */

TriPoints:
    totalpoints = (NSX+1) * (NSY+1)
    call add_begin
    call meter_begin totalpoints*2, sysnam, "Computing "totalpoints" points for 1-point polygons"
    offset=0
    rows=0
    totalpoints=0
    point=1

    do y=ymin to ymax by ymesh* tri_height
    	rows=rows+1
	    columns=0
      do x=xmin+offset to xmax+offset by xmesh

     		columns = columns + 1
        interpret ifunc

        if z<zmin then zmin=z
        if z>zmax then zmax=z     /* Just some silly stats for later */

        if (flip) then vec =x z y
        else vec = x y z
        call add_point vec
        call add_polygon point
        point=point+1
        call meter_step
      end
	    if offset=0 then offset=.5 * xmesh  /* offset alternate lines */
	    else offset=0
    end
    call meter_end
    totalpoints = point - 1
    poly=totalpoints
    call add_end()
return totalpoints
/*  */


MakeBase:
    add_begin
    point=TotalPoints-1
    say point
    TopCorner.1=1
    BotCorner.1=totalpoints
    BotCorner.2=BotCorner.1 + NSX +1
    BotCorner.3=BotCorner.2 + NSY
    BotCorner.4=BotCorner.3 + NSX +1
    z=ZMin-(xrange+yrange)/4
    y=ymin
    BotCorner.1=point+1
    do x=xmin to xmax by xmesh
      call add_point x y z
    	 point=point+1
       pointList.point = x y z
    end
    BotCorner.2=point
    x=xmin
    do y=ymin to ymax by ymesh
      call add_point x y z
    	 point=point+1
       pointList.point = x y z
    end
    y=ymax
    BotCorner.4=point+1
    do x=xmin to xmax by xmesh
      call add_point x y z
    	 point=point+1
       pointList.point = x y z
    end
    BotCorner.3=point
    x=xmax
    do y=ymin to ymax by ymesh
      call add_point x y z
    	 point=point+1
       pointList.point = x y z
    end

    call surface("AreaBottom")
    call add_quad BotCorner.1 BotCorner.2 BotCorner.3 BotCorner.4
    poly=poly+1
    call surface("AreaBase")

    do i=0 to NSX-1
      call add_quad TopCorner.1+i TopCorner.1+i+1 BotCorner.1+i+1 BotCorner.1+i
      poly=poly+1
    end

    do i=0 to NSY-1
      call add_quad TopCorner.4+i TopCorner.4+i+1 BotCorner.4+i+1 BotCorner.4+i
      poly=poly+1
    end

    do i=0 to NSX-1
      a=((TopCorner.1)+i*(NSX+1))
    	b=((TopCorner.1)+(i+1)*(NSX+1))
    	c=((BotCorner.2)+i+1)
    	d=((BotCorner.2)+i+2)
      call add_quad a b d c
      poly=poly+1
    end

    do i=0 to NSY-1
      a=((TopCorner.2)+i*(NSY+1))
    	b=((TopCorner.2)+(i+1)*(NSY+1))
    	c=((BotCorner.3)+i+1)
    	d=((BotCorner.3)+i+2)
      call add_quad a b d c
      poly=poly+1
    end

    add_end
return point
    /*  */

SYNTAX:
ERROR:
  say 'Sorry, Error #'RC' on line 'SIGL' has been detected.'
  say errortext(rc)
  t=Notify(1,'!Total Bummer Dude!','!An error has been detected.','@'ErrorText(rc),'Line 'SIGL)
  call end_all
  exit

Radius: PROCEDURE
  arg xf, yf
  return sqrt(xf*xf+yf*yf)


Sinc:   PROCEDURE EXPOSE rscale
  arg xf, yf     /* Classic the spherical Bessel f'n j0 */
  r=Radius(xf,yf)
  if (r = 0)
    zf = rscale
  else
    zf = rscale * sin(r) / r
  return zf


Wave:   PROCEDURE  EXPOSE rscale
  arg xf, yf     /* Simple wavy sheet */
  zf=rscale*sin(2*3.141592*xf/rscale)/3
  return(zf)


RadialWave:   PROCEDURE  EXPOSE rscale
  arg xf, yf     /* rings of waves sheet */
  xc=0; yc=0     /* Center coord.s */
  r=Radius(xf-xc,yf-yc)
  zf=rscale*sin(10*3.141592*r/rscale)/5
  return(zf)


Interfere:   PROCEDURE  EXPOSE rscale
  arg xf, yf     /* Interference of several (3) sources */
  r=Radius(xf,yf)      /* (0,0)  */
  d1=Radius(xf-5*rscale/10,yf+2*rscale/10) /* (5,-2) */
  d2=Radius(xf+4*rscale/10,yf-3*rscale/10) /* (-4,3) */
  zf=(rscale/7)*(sin(r)+sin(d1*1.5)+1.4*cos(d2))
  return(zf)


RandomSheet:   PROCEDURE  EXPOSE rscale
  arg xf, yf     /* Random altitudes */
  amp=(rscale/3)
  zf=randu()*amp-amp/2
  return(zf)


Polynomial:   PROCEDURE  EXPOSE rscale
  arg xf, yf
  zf = 2*xf*xf + 3*xf - 2*yf*yf + 8*yf + 2
  zf = zf / 20
  return(zf)


Gaussian:   PROCEDURE  EXPOSE rscale
  arg xf, yf
  xc=0; yc=0     /* Center coord.s */
  amp = rscale/3
  r = Radius(xf-xc,yf-yc)
  zf = amp * exp(-r*r/20)
  return(zf)