## Introduction

The planet Mercury has a highly elliptical orbit with a perihelion of about 0.31 AU and an aphelion of about 0.47 AU. This ellipse is not stationary but itself rotates about the Sun, a phenomenon known as the precession of the perihelion. A calculation carried out using Newtonian mechanics gives a value at variance with observation. The deficit is explained using General Relativity although we do not apply the relativistic correction in this article.

Just to give a flavour of the Haskell, we will have to calculate values of the infinite series of Legendre Polynomials evaluated at 0. We have

Since we are dealing with infinite series we will want to define this co-recursively. We could use the Stream package but let us stay with lists.

```
> {-# OPTIONS_GHC -Wall #-}
> {-# OPTIONS_GHC -fno-warn-name-shadowing #-}
> {-# OPTIONS_GHC -fno-warn-type-defaults #-}
> {-# OPTIONS_GHC -fno-warn-unused-do-bind #-}
>
> {-# LANGUAGE NoMonomorphismRestriction #-}
>
> module Legendre (
> legendre0s
> , main) where
>
> import Data.List
> import Text.Printf
> import Initial
>
> legendre0s :: [Rational]
> legendre0s = interleave legendre0Evens legendre0Odds
> where
> legendre0Evens = 1 : zipWith f [1..] legendre0Evens
> where f n p = negate $ p * (2 * n * (2 * n - 1)) / (2^2 * n^2)
> legendre0Odds = 0 : legendre0Odds
>
> interleave :: [a] -> [a] -> [a]
> interleave = curry $ unfoldr g
> where
> g ([], _) = Nothing
> g (x:xs, ys) = Just (x, (ys, xs))
```

And now we can calculate any number of terms we need

```
ghci> take 10 $ legendre0s
[1 % 1,0 % 1,(-1) % 2,0 % 1,3 % 8,0 % 1,(-5) % 16,0 % 1,35 % 128,0 % 1]
```

The reader wishing to skip the physics and the mathematical derivation can go straight to the section on implementation

This article calculates the precession in Haskell using Newtonian methods. Over a long enough period, the gravitational effect of each outer planet on Mercury can be considered to be the same as a ring with the same mass as the planet; in other words we assume that the mass of each planet has been smeared out over its orbit. Probably one can model Saturn’s rings using this technique but that is certainly the subject of a different blog post.

More specifically, we model the mass of the ring as being totally concentrated on one particular value of and one particular value of with total mass .

where is the Dirac delta function. Thus the density of our ring is

## Acknowledgement

This blog follows the exposition given in [@Fitz:Newtonian:Dynamics] and [@brown:SpaceTime] concretized for the precession of the perihelion of Mercury with some of the elisions expanded. More details on Legendre Polynomials can be found in [@Bowles:Legendre:Polynomials].

## Axially Symmetric Mass Distributions

We consider axially symmetric mass distributions in spherical polar co-ordinates where runs from to , (the polar angle) runs from to and (the azimuthal angle) runs from to .

For clarity we give their conversion to cartesian co-ordinates.

The volume element in spherical polar co-ordinates is given by .

The gravitational potential given by masses each of mass and at position is:

If instead of point masses, we have a mass distribution then

where is the volume element.

If the mass distribution is axially symmetric then so will the potential. In spherical polar co-ordinates:

where denotes the average over the azimuthal angle.

Expanding the middle term on the right hand size and noting that :

Writing and noting that

where are the Legendre Polynomials we see that when

Applying the Spherical Harmonic Addition Theorem (or see [@arfken]) we obtain

Similarly when we obtain

Substituting into the equation for the potential for axially symmetric mass distributions gives us

where

Note that the first integral has limits to and the second has limits to .

It is well known that the Legendre Polynomials form an orthogonal and complete set for continuous functions. Indeed

Thus we can write

Using the orthogonality condition we have

Hence

## Gravitational Potential of a Ring

We now substitute in the axially symmetric density of a ring

Substituting again

Thus for

And for

Thus at and we have

and for

Let be the mass of the Sun then the potential due to all the Sun and all planets at a distance (excluding the planet positioned at ) is

## Apsidal Angles

An apsis is the closest and furthest point that a planet reaches i.e. the perihelion and the aphelion. Without the perturbing influence of the outer planets the angle between these points, the apsidal angle, would be . In presence of the outer planets this is no longer the case.

Writing down the Lagrangian for a single planet we have

where is the total potential due to the Sun and the other planets (as calculated above). is ignorable so we have a conserved quantity . We write for which is also conserved.

Applying Lagrange’s equation for we have

Thus the radial equation of motion is

To make further progress let us take just one term for a ring outside the planet of consideration and use the trick given in [@brown:SpaceTime]. Writing for the aphelion, for the perihelion and for the major radius we have

Defining by writing

we have

giving

Using the Taylor approximation

Thus

Then since

We have

It is a nuisance to be continually writing . From now on this is denoted by . Using

We obtain

We can therefore re-write the radial equation of motion approximately as

Now let us re-write the equation of motion as a relation between and .

Thus we have

Letting we can re-write this as

This is the equation for simple harmonic motion with and since for a circular orbit we can write

and therefore the change in radians per revolution is

To convert this to arc-seconds per century we apply a conversion factor

where 414.9 is the number of orbits of Mercury per century.

## Implementation

The implementation is almost trivial given that we have previously calculated the Legendre Polynomials (evaluated at 0). First let us make the code a bit easier to read by defining arithmetic pointwise (note that for polynomials we would not want to do this).

```
> instance Num a => Num [a] where
> (*) = zipWith (*)
> (+) = zipWith (+)
> abs = error "abs makes no sense for infinite series"
> signum = error "signum makes no sense for infinite series"
> fromInteger = error "fromInteger makes no sense for infinite series"
```

Next we define our conversion function so that we can compare our results against those obtained by Le Verrier.

```
> conv :: Floating a => a -> a
> conv x = x * 414.9 * (360 / (2 * pi)) * 3600
```

The main calculation for which we can take any number of terms.

```
> perturbations :: Double -> Double -> Double -> Double -> [Double]
> perturbations mRing mSun planetR ringR =
> map ((pi * (mRing / mSun)) *) xs
> where
> xs = (map (^2) $ map fromRational legendre0s) *
> (map fromIntegral [0..]) *
> (map fromIntegral [1..]) *
> (map ((planetR / ringR)^) [1..])
```

Arbitrarily, let us take 20 terms.

```
> predict :: Double -> Double -> Double -> Double
> predict x y z = sum $
> map conv $
> take 20 $
> perturbations x sunMass y z
```

And now let us compare our calculations with Le Verrier’s.

`> main :: IO () > main = do > printf "Venus %3.1f %3.1f\n" > (280.6 :: Double) > (predict venusMass mercuryMajRad venusMajRad) > printf "Earth %3.1f %3.1f\n" > (83.6 :: Double) > (predict earthMass mercuryMajRad earthMajRad) > printf "Mars %3.1f %3.1f\n" > (2.6 :: Double) > (predict marsMass mercuryMajRad marsMajRad) > printf "Jupiter %3.1f %3.1f\n" > (152.6 :: Double) > (predict jupiterMass mercuryMajRad jupiterMajRad)`

`ghci> main Venus 280.6 286.0 Earth 83.6 95.3 Mars 2.6 2.4 Jupiter 152.6 160.1`

Not too bad.

## Appendix

Note the lectures by Fitzpatrick [@Fitz:Newtonian:Dynamics] use a different approximation for the apsidal angle

We do not derive this here but note that the expansion and approximation are not entirely straightforward and are given here for completenes. Note that the result derived this way is identical to the result obtained in the main body of the article.

The radial force is given by

We also have

Thus

Re-arranging

we note that the last two terms can be re-written with a numerator of

and a denominator which is dominated by the . Thus

Since this term is we can expand the term of interest further

## Bibliography

Arfken, George. 1985. *Mathematical Methods for Physicists*. Third.. Orlando: ap.

Bowles, Robert. “Properties of Legendre Polynomials.” http://www.ucl.ac.uk/~ucahdrb/MATHM242/OutlineCD2.pdf.

Brown, Kevin. 2013. *Physics in space and time*. lulu.com.

Fitzpatrick, Richard. 1996. “Newtonian Dynamics.” http://farside.ph.utexas.edu/teaching/336k/lectures.

You could simplify and clarify your equations if you used geometric units, such that c = G = 1. This is customary, especially when working with general relativity. If you are ever reading a book on relativity and you think the equations don’t seem to be dimensionally correct, and are missing factors of G or c, if you look more closely you’ll see they are using geometric units.

And I believe in nuclear physics the situation is similar with with c = hbar = 1. But we are functional programmers and therefore insist on type correctness

Yes, but geometric units adhere perfectly to type correctness, as well as being simpler and more clear, which is why they are routinely used by both theoreticians and programmers.

I notice that you recently posted a book reivew on Amazon asserting that an equation written in geometric units was “dimensionally incorrect”. That is untrue. You can read about geometric units in any standard book on relativity. Equations written in geometrical units are dimensionally correct. (Note that mass, energy, distance, and time all have the same units on this basis.) It’s a shame that the book’s prospects for being read are now tarnished with that incorrect claim. It would be nice if you removed that incorrect claim from your review. After all, it’s well known that you functional programmers insist on correctness!

Apologies and amended. I was following Fitzpatricks’s notes where the gravitational constant is explicit and previously following Hairer’s book where again this is explicit. However I have just looked at my copy of O’Neill’s text on Semi-Riemannian Geometry and he does indeed use geometric units. I defer to your greater knowledge. It feels like one is throwing away type information (in the type theoretic sense) by using this system – perhaps the subject matter for another blog.

For more information on geometrical units (also called geometrized units or relativistic units), you could check out the standard relativity texts by Misner, Thorne, & Wheeler, or Wald, or D’Inverno, or Rindler, etc.

I much appreciate the amended wording of your book review – that was very gracious of you. I hesitate to mention it, but the amended wording, although a big improvement, is still not actually correct, because it says the reader “will find no mention of the gravitational constant” in the calculation of Mercury’s precession, whereas that calculation actually begins with the words “We are using units such that the gravitational constant and the speed of light are both unity.” Likewise each of the other sections of the book include reminders to the reader about the use of geometrical units. I suppose what you meant is that, since the book makes use of geometrical units (as the reader is repeatedly reminded), the constants G and c do not appear explicitly in the equations – as is quite common in the relativity literature. Hopefully that’s how people will interpret your comment.