/*************************************************************************************************************************************/
/* */
/* 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.21$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 */
/* */
/*************************************************************************************************************************************/
/*===================================================================================================================================*/
/*************************************************************************************************************************************/
/* */
/* C A L C U L 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... */
INIT_ERROR;
/*..............................................................................................................................*/
Test(IZLT(ordre_l))
Bblock
PRINT_ERREUR("l'ordre 'l' est negatif, une valeur nulle est renvoyee");
Eblock
ATes
Bblock
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(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
INCR(valeur_du_polynome
,MUL3(MONX(NEGA(FU),indice_du_polynome)
,DIVI(FACT(SOUS(DOUB(ordre_l),DOUB(indice_du_polynome)))
,MUL4(MONX(FDEUX,ordre_l)
,FACT(indice_du_polynome)
,FACT(SOUS(NEUT(ordre_l),NEUT(indice_du_polynome)))
,FACT(SOUS(NEUT(ordre_l),DOUB(indice_du_polynome)))
)
)
,MONX(variable,SOUS(NEUT(ordre_l),DOUB(indice_du_polynome)))
)
);
/* Calcul (non optimise...) et cumul des differents monomes 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 */
/* */
Eblock
EDoI
Eblock
ETes
Eblock
ETes
RETU(valeur_du_polynome);
Eblock
EFonctionF
/*===================================================================================================================================*/
/*************************************************************************************************************************************/
/* */
/* C A L C U L 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