Back to index

radiance  4R0+20100331
fprism.c
Go to the documentation of this file.
00001 #ifndef lint
00002 static const char    RCSid[] = "$Id: fprism.c,v 2.7 2004/03/30 16:13:01 schorsch Exp $";
00003 #endif
00004 /* Ce programme calcule les directions et les energies des rayons lumineux
00005    resultant du passage d'un rayon au travers d'un vitrage prismatique
00006  
00007    1991, LESO - EPFL, R. Compagnon - F. Di Pasquale                        */
00008 
00009 #include "standard.h"
00010 
00011 #include "calcomp.h"
00012 #include "func.h"
00013 
00014 #ifdef NOSTRUCTASSIGN
00015 static double err = "No structure assignment!";  /* generate compiler error */
00016 #endif
00017 
00018 
00019 static double
00020 Sqrt(
00021        double x
00022 )
00023 {
00024        if (x < 0.)
00025               return(0.);
00026        return(sqrt(x));
00027 }
00028 
00029 /* definitions de macros utiles */
00030 
00031 #define ALPHA 0
00032 #define BETA 1
00033 #define GAMMA 2
00034 #define DELTA 3
00035 #define AUCUNE 4
00036 #define X(r) r.v[0]
00037 #define Y(r) r.v[1]
00038 #define Z(r) r.v[2]
00039 #define XX(v) v[0]
00040 #define YY(v) v[1]
00041 #define ZZ(v) v[2]
00042 #define alpha_beta(v_alpha,v_beta) tfm(matbt,v_alpha,v_beta)
00043 #define beta_alpha(v_beta,v_alpha) tfm(matb,v_beta,v_alpha)
00044 #define alpha_gamma(v_alpha,v_gamma) tfm(matct,v_alpha,v_gamma)
00045 #define gamma_alpha(v_gamma,v_alpha) tfm(matc,v_gamma,v_alpha)
00046 #define prob_alpha_gamma(r) (1.-prob_alpha_beta(r))
00047 #define prob_beta_gamma(r) (1.-prob_beta_alpha(r))
00048 #define prob_gamma_beta(r) (1.-prob_gamma_alpha(r))
00049 #define prob_delta_gamma(r) (1.-prob_delta_beta(r))
00050 #define prob_beta_delta(r) (prob_beta_alpha(r))
00051 #define prob_gamma_delta(r) (prob_gamma_alpha(r))
00052 #define prob_delta_beta(r) (prob_alpha_beta(r))
00053 
00054 
00055 /* Definitions des types de donnees */
00056 
00057 typedef struct { FVECT v;          /* direction */
00058                double ppar1,pper1,
00059                      ppar2,pper2;  /* polarisations */
00060                double e;           /* energie */
00061                double n;           /* milieu dans lequel on se propage */
00062                int orig,dest;             /* origine et destination          */
00063               } TRAYON;     
00064 
00065 typedef struct { double a,b,c,d;   /* longueurs caracteristiques      */
00066                double np;          /* indice de refraction            */
00067               } TPRISM;
00068 
00069 
00070 /* Definitions des variables globales */
00071 
00072 static TPRISM prism;                   /* le prisme ! */
00073 static MAT4 matb = MAT4IDENT;             /* matrices de changement de bases */
00074 static MAT4 matbt = MAT4IDENT;
00075 static MAT4 matc = MAT4IDENT;
00076 static MAT4 matct = MAT4IDENT;                   
00077 static double seuil;               /* seuil pour l'arret du trace     */
00078 static double sinus,cosinus;              /* sin et cos                  */
00079 static double rapport;                    /* rapport entre les indices des 
00080                                       milieux refracteur et incident  */
00081 static int tot_ref;                /* flag pour les surfaces 
00082                                       reflechissantes             */
00083 static double fact_ref[4]={1.0,1.0,1.0,1.0};     /* facteurs de reflexion     */ 
00084 static double tolerance;           /* degre de tol. pour les amalgames */
00085 static double tolsource;           /* degre de tol. pour les sources  */
00086 static int bidon;
00087 #define BADVAL       (-10)
00088 static long prismclock = -1;
00089 static int nosource;               /* indique que l'on ne trace pas
00090                                       en direction d'une source */
00091 static int sens;                   /* indique le sens de prop. du ray.*/  
00092 static int nbrayons;               /* indice des rayons sortants      */
00093 static TRAYON *ray;                /* tableau des rayons sortants     */
00094 static TRAYON *raytemp;                   /* variable temporaire         */
00095 
00096 
00097 static void prepare_matrices(void);
00098 static void tfm(MAT4 mat, FVECT v_old, FVECT v_new);
00099 static double prob_alpha_beta(TRAYON r);
00100 static double prob_beta_alpha(TRAYON r);
00101 static double prob_gamma_alpha(TRAYON r);
00102 static void v_par(FVECT v, FVECT v_out);
00103 static void v_per(FVECT v, FVECT v_out);
00104 static TRAYON transalphabeta(TRAYON r_initial);
00105 static TRAYON transbetaalpha(TRAYON r_initial);
00106 static TRAYON transalphagamma(TRAYON r_initial);
00107 static TRAYON transgammaalpha(TRAYON r_initial);
00108 static int compare(TRAYON r1, TRAYON r2, double marge);
00109 static void sortie(TRAYON r);
00110 static void trigo(TRAYON r);
00111 static TRAYON reflexion(TRAYON r_incident);
00112 static TRAYON transmission(TRAYON r_incident);
00113 static void trace_rayon(TRAYON r_incident);
00114 static void inverser(TRAYON *r1, TRAYON *r2);
00115 static void setprism(void);
00116 static double l_get_val(char *nm);
00117 
00118 
00119 /* Definition des routines */
00120 
00121 #define term(a,b) a/Sqrt(a*a+b*b)
00122 static void
00123 prepare_matrices(void)
00124 {
00125        /* preparation des matrices de changement de bases */
00126 
00127        matb[0][0] = matbt[0][0] = matb[1][1] = matbt[1][1] = term(prism.a,prism.d);
00128        matb[1][0] = matbt[0][1] = term(-prism.d,prism.a);
00129        matb[0][1] = matbt[1][0] = term(prism.d,prism.a);
00130        matc[0][0] = matct[0][0] = matc[1][1] = matct[1][1] = term(prism.b,prism.d);
00131        matc[1][0] = matct[0][1] = term(prism.d,prism.b);
00132        matc[0][1] = matct[1][0] = term(-prism.d,prism.b);
00133        return;
00134 }
00135 #undef term
00136 
00137 
00138 static void
00139 tfm(
00140        MAT4 mat,
00141        FVECT v_old,
00142        FVECT v_new
00143 )
00144 {
00145        /* passage d'un repere old au repere new par la matrice mat */
00146        FVECT v_temp;
00147 
00148        multv3(v_temp,v_old,mat);
00149        normalize(v_temp);
00150        VCOPY(v_new,v_temp);
00151        return;
00152 }
00153 
00154 #define A prism.a
00155 #define B prism.b
00156 #define C prism.c
00157 #define D prism.d
00158 
00159 
00160 static double
00161 prob_alpha_beta(
00162        TRAYON r
00163 )
00164 {
00165        /* calcul de la probabilite de passage de alpha a beta */
00166        double prob,test;
00167 
00168        if ( X(r) != 0. )
00169        {
00170               test = Y(r)/X(r);
00171               if ( test > B/D ) prob = 1.;
00172               else if ( test >= -A/D ) prob = (A+test*D)/(A+B);
00173               else prob = 0.;
00174        } 
00175        else prob = 0.;
00176        return prob;
00177 }
00178 
00179 
00180 static double
00181 prob_beta_alpha(
00182        TRAYON r
00183 )
00184 {
00185        /* calcul de la probabilite de passage de beta a aplha */
00186        double prob,test;
00187 
00188        if ( X(r) != 0. )
00189        {
00190               test = Y(r)/X(r);
00191               if ( test > B/D ) prob = (A+B)/(A+test*D);
00192               else if ( test >= -A/D ) prob = 1.;
00193               else prob = 0.;
00194        }
00195        else prob = 0.;
00196        return prob;
00197 }
00198 
00199 
00200 static double
00201 prob_gamma_alpha(
00202        TRAYON r
00203 )
00204 {
00205        /* calcul de la probabilite de passage de gamma a alpha */
00206        double prob,test;
00207 
00208        if ( X(r) != 0. )
00209        {
00210               test = Y(r)/X(r);
00211               if ( test > B/D ) prob = 0.;
00212               else if ( test >= -A/D ) prob = 1.;
00213               else prob = (A+B)/(B-test*D);
00214        } 
00215        else prob = 0.;
00216        return prob;
00217 }
00218 
00219 #undef A
00220 #undef B
00221 #undef C
00222 #undef D
00223 
00224 
00225 static void
00226 v_par(
00227        FVECT v,
00228        FVECT v_out
00229 )
00230 /* calcule le vecteur par au plan d'incidence lie a v */
00231 {
00232        FVECT v_temp;
00233        double det;
00234 
00235        det = Sqrt( (YY(v)*YY(v)+ZZ(v)*ZZ(v))*(YY(v)*YY(v)+ZZ(v)*ZZ(v))+
00236                      (XX(v)*XX(v)*YY(v)*YY(v))+(XX(v)*XX(v)*ZZ(v)*ZZ(v)) );
00237        XX(v_temp) = (YY(v)*YY(v)+ZZ(v)*ZZ(v))/det;
00238        YY(v_temp) = -( XX(v)*YY(v) )/det;
00239        ZZ(v_temp) = -( XX(v)*ZZ(v) )/det;
00240        VCOPY(v_out,v_temp); 
00241        return;
00242 }
00243 
00244 
00245 static void
00246 v_per(
00247        FVECT v,
00248        FVECT v_out
00249 )
00250 /* calcule le vecteur perp au plan d'incidence lie a v */
00251 {
00252        FVECT v_temp;
00253        double det;
00254 
00255        det = Sqrt( (ZZ(v)*ZZ(v)+YY(v)*YY(v)) );
00256        XX(v_temp) = 0.;
00257        YY(v_temp) = -ZZ(v)/det;
00258        ZZ(v_temp) = YY(v)/det;
00259        VCOPY(v_out,v_temp);
00260        return;
00261 }
00262 
00263 
00264 static TRAYON
00265 transalphabeta(
00266        TRAYON r_initial
00267 )
00268 /* transforme le rayon r_initial de la base associee a alpha dans
00269    la base associee a beta */ 
00270 {
00271        TRAYON r_final;
00272        FVECT vpar_temp1,vpar_temp2,vper_temp1,vper_temp2;
00273 
00274        r_final = r_initial;
00275        alpha_beta(r_initial.v,r_final.v);
00276        if ((Y(r_initial) != 0. || Z(r_initial) != 0.)&&(Y(r_final) !=0. || Z(r_final)!= 0.))
00277        {
00278               v_par(r_initial.v,vpar_temp1);
00279               alpha_beta(vpar_temp1,vpar_temp1);
00280               v_per(r_initial.v,vper_temp1);
00281               alpha_beta(vper_temp1,vper_temp1);
00282               v_par(r_final.v,vpar_temp2);
00283               v_per(r_final.v,vper_temp2); 
00284               r_final.ppar1 = (r_initial.ppar1*fdot(vpar_temp1,vpar_temp2))+
00285                      (r_initial.pper1*fdot(vper_temp1,vpar_temp2));
00286               r_final.pper1 = (r_initial.ppar1*fdot(vpar_temp1,vper_temp2))+
00287                      (r_initial.pper1*fdot(vper_temp1,vper_temp2));
00288               r_final.ppar2 = (r_initial.ppar2*fdot(vpar_temp1,vpar_temp2))+
00289                      (r_initial.pper2*fdot(vper_temp1,vpar_temp2));
00290               r_final.pper2 = (r_initial.ppar2*fdot(vpar_temp1,vper_temp2))+
00291                      (r_initial.pper2*fdot(vper_temp1,vper_temp2));
00292        }
00293        return r_final;
00294 }
00295 
00296 
00297 static TRAYON
00298 transbetaalpha(
00299        TRAYON r_initial
00300 )
00301 {
00302        /* transforme le rayon r_initial de la base associee a beta dans
00303           la base associee a alpha */ 
00304        TRAYON r_final;
00305        FVECT vpar_temp1,vpar_temp2,vper_temp1,vper_temp2;
00306 
00307        r_final = r_initial;
00308        beta_alpha(r_initial.v,r_final.v);
00309        if ((Y(r_initial) != 0. || Z(r_initial) != 0. )&&(Y(r_final) != 0. || Z(r_final)!= 0.))
00310        {
00311               v_par(r_initial.v,vpar_temp1);
00312               beta_alpha(vpar_temp1,vpar_temp1);
00313               v_per(r_initial.v,vper_temp1);
00314               beta_alpha(vper_temp1,vper_temp1);
00315               v_par(r_final.v,vpar_temp2);
00316               v_per(r_final.v,vper_temp2);
00317               r_final.ppar1 = (r_initial.ppar1*fdot(vpar_temp1,vpar_temp2))+
00318                      (r_initial.pper1*fdot(vper_temp1,vpar_temp2));
00319               r_final.pper1 = (r_initial.ppar1*fdot(vpar_temp1,vper_temp2))+
00320                      (r_initial.pper1*fdot(vper_temp1,vper_temp2));
00321               r_final.ppar2 = (r_initial.ppar2*fdot(vpar_temp1,vpar_temp2))+
00322                      (r_initial.pper2*fdot(vper_temp1,vpar_temp2));
00323               r_final.pper2 = (r_initial.ppar2*fdot(vpar_temp1,vper_temp2))+
00324                      (r_initial.pper2*fdot(vper_temp1,vper_temp2));
00325 
00326        }
00327        return r_final;
00328 }
00329 
00330  
00331 static TRAYON
00332 transalphagamma(
00333        TRAYON r_initial
00334 )
00335 /* transforme le rayon r_initial de la base associee a alpha dans
00336    la base associee a gamma */ 
00337 {
00338        TRAYON r_final;
00339        FVECT vpar_temp1,vpar_temp2,vper_temp1,vper_temp2;
00340 
00341        r_final = r_initial;
00342        alpha_gamma(r_initial.v,r_final.v);
00343        if (( Y(r_initial) != 0. || Z(r_initial) != 0. )&&(Y(r_final)!= 0. || Z(r_final) !=0.))
00344        {
00345               v_par(r_initial.v,vpar_temp1);
00346               alpha_gamma(vpar_temp1,vpar_temp1);
00347               v_per(r_initial.v,vper_temp1);
00348               alpha_gamma(vper_temp1,vper_temp1);
00349               v_par(r_final.v,vpar_temp2);
00350               v_per(r_final.v,vper_temp2); 
00351               r_final.ppar1 = (r_initial.ppar1*fdot(vpar_temp1,vpar_temp2))+
00352                      (r_initial.pper1*fdot(vper_temp1,vpar_temp2));
00353               r_final.pper1 = (r_initial.ppar1*fdot(vpar_temp1,vper_temp2))+
00354                      (r_initial.pper1*fdot(vper_temp1,vper_temp2));
00355               r_final.ppar2 = (r_initial.ppar2*fdot(vpar_temp1,vpar_temp2))+
00356                      (r_initial.pper2*fdot(vper_temp1,vpar_temp2));
00357               r_final.pper2 = (r_initial.ppar2*fdot(vpar_temp1,vper_temp2))+
00358                      (r_initial.pper2*fdot(vper_temp1,vper_temp2));
00359 
00360        }
00361        return r_final;
00362 }
00363 
00364  
00365 static TRAYON
00366 transgammaalpha(
00367        TRAYON r_initial
00368 )
00369 /* transforme le rayon r_initial de la base associee a gamma dans
00370    la base associee a alpha */ 
00371 {
00372        TRAYON r_final;
00373        FVECT vpar_temp1,vpar_temp2,vper_temp1,vper_temp2;
00374 
00375        r_final = r_initial;
00376        gamma_alpha(r_initial.v,r_final.v);
00377        if (( Y(r_initial) != 0. || Z(r_initial) != 0. )&&(Y(r_final) !=0. || Z(r_final) != 0.))
00378        {
00379               v_par(r_initial.v,vpar_temp1);
00380               gamma_alpha(vpar_temp1,vpar_temp1);
00381               v_per(r_initial.v,vper_temp1);
00382               gamma_alpha(vper_temp1,vper_temp1);
00383               v_par(r_final.v,vpar_temp2);
00384               v_per(r_final.v,vper_temp2); 
00385               r_final.ppar1 = (r_initial.ppar1*fdot(vpar_temp1,vpar_temp2))+
00386                      (r_initial.pper1*fdot(vper_temp1,vpar_temp2));
00387               r_final.pper1 = (r_initial.ppar1*fdot(vpar_temp1,vper_temp2))+
00388                      (r_initial.pper1*fdot(vper_temp1,vper_temp2));
00389               r_final.ppar2 = (r_initial.ppar2*fdot(vpar_temp1,vpar_temp2))+
00390                      (r_initial.pper2*fdot(vper_temp1,vpar_temp2));
00391               r_final.pper2 = (r_initial.ppar2*fdot(vpar_temp1,vper_temp2))+
00392                      (r_initial.pper2*fdot(vper_temp1,vper_temp2));
00393        }
00394        return r_final;
00395 }
00396  
00397 
00398 
00399 static int
00400 compare(
00401        TRAYON r1,
00402        TRAYON r2,
00403        double marge
00404 )
00405 {
00406        double arctg1, arctg2;
00407 
00408        arctg1 = atan2(Y(r1),X(r1));
00409        arctg2 = atan2(Y(r2),X(r2));
00410        if ((arctg1 - marge <= arctg2) && (arctg1 + marge >= arctg2)) return 1;
00411        else return 0;
00412 }
00413 
00414 
00415 
00416 
00417 static void
00418 sortie(
00419        TRAYON r
00420 )
00421 {
00422        int i = 0;
00423        int egalite = 0;
00424 
00425 
00426        if(r.e > seuil)
00427        {
00428               while (i < nbrayons && egalite == 0)
00429               {
00430                      raytemp = &ray[i];
00431                      egalite = compare(r,*raytemp,tolerance);
00432                      if (egalite) raytemp->e = raytemp->e + r.e;
00433                      else i = i + 1;
00434               }
00435               if (egalite == 0)
00436               {
00437                      if (nbrayons == 0) ray = (TRAYON *)calloc(1,sizeof(TRAYON));
00438                      else ray = (TRAYON *)realloc((void *)ray, (nbrayons+1)*(sizeof(TRAYON)));         
00439                      if (ray == NULL)
00440                             error(SYSTEM, "out of memory in sortie\n");
00441                      raytemp = &ray[nbrayons];
00442                      raytemp->v[0] = X(r);
00443                      raytemp->v[1] = Y(r);
00444                      raytemp->v[2] = Z(r);
00445                      raytemp->e = r.e;
00446                      nbrayons++;
00447               }
00448        }
00449        return;
00450 }
00451 
00452 
00453 static void
00454 trigo(
00455        TRAYON r
00456 )
00457 /* calcule les grandeurs trigonometriques relatives au rayon incident
00458    et le rapport entre les indices du milieu refracteur et incident   */
00459 {
00460        double det;
00461 
00462        det = Sqrt(X(r)*X(r)+Y(r)*Y(r)+Z(r)*Z(r));
00463        sinus = Sqrt(Y(r)*Y(r)+Z(r)*Z(r))/det;
00464        cosinus = Sqrt(X(r)*X(r))/det;
00465        if (r.n == 1.) rapport = prism.np * prism.np;
00466        else rapport = 1./(prism.np * prism.np);
00467        return;
00468 }
00469 
00470 
00471 static TRAYON
00472 reflexion(
00473        TRAYON r_incident
00474 )
00475 {
00476  /* calcul du rayon reflechi par une face */
00477  TRAYON r_reflechi;
00478 
00479  r_reflechi = r_incident; 
00480  trigo(r_incident);
00481  X(r_reflechi) = -X(r_incident);
00482  Y(r_reflechi) = Y(r_incident);
00483  Z(r_reflechi) = Z(r_incident);
00484  if(sinus > Sqrt(rapport) || r_incident.dest == tot_ref)
00485        {
00486         r_reflechi.ppar1 = r_incident.ppar1;
00487         r_reflechi.pper1 = r_incident.pper1;
00488          r_reflechi.ppar2 = r_incident.ppar2;
00489         r_reflechi.pper2 = r_incident.pper2;
00490         r_reflechi.e = r_incident.e * fact_ref[r_incident.dest];
00491        }
00492  else
00493        { 
00494         r_reflechi.ppar1 = r_incident.ppar1*(rapport*cosinus-Sqrt(rapport-
00495                (sinus*sinus)))/(rapport*cosinus+Sqrt(rapport-(sinus*sinus)));
00496         r_reflechi.pper1 = r_incident.pper1*(cosinus-Sqrt
00497                (rapport-(sinus*sinus)))/(cosinus+Sqrt(rapport-(sinus*sinus)));
00498         r_reflechi.ppar2 = r_incident.ppar2*(rapport*cosinus-Sqrt(rapport-
00499                (sinus*sinus)))/(rapport*cosinus+Sqrt(rapport-(sinus*sinus)));
00500         r_reflechi.pper2 = r_incident.pper2*(cosinus-Sqrt
00501                (rapport-(sinus*sinus)))/(cosinus+Sqrt(rapport-(sinus*sinus)));
00502         r_reflechi.e = r_incident.e *(((r_reflechi.ppar1*r_reflechi.ppar1+
00503          r_reflechi.pper1*r_reflechi.pper1)/(r_incident.ppar1*r_incident.ppar1+
00504         r_incident.pper1*r_incident.pper1))+((r_reflechi.ppar2*r_reflechi.ppar2
00505         +r_reflechi.pper2*r_reflechi.pper2)/(r_incident.ppar2*r_incident.ppar2
00506         +r_incident.pper2*r_incident.pper2)))/2;
00507        }
00508  
00509  /* a la sortie de cette routine r_transmis.orig et .dest ne sont pas definis!*/
00510  return r_reflechi;
00511 }
00512 
00513 
00514 static TRAYON
00515 transmission(
00516        TRAYON r_incident
00517 )
00518 {
00519  /* calcul du rayon refracte par une face */
00520  TRAYON r_transmis;
00521 
00522  r_transmis = r_incident;
00523  trigo(r_incident); 
00524  if (sinus <= Sqrt(rapport) && r_incident.dest != tot_ref) 
00525        {
00526         X(r_transmis) = (X(r_incident)/(fabs(X(r_incident))))*
00527                       (Sqrt(1.-(Y(r_incident)*Y(r_incident)+Z(r_incident)*
00528                                 Z(r_incident))/rapport));
00529         Y(r_transmis) = Y(r_incident)/Sqrt(rapport);
00530         Z(r_transmis) = Z(r_incident)/Sqrt(rapport);
00531         r_transmis.ppar1 = r_incident.ppar1*2.*Sqrt(rapport)*cosinus/
00532                         (Sqrt(rapport-sinus*sinus)+rapport*cosinus);
00533         r_transmis.pper1 = r_incident.pper1*2.*cosinus/(cosinus+Sqrt(rapport
00534                         - sinus*sinus));
00535         r_transmis.ppar2 = r_incident.ppar2*2.*Sqrt(rapport)*cosinus/
00536                         (Sqrt(rapport-sinus*sinus)+rapport*cosinus);
00537         r_transmis.pper2 = r_incident.pper2*2.*cosinus/(cosinus+Sqrt(rapport
00538                         - sinus*sinus));
00539         r_transmis.e = (r_incident.e/2)*(Sqrt(rapport-sinus*sinus)/cosinus)
00540            *(((r_transmis.ppar1*r_transmis.ppar1+r_transmis.pper1*
00541                r_transmis.pper1)
00542            /(r_incident.ppar1*r_incident.ppar1+r_incident.pper1*
00543               r_incident.pper1))+
00544            ((r_transmis.ppar2*r_transmis.ppar2+r_transmis.pper2*r_transmis.pper2)
00545           /(r_incident.ppar2*r_incident.ppar2+r_incident.pper2*r_incident.pper2)));
00546         if(r_incident.n == 1.) r_transmis.n = prism.np;
00547         else r_transmis.n = 1.;
00548        }
00549  else r_transmis.e = 0.;
00550 
00551  /* a la sortie de cette routine r_transmis.orig et .dest ne sont pas definis!*/
00552 
00553  return r_transmis;
00554 }
00555 
00556 
00557 
00558 
00559 #define ensuite(rayon,prob_passage,destination) r_suite = rayon; \
00560                              r_suite.e = prob_passage(rayon)*rayon.e; \
00561                              r_suite.dest = destination; \
00562                              if ( r_suite.e > seuil ) trace_rayon(r_suite)
00563 
00564 
00565 static void
00566 trace_rayon(
00567        TRAYON r_incident
00568 )
00569 {
00570  /* trace le rayon donne */
00571  TRAYON r_reflechi,r_transmis,r_suite;
00572 
00573  switch (r_incident.dest)
00574        {
00575         case ALPHA:
00576               if ( r_incident.orig == ALPHA )
00577                      {
00578                       r_reflechi = reflexion(r_incident);
00579                       sortie(r_reflechi);
00580 
00581                       r_transmis = transmission(r_incident);
00582                       r_transmis.orig = ALPHA;
00583 
00584                       ensuite(r_transmis,prob_alpha_beta,BETA);
00585                       ensuite(r_transmis,prob_alpha_gamma,GAMMA);
00586                      } 
00587               else
00588                      {
00589                       r_reflechi = reflexion(r_incident);
00590                       r_reflechi.orig = ALPHA;
00591                       ensuite(r_reflechi,prob_alpha_beta,BETA);
00592                       ensuite(r_reflechi,prob_alpha_gamma,GAMMA);
00593 
00594                       r_transmis = transmission(r_incident); 
00595                       sortie(r_transmis);
00596                      }
00597               break;
00598         case BETA:
00599           r_reflechi = transbetaalpha(reflexion(transalphabeta(r_incident)));
00600           r_reflechi.orig = BETA;
00601           r_transmis = transbetaalpha(transmission(transalphabeta
00602                      (r_incident)));
00603           r_transmis.orig = GAMMA; 
00604               if ( r_incident.n > 1.0 )  /* le rayon vient de l'interieur */
00605                      {
00606                       ensuite(r_reflechi,prob_beta_alpha,ALPHA);
00607                       ensuite(r_reflechi,prob_beta_gamma,GAMMA); 
00608 
00609                       ensuite(r_transmis,prob_beta_gamma,GAMMA);
00610                       ensuite(r_transmis,prob_beta_delta,DELTA);
00611                      }
00612                else  /* le rayon vient de l'exterieur */
00613                      {
00614                       ensuite(r_reflechi,prob_beta_gamma,GAMMA);
00615                       ensuite(r_reflechi,prob_beta_delta,DELTA);
00616 
00617                       ensuite(r_transmis,prob_beta_alpha,ALPHA);
00618                       ensuite(r_transmis,prob_beta_gamma,GAMMA);
00619                      }
00620                break;
00621         case GAMMA:
00622           r_reflechi = transgammaalpha(reflexion(transalphagamma(r_incident)));
00623           r_reflechi.orig = GAMMA; 
00624           r_transmis = transgammaalpha(transmission(transalphagamma
00625                       (r_incident)));             
00626           r_transmis.orig = GAMMA; 
00627               if ( r_incident.n > 1.0 )  /* le rayon vient de l'interieur */
00628                      {
00629                       ensuite(r_reflechi,prob_gamma_alpha,ALPHA);
00630                       ensuite(r_reflechi,prob_gamma_beta,BETA);
00631 
00632                       ensuite(r_transmis,prob_gamma_beta,BETA);
00633                       ensuite(r_transmis,prob_gamma_delta,DELTA);
00634                      }
00635                else  /* le rayon vient de l'exterieur */
00636                      {
00637                       ensuite(r_reflechi,prob_gamma_beta,BETA);
00638                       ensuite(r_reflechi,prob_gamma_delta,DELTA);
00639 
00640                       ensuite(r_transmis,prob_gamma_alpha,ALPHA);
00641                       ensuite(r_transmis,prob_gamma_beta,BETA);
00642                      }
00643                break;
00644         case DELTA:
00645                if ( r_incident.orig != DELTA ) sortie(r_incident);
00646                else
00647                      {
00648                       ensuite(r_incident,prob_delta_beta,BETA);
00649                       ensuite(r_incident,prob_delta_gamma,GAMMA);
00650                      }
00651                break;
00652        }
00653  return;
00654 }
00655 
00656 #undef ensuite
00657 
00658 static void
00659 inverser(
00660        TRAYON *r1,
00661        TRAYON *r2
00662 )
00663 {
00664        TRAYON temp;
00665        temp = *r1;
00666        *r1 = *r2;
00667        *r2 = temp;
00668 }
00669 
00670 
00671 
00672 static void
00673 setprism(void)
00674 {
00675  double d;
00676  TRAYON r_initial,rsource;
00677  int i,j; 
00678 
00679  prismclock = eclock;
00680 r_initial.ppar1 = r_initial.pper2 = 1.;
00681 r_initial.pper1 = r_initial.ppar2 = 0.;
00682 
00683 d = 1; prism.a = funvalue("arg", 1, &d);
00684  if(prism.a < 0.) goto badopt;
00685 d = 2; prism.b = funvalue("arg", 1, &d);
00686  if(prism.b < 0.) goto badopt;
00687 d = 3; prism.c = funvalue("arg", 1, &d);
00688  if(prism.c < 0.) goto badopt;
00689 d = 4; prism.d = funvalue("arg", 1, &d);
00690  if(prism.d < 0.) goto badopt;
00691 d = 5; prism.np = funvalue("arg", 1, &d);
00692  if(prism.np < 1.) goto badopt;
00693 d = 6; seuil = funvalue("arg", 1, &d);
00694  if (seuil < 0. || seuil >=1) goto badopt;
00695 d = 7; tot_ref = (int)(funvalue("arg", 1, &d) + .5);
00696  if (tot_ref != 1 && tot_ref != 2 && tot_ref != 4) goto badopt;
00697  if (tot_ref < 4 ) 
00698        {
00699         d = 8; fact_ref[tot_ref] = funvalue("arg", 1, &d);
00700         if (fact_ref[tot_ref] < 0. || fact_ref[tot_ref] > 1.) goto badopt;
00701        }
00702 d = 9; tolerance = funvalue("arg", 1, &d);
00703  if (tolerance <= 0.) goto badopt;
00704 d = 10; tolsource = funvalue("arg", 1, &d);
00705  if (tolsource < 0. ) goto badopt;
00706 X(r_initial) = varvalue("Dx");
00707 Y(r_initial) = varvalue("Dy");
00708 Z(r_initial) = varvalue("Dz");
00709 #ifdef DEBUG
00710 fprintf(stderr,"dx=%lf dy=%lf dz=%lf\n",X(r_initial),Y(r_initial),Z(r_initial));
00711 #endif
00712 
00713  /* initialisation */ 
00714  prepare_matrices();
00715  r_initial.e = 1.0;
00716  r_initial.n = 1.0;
00717 
00718 if(ray!=NULL) free(ray);    
00719 nbrayons = 0;
00720 /* determination de l'origine et de la destination du rayon initial */
00721 
00722 if ( X(r_initial) != 0.)
00723  {
00724   if ( X(r_initial) > 0. )
00725        {
00726         r_initial.orig = r_initial.dest = ALPHA;
00727         sens = 1;
00728        }
00729   else if ( X(r_initial) < 0. ) 
00730        {
00731         r_initial.orig = r_initial.dest = DELTA;
00732         sens = -1;
00733        }
00734  
00735   normalize(r_initial.v); 
00736 
00737   trace_rayon(r_initial);
00738 
00739   X(rsource) = varvalue("DxA");
00740   Y(rsource) = varvalue("DyA");
00741   Z(rsource) = varvalue("DzA");
00742   nosource = ( X(rsource)==0. && Y(rsource)==0. && Z(rsource)==0. );
00743   if ( !nosource )
00744        {
00745         for (j=0; j<nbrayons; j++)
00746               {
00747                if ( !compare(ray[j],rsource,tolsource) ) ray[j].e =0.;
00748               }
00749        }
00750   for  (j = 0; j < nbrayons; j++)
00751         {
00752          for (i = j+1; i < nbrayons; i++)
00753          {
00754           if (ray[j].e < ray[i].e) inverser(&ray[j],&ray[i]);
00755          }
00756        }
00757 
00758     bidon = 1;
00759   }
00760  else bidon = 0;
00761  return;
00762 
00763  /* message puis sortie si erreur dans la ligne de commande */
00764  badopt:
00765  bidon = BADVAL;
00766  return;
00767 }
00768 
00769 
00770 static double
00771 l_get_val(
00772        char *nm
00773 )
00774 {
00775  int val, dir, i, trouve, curseur;
00776  int nb;
00777  double valeur;
00778  TRAYON *rayt, raynull;
00779  
00780  if (prismclock < 0 || prismclock < eclock) setprism();
00781  if (bidon == BADVAL) {
00782        errno = EDOM;
00783        return(0.0);
00784  }
00785  val = (int)(argument(1) + .5);
00786  dir = (int)(argument(2) + .5);
00787  nb = (int)(argument(3) + .5);
00788  X(raynull) = bidon;
00789  Y(raynull) = Z(raynull) = 0.;
00790  raynull.e = 0.;
00791  trouve = curseur = 0;
00792  if ( !nosource && nb==2 ) nb=1; /* on est en train de tracer la source
00793                                  a partir de sa seconde source virtuelle */ 
00794 #ifdef DEBUG
00795  fprintf(stderr, " On considere le rayon no: %d\n", nb);
00796 #endif
00797  for(i=0; i < nbrayons &&!trouve; i++)
00798   {
00799    if(ray[i].v[0] * dir * sens >= 0.) curseur ++;
00800    if(curseur == nb)
00801    {
00802     rayt = &ray[i];
00803     trouve = 1;
00804    }
00805   }
00806  if(!trouve) rayt = &raynull; 
00807  switch(val) {
00808        case 0 : valeur = rayt->v[0];
00809                 break;
00810        case 1 : valeur = rayt->v[1];
00811                break;
00812        case 2 : valeur = rayt->v[2];
00813                break;
00814        case 3 : valeur = rayt->e;
00815                break;
00816        default : errno = EDOM; return(0.0);
00817     } 
00818 #ifdef DEBUG
00819   fprintf(stderr, "get_val( %i, %i, %i) = %lf\n",val,dir,nb,valeur);
00820 #endif
00821   return valeur;
00822 }
00823 
00824 
00825 extern void
00826 setprismfuncs(void)  /* declared in func.h */
00827 {
00828  funset("fprism_val", 3, '=', l_get_val);
00829 }