Implementation for GeographicLib::SphericalEngine class. More...
#include <GeographicLib/SphericalEngine.hpp>
#include <GeographicLib/CircularEngine.hpp>
#include <GeographicLib/Utility.hpp>
Go to the source code of this file.
Namespaces | |
GeographicLib | |
Namespace for GeographicLib. | |
Implementation for GeographicLib::SphericalEngine class.
Copyright (c) Charles Karney (2011-2017) charl and licensed under the MIT/X11 License. For more information, see es@k arney .comhttps://geographiclib.sourceforge.io/
The general sum is
V(r, theta, lambda) = sum(n = 0..N) sum(m = 0..n) q^(n+1) * (C[n,m] * cos(m*lambda) + S[n,m] * sin(m*lambda)) * P[n,m](t)
where t = cos(theta)
, q = a/r
. In addition write u = sin(theta)
.
P[n,m]
is a normalized associated Legendre function of degree n
and order m
. Here the formulas are given for full normalized functions (usually denoted Pbar
).
Rewrite outer sum
V(r, theta, lambda) = sum(m = 0..N) * P[m,m](t) * q^(m+1) * [Sc[m] * cos(m*lambda) + Ss[m] * sin(m*lambda)]
where the inner sums are
Sc[m] = sum(n = m..N) q^(n-m) * C[n,m] * P[n,m](t)/P[m,m](t) Ss[m] = sum(n = m..N) q^(n-m) * S[n,m] * P[n,m](t)/P[m,m](t)
Evaluate sums via Clenshaw method. The overall framework is similar to Deakin with the following changes:
cos(m*lambda)
and sin(m*lambda)
into the evaluation of the outer sum (rather than independently computing an array of these trigonometric terms).N
is large.For the general framework of Clenshaw, see http://mathworld.wolfram.com/ClenshawRecurrenceFormula.html
Let
S = sum(k = 0..N) c[k] * F[k](x) F[n+1](x) = alpha[n](x) * F[n](x) + beta[n](x) * F[n-1](x)
Evaluate S
with
y[N+2] = y[N+1] = 0 y[k] = alpha[k] * y[k+1] + beta[k+1] * y[k+2] + c[k] S = c[0] * F[0] + y[1] * F[1] + beta[1] * F[0] * y[2]
IF F[0](x) = 1
and beta(0,x) = 0
, then F[1](x) = alpha(0,x)
and we can continue the recursion for y[k]
until y[0]
, giving
S = y[0]
Evaluating the inner sum
l = n-m; n = l+m Sc[m] = sum(l = 0..N-m) C[l+m,m] * q^l * P[l+m,m](t)/P[m,m](t) F[l] = q^l * P[l+m,m](t)/P[m,m](t)
Holmes + Featherstone, Eq. (11), give
P[n,m] = sqrt((2*n-1)*(2*n+1)/((n-m)*(n+m))) * t * P[n-1,m] - sqrt((2*n+1)*(n+m-1)*(n-m-1)/((n-m)*(n+m)*(2*n-3))) * P[n-2,m]
thus
alpha[l] = t * q * sqrt(((2*n+1)*(2*n+3))/ ((n-m+1)*(n+m+1))) beta[l+1] = - q^2 * sqrt(((n-m+1)*(n+m+1)*(2*n+5))/ ((n-m+2)*(n+m+2)*(2*n+1)))
In this case, F[0] = 1
and beta[0] = 0
, so the Sc[m] = y[0]
.
Evaluating the outer sum
V = sum(m = 0..N) Sc[m] * q^(m+1) * cos(m*lambda) * P[m,m](t) + sum(m = 0..N) Ss[m] * q^(m+1) * cos(m*lambda) * P[m,m](t) F[m] = q^(m+1) * cos(m*lambda) * P[m,m](t) [or sin(m*lambda)]
Holmes + Featherstone, Eq. (13), give
P[m,m] = u * sqrt((2*m+1)/((m>1?2:1)*m)) * P[m-1,m-1]
also, we have
cos((m+1)*lambda) = 2*cos(lambda)*cos(m*lambda) - cos((m-1)*lambda)
thus
alpha[m] = 2*cos(lambda) * sqrt((2*m+3)/(2*(m+1))) * u * q = cos(lambda) * sqrt( 2*(2*m+3)/(m+1) ) * u * q beta[m+1] = -sqrt((2*m+3)*(2*m+5)/(4*(m+1)*(m+2))) * u^2 * q^2 * (m == 0 ? sqrt(2) : 1)
Thus
F[0] = q [or 0] F[1] = cos(lambda) * sqrt(3) * u * q^2 [or sin(lambda)] beta[1] = - sqrt(15/4) * u^2 * q^2
Here is how the various components of the gradient are computed
Differentiate wrt r
d q^(n+1) / dr = (-1/r) * (n+1) * q^(n+1)
so multiply C[n,m]
by n+1
in inner sum and multiply the sum by -1/r
.
Differentiate wrt lambda
d cos(m*lambda) = -m * sin(m*lambda) d sin(m*lambda) = m * cos(m*lambda)
so multiply terms by m
in outer sum and swap sine and cosine variables.
Differentiate wrt theta
dV/dtheta = V' = -u * dV/dt = -u * V'
here '
denotes differentiation wrt theta
.
d/dtheta (Sc[m] * P[m,m](t)) = Sc'[m] * P[m,m](t) + Sc[m] * P'[m,m](t)
Now P[m,m](t) = const * u^m
, so P'[m,m](t) = m * t/u * P[m,m](t)
, thus
d/dtheta (Sc[m] * P[m,m](t)) = (Sc'[m] + m * t/u * Sc[m]) * P[m,m](t)
Clenshaw recursion for Sc[m]
reads
y[k] = alpha[k] * y[k+1] + beta[k+1] * y[k+2] + c[k]
Substituting alpha[k] = const * t
, alpha'[k] = -u/t * alpha[k]
, beta'[k] = c'[k] = 0
gives
y'[k] = alpha[k] * y'[k+1] + beta[k+1] * y'[k+2] - u/t * alpha[k] * y[k+1]
Finally, given the derivatives of V
, we can compute the components of the gradient in spherical coordinates and transform the result into cartesian coordinates.
Definition in file SphericalEngine.cpp.