import var Plot in "Plot.t"
setscreen("nocursor")
setscreen("graphics:800;800")
setscreen("title: Bicubic Spline Plot")
setscreen("position: -300;-200");

var file:   int;					% file number
var filename: string := "teapot.data"	% our data file name containing
							% matrix information
type matrix:
	array 0..3, 0..3 of real

var xmatrix, ymatrix, zmatrix: matrix

procedure GetMatrix(var mat: matrix)
	for i: 0..3
		for j: 0..3
			get :file, mat(i,j)
		end for
	end for
end GetMatrix

%	Compute [s^3  s^2  s  1]M[t^3 t^2 t 1]+, where M is
%	BPB+ for basis B and observations P.
%	s,t are specific values in [0,1].

function poly(s,t: real, mat: matrix): real
	var T: array 0..3 of real
	var S: array 0..3 of real
	var sums: array 0..3 of real
	var sum: real

	T(3) := 1.0
	T(2) := t
	T(1) := t*t
	T(0) := T(1)*t
	S(3) := 1.0
	S(2) := s
	S(1) := s*s
	S(0) := S(1)*s
	for i: 0..3				% first evaluate MT+
		sums(i) := 0.0
		for j: 0..3
			sums(i) += T(j)*mat(j,i)
		end for
	end for
	sum := 0
	for i: 0..3					% now compute S*sums
		sum += S(i)*sums(i)
	end for
	result(sum)
end poly

function polyx(s,t: real): real
	result(poly(s,t,xmatrix))
end polyx

function polyy(s,t: real): real
	result(poly(s,t,ymatrix))
end polyy

function polyz(s,t: real): real
	result(poly(s,t,zmatrix))
end polyz

var number: int;

Plot.SetScale(100,100)
Plot.SetOffset(500,500)
Plot.SetColour(2)

open : file, filename, get
get : file, number
put "Plotting ", number," spline segments for each of x,y,z"

for i: 1..number
	GetMatrix(xmatrix)
	GetMatrix(ymatrix)
	GetMatrix(zmatrix)
	Plot.Plot3D(polyx,polyy,polyz, 0,1, 8, 0,1, 8)
end for