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 |
|
|
} |