/*************************************************************************************************************************************/
/* */
/* D E F I N I T I O N S R E L A T I V E S A U M I L I E U D E P R O P A G A T I O N : */
/* */
/* */
/* Author of '$xrk/rdn_walk.52$I' : */
/* */
/* Jean-Francois COLONNA (LACTAMME, 1998??????????). */
/* */
/*************************************************************************************************************************************/
/*===================================================================================================================================*/
/*************************************************************************************************************************************/
/* */
/* G E S T I O N D E L A P R O P A G A T I O N : */
/* */
/*************************************************************************************************************************************/
Bblock
Test(IFET(IL_FAUT(utiliser_un_milieu_de_propagation),IZNE(module_de_la_vitesse_incidente)))
Bblock
DEFV(Int,INIT(Xg,X));
DEFV(Int,INIT(Yg,Y));
DEFV(Int,INIT(Zg,Z));
/* Coordonnees du point ou evaluer le gradient : on prend {X,Y,Z} a priori, mais cela */
/* pourra etre modifie si 'IL_FAUT(adapter_les_M_nPAS)' ci-apres en prenant la position */
/* moyenne entre la position courante (a 't') et la position anticipee (a 't+dt'). */
DEFV(genere_Float,INIT(niveau_central_du_pave_du_Mgradient_local_tri_dimensionnel,FLOT__NIVEAU_UNDEF));
/* Valeur du champ "milieur de propagation" au centre du pave de calcul du gradient local. */
/* Ceci a ete introduit le 19971230143719, ainsi que 'ACCES_NIVEAUX_LOCAUX_COURANTS(...)' */
/* afin de savoir si pour une particule donnee ce niveau central change d'un instant a */
/* l'autre. On decide arbitrairement que s'il n'a pas change, il ne peut y avoir refraction. */
DEFV(deltaI_3D,pave_du_Mgradient_local_tri_dimensionnel);
/* Pave a l'interieur duquel evaluer le gradient. */
DEFV(deltaF_3D,Mgradient_local_tri_dimensionnel);
/* Gradient local (Gx,Gy,Gz) du milieu de propagation. */
DEFV(Float,INIT(module_du_Mgradient_local_tri_dimensionnel,FLOT__UNDEF));
DEFV(Float,INIT(Rho_du_Mgradient_local_tri_dimensionnel,FZERO));
/* Module |G| du gradient local (Gx,Gy,Gz) du milieu de propagation. On notera que ce */
/* module est calcule par 'longF3D(...)' ou par 'Rho_3D(...)', suivant les cas, d'ou les */
/* deux noms differents... */
DEFV(Float,INIT(Phi_du_Mgradient_local_tri_dimensionnel,FZERO));
DEFV(Float,INIT(Theta_du_Mgradient_local_tri_dimensionnel,FZERO));
/* Angles {Phi,Theta} du gradient local (Gx,Gy,Gz) du milieu de propagation. */
DEFV(deltaI_3D,delta_anticipe_dans_le_milieu_de_propagation);
/* Deplacement a priori exprime en coordonnees image... */
#define delta_X_anticipe \
ASD1(delta_anticipe_dans_le_milieu_de_propagation,dx)
#define delta_Y_anticipe \
ASD1(delta_anticipe_dans_le_milieu_de_propagation,dy)
#define delta_Z_anticipe \
ASD1(delta_anticipe_dans_le_milieu_de_propagation,dz)
INITIALISATION_ACCROISSEMENT_3D(delta_anticipe_dans_le_milieu_de_propagation
,SOUS(X_anticipe,X)
,SOUS(Y_anticipe,Y)
,SOUS(Z_anticipe,Z)
);
/* Deplacement anticipe de la particule courante (ne pouvant etre nul). */
Test(IL_FAUT(adapter_les_M_nPAS))
Bblock
EGAL(Xg,MOYE(X,X_anticipe));
EGAL(Yg,MOYE(Y,Y_anticipe));
EGAL(Zg,MOYE(Z,Z_anticipe));
/* Coordonnees du point ou evaluer le gradient : on prend la position moyenne entre la */
/* position courante (a 't') et la position anticipee (a 't+dt'). */
INITIALISATION_ACCROISSEMENT_3D(delta_anticipe_dans_le_milieu_de_propagation
,MAX2(pasX,MOIT(ABSO(delta_X_anticipe)))
,MAX2(pasY,MOIT(ABSO(delta_Y_anticipe)))
,MAX2(pasZ,MOIT(ABSO(delta_Z_anticipe)))
);
Eblock
ATes
Bblock
Eblock
ETes
#define ch_M \
milieu_de_propagation
/* ATTENTION, les trois : */
/* */
/* TRON(DIVI(delta_?_anticipe,pas?),MINIMUM_ANTICIPATION_?,MAXIMUM_ANTICIPATION_?) */
/* */
/* qui suivent sont absents de 'v $xrk/rdn_walk.41$K eM_Npas' et ont du etre ajoutes le */
/* 19981023095100 dans 'v $xrk/rdn_walk.51$K eM_Npas' car, en effet, a cause de la */
/* gravitation generalisee, les vitesses peuvent augmenter dans des proportions importantes */
/* et augmenter ainsi de facon inconsideree le volume des parallelepipedes dans lesquels */
/* le gradient est calcule. */
/* ATTENTION, le 19991224132633, 'delta_?_anticipe' a ete remplace dans 'eM_Npas?' par */
/* le maximum entre 'delta_?_anticipe' et le rayon (mis dans l'espace de visualisation) */
/* du 'sorps' courant. L'avantage de cette methode est que le gradient va etre calcule */
/* si 'IL_FAUT(adapter_les_M_nPAS)' dans une boite qui contient la particule (consideree */
/* comme non ponctuelle donc) ce qui va permettre de la "confronter" au milieu avec son */
/* encombrement reel. Ainsi, des "grosses" particules ne risquent plus de donner */
/* l'impression de penetrer dans le milieu comme cela peut se voir dans la sequence : */
/* */
/* xivPdf 9 2 / 025247_025758 */
/* */
/* en particulier sur la derniere image ('025758')... */
/* ATTENTION, le 20000328135804, les valeurs 'MAXIMUM_ANTICIPATION_?' ont ete remplacees */
/* par 'INFINI' si 'IL_FAUT(prendre_en_compte_M_l_encombrement_des_particules)'... */
#define MINIMUM_ANTICIPATION_X \
UN
#define MAXIMUM_ANTICIPATION_X \
QUATRE
#define M_NpasX \
INTE(fM_NpasX)
#define eM_NpasX \
INTE(MUL2(COND(IL_FAUT(adapter_les_M_nPAS) \
,TRON(DIVI(COND(IL_FAUT(prendre_en_compte_M_l_encombrement_des_particules) \
,MAX2(delta_X_anticipe \
,lX_PHYSIQUE_A_VISUALISATION \
(ACCES_LISTE(liste_initiale_des_RAYON_D_INTERACTION,corps)) \
) \
,delta_X_anticipe \
) \
,pasX \
) \
,MINIMUM_ANTICIPATION_X \
,COND(IL_FAUT(prendre_en_compte_M_l_encombrement_des_particules),INFINI,MAXIMUM_ANTICIPATION_X) \
) \
,MINIMUM_ANTICIPATION_X \
) \
,fM_NpasX \
) \
) \
/* ATTENTION, jusqu'au 20011015093459, par erreur, c'etait 'liste_initiale_des_RAYON' qui */ \
/* etait utilise ici a la place de 'liste_initiale_des_RAYON_D_INTERACTION'. Cela a eu des */ \
/* consequences sur pratiquement toutes les sequences generees avant cette date (voir a ce */ \
/* propos 'v _____xivPdf_14_1/$Fnota 011777_012288'), puisqu'en regle generale ces deux */ \
/* parametres different (legerement...). */
#define PreX(x) \
nPREX(x,eM_NpasX)
#define SucX(x) \
nSUCX(x,eM_NpasX)
#define MINIMUM_ANTICIPATION_Y \
UN
#define MAXIMUM_ANTICIPATION_Y \
QUATRE
#define M_NpasY \
INTE(fM_NpasY)
#define eM_NpasY \
INTE(MUL2(COND(IL_FAUT(adapter_les_M_nPAS) \
,TRON(DIVI(COND(IL_FAUT(prendre_en_compte_M_l_encombrement_des_particules) \
,MAX2(delta_Y_anticipe \
,lY_PHYSIQUE_A_VISUALISATION \
(ACCES_LISTE(liste_initiale_des_RAYON_D_INTERACTION,corps)) \
) \
,delta_Y_anticipe \
) \
,pasY \
) \
,MINIMUM_ANTICIPATION_Y \
,COND(IL_FAUT(prendre_en_compte_M_l_encombrement_des_particules),INFINI,MAXIMUM_ANTICIPATION_Y) \
) \
,MINIMUM_ANTICIPATION_Y \
) \
,fM_NpasY \
) \
) \
/* ATTENTION, jusqu'au 20011015093459, par erreur, c'etait 'liste_initiale_des_RAYON' qui */ \
/* etait utilise ici a la place de 'liste_initiale_des_RAYON_D_INTERACTION'. Cela a eu des */ \
/* consequences sur pratiquement toutes les sequences generees avant cette date (voir a ce */ \
/* propos 'v _____xivPdf_14_1/$Fnota 011777_012288'), puisqu'en regle generale ces deux */ \
/* parametres different (legerement...). */
#define PreY(y) \
nPREY(y,eM_NpasY)
#define SucY(y) \
nSUCY(y,eM_NpasY)
#define MINIMUM_ANTICIPATION_Z \
UN
#define MAXIMUM_ANTICIPATION_Z \
QUATRE
#define M_NpasZ \
INTE(fM_NpasZ)
#define eM_NpasZ \
INTE(MUL2(COND(IL_FAUT(adapter_les_M_nPAS) \
,TRON(DIVI(COND(IL_FAUT(prendre_en_compte_M_l_encombrement_des_particules) \
,MAX2(delta_Z_anticipe \
,lZ_PHYSIQUE_A_VISUALISATION \
(ACCES_LISTE(liste_initiale_des_RAYON_D_INTERACTION,corps)) \
) \
,delta_Z_anticipe \
) \
,pasZ \
) \
,MINIMUM_ANTICIPATION_Z \
,COND(IL_FAUT(prendre_en_compte_M_l_encombrement_des_particules),INFINI,MAXIMUM_ANTICIPATION_Z) \
) \
,MINIMUM_ANTICIPATION_Z \
) \
,fM_NpasZ \
) \
) \
/* ATTENTION, jusqu'au 20011015093459, par erreur, c'etait 'liste_initiale_des_RAYON' qui */ \
/* etait utilise ici a la place de 'liste_initiale_des_RAYON_D_INTERACTION'. Cela a eu des */ \
/* consequences sur pratiquement toutes les sequences generees avant cette date (voir a ce */ \
/* propos 'v _____xivPdf_14_1/$Fnota 011777_012288'), puisqu'en regle generale ces deux */ \
/* parametres different (legerement...). */
#define PreZ(z) \
nPREZ(z,eM_NpasZ)
#define SucZ(z) \
nSUCZ(z,eM_NpasZ)
#define FALOAD_POINT(champ,x,y,z) \
______NORMALISE_NIVEAU(FAload_point(champ \
,x,y,z \
,M_periodiser_X,M_periodiser_Y,M_periodiser_Z \
,M_symetriser_X,M_symetriser_Y,M_symetriser_Z \
,M_prolonger_X,M_prolonger_Y,M_prolonger_Z \
,M_niveau_hors_du_milieu_de_propagation \
) \
)
/* Pour reduire la longueur des lignes qui suivent... */
Test(IL_FAUT(affiner_le_maillage_de_calcul_du_M_gradient))
Bblock
DEFV(Logical,INIT(poursuivre_l_affinage_du_maillage_de_calcul_du_M_gradient,VRAI));
#define SURVEILLANCE_DE_L_AFFINAGE_DU_MAILLAGE_DE_CALCUL_DU_M_GRADIENT \
MILLE
DEFV(Int,INIT(surveillance_de_l_affinage_du_maillage_de_calcul_du_M_gradient
,SURVEILLANCE_DE_L_AFFINAGE_DU_MAILLAGE_DE_CALCUL_DU_M_GRADIENT
)
);
#undef SURVEILLANCE_DE_L_AFFINAGE_DU_MAILLAGE_DE_CALCUL_DU_M_GRADIENT
/* Afin d'iterer sur l'affinage... */
DEFV(Float,INIT(FXc,FLOT(X)));
DEFV(Float,INIT(FYc,FLOT(Y)));
DEFV(Float,INIT(FZc,FLOT(Z)));
/* Point (en flottant) a partir duquel on va chercher un gradient non nul, */
DEFV(Int,INIT(Xc,UNDEF));
DEFV(Int,INIT(Yc,UNDEF));
DEFV(Int,INIT(Zc,UNDEF));
/* Puis en entier... */
DEFV(Float,INIT(maximum_des_FpasXYZ,FLOT__UNDEF));
/* Plus grand deplacement en valeur absolu... */
DEFV(Float,INIT(distance_a_parcourir_a_priori_dans_le_milieu_de_propagation,FLOT__UNDEF));
/* Distance a parcourir a priori (c'est-a-dire s'il n'y a ni reflexion, ni refraction...). */
DEFV(Float,INIT(FpasX
,SOUS(X_PHYSIQUE_A_VISUALISATION(AXPB(ASD1(ACCES_VITESSE_COURANTE(corps),dx)
,DCT_EFFECTIF
,ASD1(ACCES_COORDONNEES_COURANTES(corps),x)
)
)
,X_PHYSIQUE_A_VISUALISATION(ASD1(ACCES_COORDONNEES_COURANTES(corps),x))
)
)
);
DEFV(Float,INIT(FpasY
,SOUS(Y_PHYSIQUE_A_VISUALISATION(AXPB(ASD1(ACCES_VITESSE_COURANTE(corps),dy)
,DCT_EFFECTIF
,ASD1(ACCES_COORDONNEES_COURANTES(corps),y)
)
)
,Y_PHYSIQUE_A_VISUALISATION(ASD1(ACCES_COORDONNEES_COURANTES(corps),y))
)
)
);
DEFV(Float,INIT(FpasZ
,SOUS(Z_PHYSIQUE_A_VISUALISATION(AXPB(ASD1(ACCES_VITESSE_COURANTE(corps),dz)
,DCT_EFFECTIF
,ASD1(ACCES_COORDONNEES_COURANTES(corps),z)
)
)
,Z_PHYSIQUE_A_VISUALISATION(ASD1(ACCES_COORDONNEES_COURANTES(corps),z))
)
)
);
/* Pas non entier de deplacement sur les axes dans le sens de la vitesse courante... */
/* ATTENTION, jusqu'au 20011022133229, il y avait ici 'dct' et non pas 'DCT_EFFECTIF' */
/* par erreur... */
EGAL(maximum_des_FpasXYZ,MAX3(ABSO(FpasX),ABSO(FpasY),ABSO(FpasZ)));
EGAL(FpasX,DIVZ(FpasX,maximum_des_FpasXYZ));
EGAL(FpasY,DIVZ(FpasY,maximum_des_FpasXYZ));
EGAL(FpasZ,DIVZ(FpasZ,maximum_des_FpasXYZ));
/* Pas non entier de deplacement sur les axes dans le sens de la vitesse courante... */
EGAL(distance_a_parcourir_a_priori_dans_le_milieu_de_propagation
,longI3D(delta_anticipe_dans_le_milieu_de_propagation)
);
/* Distance a parcouri a priori (c'est-a-dire s'il n'y a ni reflexion, ni refraction...). */
Tant(IL_FAUT(poursuivre_l_affinage_du_maillage_de_calcul_du_M_gradient))
Bblock
EGAL(Xc,INTE(FXc));
EGAL(Yc,INTE(FYc));
EGAL(Zc,INTE(FZc));
/* Point (en entier) ou le gradient va etre calcule... */
#define PReX(x) \
nPREX(x,M_NpasX)
#define SUcX(x) \
nSUCX(x,M_NpasX)
#define PReY(y) \
nPREY(y,M_NpasY)
#define SUcY(y) \
nSUCY(y,M_NpasY)
#define PReZ(z) \
nPREZ(z,M_NpasZ)
#define SUcZ(z) \
nSUCZ(z,M_NpasZ)
EGAL(niveau_central_du_pave_du_Mgradient_local_tri_dimensionnel
,FALOAD_POINT(ch_M,NEUT(Xc),NEUT(Yc),NEUT(Zc))
);
INITIALISATION_ACCROISSEMENT_3D(pave_du_Mgradient_local_tri_dimensionnel
,LENG(PReX(Xc),SUcX(Xc))
,LENG(PReY(Yc),SUcY(Yc))
,LENG(PReZ(Zc),SUcZ(Zc))
);
/* Definition du gradient calcule... */
INITIALISATION_ACCROISSEMENT_3D(Mgradient_local_tri_dimensionnel
,DIVI(SOUS(FALOAD_POINT(ch_M,SUcX(Xc),NEUT(Yc),NEUT(Zc))
,FALOAD_POINT(ch_M,PReX(Xc),NEUT(Yc),NEUT(Zc))
)
,_____lNORMALISE_OX(DOUB(nPAS(M_NpasX,pasX)))
)
,DIVI(SOUS(FALOAD_POINT(ch_M,NEUT(Xc),SUcY(Yc),NEUT(Zc))
,FALOAD_POINT(ch_M,NEUT(Xc),PReY(Yc),NEUT(Zc))
)
,_____lNORMALISE_OY(DOUB(nPAS(M_NpasY,pasY)))
)
,DIVI(SOUS(FALOAD_POINT(ch_M,NEUT(Xc),NEUT(Yc),SUcZ(Zc))
,FALOAD_POINT(ch_M,NEUT(Xc),NEUT(Yc),PReZ(Zc))
)
,_____lNORMALISE_OZ(DOUB(nPAS(M_NpasZ,pasZ)))
)
);
/* Calcul des trois composantes du gradient local au point {X,Y,Z} en ce qui concerne le */
/* milieu de propagation. */
#undef SUcZ
#undef PReZ
#undef SUcY
#undef PReY
#undef SUcX
#undef PReX
EGAL(module_du_Mgradient_local_tri_dimensionnel
,longF3D(Mgradient_local_tri_dimensionnel)
);
/* Calcul du module |G| du gradient local au point courant (Xc,Yc,Zc). */
Test(I3OU(IZGT(module_du_Mgradient_local_tri_dimensionnel)
,IFGT(RdisF3D(FLOT(X),FLOT(Y),FLOT(Z)
,FXc,FYc,FZc
)
,distance_a_parcourir_a_priori_dans_le_milieu_de_propagation
)
,IZLE(surveillance_de_l_affinage_du_maillage_de_calcul_du_M_gradient)
)
)
Bblock
EGAL(poursuivre_l_affinage_du_maillage_de_calcul_du_M_gradient,FAUX);
/* On s'arrete d'affiner soit sur le premier gradient non nul rencontre, soit lorsque l'on */
/* arrive au voisinage du point "anticipe" (en fonction de la vitesse courante et du 'dct'). */
Eblock
ATes
Bblock
INCR(FXc,FpasX);
INCR(FYc,FpasY);
INCR(FZc,FpasZ);
/* Deplacement le long du vecteur vitesse. */
DECR(surveillance_de_l_affinage_du_maillage_de_calcul_du_M_gradient,I);
/* Pour eviter de boucler... */
Eblock
ETes
Eblock
ETan
EGAL(Xg,Xc);
EGAL(Yg,Yc);
EGAL(Zg,Zc);
/* Coordonnees du point ou a ete evalue le gradient. */
Eblock
ATes
Bblock
EGAL(niveau_central_du_pave_du_Mgradient_local_tri_dimensionnel
,FALOAD_POINT(ch_M,NEUT(Xg),NEUT(Yg),NEUT(Zg))
);
/* Definition du gradient calcule... */
Test(IL_FAUT(calculer_un_M_gradient_tridimensionnel_simplifie))
Bblock
INITIALISATION_ACCROISSEMENT_3D(pave_du_Mgradient_local_tri_dimensionnel
,LENG(PreX(Xg),SucX(Xg))
,LENG(PreY(Yg),SucY(Yg))
,LENG(PreZ(Zg),SucZ(Zg))
);
INITIALISATION_ACCROISSEMENT_3D(Mgradient_local_tri_dimensionnel
,DIVI(SOUS(FALOAD_POINT(ch_M,SucX(Xg),NEUT(Yg),NEUT(Zg))
,FALOAD_POINT(ch_M,PreX(Xg),NEUT(Yg),NEUT(Zg))
)
,_____lNORMALISE_OX(DOUB(nPAS(eM_NpasX,pasX)))
)
,DIVI(SOUS(FALOAD_POINT(ch_M,NEUT(Xg),SucY(Yg),NEUT(Zg))
,FALOAD_POINT(ch_M,NEUT(Xg),PreY(Yg),NEUT(Zg))
)
,_____lNORMALISE_OY(DOUB(nPAS(eM_NpasY,pasY)))
)
,DIVI(SOUS(FALOAD_POINT(ch_M,NEUT(Xg),NEUT(Yg),SucZ(Zg))
,FALOAD_POINT(ch_M,NEUT(Xg),NEUT(Yg),PreZ(Zg))
)
,_____lNORMALISE_OZ(DOUB(nPAS(eM_NpasZ,pasZ)))
)
);
/* Calcul des trois composantes du gradient local au point {X,Y,Z} en ce qui concerne le */
/* milieu de propagation. Ce calcul est fait de trois calculs monodimensionnels. */
Eblock
ATes
Bblock
DEFV(Int,INIT(nombre_de_points_dans_le_pave,ZERO));
/* Nombre de points du volume dans lequel on va rechercher la moyenne des niveaux.. */
DEFV(genere_Float,INIT(niveau_moyen_dans_le_pave,FZERO));
/* Afin de calculer le niveau moyen dans le pave de calcul du gradient (dans [0,1]). */
DEFV(genere_Float,INIT(niveau_minimum_tolerable,FLOT__NIVEAU_UNDEF));
DEFV(genere_Float,INIT(niveau_maximum_tolerable,FLOT__NIVEAU_UNDEF));
/* Segment dans lequel on tolere les niveaux pour le calcul du gradient... */
DEFV(Int,INIT(nombre_de_points_du_Mgradient_local,ZERO));
/* Nombre de points du volume dans lequel on va moyenner des gradients "ponctuels". */
DEFV(Float,INIT(sphere_dans_le_pave_du_Mgradient_local_tri_dimensionnel,FLOT__UNDEF));
/* Afin de connaitre le diametre de la sphere inscrite dans le pave de calcul du gradient. */
DEFV(deltaF_3D,rayon_vecteur_courant);
DEFV(Float,INIT(module_du_rayon_vecteur_courant,FLOT__UNDEF));
/* Rayon vecteur du point courant {X,Y,Z} et son module. */
INITIALISATION_ACCROISSEMENT_3D(pave_du_Mgradient_local_tri_dimensionnel
,LENG(PREX(PreX(Xg)),SUCX(SucX(Xg)))
,LENG(PREY(PreY(Yg)),SUCY(SucY(Yg)))
,LENG(PREZ(PreZ(Zg)),SUCZ(SucZ(Zg)))
);
EGAL(sphere_dans_le_pave_du_Mgradient_local_tri_dimensionnel
,MOIT(MAX3(ASD1(pave_du_Mgradient_local_tri_dimensionnel,dx)
,ASD1(pave_du_Mgradient_local_tri_dimensionnel,dy)
,ASD1(pave_du_Mgradient_local_tri_dimensionnel,dz)
)
)
);
/* ATTENTION, jusqu'au 20000328135804 il y avait ici 'MIN3(...)', mais c'est evidemment */
/* 'MAX3(...)' qu'il faut utiliser... */
Test(IL_FAUT(M_gradient_tridimensionnel_non_simplifie_filtrer))
Bblock
begin_albumQ(DoIn,PreZ(Zg),SucZ(Zg),pasZ
,DoIn,PreY(Yg),SucY(Yg),pasY
,DoIn,PreX(Xg),SucX(Xg),pasX
)
Bblock
INITIALISATION_ACCROISSEMENT_3D(rayon_vecteur_courant
,FLOT(SOUS(X,Xg))
,FLOT(SOUS(Y,Yg))
,FLOT(SOUS(Z,Zg))
);
EGAL(module_du_rayon_vecteur_courant,longF3D(rayon_vecteur_courant));
/* Calcul du rayon vecteur du point courant {X,Y,Z} et de son module. */
Test(IFOU(IL_NE_FAUT_PAS(M_gradient_tridimensionnel_non_simplifie_spherique)
,IFET(IL_FAUT(M_gradient_tridimensionnel_non_simplifie_spherique)
,IFLE(module_du_rayon_vecteur_courant
,sphere_dans_le_pave_du_Mgradient_local_tri_dimensionnel
)
)
)
)
Bblock
INCR(niveau_moyen_dans_le_pave,FALOAD_POINT(ch_M,X,Y,Z));
/* Cumul des niveaux... */
INCR(nombre_de_points_dans_le_pave,I);
/* Comptage des points utilises... */
Eblock
ATes
Bblock
Eblock
ETes
Eblock
end_albumQ(EDoI,EDoI,EDoI)
Test(IZGT(nombre_de_points_dans_le_pave))
Bblock
EGAL(niveau_moyen_dans_le_pave
,DIVI(niveau_moyen_dans_le_pave,FLOT(nombre_de_points_dans_le_pave))
);
/* Et normalisation du cumul des niveaux... */
Eblock
ATes
Bblock
PRINT_ERREUR("le nombre de niveaux recuperes est nul");
CAL1(Prer1("rayon de la sphere = %+f\n"
,sphere_dans_le_pave_du_Mgradient_local_tri_dimensionnel
)
);
Eblock
ETes
EGAL(niveau_minimum_tolerable
,SOUS(niveau_moyen_dans_le_pave
,M_gradient_tridimensionnel_non_simplifie_tolerance
)
);
EGAL(niveau_maximum_tolerable
,ADD2(niveau_moyen_dans_le_pave
,M_gradient_tridimensionnel_non_simplifie_tolerance
)
);
/* Definition du segment d'acceptation des niveaux qui participent au calcul du gradient. */
/* */
/* ATTENTION, le 19980914135745, en preparant la sequence : */
/* */
/* xivPdf 11 2 / 029972_030483 */
/* */
/* j'ai pu constater que cela pouvait etre tres ennuyeux, meme sur des images "milieu" ne */
/* contenant que du 'NOIR' et du 'BLANC'. En effet, imaginons qu'il n'y ait, par exemple, */
/* qu'un seul point 'NOIR' ; le niveau moyen est alors tres proche de 'BLANC' et si le */
/* facteur de tolerance est inferieur a 1, le seul point 'NOIR' sera ignore lors du calcul */
/* du gradient, et celui-ci sera donc de module nul... */
Eblock
ATes
Bblock
Eblock
ETes
INITIALISATION_ACCROISSEMENT_3D(Mgradient_local_tri_dimensionnel
,FZERO
,FZERO
,FZERO
);
/* Initialisation du gradient local au point (Xg,Yg,Zg). Ceci doit etre fait quel que soit */
/* le mode 'calculer_la_moyenne_des_Mgradient_3x3x3_tri_dimensionnel' car, en effet, lorsque */
/* le milieu ne change pas il faut que le gradient soit {0,0,0}, alors que l'on ne cumulera */
/* aucun point si 'IL_NE_FAUT_PAS(calculer_la_moyenne_des_Mgradient_3x3x3_tri_dimensionnel'. */
begin_albumQ(DoIn,PreZ(Zg),SucZ(Zg),pasZ
,DoIn,PreY(Yg),SucY(Yg),pasY
,DoIn,PreX(Xg),SucX(Xg),pasX
)
/* Balayage d'un volume dont le point (Xg,Yg,Zg) est le centre. */
/* */
/* ATTENTION, le 19971229141908, afin d'eviter les refractions parasites, j'ai essaye de */
/* mettre ici : */
/* */
/* begin_albumQ(DoIn,SUCZ(PreZ(Zg)),PREZ(SucZ(Zg)),pasZ */
/* ,DoIn,SUCY(PreY(Yg)),PREY(SucY(Yg)),pasY */
/* ,DoIn,SUCX(PreX(Xg)),PREX(SucX(Xg)),pasX */
/* ) */
/* */
/* Malheureusement, cela a provoque la perte de nombreuses reflexions. Je reviens donc a */
/* la solution anterieure... */
Bblock
DEFV(deltaF_3D,Mgradient_3x3x3_tri_dimensionnel);
/* Gradient "ponctuel" au point courant. */
DEFV(Float,INIT(ponderation_du_Mgradient_3x3x3,FU));
/* Ponderation des gradients ponctuels lors du cumul du gradient local. */
INITIALISATION_ACCROISSEMENT_3D(rayon_vecteur_courant
,FLOT(SOUS(X,Xg))
,FLOT(SOUS(Y,Yg))
,FLOT(SOUS(Z,Zg))
);
EGAL(module_du_rayon_vecteur_courant,longF3D(rayon_vecteur_courant));
/* Calcul du rayon vecteur du point courant {X,Y,Z} et de son module. */
Test(IFOU(IL_NE_FAUT_PAS(M_gradient_tridimensionnel_non_simplifie_spherique)
,IFET(IL_FAUT(M_gradient_tridimensionnel_non_simplifie_spherique)
,IFLE(module_du_rayon_vecteur_courant
,sphere_dans_le_pave_du_Mgradient_local_tri_dimensionnel
)
)
)
)
Bblock
DEFV(genere_Float,INIT(niveau_sXnYnZ,FALOAD_POINT(ch_M,SUCX(X),NEUT(Y),NEUT(Z))));
DEFV(genere_Float,INIT(niveau_pXnYnZ,FALOAD_POINT(ch_M,PREX(X),NEUT(Y),NEUT(Z))));
DEFV(genere_Float,INIT(niveau_nXsYnZ,FALOAD_POINT(ch_M,NEUT(X),SUCY(Y),NEUT(Z))));
DEFV(genere_Float,INIT(niveau_nXpYnZ,FALOAD_POINT(ch_M,NEUT(X),PREY(Y),NEUT(Z))));
DEFV(genere_Float,INIT(niveau_nXnYsZ,FALOAD_POINT(ch_M,NEUT(X),NEUT(Y),SUCZ(Z))));
DEFV(genere_Float,INIT(niveau_nXnYpZ,FALOAD_POINT(ch_M,NEUT(X),NEUT(Y),PREZ(Z))));
/* Acces aux 6 points utiles au calcul du gradient tres tres local... */
Test(IFOU(IL_NE_FAUT_PAS(M_gradient_tridimensionnel_non_simplifie_filtrer)
,IFET(IL_FAUT(M_gradient_tridimensionnel_non_simplifie_filtrer)
,I3ET(IFET(IFINff(niveau_sXnYnZ
,niveau_minimum_tolerable
,niveau_maximum_tolerable
)
,IFINff(niveau_pXnYnZ
,niveau_minimum_tolerable
,niveau_maximum_tolerable
)
)
,IFET(IFINff(niveau_nXsYnZ
,niveau_minimum_tolerable
,niveau_maximum_tolerable
)
,IFINff(niveau_nXpYnZ
,niveau_minimum_tolerable
,niveau_maximum_tolerable
)
)
,IFET(IFINff(niveau_nXnYsZ
,niveau_minimum_tolerable
,niveau_maximum_tolerable
)
,IFINff(niveau_nXnYpZ
,niveau_minimum_tolerable
,niveau_maximum_tolerable
)
)
)
)
)
)
Bblock
INITIALISATION_ACCROISSEMENT_3D(Mgradient_3x3x3_tri_dimensionnel
,DIVZ(SOUS(niveau_sXnYnZ,niveau_pXnYnZ)
,_____lNORMALISE_OX(DOUB(pasX))
)
,DIVZ(SOUS(niveau_nXsYnZ,niveau_nXpYnZ)
,_____lNORMALISE_OY(DOUB(pasY))
)
,DIVZ(SOUS(niveau_nXnYsZ,niveau_nXnYpZ)
,_____lNORMALISE_OZ(DOUB(pasZ))
)
);
/* Gradient "ponctuel" au point courant {X,Y,Z}. */
Test(IL_FAUT(ponderer_un_M_gradient_tridimensionnel_non_simplifie))
Bblock
DEFV(Float,INIT(cosinus_du_rayon_vecteur_courant,FLOT__UNDEF));
/* Cosinus de l'angle entre le rayon vecteur du point courant et le vecteur vitesse. */
Test(IZNE(module_du_rayon_vecteur_courant))
Bblock
EGAL(cosinus_du_rayon_vecteur_courant
,DIVI(prdF3D(ACCES_VITESSE_COURANTE(corps)
,rayon_vecteur_courant
)
,MUL2(module_de_la_vitesse_incidente
,module_du_rayon_vecteur_courant
)
)
);
/* Cosinus de l'angle entre le rayon vecteur du point courant et le vecteur vitesse. */
EGAL(ponderation_du_Mgradient_3x3x3
,COS1(cosinus_du_rayon_vecteur_courant)
);
/* La ponderation est ainsi dans [0,1] et fonction de la distance angulaire du point courant */
/* {X,Y,Z} a l'axe du vecteur vitesse (via 'cosinus(...)'). Plus cette distance est faible, */
/* et plus la ponderation est proche de 1 ; plus cette distance est elevee, et plus la */
/* ponderation est proche de 0. */
Eblock
ATes
Bblock
Eblock
ETes
Eblock
ATes
Bblock
Eblock
ETes
Test(IL_FAUT(calculer_la_moyenne_des_Mgradient_3x3x3_tri_dimensionnel))
Bblock
INITIALISATION_ACCROISSEMENT_3D(Mgradient_local_tri_dimensionnel
,AXPB(ponderation_du_Mgradient_3x3x3
,ASD1(Mgradient_3x3x3_tri_dimensionnel,dx)
,ASD1(Mgradient_local_tri_dimensionnel,dx)
)
,AXPB(ponderation_du_Mgradient_3x3x3
,ASD1(Mgradient_3x3x3_tri_dimensionnel,dy)
,ASD1(Mgradient_local_tri_dimensionnel,dy)
)
,AXPB(ponderation_du_Mgradient_3x3x3
,ASD1(Mgradient_3x3x3_tri_dimensionnel,dz)
,ASD1(Mgradient_local_tri_dimensionnel,dz)
)
);
/* Cumul d'obtention du gradient local au point (Xg,Yg,Zg). */
INCR(nombre_de_points_du_Mgradient_local,I);
/* Comptage de tous les points utilises. */
Eblock
ATes
Bblock
DEFV(Float,INIT(Rho_du_Mgradient_3x3x3_tri_dimensionnel,FLOT__UNDEF));
DEFV(Float,INIT(Phi_du_Mgradient_3x3x3_tri_dimensionnel,FLOT__UNDEF));
DEFV(Float,INIT(Theta_du_Mgradient_3x3x3_tri_dimensionnel,FLOT__UNDEF));
/* Coordonnees spheriques {Rho,Phi,Theta} du gradient "ponctuel" au point courant. */
EGAL(Rho_du_Mgradient_3x3x3_tri_dimensionnel
,Rho_3D(ASD1(Mgradient_3x3x3_tri_dimensionnel,dx)
,ASD1(Mgradient_3x3x3_tri_dimensionnel,dy)
,ASD1(Mgradient_3x3x3_tri_dimensionnel,dz)
)
);
Test(IZNE(Rho_du_Mgradient_3x3x3_tri_dimensionnel))
Bblock
EGAL(Phi_du_Mgradient_3x3x3_tri_dimensionnel
,Phi_3D(ASD1(Mgradient_3x3x3_tri_dimensionnel,dx)
,ASD1(Mgradient_3x3x3_tri_dimensionnel,dy)
,ASD1(Mgradient_3x3x3_tri_dimensionnel,dz)
)
);
EGAL(Theta_du_Mgradient_3x3x3_tri_dimensionnel
,Theta_3D(ASD1(Mgradient_3x3x3_tri_dimensionnel,dx)
,ASD1(Mgradient_3x3x3_tri_dimensionnel,dy)
,ASD1(Mgradient_3x3x3_tri_dimensionnel,dz)
)
);
/* Coordonnees spheriques {Rho,Phi,Theta} du gradient "ponctuel" au point courant. */
INCR(Rho_du_Mgradient_local_tri_dimensionnel
,MUL2(ponderation_du_Mgradient_3x3x3
,Rho_du_Mgradient_3x3x3_tri_dimensionnel
)
);
INCR(Phi_du_Mgradient_local_tri_dimensionnel
,MUL2(ponderation_du_Mgradient_3x3x3
,Phi_du_Mgradient_3x3x3_tri_dimensionnel
)
);
INCR(Theta_du_Mgradient_local_tri_dimensionnel
,MUL2(ponderation_du_Mgradient_3x3x3
,Theta_du_Mgradient_3x3x3_tri_dimensionnel
)
);
/* Cumul d'obtention des coordonnees spheriques {Rho,Phi,Theta} du gradient local (Gx,Gy,Gz) */
/* du milieu de propagation. */
INCR(nombre_de_points_du_Mgradient_local,I);
/* Comptage des points utilises et "non vides" (c'est-a-dire dont le 'Rho' n'est pas nul. */
/* Cette technique permet, par exemple, de garantir que les mouvements confines dans des */
/* plans paralleles a {OX,OY} (a 'Z' constant) le restent... */
Eblock
ATes
Bblock
Eblock
ETes
Eblock
ETes
Eblock
ATes
Bblock
Eblock
ETes
Eblock
ATes
Bblock
Eblock
ETes
Eblock
end_albumQ(EDoI,EDoI,EDoI)
Test(IZGT(nombre_de_points_du_Mgradient_local))
Bblock
Test(IL_FAUT(calculer_la_moyenne_des_Mgradient_3x3x3_tri_dimensionnel))
Bblock
INITIALISATION_ACCROISSEMENT_3D(Mgradient_local_tri_dimensionnel
,DIVI(ASD1(Mgradient_local_tri_dimensionnel,dx)
,FLOT(nombre_de_points_du_Mgradient_local)
)
,DIVI(ASD1(Mgradient_local_tri_dimensionnel,dy)
,FLOT(nombre_de_points_du_Mgradient_local)
)
,DIVI(ASD1(Mgradient_local_tri_dimensionnel,dz)
,FLOT(nombre_de_points_du_Mgradient_local)
)
);
/* Normalisation du cumul de definition du gradient local au point (Xg,Yg,Zg). */
Eblock
ATes
Bblock
EGAL(Rho_du_Mgradient_local_tri_dimensionnel
,NEUT(DIVI(Rho_du_Mgradient_local_tri_dimensionnel
,FLOT(nombre_de_points_du_Mgradient_local)
)
)
);
EGAL(Phi_du_Mgradient_local_tri_dimensionnel
,CERC(DIVI(Phi_du_Mgradient_local_tri_dimensionnel
,FLOT(nombre_de_points_du_Mgradient_local)
)
)
);
EGAL(Theta_du_Mgradient_local_tri_dimensionnel
,CERC(DIVI(Theta_du_Mgradient_local_tri_dimensionnel
,FLOT(nombre_de_points_du_Mgradient_local)
)
)
);
/* Coordonnees spheriques {Rho,Phi,Theta} du gradient local (Gx,Gy,Gz) du milieu de */
/* propagation. */
INITIALISATION_ACCROISSEMENT_3D(Mgradient_local_tri_dimensionnel
,Xcartesienne_3D(Rho_du_Mgradient_local_tri_dimensionnel
,Phi_du_Mgradient_local_tri_dimensionnel
,Theta_du_Mgradient_local_tri_dimensionnel
)
,Ycartesienne_3D(Rho_du_Mgradient_local_tri_dimensionnel
,Phi_du_Mgradient_local_tri_dimensionnel
,Theta_du_Mgradient_local_tri_dimensionnel
)
,Zcartesienne_3D(Rho_du_Mgradient_local_tri_dimensionnel
,Phi_du_Mgradient_local_tri_dimensionnel
,Theta_du_Mgradient_local_tri_dimensionnel
)
);
/* Calcul du gradient local au point (Xg,Yg,Zg). ATTENTION : cette methode a un gros defaut. */
/* En effet, les petites imprecisions sur {Phi,Theta} peuvent provoquer des anomalies ; par */
/* exemple, dans une simulation bidimensionnelle a 'Z' constant, les particules peuvent */
/* ainsi s'echapper de leur plan 'Z' de depart. Ce probleme est evidemment cause par les */
/* passages des coordonnees cartesiennes aux coordonnees spheriques, puis retour (cela peut */
/* etre vu grace au programme 'v $xtKi/CartSph3D.01$K'). Ce probleme est corrigeable */
/* approximativement via 'seuil_de_Mgradient_local_tri_dimensionnel_?' ci-apres. */
/* */
/* On notera au passage qu'un plan de coordonnee 'Z' constante a un 'Theta' egal a pi/2. */
INITIALISATION_ACCROISSEMENT_3D(Mgradient_local_tri_dimensionnel
,TROP_PETIT(ASD1(Mgradient_local_tri_dimensionnel,dx)
,seuil_de_Mgradient_local_tri_dimensionnel_X
)
,TROP_PETIT(ASD1(Mgradient_local_tri_dimensionnel,dy)
,seuil_de_Mgradient_local_tri_dimensionnel_Y
)
,TROP_PETIT(ASD1(Mgradient_local_tri_dimensionnel,dz)
,seuil_de_Mgradient_local_tri_dimensionnel_Z
)
);
/* Afin de corriger le probleme evoque ci-dessus au-sujet des aller-retours entre les */
/* coordonnees cartesiennes et spheriques, dont le produit n'est pas la transformation */
/* unite a cause des erreurs d'arrondi. Lorsque l'une des composantes est trop petite, */
/* on considere qu'en fait elle doit etre nulle et donc on l'annule... */
Eblock
ETes
Eblock
ATes
Bblock
Test(IL_FAUT(calculer_la_moyenne_des_Mgradient_3x3x3_tri_dimensionnel))
Bblock
PRINT_ERREUR("lors du filtrage, le nombre de niveaux toleres est nul");
CAL1(Prer1("rayon de la sphere = %+f\n"
,sphere_dans_le_pave_du_Mgradient_local_tri_dimensionnel
)
);
CAL1(Prer3("niveaux = {moyen=%+f,minimum=%+f,maximum=%+f}\n"
,niveau_moyen_dans_le_pave
,niveau_minimum_tolerable
,niveau_maximum_tolerable
)
);
Eblock
ATes
Bblock
Eblock
ETes
Eblock
ETes
Eblock
ETes
Eblock
ETes
#undef FALOAD_POINT
#undef SucZ
#undef PreZ
#undef eM_NpasZ
#undef M_NpasZ
#undef MAXIMUM_ANTICIPATION_Z
#undef MINIMUM_ANTICIPATION_Z
#undef SucY
#undef PreY
#undef eM_NpasY
#undef M_NpasY
#undef MAXIMUM_ANTICIPATION_Y
#undef MINIMUM_ANTICIPATION_Y
#undef SucX
#undef PreX
#undef eM_NpasX
#undef M_NpasX
#undef MAXIMUM_ANTICIPATION_X
#undef MINIMUM_ANTICIPATION_X
#undef ch_M
#undef delta_Z_anticipe
#undef delta_Y_anticipe
#undef delta_X_anticipe
EGAL(module_du_Mgradient_local_tri_dimensionnel
,longF3D(Mgradient_local_tri_dimensionnel)
);
/* Calcul du module |G| du gradient local au point {X,Y,Z} du milieu de propagation. */
EDITION_E(BLOC(CAL2(Prin3(" MILIEU gradient={%+f,%+f,%+f}"
,ASD1(Mgradient_local_tri_dimensionnel,dx)
,ASD1(Mgradient_local_tri_dimensionnel,dy)
,ASD1(Mgradient_local_tri_dimensionnel,dz)
)
);
)
);
EDITION_E(BLOC(CAL2(Prin1(" module(gradient)=%+f"
,module_du_Mgradient_local_tri_dimensionnel
)
);
)
);
EDITION_E(BLOC(CAL2(Prin3(" au point={%d,%d,%d}"
,Xg
,Yg
,Zg
)
);
)
);
EDITION_E(BLOC(CAL2(Prin3(" dans un pave={%d,%d,%d}"
,ASD1(pave_du_Mgradient_local_tri_dimensionnel,dx)
,ASD1(pave_du_Mgradient_local_tri_dimensionnel,dy)
,ASD1(pave_du_Mgradient_local_tri_dimensionnel,dz)
)
);
)
);
/* On notera que grace a cette edition, il est possible de generer une visualisation */
/* tridimensionnelle du milieu de propagation. En effet soit 'MILIEU(n)' la suite des */
/* images definissant ce milieu. On fera alors : */
/* */
/* $xci/dilate.01$X A=MILIEU(n) dilater=FAUX $formatI | \ */
/* $xci/eor$X A1=MILIEU(n) R=CONTOUR(n) $formatI */
/* */
/* ce qui donne le 'CONTOUR' du 'MILIEU' de propagation. Puis : */
/* */
/* $xci/liste_points$X A=CONTOUR(n) eX=VRAI eY=FAUX eNIVEAU=FAUX epoints=FAUX | \ */
/* $SE -e "s/^.*=//" \ */
/* > F1(n)$COORD_X */
/* $xci/liste_points$X A=CONTOUR(n) eX=FAUX eY=VRAI eNIVEAU=FAUX epoints=FAUX | \ */
/* $SE -e "s/^.*=//" \ */
/* > F1(n)$COORD_Y */
/* */
/* ce qui donne les 'NC' coordonnees {X,Y} de 'CONTOUR(n)'. Puis on extrait 1 coordonnee */
/* sur 'N' par : */
/* */
/* $xrv/un_sur_N.01$X ATTENTION=FAUX ne=NC fichier=F1(n)$COORD_X paquets=N | \ */
/* $SE -e "s/^.*=//" \ */
/* > F2(n)$COORD_X */
/* $xrv/un_sur_N.01$X ATTENTION=FAUX ne=NC fichier=F1(n)$COORD_Y paquets=N | \ */
/* $SE -e "s/^.*=//" \ */
/* > F2(n)$COORD_Y */
/* */
/* ce qui donne un echantillonnage regulier de 'NE' coordonnees {X,Y} de 'CONTOUR(n)'. Soit */
/* alors 'Z(n)' la coordonnee 'Z' associee a la couche 'MILIEU(n)'. Il suffit de generer : */
/* */
/* $xci/valeurs_inte$X premiere=1 derniere=NE vD=Z(n) vA=Z(n) \ */
/* > F2(n)$COORD_Z */
/* */
/* On dispose alors de trois series de fichiers {F2(n)$COORD_X,F2(n)$COORD_Y,F2(n)$COORD_Z} */
/* qui definissent de facon tridimensionnelle la frontiere du milieu de propagation. On peut */
/* alors la visualiser a l'aide de : */
/* */
/* $xrv/particule.10$X */
/* */
/* mais on peut aussi injecter {F2(n)$COORD_X,F2(n)$COORD_Y,F2(n)$COORD_Z} comme conditions */
/* initiales dans : */
/* */
/* $xrk/rdn_walk.52$X */
/* */
/* et editer ici le gradient (en ne faisant qu'une seule image, visualisant donc les */
/* conditions initiales) que l'on visualisera a son tour grace a : */
/* */
/* $xrv/particule.10$X */
/* */
/* en le materialisant, par exemple, avec le RAYON des particules... */
EDITION_E(BLOC(CAL2(Prin1(" niveau central du milieu=%+f"
,niveau_central_du_pave_du_Mgradient_local_tri_dimensionnel
)
);
)
);
Test(IZGT(module_du_Mgradient_local_tri_dimensionnel))
Bblock
DEFV(Float,INIT(module_de_la_vitesse_apres_reflexion_ou_refraction,FLOT__UNDEF));
/* Module |V| du vecteur vitesse apres reflexion ou refraction. */
DEFV(deltaF_3D,vecteur_directeur_OX2);
DEFV(deltaF_3D,vecteur_directeur_OY2);
DEFV(Float,INIT(module_du_vecteur_directeur_OY2,FLOT__UNDEF));
DEFV(deltaF_3D,vecteur_directeur_OZ2);
/* Vecteurs unitaires des axes {OX2,OY2,OZ2}. */
DEFV(matrixF_3D,matrice_de_rotation_directe);
DEFV(matrixF_3D,matrice_de_rotation_inverse);
/* Definition de la matrice de passage {OX1,OY1,OZ1} --> {OX2,OY2,OZ2} et inverse. */
DEFV(deltaF_3D,vecteur_vitesse_incident);
/* Vecteur vitesse incident exprime dans le referentiel {OX2,OY2,OZ2}. */
DEFV(deltaF_3D,vecteur_vitesse_reflechi_ou_refracte);
DEFV(Float,INIT(Rho_du_vecteur_vitesse_incident,FLOT__UNDEF));
DEFV(Float,INIT(Theta_du_vecteur_vitesse_incident,FLOT__UNDEF));
/* Vecteur vitesse reflechi ou refracte exprime dans le referentiel {OX2,OY2,OZ2} et */
/* son {rho,theta} exprime dans le plan {OX2,OZ2}. */
DEFV(Float,INIT(cosinus_du_theta_incident_VG,FLOT__UNDEF));
DEFV(Float,INIT(sinus_du_theta_incident_VG,FLOT__UNDEF));
DEFV(Float,INIT(theta_incident_VG,FLOT__UNDEF));
/* Angle theta entre le vecteur vitesse incident et le gradient du milieu de propagation. */
DEFV(Float,INIT(theta_apres_reflexion_ou_refraction,FLOT__UNDEF));
/* Angle theta entre le vecteur vitesse reflechi et le gradient du milieu de propagation. */
DEFV(Float,INIT(rapport_de_l_indice_refracte_a_l_indice_incident,FU));
/* Lorsqu'il y a refraction, le module de la vitesse 'V' est divise par le rapport des */
/* indices de refraction N(refracte)/N(incident) qui est superieur a 1. Lorsqu'il y a */
/* reflexion, le module de 'V' ne change pas, d'ou cette initialisation a 1... */
DEFV(Logical,INIT(immobiliser_le_corps_courant_suite_a_un_test_probabiliste,FAUX));
/* A priori, le corps courant ne doit pas etre immobilise de facon probabiliste... */
/* Le 20020219142536 a ete introduit la possibilite d'immobiliser un corps apres une */
/* reflexion, non plus uniquement avec un decompteur ('ACCES_IMMOBILISABLES(corps)>=0'), */
/* mais, de plus, avec une probabilite donnee par '|ACCES_IMMOBILISABLES(corps)|'... */
Test(IZLT(ACCES_IMMOBILISABLES(corps)))
/* Cas ou un test d'immobilisation probabiliste est demande, mais evidemment uniquement */
/* dans le cas ou une reflexion est rencontree ci-apres... */
Bblock
DEFV(Float,INIT(tirage_aleatoire_d_immobilisation,FLOT__UNDEF));
GENERATION_D_UNE_VALEUR(tirage_aleatoire_d_immobilisation
,PROBABILITE_NULLE
,PROBABILITE_UNITE
);
/* On notera que 'tirage_aleatoire_d_immobilisation', de meme que l'indicateur */
/* 'immobiliser_le_corps_courant_suite_a_un_test_probabiliste' qui est evalue ci-apres, */
/* n'ont d'utilite que s'il y a reflexion du corps courant par la suite. Ce calcul est */
/* fait malgre tout dans tous les cas afin de simplifier les deux tests differents et */
/* exclusifs l'un de l'autre qui seront eventuellement effectues par la suite... */
Test(IFLT(tirage_aleatoire_d_immobilisation,ABSO(ACCES_IMMOBILISABLES(corps))))
Bblock
Test(IFOU(IL_NE_FAUT_PAS(tenter_de_synchroniser_les_immobilisations_sur_les_naissances)
,IFET(IL_FAUT(tenter_de_synchroniser_les_immobilisations_sur_les_naissances)
,IFLT(compteur_des_particules_immobilisees_lors_d_une_collision_avec_le_milieu
,INTE(MUL2(facteur_synchronisation_des_immobilisations_sur_les_naissances
,FLOT(compteur_des_particules_nees_apres_l_instant_initial)
)
)
)
)
)
)
Bblock
/* Une tentative de synchronisation des morts sur les naissances a ete introduite le */
/* 20020227092012. D'une part, celle-ci ne s'applique que si le test d'immobilisation */
/* est probabiliste ('IZLT(ACCES_IMMOBILISABLES(corps))'). D'autre part, elle consiste */
/* uniquement a refuser une immobilisation s'il n'y a pas au moins une naissance (au cours */
/* de l'intervalle de temps courant) qui la compense ; ainsi, il ne peut pas y avoir plus */
/* d'immobilisations (des morts suivant les options) que de naissances, mais il peut y */
/* avoir plus de naissances que d'immobilisations puisque les naissances sont controlees */
/* de facon inflexible par 'fichier_LISTE_DATE_DE_NAISSANCE'. Il conviendra donc, si cette */
/* option est activee, d'avoir des probabilites d'immobilisation suffisamment elevees */
/* de facon a ce qu'en l'absence de cette synchronisation, il y ait plus d'immobilisations */
/* que de naissances ; ainsi lorsque la synchronisation sera activee, on sera sur que les */
/* immobilisations seront "seuillees" par rapport aux naissances... */
EGAL(immobiliser_le_corps_courant_suite_a_un_test_probabiliste,VRAI);
/* Si le corps courant est reflechi ci-apres, il doit etre immobilise... */
Eblock
ATes
Bblock
Eblock
ETes
Eblock
ATes
Bblock
Eblock
ETes
Eblock
ATes
Bblock
/* Cas ou l'immobilisation eventuelle a lieu sur decomptage... */
Eblock
ETes
TRANSFERT_ACCROISSEMENT_3D(vecteur_directeur_OX2,Mgradient_local_tri_dimensionnel);
NORMALISATION_ACCROISSEMENT_3D(vecteur_directeur_OX2);
PRODUIT_VECTORIEL_ACCROISSEMENT_3D(vecteur_directeur_OY2
,ACCES_VITESSE_COURANTE(corps)
,Mgradient_local_tri_dimensionnel
);
/* ATTENTION, on utilise 'ACCES_VITESSE_COURANTE(corps)' et non pas 'vecteur_directeur_OZ2' */
/* de meme que 'Mgradient_local_tri_dimensionnel' et non pas 'vecteur_directeur_OX2' car, */
/* effet, on a besoin de 'module_du_vecteur_directeur_OY2' valant le module exact du produit */
/* vectoriel de 'V' et de 'G'... */
EGAL(module_du_vecteur_directeur_OY2
,longF3D(vecteur_directeur_OY2)
);
Test(IZGT(module_du_vecteur_directeur_OY2))
Bblock
/* Cas ou l'axe 'OY2' a pu etre defini... */
NORMALISATION_ACCROISSEMENT_3D(vecteur_directeur_OY2);
PRODUIT_VECTORIEL_ACCROISSEMENT_3D(vecteur_directeur_OZ2
,vecteur_directeur_OX2
,vecteur_directeur_OY2
);
NORMALISATION_ACCROISSEMENT_3D(vecteur_directeur_OZ2);
/* Definition du nouveau referentiel {OX2,OY2,OZ2} ou 'OX2' est dirige par le Gradient : */
/* */
/* --> --> */
/* I = G */
/* */
/* --> --> --> */
/* J = V /\ I */
/* */
/* --> --> --> */
/* K = I /\ J */
/* */
/* les trois vecteurs {I,J,K} etant evidemment normalises. On notera au passage que les */
/* vecteurs vitesses (incident et reflechi ou refracte) sont donc contenus dans le plan */
/* {OX2,OZ2} ce qui sera exploite par la suite lors d'un passage en coordonnees polaires */
/* et inversement... */
/* */
/* On appelera respectivement "Normale" et "Tangentielle" les composantes de la vitesse */
/* le long des axes 'OX2' et 'OZ2'. */
INITIALISATION_MATRICE_3D(matrice_de_rotation_directe
,ASD1(vecteur_directeur_OX2,dx)
,ASD1(vecteur_directeur_OX2,dy)
,ASD1(vecteur_directeur_OX2,dz)
,ASD1(vecteur_directeur_OY2,dx)
,ASD1(vecteur_directeur_OY2,dy)
,ASD1(vecteur_directeur_OY2,dz)
,ASD1(vecteur_directeur_OZ2,dx)
,ASD1(vecteur_directeur_OZ2,dy)
,ASD1(vecteur_directeur_OZ2,dz)
);
INITIALISATION_MATRICE_3D(matrice_de_rotation_inverse
,ASD1(vecteur_directeur_OX2,dx)
,ASD1(vecteur_directeur_OY2,dx)
,ASD1(vecteur_directeur_OZ2,dx)
,ASD1(vecteur_directeur_OX2,dy)
,ASD1(vecteur_directeur_OY2,dy)
,ASD1(vecteur_directeur_OZ2,dy)
,ASD1(vecteur_directeur_OX2,dz)
,ASD1(vecteur_directeur_OY2,dz)
,ASD1(vecteur_directeur_OZ2,dz)
);
/* Calcul de la matrice de rotation permettant de passer du referentiel {OX1,OY1,OZ1} au */
/* referentiel {OX2,OY2,OZ2} et de son inverse. */
PRODUIT_MATRICE_ACCROISSEMENT_3D(vecteur_vitesse_incident
,matrice_de_rotation_directe
,ACCES_VITESSE_COURANTE(corps)
);
/* Calcul de la vitesse incidente dans le referentiel {OX2,OY2,OZ2}. */
EGAL(cosinus_du_theta_incident_VG
,DIVI(prdF3D(ACCES_VITESSE_COURANTE(corps),Mgradient_local_tri_dimensionnel)
,MUL2(module_de_la_vitesse_incidente
,module_du_Mgradient_local_tri_dimensionnel
)
)
);
EGAL(sinus_du_theta_incident_VG
,DIVI(module_du_vecteur_directeur_OY2
,MUL2(module_de_la_vitesse_incidente
,module_du_Mgradient_local_tri_dimensionnel
)
)
);
EGAL(theta_incident_VG
,ATAN(sinus_du_theta_incident_VG
,cosinus_du_theta_incident_VG
)
);
/* Angle theta entre le vecteur vitesse incident et le gradient du milieu de propagation. */
EGAL(theta_apres_reflexion_ou_refraction,theta_incident_VG);
EGAL(rapport_de_l_indice_refracte_a_l_indice_incident,FU);
/* On fait comme si il devait n'y avoir ni reflexion, ni refraction... */
Test(IFINff(theta_incident_VG
,angle_inferieur_du_cone_de_reflexion
,angle_superieur_du_cone_de_reflexion
)
)
/* Dans ce cas le vecteur vitesse incident et le vecteur gradient ne vont pas dans la */
/* meme direction ; on decrete alors qu'il y a reflexion totale et donc l'index de */
/* refraction N(incident)/N(refracte) est superieur a 1. Il y a donc reflexion lorsque */
/* l'on passe d'un milieu de fort indice ('BLANC' par exemple) a un milieu de faible */
/* indice ('NOIR' par exemple). */
Bblock
Test(IFET(EST_VRAI(il_peut_y_avoir_reflexion)
,TOUJOURS_VRAI
)
)
Bblock
/* Le 'TOUJOURS_VRAI' est la uniquement pour assurer la symetrie de ce 'Test(...)' avec */
/* celui relatif a 'il_peut_y_avoir_refraction' qui va suivre... */
EGAL(theta_apres_reflexion_ou_refraction,SOUS(PI,theta_incident_VG));
/* Angle theta entre le vecteur vitesse reflechi et le gradient du milieu de propagation. */
/* ATTENTION, on notera que les lois de l'optique donnent : */
/* */
/* G ^ */
/* * | --* */
/* * +theta | -theta * | */
/* * | * */
/* * | * */
/* * | * */
/* * | * */
/* --------------------O------------------> */
/* | */
/* | */
/* | */
/* */
/* soit : */
/* */
/* theta(reflechi) = -theta(incident), */
/* */
/* or ici, la definition de 'theta' est differente. On considere : */
/* */
/* G ^ */
/* * | --* */
/* * |pi-theta * | */
/* * |-------* */
/* * |\ * */
/* * | \ * */
/* N(incident) * | *\ ('BLANC' par exemple) */
/* --------------------O---\--------------> */
/* N(refracte) | . \ ('NOIR' par exemple) */
/* | . \theta */
/* | .\ */
/* | . */
/* | . */
/* | . V */
/* */
/* soit : */
/* */
/* theta(reflechi) = pi-theta(incident), */
/* */
EGAL(rapport_de_l_indice_refracte_a_l_indice_incident,FU);
/* Le module de 'V' ne change pas lorsqu'il y a reflexion. Cette mise a 1 a deja eu lieu */
/* lors de la definition de 'rapport_de_l_indice_refracte_a_l_indice_incident', mais on */
/* le refait ici par symetrie avec le cas de la refraction qui va suivre... */
EGAL(ACCES_REFLEXIONS_COURANTS(corps),VRAI);
/* Et ce afin d'eviter que, par exemple, juste apres une reflexion, on considere lors du */
/* calcul du gradient qu'il y a refraction ; ceci est tout a fait possible puisqu'apres */
/* une reflexion la situation est inversee et que l'on risque donc de trouver un gradient */
/* inverse de celui qui a cause la reflexion et qui donc implique une refraction... */
INCR(ACCES_COMPTEURS_REFLEXIONS_COURANTS(corps),I);
/* Comptage des reflexions pour 'corps'. */
INCR(compteur_des_collisions_particules_milieu,I);
/* Comptage des collisions de type 'particules-milieu'. */
Test(IFOU(IZEQ(ACCES_IMMOBILISABLES(corps))
,IFET(IZLT(ACCES_IMMOBILISABLES(corps))
,IL_FAUT(immobiliser_le_corps_courant_suite_a_un_test_probabiliste)
)
)
)
/* Ce dispositif a ete introduit le 20010906164754 et passe de 'EST_VRAI(...)' a 'IZLE(...)' */
/* le 20010910092922 puisqu'a cette date il a ete transforme d'un indicateur logique en un */
/* decompteur. Le 20020219142536 a ete introduit la possibilite supplementaire de faire un */
/* test probabiliste, d'ou le passage a 'IFOU(IZEQ(...),...)'. */
Bblock
Test(EST_VRAI(faire_mourir_les_particules_immobilisables))
Bblock
EGAL(ACCES_DATES_DE_MORT(corps),temps_courant);
/* Ce dispositif a ete introduit le 20011008091649 ; il permet de choisir entre la mort */
/* ('VRAI') auquel cas la particule disparait, et l'immobilisation "inertielle" ('FAUX') */
/* auquel cas la particule est toujours visible mais ne peut plus bouger... */
Eblock
ATes
Bblock
INITIALISATION_ACCROISSEMENT_3D(vecteur_vitesse_incident
,FZERO
,FZERO
,FZERO
);
/* On immobilise le corps en lui donnant une vitesse nulle... */
EGAL(ACCES_MASSES(corps),MASSE_D_UN_CORPS_APRES_IMMOBILISATION);
/* Mais il faut de plus lui donner une masse infinie afin que lors de chocs ulterieurs */
/* avec d'autres particules (non immobilisees), celle-ci ne soit pas remise en mouvement */
/* (via la conservation de la quantite de mouvement...). */
Eblock
ETes
INCR(compteur_des_particules_immobilisees_lors_d_une_collision_avec_le_milieu,I);
/* Comptage... */
EDITION_E(BLOC(CAL2(Prin0(" IMMOBILISATION"));
)
);
Eblock
ATes
Bblock
Test(IZGT(ACCES_IMMOBILISABLES(corps)))
Bblock
DECR(ACCES_IMMOBILISABLES(corps),I);
/* Decomptage des reflexions pour 'corps' avant une eventuelle immobilisation. */
Eblock
ATes
Bblock
Eblock
ETes
EDITION_E(BLOC(CAL2(Prin2(" REFLEXION thetaI=%+f thetaR=%+f"
,theta_incident_VG
,theta_apres_reflexion_ou_refraction
)
);
)
);
Eblock
ETes
Eblock
ATes
Bblock
Eblock
ETes
EGAL(ACCES_REFRACTIONS_COURANTS(corps),FAUX);
/* Et ce afin d'eviter que, par exemple, juste apres une reflexion, on considere lors du */
/* calcul du gradient qu'il y a refraction ; ceci est tout a fait possible puisqu'apres */
/* une reflexion la situation est inversee et que l'on risque donc de trouver un gradient */
/* inverse de celui qui a cause la reflexion et qui donc implique une refraction... */
Eblock
ATes
/* Dans ce cas le vecteur vitesse incident et le vecteur gradient vont dans la meme */
/* direction ; on decrete qu'il y a refraction et donc l'index de refraction */
/* N(incident)/N(refracte) est inferieur a 1. Il y a donc refraction lorsque */
/* l'on passe d'un milieu de faible indice ('NOIR' par exemple) a un milieu de fort */
/* indice ('BLANC' par exemple). */
Bblock
Test(I3ET(EST_VRAI(il_peut_y_avoir_refraction)
,EST_FAUX(ACCES_REFLEXIONS_COURANTS(corps))
,IFNE(ACCES_NIVEAUX_LOCAUX_COURANTS(corps)
,niveau_central_du_pave_du_Mgradient_local_tri_dimensionnel
)
)
)
/* Lorsqu'il n'y a pas variation du niveau central de calcul du gradient, on decide */
/* arbitrairement qu'il ne peut y avoir de refraction. Cette methode introduite le */
/* 19971230143719 permet de lutter contre des refractions parasites que l'on decretait */
/* jusqu'a presente "en regardant" derriere les particules ; des inhomogeneites dans le */
/* milieu de propagation provoquait ce probleme. Par exemple : */
/* */
/* soit le milieu de propagation suivant (ou '000' et '255' representent respectivement */
/* le 'NOIR' et le 'BLANC') : */
/* */
/* 255 255 255 255 000 000 000 000 000 000 */
/* */
/* 255 255 255 255 000 000 000 000 000 000 */
/* */
/* 255 255 255 255 000 000 000 000 000 000 */
/* ------------------- */
/* 255 255 255|255 255 255 255 255|255 255 */
/* |+++ */
/* 255 255 255|255 255 255 255 255|255 255 */
/* | */
/* 255 255 255|255 255 255 255 255|255 255 */
/* | === */
/* 255 255 255|255 255 255 255 255|255 255 */
/* | */
/* 255 255 255|255 255 255 255 255|255 255 */
/* ------------------- */
/* 255 255 255 255 255 255 255 255 255 255 */
/* */
/* La particule est au point souligne par '===' a l'instant 't' et au point souligne par */
/* '+++' a l'instant 't+dt'. Dans les conditions de la sequence : */
/* */
/* xivPdf 11 1 / 028254_028765 */
/* */
/* il s'agit du corps '26' a la periode '41'. Le gradient du milieu de propagation y est */
/* alors calcule dans un pave 5x5x5 (materialise par un cadre carre ci-dessus). A */
/* l'instant 't' il ne voit donc pas les points de niveau '000' ; le gradient est alors */
/* d'ailleurs nul et il ne se passe rien. A l'instant 't+dt', il voit les points de niveau */
/* '000', mais malheureusement majoritairement derriere la particule (son vecteur vitesse */
/* va du point '===' au point '+++') ; par definition, cela correspond a de la refraction. */
/* Il est evident que celle-ci est evidemment parasite... */
Bblock
EGAL(rapport_de_l_indice_refracte_a_l_indice_incident
,COND(IFGE(module_du_Mgradient_local_tri_dimensionnel,FU)
,NEUT(module_du_Mgradient_local_tri_dimensionnel)
,INVE(module_du_Mgradient_local_tri_dimensionnel)
)
);
EGAL(theta_apres_reflexion_ou_refraction
,ASIX(DIVI(SINX(theta_incident_VG)
,rapport_de_l_indice_refracte_a_l_indice_incident
)
)
);
/* Modification de l'angle 'theta' suivant les lois de l'optique : */
/* */
/* sinus(incident) N(refracte) */
/* ----------------- = ------------- > 1 */
/* sinus(refracte) N(incident) */
/* */
/* (Loi de Snell). */
/* */
/* L'indice de refraction utilise est soit le module du gradient, soit son inverse, en */
/* choisissant celui qui est superieur a 1. Les lois de l'optique donnent : */
/* */
/* On notera qu'un passage "brutal" du '$NOIR' au '$BLANC', avec les parametres par defaut, */
/* correspond a un rapport 'rapport_de_l_indice_refracte_a_l_indice_incident' de : */
/* */
/* N(refracte) */
/* ------------- = 63.875000 */
/* N(incident) */
/* */
/* ce qui est enorme... */
Test(IFINff(theta_incident_VG,NEGA(PI_SUR_2),NEUT(PI_SUR_2)))
Bblock
/* Lorsque 'theta' est dans [-pi/2,+pi/2], le resultat de 'ASIX(...)' qui est donne dans */
/* [-pi/2,+pi/2] est alors conserve... */
Eblock
ATes
Bblock
EGAL(theta_apres_reflexion_ou_refraction
,SOUS(PI,theta_apres_reflexion_ou_refraction)
);
/* Lorsque 'theta' n'est pas dans [-pi/2,+pi/2], le resultat de 'ASIX(...)' qui est donne */
/* [-pi/2,+pi/2] est corrige, sachant que : */
/* */
/* arcsinus(theta) = arcsinus(pi-theta). */
/* */
Eblock
ETes
EGAL(ACCES_REFRACTIONS_COURANTS(corps),VRAI);
/* Et ce afin d'eviter que, par exemple, juste apres une reflexion, on considere lors du */
/* calcul du gradient qu'il y a refraction ; ceci est tout a fait possible puisqu'apres */
/* une reflexion la situation est inversee et que l'on risque donc de trouver un gradient */
/* inverse de celui qui a cause la reflexion et qui donc implique une refraction... */
INCR(ACCES_COMPTEURS_REFRACTIONS_COURANTS(corps),I);
/* Comptage des refractions pour 'corps'. */
EDITION_E(BLOC(CAL2(Prin3(" REFRACTION thetaI=%+f thetaR=%+f indice(R/I)=%+f"
,theta_incident_VG
,theta_apres_reflexion_ou_refraction
,rapport_de_l_indice_refracte_a_l_indice_incident
)
);
)
);
Test(IL_NE_FAUT_PAS(moduler_la_vitesse_lors_d_une_refraction))
Bblock
EGAL(rapport_de_l_indice_refracte_a_l_indice_incident,FU);
/* On annule ainsi le calcul de 'rapport_de_l_indice_refracte_a_l_indice_incident' dont la */
/* valeur a ete utile, et ce afin de ne pas modifier la vitesse de la particule, ce qui */
/* dans le cas de rapports trop grands provoque des pseudo-immobilisations... */
Eblock
ATes
Bblock
Eblock
ETes
Eblock
ATes
Bblock
Eblock
ETes
EGAL(ACCES_REFLEXIONS_COURANTS(corps),FAUX);
/* Et ce afin d'eviter que, par exemple, juste apres une reflexion, on considere lors du */
/* calcul du gradient qu'il y a refraction ; ceci est tout a fait possible puisqu'apres */
/* une reflexion la situation est inversee et que l'on risque donc de trouver un gradient */
/* inverse de celui qui a cause la reflexion et qui donc implique une refraction... */
Eblock
ETes
Test(IFNE_a_peu_pres_absolu(ASD1(vecteur_vitesse_incident,dy),FZERO,GRAND_EPSILON))
Bblock
PRINT_ERREUR("le vecteur vitesse n'est pas dans le plan [OX2,OZ2]");
CAL1(Prer3("gradient local = {%+f,%+f,%+f}\n"
,ASD1(Mgradient_local_tri_dimensionnel,dx)
,ASD1(Mgradient_local_tri_dimensionnel,dy)
,ASD1(Mgradient_local_tri_dimensionnel,dz)
)
);
CAL1(Prer3("vecteur dans [OX1,OY1,OZ1] = {%+f,%+f,%+f}\n"
,ASD1(ACCES_VITESSE_COURANTE(corps),dx)
,ASD1(ACCES_VITESSE_COURANTE(corps),dy)
,ASD1(ACCES_VITESSE_COURANTE(corps),dz)
)
);
CAL1(Prer3("matrice directe =\n {%+f,%+f,%+f}\n"
,ASD2(matrice_de_rotation_directe,cx,cx)
,ASD2(matrice_de_rotation_directe,cx,cy)
,ASD2(matrice_de_rotation_directe,cx,cz)
)
);
CAL1(Prer3("{%+f,%+f,%+f}\n"
,ASD2(matrice_de_rotation_directe,cy,cx)
,ASD2(matrice_de_rotation_directe,cy,cy)
,ASD2(matrice_de_rotation_directe,cy,cz)
)
);
CAL1(Prer3("{%+f,%+f,%+f}\n"
,ASD2(matrice_de_rotation_directe,cz,cx)
,ASD2(matrice_de_rotation_directe,cz,cy)
,ASD2(matrice_de_rotation_directe,cz,cz)
)
);
CAL1(Prer3("vecteur dans [OX2,OY2,OZ2] = {%+f,%+f,%+f}\n"
,ASD1(vecteur_vitesse_incident,dx)
,ASD1(vecteur_vitesse_incident,dy)
,ASD1(vecteur_vitesse_incident,dz)
)
);
CAL1(Prer3("matrice inverse =\n {%+f,%+f,%+f}\n"
,ASD2(matrice_de_rotation_inverse,cx,cx)
,ASD2(matrice_de_rotation_inverse,cx,cy)
,ASD2(matrice_de_rotation_inverse,cx,cz)
)
);
CAL1(Prer3("{%+f,%+f,%+f}\n"
,ASD2(matrice_de_rotation_inverse,cy,cx)
,ASD2(matrice_de_rotation_inverse,cy,cy)
,ASD2(matrice_de_rotation_inverse,cy,cz)
)
);
CAL1(Prer3("{%+f,%+f,%+f}\n"
,ASD2(matrice_de_rotation_inverse,cz,cx)
,ASD2(matrice_de_rotation_inverse,cz,cy)
,ASD2(matrice_de_rotation_inverse,cz,cz)
)
);
Eblock
ATes
Bblock
Eblock
ETes
EGAL(Rho_du_vecteur_vitesse_incident
,Rho_2D(ASD1(vecteur_vitesse_incident,dx)
,ASD1(vecteur_vitesse_incident,dz)
)
);
EGAL(Theta_du_vecteur_vitesse_incident
,Theta_2D(ASD1(vecteur_vitesse_incident,dx)
,ASD1(vecteur_vitesse_incident,dz)
)
);
/* Le vecteur vitesse incident, de meme que le vecteur vitesse reflechi ou refracte, sont */
/* dans le plan {OX2,OZ2} par definition. On peut donc y passer en coordonnees polaires. */
EGAL(Rho_du_vecteur_vitesse_incident
,DIVI(Rho_du_vecteur_vitesse_incident
,rapport_de_l_indice_refracte_a_l_indice_incident
)
);
/* Prise en compte de l'eventuel modification du module de la vitesse 'V' (dans le cas de */
/* la refraction). */
INITIALISATION_ACCROISSEMENT_3D(vecteur_vitesse_reflechi_ou_refracte
,Xcartesienne_2D(Rho_du_vecteur_vitesse_incident
,ADD2(Theta_du_vecteur_vitesse_incident
,SOUS(theta_apres_reflexion_ou_refraction
,theta_incident_VG
)
)
)
,ASD1(vecteur_vitesse_incident,dy)
,Ycartesienne_2D(Rho_du_vecteur_vitesse_incident
,ADD2(Theta_du_vecteur_vitesse_incident
,SOUS(theta_apres_reflexion_ou_refraction
,theta_incident_VG
)
)
)
);
/* Reflexion ou refraction du vecteur vitesse incident dans le referentiel {OX2,OY2,OZ2}. */
INITIALISATION_ACCROISSEMENT_3D(vecteur_vitesse_reflechi_ou_refracte
,AXPB(MUL2(facteur_vitesse_OX2_refraction_reflexion
,ACCES_FACTEUR_VITESSE_OX2_REFRACTION_REFLEXION(corps)
)
,ASD1(vecteur_vitesse_reflechi_ou_refracte,dx)
,ADD2(translation_vitesse_OX2_refraction_reflexion
,ACCES_TRANSLATION_VITESSE_OX2_REFRACTION_REFLEXION(corps)
)
)
,NEUT(ASD1(vecteur_vitesse_reflechi_ou_refracte,dy))
,AXPB(MUL2(facteur_vitesse_OZ2_refraction_reflexion
,ACCES_FACTEUR_VITESSE_OZ2_REFRACTION_REFLEXION(corps)
)
,ASD1(vecteur_vitesse_reflechi_ou_refracte,dz)
,ADD2(translation_vitesse_OZ2_refraction_reflexion
,ACCES_TRANSLATION_VITESSE_OZ2_REFRACTION_REFLEXION(corps)
)
)
);
/* Simulation eventuelle d'echange d'energie avec le milieu en agissant sur les composantes */
/* "Normale" ('OX2') et "Tangentielle" ('OZ2') du vecteur vitesse. */
PRODUIT_MATRICE_ACCROISSEMENT_3D(ACCES_VITESSE_COURANTE(corps)
,matrice_de_rotation_inverse
,vecteur_vitesse_reflechi_ou_refracte
);
/* La vitesse du corps courant est donc changee en fonction, localement, de la reflexion */
/* ou de la refraction due au milieu de propagation... */
/* */
/* On notera que j'ai tente la chose suivante : */
/* */
/* --> --> */
/* --> V /\ G */
/* U = -------------- */
/* --> --> */
/* | V /\ G | */
/* */
/* donne, une fois normalise, un vecteur unitaire orthogonal au plan forme par le vecteur */
/* vitesse 'V' et le gradient local 'G' du milieu de propagation. Appelons 'R' le vecteur */
/* vitesse Reflechi ou Refracte. On a evidemment : */
/* */
/* --> --> --> --> --> */
/* R /\ G = | R |.| G |.sin(theta). U */
/* */
/* ou theta est l'angle entre 'R' et 'G' qui est donne par les lois de l'optique (Reflexion */
/* ou Refraction). De la, on peut tirer le vecteur 'R' : */
/* */
/* --> --> --> --> -1 --> */
/* R = [ | R |.| G |.sin(theta). U ] /\ G */
/* */
/* a l'aide de l'anti-produit vectoriel ('v $ximd/operator.1$FON APvect'). Malheureusement, */
/* j'ai du y renoncer le 19971204091928 car il y a indetermination dans le calcul de */
/* l'anti-produit vectoriel comme cela est explique dans 'v $ximd/operator.1$FON APvect'. */
EGAL(module_de_la_vitesse_apres_reflexion_ou_refraction
,longF3D(ACCES_VITESSE_COURANTE(corps))
);
Test(IFET(IZNE(module_de_la_vitesse_apres_reflexion_ou_refraction)
,IFLE(module_de_la_vitesse_apres_reflexion_ou_refraction,pEPSILON)
)
)
Bblock
PRINT_ERREUR("un vecteur vitesse est trop petit, sans etre nul");
/* Cela peut se produire en particulier depuis l'introduction de l'immobilisation des */
/* particules lors de leurs collisions avec le milieu lorsque la condition */
/* 'IL_NE_FAUT_PAS(faire_mourir_les_particules_immobilisables)' ; en effet, leur masse */
/* devient quasiment infinie, ce qui donne des vitesses infinitesimales lors des chocs */
/* (voir donc 'fichier_LISTE_IMMOBILISABLE' introduit le 20010906164754). */
CAL1(Prer1("il s'agit du corps %d\n"
,corps
)
);
CAL1(Prer3("vitesse = {%+g,%+g,%+g}\n"
,ASD1(ACCES_VITESSE_COURANTE(corps),dx)
,ASD1(ACCES_VITESSE_COURANTE(corps),dy)
,ASD1(ACCES_VITESSE_COURANTE(corps),dz)
)
);
INITIALISATION_ACCROISSEMENT_3D(ACCES_VITESSE_COURANTE(corps)
,FZERO
,FZERO
,FZERO
);
/* Lorsque le module de la vitesse est trop petit, cela signifie qu'en general il y a eu */
/* trop de refraction (avec un indice trop grand). La vitesse a donc tendance a diminuer */
/* tres vite. Pour eviter des problemes ulterieurs avec les nombres flottants ("underflow"s) */
/* on annule la vitesse... */
Eblock
ATes
Bblock
Eblock
ETes
Eblock
ATes
Bblock
/* Cas ou l'axe 'OY2' n'a pu etre defini ; cela signifie que le gradient 'G' et la vitesse */
/* 'V' sont colineaires. On ne peut donc pas definir de referentiel ; on fait donc comme */
/* si le milieu etait parfaitement reflechissant par rapport a un plan orthogonal a 'G'. */
EGAL(ACCES_REFLEXIONS_COURANTS(corps),VRAI);
/* Et ce afin d'eviter que, par exemple, juste apres une reflexion, on considere lors du */
/* calcul du gradient qu'il y a refraction ; ceci est tout a fait possible puisqu'apres */
/* une reflexion la situation est inversee et que l'on risque donc de trouver un gradient */
/* inverse de celui qui a cause la reflexion et qui donc implique une refraction... */
INCR(ACCES_COMPTEURS_REFLEXIONS_COURANTS(corps),I);
/* Comptage des reflexions pour 'corps'. */
INCR(compteur_des_collisions_particules_milieu,I);
/* Comptage des collisions de type 'particules-milieu'. */
Test(IFOU(IZEQ(ACCES_IMMOBILISABLES(corps))
,IFET(IZLT(ACCES_IMMOBILISABLES(corps))
,IL_FAUT(immobiliser_le_corps_courant_suite_a_un_test_probabiliste)
)
)
)
/* Ce dispositif a ete introduit le 20010906164754 et passe de 'EST_VRAI(...)' a 'IZLE(...)' */
/* le 20010910092922 puisqu'a cette date il a ete transforme d'un indicateur logique en un */
/* decompteur. Le 20020219142536 a ete introduit la possibilite supplementaire de faire un */
/* test probabiliste, d'ou le passage a 'IFOU(IZEQ(...),...)'. */
Bblock
Test(EST_VRAI(faire_mourir_les_particules_immobilisables))
Bblock
EGAL(ACCES_DATES_DE_MORT(corps),temps_courant);
/* Ce dispositif a ete introduit le 20011008091649 ; il permet de choisir entre la mort */
/* ('VRAI') auquel cas la particule disparait, et l'immobilisation "inertielle" ('FAUX') */
/* auquel cas la particule est toujours visible mais ne peut plus bouger... */
Eblock
ATes
Bblock
INITIALISATION_ACCROISSEMENT_3D(ACCES_VITESSE_COURANTE(corps)
,FZERO
,FZERO
,FZERO
);
/* On immobilise le corps en lui donnant une vitesse nulle... */
EGAL(ACCES_MASSES(corps),MASSE_D_UN_CORPS_APRES_IMMOBILISATION);
/* Mais il faut de plus lui donner une masse infinie afin que lors de chocs ulterieurs */
/* avec d'autres particules (non immobilisees), celle-ci ne soit pas remise en mouvement */
/* (via la conservation de la quantite de mouvement...). */
Eblock
ETes
INCR(compteur_des_particules_immobilisees_lors_d_une_collision_avec_le_milieu,I);
/* Comptage... */
EDITION_E(BLOC(CAL2(Prin0(" IMMOBILISATION"));
)
);
Eblock
ATes
Bblock
Test(IZGT(ACCES_IMMOBILISABLES(corps)))
Bblock
DECR(ACCES_IMMOBILISABLES(corps),I);
/* Decomptage des reflexions pour 'corps' avant une eventuelle immobilisation. */
Eblock
ATes
Bblock
Eblock
ETes
INITIALISATION_ACCROISSEMENT_3D(ACCES_VITESSE_COURANTE(corps)
,NEGA(ASD1(ACCES_VITESSE_COURANTE(corps),dx))
,NEGA(ASD1(ACCES_VITESSE_COURANTE(corps),dy))
,NEGA(ASD1(ACCES_VITESSE_COURANTE(corps),dz))
);
/* On repart donc dans la direction inverse de la direction incidente... */
EDITION_E(BLOC(CAL2(Prin0(" REFLEXION_ABSOLUE"));
)
);
Eblock
ETes
Eblock
ETes
Eblock
ATes
Bblock
Test(IFINff(NIVEAU_LOCAL_COURANT_INITIAL
,______________NOIR_NORMALISE
,______________BLANC_NORMALISE
)
)
Bblock
PRINT_ERREUR("la valeur de 'NIVEAU_LOCAL_COURANT_INITIAL' est incompatible avec les tests");
Eblock
ATes
Bblock
Eblock
ETes
Test(I3ET(IFGT(temps_courant,instant_initial)
,IFNE(ACCES_NIVEAUX_LOCAUX_COURANTS(corps)
,NIVEAU_LOCAL_COURANT_INITIAL
)
,IFNE(niveau_central_du_pave_du_Mgradient_local_tri_dimensionnel
,ACCES_NIVEAUX_LOCAUX_COURANTS(corps)
)
)
)
/* Le test sur 'NIVEAU_LOCAL_COURANT_INITIAL' a ete introduit le 19990927173502 parce que */
/* le test sur 'instant_initial' ne suffit pas ; en effet, lorsque des particules sont */
/* initialement au repos et qu'elles ne sont mises que plus tard en mouvement, la valeur */
/* de 'ACCES_NIVEAUX_LOCAUX_COURANTS(corps)' n'est donc pas correcte dans les instants */
/* suivants. Cela s'est vu en generant la sequence : */
/* */
/* xivPdf 9 2 / 015644_016155 */
/* */
/* On notera que dans ces conditions, le test : */
/* */
/* IFGT(temps_courant,instant_initial) */
/* */
/* n'est plus utile ; je le garde malgre tout, car on ne sait jamais... */
Bblock
Test(IL_FAUT(editer_les_messages_d_erreur_du_calcul_du_gradient))
/* Test introduit le 20150208081059... */
Bblock
PRINT_ERREUR("le Gradient est nul");
CAL1(Prer3("au point........ = {%d,%d,%d}\n",Xg,Yg,Zg));
CAL1(Prer1("alors que le niveau du milieu pour le corps %d a change\n",corps));
CAL1(Prer0("(cela peut se produire si le milieu est dynamique et si son 'volume' diminue)\n"));
/* Effectivement, lors d'une reduction du "volume" du milieu, une particule peut voir */
/* changer brutalement son milieu local, sans qu'elle se soit elle-meme deplacee... */
CAL1(Prer0("(ou encore en presence d'un champ de gravitation trop fort)\n"));
CAL1(Prer0("(auquel cas des accelerations violentes peuvent provoquer)\n"));
CAL1(Prer0("(des deplacements trop importants par rapport au pas de temps)\n"));
/* Effectivement, en presence d'un champ de gravitation trop fort, les accelerations sont */
/* alors tres violentes et les vitesses peuvent donc devenir tres importantes ; certaines */
/* particules peuvent donc changer de "zones" dans le milieu sans que cela se voit (a cause */
/* d'un pas de temps trop grand et qui le sera toujours puisqu'alors les vitesses augmentent */
/* sans arret sauf collision...). */
CAL1(Prer1("niveau anterieur.......... = %g\n",ACCES_NIVEAUX_LOCAUX_COURANTS(corps)));
CAL1(Prer1("niveau courant............ = %g\n"
,niveau_central_du_pave_du_Mgradient_local_tri_dimensionnel
)
);
CAL1(Prer1("periode................... = %d\n"
,numero_de_la_periode_courante_de_la_simulation
)
);
CAL1(Prer1("temps..................... = %f\n"
,temps_courant
)
);
CAL1(Prer3("coordonnees gradient...... = {%+d,%+d,%+d}\n"
,Xg
,Yg
,Zg
)
);
CAL1(Prer6("coordonnees particule..... = {%+f,%+f,%+f} --> {%+d,%+d,%+d}\n"
,ASD1(ACCES_COORDONNEES_COURANTES(corps),x)
,ASD1(ACCES_COORDONNEES_COURANTES(corps),y)
,ASD1(ACCES_COORDONNEES_COURANTES(corps),z)
,INTE(gX_PHYSIQUE_A_VISUALISATION(ASD1(ACCES_COORDONNEES_COURANTES(corps),x),FZERO))
,INTE(gY_PHYSIQUE_A_VISUALISATION(ASD1(ACCES_COORDONNEES_COURANTES(corps),y),FZERO))
,INTE(gZ_PHYSIQUE_A_VISUALISATION(ASD1(ACCES_COORDONNEES_COURANTES(corps),z),FZERO))
)
);
CAL1(Prer3("vitesse particule......... = {%+f,%+f,%+f}\n"
,ASD1(ACCES_VITESSE_COURANTE(corps),dx)
,ASD1(ACCES_VITESSE_COURANTE(corps),dy)
,ASD1(ACCES_VITESSE_COURANTE(corps),dz)
)
);
Eblock
ATes
Bblock
Eblock
ETes
Eblock
ATes
Bblock
Eblock
ETes
EGAL(ACCES_REFLEXIONS_COURANTS(corps),FAUX);
EGAL(ACCES_REFRACTIONS_COURANTS(corps),FAUX);
/* Et ce afin d'eviter que, par exemple, juste apres une reflexion, on considere lors du */
/* calcul du gradient qu'il y a refraction ; ceci est tout a fait possible puisqu'apres */
/* une reflexion la situation est inversee et que l'on risque donc de trouver un gradient */
/* inverse de celui qui a cause la reflexion et qui donc implique une refraction... */
EDITION_E(BLOC(CAL2(Prin0(" RIEN"));
)
);
Eblock
ETes
EGAL(ACCES_NIVEAUX_LOCAUX_COURANTS(corps)
,niveau_central_du_pave_du_Mgradient_local_tri_dimensionnel
);
/* Pour tester a l'instant suivant s'il y a variation du niveau central de calcul du */
/* gradient. Lorsqu'il ne change pas, on decide arbitrairement qu'il ne peut y avoir de */
/* refraction... */
Eblock
ATes
Bblock
Eblock
ETes
Eblock