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