let progname$="wiggb"
let scrsize=.5
let sketch=1

  def rsmall=1
  def rlarge=3
  def x(u,v)= (rlarge+rsmall*cos(v))*cos(u)
  def y(u,v)= (rlarge+rsmall*cos(v))*sin(u)
  def z(u,v)= rsmall*sin(v)  + .5*sin(6*u+shift)

!!!!!!! basic  parameters !!!!!!!!!!!!!!!!!!
let a =  0          !  bounds
let b =  2*pi       !  a < u < b
let c =  0          !  c < v < d
let d =  2*pi
let n= 64      !  grid mesh to divide rectangle [a,b] x [c,d]
let m= 24
let du=(b-a)/n
let dv=(d-c)/m
call eye(6,3,2)             ! eyeball point
call lightbulb(6,6,6)       ! light source - used for shading of grid patches
let zoom = 3                ! quick way to resize things
let xmin=-1.5  * zoom       ! 2dim screen rectangle
let xmax= 1.5  * zoom       ! it lies on plane with normal set by eye
let ymin=  -1  * zoom
let ymax=   1  * zoom
let zmin=  -8  * zoom       ! only effects distance calculations
let zmax=   8  * zoom       !  used by paint storage routines
let xaxis=1.75 * zoom       ! endpoints of axis and position of labels
let yaxis=1.5  * zoom
let zaxis=1    * zoom
let numcon=16               ! number of contour planes
let maxcon=zaxis            ! largest contour plane
let mincon=-zaxis           ! smallest contour plane

! other global variables are objq,maxq,ys,zs,xpix,ypix,pix,mode$
! color names (see sub init)
!!!!!!!!!!!!!!!!!!!!! arrays
option base 0                   ! arrays start with 0
dim xscrn(2),yscrn(2),zscrn(2)  ! unit orthogonal vectors of screen plane
dim bulb(2)                     ! light source
dim index(1),scratch(1)         ! used by sort
dim object$(1)

sub redimension(maxpics)
 mat redim object$(maxpics),index(maxpics),scratch(maxpics)
end sub                      !  call redimension(1) to free up memory

!!!!!!!!!!!!!!!!!!!!!!!!
!!! 3 - dim routines !!!
!!!!!!!!!!!!!!!!!!!!!!!!
def xs(a,b,c)          ! transform into screen coordinates:  xs,ys,zs
   let ys = a*yscrn[0]+b*yscrn[1]+c*yscrn[2]    ! exploits side effect
   let zs = a*zscrn[0]+b*zscrn[1]+c*zscrn[2]    ! ys and zs are passed
   let xs = a*xscrn[0]+b*xscrn[1]+c*xscrn[2]    ! as global variables
end def

def norm(a()) = sqr(dot(a,a))

sub normalize(a())
  local norma
  let norma=norm(a())
  if norma <>0  then mat a=(1/norma)*a
end sub

sub eye(a,b,c)     ! Finds xscrn,yscrn orthogonal to (a,b,c) and such
  let xscrn[0]=-b  ! that the z-axis is straight up and down.
  let xscrn[1]= a  !  zscrn points to eye
  let xscrn[2]= 0  ! This makes z-axis straight up and down.
  let yscrn[0]= -a
  let yscrn[1]= -b
  let yscrn[2]= (a*a + b*b)/c   ! needs c not equal zero
  let zscrn[0]=a
  let zscrn[1]=b
  let zscrn[2]=c
  call normalize(xscrn)
  call normalize(yscrn)
  call normalize(zscrn)
end sub

sub lightbulb(a,b,c)        ! bulb is used to determine hue of grid patches
   let bulb(0)=a            ! if bulb perpindicular to patch use bright
   let bulb(1)=b            ! colors, otherwise use darker
   let bulb(2)=c
   call normalize(bulb)
end sub

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!! initialize graphics !!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

sub init
  set mode "vga"
  when error in
    open #1 : screen 0,scrsize,0,scrsize           ! in case #1 already open
  use
    close #1
    open #1 : screen 0,scrsize,0,scrsize
  end when
  window #1
  set window xmin,xmax,ymin,ymax
  set background color "white"
  clear
  let black=0
  let green=3
  let blue=10
  let yellow=13
  let red=14
  let white=15
  set color mix(0) 1,1,1        ! background=white
  set color mix(1) 0,1/3,0      ! dark green
  set color mix(2) 0,2/3,0
  set color mix(3) 1/3,2/3,0
  set color mix(4) 0,1,0
  set color mix(5) 1/3,1,0
  set color mix(6) 2/3,1,0      ! light green
  set color mix(7) 0,0,1/3      ! dark blue
  set color mix(8) 0,0,2/3
  set color mix(9) 0,1/3,2/3
  set color mix(10) 0,0,1
  set color mix(11) 0,1/3,1
  set color mix(12) 0,2/3,1     ! light blue
  set color mix(13) 1,1,0       ! yellow
  set color mix(14) 1,0,0       ! red
  set color mix(15) 0,0,0       ! white=black
end sub

def gethue(side,intensity)
  if side > 0 then
   let gethue=int(7+5.99*intensity)  ! blue 7-12
  else
   let gethue=int(1+5.99*intensity)  ! green 1-6
  end if
end def

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!  paint routines !!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

let xpix=640       ! two bytes each packb uses signed so clip ok
let ypix=480       ! object$(q) =  type distance pointx pointy hue
let zpix=256*256   !  where type,hue one byte, others two bytes

def type$(q)   =object$(q)[1:1]
sub settype(q,a$)                 !!!!! should be called first!!!!
  let object$(q)[1:1]=a$          ! g(rid),l(ine),t(ext),s(kip)
  if len(object$(q))<8 then let object$(q)=object$(q) & "1234567"
  ! this ensures packb and unpackb work correctly
end sub

def distance(q)=unpackb(object$(q),9,16)  ! unsigned - no clip protect
sub setdistance(q,w)
  local temp
  let temp=(w-zmin)*zpix/(zmax-zmin)  ! convert to pixel coord
  call packb(object$(q),9,16,temp)
end sub

def pointx(q)  =unpackb(object$(q),25,-16) ! signed integer
sub setpointx(q,w)
  local temp
  let temp=(w-xmin)*xpix/(xmax-xmin)  ! convert to pixel coord
  call packb(object$(q),25,-16,temp)
end sub

def pointy(q)  =unpackb(object$(q),41,-16)
sub setpointy(q,w)
  local temp, temp$
  let temp=(w-ymin)*ypix/(ymax-ymin)  ! convert to pixel coord
  call packb(object$(q),41,-16,temp)
end sub

def hue(q)     =unpackb(object$(q)[8:8],1,8)
sub sethue(q,w)
  let object$(q)[8:8]=chr$(w)
end sub

def text$(q)   =object$(q)[9:100]
sub settext(q,a$)
  let object$(q)[9:100]=a$
end sub

sub paint(minq,maxq)
 local q,xpt1,xpt2,xpt3,xpt4,ypt1,ypt2,ypt3,ypt4,indq
 set text justify "center","half"
 set window 1,xpix,1,ypix          ! fake pixel coordinates
 for q=minq to maxq
    let indq=index(q)
    select case type$(indq)
      case "g"                    ! grid
        let xpt1=pointx(indq)
        let ypt1=pointy(indq)
        let xpt2=pointx(indq+1)
        let ypt2=pointy(indq+1)
        let xpt3=pointx(indq+m+2)
        let ypt3=pointy(indq+m+2)
        let xpt4=pointx(indq+m+1)
        let ypt4=pointy(indq+m+1)
        set color hue(indq)
        plot area: xpt1,ypt1 ; xpt2,ypt2 ; xpt3,ypt3 ; xpt4,ypt4
        set color red
        plot  xpt1,ypt1 ; xpt2,ypt2 ; xpt3,ypt3 ; xpt4,ypt4 ;  xpt1,ypt1
      case "l"               ! line
        set color yellow
        plot pointx(indq),pointy(indq);pointx(indq+1),pointy(indq+1)
      case "s"               ! skip - no plotting -outer grid & end of lines
      case "t"               ! text
       set color white
       plot text, at pointx(indq),pointy(indq) : text$(indq)
    end select
!   call savescr(q)
  next q
 set window xmin,xmax,ymin,ymax  ! standard coordinates
end sub

!!!!!!!!!!!!!! line, text, and axis routines !!!!!!!!!!!!!
sub loadsegment(x1,y1,z1,x2,y2,z2,n)
  local q,dx,dy,dz,x,y,z,qstart,qfinish,temp1,temp2
  let dx=(x2-x1)/n    ! loads line segment from point1 to point2
  let dy=(y2-y1)/n    ! broken into n subsegments
  let dz=(z2-z1)/n
  let x=x1
  let y=y1
  let z=z1
  let qstart=objq
  let qfinish=objq+n
  for q=qstart to qfinish
    call settype(q,"l")        ! line
    call setpointx(q,xs(x,y,z))
    call setpointy(q,ys)
    let temp1=zs
    let temp2=xs(x+dx,y+dy,z+dz)
    let temp2=zs
    call setdistance(q,.5*(temp1+temp2)+fudge)
    let x=x+dx
    let y=y+dy
    let z=z+dz
  next q
  call settype(qfinish,"s")    ! skip - don't plot anything at end of segment
  let objq=qfinish+1
end sub

sub loadtext(a,b,c,label$)   ! loads label$ at point (a,b,c)
  call settype(objq,"t")          ! text
  call setpointx(objq,xs(a,b,c))
  call setpointy(objq,ys)
  call setdistance(objq,zs)
  call settext(objq,label$)
  let objq=objq+1
end sub

sub loadaxis(a,b,c)
  call loadsegment(-a,0,0,a,0,0,n+m)
  call loadsegment(0,-b,0,0,b,0,n+m)
  call loadsegment(0,0,-c,0,0,c,n+m)
  call loadtext(0,0,c,"Z")
  call loadtext(0,b,0,"Y")
  call loadtext(a,0,0,"X")
end sub

!!!!!!!!!!!! load surface grid !!!!!!!!!!!!!!
sub loadgrid
  local u,v,i,j,u1,u2,u3,v1,v2,v3,side,intensity,normal(2)
  let u=a
  for i=0 to n
    let v=c
    for j=0 to m
      call settype(objq,"g")
      if i=n or j=m then call settype(objq,"s")  ! don't draw last row-col
      call setpointx(objq,xs(x(u,v),y(u,v),z(u,v)))
      call setpointy(objq,ys)
      call setdistance(objq,zs)
      let u1=x(u+du,v)-x(u,v)
      let u2=y(u+du,v)-y(u,v)
      let u3=z(u+du,v)-z(u,v)
      let v1=x(u,v+dv)-x(u,v)
      let v2=y(u,v+dv)-y(u,v)
      let v3=z(u,v+dv)-z(u,v)
      let normal(0)=u2*v3-u3*v2    ! cross product
      let normal(1)=u3*v1-u1*v3
      let normal(2)=u1*v2-u2*v1
      call normalize(normal())     ! unit normal vector to patch
      let side=dot(normal,zscrn)   ! sign determines which side we see
      let intensity=abs(dot(normal,bulb))    !  intensity
      call sethue(objq,gethue(side,intensity))
      let v=v+dv
      let objq=objq+1
    next j
    let u=u+du
  next i
end sub

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!! contour routines !!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
def between(a,b,x)
   if a<=b then
     if (a<=x) and (x<=b) then let between=1 else let between=0
   else
     if (b<=x) and (x<=a) then let between=1 else let between=0
   end if
end def

sub loadcontourarray
 local i,u,v,x1,x2,x3,x4,y1,y2,y3,y4,z1,z2,z3,z4
 local mn,mx,k,ck,dck,scalar,p1,p2,q1,q2,temp
 let dck= (maxcon-mincon)/(numcon-1)
 let u=a
 for i=0 to n-1
  let v=c
  for j=0 to m-1
    let x1=x(u,v)
    let y1=y(u,v)
    let z1=z(u,v)
    let x2=x(u+du,v)
    let y2=y(u+du,v)
    let z2=z(u+du,v)
    let x3=x(u+du,v+dv)
    let y3=y(u+du,v+dv)
    let z3=z(u+du,v+dv)
    let x4=x(u,v+dv)
    let y4=y(u,v+dv)
    let z4=z(u,v+dv)
    let mx=max(z1,max(z2,max(z3,z4)))
    let mn=min(z1,min(z2,min(z3,z4)))
    for ck=mincon to maxcon step dck
      if between(mx,mn,ck)=1 then
        let p1,p2,q1,q2=0
        if between(z1,z2,ck)=1 then
            let scaler=(ck-z1)/(z2-z1)
            let p1=x1+scaler*(x2-x1)
            let p2=y1+scaler*(y2-y1)
            if      between(z2,z3,ck)=1 then
              let scaler=(ck-z2)/(z3-z2)
              let q1=x2+scaler*(x3-x2)
              let q2=y2+scaler*(y3-y2)
            else if between(z3,z4,ck)=1 then
              let scaler=(ck-z3)/(z4-z3)
              let q1=x3+scaler*(x4-x3)
              let q2=y3+scaler*(y4-y3)
            else
              let scaler=(ck-z4)/(z1-z4)
              let q1=x4+scaler*(x1-x4)
              let q2=y4+scaler*(y1-y4)
            end if
        else if between(z2,z3,ck)=1 then
            let scaler=(ck-z2)/(z3-z2)
            let p1=x2+scaler*(x3-x2)
            let p2=y2+scaler*(y3-y2)
            if      between(z3,z4,ck)=1 then
              let scaler=(ck-z3)/(z4-z3)
              let q1=x3+scaler*(x4-x3)
              let q2=y3+scaler*(y4-y3)
            else
              let scaler=(ck-z4)/(z1-z4)
              let q1=x4+scaler*(x1-x4)
              let q2=y4+scaler*(y1-y4)
            end if
        else if between(z3,z4,ck)=1 then
            let scaler=(ck-z3)/(z4-z3)
            let p1=x3+scaler*(x4-x3)
            let p2=y3+scaler*(y4-y3)
            let scaler=(ck-z4)/(z1-z4)
            let q1=x4+scaler*(x1-x4)
            let q2=y4+scaler*(y1-y4)
        end if
        !  essentially  loadsegment(p1,p2,ck,q1,q2,ck,1)
        !  with a little fudge= .25(du+dv)
        call settype(objq,"l")
        let temp=xs(x1,y1,z1)                  ! dummy call to get zs
        call setdistance(objq,zs +.25*(du+dv)) ! just in front of grid 
        call setpointx(objq,xs(p1,p2,ck))
        call setpointy(objq,ys)
         let objq=objq+1
        call settype(objq,"s")
        call setpointx(objq,xs(q1,q2,ck))
        call setpointy(objq,ys)
         let objq=objq+1
      end if
    next ck
    let v=v+dv
  next j
  let u=u+du
 next i
end sub


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!  sorting !!!!!!!!!!!!!!!
sub loadindex(index(),n)
  local i
  for i=0 to n
    let index(i)=i
  next i
end sub

sub merge(index(),p,m,q,scratch())
  local i,j,k
  let i=p
  let j=m+1
  let k=p
  do until i=m+1 or j=q+1
    if distance(index(i))< distance(index(j)) then
      let scratch(k)=index(i)
      let i=i+1
    else
      let scratch(k)=index(j)
      let j=j+1
    end if
    let k=k+1
  loop
  if i=m+1 then
     do until j=q+1
       let scratch(k)=index(j)
       let j=j+1
       let k=k+1
     loop
  else
     do until i=m+1
       let scratch(k)=index(i)
       let i=i+1
       let k=k+1
     loop
  end if
  for i=p to  q
    let index(i)=scratch(i)
  next i
end sub     ! merge

sub sort(index(),p,q,scratch())
  local m
  if (q-p)<2 then
     call merge(index,p,p,q,scratch)
  else
     let m=int(p+(q-p)/2)
     call sort(index,p,m,scratch)
     call sort(index,m+1,q,scratch)
     call merge(index,p,m,q,scratch)
  end if
end sub

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!! pbmplus save screen   !!!!!
!! globals : ppmflip, progname$!
!!  ppmnumb, #80, colormix$  !!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

dim colormix$(0 to 15)

sub getcolormix
  local point,r,g,b
  for point=0 to 15
    ask color mix (point) r,g,b
    let colormix$(point)=using$("# # #", 3*r,3*g,3*b)
  next point
end sub

sub initppm  ! save first file of transparent gifs, writes progname10.ppm
  local x,y,a$,px,py,xmin,xmax,ymin,ymax,point
  ask window xmin,xmax,ymin,ymax    ! save window coordinates
  ask pixels px,py
  set window 1,px,1,py
  call getcolormix
  open #77 : name progname$&"10.ppm", create "newold"
  erase #77
  print #77: "P3",px,py,3   ! P3=ppm, width, height, max rgb 0-3
  open #78 : name "temp0.tmp", create "newold"  ! scratch file
  erase #78
  for y= py to 1 step -1
   for x= 1 to px
     box keep x,x,y,y in a$
     let point = min(1,ord(a$[9:9]))+ 2*min(1,ord(a$[10:10]))
     let point = point + 4*min(1,ord(a$[11:11]))+8*min(1,ord(a$[12:12]))
     print #77: colormix$(point)
     print #78: colormix$(point)
   next x
  next y
  set window xmin,xmax,ymin,ymax  ! restore window coordinates
  close #77
  close #78
  open  #80 : name progname$&".scr", create "newold"
  erase #80       ! writes an animator script file
  print #80 : ">OUTPUT " & progname$&".gif"
  print #80 : ">WIDTH ";px
  print #80 : ">HEIGHT ";py
  print #80 : ">NSLOOP"
  print #80 : ">FRAME NOTHING"
  print #80 : "FILE "; progname$&"10.gif"
  print #80 : "TRANS #000000"     ! black = transparent color
  print #80 : "TIME 0"
  print #80 : "VER 0"
  print #80 : "HOR 0"
  let ppmnumb=10
  let ppmflip=0
end sub

sub moreppm  ! keeps track of transparent color, writes prognamexx.ppm
  local x,y,a$,b$,px,py,xmin,xmax,ymin,ymax,point
  ask window xmin,xmax,ymin,ymax    ! save window coordinates
  ask pixels px,py
  set window 1,px,1,py
  let ppmnumb=ppmnumb+1
  open #77 : name progname$&str$(ppmnumb)&".ppm", create "newold"
  erase #77
  print #77: "P3",px,py,3   ! P3=ppm, width, height, max rgb 0-3
  open #78 : name "temp"&str$(ppmflip)&".tmp", create "newold"  ! old scratch
  let ppmflip=1-ppmflip
  open #79 : name "temp"&str$(ppmflip)&".tmp", create "newold"  ! new scratch
  erase #79
  for y= py to 1 step -1
   for x= 1 to px
     box keep x,x,y,y in a$
     let point = min(1,ord(a$[9:9]))+ 2*min(1,ord(a$[10:10]))
     let point = point + 4*min(1,ord(a$[11:11]))+8*min(1,ord(a$[12:12]))
     let a$=colormix$(point)
     print #79: a$                      ! put screen pix into new scratch
     line input #78: b$                 ! get old scratch value
     if a$=b$ then
       print #77: "0 0 0"               ! put transparent color into file
     else
       print #77: a$                    ! put new screen value into file
     end if
   next x
  next y
  set window xmin,xmax,ymin,ymax  ! restore window coordinates
  close #77
  close #78
  close #79
  print #80 : ">FRAME NOTHING" ! animator script file
  print #80 : "FILE "; progname$&str$(ppmnumb)&".gif"
  print #80 : "TRANS"
  print #80 : "TIME"
end sub

sub closeppm  ! clean up temporary files
  when error in
   close #80
   unsave "temp0.tmp"
   unsave "temp1.tmp"
  use
  end when
end sub

sub saveppm
 if sketch=1 then exit sub
 if ppmnumb=0 then
   call initppm
 else
   call moreppm  
 end if
end sub

sub saveppmx
  ! saves cur win to pbmplus file  progname$.ppm, 480x640 1 minute 2.2MB
  local x,y,a$,px,py,xmin,xmax,ymin,ymax,point
  ask window xmin,xmax,ymin,ymax    ! save window coordinates
  ask pixels px,py
  set window 1,px,1,py
  call getcolormix
  open #77 : name ppm$&".ppm", create "newold"
  erase #77
  print #77: "P3",px,py,3   ! P3=ppm, width, height, max rgb 0-3
  for y= py to 1 step -1
   for x= 1 to px
     box keep x,x,y,y in a$
     let point = min(1,ord(a$[9:9]))+ 2*min(1,ord(a$[10:10]))
     let point = point + 4*min(1,ord(a$[11:11]))+8*min(1,ord(a$[12:12]))
     print #77: colormix$(point)
   next x
  next y
  set window xmin,xmax,ymin,ymax  ! restore window coordinates
  close #77
end sub



!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!! main program !!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

let dshift=2*pi/8
let shift=-dshift
let ppmnumb=0
do
 let shift=shift+dshift
 call redimension(4000)
 let objq=0
 call loadgrid                       ! load in grid
 let maxq=objq-1
 call loadindex(index,maxq)      !  must be done for paint even if no sort
 call sort(index,0,maxq,scratch) !  sort
 call init                       ! initiate graphics
 call paint(0,maxq)              ! do the drawing
 call saveppm
loop until shift>=2*pi
call closeppm ! clean up
set mode "history"                  ! exit program in text mode
! sound 600,1
end   !!! end of main program