Back to index

radiance  4R0+20100331
genprism.c
Go to the documentation of this file.
00001 #ifndef lint
00002 static const char    RCSid[] = "$Id: genprism.c,v 2.12 2003/11/16 10:29:38 schorsch Exp $";
00003 #endif
00004 /*
00005  *  genprism.c - generate a prism.
00006  *            2D vertices in the xy plane are given on the
00007  *            command line or from a file.  Their order together
00008  *            with the extrude direction will determine surface
00009  *            orientation.
00010  */
00011 
00012 #include  <stdio.h>
00013 #include  <string.h>
00014 #include <stdlib.h>
00015 #include  <math.h>
00016 #include  <ctype.h>
00017 
00018 #define  MAXVERT     1024          /* maximum # vertices */
00019 
00020 #define  FTINY              1e-6
00021 
00022 char  *pmtype;              /* material type */
00023 char  *pname;        /* name */
00024 
00025 double  lvect[3] = {0.0, 0.0, 1.0};
00026 int    lvdir = 1;
00027 double llen = 1.0;
00028 
00029 double  vert[MAXVERT][2];
00030 int  nverts = 0;
00031 
00032 double u[MAXVERT][2];              /* edge unit vectors */
00033 double a[MAXVERT];          /* corner trim sizes */
00034 
00035 int  do_ends = 1;           /* include end caps */
00036 int  iscomplete = 0;        /* polygon is already completed */
00037 double  crad = 0.0;         /* radius for corners (sign from lvdir) */
00038 
00039 static double  compute_rounding(void);
00040 
00041 
00042 static void
00043 readverts(fname)            /* read vertices from a file */
00044 char  *fname;
00045 {
00046        FILE  *fp;
00047 
00048        if (fname == NULL)
00049               fp = stdin;
00050        else if ((fp = fopen(fname, "r")) == NULL) {
00051               fprintf(stderr, "%s: cannot open\n", fname);
00052               exit(1);
00053        }
00054        while (fscanf(fp, "%lf %lf", &vert[nverts][0], &vert[nverts][1]) == 2)
00055               nverts++;
00056        fclose(fp);
00057 }
00058 
00059 
00060 static void
00061 side(n0, n1)                /* print single side */
00062 register int  n0, n1;
00063 {
00064        printf("\n%s polygon %s.%d\n", pmtype, pname, n0+1);
00065        printf("0\n0\n12\n");
00066        printf("\t%18.12g\t%18.12g\t%18.12g\n", vert[n0][0],
00067                      vert[n0][1], 0.0);
00068        printf("\t%18.12g\t%18.12g\t%18.12g\n", vert[n0][0]+lvect[0],
00069                      vert[n0][1]+lvect[1], lvect[2]);
00070        printf("\t%18.12g\t%18.12g\t%18.12g\n", vert[n1][0]+lvect[0],
00071                      vert[n1][1]+lvect[1], lvect[2]);
00072        printf("\t%18.12g\t%18.12g\t%18.12g\n", vert[n1][0],
00073                      vert[n1][1], 0.0);
00074 }
00075 
00076 
00077 static void
00078 rside(n0, n1)               /* print side with rounded edge */
00079 register int  n0, n1;
00080 {
00081        double  s, c, t[3];
00082 
00083                                    /* compute tanget offset vector */
00084        s = lvdir*(lvect[1]*u[n1][0] - lvect[0]*u[n1][1])/llen;
00085        if (s < -FTINY || s > FTINY) {
00086               c = sqrt(1. - s*s);
00087               t[0] = (c - 1.)*u[n1][1];
00088               t[1] = (1. - c)*u[n1][0];
00089               t[2] = s;
00090        } else
00091               t[0] = t[1] = t[2] = 0.;
00092                                    /* output side */
00093        printf("\n%s polygon %s.%d\n", pmtype, pname, n0+1);
00094        printf("0\n0\n12\n");
00095        printf("\t%18.12g\t%18.12g\t%18.12g\n",
00096                      vert[n0][0] + crad*(t[0] - a[n0]*u[n1][0]),
00097                      vert[n0][1] + crad*(t[1] - a[n0]*u[n1][1]),
00098                      crad*(t[2] + 1.));
00099        printf("\t%18.12g\t%18.12g\t%18.12g\n",
00100                      vert[n0][0] + lvect[0] + crad*(t[0] - a[n0]*u[n1][0]),
00101                      vert[n0][1] + lvect[1] + crad*(t[1] - a[n0]*u[n1][1]),
00102                      lvect[2] + crad*(t[2] - 1.));
00103        printf("\t%18.12g\t%18.12g\t%18.12g\n",
00104                      vert[n1][0] + lvect[0] + crad*(t[0] + a[n1]*u[n1][0]),
00105                      vert[n1][1] + lvect[1] + crad*(t[1] + a[n1]*u[n1][1]),
00106                      lvect[2] + crad*(t[2] - 1.));
00107        printf("\t%18.12g\t%18.12g\t%18.12g\n",
00108                      vert[n1][0] + crad*(t[0] + a[n1]*u[n1][0]),
00109                      vert[n1][1] + crad*(t[1] + a[n1]*u[n1][1]),
00110                      crad*(t[2] + 1.));
00111                                    /* output joining edge */
00112        if (lvdir*a[n1] < 0.)
00113               return;
00114        printf("\n%s cylinder %s.e%d\n", pmtype, pname, n0+1);
00115        printf("0\n0\n7\n");
00116        printf("\t%18.12g\t%18.12g\t%18.12g\n",
00117               vert[n1][0] + crad*(a[n1]*u[n1][0] - u[n1][1]),
00118               vert[n1][1] + crad*(a[n1]*u[n1][1] + u[n1][0]),
00119               crad);
00120        printf("\t%18.12g\t%18.12g\t%18.12g\n",
00121               vert[n1][0] + lvect[0] + crad*(a[n1]*u[n1][0] - u[n1][1]),
00122               vert[n1][1] + lvect[1] + crad*(a[n1]*u[n1][1] + u[n1][0]),
00123               lvect[2] - crad);
00124        printf("\t%18.12g\n", lvdir*crad);
00125 }
00126 
00127 
00128 static double
00129 compute_rounding()          /* compute vectors for rounding operations */
00130 {
00131        register int  i;
00132        register double      *v0, *v1;
00133        double  l, asum;
00134 
00135        v0 = vert[nverts-1];
00136        for (i = 0; i < nverts; i++) {            /* compute u[*] */
00137               v1 = vert[i];
00138               u[i][0] = v0[0] - v1[0];
00139               u[i][1] = v0[1] - v1[1];
00140               l = sqrt(u[i][0]*u[i][0] + u[i][1]*u[i][1]);
00141               if (l <= FTINY) {
00142                      fprintf(stderr, "Degenerate side in prism\n");
00143                      exit(1);
00144               }
00145               u[i][0] /= l;
00146               u[i][1] /= l;
00147               v0 = v1;
00148        }
00149        asum = 0.;
00150        v1 = u[0];
00151        for (i = nverts; i--; ) {          /* compute a[*] */
00152               v0 = u[i];
00153               l = v0[0]*v1[0] + v0[1]*v1[1];
00154               if (1+l <= FTINY) {
00155                      fprintf(stderr, "Overlapping sides in prism\n");
00156                      exit(1);
00157               }
00158               if (1-l <= 0.)
00159                      a[i] = 0.;
00160               else {
00161                      a[i] = sqrt((1-l)/(1+l));
00162                      asum += l = v1[0]*v0[1]-v1[1]*v0[0];
00163                      if (l < 0.)
00164                             a[i] = -a[i];
00165               }
00166               v1 = v0;
00167        }
00168        return(asum*.5);
00169 }
00170 
00171 
00172 static void
00173 printends()                 /* print ends of prism */
00174 {
00175        register int  i;
00176                                           /* bottom face */
00177        printf("\n%s polygon %s.b\n", pmtype, pname);
00178        printf("0\n0\n%d\n", nverts*3);
00179        for (i = 0; i < nverts; i++)
00180               printf("\t%18.12g\t%18.12g\t%18.12g\n", vert[i][0],
00181                             vert[i][1], 0.0);
00182                                           /* top face */
00183        printf("\n%s polygon %s.t\n", pmtype, pname);
00184        printf("0\n0\n%d\n", nverts*3);
00185        for (i = nverts; i--; )
00186               printf("\t%18.12g\t%18.12g\t%18.12g\n", vert[i][0]+lvect[0],
00187                             vert[i][1]+lvect[1], lvect[2]);
00188 }
00189 
00190 
00191 static void
00192 printrends()                /* print ends of prism with rounding */
00193 {
00194        register int  i;
00195        double c0[3], c1[3], cl[3];
00196                                           /* bottom face */
00197        printf("\n%s polygon %s.b\n", pmtype, pname);
00198        printf("0\n0\n%d\n", nverts*3);
00199        for (i = 0; i < nverts; i++) {
00200               printf("\t%18.12g\t%18.12g\t%18.12g\n",
00201                      vert[i][0] + crad*(a[i]*u[i][0] - u[i][1]),
00202                      vert[i][1] + crad*(a[i]*u[i][1] + u[i][0]),
00203                      0.0);
00204        }
00205                                           /* top face */
00206        printf("\n%s polygon %s.t\n", pmtype, pname);
00207        printf("0\n0\n%d\n", nverts*3);
00208        for (i = nverts; i--; ) {
00209               printf("\t%18.12g\t%18.12g\t%18.12g\n",
00210               vert[i][0] + lvect[0] + crad*(a[i]*u[i][0] - u[i][1]),
00211               vert[i][1] + lvect[1] + crad*(a[i]*u[i][1] + u[i][0]),
00212               lvect[2]);
00213        }
00214                                           /* bottom corners and edges */
00215        c0[0] = cl[0] = vert[nverts-1][0] +
00216                      crad*(a[nverts-1]*u[nverts-1][0] - u[nverts-1][1]);
00217        c0[1] = cl[1] = vert[nverts-1][1] +
00218                      crad*(a[nverts-1]*u[nverts-1][1] + u[nverts-1][0]);
00219        c0[2] = cl[2] = crad;
00220        for (i = 0; i < nverts; i++) {
00221               if (i < nverts-1) {
00222                      c1[0] = vert[i][0] + crad*(a[i]*u[i][0] - u[i][1]);
00223                      c1[1] = vert[i][1] + crad*(a[i]*u[i][1] + u[i][0]);
00224                      c1[2] = crad;
00225               } else {
00226                      c1[0] = cl[0]; c1[1] = cl[1]; c1[2] = cl[2];
00227               }
00228               if (lvdir*a[i] > 0.) {
00229                      printf("\n%s sphere %s.bc%d\n", pmtype, pname, i+1);
00230                      printf("0\n0\n4 %18.12g %18.12g %18.12g %18.12g\n",
00231                             c1[0], c1[1], c1[2], lvdir*crad);
00232               }
00233               printf("\n%s cylinder %s.be%d\n", pmtype, pname, i+1);
00234               printf("0\n0\n7\n");
00235               printf("\t%18.12g\t%18.12g\t%18.12g\n", c0[0], c0[1], c0[2]);
00236               printf("\t%18.12g\t%18.12g\t%18.12g\n", c1[0], c1[1], c1[2]);
00237               printf("\t%18.12g\n", lvdir*crad);
00238               c0[0] = c1[0]; c0[1] = c1[1]; c0[2] = c1[2];
00239        }
00240                                           /* top corners and edges */
00241        c0[0] = cl[0] = vert[nverts-1][0] + lvect[0] +
00242                      crad*(a[nverts-1]*u[nverts-1][0] - u[nverts-1][1]);
00243        c0[1] = cl[1] = vert[nverts-1][1] + lvect[1] +
00244                      crad*(a[nverts-1]*u[nverts-1][1] + u[nverts-1][0]);
00245        c0[2] = cl[2] = lvect[2] - crad;
00246        for (i = 0; i < nverts; i++) {
00247               if (i < nverts-1) {
00248                      c1[0] = vert[i][0] + lvect[0] +
00249                                    crad*(a[i]*u[i][0] - u[i][1]);
00250                      c1[1] = vert[i][1] + lvect[1] +
00251                                    crad*(a[i]*u[i][1] + u[i][0]);
00252                      c1[2] = lvect[2] - crad;
00253               } else {
00254                      c1[0] = cl[0]; c1[1] = cl[1]; c1[2] = cl[2];
00255               }
00256               if (lvdir*a[i] > 0.) {
00257                      printf("\n%s sphere %s.tc%d\n", pmtype, pname, i+1);
00258                      printf("0\n0\n4 %18.12g %18.12g %18.12g %18.12g\n",
00259                             c1[0], c1[1], c1[2], lvdir*crad);
00260               }
00261               printf("\n%s cylinder %s.te%d\n", pmtype, pname, i+1);
00262               printf("0\n0\n7\n");
00263               printf("\t%18.12g\t%18.12g\t%18.12g\n", c0[0], c0[1], c0[2]);
00264               printf("\t%18.12g\t%18.12g\t%18.12g\n", c1[0], c1[1], c1[2]);
00265               printf("\t%18.12g\n", lvdir*crad);
00266               c0[0] = c1[0]; c0[1] = c1[1]; c0[2] = c1[2];
00267        }
00268 }
00269 
00270 
00271 static void
00272 printsides(round)           /* print prism sides */
00273 int  round;
00274 {
00275        register int  i;
00276 
00277        for (i = 0; i < nverts-1; i++)
00278               if (round)
00279                      rside(i, i+1);
00280               else
00281                      side(i, i+1);
00282        if (!iscomplete) {
00283               if (round)
00284                      rside(nverts-1, 0);
00285               else
00286                      side(nverts-1, 0);
00287        }
00288 }
00289 
00290 
00291 static void
00292 printhead(ac, av)           /* print command header */
00293 register int  ac;
00294 register char  **av;
00295 {
00296        putchar('#');
00297        while (ac--) {
00298               putchar(' ');
00299               fputs(*av++, stdout);
00300        }
00301        putchar('\n');
00302 }
00303 
00304 
00305 int
00306 main(argc, argv)
00307 int  argc;
00308 char  **argv;
00309 {
00310        int  an;
00311        
00312        if (argc < 4)
00313               goto userr;
00314 
00315        pmtype = argv[1];
00316        pname = argv[2];
00317 
00318        if (!strcmp(argv[3], "-")) {
00319               readverts(NULL);
00320               an = 4;
00321        } else if (isdigit(argv[3][0])) {
00322               nverts = atoi(argv[3]);
00323               if (argc-3 < 2*nverts)
00324                      goto userr;
00325               for (an = 0; an < nverts; an++) {
00326                      vert[an][0] = atof(argv[2*an+4]);
00327                      vert[an][1] = atof(argv[2*an+5]);
00328               }
00329               an = 2*nverts+4;
00330        } else {
00331               readverts(argv[3]);
00332               an = 4;
00333        }
00334        if (nverts < 3) {
00335               fprintf(stderr, "%s: not enough vertices\n", argv[0]);
00336               exit(1);
00337        }
00338 
00339        for ( ; an < argc; an++) {
00340               if (argv[an][0] != '-')
00341                      goto userr;
00342               switch (argv[an][1]) {
00343               case 'l':                          /* length vector */
00344                      lvect[0] = atof(argv[++an]);
00345                      lvect[1] = atof(argv[++an]);
00346                      lvect[2] = atof(argv[++an]);
00347                      if (lvect[2] < -FTINY)
00348                             lvdir = -1;
00349                      else if (lvect[2] > FTINY)
00350                             lvdir = 1;
00351                      else {
00352                             fprintf(stderr,
00353                                    "%s: illegal extrusion vector\n",
00354                                           argv[0]);
00355                             exit(1);
00356                      }
00357                      llen = sqrt(lvect[0]*lvect[0] + lvect[1]*lvect[1] +
00358                                    lvect[2]*lvect[2]);
00359                      break;
00360               case 'r':                          /* radius */
00361                      crad = atof(argv[++an]);
00362                      break;
00363               case 'e':                          /* ends */
00364                      do_ends = !do_ends;
00365                      break;
00366               case 'c':                          /* complete */
00367                      iscomplete = !iscomplete;
00368                      break;
00369               default:
00370                      goto userr;
00371               }
00372        }
00373        if (crad > FTINY) {
00374               if (crad > lvdir*lvect[2]) {
00375                      fprintf(stderr, "%s: rounding greater than height\n",
00376                                    argv[0]);
00377                      exit(1);
00378               }
00379               crad *= lvdir;              /* simplifies formulas */
00380               compute_rounding();
00381               printhead(argc, argv);
00382               if (do_ends)
00383                      printrends();
00384               printsides(1);
00385        } else {
00386               printhead(argc, argv);
00387               if (do_ends)
00388                      printends();
00389               printsides(0);
00390        }
00391        exit(0);
00392 userr:
00393        fprintf(stderr, "Usage: %s material name ", argv[0]);
00394        fprintf(stderr, "{ - | vfile | N v1 v2 .. vN } ");
00395        fprintf(stderr, "[-l lvect][-r radius][-c][-e]\n");
00396        exit(1);
00397 }
00398