ViewVC Help

View File | Revision Log | Show Annotations | Download File | Root Listing

root/Development/ray/src/rt/fprism.c

Comparing ray/src/rt/fprism.c (file contents):
Revision 2.4 by greg, Sat Feb 22 02:07:28 2003 UTC vs.
Revision 2.7 by schorsch, Tue Mar 30 16:13:01 2004 UTC

#	Line 8 | Line 8 | static const char      RCSid[] = "$Id$";
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(x)
21	<	double x;
20	>	Sqrt(
21	>	double x
22	>	)
23		{
24		if (x < 0.)
25		return(0.);
#	Line 79 | Line 83 | static int tot_ref;                    /* flag pour les surfaces
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  */
82	–	static double Nx;
86		static int bidon;
87		#define BADVAL  (-10)
88		static long prismclock = -1;
#	Line 89 | Line 92 | 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             */
92	–	static TRAYON rtemp;                    /* variable temporaire             */
95
94	–	extern double argument();
95	–	extern double varvalue();
96	–	extern double funvalue();
97	–	extern long eclock;
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
123	<	prepare_matrices()
122	>	static void
123	>	prepare_matrices(void)
124		{
125	<	/* preparation des matrices de changement de bases */
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;
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
139	<	tfm(mat,v_old,v_new)
140	<	MAT4 mat;
141	<	FVECT v_old,v_new;
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;
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;
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
#	Line 137 | Line 158 | FVECT v_old,v_new;
158
159
160		static double
161	<	prob_alpha_beta(r)
162	<	TRAYON r;
161	>	prob_alpha_beta(
162	>	TRAYON r
163	>	)
164		{
165	<	/* calcul de la probabilite de passage de alpha a beta */
166	<	double prob,test;
165	>	/* calcul de la probabilite de passage de alpha a beta */
166	>	double prob,test;
167
168	<	if ( X(r) != 0. )
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.;
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;
175	>	else prob = 0.;
176	>	return prob;
177		}
178
179
180		static double
181	<	prob_beta_alpha(r)
182	<	TRAYON r;
181	>	prob_beta_alpha(
182	>	TRAYON r
183	>	)
184		{
185	<	/* calcul de la probabilite de passage de beta a aplha */
186	<	double prob,test;
185	>	/* calcul de la probabilite de passage de beta a aplha */
186	>	double prob,test;
187
188	<	if ( X(r) != 0. )
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.;
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;
195	>	else prob = 0.;
196	>	return prob;
197		}
198
199
200	<	double prob_gamma_alpha(r)
201	<	TRAYON r;
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;
205	>	/* calcul de la probabilite de passage de gamma a alpha */
206	>	double prob,test;
207
208	<	if ( X(r) != 0. )
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);
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;
215	>	else prob = 0.;
216	>	return prob;
217		}
218
219		#undef A
#	Line 197 | Line 222 | TRAYON r;
222		#undef D
223
224
225	<	static
226	<	v_par(v,v_out)
227	<	FVECT v,v_out;
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;
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;
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
246	<	v_per(v,v_out)
247	<	FVECT v,v_out;
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;
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;
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(r_initial)
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 */
239	–	TRAYON r_initial;
270		{
271	<	TRAYON r_final;
272	<	FVECT vpar_temp1,vpar_temp2,vper_temp1,vper_temp2;
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;
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(r_initial)
299	<	/* transforme le rayon r_initial de la base associee a beta dans
300	<	la base associee a alpha */
271	<	TRAYON r_initial;
298	>	transbetaalpha(
299	>	TRAYON r_initial
300	>	)
301		{
302	<	TRAYON r_final;
303	<	FVECT vpar_temp1,vpar_temp2,vper_temp1,vper_temp2;
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));
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;
326	>	}
327	>	return r_final;
328		}
329
330
331		static TRAYON
332	<	transalphagamma(r_initial)
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 */
304	–	TRAYON r_initial;
337		{
338	<	TRAYON r_final;
339	<	FVECT vpar_temp1,vpar_temp2,vper_temp1,vper_temp2;
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));
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;
360	>	}
361	>	return r_final;
362		}
363
364
365		static TRAYON
366	<	transgammaalpha(r_initial)
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 */
337	–	TRAYON r_initial;
371		{
372	<	TRAYON r_final;
373	<	FVECT vpar_temp1,vpar_temp2,vper_temp1,vper_temp2;
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;
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(r1,r2,marge)
401	<	TRAYON r1, r2;
402	<	double marge;
403	<
400	>	compare(
401	>	TRAYON r1,
402	>	TRAYON r2,
403	>	double marge
404	>	)
405		{
406	<	double arctg1, arctg2;
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;
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
418	<	sortie(r)
419	<	TRAYON r;
417	>	static void
418	>	sortie(
419	>	TRAYON r
420	>	)
421		{
422	<	int i = 0;
423	<	int egalite = 0;
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(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;
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
454	<	trigo(r)
455	<	TRAYON r;
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;
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(r_incident)
473	<	TRAYON r_incident;
472	>	reflexion(
473	>	TRAYON r_incident
474	>	)
475		{
476		/* calcul du rayon reflechi par une face */
477		TRAYON r_reflechi;
#	Line 475 | Line 512 | TRAYON r_incident;
512
513
514		static TRAYON
515	<	transmission(r_incident)
516	<	TRAYON r_incident;
515	>	transmission(
516	>	TRAYON r_incident
517	>	)
518		{
519		/* calcul du rayon refracte par une face */
520		TRAYON r_transmis;
#	Line 524 | Line 562 | TRAYON r_incident;
562		if ( r_suite.e > seuil ) trace_rayon(r_suite)
563
564
565	<	static
566	<	trace_rayon(r_incident)
567	<	TRAYON r_incident;
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;
#	Line 616 | Line 655 | TRAYON r_incident;
655
656		#undef ensuite
657
658	<	static
659	<	inverser(r1,r2)
660	<	TRAYON *r1,*r2;
661	<
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;
664	>	TRAYON temp;
665	>	temp = *r1;
666	>	*r1 = *r2;
667	>	*r2 = temp;
668		}
669
670
671
672	<	static
673	<	setprism()
672	>	static void
673	>	setprism(void)
674		{
675		double d;
676		TRAYON r_initial,rsource;
677	<	int i,j,k;
677	>	int i,j;
678
679		prismclock = eclock;
680		r_initial.ppar1 = r_initial.pper2 = 1.;
#	Line 728 | Line 768 | if ( X(r_initial) != 0.)
768
769
770		static double
771	<	l_get_val()
772	<
771	>	l_get_val(
772	>	char *nm
773	>	)
774		{
775		int val, dir, i, trouve, curseur;
776		int nb;
#	Line 781 | Line 822 | l_get_val()
822		}
823
824
825	<	setprismfuncs()
825	>	extern void
826	>	setprismfuncs(void)  /* declared in func.h */
827		{
828		funset("fprism_val", 3, '=', l_get_val);
829		}

Diff Legend

–	Removed lines
+	Added lines
<	Changed lines (old)
>	Changed lines (new)