DECLARE SUB SetPalette ()
DECLARE SUB PSET3D (xb%, yb%, zb%, c%)
DECLARE SUB InitControlPoints (surf() AS ANY)
DECLARE SUB curve (x00%, y00%, z00%, x01%, y01%, z01%, x02%, y02%, z02%, x03%, y03%, z03%, n%, t%)
DECLARE SUB CurveFDL2 (x00%, y00%, z00%, x01%, y01%, z01%, x02%, y02%, z02%, x03%, y03%, z03%, n%, t%)
DECLARE SUB CurveHorner (x00%, y00%, z00%, x01%, y01%, z01%, x02%, y02%, z02%, x03%, y03%, z03%, n%, tp%)
'------------------------------------------------------------------
'BEZSURF.BAS - bezier/catmull-rom surface drawing algorithm
'         demonstration by Toshi Horie (version 2.2)
' This one uses the 3 dimensional forward differencing, so it can
' handle any kind of 3D bezier surface.  It uses level of detail
' management for faster rendering.  This can be sped up immensely
' using a fast perspective projection routine instead of PSET3D.
' Copyright (c) 2000 Toshihiro Horie
' toshiman@uclink4.berkeley.edu
'------------------------------------------------------------------

DEFINT A-Z
COMMON SHARED detail
CONST FALSE = 0, TRUE = NOT FALSE

TYPE rgbtype
	red AS INTEGER
	green AS INTEGER
	blue AS INTEGER
END TYPE

TYPE Pt
	X AS INTEGER
	Y AS INTEGER
	Z AS INTEGER
END TYPE

'catmull-rom curve coefficients
CONST cc = .8:  ' c=[0-1]   c = 1-tightness
CONST tt1 = -cc
CONST tt2 = 2 - cc
CONST tt3 = cc - 2
CONST tt4 = cc
CONST tt5 = 2 * cc
CONST tt6 = cc - 3
CONST tt7 = 3 - 2 * cc
CONST tt8 = -cc
'------------------------------------------------
'These constants select the appearance and
'speed of the curve drawing algorithms.

  CONST BEZIER = TRUE
  CONST FAST = TRUE
  CONST MAXDETAIL = 60
'-------------------------------------------------
DIM surf(16, 16) AS Pt
DIM SHARED ctrl(MAXDETAIL, 4) AS Pt


SCREEN 13
CALL SetPalette
CALL InitControlPoints(surf())

' +00-----10-----20-----30-> s
' | *....`*`````*`....`*
' |
' |
' 01*``...*...``*````..*31
' |
' |
' 02*....`*.````*``````*32
' |
' |
' 03*`````*``...*.....`*33
' |
' v
' t
' STEP 1) Using the 16 control points that define the patch
'         (they are marked with the * symbol), we fill in the
'         details in the s parameter direction *only at
'         the four rows* that come from the control points.
' Note: These details will become the control points for
'       filling in surface patch details in the t direction.
'
'
' STEP 2) fill in the patch in the t direction.  The diagram shows
'         the filling of the curve at s=0.2 (close to s=0)
'         We use the detailed control points we calculated
'         in step 1 (marked with @).
' +00-----10-----20-----30-> s
' | *@...`*`````*`....`*
' |   #
' |    #
' 01*``@..*...``*````..*31
' |   #
' |   #
' 02*@..`.*.````*``````*32
' | #
' | #
' 03*@````*``...*.....`*33
' |
' v
' t

'plots 16+ patches using bezier interpolation.
'When step is 3, we get C0 continuity.
'the continuity of the patches depends on the positions of
'the control points in the surf() array.
timer1! = TIMER
IF BEZIER = TRUE THEN st = 3 ELSE st = 1
FOR si = 12 TO 0 STEP -st    'loop thru terrain patches in s direction
	detail = MAXDETAIL - si * 4
	FOR ti = 0 TO 12 STEP st 'loop thru terrain patches in t direction
		' fetch the control points from memory to registers (if they fit)
		x00% = surf(si, ti).X: x10% = surf(si + 1, ti).X: x20% = surf(si + 2, ti).X: x30% = surf(si + 3, ti).X
		y00% = surf(si, ti).Y: y10% = surf(si + 1, ti).Y: y20% = surf(si + 2, ti).Y: y30% = surf(si + 3, ti).Y
		z00% = surf(si, ti).Z: z10% = surf(si + 1, ti).Z: z20% = surf(si + 2, ti).Z: z30% = surf(si + 3, ti).Z
		'1) fill in details in s direction (t=0)
		CALL curve(x00%, y00%, z00%, x10%, y10%, z10%, x20%, y20%, z20%, x30%, y30%, z30%, detail, 0)    '1st row
	   
		x01% = surf(si, ti + 1).X: x11% = surf(si + 1, ti + 1).X: x21% = surf(si + 2, ti + 1).X: x31% = surf(si + 3, ti + 1).X
		y01% = surf(si, ti + 1).Y: y11% = surf(si + 1, ti + 1).Y: y21% = surf(si + 2, ti + 1).Y: y31% = surf(si + 3, ti + 1).Y
		z01% = surf(si, ti + 1).Z: z11% = surf(si + 1, ti + 1).Z: z21% = surf(si + 2, ti + 1).Z: z31% = surf(si + 3, ti + 1).Z
		'1) fill in details in s direction (t=1)
		CALL curve(x01%, y01%, z01%, x11%, y11%, z11%, x21%, y21%, z21%, x31%, y31%, z31%, detail, 1)    '2nd row
	   
		x02% = surf(si, ti + 2).X: x12% = surf(si + 1, ti + 2).X: x22% = surf(si + 2, ti + 2).X: x32% = surf(si + 3, ti + 2).X
		y02% = surf(si, ti + 2).Y: y12% = surf(si + 1, ti + 2).Y: y22% = surf(si + 2, ti + 2).Y: y32% = surf(si + 3, ti + 2).Y
		z02% = surf(si, ti + 2).Z: z12% = surf(si + 1, ti + 2).Z: z22% = surf(si + 2, ti + 2).Z: z32% = surf(si + 3, ti + 2).Z
		'1) fill in details in s direction (t=2)
		CALL curve(x02%, y02%, z02%, x12%, y12%, z12%, x22%, y22%, z22%, x32%, y32%, z32%, detail, 2)    '3rd row
	
		x03% = surf(si, ti + 3).X: x13% = surf(si + 1, ti + 3).X: x23% = surf(si + 2, ti + 3).X: x33% = surf(si + 3, ti + 3).X
		y03% = surf(si, ti + 3).Y: y13% = surf(si + 1, ti + 3).Y: y23% = surf(si + 2, ti + 3).Y: y33% = surf(si + 3, ti + 3).Y
		z03% = surf(si, ti + 3).Z: z13% = surf(si + 1, ti + 3).Z: z23% = surf(si + 2, ti + 3).Z: z33% = surf(si + 3, ti + 3).Z
		'1) fill in details in s direction (t=3)
		CALL curve(x03%, y03%, z03%, x13%, y13%, z13%, x23%, y23%, z23%, x33%, y33%, z33%, detail, 3)    '4th row

		'2) fill in detail in t direction and plot surface
		FOR i% = 0 TO detail
			xi0% = ctrl(i%, 0).X: yi0% = ctrl(i%, 0).Y: zi0% = ctrl(i%, 0).Z: 'row 0 (1st control point)
			xi1% = ctrl(i%, 1).X: yi1% = ctrl(i%, 1).Y: zi1% = ctrl(i%, 1).Z: 'row 1 (2nd control point)
			xi2% = ctrl(i%, 2).X: yi2% = ctrl(i%, 2).Y: zi2% = ctrl(i%, 2).Z: 'row 2 (3rd control point)
			xi3% = ctrl(i%, 3).X: yi3% = ctrl(i%, 3).Y: zi3% = ctrl(i%, 3).Z: 'row 3 (4th control point)
			'-1 = don't fill detailed control point array
			CALL curve(xi0%, yi0%, zi0%, xi1%, yi1%, zi1%, xi2%, yi2%, zi2%, xi3%, yi3%, zi3%, detail, -1)
		NEXT i%
	NEXT
NEXT
PRINT USING "##.## secs"; TIMER - timer1!
k$ = INPUT$(1)
END

SUB curve (x00%, y00%, z00%, x01%, y01%, z01%, x02%, y02%, z02%, x03%, y03%, z03%, n%, t%)
IF FAST THEN
	CALL CurveFDL2(x00%, y00%, z00%, x01%, y01%, z01%, x02%, y02%, z02%, x03%, y03%, z03%, n%, t%)
ELSE
	CALL CurveHorner(x00%, y00%, z00%, x01%, y01%, z01%, x02%, y02%, z02%, x03%, y03%, z03%, n%, t%)
END IF
END SUB

DEFLNG A-Z
SUB CurveFDL2 (x00%, y00%, z00%, x01%, y01%, z01%, x02%, y02%, z02%, x03%, y03%, z03%, n%, t%)
'params:
'  (X00,Y00,Z00) is the first control point
'  (X03,Y03,Z03) is the last control point
'  n    is the number of steps to take only 20<=n%<=64 work reliably.
'  t% is when stepping (filling in details) in  s direction, it's used as
'          the control point array index.
'----------------------------------------------------------
'* Note that this diverges from the real curve over time
'  because of accumulating roundoff error.
'* I did multiple scalings to increase accuracy
'  by minimizing the amount of sig. figs lost to truncation,
'  at the cost of four extra divisions at startup.

IF BEZIER THEN
'x coefficients
  A1 = (x03% - 3 * x02% + 3 * x01% - x00%)
  B1 = (3 * x02% - 6 * x01% + 3 * x00%)
  C1 = (3 * x01% - 3 * x00%)
  D1 = x00%

'y coefficients
  A2 = (y03% - 3 * y02% + 3 * y01% - y00%)
  B2 = (3 * y02% - 6 * y01% + 3 * y00%)
  C2 = (3 * y01% - 3 * y00%)
  D2 = y00%

'z coefficients
  A3 = (z03% - 3 * z02% + 3 * z01% - z00%)
  B3 = (3 * z02% - 6 * z01% + 3 * z00%)
  C3 = (3 * z01% - 3 * z00%)
  D3 = z00%
ELSE 'catmull-rom
 'x coefficients
  A1 = (tt1 * x00% + tt2 * x01% + tt3 * x02% + tt4 * x03%)
  B1 = (tt5 * x00% + tt6 * x01% + tt7 * x02% + tt8 * x03%)
  C1 = (-cc * x00% + cc * x02%)
  D1 = x01%
  
 'y coefficients
  A2 = (tt1 * y00% + tt2 * y01% + tt3 * y02% + tt4 * y03%)
  B2 = (tt5 * y00% + tt6 * y01% + tt7 * y02% + tt8 * y03%)
  C2 = (-cc * y00% + cc * y02%)
  D2 = y01%
  
 'z coefficients
  A3 = (tt1 * z00% + tt2 * z01% + tt3 * z02% + tt4 * z03%)
  B3 = (tt5 * z00% + tt6 * z01% + tt7 * z02% + tt8 * z03%)
  C3 = (-cc * z00% + cc * z02%)
  D3 = z01%
END IF
'note: f=delta
f = 9800& \ n%
fsqr = f * f \ 98&       'used to be \9800& but now split into \98& and \100& for more accuracy
fcube = f * fsqr \ 100&  ' divide by 9800& when PSETing
'LOCATE 10, 1: PRINT f; fsqr; fcube

'initial conditions of (0th to 3rd) derivatives for x
   X = D1 * 9800&
  dx = A1 * fcube \ 9800& + B1 * fsqr \ 100& + C1 * f
dddx = 6 * A1 * fcube \ 9800&
 ddx = dddx + B1 * 2 * fsqr \ 100&

'initial conditions of (0th to 3rd) derivatives for y
   Y = D2 * 9800&
  dy = A2 * fcube \ 9800& + B2 * fsqr \ 100& + C2 * f
dddy = 6 * A2 * fcube \ 9800&
 ddy = dddy + B2 * 2 * fsqr \ 100&

'initial conditions of (0th to 3rd) derivatives for z
   Z = D3 * 9800&
  dz = A3 * fcube \ 9800& + B3 * fsqr \ 100& + C3 * f
dddz = 6 * A3 * fcube \ 9800&
 ddz = dddz + B3 * 2 * fsqr \ 100&

i% = 0
'draw the curve, baby!
DO
	SELECT CASE t%
		CASE 0, 3
			IF BEZIER THEN CALL PSET3D(X \ 9800, Y \ 9800, Z \ 9800, 15)
			ctrl(i%, t%).X = X \ 9800
			ctrl(i%, t%).Y = Y \ 9800
			ctrl(i%, t%).Z = Z \ 9800
		CASE 1, 2
			' don't plot the middle two control points
			' because the curve doesn't go through them
			ctrl(i%, t%).X = X \ 9800
			ctrl(i%, t%).Y = Y \ 9800
			ctrl(i%, t%).Z = Z \ 9800
		CASE ELSE
			CALL PSET3D(X \ 9800, Y \ 9800, Z \ 9800, Z \ 98000)
	END SELECT
	X = X + dx: dx = dx + ddx: ddx = ddx + dddx
	Y = Y + dy: dy = dy + ddy: ddy = ddy + dddy
	Z = Z + dz: dz = dz + ddz: ddz = ddz + dddz
	i% = i% + 1
LOOP WHILE i% <= n%

END SUB

DEFSNG A-Z
SUB CurveHorner (x00%, y00%, z00%, x01%, y01%, z01%, x02%, y02%, z02%, x03%, y03%, z03%, n%, t%)
IF BEZIER THEN
'x coefficients
  A1 = x03% - 3 * x02% + 3 * x01% - x00%
  B1 = 3 * x02% - 6 * x01% + 3 * x00%
  C1 = 3 * x01% - 3 * x00%
  D1 = x00%

'y coefficients
  A2 = y03% - 3 * y02% + 3 * y01% - y00%
  B2 = 3 * y02% - 6 * y01% + 3 * y00%
  C2 = 3 * y01% - 3 * y00%
  D2 = y00%

'z coefficients
  A3 = z03% - 3 * z02% + 3 * z01% - z00%
  B3 = 3 * z02% - 6 * z01% + 3 * z00%
  C3 = 3 * z01% - 3 * z00%
  D3 = z00%
ELSE 'catmull-rom
 'x coefficients
  A1 = (tt1 * x00% + tt2 * x01% + tt3 * x02% + tt4 * x03%)
  B1 = (tt5 * x00% + tt6 * x01% + tt7 * x02% + tt8 * x03%)
  C1 = (-cc * x00% + cc * x02%)
  D1 = x01%
  
 'y coefficients
  A2 = (tt1 * y00% + tt2 * y01% + tt3 * y02% + tt4 * y03%)
  B2 = (tt5 * y00% + tt6 * y01% + tt7 * y02% + tt8 * y03%)
  C2 = (-cc * y00% + cc * y02%)
  D2 = y01%
  
 'z coefficients
  A3 = (tt1 * z00% + tt2 * z01% + tt3 * z02% + tt4 * z03%)
  B3 = (tt5 * z00% + tt6 * z01% + tt7 * z02% + tt8 * z03%)
  C3 = (-cc * z00% + cc * z02%)
  D3 = z01%
END IF

f! = 1! / n%
xb% = x00%: yb% = y00%: zb% = z00%
i% = 0: tt! = 0
DO
	SELECT CASE t%
		CASE 0, 3
			IF BEZIER THEN CALL PSET3D(xb%, yb%, zb%, 15)
			ctrl(i%, t%).X = xb%
			ctrl(i%, t%).Y = yb%
			ctrl(i%, t%).Z = zb%
		CASE 1, 2
			' don't plot the middle two control points
			' because the curve doesn't go through them
			ctrl(i%, t%).X = xb%
			ctrl(i%, t%).Y = yb%
			ctrl(i%, t%).Z = zb%
		CASE ELSE
			CALL PSET3D(xb%, yb%, zb%, zb% \ 10)
	END SELECT
	tt! = tt! + f!
	i% = i% + 1
	xb% = ((A1 * tt! + B1) * tt! + C1) * tt! + D1
	yb% = ((A2 * tt! + B2) * tt! + C2) * tt! + D2
	zb% = ((A3 * tt! + B3) * tt! + C3) * tt! + D3
LOOP WHILE i% <= n%
END SUB

DEFINT A-Z
SUB InitControlPoints (surf() AS Pt)
'initializes the 16 control points in the 4x4 grid
'and return it through the surf() array
DIM wavy(128)
FOR a = 0 TO 16
	wavy(a) = COS(a * 17 / 3.14159265359#) * 2000
NEXT a
FOR ti = 0 TO 16
	FOR si = 0 TO 16
		surf(ti, si).X = ti * 20
		surf(ti, si).Y = si * 20
		' create sinusoidal heightfield of control points
		surf(ti, si).Z = wavy(ti + si)
	NEXT si
NEXT ti

END SUB

SUB PSET3D (xb%, yb%, zb%, c)
' the 3D to 2D projection should be dependent on
' the "center" and size of the surface so we can see the
' surface from a good viewpoint.
	 dist = xb% + 20
	 PSET (160 + 40 * (yb% - 160) \ dist, 100 - zb% \ dist), c
END SUB

SUB SetPalette
' ported from the C code in the book Computer Graphics by Foley et al.
numcolors = 256
max# = 63
DIM HSVPal(numcolors) AS rgbtype
SCREEN 13
s# = 0
s# = .9: 'saturation (0 to 1) 0=greyscale, 1=marker colors
v# = .9: 'brightness (0 to 1) 0=black, 1=brightest
p# = v# * (1# - s#)
conV# = numcolors / 360#
'if H was undefined, then we would have the greyscale.
FOR H% = 0 TO 359:       'hue (0 to 360) aka rainbow spectra
	c% = INT(H% * conV#): 'index into palette array
	Hn# = H% / 60        'Hn is in the interval [0,6)
	i% = INT(Hn#)        'i=floor(Hn)
	f# = Hn# - i%        'fractional part of Hn
	q# = v# * (1# - (s# * f#))
	t# = v# * (1# - (s# * (1 - f#)))
	'PRINT i%, c%, q#, t#
	SELECT CASE i%
		CASE 0:
			HSVPal(c%).red = INT(v# * max#)
			HSVPal(c%).green = INT(t# * max#)
			HSVPal(c%).blue = INT(p# * max#)
		CASE 1:
			HSVPal(c%).red = INT(q# * max#)
			HSVPal(c%).green = INT(v# * max#)
			HSVPal(c%).blue = INT(p# * max#)
		CASE 2:
			HSVPal(c%).red = p# * max#
			HSVPal(c%).green = v# * max#
			HSVPal(c%).blue = t# * max#
		CASE 3:
			HSVPal(c%).red = INT(p# * max#)
			HSVPal(c%).green = INT(q# * max#)
			HSVPal(c%).blue = INT(v# * max#)
		CASE 4:
			HSVPal(c%).red = INT(t# * max#)
			HSVPal(c%).green = INT(p# * max#)
			HSVPal(c%).blue = INT(v# * max#)
		CASE 5:
			HSVPal(c%).red = INT(v# * max#)
			HSVPal(c%).green = INT(p# * max#)
			HSVPal(c%).blue = INT(q# * max#)
	END SELECT
NEXT H%
OUT &H3C8, 0: OUT &H3C9, 0: OUT &H3C9, 0: OUT &H3C9, 0
FOR i% = 1 TO numcolors - 1
	OUT &H3C8, i%
	OUT &H3C9, HSVPal(i%).red
	OUT &H3C9, HSVPal(i%).green
	OUT &H3C9, HSVPal(i%).blue
NEXT i%
ERASE HSVPal
END SUB