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 |
|
|
21 |
|
return(0.); |
22 |
|
return(sqrt(x)); |
23 |
|
} |
27 |
– |
#define sqrt(x) Sqrt(x) |
24 |
|
|
25 |
|
/* definitions de macros utiles */ |
26 |
|
|
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 |
|
{ |
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; |
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; |
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; |
400 |
|
if (egalite == 0) |
401 |
|
{ |
402 |
|
if (nbrayons == 0) ray = (TRAYON *)calloc(1,sizeof(TRAYON)); |
403 |
< |
else ray = (TRAYON *)realloc(ray, (nbrayons+1)*(sizeof(TRAYON))); |
403 |
> |
else ray = (TRAYON *)realloc((void *)ray, (nbrayons+1)*(sizeof(TRAYON))); |
404 |
|
if (ray == NULL) |
405 |
|
error(SYSTEM, "out of memory in sortie\n"); |
406 |
|
raytemp = &ray[nbrayons]; |
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; |
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; |
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 |
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* |
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; \ |
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; |
725 |
|
bidon = BADVAL; |
726 |
|
return; |
727 |
|
} |
728 |
+ |
|
729 |
+ |
|
730 |
+ |
static double |
731 |
+ |
l_get_val(char *nm) |
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 |
|
{ |