The F447.f program is designed to calculate the geoid undulation for EGM96. I'm trying to understand how it accomplish this calculation. I stuck in SUBROUTINE LEGFDN. I want to understand what is the formula for RLNN(N) vector. So far I'm sure only for RLNN(1)=1.

I believe RLNN(2) = SITHET*DRTS(3) = sen(theta)*sqrt(3), but I don't understanding RLNN(3) and forward.

The loop below can be translated

      DO 15  N1 = 3,M1
      N = N1-1
      N2 = 2*N
   15 RLNN(N1) = DRTS(N2+1)*DIRT(N2)*SITHET*RLNN(N1-1)

to something like this

      DO 15  N = 2,M
  15  END DO

where M=360 is the polynomial degree and THETA (I believe) is the latitude, someone recognize this recursion formula as part of some Legendre polynomial?

  • 1
    $\begingroup$ I am happy to see FORTRAN code! All the GO TOs, not so much. $\endgroup$ May 22 '19 at 19:17
  • $\begingroup$ Are those non-standard Legendre polynomials? $\endgroup$ May 22 '19 at 20:48
  • 1
    $\begingroup$ ANCIENT, ALL CAPS FORTRAN code. It's nigh impossible to read. LAWYERS USE ALL CAPS TO RENDER KEY PARTS OF A LEGAL DOCUMENT UNREADABLE. There is no reason to emulate such behavior. $\endgroup$ May 23 '19 at 12:11
  • 1
    $\begingroup$ @DavidHammen what are these "lower case characters" of which you speak? :) live.staticflickr.com/195/522484762_a2cc146805_b.jpg $\endgroup$ May 23 '19 at 15:11

This is not a recurrence for a Legendre polynomial $P_{n}(\mu)$ but an Associated Legendre Function (ALF). This1 states much more than is needed here.

ALFs are often denoted $P_{n,m}(\mu)$ in geophysics. The degree $n$ and order $m$ are non-negative integers $0, 1, 2, \dots$ with $0 \le m \le n.$ Argument $\mu$ is the sine of the latitude (measured north or south from the equator), or equivalently the cosine of the colatitude (measured southward from the north pole).

ALFs appear as factors in solid spherical harmonic formulas. Such formulas are used to represent a gravitational field. Spherical harmonics are classified as zonal, sectorial, or tesseral according to whether $m=0$ (zonal), $m=n$ (sectorial), or $0 < m < n$ (tesseral). Applying those conventions to ALFs, zonals are the same as Legendre polynomials: $P_{n,0}(\mu) = P_{n}(\mu).$ Sectorials correspond to $P_{n,n}(\mu),$ while $P_{n,m}(\mu), \, 0 < m < n$ correspond to tesserals.

An ALF might evaluate to a very large numeric value for large $n$ or $m$. Geophysicists often introduce a normalization factor which I denote $N_{n,m}$. Its formula is \begin{align*} N_{n,m} = \begin{cases} \sqrt {(2n+1) \dfrac {(n-m)!}{(n+m)!}}, & m = 0 \\ \sqrt {2 (2n+1) \dfrac {(n-m)!}{(n+m)!}}, & m > 0. \end{cases} \end{align*} A "fully normalized" ALF using $N_{n,m}$ is indicated by a bar above the $P$: \begin{align*} \overline{P}_{n,m}(\mu) = N_{n,m} \, P_{n,m}(\mu) \quad \Leftrightarrow \quad P_{n,m}(\mu) = \overline{P}_{n,m}(\mu) / N_{n,m} . \end{align*}

A sectorial recurrence using colatitude $\theta$ is \begin{align*} P_{n,n}(\cos \theta) = (2n-1) \sin \theta \; P_{n-1,n-1}(\cos \theta) \end{align*} beginning with $P_{0,0} = 1$ by definition and $P_{1,1} = \sin \theta$ by the recurrence. (Note that the ALF argument is $\cos \theta$ but that $\sin \theta$ appears as a factor in the RHS.) This recurrence can be expressed using fully normalized ALFs by substituting $P_{n,n}(\mu) = \overline{P}_{n,n}(\mu) / N_{n,n}$. Assume $n > 1$: \begin{align*} \frac {\overline{P}_{n,n}(\cos \theta)} {N_{n,n}} & = (2n-1) \sin \theta \; \frac {\overline{P}_{n-1,n-1}(\cos \theta)} {N_{n-1,n-1}} \\ \frac {\overline{P}_{n,n}(\cos \theta)} {\sqrt{2 (2n+1) \dfrac {0!} {(2n)!}}} & = (2n-1) \sin \theta \; \frac {\overline{P}_{n-1,n-1}(\cos \theta)} {\sqrt{2 (2n-1) \dfrac {0!} {(2n-2)!}}} \\ \overline{P}_{n,n}(\cos \theta) & = \sqrt{\dfrac {2 (2n+1) } {(2n)(2n-1)(2n-2)!}} (2n-1) \sqrt {\frac {(2n-2)!} {2(2n-1)} } \sin \theta \; \overline{P}_{n-1,n-1}(\cos \theta) \\ & = \sqrt{ \dfrac {2n+1} {2n}} \sin \theta \; \overline{P}_{n-1,n-1}(\cos \theta) . \end{align*} This corresponds to the OP's code fragment, with the convention that $\overline{P}_{n,n}$ is stored in RLNN(N+1).

  • $\begingroup$ Exemplary answer! It's great when someone takes the time to wade through some complicated math and bring it to a successful resolution, and do some helpful teaching in the process. +1 Any thoughts on How can I verify my reconstructed gravity field of Ceres from spherical harmonics? they would be $\endgroup$
    – uhoh
    May 26 '19 at 21:16
  • $\begingroup$ So the ALF are simplified forms of Legendre polynomials and can be used instead of them in the appropriated framework? $\endgroup$ May 27 '19 at 14:43
  • $\begingroup$ I wouldn't say ALFs are simplified; if anything, they are more complicated because they are more general. Legendre polynomials (the m=0 case of ALFs) are associated only with zonal terms. ALFs are associated with sectorial and tesseral terms. All three sets are needed to represent a gravitational field using solid spherical harmonics. $\endgroup$
    – Bob Werner
    May 27 '19 at 15:56
  • $\begingroup$ @uhoh I agree with ChrisR's answer to your linked question -- validation is tough without someone else's results. What I could show is recurrences for Pnm, fullly normalized and not. Would that help? $\endgroup$
    – Bob Werner
    May 27 '19 at 23:23
  • $\begingroup$ @BobWerner Hmm... At the time, the only "someone else's result" I'd found was the Bouguer anomaly map I'd included in the question. That my maps truncated at n=4 and 5 don't look even a little bit like it may mean simply that I didn't know what I was doing. Why don't I have another look at each 1) my question 2) the literature and 3) ChrisR's answer first. In the mean time, I wonder if "recurrences for Pnm, fullly normalized and not" might be helpful for How are the coefficients in the EGM96 model normalized? $\endgroup$
    – uhoh
    May 27 '19 at 23:49

(This answer continues the one above.)

A table of $P_{n,m}(\mu)$ values can be envisioned as a square array. Rows down the page are indexed by degree $n$, and columns to the right by order $m$.

Due to the underlying math, the only sensible $P_{n,m}(\mu)$ values occur for $0 \le m \le n$. The effect is that the table is lower-triangular instead of square; $P_{n,m}(\mu)$ do not exist above the $m=n$ diagonal.

Values along the $m=n$ diagonal correspond to the sectorial terms described in the preceding answer. A simple recurrence provides the value of $P_{n,n}(\mu)$ from the preceding value $P_{n-1,n-1}(\mu)$.

A different recurrence provides values down a column. This "vertical" recurrence uses the two immediately preceding values higher in the column. For example in column $m=3$, the value of $P_{6,3}(\mu)$ can be computed from $P_{5,3}(\mu)$ and $P_{4,3}(\mu)$.

The vertical recurrence using colatitude $\theta$ is: \begin{align*} (n-m) P_{n,m}(\cos\theta) = (2n-1) \cos\theta \; P_{n-1,m} (\sin\theta) - (n-m-1) P_{n-2,m}(\cos\theta) \end{align*} This recurrence using unnormalized $P_{n,m}(\cos\theta)$ can be converted to fully normalized $\overline{P}_{n,m}(\cos\theta)$ in the same way as used above for converting the sectorial recurrence: $\overline{P}_{n,m}(\cos\theta) / N_{n,m}$ is substituted for $P_{n,m}(\cos\theta)$. \begin{align*} \overline{P}_{n,m}(\cos\theta) & = \frac {N_{n,m}} {n-m} \left[ \frac {2n-1} {N_{n-1,m}} \cos\theta \; \overline{P}_{n-1,m} (\cos\theta) - \frac {n+m-1} {N_{n-2,m}} \overline{P}_{n-2,m}(\cos\theta) \right] \\ & = \frac {\sqrt{(2-\delta_{0,m}) (2n+1) \frac{(n-m)!}{(n+m)!}}} {n-m} \left[ \begin{array}{l} \dfrac {2n-1} {\sqrt{(2-\delta_{0,m}) (2n-1) \frac{(n-m-1)!} {(n+m-1)!}}} \cos\theta \; \overline{P}_{n-1,m} (\cos\theta) \\ \quad - \dfrac {n+m-1} {\sqrt{(2-\delta_{0,m}) (2n-3) \frac{(n-m-2)!} {(n+m-2)!}}} \overline{P}_{n-2,m}(\cos\theta) \end{array} \right] \\ & = \frac {\sqrt{(2n+1) \frac{(n-m)(n-m-1)(n-m-2)!}{(n+m)(n+m-1)(n+m-2)!}}} {n-m} \left[ \begin{array}{l} \dfrac {2n-1} {\sqrt{(2n-1) \frac{(n-m-1)(n-m-2)!} {(n+m-1)(n+m-2)!}}} \cos\theta \; \overline{P}_{n-1,m} (\cos\theta) \\ \quad - \dfrac {n+m-1} {\sqrt{(2n-3) \frac{(n-m-2)!} {(n+m-2)!}}} \overline{P}_{n-2,m}(\cos\theta) \end{array} \right] \\ & = \frac {\sqrt{(2n+1) \frac{(n-m)(n-m-1)}{(n+m)(n+m-1)}}} {n-m} \left[ \begin{array}{l} \dfrac {2n-1} {\sqrt{(2n-1) \frac{(n-m-1)} {(n+m-1)}}} \cos\theta \; \overline{P}_{n-1,m} (\cos\theta) \\ \quad - \dfrac {n+m-1} {\sqrt{2n-3}} \overline{P}_{n-2,m}(\cos\theta) \end{array} \right] \\ & = \sqrt{\frac{2n+1} {(n-m)(n+m)}} \sqrt{\frac{n-m-1}{n+m-1}} \left[ \begin{array}{l} \sqrt{(2n-1) \frac{n+m-1} {n-m-1}} \cos\theta \; \overline{P}_{n-1,m} (\cos\theta) \\ \quad - \dfrac {n+m-1} {\sqrt{2n-3}} \overline{P}_{n-2,m}(\cos\theta) \end{array} \right] \\ & = \sqrt{\frac{2n+1} {(n-m)(n+m)}} \left[ \begin{array}{l} \sqrt{2n-1} \cos\theta \; \overline{P}_{n-1,m} (\cos\theta) \\ \quad - \sqrt{ \dfrac {(n-m-1)(n+m-1)} {2n-3} } \overline{P}_{n-2,m}(\cos\theta) \end{array} \right] . \end{align*} This is the vertical recurrence for fully normalized $\overline{P}_{n,m}(\cos\theta)$.

One wrinkle is that when computing the subdiagonal term $\overline{P}_{n+1,n}(\cos\theta)$, there is only one preceding term (the diagonal $\overline{P}_{n,n}(\cos\theta)$). The second preceding term (above the diagonal) does not exist. In this special case, only the first term within the brackets is needed.

F447.f program

The subroutine LEGFDN saves values of $\sqrt {i}$ in array DRTS(i) and the reciprocal $1 / \sqrt {i}$ in DIRT(i). The subroutine's code

      RLEG(N1) = DRTS(N2+1)*DIRT(N+M)*DIRT(N-M)*(DRTS(N2-1)*COTHET*  
     2 RLEG(N1-1)-DRTS(N+M-1)*DRTS(N-M-1)*DIRT(N2-3)*RLEG(N1-2))

near the bottom of DO-loop 30 is equivalent to

      RLEG(N1) = SQRT(2*N+1)/SQRT(N+M)/SQRT(N-M)*
     1 ( SQRT(2*N-1)*COS(THETA)*RLEG(N1-1)
     2 - SQRT(N+M-1)*SQRT(N-M-1)/SQRT(2*N-3)*RLEG(N1-2)
     3 )

This corresponds exactly to the vertical recurrence. Evidently $\overline{P}_{n,m}(\cos\theta)$ is stored in one-dimensional array RLEG(n). (I've not studied the code to understand how it handles the various columns.)

One conclusion is that the normalization formula used in the F447.f program is the $N_{n,m}$ formula given in the preceding answer.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.