// File: SPLINE3.BCPL
// P.Baudelaire
// December 5, 1977  4:49 PM

// The procedure ParametricSpline implements the algorithms
// for local cubic splines described in
//		"Spline Curve Techniques"
//		by P.Baudelaire, R.Flegal, & R.Sproull
//		Xerox Internal Report  (1977)

// Uses MICROCODED floating point routines

// outgoing procedures:

external [
	ParametricSpline
	PSerror
	]

// outgoing statics:

external [
	PSzone
	]

static [
	PSzone=0		// storage zone
	]

// incoming procedures:

external [
	FLD; FAD; FML; FST	// FLOAT (Alto floating point package)
	FLDI; FSB; FDV
	FCM; FNEG; FTR
	FPSetup

	Allocate		// Alto SYSTEM
	Free
	MoveBlock
	Zero
	]

// incoming statics:

external [
	FPwork			// FLOAT  (floating point accumulator)
	]

// local definitions:

manifest [
	naturalSpline=0
	periodicSpline=1
	// floating point registers: 0 to 4
	ac0=0; ac1=1; ac2=2; ac3=3; ac4=4
	// constants:
	zero=5; one=6; two=7; six=8
	numFPacs=9
	// local spline type
	CatmullRom=1
	FourPoints=2
	Bspline1=3
	Bspline2=4
	Bspline3=5
	Hermit=6
	]

structure PSVEC [
	FPworkSave word
	FPworkNew word
	fpx word
	fpy word
	a word
	]

manifest lPSVEC=size PSVEC/16

// local statics:

static [
	// working storage
	PSvec=0
	// various parameters
	defaultPStype=Bspline3
	PScoefNums=0
	PScoefDens=0
	PScoefEndNum=1
	PScoefEndDen=1
	HermitCoef=3
	]




let ParametricSpline(n, px, py, p1x, p2x, p3x, p1y, p2y, p3y, endCondition, PStype;
numargs nargs) = valof [

	// default arguments, get storage, check various things
	let tempVec= vec lPSVEC
	if PSinit(tempVec) eq 0 resultis 0


	PScoefNums= table [ 1; 1; 1; 1; 1; 1 ]
	PScoefDens= table [ 2; 6; 1; 2; 6; 2 ]
	let convertXY= n ls 0
	if convertXY then n=-n
	let x=PSallocate(lv(PSvec>>PSVEC.fpx), 2*n+4)
	let y=PSallocate(lv(PSvec>>PSVEC.fpy), 2*n+4)
	if (x eq 0) % (y eq 0) resultis 0

	test convertXY
	ifso [
		// convert coordinates from integer to floating point
		for i=0 to n-1 do [
			FST(FLDI(ac1, x!i), px+2*(i+1))
			FST(FLDI(ac1, y!i), py+2*(i+1))
			]
		]
	ifnot [
		MoveBlock(x+2, px, 2*n)
		MoveBlock(y+2, py, 2*n)
		]

	switchon nargs into [
		case 9:
			endCondition=naturalSpline
		case 10:
			PStype=defaultPStype
		case 11:
			if endCondition ne naturalSpline &
			   endCondition ne periodicSpline resultis PSquit(PSerror(3))
			if n ls 3 then endCondition=naturalSpline
			if PStype eq 0 then PStype=defaultPStype
			endcase
		default:
			resultis PSquit(PSerror(4))
		]

	test endCondition eq periodicSpline
	ifso [
		if ((FCM(FLD(ac1, x+2), x+2*n) ne 0) %
		   (FCM(FLD(ac2, y+2), y+2*n) ne 0)) resultis PSquit(PSerror(2))
		MoveBlock(x, x+2*(n-1), 2)
		MoveBlock(y, y+2*(n-1), 2)
		MoveBlock(x+2*(n+1), x+4, 2)
		MoveBlock(y+2*(n+1), y+4, 2)
		]
	ifnot [
		// ac3=b=PScoefEndNum/PScoefEndDen
		FDV(FLDI(ac3, PScoefEndNum), FLDI(ac2, PScoefEndDen))
		// x(0)=x(1)+b*(x(1)-x(2))
		FST(FAD(FML(FSB(FLD(ac1, x+2), x+4), ac3), x+2), x)
		FST(FAD(FML(FSB(FLD(ac1, y+2), y+4), ac3), y+2), y)
		// x(n+1)=x(n)+b*(x(n)-x(n-1))
		FST(FAD(FML(FSB(FLD(ac1, x+2*n), x+2*(n-1)), ac3), x+2*n), x+2*(n+1))
		FST(FAD(FML(FSB(FLD(ac1, y+2*n), y+2*(n-1)), ac3), y+2*n), y+2*(n+1))
		]

	let a=PSallocate(lv(PSvec>>PSVEC.a), 32)
	if (a eq 0) resultis 0
	Zero(a, 32)
	let c=nil
	let saveXY=false
	FDV(FLDI(ac2, PScoefNums!(PStype-1)), FLDI(ac3, PScoefDens!(PStype-1)))
	test PStype eq CatmullRom
	ifnot [
		switchon PStype into [
		case FourPoints:
			c=table [
				-1;  3; -3;  1;
				 3; -6;  3;  0;
				-2; -3;  6; -1;
				 0;  6;  0;  0 ]
			endcase
		case Bspline1:
			c=table [
				 0;  0;  0;  0;
				 0;  0;  0;  0;
				 0; -1;  1;  0;
				 0;  1;  0;  0 ]
			endcase
		case Bspline2:
			c=table [
				 0;  0;  0;  0;
				 1; -2;  1;  0;
				-2;  2;  0;  0;
				 1;  1;  0;  0 ]
			saveXY=true
			endcase
		case Bspline3:
			c=table [
				-1;  3; -3;  1;
				 3; -6;  3;  0;
				-3;  0;  3;  0;
				 1;  4;  1;  0 ]
			saveXY=true
			endcase
		case Hermit:
			c=table [
				 0;  0;  0;  0;
				 0;  0;  0;  0;
				 0;  0;  0;  0;
				 1;  1;  0;  0 ]
			c!0=2-HermitCoef
			c!2=HermitCoef-2
			c!4=2*HermitCoef-3
			c!5=-HermitCoef
			c!6=3-HermitCoef
			c!8=-HermitCoef
			c!9=HermitCoef
			saveXY=true
			endcase
			]
		for i=0 to 15 do [
			FST(FML(FLDI(ac1, c!i), ac2), a+2*i)
			]
		]
	ifso [
		// case CatmullRom:
		FST(FST(ac2, a+6), a+20)
		FST(FST(FST(FNEG(FLD(ac1, ac2)), a), a+14), a+16)
		FST(FNEG(FST(FSB(FLDI(ac1, 2), ac2), a+2)), a+4)
		FST(FNEG(FSB(FLDI(ac1, 3), ac2)), a+10)
		FST(FSB(FLDI(ac1, 3), FST(FAD(ac2, ac2), a+8)), a+12)
		FST(one, a+26)
		]
	for i=0 to 3 do [
		FST(FML(FLD(ac1, a+2*i), six), a+2*i)
		FST(FML(FLD(ac1, a+2*(i+4)), two), a+2*(i+4))
		]

	//do everything twice to get x(t) and y(t).
	let p1,p2,p3,p,s=nil,nil,nil,nil,nil
	for t=1 to 2 do [
		test (t eq 1)
		   ifso [ s=px; p=x; p1=p1x; p2=p2x; p3=p3x ]
		   ifnot [ s=py; p=y; p1=p1y; p2=p2y; p3=p3y ]
		for i=0 to n-1 do [
			FLDI(ac1, 0); FLDI(ac2, 0)
			FLDI(ac3, 0); FLDI(ac4, 0) 
			for j=0 to 3 do [
				FAD(ac4, FML(FLD(ac0, a+2*j), p+2*j))
				FAD(ac3, FML(FLD(ac0, a+2*(j+4)), p+2*j))
				FAD(ac2, FML(FLD(ac0, a+2*(j+8)), p+2*j))
				FAD(ac1, FML(FLD(ac0, a+2*(j+12)), p+2*j))
				] 
			FST(ac4, p3+2*i)
			FST(ac3, p2+2*i)
			FST(ac2, p1+2*i)
			if saveXY then test convertXY
				ifso s!i=FTR(ac1)
				ifnot FST(ac1, s+2*i) 
			p=p+2
			]
		]
	resultis PSquit(true)
	]


and PSinit(psv) = valof [
	if PSzone eq 0 resultis PSerror(0)
	PSvec=psv
	Zero(PSvec, lPSVEC)
	// new floating point work area
	let FPworkNew= Allocate(PSzone, 4*numFPacs+1)
	if FPworkNew eq 0 resultis PSerror(1)
	PSvec>>PSVEC.FPworkSave=FPwork
	PSvec>>PSVEC.FPworkNew=FPworkNew
	FPworkNew!0=numFPacs
	FPSetup(FPworkNew)
	FLDI(zero,0); FLDI(one, 1); FLDI(two,2); FLDI(six,6)
	resultis true
	]



and PSallocate(location, m) = valof [
	let b=Allocate(PSzone, m)
	if b eq 0 resultis PSquit(PSerror(1))
	@location=b
	resultis b
	]



and PSerror(errorCode,a1,a2,a3,a4) = valof [
	(table[#77403; #1401]) ("PS.ERRORS", lv errorCode)
	resultis 0
	]


and PSquit(result) = valof [
	if PSvec eq 0 resultis result
	FPSetup(PSvec>>PSVEC.FPworkSave)
	PSvec>>PSVEC.FPworkSave=0
	for i=0 to lPSVEC-1 do if PSvec!i ne 0 then Free(PSzone, PSvec!i)
	PSvec=0
	resultis result
	]