ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/rt/fprism.c
Revision: 2.7
Committed: Tue Mar 30 16:13:01 2004 UTC (20 years ago) by schorsch
Content type: text/plain
Branch: MAIN
CVS Tags: rad3R7P2, rad3R7P1, rad4R0, rad3R6, rad3R6P1, rad3R8, rad3R9
Changes since 2.6: +299 -257 lines
Log Message:
Continued ANSIfication. There are only bits and pieces left now.

File Contents

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