ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/rt/fprism.c
(Generate patch)

Comparing ray/src/rt/fprism.c (file contents):
Revision 2.1 by greg, Wed Sep 29 10:40:31 1993 UTC vs.
Revision 2.8 by greg, Fri Feb 18 00:40:25 2011 UTC

# Line 1 | Line 1
1 – /* Copyright (c) 1993 Regents of the University of California */
2 –
1   #ifndef lint
2 < static char SCCSid[] = "$SunId$ LBL";
2 > static const char       RCSid[] = "$Id$";
3   #endif
6 –
4   /* Ce programme calcule les directions et les energies des rayons lumineux
5     resultant du passage d'un rayon au travers d'un vitrage prismatique
6  
# Line 11 | Line 8 | static char SCCSid[] = "$SunId$ LBL";
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(x)
22 < double x;
21 > Sqrt(
22 >        double x
23 > )
24   {
25          if (x < 0.)
26                  return(0.);
27          return(sqrt(x));
28   }
27 – #define sqrt(x) Sqrt(x)
29  
30   /* definitions de macros utiles */
31  
# Line 83 | Line 84 | static int tot_ref;                    /* flag pour les surfaces
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  */
86 – static double Nx;
87   static int bidon;
88   #define BADVAL  (-10)
89   static long prismclock = -1;
# Line 93 | Line 93 | 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 – static TRAYON rtemp;                    /* variable temporaire             */
96  
98 – extern double argument();
99 – extern double varvalue();
100 – extern double funvalue();
101 – extern long eclock;
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
124 < prepare_matrices()
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 */
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;
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
140 < tfm(mat,v_old,v_new)
141 < MAT4 mat;
142 < FVECT v_old,v_new;
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;
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;
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
# Line 141 | Line 159 | FVECT v_old,v_new;
159  
160  
161   static double
162 < prob_alpha_beta(r)
163 < TRAYON r;
162 > prob_alpha_beta(
163 >        TRAYON r
164 > )
165   {
166 < /* calcul de la probabilite de passage de alpha a beta */
167 < double prob,test;
166 >        /* calcul de la probabilite de passage de alpha a beta */
167 >        double prob,test;
168  
169 < if ( X(r) != 0. )
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.;
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;
176 >        else prob = 0.;
177 >        return prob;
178   }
179  
180  
181   static double
182 < prob_beta_alpha(r)
183 < TRAYON r;
182 > prob_beta_alpha(
183 >        TRAYON r
184 > )
185   {
186 < /* calcul de la probabilite de passage de beta a aplha */
187 < double prob,test;
186 >        /* calcul de la probabilite de passage de beta a aplha */
187 >        double prob,test;
188  
189 < if ( X(r) != 0. )
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.;
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;
196 >        else prob = 0.;
197 >        return prob;
198   }
199  
200  
201 < double prob_gamma_alpha(r)
202 < TRAYON r;
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;
206 >        /* calcul de la probabilite de passage de gamma a alpha */
207 >        double prob,test;
208  
209 < if ( X(r) != 0. )
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);
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;
216 >        else prob = 0.;
217 >        return prob;
218   }
219  
220   #undef A
# Line 201 | Line 223 | TRAYON r;
223   #undef D
224  
225  
226 < static
227 < v_par(v,v_out)
228 < FVECT v,v_out;
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;
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;
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
247 < v_per(v,v_out)
248 < FVECT v,v_out;
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;
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;
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(r_initial)
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 */
243 – TRAYON r_initial;
271   {
272 < TRAYON r_final;
273 < FVECT vpar_temp1,vpar_temp2,vper_temp1,vper_temp2;
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;
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(r_initial)
300 < /* transforme le rayon r_initial de la base associee a beta dans
301 <   la base associee a alpha */
275 < TRAYON r_initial;
299 > transbetaalpha(
300 >        TRAYON r_initial
301 > )
302   {
303 < TRAYON r_final;
304 < FVECT vpar_temp1,vpar_temp2,vper_temp1,vper_temp2;
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));
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;
327 >        }
328 >        return r_final;
329   }
330  
331  
332   static TRAYON
333 < transalphagamma(r_initial)
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 */
308 – TRAYON r_initial;
338   {
339 < TRAYON r_final;
340 < FVECT vpar_temp1,vpar_temp2,vper_temp1,vper_temp2;
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));
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;
361 >        }
362 >        return r_final;
363   }
364  
365  
366   static TRAYON
367 < transgammaalpha(r_initial)
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 */
341 – TRAYON r_initial;
372   {
373 < TRAYON r_final;
374 < FVECT vpar_temp1,vpar_temp2,vper_temp1,vper_temp2;
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;
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 <
401 < static
402 < sortie(r)
403 < TRAYON r;
400 > static int
401 > compare(
402 >        TRAYON r1,
403 >        TRAYON r2,
404 >        double marge
405 > )
406   {
407 < int i = 0;
376 < int egalite = 0;
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 < if(r.e > seuil)
415 >
416 >
417 >
418 > static void
419 > sortie(
420 >        TRAYON r
421 > )
422   {
423 < while (i < nbrayons && egalite == 0)
424 <  {
425 <   raytemp = &ray[i];
426 <   egalite = compare(r,*raytemp,tolerance);
427 <   if (egalite) raytemp->e = raytemp->e + r.e;
428 <   else i = i + 1;
429 <  }
430 < if (egalite == 0)
431 <  {
432 <   if (nbrayons == 0) ray = (TRAYON *)calloc(1,sizeof(TRAYON));
433 <   else ray = (TRAYON *)realloc(ray, (nbrayons+1)*(sizeof(TRAYON)));        
434 <   if (ray == NULL)
435 <     error(SYSTEM, "out of memory in sortie\n");
436 <   raytemp = &ray[nbrayons];
437 <   raytemp->v[0] = X(r);
438 <   raytemp->v[1] = Y(r);
439 <   raytemp->v[2] = Z(r);
440 <   raytemp->e = r.e;
441 <   nbrayons++;
442 <  }
443 < }
444 < return;
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
455 < trigo(r)
456 < TRAYON r;
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;
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(r_incident)
474 < TRAYON r_incident;
473 > reflexion(
474 >        TRAYON r_incident
475 > )
476   {
477   /* calcul du rayon reflechi par une face */
478   TRAYON r_reflechi;
# Line 432 | Line 482 | TRAYON 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)
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;
# Line 442 | Line 492 | TRAYON r_incident;
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)));
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
# Line 463 | Line 513 | TRAYON r_incident;
513  
514  
515   static TRAYON
516 < transmission(r_incident)
517 < TRAYON r_incident;
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)
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)*
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
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
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)
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*
# Line 506 | Line 557 | TRAYON r_incident;
557  
558  
559  
509 – static int
510 – compare(r1,r2,marge)
511 – TRAYON r1, r2;
512 – double marge;
513 –
514 – {
515 – double arctg1, arctg2;
516 –
517 – arctg1 = atan2(Y(r1),X(r1));
518 – arctg2 = atan2(Y(r2),X(r2));
519 – if ((arctg1 - marge <= arctg2) && (arctg1 + marge >= arctg2)) return 1;
520 – else return 0;
521 – }
522 –
523 –
524 –
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
567 < trace_rayon(r_incident)
568 < TRAYON r_incident;
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;
# Line 620 | Line 656 | TRAYON r_incident;
656  
657   #undef ensuite
658  
659 < static
660 < inverser(r1,r2)
661 < TRAYON *r1,*r2;
662 <
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;
665 >        TRAYON temp;
666 >        temp = *r1;
667 >        *r1 = *r2;
668 >        *r2 = temp;
669   }
670  
671  
672  
673 < static double
674 < l_get_val()
638 <
673 > static void
674 > setprism(void)
675   {
640 – int val, dir, i, trouve, curseur;
641 – int nb;
642 – double valeur;
643 – TRAYON *rayt, raynull;
644 –
645 – if (prismclock < 0 || prismclock < eclock) setprism();
646 – if (bidon == BADVAL) {
647 –        errno = EDOM;
648 –        return(0.0);
649 – }
650 – val = (int)(argument(1) + .5);
651 – dir = (int)(argument(2) + .5);
652 – nb = (int)(argument(3) + .5);
653 – X(raynull) = bidon;
654 – Y(raynull) = Z(raynull) = 0.;
655 – raynull.e = 0.;
656 – trouve = curseur = 0;
657 – if ( !nosource && nb==2 ) nb=1; /* on est en train de tracer la source
658 –                                     a partir de sa seconde source virtuelle */
659 – #ifdef DEBUG
660 – fprintf(stderr, " On considere le rayon no: %d\n", nb);
661 – #endif
662 – for(i=0; i < nbrayons &&!trouve; i++)
663 –  {
664 –   if(ray[i].v[0] * dir * sens >= 0.) curseur ++;
665 –   if(curseur == nb)
666 –   {
667 –    rayt = &ray[i];
668 –    trouve = 1;
669 –   }
670 –  }
671 – if(!trouve) rayt = &raynull;
672 – switch(val) {
673 –        case 0 : valeur = rayt->v[0];
674 –                 break;
675 –        case 1 : valeur = rayt->v[1];
676 –                 break;
677 –        case 2 : valeur = rayt->v[2];
678 –                 break;
679 –        case 3 : valeur = rayt->e;
680 –                 break;
681 –        default : errno = EDOM; return(0.0);
682 –    }
683 – #ifdef DEBUG
684 –  fprintf(stderr, "get_val( %i, %i, %i) = %lf\n",val,dir,nb,valeur);
685 – #endif
686 –  return valeur;
687 – }
688 –
689 –
690 – static
691 – setprism()
692 –
693 – {
676   double d;
677   TRAYON r_initial,rsource;
678 < int i,j,k;
678 > int i,j;
679  
680   prismclock = eclock;
681   r_initial.ppar1 = r_initial.pper2 = 1.;
# Line 785 | Line 767 | if ( X(r_initial) != 0.)
767   return;
768   }
769  
770 < setprismfuncs()
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, 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 > extern void
827 > setprismfuncs(void)  /* declared in func.h */
828   {
829   funset("fprism_val", 3, '=', l_get_val);
830   }

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines