// File: SPLINE2.BCPL
// P.Baudelaire
// December 5, 1977  4:45 PM

// The algorithms implemented by these procedures are 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
	CubicSpline
	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; FSN
	FPSetup

	Allocate		// Alto SYSTEM
	Free
	Zero
	]

// incoming statics:

external [
	FPwork			// microFLOAT  (floating point registers)
	]

// 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
	h 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, d1x, d2x, d3x, d1y, d2y, d3y, splineType, nil;
numargs nargs) = valof [

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

	if n ls 0 then [
		n=-n
		if ConvertToFP(n, lv x, lv y) eq 0 resultis 0
		]

	switchon nargs into [
		case 9:
			splineType=naturalSpline
		case 10:
		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))
		]

	// compute parametrization (polygonal line approximation)
	let h=PSallocate(lv(PSvec>>PSVEC.h), 2*n)
	if h eq 0 resultis 0
	for i=0 to n-2 do [
		FSB(FLD(ac2, x+2*(i+1)), x+2*i)
		if FSN(ac2) eq -1 then FNEG(ac2)
		FSB(FLD(ac3, y+2*(i+1)), y+2*i)
		if FSN(ac3) eq -1 then FNEG(ac3)
		FDV((FCM(ac2, ac3) eq 1 ? ac3, ac2), two)
		FST(FAD(ac2, ac3), h+2*i)
		]

	// now, compute the cubic splines x(h) & y(h)
	let done=ComputeCubicSpline(n, h, x, d1x, d2x, d3x, splineType)
	if done then
		done=ComputeCubicSpline(n, h, y, d1y, d2y, d3y, splineType)

	// reparametrize to the interval [0, 1]
	if done then for i=0 to n-2 do [
		FLD(ac3, FLD(ac1, h+2*i))
		FST(FML(FLD(ac2, d1x+2*i), ac1), d1x+2*i)
		FST(FML(FLD(ac2, d1y+2*i), ac1), d1y+2*i)
		FML(ac1, ac3)
		FST(FML(FLD(ac2, d2x+2*i), ac1), d2x+2*i)
		FST(FML(FLD(ac2, d2y+2*i), ac1), d2y+2*i)
		FML(ac1, ac3)
		FST(FML(FLD(ac2, d3x+2*i), ac1), d3x+2*i)
		FST(FML(FLD(ac2, d3y+2*i), ac1), d3y+2*i)
		]

	resultis PSquit(done)
	]



and CubicSpline(n,x,y,d1y,d2y,d3y,splineType; numargs nargs) = valof [

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

	if n ls 0 then [
		n=-n
		if ConvertToFP(n, lv x, lv y) eq 0 resultis 0
		]

	switchon nargs into [
		case 6:
			splineType=naturalSpline
		case 7:
			if splineType ne naturalSpline &
			   splineType ne periodicSpline resultis PSquit(PSerror(3))
			if n ls 3 then splineType=naturalSpline
			endcase
		default:
			resultis PSquit(PSerror(4))
		]

	let h=PSallocate(lv(PSvec>>PSVEC.h), 2*n)
	if h eq 0 resultis 0
	for i=0 to n-2 do FST(FSB(FLD(ac1, x+2*(i+1)), x+2*i), h+2*i)

	// now, compute
	let done=ComputeCubicSpline(n, h, y, d1y, d2y, d3y, splineType)
	if done eq 0 resultis 0
	resultis PSquit()
	]



and ComputeCubicSpline(n, h, y, d1, d2, d3, splineType) = valof [
	// IMPORTANT: assume all arguments are right !!!!
	// ONLY called by CubicSpline & ParametricSpline
	let c,r,s=0,0,0
	if splineType eq periodicSpline then [
		if 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

	// check overlapping knots!
	for i=0 to n-1 do
		if FCM(zero, h+2*i) eq 0 then FST(one, h+2*i)

	test splineType eq naturalSpline
	ifso [
		// a(0)=2*(h(0)+h(1))
		FST(FML(FAD(FLD(ac1, h), h+2), two), a)
		// a(i)=2*(h(i)+h(i+1))-h(i)*h(i)/a(i-1)  {i=1,2,...,n-3}
		for i=1 to n-3 do [
			FML(FAD(FLD(ac4, h+2*i), h+2*(i+1)), two)
			FDV(FML(FLD(ac2,h+2*i), ac2), ac1)
			FLD(ac1, FST(FSB(ac4, ac2), a+2*i))
			]
		]

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

		// c(0)=h(n-2)
		FST(FLD(ac1, h+2*(n-2)), c)
		// c(i)=-h(i-1)*c(i-1)/a(i-1)   {i=1,2,...,n-3}
		// Notice: c(i-1) is ac1
		for i=1 to n-3 do
			FST(FNEG(FDV(FML(ac1, h+2*(i-1)), a+2*(i-1))), c+2*i)
		]

	computebc:
	if n ge 3 then test splineType eq naturalSpline
	ifso [
		// first: D(i)=(y(i+1)-y(i))/h(i)   {i=0,1,...,n-2}
		// and: b(i)=6*(D(i+1)-D(i))   {i=0,1,...,n-3}
		// Notice: D(i) is ac1, D(i+1) is ac2
		FDV(FSB(FLD(ac1, y+2), y), h)
		for i=0 to n-3 do [
			FDV(FSB(FLD(ac2, y+2*(i+2)), y+2*(i+1)), h+2*(i+1))
			FST(FML(FSB(FLD(ac3, ac2), ac1), six), b+2*i)
			FLD(ac1, ac2)
			]
		// then: b(i)=b(i)-h(i)*b(i-1)/a(i-1)  {i=1,2,...,n-3}
		// Notice: b(i-1) is ac1
		FLD(ac1, b)
		for i=1 to n-3 do [
			FDV(FML(FLD(ac2, ac1), h+2*i), a+2*(i-1))
			FST(FSB(FLD(ac1, b+2*i), ac2), b+i*2)
			]
		]

	ifnot [
		// first: D(i)=(y(i+1)-y(i))/h(i)   {i=0,1,...,n-2}
		// and: b(i)=6*(D(i)-D(i-1))   {i=1,2,...,n-3}
		// and: b(0)= 6*(D(0)-D(n-2))
		// Notice: D(i-1) is ac1, D(i) is ac2
		FDV(FSB(FLD(ac1, y), y+2*(n-2)), h+2*(n-2))
		for i=0 to n-3 do [
			FDV(FSB(FLD(ac2, y+2*(i+1)), y+2*i), h+2*i)
			FST(FML(FSB(FLD(ac3, ac2), ac1), six), b+2*i)
			FLD(ac1, ac2)
			]
		// then: b(i)=b(i)-h(i-1)*b(i-1)/a(i-1)  {i=1,2,...,n-3}
		// Notice: b(i-1) is ac1
		FLD(ac1, b)
		for i=1 to n-3 do [
			FDV(FML(FLD(ac2, ac1), h+2*(i-1)), a+2*(i-1))
			FST(FSB(FLD(ac1, b+2*i), ac2), b+i*2)
			]

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

	computed2:
	// COMPUTE SECOND DERIVATIVES
	test splineType eq naturalSpline
	ifso [
		// d2(0)=d2(n-1)=0
		FST(zero,d2); FST(zero,d2+2*(n-1))
		// d2(i)=(b(i-1)-h(i)*d2(i+1))/a(i-1)  {i=n-2,...,2,1}
		for i=n-2 to 1 by -1 do [
			FML(FLD(ac2, d2+2*(i+1)), h+2*i)
			FST(FDV(FSB(FLD(ac1, b+2*(i-1)), ac2), a+2*(i-1)), d2+2*i)
			]
		]

	ifnot [
		// ac1=(6*(D(n-2)-D(n-3))-h(n-2)*s(0)-h(n-3)*s(n-3))
		FDV(FSB(FLD(ac1, y+2*(n-1)), y+2*(n-2)), h+2*(n-2))
		FSB(ac1, FDV(FSB(FLD(ac2, y+2*(n-2)), y+2*(n-3)), h+2*(n-3)))
		FSB(FML(ac1, six), FML(FLD(ac3, s), h+2*(n-2)))
		FSB(ac1, FML(FLD(ac3, s+2*(n-3)), h+2*(n-3)))

		// ac2=(h(n-2)*r(0)+h(n-3)*r(n-3)+2*(h(n-3)+h(n-2)))
		FAD(FML(FLD(ac2, r), h+2*(n-2)), FML(FLD(ac3, r+2*(n-3)), h+2*(n-3)))
		FAD(ac2, FML(FAD(FLD(ac3, h+2*(n-3)), h+2*(n-2)), two))

		// ac1=d2(n-2)=ac1/ac2
		FST(FDV(ac1, ac2), d2+2*(n-2))

		// d2(i)=r(i)*d2(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), d2+2*i)

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

	computed1d3:
	// COMPUTE FIRST & THIRD DERIVATIVES
	// d1(i)=(y(i+1)-y(i))/h(i) - h(i)*(2*d2(i)+d2(i+1))/6
	// d3(i)=(d2(i+1)-d2(i))/h(i)   {i=0,1,2,...,n-2}
	for i=0 to n-2 do [
		FDV(FSB(FLD(ac1, y+2*(i+1)), y+2*i), h+2*i)
		FAD(FML(FLD(ac2, d2+2*i), two), d2+2*(i+1))
		FST(FSB(ac1, FDV(FML(ac2, h+2*i), six)), d1+2*i)
		FST(FDV(FSB(FLD(ac1, d2+2*(i+1)), d2+2*i), h+2*i), d3+2*i)
		]

	resultis true
	]


and ConvertToFP(n, intXloc, intYloc) = valof [
	// convert coordinates from integer to floating point
	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 PSquit(PSerror(1))
	for i=0 to n-1 do [
		FST(FLDI(ac1, (@intXloc)!i), fpx+2*i)
		FST(FLDI(ac1, (@intYloc)!i), fpy+2*i)
		]
	@intXloc=fpx
	@intYloc=fpy
	resultis true
	]


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


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 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
	]