/*************************************************************************************************************************************/
/* */
/* E T U D E D E L ' E R G O D I C I T E D E L A D I F F U S I O N D ' U N E P A R T I C U L E : */
/* */
/* */
/* .. */
/* .: */
/* . . .. */
/* . .. .. .: */
/* .. : : . .. . : .. */
/* ... . o . + . :. : :- . . */
/* o. .oo. *- o o..-o. ..+.-. .. */
/* -o:..:-+.. -+..: :: +--. o .. .... */
/* : .*::.-:.- .++o: +. */
/* . :. -#+*-:. o */
/* .::o:. .+ oo +*-. */
/* ...-o#:. :-. :#+.. */
/* o#++ o o#: */
/* .* : o .-::. */
/* ++ .o+#*. +.. */
/* - ooo** -o. */
/* * . */
/* .o o */
/* *. .. . */
/* +...+o-* *- o-+. */
/* ::... +o-o .. - */
/* . :..: +.:: .. */
/* :+ .. :+ + */
/* + o + -.:+o-. */
/* .. :* ...- */
/* o.* .. */
/* ++: o */
/* .. . */
/* : */
/* */
/* */
/* Author of '$xrk/ergodique.11$K' : */
/* */
/* Jean-Francois COLONNA (LACTAMME, 1993??????????). */
/* */
/*************************************************************************************************************************************/
/*===================================================================================================================================*/
/*************************************************************************************************************************************/
/* */
/* I N T E R F A C E ' listG ' : */
/* */
/* */
/* :Debut_listG: */
/* :Fin_listG: */
/* */
/*************************************************************************************************************************************/
/*===================================================================================================================================*/
/*************************************************************************************************************************************/
/* */
/* D I R E C T I V E S S P E C I F I Q U E S D E C O M P I L A T I O N : */
/* */
/*************************************************************************************************************************************/
@define PRAGMA_CL_____MODULE_NON_OPTIMISABLE
/*===================================================================================================================================*/
/*************************************************************************************************************************************/
/* */
/* F I C H I E R S D ' I N C L U D E S : */
/* */
/*************************************************************************************************************************************/
#include INCLUDES_BASE
/*===================================================================================================================================*/
/*************************************************************************************************************************************/
/* */
/* D E F I N I T I O N S D E B A S E E T U N I V E R S E L L E S : */
/* */
/*************************************************************************************************************************************/
#include xrk/attractor.11.I"
/*===================================================================================================================================*/
/*************************************************************************************************************************************/
/* */
/* 3 */
/* D E F I N I T I O N D E L ' E S P A C E P H Y S I Q U E D A N S R ( D E B U T ) : */
/* */
/* */
/* Nota : */
/* */
/* Les extrema des coordonnees {x,y,z} */
/* ainsi que ceux de leurs differentielles */
/* {dx,dy,dz} sont fixees un peu arbitrairement */
/* et sans etre parametrees. */
/* */
/* */
/*************************************************************************************************************************************/
#define hXmin_ESPACE \
PARE(-20.0)
#define hYmin_ESPACE \
PARE(-20.0)
#define hZmin_ESPACE \
PARE(-20.0)
/* Definition du "coin" inferieur-gauche-arriere de l'espace physique. */
#define hXmax_ESPACE \
PARE(20.0)
#define hYmax_ESPACE \
PARE(20.0)
#define hZmax_ESPACE \
PARE(20.0)
/* Definition du "coin" superieur-droite-avant de l'espace physique. */
/*===================================================================================================================================*/
/*************************************************************************************************************************************/
/* */
/* 3 */
/* D E F I N I T I O N D E L ' E S P A C E P H Y S I Q U E D A N S R ( D E B U T ) : */
/* */
/*************************************************************************************************************************************/
#include xrk/attractor.12.I"
#define dXmin_ESPACE \
PARE(-2.0)
#define dYmin_ESPACE \
PARE(-2.0)
#define dZmin_ESPACE \
PARE(-2.0)
/* Definition des minima des differentielles {dx,dy,dz}. */
#define dXmax_ESPACE \
PARE(2.0)
#define dYmax_ESPACE \
PARE(2.0)
#define dZmax_ESPACE \
PARE(2.0)
/* Definition des maxima des differentielles {dx,dy,dz}. */
#include xrk/attractor.1D.I"
/* Formules de renormalisation des differentielles dans [0,1] ; elles sont utilisees lorsque */
/* la production d'images en couleurs est demandee (voir 'visualiser_en_RVB'). */
/*===================================================================================================================================*/
/*************************************************************************************************************************************/
/* */
/* D E F I N I T I O N D E S D I F F E R E N T S E S P A C E S E T D E L ' E F F E T D E B R U M E : */
/* */
/*************************************************************************************************************************************/
#include xrk/attractor.13.I"
/*===================================================================================================================================*/
/*************************************************************************************************************************************/
/* */
/* A I D E A U C A D R A G E D E S I M A G E S : */
/* */
/*************************************************************************************************************************************/
#include xrk/attractor.1C.I"
DONNEES_DE_RECHERCHE_DES_EXTREMA_DES_COORDONNEES_ET_DES_DERIVEES
/* Definition des extrema des coordonnees et des derivees. On notera bien l'absence de */
/* point-virgule apres 'DONNEES_DE_RECHERCHE_DES_EXTREMA_DES_COORDONNEES_ET_DES_DERIVEES'. */
/*===================================================================================================================================*/
/*************************************************************************************************************************************/
/* */
/* G E N E R A T I O N D E S I M A G E S : */
/* */
/*************************************************************************************************************************************/
#include xrv/champs_5.14.I"
/*===================================================================================================================================*/
/*************************************************************************************************************************************/
/* */
/* D E F I N I T I O N S G E N E R A L E S R E L A T I V E S A L A V I S U A L I S A T I O N : */
/* */
/*************************************************************************************************************************************/
#define DCT \
GRO1(FRA8(FU))
DEFV(Local,DEFV(Float,INIT(dct,DCT)));
/* Definition de 'dt'. On notera la valeur 0.125 utilisee (apres que 0.2 l'ait ete) ; cette */
#include xrk/attractor.14.I"
/*===================================================================================================================================*/
/*************************************************************************************************************************************/
/* */
/* F O N C T I O N D E M E M O R I S A T I O N D U P O I N T C O U R A N T : */
/* */
/*************************************************************************************************************************************/
#include xrk/attractor.16.I"
#define RAYON_DE_VISUALISATION \
GRO1(FRA10(FRA10(mhXYZlongueur_ESPACE)))
DEFV(Local,DEFV(Float,INIT(rayon_de_visualisation,RAYON_DE_VISUALISATION)));
/* Rayon du disque materialisant une iteration. Il fut exprime longtemps sous la */
/* forme : */
/* */
/* GRO4(FRA10(FU)) */
/* */
BFonctionI
DEFV(Local,DEFV(FonctionI,memorisation_1_point_07(AXf,AYf,AZf,AdXf,AdYf,AdZf,numero_de_l_iteration_courante)))
DEFV(Argument,DEFV(Float,AXf));
DEFV(Argument,DEFV(Float,AYf));
DEFV(Argument,DEFV(Float,AZf));
/* Definition de la position {x,y,z} de l'iteration courante. */
DEFV(Argument,DEFV(Float,AdXf));
DEFV(Argument,DEFV(Float,AdYf));
DEFV(Argument,DEFV(Float,AdZf));
/* Definition des differentielles {dx,dy,dz} de la position de l'iteration courante. */
DEFV(Argument,DEFV(Int,numero_de_l_iteration_courante));
/* Numero de l'iteration courante afin d'attenuer eventuellement la luminance des points */
/* materialisant chaque iteration en fonction de leur numero (les premieres iterations etant */
/* plus sombres, et les dernieres etant plus lumineuses). */
/*-----------------------------------------------------------------------------------------------------------------------------------*/
Bblock
#include xrk/attractor.15.I"
INIT_ERROR;
/*..............................................................................................................................*/
MEMORISATION_DU_POINT_COURANT(X_DERIVEE_DANS_01(AdXf)
,Y_DERIVEE_DANS_01(AdYf)
,Z_DERIVEE_DANS_01(AdZf)
);
/* Memorisation du point courant en Noir et Blanc ou en Couleurs, mais uniquement s'il est */
/* visible en fonction des conditions de visualisation... */
RETU_ERROR;
Eblock
EFonctionI
/*===================================================================================================================================*/
/*************************************************************************************************************************************/
/* */
/* D E F I N I T I O N D U S Y S T E M E D E D I F F U S I O N : */
/* */
/* */
/* Definition (Claude Bardos & Francois Golse) : */
/* */
/* Soit une particule P de coordonnees */
/* {x,y,z}. Son mouvement se fait dans l'espace */
/* tridimensionnel entre deux plans paralleles */
/* sur lesquels elle rebondit. Il est calcule */
/* iterativement connaissant les derivees en 't' */
/* des trois coordonnees : */
/* */
/* dx */
/* ---- = F (u,v) */
/* dt x */
/* */
/* dy */
/* ---- = F (u,v) */
/* dt y */
/* */
/* dz */
/* ---- = constante (en valeur absolue) */
/* dt */
/* */
/* ou : */
/* */
/* u E [0,2.pi] */
/* v E [0,2.pi] */
/* */
/* et : */
/* */
/* 'F (u,v)' et 'F (u,v)' etant deux fonctions quelconques de moyenne nulle : */
/* x y */
/* */
/* */
/* <F (u,v)> = 0 */
/* x */
/* */
/* <F (u,v)> = 0 */
/* y */
/* */
/* par exemple : */
/* */
/* */
/* F (u,v) = cos(u) */
/* x */
/* */
/* F (u,v) = cos(v) */
/* y */
/* */
/* */
/* */
/* Les donnees initiales du probleme */
/* sont : */
/* */
/* (x ,y ,z ), */
/* 0 0 0 */
/* */
/* (u ,v ). */
/* 0 0 */
/* */
/* */
/* La reflexion sur l'un des deux plans paralleles */
/* (que l'on choisit perpendiculaires a l'axe 'OZ' */
/* a cause des definitions des trois derivees, et */
/* particulierement de la constance de la derivee de */
/* la coordonnee 'z') se fait en changeant le couple */
/* (u,v) courant suivant la loi ("chat d'Arnold") : */
/* */
/* */
/* | u | | 2 1 || u | */
/* | | = | || | (modulo 2.pi) */
/* | v | | 1 1 || v | */
/* */
/* [apres] [avant] */
/* */
/* associee a : */
/* */
/* dz dz */
/* ---- = - ---- */
/* dt dt */
/* */
/* */
/*************************************************************************************************************************************/
#include xrk/attractor.17.I"
dfTRANSFORMAT_31(liste_NOMBRE_D_ITERATIONS,fichier_NOMBRE_D_ITERATIONS,NOMBRE_D_ITERATIONS_IMPLICITE,NOMBRE_D_ITERATIONS)
/* Definition du fichier des nombres d'iterations. */
#define sNOMBRE_D_ITERATIONS(numero_de_la_periode) \
INTE(sTRANSFORMAT_31(numero_de_la_periode,liste_NOMBRE_D_ITERATIONS)) \
/* Formule generale definissant les variations du nombre d'iterations. */
#define CU0 \
GRO1(FRA1(FU))
#define CV0 \
GRO1(FRA1(FU))
DEFV(Local,DEFV(Float,INIT(cu,CU0)));
DEFV(Local,DEFV(Float,INIT(cv,CV0)));
/* Definition des deux angles (u,v), */
#define COMPOSANTE_EN_Z_DE_LA_VITESSE \
FU
#define DCZ \
COMPOSANTE_EN_Z_DE_LA_VITESSE
#define CX0 \
FZERO
#define CY0 \
FZERO
#define CZ0 \
FZERO
/* Nouvelle valeur est destinee a limiter un probleme vu dans la sequence : */
/* */
/* xivPdf 2 1 / 017100_017227 */
/* */
/* pour laquelle la transition entre l'avant-derniere et la derniere images est bizarre ; en */
/* particulier une des "branches" en bas a droite change de couleur (du bleu a l'orange). */
/* Une partie du probleme vient de l'incrementation/decrementation de la coordonnee 'z' */
/* (voir la definition de la fonction 'Fz(...)') ce qui a ete mise en evidence grace au */
/* programme 'v $xtc/increment.01$c' ; ce probleme a disparu pour un pas 'dct' inverse d'une */
/* puissance de 2. Malheureusement ce n'est pas le seul probleme, et il reste surement */
/* quelque chose qui releve du meme principe. Ca y est j'ai compris : */
/* */
/* Soient par exemple les trois points suivants (dans l'espace physique) : */
/* */
/* P1=(X,Y,Z)=(7.3155693120848868,8.2766071523526143,-0.8750000000000000) */
/* P2=(X,Y,Z)=(7.4171054615643000,8.3137138067576046,-1.0000000000000000) */
/* P3=(X,Y,Z)=(7.3158618945904497,8.2742346945782863,-0.8750000000000000) */
/* */
/* . */
/* Y /|\ */
/* | */
/* | P2 */
/* | */
/* | */
/* | */
/* | P1 */
/* | */
/* | P3 */
/* | */
/* | */
/* O---------------------------------------> */
/* */
/* X */
/* */
/* Les points 'P1' et 'P3' sont donc dans le meme plan orthogonal a l'axe 'OZ' physique, */
/* ce que l'on va noter : */
/* */
/* P1=P3 */
/* */
/* */
/* Faisons maintenant les deux rotations suivantes d'angle 'ANGLE' autour de l'axe 'OX' : */
/* */
/* 1-ANGLE=pi-epsilon (=3.1) : */
/* */
/* Pm1=(X,Y,Z)=(0.7438523104028295,0.2255644649570133,0.0406130021528017) */
/* Pm2=(X,Y,Z)=(0.7472368487188100,0.2245018989596498,0.0448274959213018) */
/* Pm3=(X,Y,Z)=(0.7438620631530150,0.2256434784888424,0.0406097138739400) */
/* */
/* et le point 'Pm1' est devant le point 'Pm3', soit : */
/* */
/* Pm1>Pm3 */
/* */
/* (l'axe 'OZ' va d'arriere en avant) */
/* */
/* */
/* 2-ANGLE=pi+epsilon (=3.2) : */
/* */
/* Pp1=(X,Y,Z)=(0.7438523104028295,0.2228809647667255,0.0130122691938415) */
/* Pp2=(X,Y,Z)=(0.7472368487188100,0.2214029598611198,0.0170996284541076) */
/* Pp3=(X,Y,Z)=(0.7438620631530150,0.2229599118401223,0.0130168855335212) */
/* */
/* et le point 'Pp1' est derriere le point 'Pp3', soit : */
/* */
/* Pp1<Pp3 */
/* */
/* */
/* Ainsi, les points 'P1' et 'P3' n'ayant pas la meme couleur, il y a, lors du passage */
/* de l'angle de rotation par la valeur 'pi' un basculement des couleurs... */
/*===================================================================================================================================*/
/*************************************************************************************************************************************/
/* */
/* F O N C T I O N S D E V I S U A L I S A T I O N E T D ' I N T E R P O L A T I O N : */
/* */
/*************************************************************************************************************************************/
#include xrv/particule.31.I"
dfTRANSFORMAT_31(liste_VARIABLE_cu0,fichier_VARIABLE_cu0,VARIABLE_cu0_IMPLICITE,CU0)
dfTRANSFORMAT_31(liste_VARIABLE_cv0,fichier_VARIABLE_cv0,VARIABLE_cv0_IMPLICITE,CV0)
/* Definition des fichiers des valeurs initiales des deux angles (u0,v0). */
#define sVARIABLE_cu0(numero_de_la_periode) \
FLOT(sTRANSFORMAT_31(numero_de_la_periode,liste_VARIABLE_cu0))
#define sVARIABLE_cv0(numero_de_la_periode) \
FLOT(sTRANSFORMAT_31(numero_de_la_periode,liste_VARIABLE_cv0))
/* Formule generale definissant les variations des valeurs initiales des deux angles. */
dfTRANSFORMAT_31(liste_VARIABLE_dcz,fichier_VARIABLE_dcz,VARIABLE_dcz_IMPLICITE,DCZ)
/* Definition du fichier de la valeur initiale de la composante en 'Z' de la vitesse. */
#define sVARIABLE_dcz(numero_de_la_periode) \
FLOT(sTRANSFORMAT_31(numero_de_la_periode,liste_VARIABLE_dcz))
/* Formule generale definissant la variation de la valeur initiale de la composante en 'Z' */
/* de la vitesse. */
dfTRANSFORMAT_31(liste_VARIABLE_cx0,fichier_VARIABLE_cx0,VARIABLE_cx0_IMPLICITE,CX0)
dfTRANSFORMAT_31(liste_VARIABLE_cy0,fichier_VARIABLE_cy0,VARIABLE_cy0_IMPLICITE,CY0)
dfTRANSFORMAT_31(liste_VARIABLE_cz0,fichier_VARIABLE_cz0,VARIABLE_cz0_IMPLICITE,CZ0)
/* Definition des fichiers des valeurs initiales des trois variables (x0,y0,z0). */
#define sVARIABLE_cx0(numero_de_la_periode) \
FLOT(sTRANSFORMAT_31(numero_de_la_periode,liste_VARIABLE_cx0))
#define sVARIABLE_cy0(numero_de_la_periode) \
FLOT(sTRANSFORMAT_31(numero_de_la_periode,liste_VARIABLE_cy0))
#define sVARIABLE_cz0(numero_de_la_periode) \
FLOT(sTRANSFORMAT_31(numero_de_la_periode,liste_VARIABLE_cz0))
/* Formule generale definissant les variations des valeurs initiales des trois variables. */
dfTRANSFORMAT_31(liste_PAS_DE_TEMPS_dct,fichier_PAS_DE_TEMPS_dct,PAS_DE_TEMPS_dct_IMPLICITE,DCT)
/* Definition du fichier des pas de temps. */
#define sPAS_DE_TEMPS_dct(numero_de_la_periode) \
FLOT(sTRANSFORMAT_31(numero_de_la_periode,liste_PAS_DE_TEMPS_dct))
/* Formule generale definissant les variations du pas de temps. */
/*===================================================================================================================================*/
/*************************************************************************************************************************************/
/* */
/* M I S E E N P L A C E D ' U N S Y S T E M E " N O N D I F F U S I F " : */
/* */
/* */
/* Definition : */
/* */
/* Au lieu de prendre les fonctions 'Fx' et */
/* 'Fy' definies precedemment, on utilise les */
/* fonction 'Gx' et 'Gy' definies ainsi : */
/* */
/* dx */
/* ---- = G (u,v) */
/* dt x */
/* */
/* dy */
/* ---- = G (u,v) */
/* dt y */
/* */
/* ou : */
/* */
/* G (u,v) = F (T(u,v)) - F (u,v) */
/* x x x */
/* */
/* G (u,v) = F (T(u,v)) - F (u,v) */
/* y y y */
/* */
/* mais en fait pour simplifier les choses, nous allons utiliser : */
/* */
/* -1 */
/* G (u,v) = F (u,v) - F (T (u,v)) */
/* x x x */
/* */
/* -1 */
/* G (u,v) = F (u,v) - F (T (u,v)) */
/* y y y */
/* */
/* (ou 'T' represente la matrice d'Arnold) */
/* */
/* Dans ces conditions, les particules devraient */
/* diffuser beaucoup moins qu'avec les fonctions */
/* 'Fx' et 'Fy'... */
/* */
/* */
/*************************************************************************************************************************************/
#define UN_SYSTEME_DIFFUSIF \
VRAI
DEFV(Local,DEFV(Logical,INIT(un_systeme_diffusif,UN_SYSTEME_DIFFUSIF)));
/* Indique si l'on utilise les fonctions [Fx,Fy] ('VRAI') ou les fonctions [Gx,Gy] ('FAUX'). */
DEFV(Local,DEFV(Float,INIT(cu_avant_la_reflexion,FLOT__UNDEF)));
DEFV(Local,DEFV(Float,INIT(cv_avant_la_reflexion,FLOT__UNDEF)));
/* Memorisation des angles 'u' et 'v' avant la precedente reflexion (ces deux valeurs sont */
/* remises a 0 lors de toute reinitialisation...). */
/*===================================================================================================================================*/
/*************************************************************************************************************************************/
/* */
/* D E F I N I T I O N D E S T R O I S F O N C T I O N S ' F ' : */
/* */
/*************************************************************************************************************************************/
#define Fx(u,v) \
COSX(u) \
/* Definition de la fonction F (u,v). */ \
/* x */
#define Fy(u,v) \
COSX(v) \
/* Definition de la fonction F (u,v). */ \
/* y */
#define Fz(u,v) \
composante_en_Z_de_la_vitesse \
/* Definition de la fonction F (u,v). */ \
/* z */
DEFV(Local,DEFV(Float,INIT(composante_en_Z_de_la_vitesse,COMPOSANTE_EN_Z_DE_LA_VITESSE)));
/* Definition de la composante en 'z' de la vitesse. */
/*===================================================================================================================================*/
/*************************************************************************************************************************************/
/* */
/* D E F I N I T I O N D E S D E U X P L A N S P A R A L L E L E S : */
/* */
/*************************************************************************************************************************************/
#define Z_DU_PLAN_1 \
NEGA(FU)
DEFV(Local,DEFV(Float,INIT(z_du_plan_1,Z_DU_PLAN_1)));
/* Definition du premier orthogonal a l'axe 'OZ', */
#define Z_DU_PLAN_2 \
NEUT(FU)
DEFV(Local,DEFV(Float,INIT(z_du_plan_2,Z_DU_PLAN_2)));
/* Definition du second orthogonal a l'axe 'OZ'. */
/*===================================================================================================================================*/
/*************************************************************************************************************************************/
/* */
/* D E F I N I T I O N D E L A R E F L E X I O N S U R L E S D E U X P L A N S P A R A L L E L E S : */
/* */
/*************************************************************************************************************************************/
#define REFLEXION_11 \
FDEUX
#define REFLEXION_12 \
FU
#define REFLEXION_21 \
FU
#define REFLEXION_22 \
FU
/*===================================================================================================================================*/
/*************************************************************************************************************************************/
/* */
/* D E F I N I T I O N D E L ' I N T E G R A T I O N D U S Y S T E M E */
/* D ' E Q U A T I O N S D I F F E R E N T I E L L E S : */
/* */
/*************************************************************************************************************************************/
#include xrk/integr.1B.vv.I"
/*===================================================================================================================================*/
/*************************************************************************************************************************************/
/* */
/* D E F I N I T I O N D E S I N I T I A L I S A T I O N S : */
/* */
/*************************************************************************************************************************************/
#include xrk/attractor.18.I"
/*===================================================================================================================================*/
/*************************************************************************************************************************************/
/* */
/* E T U D E D E L ' E R G O D I C I T E D E L A D I F F U S I O N D ' U N E P A R T I C U L E : */
/* */
/*************************************************************************************************************************************/
BCommande(nombre_d_arguments,arguments)
/*-----------------------------------------------------------------------------------------------------------------------------------*/
Bblock
/*..............................................................................................................................*/
INITIALISATIONS_GENERALES;
/* Initialisations generales faites au tout debut... */
iTRANSFORMAT_31(liste_VARIABLE_cu0,VARIABLE_cu0_IMPLICITE);
iTRANSFORMAT_31(liste_VARIABLE_cv0,VARIABLE_cv0_IMPLICITE);
/* Initialisation des valeurs initiales des deux angles (u0,v0). */
iTRANSFORMAT_31(liste_VARIABLE_dcz,VARIABLE_dcz_IMPLICITE);
/* Initialisation de la valeur implicite de la valeur initiale de la composante en 'Z' de */
/* la vitesse. */
iTRANSFORMAT_31(liste_VARIABLE_cx0,VARIABLE_cx0_IMPLICITE);
iTRANSFORMAT_31(liste_VARIABLE_cy0,VARIABLE_cy0_IMPLICITE);
iTRANSFORMAT_31(liste_VARIABLE_cz0,VARIABLE_cz0_IMPLICITE);
/* Initialisation des valeurs initiales des trois variables (x0,y0,z0). */
iTRANSFORMAT_31(liste_PAS_DE_TEMPS_dct,PAS_DE_TEMPS_dct_IMPLICITE);
/* Initialisation du pas de temps. */
iTRANSFORMAT_31(liste_NOMBRE_D_ITERATIONS,NOMBRE_D_ITERATIONS_IMPLICITE);
/* Initialisation du nombre d'iterations. */
#include xrv/champs_5.1A.I"
GET_ARGUMENTSv(nombre_d_arguments
,BLOC(PROCESS_ARGUMENTS_GEOMETRIQUES;
PROCESS_ARGUMENT_FICHIER("VARIABLE_cu0="
,fichier_VARIABLE_cu0
,liste_VARIABLE_cu0
,VARIABLE_cu0_IMPLICITE
,gTRANSFORMAT_31
);
PROCESS_ARGUMENT_FICHIER("VARIABLE_cv0="
,fichier_VARIABLE_cv0
,liste_VARIABLE_cv0
,VARIABLE_cv0_IMPLICITE
,gTRANSFORMAT_31
);
PROCESS_ARGUMENT_FICHIER("VARIABLE_dcz="
,fichier_VARIABLE_dcz
,liste_VARIABLE_dcz
,VARIABLE_dcz_IMPLICITE
,gTRANSFORMAT_31
);
PROCESS_ARGUMENT_FICHIER("VARIABLE_cx0="
,fichier_VARIABLE_cx0
,liste_VARIABLE_cx0
,VARIABLE_cx0_IMPLICITE
,gTRANSFORMAT_31
);
PROCESS_ARGUMENT_FICHIER("VARIABLE_cy0="
,fichier_VARIABLE_cy0
,liste_VARIABLE_cy0
,VARIABLE_cy0_IMPLICITE
,gTRANSFORMAT_31
);
PROCESS_ARGUMENT_FICHIER("VARIABLE_cz0="
,fichier_VARIABLE_cz0
,liste_VARIABLE_cz0
,VARIABLE_cz0_IMPLICITE
,gTRANSFORMAT_31
);
PROCESS_ARGUMENT_FICHIER("PAS_DE_TEMPS_dct="
,fichier_PAS_DE_TEMPS_dct
,liste_PAS_DE_TEMPS_dct
,PAS_DE_TEMPS_dct_IMPLICITE
,gTRANSFORMAT_31
);
PROCESS_ARGUMENT_FICHIER("NOMBRE_D_ITERATIONS="
,fichier_NOMBRE_D_ITERATIONS
,liste_NOMBRE_D_ITERATIONS
,NOMBRE_D_ITERATIONS_IMPLICITE
,gTRANSFORMAT_31
);
PROCESS_ARGUMENTS_DE_VISUALISATION;
PROCESS_ARGUMENTS_DE_VISUALISATION_DES_AXES_DE_COORDONNEES;
GET_ARGUMENT_I("n=""iterations=",nombre_d_iterations);
GET_ARGUMENT_L("diffusif=",un_systeme_diffusif);
GET_ARGUMENT_F("vz=",composante_en_Z_de_la_vitesse);
GET_ARGUMENT_F("z1=",z_du_plan_1);
GET_ARGUMENT_F("z2=",z_du_plan_2);
GET_ARGUMENT_F("dt=""dct=",dct);
)
);
#include xrv/champs_5.19.I"
/* Pour eviter le message : */
/* */
/* Static function is not referenced. */
/* */
/* sur 'SYSTEME_ES9000_AIX_CC'... */
#include xrk/attractor.19.I"
/* Validations et definition de l'espace physique. */
Komp(numero_de_la_periode_courante_de_la_simulation,nombre_de_periodes_de_la_simulation)
Bblock
RE_INITIALISATION_DE_L_HORLOGE;
INITIALISATIONS_RELATIVES_A_CHAQUE_NOUVELLE_IMAGE(numero_de_la_periode_courante);
/* Initialisations necessaires avant le calcul et la generation de chaque nouvelle image. */
EGAL(cu,sVARIABLE_cu0(numero_de_la_periode_courante));
EGAL(cv,sVARIABLE_cv0(numero_de_la_periode_courante));
/* Calcul des valeurs initiales des deux angles (u0,v0). */
Test(IL_NE_FAUT_PAS(un_systeme_diffusif))
Bblock
EGAL(cu_avant_la_reflexion,FZERO);
EGAL(cv_avant_la_reflexion,FZERO);
/* Initialisation des valeurs des fonctions 'Fx' et 'Fy' a l'instant precedent... */
Eblock
ATes
Bblock
Eblock
ETes
EGAL(composante_en_Z_de_la_vitesse,sVARIABLE_dcz(numero_de_la_periode_courante));
/* Calcul de la valeur initiale de la composante en 'Z' de la vitesse. */
EGAL(cx,sVARIABLE_cx0(numero_de_la_periode_courante));
EGAL(cy,sVARIABLE_cy0(numero_de_la_periode_courante));
EGAL(cz,sVARIABLE_cz0(numero_de_la_periode_courante));
/* Calcul des valeurs initiales des trois variables (x0,y0,z0). */
Test(NINCff(cz,z_du_plan_1,z_du_plan_2))
Bblock
PRINT_ERREUR("la coordonnee 'z' initiale n'est pas situee entre les deux plans paralleles");
CAL1(Prer3("z=%g plans=(%g,%g)\n",cz,z_du_plan_1,z_du_plan_2));
Eblock
ATes
Bblock
Eblock
ETes
EGAL(dcx,FZERO);
EGAL(dcy,FZERO);
EGAL(dcz,FZERO);
/* Initialisation des differentielles pour la premiere visualisation si celle-ci a lieu */
/* en couleurs. On notera que 'FZERO' est la valeur la plus "logique"... */
vTRANSFORMAT_31(nombre_d_iterations,sNOMBRE_D_ITERATIONS,numero_de_la_periode_courante,fichier_NOMBRE_D_ITERATIONS);
/* Calcul du nombre d'iterations lorsqu'il est variable. */
vTRANSFORMAT_31(dct,sPAS_DE_TEMPS_dct,numero_de_la_periode_courante,fichier_PAS_DE_TEMPS_dct);
/* Calcul du pas de temps. */
Test(IFLT(DIVI(SOUA(z_du_plan_1,z_du_plan_2)
,ABSO(composante_en_Z_de_la_vitesse)
)
,GRO1(GRO10(dct))
)
)
Bblock
PRINT_ERREUR("le pas de temps est trop grand, ce qui creera des artefacts au moment de la reflexion");
CAL1(Prer2("temps de vol du plan 1 au plan 2=%g dct=%g\n"
,DIVI(SOUA(z_du_plan_1,z_du_plan_2)
,ABSO(composante_en_Z_de_la_vitesse)
)
,dct
)
);
/* Voir a ce propos l'instruction : */
/* */
/* EGAL(cz,TRON(cz,...)); */
/* */
/* ci-apres et les commentaires associes... */
Eblock
ATes
Bblock
Eblock
ETes
Komp(numero_de_l_iteration_courante,nombre_d_iterations)
Bblock
RECHERCHE_DES_EXTREMA_DES_COORDONNEES_ET_DES_DERIVEES;
/* On notera que cette recherche n'est pas conditionnee par 'editer_les_extrema', car les */
/* extrema pourraient etre utilises pour la visualisation... */
Test(IFOU(IL_FAUT(visualiser_le_fantome)
,IFGE(numero_de_l_iteration_courante,PREMIERE_ITERATION_VISUALISEE)
)
)
Bblock
CALS(memorisation_1_point_07(SOUS(cx,Xcentre_ESPACE)
,SOUS(cy,Ycentre_ESPACE)
,SOUS(cz,Zcentre_ESPACE)
,dcx
,dcy
,dcz
,numero_de_l_iteration_courante
)
);
/* Memorisation de l'iteration courante... */
Eblock
ATes
Bblock
Eblock
ETes
INTEGRE(Icx,cx,dcx
,COND(IL_FAUT(un_systeme_diffusif)
,Fx(cu,cv)
,SOUS(Fx(cu,cv),Fx(cu_avant_la_reflexion,cv_avant_la_reflexion))
)
,dct
);
INTEGRE(Icy,cy,dcy
,COND(IL_FAUT(un_systeme_diffusif)
,Fy(cu,cv)
,SOUS(Fy(cu,cv),Fy(cu_avant_la_reflexion,cv_avant_la_reflexion))
)
,dct
);
INTEGRE(Icz,cz,dcz
,Fz(cu,cv)
,dct
);
MISE_A_JOUR(cx,Icx,cy,Icy,cz,Icz);
/* Integration de la trajectoire. On notera qu'il est possible d'observer des retours en */
/* arriere presque parfaits ; par exemple avec : */
/* */
/* (u ,v )=(1.0,1.0) */
/* 0 0 */
/* */
/* et : */
/* */
/* dct=0.2 */
/* */
/* on observe apres la dix-neuvieme reflexion : */
/* */
/* (u,v)=(3.351,5.768) ==> {dx,dy,dz}=(-0.9781,+0.8702,+1) */
/* */
/* et lors de la reflexion suivante : */
/* */
/* (u,v)=(6.188,2.836) ==> {dx,dy,dz}=(+0.9954,-0.9537,-1) */
/* */
/* ou l'on voit les signes des trois derivees s'inverser, alors que leurs valeurs absolues */
/* ne changent pratiquement pas... */
Test(NINCoo(cz,z_du_plan_1,z_du_plan_2))
Bblock
/* Cas ou la particule tente de franchir l'un des deux plans, il y a reflexion : */
DEFV(Float,INIT(cu_manoeuvre,FLOT__UNDEF));
DEFV(Float,INIT(cv_manoeuvre,FLOT__UNDEF));
/* Deux valeurs de manoeuvre pour faire le produit matriciel de reflexion... */
Test(IL_NE_FAUT_PAS(un_systeme_diffusif))
Bblock
INIT(cu_avant_la_reflexion,cu);
INIT(cv_avant_la_reflexion,cv);
/* Memorisation des angles 'u' et 'v' avant la reflexion... */
Eblock
ATes
Bblock
Eblock
ETes
EGAL(cu_manoeuvre,LIZ2(REFLEXION_11,cu,REFLEXION_12,cv));
EGAL(cv_manoeuvre,LIZ2(REFLEXION_21,cu,REFLEXION_22,cv));
EGAL(cu,CERC(cu_manoeuvre));
EGAL(cv,CERC(cv_manoeuvre));
/* Reflexion par modification des deux angles (u,v) suivant : */
/* */
/* | u | | 2 1 || u | */
/* | | = | || | (modulo 2.pi) */
/* | v | | 1 1 || v | */
/* */
/* [apres] [avant] */
/* */
/* On pourrait croire que le modulo '2.pi' n'est pas necessaire puisque les angles 'u' */
/* et 'v' sont utilises via des lignes trigonometriques, mais en fait, il n'en n'est rien */
/* car en effet, on est en presence ci-dessus approximativement d'une exponentielle en base */
/* 2 lors du calcul de 'u', or il est connu que cela va tres vite, d'ou des risques de */
/* debordement... */
EGAL(composante_en_Z_de_la_vitesse,NEGA(composante_en_Z_de_la_vitesse));
/* Et bien sur, il faut aussi inverser la composante en 'Z' de la vitesse : */
/* */
/* dz dz */
/* ---- = - ---- */
/* dt dt */
/* */
EGAL(cz,TROQ(cz,z_du_plan_1,z_du_plan_2));
/* Malheureusement, il faut tricher un peu en ramenant la coordonnee 'z' sur le plan */
/* adequat car en effet, si l'on restait la ou l'on est, et surtout a l'exterieur des */
/* des deux plans, a l'iteration suivante on pourrait, en cas de deplacement reduit, */
/* rester a l'exterieur des deux plans. Ceci explique que l'on valide le pas de temps */
/* utilise, car s'il est trop grand, cette operation 'TROQ(...)' peut modifier de facon */
/* considerable la trajectoire de la particule... */
Eblock
ATes
Bblock
/* Cas ou il n'y a pas reflexion, la particule continue a se deplacer en ligne droite... */
Eblock
ETes
INCREMENTATION_DE_L_HORLOGE(dct);
/* Simulation du temps de la simulation... */
Eblock
EKom
#include xrk/attractor.1A.I"
VISUALISATION_DES_AXES_DE_COORDONNEES;
/* Visualisation si necessaire des trois axes de coordonnees. */
GENERATION_D_UNE_IMAGE_ET_PASSAGE_A_LA_SUIVANTE(BLOC(VIDE;));
/* Generation de l'image courante... */
Eblock
EKom
EDITION_DES_EXTREMA_DES_COORDONNEES_ET_DES_DERIVEES;
/* Edition facultative des extrema des coordonnees et des derivees. */
RETU_Commande;
Eblock
ECommande