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, 8 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

# Content
1 #ifndef lint
2 static const char RCSid[] = "$Id: fprism.c,v 2.9 2013/08/07 05:10:09 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 "ray.h"
12 #include "calcomp.h"
13 #include "func.h"
14
15 #ifdef NOSTRUCTASSIGN
16 static double err = "No structure assignment!"; /* generate compiler error */
17 #endif
18
19
20 static double
21 Sqrt(
22 double x
23 )
24 {
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
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
119
120 /* Definition des routines */
121
122 #define term(a,b) a/Sqrt(a*a+b*b)
123 static void
124 prepare_matrices(void)
125 {
126 /* preparation des matrices de changement de bases */
127
128 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 }
136 #undef term
137
138
139 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 }
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 prob_alpha_beta(
163 TRAYON r
164 )
165 {
166 /* calcul de la probabilite de passage de alpha a beta */
167 double prob,test;
168
169 if ( X(r) != 0. )
170 {
171 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 }
176 else prob = 0.;
177 return prob;
178 }
179
180
181 static double
182 prob_beta_alpha(
183 TRAYON r
184 )
185 {
186 /* calcul de la probabilite de passage de beta a aplha */
187 double prob,test;
188
189 if ( X(r) != 0. )
190 {
191 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 }
196 else prob = 0.;
197 return prob;
198 }
199
200
201 static double
202 prob_gamma_alpha(
203 TRAYON r
204 )
205 {
206 /* calcul de la probabilite de passage de gamma a alpha */
207 double prob,test;
208
209 if ( X(r) != 0. )
210 {
211 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 }
216 else prob = 0.;
217 return prob;
218 }
219
220 #undef A
221 #undef B
222 #undef C
223 #undef D
224
225
226 static void
227 v_par(
228 FVECT v,
229 FVECT v_out
230 )
231 /* calcule le vecteur par au plan d'incidence lie a v */
232 {
233 FVECT v_temp;
234 double det;
235
236 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 }
244
245
246 static void
247 v_per(
248 FVECT v,
249 FVECT v_out
250 )
251 /* calcule le vecteur perp au plan d'incidence lie a v */
252 {
253 FVECT v_temp;
254 double det;
255
256 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 }
263
264
265 static TRAYON
266 transalphabeta(
267 TRAYON r_initial
268 )
269 /* transforme le rayon r_initial de la base associee a alpha dans
270 la base associee a beta */
271 {
272 TRAYON r_final;
273 FVECT vpar_temp1,vpar_temp2,vper_temp1,vper_temp2;
274
275 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 }
296
297
298 static TRAYON
299 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
327 }
328 return r_final;
329 }
330
331
332 static TRAYON
333 transalphagamma(
334 TRAYON r_initial
335 )
336 /* transforme le rayon r_initial de la base associee a alpha dans
337 la base associee a gamma */
338 {
339 TRAYON r_final;
340 FVECT vpar_temp1,vpar_temp2,vper_temp1,vper_temp2;
341
342 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
361 }
362 return r_final;
363 }
364
365
366 static TRAYON
367 transgammaalpha(
368 TRAYON r_initial
369 )
370 /* transforme le rayon r_initial de la base associee a gamma dans
371 la base associee a alpha */
372 {
373 TRAYON r_final;
374 FVECT vpar_temp1,vpar_temp2,vper_temp1,vper_temp2;
375
376 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 }
397
398
399
400 static int
401 compare(
402 TRAYON r1,
403 TRAYON r2,
404 double marge
405 )
406 {
407 double arctg1, arctg2;
408
409 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 }
414
415
416
417
418 static void
419 sortie(
420 TRAYON r
421 )
422 {
423 int i = 0;
424 int egalite = 0;
425
426
427 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 }
452
453
454 static void
455 trigo(
456 TRAYON r
457 )
458 /* calcule les grandeurs trigonometriques relatives au rayon incident
459 et le rapport entre les indices du milieu refracteur et incident */
460 {
461 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 }
470
471
472 static TRAYON
473 reflexion(
474 TRAYON r_incident
475 )
476 {
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 if(sinus > Sqrt(rapport) || r_incident.dest == tot_ref)
486 {
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 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 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 transmission(
517 TRAYON r_incident
518 )
519 {
520 /* calcul du rayon refracte par une face */
521 TRAYON r_transmis;
522
523 r_transmis = r_incident;
524 trigo(r_incident);
525 if (sinus <= Sqrt(rapport) && r_incident.dest != tot_ref)
526 {
527 X(r_transmis) = (X(r_incident)/(fabs(X(r_incident))))*
528 (Sqrt(1.-(Y(r_incident)*Y(r_incident)+Z(r_incident)*
529 Z(r_incident))/rapport));
530 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 - sinus*sinus));
536 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 - sinus*sinus));
540 r_transmis.e = (r_incident.e/2)*(Sqrt(rapport-sinus*sinus)/cosinus)
541 *(((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 static void
567 trace_rayon(
568 TRAYON r_incident
569 )
570 {
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 static void
660 inverser(
661 TRAYON *r1,
662 TRAYON *r2
663 )
664 {
665 TRAYON temp;
666 temp = *r1;
667 *r1 = *r2;
668 *r2 = temp;
669 }
670
671
672
673 static void
674 setprism(void)
675 {
676 double d;
677 TRAYON r_initial,rsource;
678 int i,j;
679
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
770
771 static double
772 l_get_val(
773 char *nm
774 )
775 {
776 int val, dir, i, trouve, curseur;
777 int nb;
778 double valeur;
779 TRAYON *rayt=NULL, raynull;
780
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
826 void
827 setprismfuncs(void) /* declared in func.h */
828 {
829 funset("fprism_val", 3, '=', l_get_val);
830 }