src/gen/gensurf.c

/* Copyright (c) 1989 Regents of the University of California */

#ifndef lint
static char SCCSid[] = "$SunId$ LBL";
#endif

/*
 *  gensurf.c - program to generate functional surfaces
 *
 *      Parametric functions x(s,t), y(s,t) and z(s,t)
 *  specify the surface, which is tesselated into an m by n
 *  array of paired triangles.
 *      The surface normal is defined by the right hand
 *  rule applied to (s,t).
 *
 *      4/3/87
 */

#include  "standard.h"

#define  XNAME          "X_"                    /* x function name */
#define  YNAME          "Y_"                    /* y function name */
#define  ZNAME          "Z_"                    /* z function name */

#define  ABS(x)         ((x)>=0 ? (x) : -(x))

#define  pvect(p)       printf(vformat, (p)[0], (p)[1], (p)[2])

char  vformat[] = "%15.9g %15.9g %15.9g\n";
char  tsargs[] = "4 surf_dx surf_dy surf_dz surf.cal\n";
char  texname[] = "Phong";

int  smooth = 0;                /* apply smoothing? */

char  *modname, *surfname;

double  funvalue(), l_hermite(), argument();

typedef struct {
        FVECT  p;       /* vertex position */
        FVECT  n;       /* average normal */
} POINT;


main(argc, argv)
int  argc;
char  *argv[];
{
        POINT  *row0, *row1, *row2, *rp;
        int  i, j, m, n;
        char  stmp[256];

        varset("PI", PI);
        funset("hermite", 5, l_hermite);

        if (argc < 8)
                goto userror;

        for (i = 8; i < argc; i++)
                if (!strcmp(argv[i], "-e"))
                        scompile(NULL, argv[++i]);
                else if (!strcmp(argv[i], "-f"))
                        fcompile(argv[++i]);
                else if (!strcmp(argv[i], "-s"))
                        smooth++;
                else
                        goto userror;

        modname = argv[1];
        surfname = argv[2];
        sprintf(stmp, "%s(s,t)=%s;", XNAME, argv[3]);
        scompile(NULL, stmp);
        sprintf(stmp, "%s(s,t)=%s;", YNAME, argv[4]);
        scompile(NULL, stmp);
        sprintf(stmp, "%s(s,t)=%s;", ZNAME, argv[5]);
        scompile(NULL, stmp);
        m = atoi(argv[6]);
        n = atoi(argv[7]);
        if (m <= 0 || n <= 0)
                goto userror;

        row0 = (POINT *)malloc((n+3)*sizeof(POINT));
        row1 = (POINT *)malloc((n+3)*sizeof(POINT));
        row2 = (POINT *)malloc((n+3)*sizeof(POINT));
        if (row0 == NULL || row1 == NULL || row2 == NULL) {
                fprintf(stderr, "%s: out of memory\n", argv[0]);
                quit(1);
        }
        row0++; row1++; row2++;
                                                /* print header */
        printhead(argc, argv);
                                                /* initialize */
        comprow(-1.0/m, row0, n);
        comprow(0.0, row1, n);
        comprow(1.0/m, row2, n);
        compnorms(row0, row1, row2, n);
                                                /* for each row */
        for (i = 0; i < m; i++) {
                                                /* compute next row */
                rp = row0;
                row0 = row1;
                row1 = row2;
                row2 = rp;
                comprow((double)(i+2)/m, row2, n);
                compnorms(row0, row1, row2, n);

                for (j = 0; j < n; j++) {
                                                        /* put polygons */
                        if ((i+j) & 1)
                                putsquare(&row0[j], &row1[j],
                                                &row0[j+1], &row1[j+1]);
                        else
                                putsquare(&row1[j], &row1[j+1],
                                                &row0[j], &row0[j+1]);
                }
        }

        quit(0);

userror:
        fprintf(stderr, "Usage: %s material name ", argv[0]);
        fprintf(stderr, "x(s,t) y(s,t) z(s,t) m n [-s][-e expr][-f file]\n");
        quit(1);
}


putsquare(p0, p1, p2, p3)               /* put out a square */
POINT  *p0, *p1, *p2, *p3;
{
        static int  nout = 0;
        FVECT  norm[4];
        int  axis;
        FVECT  v1, v2, vc1, vc2;
        int  ok1, ok2;
                                        /* compute exact normals */
        fvsum(v1, p1->p, p0->p, -1.0);
        fvsum(v2, p2->p, p0->p, -1.0);
        fcross(vc1, v1, v2);
        ok1 = normalize(vc1) != 0.0;
        fvsum(v1, p2->p, p3->p, -1.0);
        fvsum(v2, p1->p, p3->p, -1.0);
        fcross(vc2, v1, v2);
        ok2 = normalize(vc2) != 0.0;
        if (!(ok1 | ok2))
                return;
                                        /* compute normal interpolation */
        axis = norminterp(norm, p0, p1, p2, p3);

                                        /* put out quadrilateral? */
        if (ok1 & ok2 && fdot(vc1,vc2) >= 1.0-FTINY*FTINY) {
                printf("\n%s ", modname);
                if (axis != -1) {
                        printf("texfunc %s\n", texname);
                        printf(tsargs);
                        printf("0\n13\t%d\n", axis);
                        pvect(norm[0]);
                        pvect(norm[1]);
                        pvect(norm[2]);
                        fvsum(v1, norm[3], vc1, -0.5);
                        fvsum(v1, v1, vc2, -0.5);
                        pvect(v1);
                        printf("\n%s ", texname);
                }
                printf("polygon %s.%d\n", surfname, ++nout);
                printf("0\n0\n12\n");
                pvect(p0->p);
                pvect(p1->p);
                pvect(p3->p);
                pvect(p2->p);
                return;
        }
                                        /* put out triangles? */
        if (ok1) {
                printf("\n%s ", modname);
                if (axis != -1) {
                        printf("texfunc %s\n", texname);
                        printf(tsargs);
                        printf("0\n13\t%d\n", axis);
                        pvect(norm[0]);
                        pvect(norm[1]);
                        pvect(norm[2]);
                        fvsum(v1, norm[3], vc1, -1.0);
                        pvect(v1);
                        printf("\n%s ", texname);
                }
                printf("polygon %s.%d\n", surfname, ++nout);
                printf("0\n0\n9\n");
                pvect(p0->p);
                pvect(p1->p);
                pvect(p2->p);
        }
        if (ok2) {
                printf("\n%s ", modname);
                if (axis != -1) {
                        printf("texfunc %s\n", texname);
                        printf(tsargs);
                        printf("0\n13\t%d\n", axis);
                        pvect(norm[0]);
                        pvect(norm[1]);
                        pvect(norm[2]);
                        fvsum(v2, norm[3], vc2, -1.0);
                        pvect(v2);
                        printf("\n%s ", texname);
                }
                printf("polygon %s.%d\n", surfname, ++nout);
                printf("0\n0\n9\n");
                pvect(p2->p);
                pvect(p1->p);
                pvect(p3->p);
        }
}


comprow(s, row, siz)                    /* compute row of values */
double  s;
register POINT  *row;
int  siz;
{
        double  st[2];
        register int  i;
                                        /* compute one past each end */
        st[0] = s;
        for (i = -1; i <= siz+1; i++) {
                st[1] = (double)i/siz;
                row[i].p[0] = funvalue(XNAME, 2, st);
                row[i].p[1] = funvalue(YNAME, 2, st);
                row[i].p[2] = funvalue(ZNAME, 2, st);
        }
}


compnorms(r0, r1, r2, siz)              /* compute row of averaged normals */
register POINT  *r0, *r1, *r2;
int  siz;
{
        FVECT  v1, v2, vc;
        register int  i;

        if (!smooth)                    /* not needed if no smoothing */
                return;
                                        /* compute middle points */
        while (siz-- >= 0) {
                fvsum(v1, r2[0].p, r1[0].p, -1.0);
                fvsum(v2, r1[1].p, r1[0].p, -1.0);
                fcross(r1[0].n, v1, v2);
                fvsum(v1, r0[0].p, r1[0].p, -1.0);
                fcross(vc, v2, v1);
                fvsum(r1[0].n, r1[0].n, vc, 1.0);
                fvsum(v2, r1[-1].p, r1[0].p, -1.0);
                fcross(vc, v1, v2);
                fvsum(r1[0].n, r1[0].n, vc, 1.0);
                fvsum(v1, r2[0].p, r1[0].p, -1.0);
                fcross(vc, v2, v1);
                fvsum(r1[0].n, r1[0].n, vc, 1.0);
                normalize(r1[0].n);
                r0++; r1++; r2++;
        }
}


int
norminterp(resmat, p0, p1, p2, p3)      /* compute normal interpolation */
register FVECT  resmat[4];
POINT  *p0, *p1, *p2, *p3;
{
#define u  ((ax+1)%3)
#define v  ((ax+2)%3)

        register int  ax;
        double  eqnmat[4][4];
        FVECT  v1;
        register int  i, j;

        if (!smooth)                    /* no interpolation if no smoothing */
                return(-1);
                                        /* find dominant axis */
        VCOPY(v1, p0->n);
        fvsum(v1, v1, p1->n, 1.0);
        fvsum(v1, v1, p2->n, 1.0);
        fvsum(v1, v1, p3->n, 1.0);
        ax = ABS(v1[0]) > ABS(v1[1]) ? 0 : 1;
        ax = ABS(v1[ax]) > ABS(v1[2]) ? ax : 2;
                                        /* assign equation matrix */
        eqnmat[0][0] = p0->p[u]*p0->p[v];
        eqnmat[0][1] = p0->p[u];
        eqnmat[0][2] = p0->p[v];
        eqnmat[0][3] = 1.0;
        eqnmat[1][0] = p1->p[u]*p1->p[v];
        eqnmat[1][1] = p1->p[u];
        eqnmat[1][2] = p1->p[v];
        eqnmat[1][3] = 1.0;
        eqnmat[2][0] = p2->p[u]*p2->p[v];
        eqnmat[2][1] = p2->p[u];
        eqnmat[2][2] = p2->p[v];
        eqnmat[2][3] = 1.0;
        eqnmat[3][0] = p3->p[u]*p3->p[v];
        eqnmat[3][1] = p3->p[u];
        eqnmat[3][2] = p3->p[v];
        eqnmat[3][3] = 1.0;
                                        /* invert matrix (solve system) */
        if (!invmat(eqnmat, eqnmat))
                return(-1);                     /* no solution */
                                        /* compute result matrix */
        for (j = 0; j < 4; j++)
                for (i = 0; i < 3; i++)
                        resmat[j][i] =  eqnmat[j][0]*p0->n[i] +
                                        eqnmat[j][1]*p1->n[i] +
                                        eqnmat[j][2]*p2->n[i] +
                                        eqnmat[j][3]*p3->n[i];
        return(ax);

#undef u
#undef v
}


/*
 * invmat - computes the inverse of mat into inverse.  Returns 1
 * if there exists an inverse, 0 otherwise.  It uses Gaussian Elimination
 * method.
 */

invmat(inverse,mat)
double mat[4][4],inverse[4][4];
{
#define SWAP(a,b,t) (t=a,a=b,b=t)

        double  m4tmp[4][4];
        register int i,j,k;
        register double temp;

        bcopy((char *)mat, (char *)m4tmp, sizeof(m4tmp));
                                        /* set inverse to identity */
        for (i = 0; i < 4; i++)
                for (j = 0; j < 4; j++)
                        inverse[i][j] = i==j ? 1.0 : 0.0;

        for(i = 0; i < 4; i++) {
                /* Look for raw with largest pivot and swap raws */
                temp = FTINY; j = -1;
                for(k = i; k < 4; k++)
                        if(ABS(m4tmp[k][i]) > temp) {
                                temp = ABS(m4tmp[k][i]);
                                j = k;
                                }
                if(j == -1)     /* No replacing raw -> no inverse */
                        return(0);
                if (j != i)
                        for(k = 0; k < 4; k++) {
                                SWAP(m4tmp[i][k],m4tmp[j][k],temp);
                                SWAP(inverse[i][k],inverse[j][k],temp);
                                }

                temp = m4tmp[i][i];
                for(k = 0; k < 4; k++) {
                        m4tmp[i][k] /= temp;
                        inverse[i][k] /= temp;
                        }
                for(j = 0; j < 4; j++) {
                        if(j != i) {
                                temp = m4tmp[j][i];
                                for(k = 0; k < 4; k++) {
                                        m4tmp[j][k] -= m4tmp[i][k]*temp;
                                        inverse[j][k] -= inverse[i][k]*temp;
                                        }
                                }
                        }
                }
        return(1);

#undef SWAP
}


eputs(msg)
char  *msg;
{
        fputs(msg, stderr);
}


wputs(msg)
char  *msg;
{
        eputs(msg);
}


quit(code)
{
        exit(code);
}


printhead(ac, av)               /* print command header */
register int  ac;
register char  **av;
{
        putchar('#');
        while (ac--) {
                putchar(' ');
                fputs(*av++, stdout);
        }
        putchar('\n');
}


double
l_hermite()                     
{
        double  t;
        
        t = argument(5);
        return( argument(1)*((2.0*t-3.0)*t*t+1.0) +
                argument(2)*(-2.0*t+3.0)*t*t +
                argument(3)*((t-2.0)*t+1.0)*t +
                argument(4)*(t-1.0)*t*t );
}
Revision:	1.5
Committed:	Fri Jan 19 00:03:08 1990 UTC (34 years, 3 months ago) by greg
Content type:	text/plain
Branch:	MAIN
Changes since 1.4:	+2 -7 lines
Log Message:	improved portability of bcopy()
#	User	Rev	Content
1	greg	1.2	/* Copyright (c) 1989 Regents of the University of California */
2	greg	1.1
3			#ifndef lint
4			static char SCCSid[] = "$SunId$ LBL";
5			#endif
6	greg	1.2
7			/*
8	greg	1.1	* gensurf.c - program to generate functional surfaces
9			*
10			* Parametric functions x(s,t), y(s,t) and z(s,t)
11			* specify the surface, which is tesselated into an m by n
12			* array of paired triangles.
13			* The surface normal is defined by the right hand
14			* rule applied to (s,t).
15			*
16			* 4/3/87
17			*/
18
19	greg	1.5	#include "standard.h"
20	greg	1.1
21			#define XNAME "X_" /* x function name */
22			#define YNAME "Y_" /* y function name */
23			#define ZNAME "Z_" /* z function name */
24
25	greg	1.4	#define ABS(x) ((x)>=0 ? (x) : -(x))
26
27	greg	1.3	#define pvect(p) printf(vformat, (p)[0], (p)[1], (p)[2])
28	greg	1.1
29			char vformat[] = "%15.9g %15.9g %15.9g\n";
30	greg	1.3	char tsargs[] = "4 surf_dx surf_dy surf_dz surf.cal\n";
31			char texname[] = "Phong";
32	greg	1.1
33	greg	1.3	int smooth = 0; /* apply smoothing? */
34	greg	1.1
35	greg	1.3	char modname, surfname;
36	greg	1.1
37	greg	1.4	double funvalue(), l_hermite(), argument();
38	greg	1.3
39			typedef struct {
40			FVECT p; /* vertex position */
41			FVECT n; /* average normal */
42			} POINT;
43
44
45	greg	1.1	main(argc, argv)
46			int argc;
47			char *argv[];
48			{
49	greg	1.3	POINT row0, row1, row2, rp;
50	greg	1.1	int i, j, m, n;
51			char stmp[256];
52
53			varset("PI", PI);
54			funset("hermite", 5, l_hermite);
55
56			if (argc < 8)
57			goto userror;
58
59			for (i = 8; i < argc; i++)
60			if (!strcmp(argv[i], "-e"))
61			scompile(NULL, argv[++i]);
62			else if (!strcmp(argv[i], "-f"))
63			fcompile(argv[++i]);
64	greg	1.3	else if (!strcmp(argv[i], "-s"))
65			smooth++;
66	greg	1.1	else
67			goto userror;
68
69	greg	1.3	modname = argv[1];
70			surfname = argv[2];
71	greg	1.1	sprintf(stmp, "%s(s,t)=%s;", XNAME, argv[3]);
72			scompile(NULL, stmp);
73			sprintf(stmp, "%s(s,t)=%s;", YNAME, argv[4]);
74			scompile(NULL, stmp);
75			sprintf(stmp, "%s(s,t)=%s;", ZNAME, argv[5]);
76			scompile(NULL, stmp);
77			m = atoi(argv[6]);
78			n = atoi(argv[7]);
79			if (m <= 0 \|\| n <= 0)
80			goto userror;
81
82	greg	1.4	row0 = (POINT )malloc((n+3)sizeof(POINT));
83			row1 = (POINT )malloc((n+3)sizeof(POINT));
84			row2 = (POINT )malloc((n+3)sizeof(POINT));
85	greg	1.3	if (row0 == NULL \|\| row1 == NULL \|\| row2 == NULL) {
86	greg	1.1	fprintf(stderr, "%s: out of memory\n", argv[0]);
87			quit(1);
88			}
89	greg	1.4	row0++; row1++; row2++;
90	greg	1.3	/* print header */
91	greg	1.1	printhead(argc, argv);
92	greg	1.4	/* initialize */
93			comprow(-1.0/m, row0, n);
94	greg	1.3	comprow(0.0, row1, n);
95			comprow(1.0/m, row2, n);
96	greg	1.4	compnorms(row0, row1, row2, n);
97	greg	1.3	/* for each row */
98	greg	1.1	for (i = 0; i < m; i++) {
99			/* compute next row */
100	greg	1.3	rp = row0;
101	greg	1.1	row0 = row1;
102	greg	1.3	row1 = row2;
103			row2 = rp;
104	greg	1.4	comprow((double)(i+2)/m, row2, n);
105			compnorms(row0, row1, row2, n);
106	greg	1.1
107			for (j = 0; j < n; j++) {
108	greg	1.3	/* put polygons */
109			if ((i+j) & 1)
110			putsquare(&row0[j], &row1[j],
111			&row0[j+1], &row1[j+1]);
112			else
113			putsquare(&row1[j], &row1[j+1],
114			&row0[j], &row0[j+1]);
115	greg	1.1	}
116			}
117
118			quit(0);
119
120			userror:
121			fprintf(stderr, "Usage: %s material name ", argv[0]);
122	greg	1.3	fprintf(stderr, "x(s,t) y(s,t) z(s,t) m n [-s][-e expr][-f file]\n");
123	greg	1.1	quit(1);
124			}
125
126
127	greg	1.3	putsquare(p0, p1, p2, p3) /* put out a square */
128			POINT p0, p1, p2, p3;
129			{
130			static int nout = 0;
131			FVECT norm[4];
132			int axis;
133			FVECT v1, v2, vc1, vc2;
134			int ok1, ok2;
135			/* compute exact normals */
136			fvsum(v1, p1->p, p0->p, -1.0);
137			fvsum(v2, p2->p, p0->p, -1.0);
138			fcross(vc1, v1, v2);
139			ok1 = normalize(vc1) != 0.0;
140			fvsum(v1, p2->p, p3->p, -1.0);
141			fvsum(v2, p1->p, p3->p, -1.0);
142			fcross(vc2, v1, v2);
143			ok2 = normalize(vc2) != 0.0;
144			if (!(ok1 \| ok2))
145			return;
146			/* compute normal interpolation */
147			axis = norminterp(norm, p0, p1, p2, p3);
148
149			/* put out quadrilateral? */
150			if (ok1 & ok2 && fdot(vc1,vc2) >= 1.0-FTINY*FTINY) {
151			printf("\n%s ", modname);
152			if (axis != -1) {
153			printf("texfunc %s\n", texname);
154			printf(tsargs);
155			printf("0\n13\t%d\n", axis);
156			pvect(norm[0]);
157			pvect(norm[1]);
158			pvect(norm[2]);
159			fvsum(v1, norm[3], vc1, -0.5);
160			fvsum(v1, v1, vc2, -0.5);
161			pvect(v1);
162			printf("\n%s ", texname);
163			}
164			printf("polygon %s.%d\n", surfname, ++nout);
165			printf("0\n0\n12\n");
166			pvect(p0->p);
167			pvect(p1->p);
168			pvect(p3->p);
169			pvect(p2->p);
170			return;
171			}
172			/* put out triangles? */
173			if (ok1) {
174			printf("\n%s ", modname);
175			if (axis != -1) {
176			printf("texfunc %s\n", texname);
177			printf(tsargs);
178			printf("0\n13\t%d\n", axis);
179			pvect(norm[0]);
180			pvect(norm[1]);
181			pvect(norm[2]);
182			fvsum(v1, norm[3], vc1, -1.0);
183			pvect(v1);
184			printf("\n%s ", texname);
185			}
186			printf("polygon %s.%d\n", surfname, ++nout);
187			printf("0\n0\n9\n");
188			pvect(p0->p);
189			pvect(p1->p);
190			pvect(p2->p);
191			}
192			if (ok2) {
193			printf("\n%s ", modname);
194			if (axis != -1) {
195			printf("texfunc %s\n", texname);
196			printf(tsargs);
197			printf("0\n13\t%d\n", axis);
198			pvect(norm[0]);
199			pvect(norm[1]);
200			pvect(norm[2]);
201			fvsum(v2, norm[3], vc2, -1.0);
202			pvect(v2);
203			printf("\n%s ", texname);
204			}
205			printf("polygon %s.%d\n", surfname, ++nout);
206			printf("0\n0\n9\n");
207			pvect(p2->p);
208			pvect(p1->p);
209			pvect(p3->p);
210			}
211			}
212
213
214	greg	1.1	comprow(s, row, siz) /* compute row of values */
215			double s;
216	greg	1.3	register POINT *row;
217	greg	1.1	int siz;
218			{
219	greg	1.4	double st[2];
220			register int i;
221			/* compute one past each end */
222	greg	1.1	st[0] = s;
223	greg	1.4	for (i = -1; i <= siz+1; i++) {
224			st[1] = (double)i/siz;
225			row[i].p[0] = funvalue(XNAME, 2, st);
226			row[i].p[1] = funvalue(YNAME, 2, st);
227			row[i].p[2] = funvalue(ZNAME, 2, st);
228	greg	1.1	}
229	greg	1.3	}
230
231
232			compnorms(r0, r1, r2, siz) /* compute row of averaged normals */
233			register POINT r0, r1, *r2;
234			int siz;
235			{
236			FVECT v1, v2, vc;
237	greg	1.4	register int i;
238	greg	1.3
239			if (!smooth) /* not needed if no smoothing */
240			return;
241			/* compute middle points */
242	greg	1.4	while (siz-- >= 0) {
243	greg	1.3	fvsum(v1, r2[0].p, r1[0].p, -1.0);
244			fvsum(v2, r1[1].p, r1[0].p, -1.0);
245			fcross(r1[0].n, v1, v2);
246			fvsum(v1, r0[0].p, r1[0].p, -1.0);
247			fcross(vc, v2, v1);
248			fvsum(r1[0].n, r1[0].n, vc, 1.0);
249			fvsum(v2, r1[-1].p, r1[0].p, -1.0);
250			fcross(vc, v1, v2);
251			fvsum(r1[0].n, r1[0].n, vc, 1.0);
252			fvsum(v1, r2[0].p, r1[0].p, -1.0);
253			fcross(vc, v2, v1);
254			fvsum(r1[0].n, r1[0].n, vc, 1.0);
255			normalize(r1[0].n);
256			r0++; r1++; r2++;
257			}
258			}
259
260
261			int
262			norminterp(resmat, p0, p1, p2, p3) /* compute normal interpolation */
263			register FVECT resmat[4];
264			POINT p0, p1, p2, p3;
265			{
266			#define u ((ax+1)%3)
267			#define v ((ax+2)%3)
268
269			register int ax;
270	greg	1.4	double eqnmat[4][4];
271	greg	1.3	FVECT v1;
272			register int i, j;
273
274			if (!smooth) /* no interpolation if no smoothing */
275			return(-1);
276			/* find dominant axis */
277			VCOPY(v1, p0->n);
278			fvsum(v1, v1, p1->n, 1.0);
279			fvsum(v1, v1, p2->n, 1.0);
280			fvsum(v1, v1, p3->n, 1.0);
281	greg	1.4	ax = ABS(v1[0]) > ABS(v1[1]) ? 0 : 1;
282			ax = ABS(v1[ax]) > ABS(v1[2]) ? ax : 2;
283	greg	1.3	/* assign equation matrix */
284			eqnmat[0][0] = p0->p[u]*p0->p[v];
285			eqnmat[0][1] = p0->p[u];
286			eqnmat[0][2] = p0->p[v];
287			eqnmat[0][3] = 1.0;
288			eqnmat[1][0] = p1->p[u]*p1->p[v];
289			eqnmat[1][1] = p1->p[u];
290			eqnmat[1][2] = p1->p[v];
291			eqnmat[1][3] = 1.0;
292			eqnmat[2][0] = p2->p[u]*p2->p[v];
293			eqnmat[2][1] = p2->p[u];
294			eqnmat[2][2] = p2->p[v];
295			eqnmat[2][3] = 1.0;
296			eqnmat[3][0] = p3->p[u]*p3->p[v];
297			eqnmat[3][1] = p3->p[u];
298			eqnmat[3][2] = p3->p[v];
299			eqnmat[3][3] = 1.0;
300			/* invert matrix (solve system) */
301	greg	1.4	if (!invmat(eqnmat, eqnmat))
302	greg	1.3	return(-1); /* no solution */
303			/* compute result matrix */
304			for (j = 0; j < 4; j++)
305			for (i = 0; i < 3; i++)
306	greg	1.4	resmat[j][i] = eqnmat[j][0]*p0->n[i] +
307			eqnmat[j][1]*p1->n[i] +
308			eqnmat[j][2]*p2->n[i] +
309			eqnmat[j][3]*p3->n[i];
310	greg	1.3	return(ax);
311
312			#undef u
313			#undef v
314			}
315
316
317			/*
318			* invmat - computes the inverse of mat into inverse. Returns 1
319			* if there exists an inverse, 0 otherwise. It uses Gaussian Elimination
320			* method.
321			*/
322
323			invmat(inverse,mat)
324			double mat[4][4],inverse[4][4];
325			{
326			#define SWAP(a,b,t) (t=a,a=b,b=t)
327
328	greg	1.4	double m4tmp[4][4];
329	greg	1.3	register int i,j,k;
330			register double temp;
331
332	greg	1.5	bcopy((char )mat, (char )m4tmp, sizeof(m4tmp));
333	greg	1.4	/* set inverse to identity */
334			for (i = 0; i < 4; i++)
335			for (j = 0; j < 4; j++)
336			inverse[i][j] = i==j ? 1.0 : 0.0;
337	greg	1.3
338			for(i = 0; i < 4; i++) {
339	greg	1.4	/* Look for raw with largest pivot and swap raws */
340			temp = FTINY; j = -1;
341			for(k = i; k < 4; k++)
342			if(ABS(m4tmp[k][i]) > temp) {
343			temp = ABS(m4tmp[k][i]);
344			j = k;
345			}
346			if(j == -1) /* No replacing raw -> no inverse */
347			return(0);
348			if (j != i)
349			for(k = 0; k < 4; k++) {
350			SWAP(m4tmp[i][k],m4tmp[j][k],temp);
351			SWAP(inverse[i][k],inverse[j][k],temp);
352			}
353	greg	1.3
354			temp = m4tmp[i][i];
355			for(k = 0; k < 4; k++) {
356			m4tmp[i][k] /= temp;
357			inverse[i][k] /= temp;
358			}
359			for(j = 0; j < 4; j++) {
360			if(j != i) {
361			temp = m4tmp[j][i];
362			for(k = 0; k < 4; k++) {
363			m4tmp[j][k] -= m4tmp[i][k]*temp;
364			inverse[j][k] -= inverse[i][k]*temp;
365			}
366			}
367			}
368			}
369			return(1);
370	greg	1.4
371	greg	1.3	#undef SWAP
372	greg	1.1	}
373
374
375			eputs(msg)
376			char *msg;
377			{
378			fputs(msg, stderr);
379			}
380
381
382			wputs(msg)
383			char *msg;
384			{
385			eputs(msg);
386			}
387
388
389			quit(code)
390			{
391			exit(code);
392			}
393
394
395			printhead(ac, av) /* print command header */
396			register int ac;
397			register char **av;
398			{
399			putchar('#');
400			while (ac--) {
401			putchar(' ');
402			fputs(*av++, stdout);
403			}
404			putchar('\n');
405			}
406
407
408			double
409			l_hermite()
410			{
411			double t;
412
413			t = argument(5);
414			return( argument(1)((2.0t-3.0)tt+1.0) +
415			argument(2)(-2.0t+3.0)tt +
416			argument(3)((t-2.0)t+1.0)*t +
417			argument(4)(t-1.0)t*t );
418			}