// File: SPLINE1.BCPL
// P.Baudelaire & R.Flegal
// December 5, 1977  4:45 PM

// The procedure ParametricSpline implements the algorithm (1.2.7) described in
//		"Spline Curve Techniques"
//		by P.Baudelaire, R.Flegal, & R.Sproull
//		Xerox Internal Report  (May 1977)

//  Uses MICROCODE floating point routines

// outgoing procedures:

external [
	ParametricSpline
	PSerror
	]

// outgoing statics:

external [
	PSzone
	]

static [
	PSzone=0		// storage zone
	]

// incoming procedures:

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

	Allocate		// Alto SYSTEM
	Free
	Zero
	]

// incoming statics:

external [
	FPwork		// microFLOAT (Alto floating point package)
	]

// local definitions:

manifest [
	naturalSpline=0
	periodicSpline=1
	// floating point registers: 1 to 4
	ac1=1; ac2=2; ac3=3; ac4=4
	// constants:
	zero=5; one=6; two=7; six=8
	numFPacs=9
	]

structure PSVEC [
	FPworkSave word
	FPworkNew word
	fpx word
	fpy word
	a word
	b word
	c word
	r word
	s word
	]

manifest lPSVEC=size PSVEC/16

// local statics:

static [
	PSvec=0
	]



let ParametricSpline(n,x,y,p1x,p2x,p3x,p1y,p2y,p3y,splineType,w; numargs nargs) = valof [

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

	let p1,p2,p3,p=nil,nil,nil,nil
	let c,r,s=0,0,0

	if n ls 0 then [
		// convert coordinates from integer to floating point
		n=-n
		let fpx=PSallocate(lv(PSvec>>PSVEC.fpx), 2*n)
		let fpy=PSallocate(lv(PSvec>>PSVEC.fpy), 2*n)
		if (fpx eq 0) % (fpy eq 0) resultis 0
		for i=0 to n-1 do [
			FST(FLDI(ac1, x!i), fpx+2*i)
			FST(FLDI(ac1, y!i), fpy+2*i)
			]
		x=fpx
		y=fpy
		]

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

	if splineType eq periodicSpline then [
		if ((FCM(FLD(ac1,x), x+2*(n-1)) ne 0) %
		   (FCM(FLD(ac2,y), y+2*(n-1)) ne 0)) resultis PSquit(PSerror(2))
		c=PSallocate(lv(PSvec>>PSVEC.c), 2*n)
		r=PSallocate(lv(PSvec>>PSVEC.r), 2*n)
		s=PSallocate(lv(PSvec>>PSVEC.s), 2*n)
		if (c eq 0) % (r eq 0) % (s eq 0) resultis 0
		]

	let a=PSallocate(lv(PSvec>>PSVEC.a), 2*n)
	let b=PSallocate(lv(PSvec>>PSVEC.b), 2*n)
	if (a eq 0) % (b eq 0) resultis 0

	// a(0)=w(0)
	FST(FLDI(ac1, (w ? (w!1)+4,  4)), a)

	//a(i)=w(i)-1/a(i-1)  {i=1,2,...,n-3}
	//w(i) defaults to 4.   {1=0,1,...,n-3}
	for i=1 to n-3 do [
		FST( FSB (FLDI(ac4,(w ? (w!(i+1))+4, 4)), FDV(FLDI(ac2,1), ac1)), a+i*2)
		FLD(ac1,ac4)
		]

	if splineType eq periodicSpline then [
		// c(0)=1
		FST(one, c)
		// c(i)=-c(i-1)/a(i-1)   {i=1,2,...,n-3}
		for i=1 to n-3 do
			FST(FNEG(FDV(FLD(ac1, c+2*(i-1)), a+2*(i-1))), c+2*i)
		]

	//do everything twice to get x(t)  and y(t).
	for t=1 to 2 do [
		test ( t eq 1 )
		   ifso [ p=x; p1=p1x; p2=p2x; p3=p3x ]
		   ifnot [ p=y; p1=p1y; p2=p2y; p3=p3y ]

		computebc:
		if n ge 3 then test splineType eq naturalSpline
		ifso [
			//b(0)=6*(p(2)-2*p(1)+p(0))
			FST(FML(FAD(FSB(FSB(FLD(ac1, p+4), p+2), p+2), p), six), b)

			//b(i)=6*(p(i+2)-2*p(i+1)+p(i))-b(i-1)/a(i-1)  {i=1,2,...,n-3}
			for i=1 to n-3 do [
				FML(FLD(ac2, p+2*(i+1)), two)
				FML(FAD(FSB(FLD(ac1, p+2*(i+2)), ac2), p+2*i), six)
				FST(FSB(ac1, FDV(FLD(ac2, b+2*(i-1)), a+(i-1)*2)), b+i*2)
				]
			]
		ifnot [
			// b(0)=6*(p(1)-2*p(0)+p(n-2))
			FST(FML(FAD(FSB(FSB(FLD(ac1, p+2), p), p), p+2*(n-2)), six), b)

			// b(i)=6*(p(i+1)-2*p(i)+p(i-1))-b(i-1)/a(i-1)  {i=1,2,...,n-3}
			for i=1 to n-3 do [
				FML(FLD(ac2, p+2*i), two)
				FML(FAD(FSB(FLD(ac1, p+2*(i+1)), ac2), p+2*(i-1)), six)
				FST(FSB(ac1, FDV(FLD(ac2, b+2*(i-1)), a+(i-1)*2)), b+i*2)
				]

			// r(n-2)=1  and  s(n-2)=0
			FST(one, r+2*(n-2))
			FST(zero, s+2*(n-2), 0)
			// r(i)=-(r(i+1)+c(i))/a(i)   {i=n-3,...,1,0}
			// s(i)=(b(i)-s(i+1))/a(i)   {i=n-3,...,1,0}
			for i=n-3 to 0 by -1 do [
				FST(FDV(FNEG(FAD(FLD(ac1, r+2*(i+1)), c+2*i)), a+2*i), r+2*i)
				FST(FDV(FSB(FLD(ac1, b+2*i), s+2*(i+1)), a+2*i), s+2*i)
				]
			]

		computep2:
		// COMPUTE SECOND DERIVATIVES
		test splineType eq naturalSpline
		ifso [
			// p2(0)=p2(n-1)=0
			FST(zero, p2); FST(zero, p2+2*(n-1))
			// p2(i)=(b(i-1)-p2(i+1))/a(i-1)  {i= n-2,...,2,1}
			for i=n-2 to 1 by -1 do
				FST(FDV(FSB(FLD(ac1, b+2*(i-1)), p2+(i+1)*2), a+(i-1)*2), p2+2*i)
			]
		ifnot [
			// D2=p(n-1)-2*p(n-2)+p(n-3)
			// ac1=p2(n-2)=(6*D2-s(0)-s(n-3))/(r(0)+r(n-3)+4)
			FAD(FAD(FLD(ac2, r), r+2*(n-3)), FLDI(ac1, 4))
			FAD(FSB(FLD(ac1, p+2*(n-1)), FML(FLD(ac3, p+2*(n-2)), two)), p+2*(n-3))
			FST(FDV(FSB(FSB(FML(ac1, six), s), s+2*(n-3)), ac2), p2+2*(n-2))

			// p2(i)=r(i)*p2(n-2) + s(i)   {i=0,1,2,...,n-3}
			for i=0 to n-3 do
				FST(FAD(FML(FLD(ac2, ac1), r+2*i), s+2*i), p2+2*i)

			// p2(n-1)=p2(0)
			FST(FLD(ac1, p2), p2+2*(n-1))
			]

		computep1p3:
		// COMPUTE FIRST & THIRD DERIVATIVES
		// p1(i)=p(i+1)-p(i)-(2*p2(i)+p2(i+1))/6
		// p3(i)=p2(i+1)-p2(i)  {i=0,1,2,...,n-2}
		for i=0 to n-2 do [
			FSB(FLD(ac1, p+2*(i+1)), p+2*i)
			FAD(FML(FLD(ac2, p2+2*i), two), p2+(i+1)*2)
			FST(FSB(ac1, FDV(ac2, six)), p1+i*2)
			FST(FSB(FLD(ac1, p2+(i+1)*2), p2+i*2), p3+i*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
	]