ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/rt/fprism.c
(Generate patch)

Comparing ray/src/rt/fprism.c (file contents):
Revision 2.1 by greg, Wed Sep 29 10:40:31 1993 UTC vs.
Revision 2.4 by greg, Sat Feb 22 02:07:28 2003 UTC

# Line 1 | Line 1
1 – /* Copyright (c) 1993 Regents of the University of California */
2 –
1   #ifndef lint
2 < static char SCCSid[] = "$SunId$ LBL";
2 > static const char       RCSid[] = "$Id$";
3   #endif
6 –
4   /* Ce programme calcule les directions et les energies des rayons lumineux
5     resultant du passage d'un rayon au travers d'un vitrage prismatique
6  
# Line 24 | Line 21 | double x;
21                  return(0.);
22          return(sqrt(x));
23   }
27 – #define sqrt(x) Sqrt(x)
24  
25   /* definitions de macros utiles */
26  
# Line 103 | Line 99 | extern long eclock;
99  
100   /* Definition des routines */
101  
102 < #define term(a,b) a/sqrt(a*a+b*b)
102 > #define term(a,b) a/Sqrt(a*a+b*b)
103   static
104   prepare_matrices()
105   {
# Line 209 | Line 205 | FVECT v,v_out;
205   FVECT v_temp;
206   double det;
207  
208 < det = sqrt( (YY(v)*YY(v)+ZZ(v)*ZZ(v))*(YY(v)*YY(v)+ZZ(v)*ZZ(v))+
208 > det = Sqrt( (YY(v)*YY(v)+ZZ(v)*ZZ(v))*(YY(v)*YY(v)+ZZ(v)*ZZ(v))+
209           (XX(v)*XX(v)*YY(v)*YY(v))+(XX(v)*XX(v)*ZZ(v)*ZZ(v)) );
210   XX(v_temp) = (YY(v)*YY(v)+ZZ(v)*ZZ(v))/det;
211   YY(v_temp) = -( XX(v)*YY(v) )/det;
# Line 227 | Line 223 | FVECT v,v_out;
223   FVECT v_temp;
224   double det;
225  
226 < det = sqrt( (ZZ(v)*ZZ(v)+YY(v)*YY(v)) );
226 > det = Sqrt( (ZZ(v)*ZZ(v)+YY(v)*YY(v)) );
227   XX(v_temp) = 0.;
228   YY(v_temp) = -ZZ(v)/det;
229   ZZ(v_temp) = YY(v)/det;
# Line 367 | Line 363 | TRAYON r_initial;
363  
364  
365  
366 + static int
367 + compare(r1,r2,marge)
368 + TRAYON r1, r2;
369 + double marge;
370  
371 + {
372 + double arctg1, arctg2;
373 +
374 + arctg1 = atan2(Y(r1),X(r1));
375 + arctg2 = atan2(Y(r2),X(r2));
376 + if ((arctg1 - marge <= arctg2) && (arctg1 + marge >= arctg2)) return 1;
377 + else return 0;
378 + }
379 +
380 +
381 +
382 +
383   static
384   sortie(r)
385   TRAYON r;
# Line 411 | Line 423 | TRAYON r;
423   {
424   double det;
425  
426 < det = sqrt(X(r)*X(r)+Y(r)*Y(r)+Z(r)*Z(r));
427 < sinus = sqrt(Y(r)*Y(r)+Z(r)*Z(r))/det;
428 < cosinus = sqrt(X(r)*X(r))/det;
426 > det = Sqrt(X(r)*X(r)+Y(r)*Y(r)+Z(r)*Z(r));
427 > sinus = Sqrt(Y(r)*Y(r)+Z(r)*Z(r))/det;
428 > cosinus = Sqrt(X(r)*X(r))/det;
429   if (r.n == 1.) rapport = prism.np * prism.np;
430   else rapport = 1./(prism.np * prism.np);
431   return;
# Line 432 | Line 444 | TRAYON r_incident;
444   X(r_reflechi) = -X(r_incident);
445   Y(r_reflechi) = Y(r_incident);
446   Z(r_reflechi) = Z(r_incident);
447 < if(sinus > sqrt(rapport) || r_incident.dest == tot_ref)
447 > if(sinus > Sqrt(rapport) || r_incident.dest == tot_ref)
448          {
449           r_reflechi.ppar1 = r_incident.ppar1;
450           r_reflechi.pper1 = r_incident.pper1;
# Line 442 | Line 454 | TRAYON r_incident;
454          }
455   else
456          {
457 <         r_reflechi.ppar1 = r_incident.ppar1*(rapport*cosinus-sqrt(rapport-
458 <                 (sinus*sinus)))/(rapport*cosinus+sqrt(rapport-(sinus*sinus)));
459 <         r_reflechi.pper1 = r_incident.pper1*(cosinus-sqrt
460 <                (rapport-(sinus*sinus)))/(cosinus+sqrt(rapport-(sinus*sinus)));
461 <         r_reflechi.ppar2 = r_incident.ppar2*(rapport*cosinus-sqrt(rapport-
462 <                 (sinus*sinus)))/(rapport*cosinus+sqrt(rapport-(sinus*sinus)));
463 <         r_reflechi.pper2 = r_incident.pper2*(cosinus-sqrt
464 <                (rapport-(sinus*sinus)))/(cosinus+sqrt(rapport-(sinus*sinus)));
457 >         r_reflechi.ppar1 = r_incident.ppar1*(rapport*cosinus-Sqrt(rapport-
458 >                 (sinus*sinus)))/(rapport*cosinus+Sqrt(rapport-(sinus*sinus)));
459 >         r_reflechi.pper1 = r_incident.pper1*(cosinus-Sqrt
460 >                (rapport-(sinus*sinus)))/(cosinus+Sqrt(rapport-(sinus*sinus)));
461 >         r_reflechi.ppar2 = r_incident.ppar2*(rapport*cosinus-Sqrt(rapport-
462 >                 (sinus*sinus)))/(rapport*cosinus+Sqrt(rapport-(sinus*sinus)));
463 >         r_reflechi.pper2 = r_incident.pper2*(cosinus-Sqrt
464 >                (rapport-(sinus*sinus)))/(cosinus+Sqrt(rapport-(sinus*sinus)));
465           r_reflechi.e = r_incident.e *(((r_reflechi.ppar1*r_reflechi.ppar1+
466           r_reflechi.pper1*r_reflechi.pper1)/(r_incident.ppar1*r_incident.ppar1+
467           r_incident.pper1*r_incident.pper1))+((r_reflechi.ppar2*r_reflechi.ppar2
# Line 471 | Line 483 | TRAYON r_incident;
483  
484   r_transmis = r_incident;
485   trigo(r_incident);
486 < if (sinus <= sqrt(rapport) && r_incident.dest != tot_ref)
486 > if (sinus <= Sqrt(rapport) && r_incident.dest != tot_ref)
487          {
488           X(r_transmis) = (X(r_incident)/(fabs(X(r_incident))))*
489 <                         (sqrt(1.-(Y(r_incident)*Y(r_incident)+Z(r_incident)*
489 >                         (Sqrt(1.-(Y(r_incident)*Y(r_incident)+Z(r_incident)*
490                                     Z(r_incident))/rapport));
491 <         Y(r_transmis) = Y(r_incident)/sqrt(rapport);
492 <         Z(r_transmis) = Z(r_incident)/sqrt(rapport);
493 <         r_transmis.ppar1 = r_incident.ppar1*2.*sqrt(rapport)*cosinus/
494 <                           (sqrt(rapport-sinus*sinus)+rapport*cosinus);
495 <         r_transmis.pper1 = r_incident.pper1*2.*cosinus/(cosinus+sqrt(rapport
491 >         Y(r_transmis) = Y(r_incident)/Sqrt(rapport);
492 >         Z(r_transmis) = Z(r_incident)/Sqrt(rapport);
493 >         r_transmis.ppar1 = r_incident.ppar1*2.*Sqrt(rapport)*cosinus/
494 >                           (Sqrt(rapport-sinus*sinus)+rapport*cosinus);
495 >         r_transmis.pper1 = r_incident.pper1*2.*cosinus/(cosinus+Sqrt(rapport
496                             - sinus*sinus));
497 <         r_transmis.ppar2 = r_incident.ppar2*2.*sqrt(rapport)*cosinus/
498 <                           (sqrt(rapport-sinus*sinus)+rapport*cosinus);
499 <         r_transmis.pper2 = r_incident.pper2*2.*cosinus/(cosinus+sqrt(rapport
497 >         r_transmis.ppar2 = r_incident.ppar2*2.*Sqrt(rapport)*cosinus/
498 >                           (Sqrt(rapport-sinus*sinus)+rapport*cosinus);
499 >         r_transmis.pper2 = r_incident.pper2*2.*cosinus/(cosinus+Sqrt(rapport
500                             - sinus*sinus));
501 <         r_transmis.e = (r_incident.e/2)*(sqrt(rapport-sinus*sinus)/cosinus)
501 >         r_transmis.e = (r_incident.e/2)*(Sqrt(rapport-sinus*sinus)/cosinus)
502              *(((r_transmis.ppar1*r_transmis.ppar1+r_transmis.pper1*
503                  r_transmis.pper1)
504              /(r_incident.ppar1*r_incident.ppar1+r_incident.pper1*
# Line 506 | Line 518 | TRAYON r_incident;
518  
519  
520  
509 – static int
510 – compare(r1,r2,marge)
511 – TRAYON r1, r2;
512 – double marge;
513 –
514 – {
515 – double arctg1, arctg2;
516 –
517 – arctg1 = atan2(Y(r1),X(r1));
518 – arctg2 = atan2(Y(r2),X(r2));
519 – if ((arctg1 - marge <= arctg2) && (arctg1 + marge >= arctg2)) return 1;
520 – else return 0;
521 – }
522 –
523 –
524 –
521   #define ensuite(rayon,prob_passage,destination) r_suite = rayon; \
522                                   r_suite.e = prob_passage(rayon)*rayon.e; \
523                                   r_suite.dest = destination; \
# Line 633 | Line 629 | TRAYON *r1,*r2;
629  
630  
631  
636 – static double
637 – l_get_val()
638 –
639 – {
640 – int val, dir, i, trouve, curseur;
641 – int nb;
642 – double valeur;
643 – TRAYON *rayt, raynull;
644 –
645 – if (prismclock < 0 || prismclock < eclock) setprism();
646 – if (bidon == BADVAL) {
647 –        errno = EDOM;
648 –        return(0.0);
649 – }
650 – val = (int)(argument(1) + .5);
651 – dir = (int)(argument(2) + .5);
652 – nb = (int)(argument(3) + .5);
653 – X(raynull) = bidon;
654 – Y(raynull) = Z(raynull) = 0.;
655 – raynull.e = 0.;
656 – trouve = curseur = 0;
657 – if ( !nosource && nb==2 ) nb=1; /* on est en train de tracer la source
658 –                                     a partir de sa seconde source virtuelle */
659 – #ifdef DEBUG
660 – fprintf(stderr, " On considere le rayon no: %d\n", nb);
661 – #endif
662 – for(i=0; i < nbrayons &&!trouve; i++)
663 –  {
664 –   if(ray[i].v[0] * dir * sens >= 0.) curseur ++;
665 –   if(curseur == nb)
666 –   {
667 –    rayt = &ray[i];
668 –    trouve = 1;
669 –   }
670 –  }
671 – if(!trouve) rayt = &raynull;
672 – switch(val) {
673 –        case 0 : valeur = rayt->v[0];
674 –                 break;
675 –        case 1 : valeur = rayt->v[1];
676 –                 break;
677 –        case 2 : valeur = rayt->v[2];
678 –                 break;
679 –        case 3 : valeur = rayt->e;
680 –                 break;
681 –        default : errno = EDOM; return(0.0);
682 –    }
683 – #ifdef DEBUG
684 –  fprintf(stderr, "get_val( %i, %i, %i) = %lf\n",val,dir,nb,valeur);
685 – #endif
686 –  return valeur;
687 – }
688 –
689 –
632   static
633   setprism()
692 –
634   {
635   double d;
636   TRAYON r_initial,rsource;
# Line 784 | Line 725 | if ( X(r_initial) != 0.)
725   bidon = BADVAL;
726   return;
727   }
728 +
729 +
730 + static double
731 + l_get_val()
732 +
733 + {
734 + int val, dir, i, trouve, curseur;
735 + int nb;
736 + double valeur;
737 + TRAYON *rayt, raynull;
738 +
739 + if (prismclock < 0 || prismclock < eclock) setprism();
740 + if (bidon == BADVAL) {
741 +        errno = EDOM;
742 +        return(0.0);
743 + }
744 + val = (int)(argument(1) + .5);
745 + dir = (int)(argument(2) + .5);
746 + nb = (int)(argument(3) + .5);
747 + X(raynull) = bidon;
748 + Y(raynull) = Z(raynull) = 0.;
749 + raynull.e = 0.;
750 + trouve = curseur = 0;
751 + if ( !nosource && nb==2 ) nb=1; /* on est en train de tracer la source
752 +                                     a partir de sa seconde source virtuelle */
753 + #ifdef DEBUG
754 + fprintf(stderr, " On considere le rayon no: %d\n", nb);
755 + #endif
756 + for(i=0; i < nbrayons &&!trouve; i++)
757 +  {
758 +   if(ray[i].v[0] * dir * sens >= 0.) curseur ++;
759 +   if(curseur == nb)
760 +   {
761 +    rayt = &ray[i];
762 +    trouve = 1;
763 +   }
764 +  }
765 + if(!trouve) rayt = &raynull;
766 + switch(val) {
767 +        case 0 : valeur = rayt->v[0];
768 +                 break;
769 +        case 1 : valeur = rayt->v[1];
770 +                 break;
771 +        case 2 : valeur = rayt->v[2];
772 +                 break;
773 +        case 3 : valeur = rayt->e;
774 +                 break;
775 +        default : errno = EDOM; return(0.0);
776 +    }
777 + #ifdef DEBUG
778 +  fprintf(stderr, "get_val( %i, %i, %i) = %lf\n",val,dir,nb,valeur);
779 + #endif
780 +  return valeur;
781 + }
782 +
783  
784   setprismfuncs()
785   {

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines