ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/rt/fprism.c
Revision: 2.10
Committed: Tue Jul 8 18:25:00 2014 UTC (9 years, 9 months ago) by greg
Content type: text/plain
Branch: MAIN
CVS Tags: rad5R4, rad5R2, rad4R2P2, rad5R0, rad5R1, rad4R2, rad4R2P1, rad5R3, HEAD
Changes since 2.9: +2 -2 lines
Log Message:
Eliminated unnecessary "extern" and "register" modifiers

File Contents

# User Rev Content
1 greg 2.1 #ifndef lint
2 greg 2.10 static const char RCSid[] = "$Id: fprism.c,v 2.9 2013/08/07 05:10:09 greg Exp $";
3 greg 2.1 #endif
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    
7     1991, LESO - EPFL, R. Compagnon - F. Di Pasquale */
8    
9     #include "standard.h"
10    
11 greg 2.8 #include "ray.h"
12 schorsch 2.7 #include "calcomp.h"
13     #include "func.h"
14    
15 greg 2.1 #ifdef NOSTRUCTASSIGN
16     static double err = "No structure assignment!"; /* generate compiler error */
17     #endif
18    
19    
20     static double
21 schorsch 2.7 Sqrt(
22     double x
23     )
24 greg 2.1 {
25     if (x < 0.)
26     return(0.);
27     return(sqrt(x));
28     }
29    
30     /* definitions de macros utiles */
31    
32     #define ALPHA 0
33     #define BETA 1
34     #define GAMMA 2
35     #define DELTA 3
36     #define AUCUNE 4
37     #define X(r) r.v[0]
38     #define Y(r) r.v[1]
39     #define Z(r) r.v[2]
40     #define XX(v) v[0]
41     #define YY(v) v[1]
42     #define ZZ(v) v[2]
43     #define alpha_beta(v_alpha,v_beta) tfm(matbt,v_alpha,v_beta)
44     #define beta_alpha(v_beta,v_alpha) tfm(matb,v_beta,v_alpha)
45     #define alpha_gamma(v_alpha,v_gamma) tfm(matct,v_alpha,v_gamma)
46     #define gamma_alpha(v_gamma,v_alpha) tfm(matc,v_gamma,v_alpha)
47     #define prob_alpha_gamma(r) (1.-prob_alpha_beta(r))
48     #define prob_beta_gamma(r) (1.-prob_beta_alpha(r))
49     #define prob_gamma_beta(r) (1.-prob_gamma_alpha(r))
50     #define prob_delta_gamma(r) (1.-prob_delta_beta(r))
51     #define prob_beta_delta(r) (prob_beta_alpha(r))
52     #define prob_gamma_delta(r) (prob_gamma_alpha(r))
53     #define prob_delta_beta(r) (prob_alpha_beta(r))
54    
55    
56     /* Definitions des types de donnees */
57    
58     typedef struct { FVECT v; /* direction */
59     double ppar1,pper1,
60     ppar2,pper2; /* polarisations */
61     double e; /* energie */
62     double n; /* milieu dans lequel on se propage */
63     int orig,dest; /* origine et destination */
64     } TRAYON;
65    
66     typedef struct { double a,b,c,d; /* longueurs caracteristiques */
67     double np; /* indice de refraction */
68     } TPRISM;
69    
70    
71     /* Definitions des variables globales */
72    
73     static TPRISM prism; /* le prisme ! */
74     static MAT4 matb = MAT4IDENT; /* matrices de changement de bases */
75     static MAT4 matbt = MAT4IDENT;
76     static MAT4 matc = MAT4IDENT;
77     static MAT4 matct = MAT4IDENT;
78     static double seuil; /* seuil pour l'arret du trace */
79     static double sinus,cosinus; /* sin et cos */
80     static double rapport; /* rapport entre les indices des
81     milieux refracteur et incident */
82     static int tot_ref; /* flag pour les surfaces
83     reflechissantes */
84     static double fact_ref[4]={1.0,1.0,1.0,1.0}; /* facteurs de reflexion */
85     static double tolerance; /* degre de tol. pour les amalgames */
86     static double tolsource; /* degre de tol. pour les sources */
87     static int bidon;
88     #define BADVAL (-10)
89     static long prismclock = -1;
90     static int nosource; /* indique que l'on ne trace pas
91     en direction d'une source */
92     static int sens; /* indique le sens de prop. du ray.*/
93     static int nbrayons; /* indice des rayons sortants */
94     static TRAYON *ray; /* tableau des rayons sortants */
95     static TRAYON *raytemp; /* variable temporaire */
96    
97 schorsch 2.7
98     static void prepare_matrices(void);
99     static void tfm(MAT4 mat, FVECT v_old, FVECT v_new);
100     static double prob_alpha_beta(TRAYON r);
101     static double prob_beta_alpha(TRAYON r);
102     static double prob_gamma_alpha(TRAYON r);
103     static void v_par(FVECT v, FVECT v_out);
104     static void v_per(FVECT v, FVECT v_out);
105     static TRAYON transalphabeta(TRAYON r_initial);
106     static TRAYON transbetaalpha(TRAYON r_initial);
107     static TRAYON transalphagamma(TRAYON r_initial);
108     static TRAYON transgammaalpha(TRAYON r_initial);
109     static int compare(TRAYON r1, TRAYON r2, double marge);
110     static void sortie(TRAYON r);
111     static void trigo(TRAYON r);
112     static TRAYON reflexion(TRAYON r_incident);
113     static TRAYON transmission(TRAYON r_incident);
114     static void trace_rayon(TRAYON r_incident);
115     static void inverser(TRAYON *r1, TRAYON *r2);
116     static void setprism(void);
117     static double l_get_val(char *nm);
118 greg 2.1
119    
120     /* Definition des routines */
121    
122 greg 2.3 #define term(a,b) a/Sqrt(a*a+b*b)
123 schorsch 2.7 static void
124     prepare_matrices(void)
125 greg 2.1 {
126 schorsch 2.7 /* preparation des matrices de changement de bases */
127 greg 2.1
128 schorsch 2.7 matb[0][0] = matbt[0][0] = matb[1][1] = matbt[1][1] = term(prism.a,prism.d);
129     matb[1][0] = matbt[0][1] = term(-prism.d,prism.a);
130     matb[0][1] = matbt[1][0] = term(prism.d,prism.a);
131     matc[0][0] = matct[0][0] = matc[1][1] = matct[1][1] = term(prism.b,prism.d);
132     matc[1][0] = matct[0][1] = term(prism.d,prism.b);
133     matc[0][1] = matct[1][0] = term(-prism.d,prism.b);
134     return;
135 greg 2.1 }
136     #undef term
137    
138    
139 schorsch 2.7 static void
140     tfm(
141     MAT4 mat,
142     FVECT v_old,
143     FVECT v_new
144     )
145     {
146     /* passage d'un repere old au repere new par la matrice mat */
147     FVECT v_temp;
148    
149     multv3(v_temp,v_old,mat);
150     normalize(v_temp);
151     VCOPY(v_new,v_temp);
152     return;
153 greg 2.1 }
154    
155     #define A prism.a
156     #define B prism.b
157     #define C prism.c
158     #define D prism.d
159    
160    
161     static double
162 schorsch 2.7 prob_alpha_beta(
163     TRAYON r
164     )
165 greg 2.1 {
166 schorsch 2.7 /* calcul de la probabilite de passage de alpha a beta */
167     double prob,test;
168 greg 2.1
169 schorsch 2.7 if ( X(r) != 0. )
170 greg 2.1 {
171 schorsch 2.7 test = Y(r)/X(r);
172     if ( test > B/D ) prob = 1.;
173     else if ( test >= -A/D ) prob = (A+test*D)/(A+B);
174     else prob = 0.;
175 greg 2.1 }
176 schorsch 2.7 else prob = 0.;
177     return prob;
178 greg 2.1 }
179    
180    
181     static double
182 schorsch 2.7 prob_beta_alpha(
183     TRAYON r
184     )
185 greg 2.1 {
186 schorsch 2.7 /* calcul de la probabilite de passage de beta a aplha */
187     double prob,test;
188 greg 2.1
189 schorsch 2.7 if ( X(r) != 0. )
190 greg 2.1 {
191 schorsch 2.7 test = Y(r)/X(r);
192     if ( test > B/D ) prob = (A+B)/(A+test*D);
193     else if ( test >= -A/D ) prob = 1.;
194     else prob = 0.;
195 greg 2.1 }
196 schorsch 2.7 else prob = 0.;
197     return prob;
198 greg 2.1 }
199    
200    
201 schorsch 2.7 static double
202     prob_gamma_alpha(
203     TRAYON r
204     )
205 greg 2.1 {
206 schorsch 2.7 /* calcul de la probabilite de passage de gamma a alpha */
207     double prob,test;
208 greg 2.1
209 schorsch 2.7 if ( X(r) != 0. )
210 greg 2.1 {
211 schorsch 2.7 test = Y(r)/X(r);
212     if ( test > B/D ) prob = 0.;
213     else if ( test >= -A/D ) prob = 1.;
214     else prob = (A+B)/(B-test*D);
215 greg 2.1 }
216 schorsch 2.7 else prob = 0.;
217     return prob;
218 greg 2.1 }
219    
220     #undef A
221     #undef B
222     #undef C
223     #undef D
224    
225    
226 schorsch 2.7 static void
227     v_par(
228     FVECT v,
229     FVECT v_out
230     )
231 greg 2.1 /* calcule le vecteur par au plan d'incidence lie a v */
232     {
233 schorsch 2.7 FVECT v_temp;
234     double det;
235 greg 2.1
236 schorsch 2.7 det = Sqrt( (YY(v)*YY(v)+ZZ(v)*ZZ(v))*(YY(v)*YY(v)+ZZ(v)*ZZ(v))+
237     (XX(v)*XX(v)*YY(v)*YY(v))+(XX(v)*XX(v)*ZZ(v)*ZZ(v)) );
238     XX(v_temp) = (YY(v)*YY(v)+ZZ(v)*ZZ(v))/det;
239     YY(v_temp) = -( XX(v)*YY(v) )/det;
240     ZZ(v_temp) = -( XX(v)*ZZ(v) )/det;
241     VCOPY(v_out,v_temp);
242     return;
243 greg 2.1 }
244    
245    
246 schorsch 2.7 static void
247     v_per(
248     FVECT v,
249     FVECT v_out
250     )
251 greg 2.1 /* calcule le vecteur perp au plan d'incidence lie a v */
252     {
253 schorsch 2.7 FVECT v_temp;
254     double det;
255 greg 2.1
256 schorsch 2.7 det = Sqrt( (ZZ(v)*ZZ(v)+YY(v)*YY(v)) );
257     XX(v_temp) = 0.;
258     YY(v_temp) = -ZZ(v)/det;
259     ZZ(v_temp) = YY(v)/det;
260     VCOPY(v_out,v_temp);
261     return;
262 greg 2.1 }
263    
264    
265     static TRAYON
266 schorsch 2.7 transalphabeta(
267     TRAYON r_initial
268     )
269 greg 2.1 /* transforme le rayon r_initial de la base associee a alpha dans
270     la base associee a beta */
271     {
272 schorsch 2.7 TRAYON r_final;
273     FVECT vpar_temp1,vpar_temp2,vper_temp1,vper_temp2;
274 greg 2.1
275 schorsch 2.7 r_final = r_initial;
276     alpha_beta(r_initial.v,r_final.v);
277     if ((Y(r_initial) != 0. || Z(r_initial) != 0.)&&(Y(r_final) !=0. || Z(r_final)!= 0.))
278     {
279     v_par(r_initial.v,vpar_temp1);
280     alpha_beta(vpar_temp1,vpar_temp1);
281     v_per(r_initial.v,vper_temp1);
282     alpha_beta(vper_temp1,vper_temp1);
283     v_par(r_final.v,vpar_temp2);
284     v_per(r_final.v,vper_temp2);
285     r_final.ppar1 = (r_initial.ppar1*fdot(vpar_temp1,vpar_temp2))+
286     (r_initial.pper1*fdot(vper_temp1,vpar_temp2));
287     r_final.pper1 = (r_initial.ppar1*fdot(vpar_temp1,vper_temp2))+
288     (r_initial.pper1*fdot(vper_temp1,vper_temp2));
289     r_final.ppar2 = (r_initial.ppar2*fdot(vpar_temp1,vpar_temp2))+
290     (r_initial.pper2*fdot(vper_temp1,vpar_temp2));
291     r_final.pper2 = (r_initial.ppar2*fdot(vpar_temp1,vper_temp2))+
292     (r_initial.pper2*fdot(vper_temp1,vper_temp2));
293     }
294     return r_final;
295 greg 2.1 }
296    
297    
298     static TRAYON
299 schorsch 2.7 transbetaalpha(
300     TRAYON r_initial
301     )
302     {
303     /* transforme le rayon r_initial de la base associee a beta dans
304     la base associee a alpha */
305     TRAYON r_final;
306     FVECT vpar_temp1,vpar_temp2,vper_temp1,vper_temp2;
307    
308     r_final = r_initial;
309     beta_alpha(r_initial.v,r_final.v);
310     if ((Y(r_initial) != 0. || Z(r_initial) != 0. )&&(Y(r_final) != 0. || Z(r_final)!= 0.))
311     {
312     v_par(r_initial.v,vpar_temp1);
313     beta_alpha(vpar_temp1,vpar_temp1);
314     v_per(r_initial.v,vper_temp1);
315     beta_alpha(vper_temp1,vper_temp1);
316     v_par(r_final.v,vpar_temp2);
317     v_per(r_final.v,vper_temp2);
318     r_final.ppar1 = (r_initial.ppar1*fdot(vpar_temp1,vpar_temp2))+
319     (r_initial.pper1*fdot(vper_temp1,vpar_temp2));
320     r_final.pper1 = (r_initial.ppar1*fdot(vpar_temp1,vper_temp2))+
321     (r_initial.pper1*fdot(vper_temp1,vper_temp2));
322     r_final.ppar2 = (r_initial.ppar2*fdot(vpar_temp1,vpar_temp2))+
323     (r_initial.pper2*fdot(vper_temp1,vpar_temp2));
324     r_final.pper2 = (r_initial.ppar2*fdot(vpar_temp1,vper_temp2))+
325     (r_initial.pper2*fdot(vper_temp1,vper_temp2));
326 greg 2.1
327 schorsch 2.7 }
328     return r_final;
329 greg 2.1 }
330    
331    
332     static TRAYON
333 schorsch 2.7 transalphagamma(
334     TRAYON r_initial
335     )
336 greg 2.1 /* transforme le rayon r_initial de la base associee a alpha dans
337     la base associee a gamma */
338     {
339 schorsch 2.7 TRAYON r_final;
340     FVECT vpar_temp1,vpar_temp2,vper_temp1,vper_temp2;
341 greg 2.1
342 schorsch 2.7 r_final = r_initial;
343     alpha_gamma(r_initial.v,r_final.v);
344     if (( Y(r_initial) != 0. || Z(r_initial) != 0. )&&(Y(r_final)!= 0. || Z(r_final) !=0.))
345     {
346     v_par(r_initial.v,vpar_temp1);
347     alpha_gamma(vpar_temp1,vpar_temp1);
348     v_per(r_initial.v,vper_temp1);
349     alpha_gamma(vper_temp1,vper_temp1);
350     v_par(r_final.v,vpar_temp2);
351     v_per(r_final.v,vper_temp2);
352     r_final.ppar1 = (r_initial.ppar1*fdot(vpar_temp1,vpar_temp2))+
353     (r_initial.pper1*fdot(vper_temp1,vpar_temp2));
354     r_final.pper1 = (r_initial.ppar1*fdot(vpar_temp1,vper_temp2))+
355     (r_initial.pper1*fdot(vper_temp1,vper_temp2));
356     r_final.ppar2 = (r_initial.ppar2*fdot(vpar_temp1,vpar_temp2))+
357     (r_initial.pper2*fdot(vper_temp1,vpar_temp2));
358     r_final.pper2 = (r_initial.ppar2*fdot(vpar_temp1,vper_temp2))+
359     (r_initial.pper2*fdot(vper_temp1,vper_temp2));
360 greg 2.1
361 schorsch 2.7 }
362     return r_final;
363 greg 2.1 }
364    
365    
366     static TRAYON
367 schorsch 2.7 transgammaalpha(
368     TRAYON r_initial
369     )
370 greg 2.1 /* transforme le rayon r_initial de la base associee a gamma dans
371     la base associee a alpha */
372     {
373 schorsch 2.7 TRAYON r_final;
374     FVECT vpar_temp1,vpar_temp2,vper_temp1,vper_temp2;
375 greg 2.1
376 schorsch 2.7 r_final = r_initial;
377     gamma_alpha(r_initial.v,r_final.v);
378     if (( Y(r_initial) != 0. || Z(r_initial) != 0. )&&(Y(r_final) !=0. || Z(r_final) != 0.))
379     {
380     v_par(r_initial.v,vpar_temp1);
381     gamma_alpha(vpar_temp1,vpar_temp1);
382     v_per(r_initial.v,vper_temp1);
383     gamma_alpha(vper_temp1,vper_temp1);
384     v_par(r_final.v,vpar_temp2);
385     v_per(r_final.v,vper_temp2);
386     r_final.ppar1 = (r_initial.ppar1*fdot(vpar_temp1,vpar_temp2))+
387     (r_initial.pper1*fdot(vper_temp1,vpar_temp2));
388     r_final.pper1 = (r_initial.ppar1*fdot(vpar_temp1,vper_temp2))+
389     (r_initial.pper1*fdot(vper_temp1,vper_temp2));
390     r_final.ppar2 = (r_initial.ppar2*fdot(vpar_temp1,vpar_temp2))+
391     (r_initial.pper2*fdot(vper_temp1,vpar_temp2));
392     r_final.pper2 = (r_initial.ppar2*fdot(vpar_temp1,vper_temp2))+
393     (r_initial.pper2*fdot(vper_temp1,vper_temp2));
394     }
395     return r_final;
396 greg 2.1 }
397    
398    
399    
400 greg 2.2 static int
401 schorsch 2.7 compare(
402     TRAYON r1,
403     TRAYON r2,
404     double marge
405     )
406 greg 2.2 {
407 schorsch 2.7 double arctg1, arctg2;
408 greg 2.2
409 schorsch 2.7 arctg1 = atan2(Y(r1),X(r1));
410     arctg2 = atan2(Y(r2),X(r2));
411     if ((arctg1 - marge <= arctg2) && (arctg1 + marge >= arctg2)) return 1;
412     else return 0;
413 greg 2.2 }
414    
415    
416    
417    
418 schorsch 2.7 static void
419     sortie(
420     TRAYON r
421     )
422 greg 2.1 {
423 schorsch 2.7 int i = 0;
424     int egalite = 0;
425 greg 2.1
426    
427 schorsch 2.7 if(r.e > seuil)
428     {
429     while (i < nbrayons && egalite == 0)
430     {
431     raytemp = &ray[i];
432     egalite = compare(r,*raytemp,tolerance);
433     if (egalite) raytemp->e = raytemp->e + r.e;
434     else i = i + 1;
435     }
436     if (egalite == 0)
437     {
438     if (nbrayons == 0) ray = (TRAYON *)calloc(1,sizeof(TRAYON));
439     else ray = (TRAYON *)realloc((void *)ray, (nbrayons+1)*(sizeof(TRAYON)));
440     if (ray == NULL)
441     error(SYSTEM, "out of memory in sortie\n");
442     raytemp = &ray[nbrayons];
443     raytemp->v[0] = X(r);
444     raytemp->v[1] = Y(r);
445     raytemp->v[2] = Z(r);
446     raytemp->e = r.e;
447     nbrayons++;
448     }
449     }
450     return;
451 greg 2.1 }
452    
453    
454 schorsch 2.7 static void
455     trigo(
456     TRAYON r
457     )
458 greg 2.1 /* calcule les grandeurs trigonometriques relatives au rayon incident
459     et le rapport entre les indices du milieu refracteur et incident */
460     {
461 schorsch 2.7 double det;
462    
463     det = Sqrt(X(r)*X(r)+Y(r)*Y(r)+Z(r)*Z(r));
464     sinus = Sqrt(Y(r)*Y(r)+Z(r)*Z(r))/det;
465     cosinus = Sqrt(X(r)*X(r))/det;
466     if (r.n == 1.) rapport = prism.np * prism.np;
467     else rapport = 1./(prism.np * prism.np);
468     return;
469 greg 2.1 }
470    
471    
472     static TRAYON
473 schorsch 2.7 reflexion(
474     TRAYON r_incident
475     )
476 greg 2.1 {
477     /* calcul du rayon reflechi par une face */
478     TRAYON r_reflechi;
479    
480     r_reflechi = r_incident;
481     trigo(r_incident);
482     X(r_reflechi) = -X(r_incident);
483     Y(r_reflechi) = Y(r_incident);
484     Z(r_reflechi) = Z(r_incident);
485 greg 2.3 if(sinus > Sqrt(rapport) || r_incident.dest == tot_ref)
486 greg 2.1 {
487     r_reflechi.ppar1 = r_incident.ppar1;
488     r_reflechi.pper1 = r_incident.pper1;
489     r_reflechi.ppar2 = r_incident.ppar2;
490     r_reflechi.pper2 = r_incident.pper2;
491     r_reflechi.e = r_incident.e * fact_ref[r_incident.dest];
492     }
493     else
494     {
495 greg 2.3 r_reflechi.ppar1 = r_incident.ppar1*(rapport*cosinus-Sqrt(rapport-
496     (sinus*sinus)))/(rapport*cosinus+Sqrt(rapport-(sinus*sinus)));
497     r_reflechi.pper1 = r_incident.pper1*(cosinus-Sqrt
498     (rapport-(sinus*sinus)))/(cosinus+Sqrt(rapport-(sinus*sinus)));
499     r_reflechi.ppar2 = r_incident.ppar2*(rapport*cosinus-Sqrt(rapport-
500     (sinus*sinus)))/(rapport*cosinus+Sqrt(rapport-(sinus*sinus)));
501     r_reflechi.pper2 = r_incident.pper2*(cosinus-Sqrt
502     (rapport-(sinus*sinus)))/(cosinus+Sqrt(rapport-(sinus*sinus)));
503 greg 2.1 r_reflechi.e = r_incident.e *(((r_reflechi.ppar1*r_reflechi.ppar1+
504     r_reflechi.pper1*r_reflechi.pper1)/(r_incident.ppar1*r_incident.ppar1+
505     r_incident.pper1*r_incident.pper1))+((r_reflechi.ppar2*r_reflechi.ppar2
506     +r_reflechi.pper2*r_reflechi.pper2)/(r_incident.ppar2*r_incident.ppar2
507     +r_incident.pper2*r_incident.pper2)))/2;
508     }
509    
510     /* a la sortie de cette routine r_transmis.orig et .dest ne sont pas definis!*/
511     return r_reflechi;
512     }
513    
514    
515     static TRAYON
516 schorsch 2.7 transmission(
517     TRAYON r_incident
518     )
519 greg 2.1 {
520     /* calcul du rayon refracte par une face */
521     TRAYON r_transmis;
522    
523     r_transmis = r_incident;
524     trigo(r_incident);
525 greg 2.3 if (sinus <= Sqrt(rapport) && r_incident.dest != tot_ref)
526 greg 2.1 {
527     X(r_transmis) = (X(r_incident)/(fabs(X(r_incident))))*
528 greg 2.3 (Sqrt(1.-(Y(r_incident)*Y(r_incident)+Z(r_incident)*
529 greg 2.1 Z(r_incident))/rapport));
530 greg 2.3 Y(r_transmis) = Y(r_incident)/Sqrt(rapport);
531     Z(r_transmis) = Z(r_incident)/Sqrt(rapport);
532     r_transmis.ppar1 = r_incident.ppar1*2.*Sqrt(rapport)*cosinus/
533     (Sqrt(rapport-sinus*sinus)+rapport*cosinus);
534     r_transmis.pper1 = r_incident.pper1*2.*cosinus/(cosinus+Sqrt(rapport
535 greg 2.1 - sinus*sinus));
536 greg 2.3 r_transmis.ppar2 = r_incident.ppar2*2.*Sqrt(rapport)*cosinus/
537     (Sqrt(rapport-sinus*sinus)+rapport*cosinus);
538     r_transmis.pper2 = r_incident.pper2*2.*cosinus/(cosinus+Sqrt(rapport
539 greg 2.1 - sinus*sinus));
540 greg 2.3 r_transmis.e = (r_incident.e/2)*(Sqrt(rapport-sinus*sinus)/cosinus)
541 greg 2.1 *(((r_transmis.ppar1*r_transmis.ppar1+r_transmis.pper1*
542     r_transmis.pper1)
543     /(r_incident.ppar1*r_incident.ppar1+r_incident.pper1*
544     r_incident.pper1))+
545     ((r_transmis.ppar2*r_transmis.ppar2+r_transmis.pper2*r_transmis.pper2)
546     /(r_incident.ppar2*r_incident.ppar2+r_incident.pper2*r_incident.pper2)));
547     if(r_incident.n == 1.) r_transmis.n = prism.np;
548     else r_transmis.n = 1.;
549     }
550     else r_transmis.e = 0.;
551    
552     /* a la sortie de cette routine r_transmis.orig et .dest ne sont pas definis!*/
553    
554     return r_transmis;
555     }
556    
557    
558    
559    
560     #define ensuite(rayon,prob_passage,destination) r_suite = rayon; \
561     r_suite.e = prob_passage(rayon)*rayon.e; \
562     r_suite.dest = destination; \
563     if ( r_suite.e > seuil ) trace_rayon(r_suite)
564    
565    
566 schorsch 2.7 static void
567     trace_rayon(
568     TRAYON r_incident
569     )
570 greg 2.1 {
571     /* trace le rayon donne */
572     TRAYON r_reflechi,r_transmis,r_suite;
573    
574     switch (r_incident.dest)
575     {
576     case ALPHA:
577     if ( r_incident.orig == ALPHA )
578     {
579     r_reflechi = reflexion(r_incident);
580     sortie(r_reflechi);
581    
582     r_transmis = transmission(r_incident);
583     r_transmis.orig = ALPHA;
584    
585     ensuite(r_transmis,prob_alpha_beta,BETA);
586     ensuite(r_transmis,prob_alpha_gamma,GAMMA);
587     }
588     else
589     {
590     r_reflechi = reflexion(r_incident);
591     r_reflechi.orig = ALPHA;
592     ensuite(r_reflechi,prob_alpha_beta,BETA);
593     ensuite(r_reflechi,prob_alpha_gamma,GAMMA);
594    
595     r_transmis = transmission(r_incident);
596     sortie(r_transmis);
597     }
598     break;
599     case BETA:
600     r_reflechi = transbetaalpha(reflexion(transalphabeta(r_incident)));
601     r_reflechi.orig = BETA;
602     r_transmis = transbetaalpha(transmission(transalphabeta
603     (r_incident)));
604     r_transmis.orig = GAMMA;
605     if ( r_incident.n > 1.0 ) /* le rayon vient de l'interieur */
606     {
607     ensuite(r_reflechi,prob_beta_alpha,ALPHA);
608     ensuite(r_reflechi,prob_beta_gamma,GAMMA);
609    
610     ensuite(r_transmis,prob_beta_gamma,GAMMA);
611     ensuite(r_transmis,prob_beta_delta,DELTA);
612     }
613     else /* le rayon vient de l'exterieur */
614     {
615     ensuite(r_reflechi,prob_beta_gamma,GAMMA);
616     ensuite(r_reflechi,prob_beta_delta,DELTA);
617    
618     ensuite(r_transmis,prob_beta_alpha,ALPHA);
619     ensuite(r_transmis,prob_beta_gamma,GAMMA);
620     }
621     break;
622     case GAMMA:
623     r_reflechi = transgammaalpha(reflexion(transalphagamma(r_incident)));
624     r_reflechi.orig = GAMMA;
625     r_transmis = transgammaalpha(transmission(transalphagamma
626     (r_incident)));
627     r_transmis.orig = GAMMA;
628     if ( r_incident.n > 1.0 ) /* le rayon vient de l'interieur */
629     {
630     ensuite(r_reflechi,prob_gamma_alpha,ALPHA);
631     ensuite(r_reflechi,prob_gamma_beta,BETA);
632    
633     ensuite(r_transmis,prob_gamma_beta,BETA);
634     ensuite(r_transmis,prob_gamma_delta,DELTA);
635     }
636     else /* le rayon vient de l'exterieur */
637     {
638     ensuite(r_reflechi,prob_gamma_beta,BETA);
639     ensuite(r_reflechi,prob_gamma_delta,DELTA);
640    
641     ensuite(r_transmis,prob_gamma_alpha,ALPHA);
642     ensuite(r_transmis,prob_gamma_beta,BETA);
643     }
644     break;
645     case DELTA:
646     if ( r_incident.orig != DELTA ) sortie(r_incident);
647     else
648     {
649     ensuite(r_incident,prob_delta_beta,BETA);
650     ensuite(r_incident,prob_delta_gamma,GAMMA);
651     }
652     break;
653     }
654     return;
655     }
656    
657     #undef ensuite
658    
659 schorsch 2.7 static void
660     inverser(
661     TRAYON *r1,
662     TRAYON *r2
663     )
664 greg 2.1 {
665 schorsch 2.7 TRAYON temp;
666     temp = *r1;
667     *r1 = *r2;
668     *r2 = temp;
669 greg 2.1 }
670    
671    
672    
673 schorsch 2.7 static void
674     setprism(void)
675 greg 2.1 {
676     double d;
677     TRAYON r_initial,rsource;
678 schorsch 2.7 int i,j;
679 greg 2.1
680     prismclock = eclock;
681     r_initial.ppar1 = r_initial.pper2 = 1.;
682     r_initial.pper1 = r_initial.ppar2 = 0.;
683    
684     d = 1; prism.a = funvalue("arg", 1, &d);
685     if(prism.a < 0.) goto badopt;
686     d = 2; prism.b = funvalue("arg", 1, &d);
687     if(prism.b < 0.) goto badopt;
688     d = 3; prism.c = funvalue("arg", 1, &d);
689     if(prism.c < 0.) goto badopt;
690     d = 4; prism.d = funvalue("arg", 1, &d);
691     if(prism.d < 0.) goto badopt;
692     d = 5; prism.np = funvalue("arg", 1, &d);
693     if(prism.np < 1.) goto badopt;
694     d = 6; seuil = funvalue("arg", 1, &d);
695     if (seuil < 0. || seuil >=1) goto badopt;
696     d = 7; tot_ref = (int)(funvalue("arg", 1, &d) + .5);
697     if (tot_ref != 1 && tot_ref != 2 && tot_ref != 4) goto badopt;
698     if (tot_ref < 4 )
699     {
700     d = 8; fact_ref[tot_ref] = funvalue("arg", 1, &d);
701     if (fact_ref[tot_ref] < 0. || fact_ref[tot_ref] > 1.) goto badopt;
702     }
703     d = 9; tolerance = funvalue("arg", 1, &d);
704     if (tolerance <= 0.) goto badopt;
705     d = 10; tolsource = funvalue("arg", 1, &d);
706     if (tolsource < 0. ) goto badopt;
707     X(r_initial) = varvalue("Dx");
708     Y(r_initial) = varvalue("Dy");
709     Z(r_initial) = varvalue("Dz");
710     #ifdef DEBUG
711     fprintf(stderr,"dx=%lf dy=%lf dz=%lf\n",X(r_initial),Y(r_initial),Z(r_initial));
712     #endif
713    
714     /* initialisation */
715     prepare_matrices();
716     r_initial.e = 1.0;
717     r_initial.n = 1.0;
718    
719     if(ray!=NULL) free(ray);
720     nbrayons = 0;
721     /* determination de l'origine et de la destination du rayon initial */
722    
723     if ( X(r_initial) != 0.)
724     {
725     if ( X(r_initial) > 0. )
726     {
727     r_initial.orig = r_initial.dest = ALPHA;
728     sens = 1;
729     }
730     else if ( X(r_initial) < 0. )
731     {
732     r_initial.orig = r_initial.dest = DELTA;
733     sens = -1;
734     }
735    
736     normalize(r_initial.v);
737    
738     trace_rayon(r_initial);
739    
740     X(rsource) = varvalue("DxA");
741     Y(rsource) = varvalue("DyA");
742     Z(rsource) = varvalue("DzA");
743     nosource = ( X(rsource)==0. && Y(rsource)==0. && Z(rsource)==0. );
744     if ( !nosource )
745     {
746     for (j=0; j<nbrayons; j++)
747     {
748     if ( !compare(ray[j],rsource,tolsource) ) ray[j].e =0.;
749     }
750     }
751     for (j = 0; j < nbrayons; j++)
752     {
753     for (i = j+1; i < nbrayons; i++)
754     {
755     if (ray[j].e < ray[i].e) inverser(&ray[j],&ray[i]);
756     }
757     }
758    
759     bidon = 1;
760     }
761     else bidon = 0;
762     return;
763    
764     /* message puis sortie si erreur dans la ligne de commande */
765     badopt:
766     bidon = BADVAL;
767     return;
768     }
769 greg 2.2
770    
771     static double
772 schorsch 2.7 l_get_val(
773     char *nm
774     )
775 greg 2.2 {
776     int val, dir, i, trouve, curseur;
777     int nb;
778     double valeur;
779 greg 2.9 TRAYON *rayt=NULL, raynull;
780 greg 2.2
781     if (prismclock < 0 || prismclock < eclock) setprism();
782     if (bidon == BADVAL) {
783     errno = EDOM;
784     return(0.0);
785     }
786     val = (int)(argument(1) + .5);
787     dir = (int)(argument(2) + .5);
788     nb = (int)(argument(3) + .5);
789     X(raynull) = bidon;
790     Y(raynull) = Z(raynull) = 0.;
791     raynull.e = 0.;
792     trouve = curseur = 0;
793     if ( !nosource && nb==2 ) nb=1; /* on est en train de tracer la source
794     a partir de sa seconde source virtuelle */
795     #ifdef DEBUG
796     fprintf(stderr, " On considere le rayon no: %d\n", nb);
797     #endif
798     for(i=0; i < nbrayons &&!trouve; i++)
799     {
800     if(ray[i].v[0] * dir * sens >= 0.) curseur ++;
801     if(curseur == nb)
802     {
803     rayt = &ray[i];
804     trouve = 1;
805     }
806     }
807     if(!trouve) rayt = &raynull;
808     switch(val) {
809     case 0 : valeur = rayt->v[0];
810     break;
811     case 1 : valeur = rayt->v[1];
812     break;
813     case 2 : valeur = rayt->v[2];
814     break;
815     case 3 : valeur = rayt->e;
816     break;
817     default : errno = EDOM; return(0.0);
818     }
819     #ifdef DEBUG
820     fprintf(stderr, "get_val( %i, %i, %i) = %lf\n",val,dir,nb,valeur);
821     #endif
822     return valeur;
823     }
824    
825 greg 2.1
826 greg 2.10 void
827 schorsch 2.7 setprismfuncs(void) /* declared in func.h */
828 greg 2.1 {
829     funset("fprism_val", 3, '=', l_get_val);
830     }