// April 14, 1978  5:49 PM				*** overlay C ***


get "zpDefs.bcpl"

get "zpPressDf.bcpl"

// outgoing procedures:

external [
	WriteDLspline
	]

// incoming procedures:

external [
	MoveBlock			// SYSTEM
	SetBlock
	Gets

	ParametricSpline	// SPLINE

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

	SquareRoot			// ZPUTIL
	typeForm

	obtainBlock		// ZPBLOCK
	putBlock
	giveUp

	ScaleKnotTable		// ZPPIECE
	GetThickness
	ComputeNormal
	MakeCubic

	PutObjectInELtable	// ZPPRESS
	PrintFpValues

	WritePieceObject			// ZPOBJECT2
	]

// incoming statics:

external [
	keys			// SYSTEM

	@computeRange		// ZPMAKE
	@denRate
	@numRate
	@pseudoCyclic
	@splineTension

	@showValues		// ZPPRESS
	]


// local statics

static [
	@prec=20
	]


// local floating point registers

manifest [
	t0=0; t1=1; t2=2; t3=3
	]


let WriteDLspline(pressFile, splinePointer, sLeft, sBottom) = valof [WriteDLspline
	let thickness=GetThickness(splinePointer)
	let objectWordCount=0
	let procedureName="[WriteDLspline]"
	let brushShape=splinePointer>>SPLINE.shape
	let n=splinePointer>>SPLINE.nKnots
	let dn=2*n
	// 4 extra knots for cyclic curve
	let cyclic=splinePointer>>SPLINE.cyclic
	let xn=(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 for cyclic curve
		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 [
		//caution: >4 are global floating point registers
		X=5; Y=6
		d1X=7; d2X=8; d3X=9
		d1Y=10; d2Y=11; d3Y=12
		d1TX=13; d2TX=14
		d1TY=15; d2TY=16
		deltaT=17
		]

	// constants
	let two=vec 2; FLDI(t0,2); FST(t0,two)
	let three=vec 2; FLDI(t0,3); FST(t0,three)
	let six=vec 2; FLDI(t0,6); FST(t0,six)
	let epsilon=vec 2
	FLDI(t0,prec); FLDI(t1,100); FDV(t0,t1); FST(t0,epsilon)

	// scale up
	ScaleKnotTable(kTable, kTable+dnt, nt, scaleFactor, sLeft, sBottom)

	let wTable=0
	if splineTension then [
		wTable=obtainBlock(nt)
		if wTable ne 0 then [
			SetBlock(wTable, 0)
			typeForm(0,"Input spline weights now @", 8, wTable, 1, $*N)
			Gets(keys)
			]
		]

	// spline computation:
	// caution: ParametricSpline(n,X,Y,X',X'',X''',Y',Y'',Y''',type,W)
	let psDone=ParametricSpline(nt, kTable, kTable+dnt, dTable, dTable+dnt, dTable+2*dnt, 
		dTable+3*dnt, dTable+4*dnt, dTable+5*dnt,
		((cyclic & not pseudoCyclic) ? periodicSpline, naturalSpline),
		wTable)

	putBlock(wTable)

	unless psDone resultis giveUp(procedureName, dTable, kTable)

	for k=0 to n-2 do [
	  // knot k - floating point
	  let kX,kY=xTable+2*k,yTable+2*k
	  // stepping parameters
	  FLDI(t2, scaleFactor)
	  FLD(t0,kX); FSB(t0,kX+2); if FSN(t0) eq -1 then FNEG(t0); FDV(t0, t2)
	  FLD(t1,kY); FSB(t1,kY+2); if FSN(t1) eq -1 then FNEG(t1); FDV(t1, t2)
	  let m=(((FCM(t0,t1) eq 1) ? FTR(t0), FTR(t1))*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; FLDI(t0,ni); FST(t0,fni)
	  FLDI(t1,1)
	  if m ne 0 then [ FLDI(t0,m); FDV(t1,t0) ]
	  let delta=vec 2;  FST(t1,delta)
	  let delta2=vec 2; FML(t1,delta); FST(t1,delta2)
	  let delta3=vec 2; FML(t1,delta); FST(t1,delta3)
	  // derivatives @ knot - 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
	  // start of first piece
	  let xStart=vec 2; FLD(t0,kX); FST(t0,xStart)
	  let yStart=vec 2; FLD(t0,kY); FST(t0,yStart)
	  let txStart=vec 2; FLD(t0,kX1); FST(t0,txStart)
	  let tyStart=vec 2; FLD(t1,kY1); FST(t1,tyStart)
	  let tzStart=vec 2; FML(t0,txStart); FML(t1,tyStart); FAD(t0,t1); FST(t0,tzStart)
	  FLDI(deltaT,0)
	  let xEnd=vec 2; let yEnd=vec 2; let txEnd=vec 2; let tyEnd=vec 2

	  // derivatives @ start 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
	  // forward differences - double precision fixed point -
	  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
	  let dpTX=vec 2; let dpd1TX=vec 2; let dpd2TX=vec 2
	  let dpTY=vec 2; let dpd1TY=vec 2; let dpd2TY=vec 2
	  // 3rd derivatives & differences are constant
	  FLD(d3X,delta3); FML(d3X,X3); FSTDP(d3X,dpd3X)
	  FLD(d3Y,delta3); FML(d3Y,Y3); FSTDP(d3Y,dpd3Y)
	  FLD(d2TX,delta2); FML(d2TX,X3); FSTDP(d2TX,dpd2TX)
	  FLD(d2TY,delta2); FML(d2TY,Y3); FSTDP(d2TY,dpd2TY)
	  // forward difference computation in subinterval
	  for i=1 to ni do [
	    // initial values of forward differences (X & Y)
	    //		(floating point)
	    FLD(d2X,delta2); FML(d2X,X2)
	    FLD(d2Y,delta2); FML(d2Y,Y2)
	    FLD(d1X,delta);  FML(d1X,X1)
	    FLD(t0,d2X);   FDV(t0,two);    FAD(d1X,t0)
	    FLD(t0,d3X);   FDV(t0,six);    FAD(d1X,t0)
	    FLD(d1Y,delta);  FML(d1Y,Y1)
	    FLD(t0,d2Y);   FDV(t0,two);    FAD(d1Y,t0)
	    FLD(t0,d3Y);   FDV(t0,six);    FAD(d1Y,t0)
	    FAD(d2X,d3X);    FAD(d2Y,d3Y)
	    //		(double precision)
	    FSTDP(X,dpX); FSTDP(d1X,dpd1X); FSTDP(d2X,dpd2X)
	    FSTDP(Y,dpY); FSTDP(d1Y,dpd1Y); FSTDP(d2Y,dpd2Y)
	    // initial values of forward differences (TX & TY)
	    //		(floating point)
	    FLD(d1TX,delta); FML(d1TX,X2);
	    FLD(t0,d2TX); FDV(t0,two); FAD(d1TX,t0)
	    FLD(d1TY,delta); FML(d1TY,Y2);
	    FLD(t0,d2TY); FDV(t0,two); FAD(d1TY,t0)
	    //		(double precision)
	    FLD(t0,X1); FSTDP(t0,dpTX); FSTDP(d1TX,dpd1TX)
	    FLD(t0,Y1); FSTDP(t0,dpTY); FSTDP(d1TY,dpd1TY)
	    //compute
	    for j=1 to r do [
		DPAD(dpX,dpd1X); DPAD(dpY,dpd1Y)
		DPAD(dpd1X,dpd2X); DPAD(dpd1Y,dpd2Y)
		DPAD(dpd2X,dpd3X); DPAD(dpd2Y,dpd3Y)
		// check variation in first derivative
		DPAD(dpTX,dpd1TX); DPAD(dpTY,dpd1TY)
		DPAD(dpd1TX,dpd2TX); DPAD(dpd1TY,dpd2TY)
		FAD(deltaT,delta)

		FLDDP(t0,dpTX); FML(t0,tyStart)
		FLDDP(t1,dpTY); FML(t1,txStart)
		FSB(t0,t1); FDV(t0,tzStart)	

		if FSN(t0) eq -1 then FNEG(t0)
		if FCM(t0,epsilon) eq 1 then [
//			if showValues then [
//				let temp=vec 2
//				FST(t0, temp)
//				PrintFpValues("WDLspline: " , temp)
//				]
			FLDDP(t0,dpX); FST(t0,xEnd)
			FLDDP(t0,dpY); FST(t0,yEnd)
			FLDDP(t0,dpTX); FML(t0,deltaT); FST(t0,txEnd)
			FLDDP(t0,dpTY); FML(t0,deltaT); FST(t0,tyEnd)
			FLD(t0,txStart); FML(t0,deltaT); FST(t0,txStart)
			FLD(t0,tyStart); FML(t0,deltaT); FST(t0,tyStart)
			let PECstart= fPEC
			if not cyclic & (objectWordCount eq 0) then
				PECstart=selecton brushShape into [
					case sBrush: sPEC;
					case rBrush: rPEC;
					default: fPEC ]
			objectWordCount=objectWordCount +
				WritePieceObject(pressFile, thickness, xStart, yStart,
					txStart, tyStart, xEnd, yEnd, txEnd, tyEnd,
					PECstart, fPEC)
			FLDDP(t0,dpX); FST(t0,xStart)
			FLDDP(t1,dpY); FST(t1,yStart)
			FLDDP(t0,dpTX); FST(t0,txStart); FML(t0,txStart)
			FLDDP(t1,dpTY); FST(t1,tyStart); FML(t1,tyStart)
			FAD(t0,t1); FST(t0,tzStart)
			FLDI(deltaT,0)
			]
		]
	    if i ne ni then [
		//starting point of next interval
		FLDI(t1,i); FDV(t1,fni)
		FLD(t2,t1); FDV(t2,two)
		FLD(t3,t1); FDV(t3,three)
		FLD(t0,kX3); FML(t0,t1); FAD(t0,kX2)
		FST(t0,X2)
		FLD(t0,kY3); FML(t0,t1); FAD(t0,kY2)
		FST(t0,Y2)
		FLD(t0,kX3); FML(t0,t2); FAD(t0,kX2); FML(t0,t1); FAD(t0,kX1)
		FST(t0,X1)
		FLD(t0,kY3); FML(t0,t2); FAD(t0,kY2); FML(t0,t1); FAD(t0,kY1)
		FST(t0,Y1)
		FLD(X,kX3); FML(X,t3); FAD(X,kX2)
		FML(X,t2); FAD(X,kX1); FML(X,t1); FAD(X,kX)
		FLD(Y,kY3); FML(Y,t3); FAD(Y,kY2)
		FML(Y,t2); FAD(Y,kY1); FML(Y,t1); FAD(Y,kY)
 		] 
	    ]

	// last piece for this section
	FLD(t0,kX+2); FST(t0,xEnd)
	FLD(t0,kY+2); FST(t0,yEnd)
	FLD(t0,kX3); FDV(t0,two); FAD(t0,kX2); FAD(t0,kX1)
	FML(t0,deltaT); FST(t0,txEnd)
	FLD(t0,kY3); FDV(t0,two); FAD(t0,kY2); FAD(t0,kY1)
	FML(t0,deltaT); FST(t0,tyEnd)
	FLD(t0,txStart); FML(t0,deltaT); FST(t0,txStart)
	FLD(t0,tyStart); FML(t0,deltaT); FST(t0,tyStart)
	let PECstart= fPEC
	if not cyclic & (objectWordCount eq 0) then
		PECstart=selecton brushShape into [
			case sBrush: sPEC;
			case rBrush: rPEC;
			default: fPEC ]
	let PECend=fPEC
	if not cyclic & (k eq (n-2)) then PECend=selecton brushShape into [
				case sBrush: sPEC;
				case rBrush: rPEC;
				default: fPEC ]
	objectWordCount=objectWordCount + WritePieceObject(pressFile, thickness,
		xStart, yStart, txStart, tyStart, xEnd, yEnd, txEnd, tyEnd, PECstart, PECend)
	typeForm(1,$,) 
	]

	//return storage
	putBlock(dTable)
	putBlock(kTable)
	resultis objectWordCount
	]WriteDLspline