/*************************************************************************************************************************************/
/* */
/* D E F I N I T I O N S D E S F O N C T I O N S N E C E S S A I R E S */
/* A L ' E T U D E D E L ' A T O M E D ' H Y D R O G E N E : */
/* */
/* */
/* Author of '$xrq/Legendre.31$I' : */
/* */
/* Jean-Francois COLONNA (LACTAMME, 1993??????????). */
/* */
/*************************************************************************************************************************************/
/*===================================================================================================================================*/
/*************************************************************************************************************************************/
/* */
/* P O L Y N O M E S D E L E G E N D R E : */
/* */
/* */
/* definition du polynome de Legendre d'ordre (l,m) : */
/* */
/* m */
/* --- */
/* [ 2] 2 l+m l */
/* m [1 - x ] d [ 2 ] */
/* P (x) = -------------.-------[x - 1] */
/* l l l+m */
/* l!.2 dx */
/* */
/* avec : */
/* */
/* l = 0,1,2,...,+infini */
/* m = 0,1,2,...,l */
/* */
/* soit : */
/* */
/* 0 <= m <= l */
/* */
/* */
/* formule de recurrence : */
/* */
/* _______ */
/* m-1 m-1 / 2 m */
/* [l-(m-1)].[(l+1)-(m-1)].P (x) - [l+(m-1)].[(l+1)+(m-1)].P (x) = (2.l+1).(\/ 1 - x ).P (x) */
/* l+1 l-1 l */
/* */
/* ou : */
/* */
/* _______ */
/* m-1 m-1 / 2 m */
/* (l-m+1).(l-m+2).P (x) - (l+m-1).(l+m).P (x) = (2.l+1).(\/ 1 - x ).P (x) */
/* l+1 l-1 l */
/* */
/* Cette formule va donc permettre de calculer */
/* recursivement P(l,m,x) a partir des valeurs de */
/* P(l+1,m-1,x) et P(l-1,m-1,x) jusqu'a tomber sur */
/* des polynomes du type P(0,n,x) qui sont definis */
/* explicitement. Une petite difficulte surgit dans */
/* le cas ou : */
/* */
/* _______ */
/* / 2 */
/* \/ 1 - x = 0 */
/* */
/* Pour eviter une division par zero lors du calcul */
/* de P(l,m,x) a partir de P(l+1,m-1,x) et P(l-1,m-1,x) */
/* on substitue donc a 0 un petit epsilon en esperant */
/* que par continuite tout marche bien... */
/* */
/* */
/* definition du polynome de Legendre d'ordre (0,m) : */
/* */
/* 0 */
/* P (x) = 1 */
/* 0 */
/* */
/* l */
/* k=E(---) */
/* 2 */
/* __________ */
/* \ */
/* 0 \ k (2.l-2.k)! l-2.k */
/* P (x) = / (-1) .-----------------------.x */
/* l /_________ l */
/* 2 .k!.(l-k).!(l-2.k)! */
/* k=0 */
/* */
/*************************************************************************************************************************************/
/*===================================================================================================================================*/
/*************************************************************************************************************************************/
/* */
/* T A B L E D E S F A C T O R I E L L E S P R E - C A L C U L E E S : */
/* */
/*************************************************************************************************************************************/
#define PREMIERE_FACTORIELLE \
ZERO
#define DERNIERE_FACTORIELLE \
PRED(ENTIER_DE_DEBORDEMENT_DE_FACT_Float)
#define NOMBRE_DE_FACTORIELLES \
LENG(PREMIERE_FACTORIELLE,DERNIERE_FACTORIELLE)
/* Definition de la premiere et de la derniere entree de la liste des factorielles que */
/* l'on pre-calcule. ATTENTION, il est souhaitable que 'DERNIERE_FACTORIELLE' n'excede */
/* pas la "capacite" de la fonction 'FACT(...)' (voir 'v $ximf/produits$FON')... */
#define oFACT(n) \
ITb1(table_des_factorielles,INDX(n,PREMIERE_FACTORIELLE)) \
/* Fonction d'acces a la factorielle du nombre 'n'. */
#define INITIALISATION_DE_LA_TABLE_DES_FACTORIELLES \
Bblock \
Test(EST_INVALIDE(etat_de_la_table_des_factorielles)) \
Bblock \
DEFV(Int,INIT(nombre_courant,UNDEF)); \
/* Nombre courant dont on veut la factorielle... */ \
DoIn(nombre_courant,PREMIERE_FACTORIELLE,DERNIERE_FACTORIELLE,I) \
Bblock \
EGAL(oFACT(nombre_courant),FACT(nombre_courant)); \
/* Calcul de la factorielle du nombre courant... */ \
Eblock \
EDoI \
Eblock \
ATes \
Bblock \
Eblock \
ETes \
\
EGAL(etat_de_la_table_des_factorielles,VALIDE); \
/* Et maintenant, on sait que la table est initialisee... */ \
Eblock \
/* Initialisation, lorsqu'elle est necessaire de la table des factorielles pre-calculees... */
DEFV(Local,DEFV(Logical,INIT(etat_de_la_table_des_factorielles,INVALIDE)));
/* Indicateur precisant si la table est initialisee ('VALIDE') ou pas ('INVALIDE'). */
DEFV(Local,DEFV(Float,DTb1(table_des_factorielles,NOMBRE_DE_FACTORIELLES)));
/* Table des factorielles pre-calculees... */
/*===================================================================================================================================*/
/*************************************************************************************************************************************/
/* */
/* C A L C U L O P T I M I S E D ' U N P O L Y N O M E D E L E G E N D R E D E T Y P E P(l,0,x) : */
/* */
/*************************************************************************************************************************************/
BFonctionF
DEFV(Local,DEFV(FonctionF,polynome_de_Legendre_l_0(ordre_l,variable)))
DEFV(Argument,DEFV(Int,ordre_l));
/* Ordre 'l' du polynome de Legendre. */
DEFV(Argument,DEFV(Float,variable));
/* Variable pour laquelle evaluer le polynome de Legendre. */
/*-----------------------------------------------------------------------------------------------------------------------------------*/
Bblock
DEFV(Float,INIT(valeur_du_polynome,FZERO));
/* Valeur courante du polynome lors du processus iteratif. Cette valeur convient tout a */
/* fait a la methode de Horner qui sera utilisee par la suite. */
INIT_ERROR;
/*..............................................................................................................................*/
Test(IZLT(ordre_l))
Bblock
PRINT_ERREUR("l'ordre 'l' est negatif, une valeur nulle est renvoyee");
Eblock
ATes
Bblock
Test(IFGT(DOUB(ordre_l),DERNIERE_FACTORIELLE))
Bblock
PRINT_ERREUR("la table de pre-calcul des factorielles n'est pas suffisamment longue");
CAL1(Prer1("l=%d\n",ordre_l));
CAL1(Prer2("%d>%d\n",DOUB(ordre_l),DERNIERE_FACTORIELLE));
Eblock
ATes
Bblock
Eblock
ETes
INITIALISATION_DE_LA_TABLE_DES_FACTORIELLES;
/* Initialisation, lorsqu'elle est necessaire de la table des factorielles pre-calculees... */
Test(IZEQ(ordre_l))
Bblock
EGAL(valeur_du_polynome,FU);
/* Le cas 'P(0,0,x)' est traite a part : */
/* */
/* 0 */
/* P (x) = 1 */
/* 0 */
/* */
Eblock
ATes
Bblock
DEFV(Float,INIT(variable_au_carre,EXP2(variable)));
/* Pre-calcul de la variable au carre qui ne varie pas ou cours de l'iteration... */
DEFV(Float,INIT(signe_du_monome_courant,SIGNE_PLUS));
/* Signe du monome courant qui s'alterne a chaque iteration... */
DEFV(Int,INIT(indice_du_polynome,UNDEF));
/* Indice 'k' de calcul du polynome... */
DoIn(indice_du_polynome,ZERO,MOIT(ordre_l),I)
/* On notera que 'MOIT(ordre_l)' prend la moitie de l'ordre 'l' par defaut... */
Bblock
EGAL(valeur_du_polynome
,AXPB(valeur_du_polynome
,variable_au_carre
,MUL2(signe_du_monome_courant
,DIVI(oFACT(SOUS(DOUB(ordre_l),DOUB(indice_du_polynome)))
,MUL3(oFACT(indice_du_polynome)
,oFACT(SOUS(NEUT(ordre_l),NEUT(indice_du_polynome)))
,oFACT(SOUS(NEUT(ordre_l),DOUB(indice_du_polynome)))
)
)
)
)
);
/* Calcul optimise a l'aide de la methode de Horner du polynome de Legendre : */
/* */
/* l */
/* k=E(---) */
/* 2 */
/* __________ */
/* \ */
/* 0 \ k (2.l-2.k)! l-2.k */
/* P (x) = / (-1) .-----------------------.x */
/* l /_________ l */
/* 2 .k!.(l-k).!(l-2.k)! */
/* k=0 */
/* */
/* ou encore : */
/* */
/* l */
/* k=E(---) */
/* 2 */
/* __________ */
/* \ */
/* 0 1 \ k (2.l-2.k)! l-2.k */
/* P (x) = ----. / (-1) .--------------------.x */
/* l l /_________ */
/* 2 k!.(l-k).!(l-2.k)! */
/* k=0 */
/* */
EGAL(signe_du_monome_courant,NEGA(signe_du_monome_courant));
/* Inversion du signe du monome suivant... */
Eblock
EDoI
Test(EST_IMPAIR(ordre_l))
Bblock
EGAL(valeur_du_polynome
,MUL2(valeur_du_polynome
,variable
)
);
/* Lorsque le degre du polynome de Legendre est impair, il faut compenser le fait que la */
/* methode Horner ci-dessus conduit a un polynome dont toutes les puissances de la variable */
/* sont paires... */
Eblock
ATes
Bblock
Eblock
ETes
EGAL(valeur_du_polynome,DIVI(valeur_du_polynome,MONX(FDEUX,ordre_l)));
/* Puis division par 2 a la puissance 'l' qui ne change pas a l'interieur de la boucle de */
/* calcul du polynome d'ordre 'l'... */
Eblock
ETes
Eblock
ETes
RETU(valeur_du_polynome);
Eblock
EFonctionF
/*===================================================================================================================================*/
/*************************************************************************************************************************************/
/* */
/* C A L C U L O P T I M I S E D ' U N P O L Y N O M E D E L E G E N D R E D E T Y P E P(l,m,x) : */
/* */
/*************************************************************************************************************************************/
BFonctionF
DEFV(Local,DEFV(FonctionF,polynome_de_Legendre_l_m(ordre_l,ordre_m,variable)))
DEFV(Argument,DEFV(Int,ordre_l));
/* Ordre 'l' du polynome de Legendre. */
DEFV(Argument,DEFV(Int,ordre_m));
/* Ordre 'm' du polynome de Legendre. */
DEFV(Argument,DEFV(Float,variable));
/* Variable pour laquelle evaluer le polynome de Legendre. */
/*-----------------------------------------------------------------------------------------------------------------------------------*/
Bblock
DEFV(Float,INIT(valeur_du_polynome,FZERO));
/* Valeur courante du polynome lors du processus iteratif... */
INIT_ERROR;
/*..............................................................................................................................*/
Test(IFGT(ordre_m,ordre_l))
Bblock
PRINT_ERREUR("l'ordre 'm' est superieur a l'ordre 'l', une valeur nulle est renvoyee");
Eblock
ATes
Bblock
Test(IFOU(IZLT(ordre_l),IZLT(ordre_m)))
Bblock
PRINT_ERREUR("l'ordre 'l' et/ou l'ordre 'm' sont negatifs, une valeur nulle est renvoyee");
Eblock
ATes
Bblock
Test(IZGT(ordre_m))
Bblock
DEFV(Float,INIT(expression_pouvant_etre_negative_ou_nulle,SOUS(FU,EXP2(variable))));
/* Introduction d'une valeur intermediaire correspondant a une expression pouvant etre */
/* negative ou nulle... */
Test(IZGE(expression_pouvant_etre_negative_ou_nulle))
Bblock
EGAL(expression_pouvant_etre_negative_ou_nulle
,MAX2(expression_pouvant_etre_negative_ou_nulle
,pEPSILON
)
);
/* Lorsque le diviseur se trouve etre nul, on lui substitue une toute petite valeur... */
EGAL(valeur_du_polynome
,DIVI(SOUS(MUL3(SOUS(NEUT(ordre_l),PRED(ordre_m))
,SOUS(SUCC(ordre_l),PRED(ordre_m))
,polynome_de_Legendre_l_m(SUCC(ordre_l),PRED(ordre_m),variable)
)
,MUL3(ADD2(NEUT(ordre_l),PRED(ordre_m))
,ADD2(SUCC(ordre_l),PRED(ordre_m))
,polynome_de_Legendre_l_m(PRED(ordre_l),PRED(ordre_m),variable)
)
)
,MUL2(DOUP(ordre_l)
,RACX(expression_pouvant_etre_negative_ou_nulle)
)
)
);
/* Calcul recurrent (non optimise...) du polynome de Legendre : */
/* */
/* _______ */
/* / 2 m */
/* (2.l+1).(\/ 1 - x ).P (x) = */
/* l */
/* */
/* m-1 m-1 */
/* [l-(m-1)].[(l+1)-(m-1)].P (x) - [l+(m-1)].[(l+1)+(m-1)].P (x) */
/* l+1 l-1 */
/* */
Eblock
ATes
Bblock
PRINT_ERREUR("la recurrence ne peut avoir lieu, une valeur nulle est renvoyee");
Eblock
ETes
Eblock
ATes
Bblock
EGAL(valeur_du_polynome
,polynome_de_Legendre_l_0(ordre_l,variable)
);
/* La recursivite s'arrete lorsque l'ordre 'm' est nul... */
Eblock
ETes
Eblock
ETes
Eblock
ETes
RETU(valeur_du_polynome);
Eblock
EFonctionF
#undef INITIALISATION_DE_LA_TABLE_DES_FACTORIELLES
#undef oFACT
#undef NOMBRE_DE_FACTORIELLES
#undef DERNIERE_FACTORIELLE
#undef PREMIERE_FACTORIELLE