src/gen/gensurf.c

/*

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