SphericalEngine.cpp File Reference

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

namespace  GeographicLib
 

Namespace for GeographicLib.



Detailed Description

Implementation for GeographicLib::SphericalEngine class.

Copyright (c) Charles Karney (2011) <charles@karney.com> and licensed under the MIT/X11 License. For more information, see http://geographiclib.sourceforge.net/

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:

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.

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines

Generated on 6 Oct 2014 for GeographicLib by  doxygen 1.6.1