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  <stdio.h>

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

#define  PI             3.14159265358979323846

#define  FTINY          1e-7

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

char  vformat[] = "%15.9g %15.9g %15.9g\n";

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


main(argc, argv)
int  argc;
char  *argv[];
{
        static double  *xyz[4];
        double  *row0, *row1, *dp;
        double  v1[3], v2[3], vc1[3], vc2[3];
        double  a1, a2;
        int  i, j, m, n;
        char  stmp[256];
        double  d;
        register int  k;

        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
                        goto userror;

        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 = (double *)malloc((n+1)*3*sizeof(double));
        row1 = (double *)malloc((n+1)*3*sizeof(double));
        if (row0 == NULL || row1 == NULL) {
                fprintf(stderr, "%s: out of memory\n", argv[0]);
                quit(1);
        }

        printhead(argc, argv);

        comprow(0.0, row1, n);                  /* compute zeroeth row */

        for (i = 0; i < m; i++) {
                                                /* compute next row */
                dp = row0;
                row0 = row1;
                row1 = dp;
                comprow((double)(i+1)/m, row1, n);

                for (j = 0; j < n; j++) {
                                                        /* get vertices */
                        xyz[0] = row0 + 3*j;
                        xyz[1] = row1 + 3*j;
                        xyz[2] = xyz[0] + 3;
                        xyz[3] = xyz[1] + 3;
                                                        /* rotate vertices */
                        if (dist2(xyz[0],xyz[3]) < dist2(xyz[1],xyz[2])-FTINY) {
                                dp = xyz[0];
                                xyz[0] = xyz[1];
                                xyz[1] = xyz[3];
                                xyz[3] = xyz[2];
                                xyz[2] = dp;
                        }
                                                        /* get normals */
                        for (k = 0; k < 3; k++) {
                                v1[k] = xyz[1][k] - xyz[0][k];
                                v2[k] = xyz[2][k] - xyz[0][k];
                        }
                        fcross(vc1, v1, v2);
                        a1 = fdot(vc1, vc1);
                        for (k = 0; k < 3; k++) {
                                v1[k] = xyz[2][k] - xyz[3][k];
                                v2[k] = xyz[1][k] - xyz[3][k];
                        }
                        fcross(vc2, v1, v2);
                        a2 = fdot(vc2, vc2);
                                                        /* check coplanar */
                        if (a1 > FTINY*FTINY && a2 > FTINY*FTINY) {
                                d = fdot(vc1, vc2);
                                if (d*d/a1/a2 >= 1.0-FTINY*FTINY) {
                                        if (d > 0.0) {  /* coplanar */
                                                printf(
                                                "\n%s polygon %s.%d.%d\n",
                                                argv[1], argv[2], i+1, j+1);
                                                printf("0\n0\n12\n");
                                                vertex(xyz[0]);
                                                vertex(xyz[1]);
                                                vertex(xyz[3]);
                                                vertex(xyz[2]);
                                        }               /* else overlapped */
                                        continue;
                                }                       /* else bent */
                        }
                                                        /* check triangles */
                        if (a1 > FTINY*FTINY) {
                                printf("\n%s polygon %s.%da%d\n",
                                        argv[1], argv[2], i+1, j+1);
                                printf("0\n0\n9\n");
                                vertex(xyz[0]);
                                vertex(xyz[1]);
                                vertex(xyz[2]);
                        }
                        if (a2 > FTINY*FTINY) {
                                printf("\n%s polygon %s.%db%d\n",
                                        argv[1], argv[2], i+1, j+1);
                                printf("0\n0\n9\n");
                                vertex(xyz[2]);
                                vertex(xyz[1]);
                                vertex(xyz[3]);
                        }
                }
        }

        quit(0);

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


comprow(s, row, siz)                    /* compute row of values */
double  s;
register double  *row;
int  siz;
{
        double  st[2], step;

        st[0] = s;
        st[1] = 0.0;
        step = 1.0 / siz;
        while (siz-- >= 0) {
                *row++ = funvalue(XNAME, 2, st);
                *row++ = funvalue(YNAME, 2, st);
                *row++ = funvalue(ZNAME, 2, st);
                st[1] += step;
        }
}


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.2
Committed:	Sun Sep 10 16:45:16 1989 UTC (36 years, 1 month ago) by greg
Content type:	text/plain
Branch:	MAIN
Changes since 1.1:	+3 -1 lines
Log Message:	fixed location of SCCSid
#	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			#include <stdio.h>
20
21			#define XNAME "X_" /* x function name */
22			#define YNAME "Y_" /* y function name */
23			#define ZNAME "Z_" /* z function name */
24
25			#define PI 3.14159265358979323846
26
27			#define FTINY 1e-7
28
29			#define vertex(p) printf(vformat, (p)[0], (p)[1], (p)[2])
30
31			char vformat[] = "%15.9g %15.9g %15.9g\n";
32
33			double funvalue(), dist2(), fdot(), l_hermite(), argument();
34
35
36			main(argc, argv)
37			int argc;
38			char *argv[];
39			{
40			static double *xyz[4];
41			double row0, row1, *dp;
42			double v1[3], v2[3], vc1[3], vc2[3];
43			double a1, a2;
44			int i, j, m, n;
45			char stmp[256];
46			double d;
47			register int k;
48
49			varset("PI", PI);
50			funset("hermite", 5, l_hermite);
51
52			if (argc < 8)
53			goto userror;
54
55			for (i = 8; i < argc; i++)
56			if (!strcmp(argv[i], "-e"))
57			scompile(NULL, argv[++i]);
58			else if (!strcmp(argv[i], "-f"))
59			fcompile(argv[++i]);
60			else
61			goto userror;
62
63			sprintf(stmp, "%s(s,t)=%s;", XNAME, argv[3]);
64			scompile(NULL, stmp);
65			sprintf(stmp, "%s(s,t)=%s;", YNAME, argv[4]);
66			scompile(NULL, stmp);
67			sprintf(stmp, "%s(s,t)=%s;", ZNAME, argv[5]);
68			scompile(NULL, stmp);
69			m = atoi(argv[6]);
70			n = atoi(argv[7]);
71			if (m <= 0 \|\| n <= 0)
72			goto userror;
73
74			row0 = (double )malloc((n+1)3*sizeof(double));
75			row1 = (double )malloc((n+1)3*sizeof(double));
76			if (row0 == NULL \|\| row1 == NULL) {
77			fprintf(stderr, "%s: out of memory\n", argv[0]);
78			quit(1);
79			}
80
81			printhead(argc, argv);
82
83			comprow(0.0, row1, n); /* compute zeroeth row */
84
85			for (i = 0; i < m; i++) {
86			/* compute next row */
87			dp = row0;
88			row0 = row1;
89			row1 = dp;
90			comprow((double)(i+1)/m, row1, n);
91
92			for (j = 0; j < n; j++) {
93			/* get vertices */
94			xyz[0] = row0 + 3*j;
95			xyz[1] = row1 + 3*j;
96			xyz[2] = xyz[0] + 3;
97			xyz[3] = xyz[1] + 3;
98			/* rotate vertices */
99			if (dist2(xyz[0],xyz[3]) < dist2(xyz[1],xyz[2])-FTINY) {
100			dp = xyz[0];
101			xyz[0] = xyz[1];
102			xyz[1] = xyz[3];
103			xyz[3] = xyz[2];
104			xyz[2] = dp;
105			}
106			/* get normals */
107			for (k = 0; k < 3; k++) {
108			v1[k] = xyz[1][k] - xyz[0][k];
109			v2[k] = xyz[2][k] - xyz[0][k];
110			}
111			fcross(vc1, v1, v2);
112			a1 = fdot(vc1, vc1);
113			for (k = 0; k < 3; k++) {
114			v1[k] = xyz[2][k] - xyz[3][k];
115			v2[k] = xyz[1][k] - xyz[3][k];
116			}
117			fcross(vc2, v1, v2);
118			a2 = fdot(vc2, vc2);
119			/* check coplanar */
120			if (a1 > FTINYFTINY && a2 > FTINYFTINY) {
121			d = fdot(vc1, vc2);
122			if (dd/a1/a2 >= 1.0-FTINYFTINY) {
123			if (d > 0.0) { /* coplanar */
124			printf(
125			"\n%s polygon %s.%d.%d\n",
126			argv[1], argv[2], i+1, j+1);
127			printf("0\n0\n12\n");
128			vertex(xyz[0]);
129			vertex(xyz[1]);
130			vertex(xyz[3]);
131			vertex(xyz[2]);
132			} /* else overlapped */
133			continue;
134			} /* else bent */
135			}
136			/* check triangles */
137			if (a1 > FTINY*FTINY) {
138			printf("\n%s polygon %s.%da%d\n",
139			argv[1], argv[2], i+1, j+1);
140			printf("0\n0\n9\n");
141			vertex(xyz[0]);
142			vertex(xyz[1]);
143			vertex(xyz[2]);
144			}
145			if (a2 > FTINY*FTINY) {
146			printf("\n%s polygon %s.%db%d\n",
147			argv[1], argv[2], i+1, j+1);
148			printf("0\n0\n9\n");
149			vertex(xyz[2]);
150			vertex(xyz[1]);
151			vertex(xyz[3]);
152			}
153			}
154			}
155
156			quit(0);
157
158			userror:
159			fprintf(stderr, "Usage: %s material name ", argv[0]);
160			fprintf(stderr, "x(s,t) y(s,t) z(s,t) m n [-e expr] [-f file]\n");
161			quit(1);
162			}
163
164
165			comprow(s, row, siz) /* compute row of values */
166			double s;
167			register double *row;
168			int siz;
169			{
170			double st[2], step;
171
172			st[0] = s;
173			st[1] = 0.0;
174			step = 1.0 / siz;
175			while (siz-- >= 0) {
176			*row++ = funvalue(XNAME, 2, st);
177			*row++ = funvalue(YNAME, 2, st);
178			*row++ = funvalue(ZNAME, 2, st);
179			st[1] += step;
180			}
181			}
182
183
184			eputs(msg)
185			char *msg;
186			{
187			fputs(msg, stderr);
188			}
189
190
191			wputs(msg)
192			char *msg;
193			{
194			eputs(msg);
195			}
196
197
198			quit(code)
199			{
200			exit(code);
201			}
202
203
204			printhead(ac, av) /* print command header */
205			register int ac;
206			register char **av;
207			{
208			putchar('#');
209			while (ac--) {
210			putchar(' ');
211			fputs(*av++, stdout);
212			}
213			putchar('\n');
214			}
215
216
217			double
218			l_hermite()
219			{
220			double t;
221
222			t = argument(5);
223			return( argument(1)((2.0t-3.0)tt+1.0) +
224			argument(2)(-2.0t+3.0)tt +
225			argument(3)((t-2.0)t+1.0)*t +
226			argument(4)(t-1.0)t*t );
227			}