## "Fossies" - the Fresh Open Source Software Archive

### Source code changes of the file "Numeric/Legendre.cpp" betweengetdp-3.4.0-source.tgz and getdp-3.5.0-source.tgz

About: GetDP is a general finite element solver using mixed elements to discretize de Rham-type complexes in one, two and three dimensions.

Legendre.cpp  (getdp-3.4.0-source.tgz):Legendre.cpp  (getdp-3.5.0-source.tgz)
// GetDP - Copyright (C) 1997-202 P. Dular and C. Geuzaine, University of Liege // GetDP - Copyright (C) 1997-2022 P. Dular and C. Geuzaine, University of Liege
// //
// issues on https://gitlab.onelab.info/getdp/getdp/issues. // issues on https://gitlab.onelab.info/getdp/getdp/issues.
// //
// Contributor(s): // Contributor(s):
// Ruth Sabariego // Ruth Sabariego
// //
#include <math.h> #include <math.h>
#include <stdio.h> #include <stdio.h>
#include "Message.h" #include "Message.h"
#define THESIGN(a) ? 1 : -1) #define THESIGN(a) ((a) >= 0 ? 1 : -1)
#define THEABS(a) ? a : -a) #define THEABS(a) ((a) >= 0 ? a : -a)
#define ONE_OVER_FOUR_PI 7.9577471545947668E-02 #define ONE_OVER_FOUR_PI 7.9577471545947668E-02
double Factorial(double n) double Factorial(double n)
{ {
/* FACTORIAL(n) is the product of all the integers from 1 to n */ /* FACTORIAL(n) is the product of all the integers from 1 to n */
double F; double F;
< if(n < 0) {
Message::Error("Factorial(n): n must be a positive Message::Error("Factorial(n): n must be a positive integer");
return 0; return 0;
} }
== return if(n == 0) return 1.;
<= return if(n <= 2) return n;
F = F = n;
> F *= } while(n > 2) {
n--;
F *= n;
}
return F; return F;
} }
double BinomialCoef(double n, double m) double BinomialCoef(double n, double m)
{ {
/* Binomial Coefficients: (n m) Computes de number of ways of /* Binomial Coefficients: (n m) Computes de number of ways of
choosing m objects from a collection of n distinct objects */ choosing m objects from a collection of n distinct objects */
int int i;
double B = double B = 1.;
|| return if(m == 0 || n == m) return 1.;
for(i = i > i--) for(i = (int)n; i > m; i--) B *= (double)i / (double)(i - m);
B *=
return B; return B;
} }
double Legendre(int l, int m, double x) double Legendre(int l, int m, double x)
{ {
/* Computes the associated Legendre polynomial P_l^m(x). Here the /* Computes the associated Legendre polynomial P_l^m(x). Here the
degree l and the order m are the integers satisfying -l<=m<=l, degree l and the order m are the integers satisfying -l<=m<=l,
while x lies in the range -1<=x<=1 */ while x lies in the range -1<=x<=1 */
double fact, pll; double fact, pll = 0., pmm, pmmp1, somx2, Cte;
int i, ll; int i, ll;
if{ if(THEABS(m) > l || fabs(x) > 1.) {
Message::Error("Bad arguments for Legendre: P_l^m(x) with -l<=m<=l (int)," Message::Error("Bad arguments for Legendre: P_l^m(x) with -l<=m<=l (int),"
" -1<=x<=1 l = %d m = %d x = %.8g", l, m, x); " -1<=x<=1 l = %d m = %d x = %.8g",
l, m, x);
return 0.; return 0.;
} }
Cte = (m > 0) ? 1. : / Cte = (m > 0) ?
* 1. :
m = Factorial((double)(l - THEABS(m))) /
Factorial((double)(l + THEABS(m))) * pow(-1., (double)THEABS(m));
pmm = m = THEABS(m);
> 0) { pmm = 1.;
somx2 =
fact = if(m > 0) {
somx2 = sqrt((1. - x) * (1. + x));
pmm *= fact = 1.;
fact += for(i = 1; i <= m; i++) {
pmm *= -fact * somx2;
fact += 2.;
} }
} }
if(l == m) { return Cte * pmm; }
return
}
else { else {
pmmp1 = x * pmmp1 = x * (2 * m + 1) * pmm;
if(l == (m + 1)) { return Cte * pmmp1; }
return
}
else { else {
{ for(ll = (m + 2); ll <= l; ll++) {
pll = pll = (x * (2 * ll - 1) * pmmp1 - (ll + m - 1) * pmm) / (ll - m);
pmm = pmm = pmmp1;
pmmp1 = pmmp1 = pll;
} }
return Cte; return Cte * pll;
} }
} }
} }
void LegendreRecursive(int l, int m, double x, double P[]) void LegendreRecursive(int l, int m, double x, double P[])
{ {
/* Computes recursively a (l+1)-terms sequence of the associated /* Computes recursively a (l+1)-terms sequence of the associated
Legendre polynomial P_l^m(x). Legendre polynomial P_l^m(x).
l and m are the integers satisfying 0<=m<=l l and m are the integers satisfying 0<=m<=l
x lies in the range -1<=x<=1 x lies in the range -1<=x<=1
l = maximum order considered, m = invariable */ l = maximum order considered, m = invariable */
int int il;
double Pl_m, double Pl_m, Plm1_m;
P[0] = Plm1_m = Legendre(0, m, P[0] = Plm1_m = Legendre(0, m, x);
P[1] = Pl_m = Legendre(1, m, P[1] = Pl_m = Legendre(1, m, x);
if(l >= 2)
for(il = il < for(il = 1; il < l; il++) {
= + P[il + 1] = (2 * il + 1) * x * Pl_m / (il - m + 1) +
Plm1_m = (il + m) * Plm1_m / (m - il - 1);
Pl_m = Plm1_m = Pl_m;
Pl_m = P[il + 1];
} }
} }
void LegendreRecursiveM(int l, double x, double P[]) void LegendreRecursiveM(int l, double x, double P[])
{ {
/* Computes recursively a (l+1)-terms sequence of the associated /* Computes recursively a (l+1)-terms sequence of the associated
Legendre polynomial P_l^m(x). Legendre polynomial P_l^m(x).
x lies in the range -1<=x<=1, l = invariable, -l<=m<=l */ x lies in the range -1<=x<=1, l = invariable, -l<=m<=l */
int int m;
double Pl_m, double Pl_m, Plm1_m;
== 1.) if(fabs(x) == 1.)
for(m = m <= for(m = -l; m <= l; m++)
= ? : P[l + m] = (m == 0) ? pow(THESIGN(x), (double)l) : 0.;
else {
if(l == 0) {
P[0] = Legendre(0, 0, P[0] = Legendre(0, 0, x);
return; return;
} }
P[0] = Plm1_m = Legendre(l, -l, P[0] = Plm1_m = Legendre(l, -l, x);
P[1] = Pl_m = Legendre(l, P[1] = Pl_m = Legendre(l, -l + 1, x);
>= 1) if(l >= 1)
for(m = m < m for(m = -l + 1; m < l; m++) {
= + P[l + m + 1] = -2 * m * x * Pl_m / sqrt(1 - x * x) +
Plm1_m = (m * (m - 1) - l * (l + 1)) * Plm1_m;
Pl_m = Plm1_m = Pl_m;
Pl_m = P[l + m + 1];
} }
else else
return;
} }
} }
double dLegendre(int l, int m, double x) double dLegendre(int l, int m, double x)
{ {
/* Computes the derivative of the associated Legendre polynomial /* Computes the derivative of the associated Legendre polynomial
P_l^m(x) */ P_l^m(x) */
double dpl; double dpl;
if{ if(THEABS(m) > l || fabs(x) > 1.) {
Message::Error("Bad arguments for dLegendre: -l<=m<=l (integers), -1<=x<=1." Message::Error("Bad arguments for dLegendre: -l<=m<=l (integers), -1<=x<=1."
" Current values: l %d m %d x %.8g", l, m, " Current values: l %d m %d x %.8g",
l, m, x);
return 0.; return 0.;
} }
dpl = 0.; if(fabs(x) == 1.)
dpl = 0.;
else else
dpl = ? 0. : dpl = ((l + m) * (l - m + 1) * sqrt(1 - x * x) *
Legendre(l, x)) + ((THEABS((m - 1)) > l) ? 0. : Legendre(l, m - 1, x)) +
m * x * Legendre(l, m, x)) /
(1 - x * x);
return dpl; return dpl;
} }
double dLegendreFinDif(int l, int m, double x) double dLegendreFinDif(int l, int m, double x)
{ {
/* Computes the derivative of the associated Legendre polynomial /* Computes the derivative of the associated Legendre polynomial
P_l^m(x) using Finite Differences: f'(x) = (f(x+\delta P_l^m(x) using Finite Differences: f'(x) = (f(x+\delta
x)-f(x-\delta x))/(2 \delta) */ x)-f(x-\delta x))/(2 \delta) */
double dpl, delta = 1e-6; double dpl, delta = 1e-6;
if{ if(THEABS(m) > l || fabs(x) > 1.) {
Message::Error("Bad arguments for dLegendreFinDif: -l<=m<=l (integers), " Message::Error("Bad arguments for dLegendreFinDif: -l<=m<=l (integers), "
"-1<=x<=1. Current values: l %d m %d x %.8g", l, m, "-1<=x<=1. Current values: l %d m %d x %.8g",
l, m, x);
return 0.; return 0.;
} }
dpl = delta); dpl = (Legendre(l, m, x + delta) - Legendre(l, m, x - delta)) / (2 * delta);
return dpl; return dpl;
} }
void SphericalHarmonics(int l, int m, double Theta, double Phi, double Yl_m[]) void SphericalHarmonics(int l, int m, double Theta, double Phi, double Yl_m[])
{ {
int int am;
double cn, Pl_m, F, double cn, Pl_m, F, cRe;
cn = /* Normalization Factor */ cn = sqrt((2 * l + 1) * ONE_OVER_FOUR_PI); /* Normalization Factor */
am = am = THESIGN(m) * m;
F; F = sqrt(Factorial((double)(l - am)) / Factorial((double)(l + am)));
Pl_m = Legendre(l, am, cos(Theta)); Pl_m = Legendre(l, am, cos(Theta));
cRe = cn * F * Pl_m; cRe = cn * F * Pl_m;
Yl_m[0] = Yl_m[0] = cRe * cos(m * Phi);
Yl_m[1] = Yl_m[1] = cRe * sin(m * Phi);
} }
End of changes. 38 change blocks.
87 lines changed or deleted 95 lines changed or added