ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/gen/genprism.c
Revision: 2.9
Committed: Sun Jun 8 12:03:09 2003 UTC (20 years, 9 months ago) by schorsch
Content type: text/plain
Branch: MAIN
Changes since 2.8: +171 -165 lines
Log Message:
Reduced compile warnings/errors on Windows.

File Contents

# Content
1 #ifndef lint
2 static const char RCSid[] = "$Id: genprism.c,v 2.8 2003/02/22 02:07:23 greg Exp $";
3 #endif
4 /*
5 * genprism.c - generate a prism.
6 * 2D vertices in the xy plane are given on the
7 * command line or from a file. Their order together
8 * with the extrude direction will determine surface
9 * orientation.
10 */
11
12 #include <stdio.h>
13 #include <string.h>
14 #include <stdlib.h>
15 #include <math.h>
16 #include <ctype.h>
17
18 #define MAXVERT 1024 /* maximum # vertices */
19
20 #define FTINY 1e-6
21
22 char *pmtype; /* material type */
23 char *pname; /* name */
24
25 double lvect[3] = {0.0, 0.0, 1.0};
26 int lvdir = 1;
27 double llen = 1.0;
28
29 double vert[MAXVERT][2];
30 int nverts = 0;
31
32 double u[MAXVERT][2]; /* edge unit vectors */
33 double a[MAXVERT]; /* corner trim sizes */
34
35 int do_ends = 1; /* include end caps */
36 int iscomplete = 0; /* polygon is already completed */
37 double crad = 0.0; /* radius for corners (sign from lvdir) */
38
39 extern double compute_rounding();
40
41
42 static void
43 readverts(fname) /* read vertices from a file */
44 char *fname;
45 {
46 FILE *fp;
47
48 if (fname == NULL)
49 fp = stdin;
50 else if ((fp = fopen(fname, "r")) == NULL) {
51 fprintf(stderr, "%s: cannot open\n", fname);
52 exit(1);
53 }
54 while (fscanf(fp, "%lf %lf", &vert[nverts][0], &vert[nverts][1]) == 2)
55 nverts++;
56 fclose(fp);
57 }
58
59
60 static void
61 side(n0, n1) /* print single side */
62 register int n0, n1;
63 {
64 printf("\n%s polygon %s.%d\n", pmtype, pname, n0+1);
65 printf("0\n0\n12\n");
66 printf("\t%18.12g\t%18.12g\t%18.12g\n", vert[n0][0],
67 vert[n0][1], 0.0);
68 printf("\t%18.12g\t%18.12g\t%18.12g\n", vert[n0][0]+lvect[0],
69 vert[n0][1]+lvect[1], lvect[2]);
70 printf("\t%18.12g\t%18.12g\t%18.12g\n", vert[n1][0]+lvect[0],
71 vert[n1][1]+lvect[1], lvect[2]);
72 printf("\t%18.12g\t%18.12g\t%18.12g\n", vert[n1][0],
73 vert[n1][1], 0.0);
74 }
75
76
77 static void
78 rside(n0, n1) /* print side with rounded edge */
79 register int n0, n1;
80 {
81 double s, c, t[3];
82
83 /* compute tanget offset vector */
84 s = lvdir*(lvect[1]*u[n1][0] - lvect[0]*u[n1][1])/llen;
85 if (s < -FTINY || s > FTINY) {
86 c = sqrt(1. - s*s);
87 t[0] = (c - 1.)*u[n1][1];
88 t[1] = (1. - c)*u[n1][0];
89 t[2] = s;
90 } else
91 t[0] = t[1] = t[2] = 0.;
92 /* output side */
93 printf("\n%s polygon %s.%d\n", pmtype, pname, n0+1);
94 printf("0\n0\n12\n");
95 printf("\t%18.12g\t%18.12g\t%18.12g\n",
96 vert[n0][0] + crad*(t[0] - a[n0]*u[n1][0]),
97 vert[n0][1] + crad*(t[1] - a[n0]*u[n1][1]),
98 crad*(t[2] + 1.));
99 printf("\t%18.12g\t%18.12g\t%18.12g\n",
100 vert[n0][0] + lvect[0] + crad*(t[0] - a[n0]*u[n1][0]),
101 vert[n0][1] + lvect[1] + crad*(t[1] - a[n0]*u[n1][1]),
102 lvect[2] + crad*(t[2] - 1.));
103 printf("\t%18.12g\t%18.12g\t%18.12g\n",
104 vert[n1][0] + lvect[0] + crad*(t[0] + a[n1]*u[n1][0]),
105 vert[n1][1] + lvect[1] + crad*(t[1] + a[n1]*u[n1][1]),
106 lvect[2] + crad*(t[2] - 1.));
107 printf("\t%18.12g\t%18.12g\t%18.12g\n",
108 vert[n1][0] + crad*(t[0] + a[n1]*u[n1][0]),
109 vert[n1][1] + crad*(t[1] + a[n1]*u[n1][1]),
110 crad*(t[2] + 1.));
111 /* output joining edge */
112 if (lvdir*a[n1] < 0.)
113 return;
114 printf("\n%s cylinder %s.e%d\n", pmtype, pname, n0+1);
115 printf("0\n0\n7\n");
116 printf("\t%18.12g\t%18.12g\t%18.12g\n",
117 vert[n1][0] + crad*(a[n1]*u[n1][0] - u[n1][1]),
118 vert[n1][1] + crad*(a[n1]*u[n1][1] + u[n1][0]),
119 crad);
120 printf("\t%18.12g\t%18.12g\t%18.12g\n",
121 vert[n1][0] + lvect[0] + crad*(a[n1]*u[n1][0] - u[n1][1]),
122 vert[n1][1] + lvect[1] + crad*(a[n1]*u[n1][1] + u[n1][0]),
123 lvect[2] - crad);
124 printf("\t%18.12g\n", lvdir*crad);
125 }
126
127
128 static double
129 compute_rounding() /* compute vectors for rounding operations */
130 {
131 register int i;
132 register double *v0, *v1;
133 double l, asum;
134
135 v0 = vert[nverts-1];
136 for (i = 0; i < nverts; i++) { /* compute u[*] */
137 v1 = vert[i];
138 u[i][0] = v0[0] - v1[0];
139 u[i][1] = v0[1] - v1[1];
140 l = sqrt(u[i][0]*u[i][0] + u[i][1]*u[i][1]);
141 if (l <= FTINY) {
142 fprintf(stderr, "Degenerate side in prism\n");
143 exit(1);
144 }
145 u[i][0] /= l;
146 u[i][1] /= l;
147 v0 = v1;
148 }
149 asum = 0.;
150 v1 = u[0];
151 for (i = nverts; i--; ) { /* compute a[*] */
152 v0 = u[i];
153 l = v0[0]*v1[0] + v0[1]*v1[1];
154 if (1+l <= FTINY) {
155 fprintf(stderr, "Overlapping sides in prism\n");
156 exit(1);
157 }
158 if (1-l <= 0.)
159 a[i] = 0.;
160 else {
161 a[i] = sqrt((1-l)/(1+l));
162 asum += l = v1[0]*v0[1]-v1[1]*v0[0];
163 if (l < 0.)
164 a[i] = -a[i];
165 }
166 v1 = v0;
167 }
168 return(asum*.5);
169 }
170
171
172 static void
173 printends() /* print ends of prism */
174 {
175 register int i;
176 /* bottom face */
177 printf("\n%s polygon %s.b\n", pmtype, pname);
178 printf("0\n0\n%d\n", nverts*3);
179 for (i = 0; i < nverts; i++)
180 printf("\t%18.12g\t%18.12g\t%18.12g\n", vert[i][0],
181 vert[i][1], 0.0);
182 /* top face */
183 printf("\n%s polygon %s.t\n", pmtype, pname);
184 printf("0\n0\n%d\n", nverts*3);
185 for (i = nverts; i--; )
186 printf("\t%18.12g\t%18.12g\t%18.12g\n", vert[i][0]+lvect[0],
187 vert[i][1]+lvect[1], lvect[2]);
188 }
189
190
191 static void
192 printrends() /* print ends of prism with rounding */
193 {
194 register int i;
195 double c0[3], c1[3], cl[3];
196 /* bottom face */
197 printf("\n%s polygon %s.b\n", pmtype, pname);
198 printf("0\n0\n%d\n", nverts*3);
199 for (i = 0; i < nverts; i++) {
200 printf("\t%18.12g\t%18.12g\t%18.12g\n",
201 vert[i][0] + crad*(a[i]*u[i][0] - u[i][1]),
202 vert[i][1] + crad*(a[i]*u[i][1] + u[i][0]),
203 0.0);
204 }
205 /* top face */
206 printf("\n%s polygon %s.t\n", pmtype, pname);
207 printf("0\n0\n%d\n", nverts*3);
208 for (i = nverts; i--; ) {
209 printf("\t%18.12g\t%18.12g\t%18.12g\n",
210 vert[i][0] + lvect[0] + crad*(a[i]*u[i][0] - u[i][1]),
211 vert[i][1] + lvect[1] + crad*(a[i]*u[i][1] + u[i][0]),
212 lvect[2]);
213 }
214 /* bottom corners and edges */
215 c0[0] = cl[0] = vert[nverts-1][0] +
216 crad*(a[nverts-1]*u[nverts-1][0] - u[nverts-1][1]);
217 c0[1] = cl[1] = vert[nverts-1][1] +
218 crad*(a[nverts-1]*u[nverts-1][1] + u[nverts-1][0]);
219 c0[2] = cl[2] = crad;
220 for (i = 0; i < nverts; i++) {
221 if (i < nverts-1) {
222 c1[0] = vert[i][0] + crad*(a[i]*u[i][0] - u[i][1]);
223 c1[1] = vert[i][1] + crad*(a[i]*u[i][1] + u[i][0]);
224 c1[2] = crad;
225 } else {
226 c1[0] = cl[0]; c1[1] = cl[1]; c1[2] = cl[2];
227 }
228 if (lvdir*a[i] > 0.) {
229 printf("\n%s sphere %s.bc%d\n", pmtype, pname, i+1);
230 printf("0\n0\n4 %18.12g %18.12g %18.12g %18.12g\n",
231 c1[0], c1[1], c1[2], lvdir*crad);
232 }
233 printf("\n%s cylinder %s.be%d\n", pmtype, pname, i+1);
234 printf("0\n0\n7\n");
235 printf("\t%18.12g\t%18.12g\t%18.12g\n", c0[0], c0[1], c0[2]);
236 printf("\t%18.12g\t%18.12g\t%18.12g\n", c1[0], c1[1], c1[2]);
237 printf("\t%18.12g\n", lvdir*crad);
238 c0[0] = c1[0]; c0[1] = c1[1]; c0[2] = c1[2];
239 }
240 /* top corners and edges */
241 c0[0] = cl[0] = vert[nverts-1][0] + lvect[0] +
242 crad*(a[nverts-1]*u[nverts-1][0] - u[nverts-1][1]);
243 c0[1] = cl[1] = vert[nverts-1][1] + lvect[1] +
244 crad*(a[nverts-1]*u[nverts-1][1] + u[nverts-1][0]);
245 c0[2] = cl[2] = lvect[2] - crad;
246 for (i = 0; i < nverts; i++) {
247 if (i < nverts-1) {
248 c1[0] = vert[i][0] + lvect[0] +
249 crad*(a[i]*u[i][0] - u[i][1]);
250 c1[1] = vert[i][1] + lvect[1] +
251 crad*(a[i]*u[i][1] + u[i][0]);
252 c1[2] = lvect[2] - crad;
253 } else {
254 c1[0] = cl[0]; c1[1] = cl[1]; c1[2] = cl[2];
255 }
256 if (lvdir*a[i] > 0.) {
257 printf("\n%s sphere %s.tc%d\n", pmtype, pname, i+1);
258 printf("0\n0\n4 %18.12g %18.12g %18.12g %18.12g\n",
259 c1[0], c1[1], c1[2], lvdir*crad);
260 }
261 printf("\n%s cylinder %s.te%d\n", pmtype, pname, i+1);
262 printf("0\n0\n7\n");
263 printf("\t%18.12g\t%18.12g\t%18.12g\n", c0[0], c0[1], c0[2]);
264 printf("\t%18.12g\t%18.12g\t%18.12g\n", c1[0], c1[1], c1[2]);
265 printf("\t%18.12g\n", lvdir*crad);
266 c0[0] = c1[0]; c0[1] = c1[1]; c0[2] = c1[2];
267 }
268 }
269
270
271 static void
272 printsides(round) /* print prism sides */
273 int round;
274 {
275 register int i;
276
277 for (i = 0; i < nverts-1; i++)
278 if (round)
279 rside(i, i+1);
280 else
281 side(i, i+1);
282 if (!iscomplete)
283 if (round)
284 rside(nverts-1, 0);
285 else
286 side(nverts-1, 0);
287 }
288
289
290 static void
291 printhead(ac, av) /* print command header */
292 register int ac;
293 register char **av;
294 {
295 putchar('#');
296 while (ac--) {
297 putchar(' ');
298 fputs(*av++, stdout);
299 }
300 putchar('\n');
301 }
302
303
304 main(argc, argv)
305 int argc;
306 char **argv;
307 {
308 int an;
309
310 if (argc < 4)
311 goto userr;
312
313 pmtype = argv[1];
314 pname = argv[2];
315
316 if (!strcmp(argv[3], "-")) {
317 readverts(NULL);
318 an = 4;
319 } else if (isdigit(argv[3][0])) {
320 nverts = atoi(argv[3]);
321 if (argc-3 < 2*nverts)
322 goto userr;
323 for (an = 0; an < nverts; an++) {
324 vert[an][0] = atof(argv[2*an+4]);
325 vert[an][1] = atof(argv[2*an+5]);
326 }
327 an = 2*nverts+4;
328 } else {
329 readverts(argv[3]);
330 an = 4;
331 }
332 if (nverts < 3) {
333 fprintf(stderr, "%s: not enough vertices\n", argv[0]);
334 exit(1);
335 }
336
337 for ( ; an < argc; an++) {
338 if (argv[an][0] != '-')
339 goto userr;
340 switch (argv[an][1]) {
341 case 'l': /* length vector */
342 lvect[0] = atof(argv[++an]);
343 lvect[1] = atof(argv[++an]);
344 lvect[2] = atof(argv[++an]);
345 if (lvect[2] < -FTINY)
346 lvdir = -1;
347 else if (lvect[2] > FTINY)
348 lvdir = 1;
349 else {
350 fprintf(stderr,
351 "%s: illegal extrusion vector\n",
352 argv[0]);
353 exit(1);
354 }
355 llen = sqrt(lvect[0]*lvect[0] + lvect[1]*lvect[1] +
356 lvect[2]*lvect[2]);
357 break;
358 case 'r': /* radius */
359 crad = atof(argv[++an]);
360 break;
361 case 'e': /* ends */
362 do_ends = !do_ends;
363 break;
364 case 'c': /* complete */
365 iscomplete = !iscomplete;
366 break;
367 default:
368 goto userr;
369 }
370 }
371 if (crad > FTINY) {
372 if (crad > lvdir*lvect[2]) {
373 fprintf(stderr, "%s: rounding greater than height\n",
374 argv[0]);
375 exit(1);
376 }
377 crad *= lvdir; /* simplifies formulas */
378 compute_rounding();
379 printhead(argc, argv);
380 if (do_ends)
381 printrends();
382 printsides(1);
383 } else {
384 printhead(argc, argv);
385 if (do_ends)
386 printends();
387 printsides(0);
388 }
389 exit(0);
390 userr:
391 fprintf(stderr, "Usage: %s material name ", argv[0]);
392 fprintf(stderr, "{ - | vfile | N v1 v2 .. vN } ");
393 fprintf(stderr, "[-l lvect][-r radius][-c][-e]\n");
394 exit(1);
395 }
396