/*
 * Header Begins--Do not remove################################################
 *
 *	C-code to exercise the plotting package by producing a teapot.
 *
 *	(c) Eugene Fiume, University of Toronto, 1995.
 *
 *	You are permitted to use this code for noncommercial, educational
 *	purposes only.  This header must not be removed from this file.
 *
 * Header Ends--Do not remove##################################################
 */

#include <stdio.h>
#include "Plot.h"
#define PI	3.141592654
#define PI4	(PI/4.0)
#define PI3	(PI/3.0)

typedef double matrix[4][4];
typedef double vector[4];

matrix xmatrix, ymatrix, zmatrix;


void GetMatrix(mat)
	matrix mat;
{
	int i,j;

	for (i=0; i<4; i++) {
	    for (j=0; j<4; j++) {
		scanf("%lf", &(mat[i][j]));
	    }
	}

} /* GetMatrix */

void PutMatrix(mat)
	matrix mat;
{
	int i,j;

	for (i=0; i<4; i++)
		printf("%f %f %f %f\n",mat[i][0],mat[i][1],mat[i][2],mat[i][3]);
	putchar('\n');

} /* PutMatrix */

/*
 *	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].
 */
double poly(s, t, mat)
	double s, t;
	matrix mat;
{
	vector T, S, sums;
	double sum;
	int i,j;

	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; i<4; i++) {			/* first evaluate MT+ */
		sums[i] = 0.0;
	    	for (j=0; j<4; j++)
			sums[i] += T[j]*mat[j][i];
	}

	sum = 0.0;
	for (i=0; i<4; i++)			/* now compute S*sums */
		sum += S[i]*sums[i];

	return(sum);
} /* poly */


double polyx(s,t)
	double s,t;
{
	return(poly(s,t,xmatrix));
} /* polyx */

double polyy(s,t)
	double s,t;
{
	return(poly(s,t,ymatrix));
} /* polyy */

double polyz(s,t)
	double s,t;
{
	return(poly(s,t,zmatrix));
} /* polyz */



main(argc, argv)
	int argc;
	char *argv[];
{
	int i, n;

	printf("Plotting 3D Parametric Polynomials\n");
	XOpen(800,800, argc, argv);

	printf("Hit a mouse button to start\n");
	GetEvent();

	SetScale(100.0,-100.0);
	SetOffset(500,500);
	SetColour(5);

	scanf("%d", &n);
	printf("Plotting %d spline segments for each of x,y,z\n", n);
	for (i=0; i<n; i++) {
		GetMatrix(xmatrix);
		GetMatrix(ymatrix);
		GetMatrix(zmatrix);
		Plot3D(polyx,polyy,polyz, 0.0,1.0, 8, 0.0,1.0, 8);
	}

/*
 *	Let's do it again for good measure.
 */
	rewind(stdin);
	SetScale(50.0,-50.0);
	SetOffset(-200,-200);
	SetColour(4);
	SetRotation(PI4, PI3, 0.0);

	scanf("%d", &n);
	for (i=0; i<n; i++) {
		GetMatrix(xmatrix);
		GetMatrix(ymatrix);
		GetMatrix(zmatrix);
		Plot3D(polyx,polyy,polyz, 0.0,1.0, 6, 0.0,1.0, 6);
	}
	printf("Hit a mouse button to end\n");
	GetEvent();
} /* main */