/*************************************************************************************************************************************/
/* */
/* M E T H O D E S D ' I N T E G R A T I O N N U M E R I Q U E P O U R D E S */
/* S Y S T E M E S 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 D U S E C O N D O R D R E : */
/* */
/* */
/* Definition : */
/* */
/* Soit le systeme d'equations differentielles */
/* suivant : */
/* */
/* 2 */
/* d x dx */
/* F (x,y,z,t).----- + F (x,y,z,t).---- + F (x,y,z,t).x = F (t) */
/* 3x 2 2x dt 1x 0x */
/* dt */
/* */
/* 2 */
/* d y dy */
/* F (x,y,z,t).----- + F (x,y,z,t).---- + F (x,y,z,t).y = F (t) */
/* 3y 2 2y dt 1y 0y */
/* dt */
/* */
/* 2 */
/* d z dz */
/* F (x,y,z,t).----- + F (x,y,z,t).---- + F (x,y,z,t).z = F (t) */
/* 3z 2 2z dt 1z 0z */
/* dt */
/* */
/* */
/* */
/* Resolution (voir 'METHODE DE CALCUL NUMERIQUE' de JP. Nougier, page 214) */
/* par une methode explicite du deuxieme ordre : */
/* */
/* posons : */
/* */
/* h = dt */
/* */
/* et : */
/* */
/* F (x ,y ,z ,n) F (x ,y ,z ,n) */
/* 3? n n n 2? n n n */
/* A = ----------------- - ----------------- */
/* ? 2 2.h */
/* h */
/* */
/* 2.F (x ,y ,z ,n) */
/* 3? n n n */
/* B = F (x ,y ,z ,n) - ------------------- */
/* ? 1? n n n 2 */
/* h */
/* */
/* F (x ,y ,z ,n) F (x ,y ,z ,n) */
/* 3? n n n 2? n n n */
/* C = ----------------- + ----------------- */
/* ? 2 2.h */
/* h */
/* */
/* ou "?" designe indifferemment "x", "y" ou "z". */
/* */
/* On a alors : */
/* */
/* d? */
/* (F (1) + (3.C - A ).? + 2.C .----(0).h */
/* 0? ? ? 0 ? dt */
/* ? = ------------------------------------------- */
/* 1 B + 4.C */
/* ? ? */
/* */
/* puis : */
/* */
/* (F (n) - B .? - A .? ) */
/* 0? ? n ? n-1 */
/* ? = ---------------------------- */
/* n+1 C */
/* ? */
/* */
/* pour : */
/* */
/* n > 0 */
/* */
/* On notera que la reference : */
/* */
/* F (1) */
/* 0? */
/* */
/* devrait etre d'une maniere tres generale : */
/* */
/* F (x ,y ,z ,1) */
/* 0? 1 1 1 */
/* */
/* dans l'expression '?1' ci-dessus. Or clairement */
/* cela donnerait une definition recursive (puisque */
/* ce sont les variables '?1' que l'on est precisemment */
/* en train d'evaluer. C'est pourquoi on ne fait dependre */
/* les fonctions 'F0?' que de la variable 't'... */
/* */
/* */
/* Author of '$xrk/integr.2B$vv$I' : */
/* */
/* Jean-Francois Colonna (LACTAMME, 1994??????????). */
/* */
/*************************************************************************************************************************************/
/*===================================================================================================================================*/
/*************************************************************************************************************************************/
/* */
/* D O N N E E S U T I L E S : */
/* */
/*************************************************************************************************************************************/
DEFV(Local,DEFV(Float,INIT(cx0,FLOT__UNDEF)));
DEFV(Local,DEFV(Float,INIT(dcx0,FLOT__UNDEF)));
DEFV(Local,DEFV(Float,INIT(cx,FLOT__UNDEF)));
DEFV(Local,DEFV(Float,INIT(cx_1,FLOT__UNDEF)));
/* Definition de 'x' et 'dx', ainsi que de leurs valeurs initiales (0) et anterieures. */
DEFV(Local,DEFV(Float,INIT(cy0,FLOT__UNDEF)));
DEFV(Local,DEFV(Float,INIT(dcy0,FLOT__UNDEF)));
DEFV(Local,DEFV(Float,INIT(cy,FLOT__UNDEF)));
DEFV(Local,DEFV(Float,INIT(cy_1,FLOT__UNDEF)));
/* Definition de 'y' et 'dy', ainsi que de leurs valeurs initiales (0) et anterieures. */
DEFV(Local,DEFV(Float,INIT(cz0,FLOT__UNDEF)));
DEFV(Local,DEFV(Float,INIT(dcz0,FLOT__UNDEF)));
DEFV(Local,DEFV(Float,INIT(cz,FLOT__UNDEF)));
DEFV(Local,DEFV(Float,INIT(cz_1,FLOT__UNDEF)));
/* Definition de 'z' et 'dz', ainsi que de leurs valeurs initiales (0) et anterieures. */
DEFV(Local,DEFV(Logical,INIT(il_s_agit_du_premier_pas_de_temps,VRAI)));
DEFV(Local,DEFV(Logical,INIT(il_s_agit_du_second_pas_de_temps,FAUX)));
/* Afin de traiter differement le premier et le second pas de temps... */
/*===================================================================================================================================*/
/*************************************************************************************************************************************/
/* */
/* F O R M U L E S D E B A S E : */
/* */
/*************************************************************************************************************************************/
#define FONCTION_A(F3,F2,F1,F0,h) \
SOUS(DIVI(F3,EXP2(h)),DIVI(F2,DOUB(h))) \
/* Calcul des fonction 'A?'. */
#define FONCTION_B(F3,F2,F1,F0,h) \
SOUS(F1,DIVI(DOUB(F3),EXP2(h))) \
/* Calcul des fonction 'B?'. */
#define FONCTION_C(F3,F2,F1,F0,h) \
ADD2(DIVI(F3,EXP2(h)),DIVI(F2,DOUB(h))) \
/* Calcul des fonction 'C?'. */
#define MISE_A_JOUR(v12,v11,v22,v21,v32,v31) \
Bblock \
EGAL(v12,v11); \
EGAL(v22,v21); \
EGAL(v32,v31); \
Eblock \
/* Procedure de mise a jour finale des variables... */
/*===================================================================================================================================*/
/*************************************************************************************************************************************/
/* */
/* F O R M U L E S G E N E R A L E S D ' I N T E G R A T I O N : */
/* */
/*************************************************************************************************************************************/
#define INITIALISATION_DE_L_INTEGRATION_D_UN_SYSTEME_D_EQUATIONS_DIFFERENTIELLES_O2 \
Bblock \
EGAL(il_s_agit_du_premier_pas_de_temps,VRAI); \
EGAL(il_s_agit_du_second_pas_de_temps,FAUX); \
/* Afin de traiter differement le premier et le second pas de temps... */ \
Eblock \
/* Initialisation de l'integration du systeme d'equations differentielles du second ordre... */
#define INITIALISATION_POUR_UN_POINT_DE_L_INTEGRATION_D_UN_SYSTEME_D_EQUATIONS_DIFFERENTIELLES_O2(x0,y0,z0,dx0,dy0,dz0) \
Bblock \
EGAL(cx0,x0); \
EGAL(dcx0,dx0); \
/* Initialisation de la variable 'x'. */ \
EGAL(cy0,y0); \
EGAL(dcy0,dy0); \
/* Initialisation de la variable 'y'. */ \
EGAL(cz0,z0); \
EGAL(dcz0,dz0); \
/* Initialisation de la variable 'z'. */ \
Eblock \
/* Initialisation de l'integration du systeme d'equations differentielles du second ordre... */
#define DEPILER_UN_POINT_LORS_DE_L_INTEGRATION_O2(Pcx0,Pdcx0,Pcx,Pdcx,Pcx_1,Pcy0,Pdcy0,Pcy,Pdcy,Pcy_1,Pcz0,Pdcz0,Pcz,Pdcz,Pcz_1) \
Bblock \
EGAL(cx0,Pcx0); \
EGAL(dcx0,Pdcx0); \
EGAL(cx,Pcx); \
EGAL(dcx,Pdcx); \
EGAL(cx_1,Pcx_1); \
/* Definition de 'x' et 'dx', ainsi que de leurs valeurs initiales (0) et anterieures. */ \
EGAL(cy0,Pcy0); \
EGAL(dcy0,Pdcy0); \
EGAL(cy,Pcy); \
EGAL(dcy,Pdcy); \
EGAL(cy_1,Pcy_1); \
/* Definition de 'y' et 'dy', ainsi que de leurs valeurs initiales (0) et anterieures. */ \
EGAL(cz0,Pcz0); \
EGAL(dcz0,Pdcz0); \
EGAL(cz,Pcz); \
EGAL(dcz,Pdcz); \
EGAL(cz_1,Pcz_1); \
/* Definition de 'z' et 'dz', ainsi que de leurs valeurs initiales (0) et anterieures. */ \
Eblock \
/* Mise en place d'un point particulier avant l'integration. */
#define EMPILER_UN_POINT_LORS_DE_L_INTEGRATION_O2(Pcx0,Pdcx0,Pcx,Pdcx,Pcx_1,Pcy0,Pdcy0,Pcy,Pdcy,Pcy_1,Pcz0,Pdcz0,Pcz,Pdcz,Pcz_1) \
Bblock \
EGAL(Pcx0,cx0); \
EGAL(Pdcx0,dcx0); \
EGAL(Pcx,cx); \
EGAL(Pdcx,dcx); \
EGAL(Pcx_1,cx_1); \
/* Definition de 'x' et 'dx', ainsi que de leurs valeurs initiales (0) et anterieures. */ \
EGAL(Pcy0,cy0); \
EGAL(Pdcy0,dcy0); \
EGAL(Pcy,cy); \
EGAL(Pdcy,dcy); \
EGAL(Pcy_1,cy_1); \
/* Definition de 'y' et 'dy', ainsi que de leurs valeurs initiales (0) et anterieures. */ \
EGAL(Pcz0,cz0); \
EGAL(Pdcz0,dcz0); \
EGAL(Pcz,cz); \
EGAL(Pdcz,dcz); \
EGAL(Pcz_1,cz_1); \
/* Definition de 'z' et 'dz', ainsi que de leurs valeurs initiales (0) et anterieures. */ \
Eblock \
/* Rangement d'un point particulier apres l'integration. */
#define INTEGRATION_POUR_UN_POINT_D_UN_SYSTEME_D_EQUATIONS_DIFFERENTIELLES_O2(t,h) \
Bblock \
DEFV(Float,INIT(Icx,FLOT__UNDEF)); \
DEFV(Float,INIT(Icy,FLOT__UNDEF)); \
DEFV(Float,INIT(Icz,FLOT__UNDEF)); \
/* Variables intermediaires destinees a contenir les valeurs de courante de {cx,cy,cz} */ \
/* tant que les valeurs anciennes sont encore utiles... */ \
\
Test(EST_VRAI(il_s_agit_du_premier_pas_de_temps)) \
Bblock \
EGAL(cx,cx0); \
EGAL(cx_1,cx0); \
EGAL(cy,cy0); \
EGAL(cy_1,cy0); \
EGAL(cz,cz0); \
EGAL(cz_1,cz0); \
/* Pour le premier pas de temps, le calcul est on ne peut plus simple... */ \
/* */ \
/* ATTENTION, il y a eu pendant longtemps : */ \
/* */ \
/* EGAL(cx,cx0); */ \
/* EGAL(cx_1,FLOT__UNDEF); */ \
/* EGAL(cy,cy0); */ \
/* EGAL(cy_1,FLOT__UNDEF); */ \
/* EGAL(cz,cz0); */ \
/* EGAL(cz_1,FLOT__UNDEF); */ \
/* */ \
/* Or l'evaluation des fonction 'F??(...)' peut demander les coordonnees courantes comme */ \
/* les coordonnees precedentes (voir par exemple 'F1?(...)' dans 'v $xrr/N_corps.11$K'). */ \
/* D'ou cette initialisation telle que : */ \
/* */ \
/* {cx,cy,cz} = {cx_1,cy_1,cz_1} */ \
/* ainsi, deux corps differents n'auront pas les memes coordonnees precedentes a l'etat */ \
/* initial (sauf bien sur si leurs positions initiales sont identiques...). */ \
Eblock \
ATes \
Bblock \
DEFV(Float,INIT(fonction_Ax,FONCTION_A(F3x(cx,cy,cz,t),F2x(cx,cy,cz,t),F1x(cx,cy,cz,t),F0x(t),h))); \
DEFV(Float,INIT(fonction_Ay,FONCTION_A(F3y(cx,cy,cz,t),F2y(cx,cy,cz,t),F1y(cx,cy,cz,t),F0y(t),h))); \
DEFV(Float,INIT(fonction_Az,FONCTION_A(F3z(cx,cy,cz,t),F2z(cx,cy,cz,t),F1z(cx,cy,cz,t),F0z(t),h))); \
/* Definition des fonction 'A?'. */ \
DEFV(Float,INIT(fonction_Bx,FONCTION_B(F3x(cx,cy,cz,t),F2x(cx,cy,cz,t),F1x(cx,cy,cz,t),F0x(t),h))); \
DEFV(Float,INIT(fonction_By,FONCTION_B(F3y(cx,cy,cz,t),F2y(cx,cy,cz,t),F1y(cx,cy,cz,t),F0y(t),h))); \
DEFV(Float,INIT(fonction_Bz,FONCTION_B(F3z(cx,cy,cz,t),F2z(cx,cy,cz,t),F1z(cx,cy,cz,t),F0z(t),h))); \
/* Definition des fonction 'B?'. */ \
DEFV(Float,INIT(fonction_Cx,FONCTION_C(F3x(cx,cy,cz,t),F2x(cx,cy,cz,t),F1x(cx,cy,cz,t),F0x(t),h))); \
DEFV(Float,INIT(fonction_Cy,FONCTION_C(F3y(cx,cy,cz,t),F2y(cx,cy,cz,t),F1y(cx,cy,cz,t),F0y(t),h))); \
DEFV(Float,INIT(fonction_Cz,FONCTION_C(F3z(cx,cy,cz,t),F2z(cx,cy,cz,t),F1z(cx,cy,cz,t),F0z(t),h))); \
/* Definition des fonction 'C?'. */ \
\
Test(EST_VRAI(il_s_agit_du_second_pas_de_temps)) \
Bblock \
DEFV(Float,INIT(fonction_Bx_p_4Cx,ADD2(GRO1(fonction_Bx),GRO4(fonction_Cx)))); \
DEFV(Float,INIT(fonction_By_p_4Cy,ADD2(GRO1(fonction_By),GRO4(fonction_Cy)))); \
DEFV(Float,INIT(fonction_Bz_p_4Cz,ADD2(GRO1(fonction_Bz),GRO4(fonction_Cz)))); \
/* Definition des fonction 'B? + 4.C?'. */ \
\
Test(I3ET(IZNE(fonction_Bx_p_4Cx),IZNE(fonction_By_p_4Cy),IZNE(fonction_Bz_p_4Cz))) \
Bblock \
EGAL(Icx \
,DIVI(ADD3(F0x(t) \
,MUL2(SOUS(GRO3(fonction_Cx),GRO1(fonction_Ax)),cx0) \
,MUL3(GRO2(fonction_Cx),dcx0,h) \
) \
,fonction_Bx_p_4Cx \
) \
); \
EGAL(Icy \
,DIVI(ADD3(F0y(t) \
,MUL2(SOUS(GRO3(fonction_Cy),GRO1(fonction_Ay)),cy0) \
,MUL3(GRO2(fonction_Cy),dcy0,h) \
) \
,fonction_By_p_4Cy \
) \
); \
EGAL(Icz \
,DIVI(ADD3(F0z(t) \
,MUL2(SOUS(GRO3(fonction_Cz),GRO1(fonction_Az)),cz0) \
,MUL3(GRO2(fonction_Cz),dcz0,h) \
) \
,fonction_Bz_p_4Cz \
) \
); \
/* Pour le second pas de temps, le calcul est beaucoup plus complique... */ \
Eblock \
ATes \
Bblock \
PRINT_ERREUR("au moins une des fonctions 'Bx_p_4Cx', 'By_p_4Cy' ou 'Bz_p_4Cz' est nulle"); \
CAL1(Prer1("Bx_p_4Cx=%f\n",fonction_Bx_p_4Cx)); \
CAL1(Prer1("By_p_4Cy=%f\n",fonction_By_p_4Cy)); \
CAL1(Prer1("Bz_p_4Cz=%f\n",fonction_Bz_p_4Cz)); \
Eblock \
ETes \
Eblock \
ATes \
Bblock \
Test(I3ET(IZNE(fonction_Cx),IZNE(fonction_Cy),IZNE(fonction_Cz))) \
Bblock \
EGAL(Icx \
,DIVI(ADD3(F0x(t) \
,NEGA(MUL2(fonction_Bx,cx)) \
,NEGA(MUL2(fonction_Ax,cx_1)) \
) \
,GRO1(fonction_Cx) \
) \
); \
EGAL(Icy \
,DIVI(ADD3(F0y(t) \
,NEGA(MUL2(fonction_By,cy)) \
,NEGA(MUL2(fonction_Ay,cy_1)) \
) \
,GRO1(fonction_Cy) \
) \
); \
EGAL(Icz \
,DIVI(ADD3(F0z(t) \
,NEGA(MUL2(fonction_Bz,cz)) \
,NEGA(MUL2(fonction_Az,cz_1)) \
) \
,GRO1(fonction_Cz) \
) \
); \
/* Cas des pas de temps posterieurs au second... */ \
Eblock \
ATes \
Bblock \
PRINT_ERREUR("au moins une des fonctions 'Cx', 'Cy' ou 'Cz' est nulle"); \
CAL1(Prer1("Cx=%f\n",fonction_Cx)); \
CAL1(Prer1("Cy=%f\n",fonction_Cy)); \
CAL1(Prer1("Cz=%f\n",fonction_Cz)); \
Eblock \
ETes \
Eblock \
ETes \
\
EGAL(cx_1,cx); \
EGAL(cy_1,cy); \
EGAL(cz_1,cz); \
/* Memorisation du pas de temps precedent... */ \
MISE_A_JOUR(cx,Icx,cy,Icy,cz,Icz); \
/* Mise a jour finale des variables... */ \
Eblock \
ETes \
\
Test(EST_VRAI(il_s_agit_du_premier_pas_de_temps)) \
Bblock \
EGAL(dcx,dcx0); \
EGAL(dcy,dcy0); \
EGAL(dcz,dcz0); \
/* La vitesse courante est la vitesse initiale... */ \
Eblock \
ATes \
Bblock \
EGAL(dcx,DERIVATION_PARTIELLE(cx_1,cx,h)); \
EGAL(dcy,DERIVATION_PARTIELLE(cy_1,cy,h)); \
EGAL(dcz,DERIVATION_PARTIELLE(cz_1,cz,h)); \
/* La vitesse courante est approximee par une difference sur les coordonnees... */ \
Eblock \
ETes \
\
Eblock \
/* Integration du systeme d'equations differentielles du second ordre... */
#define GESTION_DE_L_INTEGRATION_D_UN_SYSTEME_D_EQUATIONS_DIFFERENTIELLES_O2 \
Bblock \
Test(EST_VRAI(il_s_agit_du_premier_pas_de_temps)) \
Bblock \
EGAL(il_s_agit_du_premier_pas_de_temps,FAUX); \
EGAL(il_s_agit_du_second_pas_de_temps,VRAI); \
/* Afin de traiter a part le second pas de temps... */ \
Eblock \
ATes \
Bblock \
Test(EST_VRAI(il_s_agit_du_second_pas_de_temps)) \
Bblock \
EGAL(il_s_agit_du_second_pas_de_temps,FAUX); \
/* Le premier et le second pas de temps on ete traites... */ \
Eblock \
ATes \
Bblock \
Eblock \
ETes \
Eblock \
ETes \
\
Eblock \
/* Integration du systeme d'equations differentielles du second ordre... */
/*===================================================================================================================================*/
/*************************************************************************************************************************************/
/* */
/* P O U R D E S R A I S O N S D E C O M P A T I B I L I T E : */
/* */
/*************************************************************************************************************************************/
DEFV(Local,DEFV(Float,INIT(dcx,FLOT__UNDEF)));
DEFV(Local,DEFV(Float,INIT(dcy,FLOT__UNDEF)));
DEFV(Local,DEFV(Float,INIT(dcz,FLOT__UNDEF)));
/* Afin d'assurer la compatibilite avec '$xrk/integr.1B$vv$I'. Malgre tout, le 1996021500 */
/* elles deviennent utiles car elles contiennent le vecteur vitesse courant approxime en */
/* sortie de 'INTEGRATION_POUR_UN_POINT_D_UN_SYSTEME_D_EQUATIONS_DIFFERENTIELLES_O2(...)'. */
/* Elles sont empilees par 'EMPILER_UN_POINT_LORS_DE_L_INTEGRATION_O2(...)' et bien entendu */
/* depilees par 'DEPILER_UN_POINT_LORS_DE_L_INTEGRATION_O2(...)' meme si elles sont inutiles */
/* en realite au processus d'integration ; en fait elles ne sont utiles que si l'on souhaite */
/* editer les coordonnees et les vitesses (voir 'v $xrr/N_corps.11$K'). */