ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/src/gen/genprism.c
Revision: 2.15
Committed: Thu Aug 14 19:32:00 2025 UTC (7 weeks ago) by greg
Content type: text/plain
Branch: MAIN
CVS Tags: HEAD
Changes since 2.14: +32 -19 lines
Log Message:
fix(genprism): Added better self-checks

File Contents

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