// March 31, 1978  2:02 PM				*** RESIDENT ***

get "zpDefs.bcpl"


// outgoing procedures:

external [
	newSplineID
	splineType
	makeSpline
	remakeSpline
	computeSpline
	]
  

// outgoing statics:

external [
	@computeRange
	@numRate
	@denRate
	@pseudoCyclic
	@splineTension
	]

static [
	// spline computation parameters
	@computeRange=64
	@numRate=3
	@denRate=2
	// flag for pseudo cyclic hack vs periodic splines
	@pseudoCyclic=false
	@splineTension=false
	]


// incoming procedures:

external [
	MoveBlock		// SYSTEM
	SetBlock
	Zero
	Gets

	typeForm		// ZPUTIL
	giveUp

	ParametricSpline	// PSPLINE

	FLD; FST; FTR; FLDI	// FLOAT
	FNEG; FAD; FSB; FML
	FDV; FCM; FSN
	FSTDP; FLDDP
	DPAD; DPSB

	obtainBlock		// ZPBLOCK
	putBlock

	checkPoint		// ZPDRAW
	startCode
	endCode
	checkSplineID
	curve
	splineColorSymbol
	giveMeXY
	]


// incoming statics:

external [
	keys			// SYSTEM

	@splineTable		// ZPINIT
	@maxSplineID

	// globals for encoding:
	@newX			// ZPDRAW 
	@newY
	]

// local statics:

static [
	@debugSpline=0
	]

// local definitions





//*****************************************************************
// First: make spline block
//*****************************************************************

let makeSpline(n, xTable, yTable, brush, color, cyclic; numargs m) = valof [makeSpline
	// xTable & yTable are FLOATING POINT
	//returns id of new spline, 0 if error
	unless n gr 0 resultis 0
	if m ls 6 then cyclic=false
	// check for false curves that actually are lines or dots
	let minVal=vec 2
	let maxVal=vec 2
	if n gr 2 & TableConstant(n, yTable) then [
		TableMax(n, xTable, maxVal)
		TableMin(n, xTable, minVal)
		n=2
		FST(FLD(0, minVal), xTable)
		FST(FLD(0, maxVal), xTable+2)
		]
	if n gr 2 & TableConstant(n, xTable) then [
		TableMax(n, yTable, maxVal)
		TableMin(n, yTable, minVal)
		n=2
		FST(FLD(0, minVal), yTable)
		FST(FLD(0, maxVal), yTable+2)
		]
	if n eq 2 & TableConstant(2, xTable) & TableConstant(2, yTable) then n=1
	// now make the spline
	let splinePointer=obtainBlock(SPLINEknotBase+4*n)
	unless splinePointer resultis giveUp("[makeSpline]")
	Zero(splinePointer, SPLINEknotBase+4*n)
	//make spline block
	let knotTable=splinePointer+SPLINEknotBase
	let dn=2*n
	MoveBlock(knotTable, xTable, dn)
	MoveBlock(knotTable+dn, yTable, dn)
	splinePointer>>SPLINE.nKnots=n
	splinePointer>>SPLINE.brush=brush
	splinePointer>>SPLINE.color=color
	splinePointer>>SPLINE.cyclic=cyclic
	splineType(splinePointer)
	giveMeXY(splinePointer, lv splinePointer>>SPLINE.xColor,
			lv splinePointer>>SPLINE.yColor)

	let id=remakeSpline(splinePointer)
	unless id then putBlock(splinePointer)
	resultis id
	]makeSpline



and splineType(splinePointer) be [splineType
	let n=splinePointer>>SPLINE.nKnots
	let xTable=splinePointer+SPLINEknotBase
	let yTable=xTable+2*n
	let x1=FTR(FLD(0, xTable))
	let y1=FTR(FLD(0, yTable))
	let x2, y2=nil, nil
	test n eq 1
	ifso [ x2=x1; y2=y1 ]
	ifnot [
		x2=FTR(FLD(0, xTable+2))
		y2=FTR(FLD(0, yTable+2))
		]
	splinePointer>>SPLINE.left= (x1 ls x2) ? x1, x2 
	splinePointer>>SPLINE.right= (x1 gr x2) ? x1, x2 
	splinePointer>>SPLINE.top= (y1 gr y2) ? y1, y2 
	splinePointer>>SPLINE.bottom= (y1 ls y2) ? y1, y2 
	splinePointer>>SPLINE.type = selecton n  into [
	case 1: dotSpline;
	case 2: ((y1 eq y2) ? horSpline, ((x1 eq x2) ? verSpline, regSpline));
	default: regSpline ]
	]splineType



//*****************************************************************
// Second: get spline id
//*****************************************************************


and remakeSpline(splinePointer) = valof [remakeSpline
	let id=newSplineID()
	unless id resultis 0
	test computeSpline(splinePointer)
	ifso [
		splineTable!0=splineTable!0+1
		splineTable!id=splinePointer
		resultis id
		]
	ifnot resultis 0
	]remakeSpline



and newSplineID() = valof [newSplineID
	for id=1 to maxSplineID do unless splineTable!id resultis id
	typeForm(0, "Sorry, no room for more than ", 10, maxSplineID, 0, "curves*N", 
	0, "To get a larger work space for curves, start DRAW with switch /S (e.g.: DRAW ", 10, 2*maxSplineID, 0, "/S )*N")
	resultis 0
	]newSplineID




//*****************************************************************
// Third: compute [ & generate chain encoding ]
//*****************************************************************



and computeSpline(splinePointer) = valof [
	let r=((splinePointer>>SPLINE.type eq regSpline) ?
		computeRegularSpline(splinePointer), computeFalseSpline(splinePointer))
	splineColorSymbol(splinePointer)
	resultis r
	]



and computeFalseSpline(splinePointer) = valof [computeFalseSpline
	// i.e. single point, or horizontal or vertical line
	splinePointer>>SPLINE.chain=0
	splinePointer>>SPLINE.nBeads=0
	splinePointer>>SPLINE.xStart=splinePointer>>SPLINE.left
	splinePointer>>SPLINE.yStart=splinePointer>>SPLINE.top
	curve(splinePointer, drawMode)
	resultis true
	]computeFalseSpline



and computeRegularSpline(splinePointer) = valof [computeRegularSpline
	// returns 0 if no storage

	let procedureName="[computeRegularSpline]"
	let n=splinePointer>>SPLINE.nKnots
	let dn=2*n
	// 4 extra knots for pseudo-cyclic curve
	let xn=(splinePointer>>SPLINE.cyclic) & pseudoCyclic ? 4, 0
	let nt=n+xn
	let dnt=2*nt

	// table for derivatives
	let dTable=obtainBlock(nt*12)
	unless dTable resultis giveUp(procedureName)
	let knotTable=splinePointer+SPLINEknotBase
	let kTable=obtainBlock(2*dnt)
	unless kTable resultis giveUp(procedureName, dTable)
	let xTable=kTable+xn
	let yTable=xTable+dnt
	MoveBlock(xTable, knotTable, dn)
	MoveBlock(yTable, knotTable+dn, dn)
	if xn ne 0 then [
		// extra knots
		MoveBlock(kTable, xTable+dn-6, 4)
		MoveBlock(xTable+dn, xTable+2, 4)
		MoveBlock(yTable-4, yTable+dn-6, 4)
		MoveBlock(yTable+dn, yTable+2, 4)
		]

	manifest [
		//local floating Point registers
		t0=0; t1=1; t2=2; t3=3
		//caution: global floating Point registers
		X=8; Y=9
		d1X=10; d2X=11; d3X=12
		d1Y=13; d2Y=14; d3Y=15
		]

	// get storage for chain encoding
	// (try to estimate block size...)
	let s=1
	let kx=FTR(FLD(t0,xTable))
	let ky=FTR(FLD(t1,yTable))
	for i=1 to n-1 do [
		let kx1=FTR(FLD(t0,xTable+2*i))
		let ky1=FTR(FLD(t1,yTable+2*i))
		s=s+1+max(abs(kx-kx1),abs(ky-ky1))/16
		kx,ky=kx1,ky1
		]
	let nBeads=(((s*4)/3)/BEADsize)+1
	let maxByteCount=16*nBeads
	let chainSize=nBeads*(BEADsize+2)
	let chainBlock=obtainBlock(chainSize+maxByteCount/2)
	unless chainBlock resultis giveUp(procedureName, dTable, (xn ? kTable, 0))

	let wTable=0
	if splineTension then [
		wTable=obtainBlock(nt)
		if wTable ne 0 then [
			SetBlock(wTable, 0, nt)
			typeForm(0,"Input spline weigths now @", 8, wTable, 1, $*N)
			Gets(keys)
			]
		]
	
	// spline computation:
	// caution: ParametricSpline(n,X,Y,X',X'',X''',Y',Y'',Y''')
	let psDone=ParametricSpline(nt, kTable, kTable+dnt, dTable, dTable+dnt, dTable+2*dnt, 
		dTable+3*dnt, dTable+4*dnt, dTable+5*dnt,
		((splinePointer>>SPLINE.cyclic & not pseudoCyclic) ? periodicSpline, naturalSpline),
		wTable)

	putBlock(wTable)
	unless psDone resultis giveUp(procedureName, dTable, chainBlock, (xn ? kTable, 0))

	newX=FTR(FLD(t0,xTable))
	newY=FTR(FLD(t1,yTable))
	splinePointer>>SPLINE.chain=chainBlock
	chainBlock!(chainSize-1)=maxByteCount
	splinePointer>>SPLINE.nBeads=nBeads
	splinePointer>>SPLINE.xStart=newX
	splinePointer>>SPLINE.yStart=newY
	startCode(splinePointer,newX,newY)

	let two=vec 2;   FST(FLDI(t0,2),two)
	let three=vec 2; FST(FLDI(t0,3),three)
	let six=vec 2;   FST(FLDI(t0,6),six)

	for k=0 to n-2 do [
	  //knot k
	  let kX,kY=xTable+2*k,yTable+2*k
	  //stepping parameters
	  let m=max(abs(FTR(FLD(t0,kX))-FTR(FLD(t1,kX+2))),
				abs(FTR(FLD(t2,kY))-FTR(FLD(t3,kY+2))))
	  m=(m*numRate)/denRate
	  let ni,r=1,m
	  if m gr computeRange then [ r=computeRange; ni=m/r; m=ni*r ]
	  //constants
	  let fni=vec 2; FST(FLDI(t0,ni), fni)
	  test m eq 0
	  ifso FLDI(t1,0)
	  ifnot FDV(FLDI(t1,1), FLDI(t0,m))
	  let delta=vec 2;  FST(t1,delta)
	  let delta2=vec 2; FST(FML(t1,delta), delta2)
	  let delta3=vec 2; FST(FML(t1,delta), delta3)
	  //derivatives - floating point -
	  let kX1=dTable+xn+2*k
	  let kX2=kX1+dnt
	  let kX3=kX2+dnt
	  let kY1=kX3+dnt
	  let kY2=kY1+dnt
	  let kY3=kY2+dnt
	  // show values
	  if debugSpline then [
		typeForm(0, "m=", 10, m, 0, ", delta=", -2, delta, 0, ",  X, X', X'', X''' =*N", -2, kX)
		for i=0 to 2 do typeForm(1, $|, -2, kX1+dnt*i)
		typeForm(1, $*N, -2, kY)
		for i=0 to 2 do typeForm(1, $|, -2, kY1+dnt*i)
		Gets(keys)
		typeForm(1, $*N)
		]
	  //starting derivatives of subintervals - floating point -
	  let X1=vec 2; let X2=vec 2; let X3=vec 2
	  let Y1=vec 2; let Y2=vec 2; let Y3=vec 2
	  //start at knot k
	  FLD(X,kX);  FLD(Y,kY)
	  X1!0=kX1!0; X1!1=kX1!1
	  Y1!0=kY1!0; Y1!1=kY1!1
	  X2!0=kX2!0; X2!1=kX2!1
	  Y2!0=kY2!0; Y2!1=kY2!1
	  X3!0=kX3!0; X3!1=kX3!1
	  Y3!0=kY3!0; Y3!1=kY3!1
	  //3rd derivatives & differences are constant
	  FML(FLD(d3X,delta3), X3)
	  FML(FLD(d3Y,delta3), Y3)
	  //forward difference computation in subinterval
	  for i=1 to ni do [
	    //floating point computation of initial values (blahh!)
	    FML(FLD(d2X,delta2),X2)
	    FML(FLD(d1X,delta),X1)
	    FAD(d1X,FDV(FLD(t0,d2X),two))
	    FAD(d1X,FDV(FLD(t0,d3X),six))
	    FML(FLD(d2Y,delta2),Y2)
	    FML(FLD(d1Y,delta),Y1)
	    FAD(d1Y,FDV(FLD(t0,d2Y),two))
	    FAD(d1Y,FDV(FLD(t0,d3Y),six))
	    FAD(d2X,d3X);    FAD(d2Y,d3Y)
	    //double precision fixed point variables
	    let dpX=vec 2; let dpd1X=vec 2; let dpd2X=vec 2; let dpd3X=vec 2
	    let dpY=vec 2; let dpd1Y=vec 2; let dpd2Y=vec 2; let dpd3Y=vec 2
	    //initial values
	    FSTDP(X,dpX); FSTDP(d1X,dpd1X); FSTDP(d2X,dpd2X); FSTDP(d3X,dpd3X)
	    FSTDP(Y,dpY); FSTDP(d1Y,dpd1Y); FSTDP(d2Y,dpd2Y); FSTDP(d3Y,dpd3Y)
	    //compute
	    for j=1 to r do [
		newX=DPAD(dpX,dpd1X)
		newY=DPAD(dpY,dpd1Y)
		DPAD(dpd1X,dpd2X); DPAD(dpd1Y,dpd2Y)
		DPAD(dpd2X,dpd3X); DPAD(dpd2Y,dpd3Y)
		checkPoint()
		]
	    //now, next interval or next knot
	    test i eq ni
	    ifso [
		//next knot
		newX=FTR(FLD(t0,kX+2))
		newY=FTR(FLD(t0,kY+2))
		]
	    ifnot [
		//starting point of next interval
		FDV(FLDI(t1,i),fni)
		FDV(FLD(t2,t1),two)
		FDV(FLD(t3,t1),three)
		FST(FAD(FML(FLD(t0,kX3),t1),kX2),X2)
		FST(FAD(FML(FLD(t0,kY3),t1),kY2),Y2)
		FST(FAD(FML(FAD(FML(FLD(t0,kX3),t2),kX2),t1),kX1),X1)
		FST(FAD(FML(FAD(FML(FLD(t0,kY3),t2),kY2),t1),kY1),Y1)
		newX=FTR(FAD(FML(FAD(FML(FAD(FML(FLD(X,kX3),t3),kX2),t2),kX1),t1),kX))
		newY=FTR(FAD(FML(FAD(FML(FAD(FML(FLD(Y,kY3),t3),kY2),t2),kY1),t1),kY))
 		] 
	    checkPoint()
	    ]
	  ]

	//end of chain encoding
	endCode(splinePointer)

	//return storage [ chain block is trimmed by endCode() ]
	putBlock(dTable)
	putBlock(kTable)

	resultis true
	]computeRegularSpline


//*****************************************************************
// miscellaneous useful stuff
//*****************************************************************

and abs(i) = ((i gr 0) ? i, -i)

and max(i1,i2) = ((i1 ge i2) ? i1, i2)

and TableMax(n, fpTable, fpMax) be [
	manifest [ t0=0 ]
	FST(FLD(t0, fpTable), fpMax)
	for k=1 to n-1 do if FCM(t0, fpTable+2*k) eq 1 then
				FST(FLD(t0, fpTable+2*k), fpMax)
	]

and TableMin(n, fpTable, fpMin) be [
	manifest [ t0=0 ]
	FST(FLD(t0, fpTable), fpMin)
	for k=1 to n-1 do if FCM(t0, fpTable+2*k) eq -1 then
				FST(FLD(t0, fpTable+2*k), fpMin)
	]

and TableConstant(n, fpTable) = valof [
	manifest [ t0=0 ]
	FLD(t0, fpTable)
	for k=1 to n-1 do if FCM(t0, fpTable+2*k) ne 0 then resultis false
	resultis true
	]