Legendre Functions, Polynomials, and Associated Legendre Polynomials
These polynomials which has used in physic, first discovered
by French mathematics professor, Adrein-Marie Legendre (1752-1833)[1]. Here we
will deal with their application in physical geodesy. Legendre was involved in
measuring distance between Paris Observatory and British Royal observatory by
triangulation method ( French-Anglo Survey).
Legendre Polynomial
A compact definition is given by Rodrigues’ formula
The following maxima script will give Legendre polynomials
up to n = 15.
s: 1 -
x^2;
N:15;
Pl:makelist(factor(diff(((-1)^n)*s^n,x,n)/((2^n)*n!)),n,0,N);
for
p:1 thru length(Pl) do display(Pl[p]);
plot2d([Pl[2],Pl[3],Pl[4],Pl[5],Pl[6],Pl[7]],[x,-1.0,1.0],[y,-1.5,1.5]);
The graph of P2 to P7 will be:
And 16 first polynomials are:
Another way to get these functions, is python:
from scipy.special import legendre
for n in range(10):
Pn =
legendre(n)
print(Pn)
the result will be:
P0 = 1
P1 = X
P2 = 1.5 x2 - 0.5
P3 = 2.5 x3 - 1.5 x
P4 = 4.375 x4 + 4.857e-16 x3
- 3.75 x2 + 2.429e-16 x + 0.375
P5 = 7.875 x5 - 8.75 x3
- 4.372e-16 x2 + 1.875 x
P6 = 14.44 x6 - 19.69 x4
+ 1.603e-15 x3 + 6.562 x2 - 0.3125
P7 = 26.81 x7 + 5.954e-15 x6
- 43.31 x5 + 2.977e-15 x4 + 19.69 x3 -
1.116e-15 x2 - 2.187 x
P8 = 50.27 x8 - 5.581e-15 x7
- 93.84 x6 - 2.233e-14 x5 + 54.14 x4 -
5.581e-15 x3 - 9.844 x2 - 4.361e-17 x + 0.2734
P9 = 94.96 x9 - 2.109e-14 x8
- 201.1 x7 + 6.326e-14 x6 + 140.8 x5 +
1.581e-14 x4 - 36.09 x3 - 1.977e-15 x2 + 2.461
x
Each of these functions can be generated by previous ones.
This property makes them suitable for computer programming by iteration,
although it can be done by recursion too, but recursion method is very slow. In
physical geodesy, to model Earth gravity, these functions are seen in solving
Laplace equation in spherical coordinates. The following is equation in book “Physical Geodesy” by Martin Vermeer.
These are computer
programs to get numerical values of Pn(x) for x between -1 and 1 (-1<=x<=1).
C# iterative
public static double
Legendre(int n, double x)
{
double LFS0 = 1.0;
double LFS1 = x;
double Lf = 0.0;
if(n == 0)
{
Lf = LFS0;
}
else if(n == 1)
{
Lf = LFS1;
}
else if (n >= 1)
{
for (int i = 2; i <=n; i++)
{
Lf =
(-(i - 1) * LFS0 + (2 * i - 1) * x * LFS1) / i;
LFS0 =
LFS1;
LFS1 =
Lf;
}
}
return Lf;
}
Python iterative
def
Legendre(x,n):
LFS0 = 1.0
LFS1 = x
Legen = 0.0
if(n==0) :
return LFS0
elif(n==1) :
return LFS1
elif(n>1):
for I in range(2,n+1):
Legen = (-(I-1)*LFS0 + (2*I-1)*x*LFS1)/(I)
LFS0 = LFS1
LFS1 = Legen
return Legen
Fortran iterative
function
Legendre(n,x)
implicit none
integer :: n, I
real(kind=8) :: x , Legendre
real(kind=8),dimension(0:1) :: LFS
LFS(0) = 1.0
LFS(1) = x
if(n==0) then
Legendre = LFS(0)
elseif(n==1) then
Legendre = LFS(1)
elseif(n>=2) then
do I = 2,n
Legendre = (-(I-1)*LFS(0) +
(2*I-1)*x*LFS(1))/I
LFS(0) = LFS(1)
LFS(1) = Legendre
end do
end if
end function
The following are recursive procedures, which are slow
compared to iterative ones , but code is more compact.
C# recursive
public static double LegendreR(int n, double x)
{
double Pn = 0.0;
if (n == 0)
{
return 1.0;
}
else if (n == 1)
{
return x;
}
else
{
Pn =
(-(n-1)*LegendreR(n-2,x) + (2*n-1)*x*LegendreR(n-1,x))/n;
}
return Pn;
}
Python recursive
def
legendR(x,n) :
if(n==0):
return 1.0
elif(n==1):
return x
else:
return (-(n-1)*legendf(x,n-2) +
(2*n-1)*x*legendf(x,n-1))/n
Fortran recursive
recursive
function Legendref(n,x) result(Pn)
integer :: n
real(kind=8):: Pn , x
if (n == 0) then
Pn = 1.0
else if (n == 1) then
Pn = x
else
Pn = (-(n-1)*Legendref(n-2,x) +
(2*n-1)*x*Legendref(n-1,x))/n
end if
end function Legendref
Associated Legendre Polynomials
In compact form:
Maxima script will give equations
s:1 - x^2;
N:15;
Pl:makelist(factor(diff(((-1)^n)*s^n,x,n)/((2^n)*n!)),n,0,N);
Plm:makelist(makelist(factor(((-1)^m)*(s^(m/2))*diff(Pl[i],x,m)),m,0,i-1),i,length(Pl));
for i:1 thru length(Plm) do
for j:1 thru
length(Plm[i]) do display(Plm[i][j]) ;
plot2d([Plm[4][1],Plm[4][2],Plm[4][3],Plm[4][4],Plm[4][4]],[x,-1.0,1.0],[y,-15,15]);
plot2d([Plm[6][2],Plm[5][2],Plm[4][2],Plm[3][2],Plm[2][2]],[x,-1.0,1.0],[y,-3.5,3]);
plot2d([Plm[7][3],Plm[6][3],Plm[5][3],Plm[4][3],Plm[3][3]],[x,-1.0,1.0],[y,-15,22]);
Associated Legendre polynomials up to Pl4 shown here.
By increasing l and m Associated Legendre
polynomials, become very large and unstable, therefore normalized associated
polynomials are used. Fully normalized Legendre polynomials
The above relations are not useful for computer programming,
because large number of factorials, and Plm instead the
following iterative, equations are used.
Numerical values of series of polynomials for every degree(l).
C#
public static double[]
FulNormALegenSerie(int n, double x)
{
double[] Pnm = new double[n + 1];
double[] BLFS = new double[n + 1];
double[] ALFS = new double[n + 1];
double s = Math.Sqrt(1.0 - x * x);
if (n == 0)
{
Pnm[0] = 1.0;
return Pnm;
}
if (n == 1)
{
Pnm[0] = x *
Math.Sqrt(3.0);
Pnm[1] = s *
Math.Sqrt(3.0);
return Pnm;
}
else if(n == 2)
{
Pnm[0] =
Math.Sqrt(5.0) * LegendreF(2, x);
Pnm[1] = 3.0
* s * x * Math.Sqrt(10.0 / 6.0);
Pnm[2] = 3.0
* (1.0 - x * x) * Math.Sqrt(10.0 / 24.0);
return Pnm;
}
else
{
BLFS[0] =
x * Math.Sqrt(3.0);
BLFS[1] =
s * Math.Sqrt(3.0);
ALFS[0] =
Math.Sqrt(5.0) * LegendreF(2, x);
ALFS[1] = 3.0 * s * x * Math.Sqrt(10.0 /
6.0);
ALFS[2] =
3.0 * (1.0 - x * x) * Math.Sqrt(10.0 / 24.0);
Pnm[0] =
Math.Sqrt(2.0 * n + 1.0) * LegendreF(n, x);
int l = 3;
while (l <= n)
{
int k = 1;
while (k <= l)
{
if (k < l - 1)
{
double Anm = Math.Sqrt((2.0 * l - 1.0) *
(2.0 * l + 1.0)/
((l - k) * (l + k)));
double Bnm = -Math.Sqrt(((l - 1.0) * (l -
1.0) - k * k)
/ (4.0 * (l - 1.0) * (l - 1.0) - 1.0));
Pnm[k] = Anm * (x * ALFS[k] + Bnm * BLFS[k]);
}
else
if(k >= l - 1)
{
if (k == l - 1)
{
Pnm[k]= Math.Sqrt(2.0 * k + 3.0) * x * ALFS[l - 1];
}
else if(k == l)
{
Pnm[k] = s * Math.Sqrt(1.0 + 1.0 / (2.0 * k)) * ALFS[l - 1];
}
}
k =
k + 1;
}
Array.Copy(ALFS, BLFS, ALFS.Length);
Array.Copy(Pnm, ALFS, Pnm.Length);
l = l +
1;
}
return Pnm;
}
Phyton
def
FulNormALegenSerie(x,n):
s = math.sqrt(1.0-x*x)
Pnm = []
BLFS = []
ALFS = []
BLFS.append(x*math.sqrt(3.0))
BLFS.append(s*math.sqrt(3.0))
ALFS.append(math.sqrt(5.0)*Legendre(x,2))
ALFS.append(3.0*s*x*math.sqrt(10.0/6.0))
ALFS.append(3.0*(1.0-x*x)*math.sqrt(10.0/24.0))
if(n==0) :
Pnm.append(1.0)
elif(n == 1) :
Pnm = BLFS
elif(n==2):
Pnm= ALFS
elif(n>2):
l = 3
while(l<= n):
Pnm = []
Pnm.append(math.sqrt(2.0*(l)+1.0)*Legendre(x,l))
k = 1
while(k <= l):
if(k < l-1) :
Anm =
math.sqrt(((2*l-1)*(2*l+1))/((l-k)*(l+k)))
Bnm =
-math.sqrt((((l-1)*(l-1)-k*k))/(4*(l-1)*(l-1)-1))
Pnm.append(Anm*(x*ALFS[k] + Bnm*BLFS[k]))
elif(k >= l-1) :
if(k == l-1):
Pnm.append(math.sqrt(2.0*(k)+3.0)*x*ALFS[l-1])
elif(k == l)
:
Pnm.append(s*math.sqrt(1.0+1.0/(2.0*k))*ALFS[l-1])
k = k + 1
BLFS = ALFS
ALFS = Pnm
l = l + 1
return Pnm
Fortran
subroutine
FulNormALegenSerie(n,x,Pnm)
implicit none
integer :: n , l , k
real(kind=8) :: x , s , Anm , Bnm
real(kind=8) :: Legendre, ld
real(kind=8),dimension(0:n+1) :: ALFS,BLFS,Pnm
s = dsqrt(1.0 - x*x)
Pnm = 0
BLFS(1) = s*dsqrt(3.D0)
ALFS(1) = 3.D0*s*x*dsqrt(5.D0/3.D0)
ALFS(2) =
3.D0*(1.D0-x*x)*dsqrt(5.D0/12.D0)
if(n == 1) then
Pnm = BLFS
elseif(n==2) then
Pnm = ALFS
else
l = 3
do while(l<= n)
k = 1
do while(k <= l)
if(k < l-1) then
Anm =
dsqrt(real((2*l-1)*(2*l+1),kind=8) &
/real((l-k)*(l+k),kind=8))
Bnm =
-dsqrt(real(((l-1)*(l-1)-k*k),kind=8) &
/real(4*(l-1)*(l-1)-1,kind=8))
Pnm(k) =
Anm*(x*ALFS(K) + Bnm*BLFS(K))
elseif(k >= l-1)
then
if(k == l-1) then
Pnm(k) =
dsqrt(2.D0*real(k,kind=8)+3.D0)*x*ALFS(l-1)
elseif(k == l)
then
ld = real(k,kind=8)
Pnm(k) =
s*dsqrt(1.D0+1.D0/(2.D0*ld))*ALFS(l-1)
endif
end if
k = k + 1
end do
BLFS = ALFS
ALFS = Pnm
l = l + 1
end do
end if
Pnm(0) =
dsqrt(2D0*real(n,kind=8)+1.D0)*Legendre(n,x)
end
subroutine
Testing
Normal potential of the Erath gravity field is[2]:
By derivation by r in spherical coordinates will give
approximate normal gravity.
GRS 80 constants
a = 6378137 m semi-major axis
GM0 = 3.98600.1014 geocentric
gravitational constant
w = 7.292115.10-5 angular velocity
In reference[2] r is:
In Geometric Geodesy [3]
We will set up five equation to find Cl,0
from reference or calculated from another normal equation, to check validity of
Legendre polynomials.
Calculated coefficients are:
C
|
0
|
2
|
4
|
6
|
8
|
r0
|
0.999993551446001
|
-0.00108416191395218
|
0.000004571732568181
|
-0.000000017043632669
|
0.000000000081735028
|
r1
|
1.00000556957139
|
-0.00108129896967865
|
0.000000471322575981
|
-0.000000000922562314
|
0.000000000010459033
|
Inserting coefficients in above gravity formula, we will get
gravity at latitudes.
Table of normal gravity at different latitudes in mGal
f
|
g0
|
g1
|
g adapted
|
g normal
|
0.00
|
978032.67715000
|
978032.67715000
|
978032.67715000
|
978032.67715000
|
10.00
|
978188.38359439
|
978188.38359288
|
978188.38360804
|
978188.38360804
|
20.00
|
978636.95381125
|
978636.95380894
|
978636.95383748
|
978636.95383748
|
30.00
|
979324.87035723
|
979324.87035723
|
979324.87035723
|
979324.87035723
|
40.00
|
980169.82967808
|
980169.82967888
|
980169.82963497
|
980169.82963497
|
50.00
|
981070.35687238
|
981070.35687157
|
981070.35682771
|
981070.35682771
|
60.00
|
981917.83849727
|
981917.83849727
|
981917.83849727
|
981917.83849727
|
70.00
|
982609.61953273
|
982609.61953507
|
982609.61956225
|
982609.61956225
|
80.00
|
983061.58822013
|
983061.58822166
|
983061.58823246
|
983061.58823246
|
90.00
|
983218.63685000
|
983218.63685000
|
983218.63684389
|
983218.63684389
|
For testing fully
normalized associated Legendre polynomials, we have used gravity anomaly of
following formula from “Definition
of Functions of the Geopotential and Their Calculation from Spherical Harmonics
Models”[4].
R = a and Plm is fully normalized Legendre
polynomial. Coefficients C and S are taken from file “ dV_ELL_Earth2014.gfc”[5]'' . Whether this assumption is correct or not, the results are
not as perfect as expected.
Gravity at Longitude and Latitude in mGal
Longitude
|
Latitude
|
g from file
|
calculated g
|
46
|
45
|
980915.91755
|
980915.05329
|
47
|
45
|
980908.90583
|
980911.38617
|
48
|
45
|
980909.65238
|
980911.22103
|
Software
All programs are in GitHub repositoryhttps://github.com/AminKH/Legendre-Polynomials
References
Identites and properties for associated Legendre functionsAssociated Legendre Polynomials and Spherical Harmonics Computation for Chemistry Applications
Taweetham Limpanuparb, Josh Milthorpey October 8, 2014
Geodetic Reference System 1980 by H. Moritz
http://mitgcm.org/~mlosch/geoidcookbook/node11.html
[1]
Adrien-Marie Legendre https://en.wikipedia.org/wiki/Adrien-Marie_Legendre
[2]
Physical Geodesy, Nico Sneeuw, Institute of Geodesy Uinversitat Stuttgart
[3]
Geometric Geodesy, Richard Rapp, The Ohio State University April 1991
No comments:
Post a Comment