Seems like every couple years I end up writing a program to calculate planetary orbits. This time I chose to use Haskell.

Paul Schlyter has written a really excellent article about the matter that is available on their homepage. In order to be able to follow along, it might be a good idea to read up about orbital elements on Wikipedia. I’ll try to explain them briefly as they’re encountered, but the Wikipedia article will still be helpful.

Since the calculations will concern various angles and coordinate systems and I know from experience that they’re tricky to keep in order, we’ll start by defining types to represent them. Following types are ones used to define orbit and they’re longitude of ascending node, inclination to ecliptic, argument of perihelion (or periapsis in general), semimajor axis, eccentricity and mean anomaly. All angles are in radians. Type parameters *a* and *b* are used to fill in which body is orbiting around which one. This way it’s easier to keep track of values and not accidentally mix them.

data LongAscNode a b = LongAscNode a b Float deriving Show data InclToEcl a b = InclToEcl a b Float deriving Show data ArgPeri a b = ArgPeri a b Float deriving Show data SemiMajor a b = SemiMajor a b Float deriving Show data Ecc a b = Ecc a b Float deriving Show data MeanAno a b = MeanAno a b Float deriving Show

Orbital elements depend on time. Trajectories of planets change slowly over time and for this we introduce way to calculate them. First representation of time and way to calculate how many days has passed since 31st of December, 1999 (taking leap years into account and factoring in time of day):

data Day d = Day Float deriving Show day :: Int -> Int -> Int -> Float -> Day Float day year month day time = Day $ fromIntegral d' + t' where d' = (367 * year - 7 * (year + (month + 9) `quot` 12 ) `quot` 4 + 275 * month `quot` 9 + day - 730530) t' = time / 24.0

Planets and sun are simple types. Advantage of having them separate like is that it’s possible to add new ones without modifying existing code. I was thinking if I should have made Bodies algebraic data type and have sun and planets be part of it, but decided against it. It might be interesting excercise to change this aspect and see how rest of the solution changes.

data Sun = Sun deriving Show data Earth = Earth deriving Show data Mercury = Mercury deriving Show

So, there needs to be a way of calculating orbital elements for given time. This can be done by defining type class *Orbit* with required functions and making instances for each planet:

class Orbit body center where n :: body -> center -> Day d -> LongAscNode body center i :: body -> center -> Day d -> InclToEcl body center w :: body -> center -> Day d -> ArgPeri body center a :: body -> center -> Day d -> SemiMajor body center e :: body -> center -> Day d -> Ecc body center m :: body -> center -> Day d -> MeanAno body center instance Orbit Sun Earth where n Sun Earth (Day _) = LongAscNode Sun Earth 0.0 i Sun Earth (Day _) = InclToEcl Sun Earth 0.0 w Sun Earth (Day d) = ArgPeri Sun Earth $ degToRad $ clamp $ 282.9404 + 4.70935E-5 * d a Sun Earth (Day _) = SemiMajor Sun Earth 1.0 e Sun Earth (Day d) = Ecc Sun Earth $ 0.016709 - 1.151E-9 * d m Sun Earth (Day d) = MeanAno Sun Earth $ degToRad $ clamp $ 356.0470 + 0.9856002585 * d instance Orbit Mercury Sun where n Mercury Sun (Day d) = LongAscNode Mercury Sun $ degToRad $ clamp $ 48.3313 + 3.24587E-5 * d i Mercury Sun (Day d) = InclToEcl Mercury Sun $ degToRad $ clamp $ 7.0047 + 5.00E-8 * d w Mercury Sun (Day d) = ArgPeri Mercury Sun $ degToRad $ clamp $ 29.1241 + 1.01444E-5 * d a Mercury Sun (Day _) = SemiMajor Mercury Sun 0.387098 e Mercury Sun (Day d) = Ecc Mercury Sun $ 0.205635 + 5.59E-10 * d m Mercury Sun (Day d) = MeanAno Mercury Sun $ degToRad $ clamp $ 168.6562 + 4.0923344368 * d

The tutorial I’m using has Sun orbiting around the Earth (presumably to simplify some calculations as you’re usually looking at things from Earth). This is the reason why we have *instance Orbit Sun Earth where* here.

As mentioned earlier, all angles are in radians and they’re converted back and forth with *degToRad* and *radToDeg* functions. Moreover, over time some of the values might go over 360 degrees or under 0 degrees and for those occasions we introduce *clamp* function that keeps angles within defined range:

radToDeg :: Floating a => a -> a radToDeg angle = 180 / pi * angle degToRad :: Floating a => a -> a degToRad angle = angle * pi / 180 clamp :: Float -> Float clamp angle | angle 360 = clamp $ angle - 360 | otherwise = angle

While it would be possible to use just these to check values, I thought that it would be nicer to have a system that tries to make sure that I don’t accidentally mix degrees and radians. Since all calculations are done in radians anyway, I decided to keep that as default. For viewing values in degrees, new type *Degrees a* is introduced. This can wrap any angle (or rather anything) and show it in degrees. By adding type class *Degreeable* (something that can be turned into degrees) and defining instances for elements we want to support conversion to degrees, we get nice system:

data Degrees a = Degrees a deriving Show class Degreeable element where toDeg :: element -> Degrees element instance Degreeable (ArgPeri a b) where toDeg (ArgPeri a b angle) = Degrees $ ArgPeri a b $ radToDeg angle instance Degreeable (LongAscNode a b) where toDeg (LongAscNode a b angle) = Degrees $ LongAscNode a b $ radToDeg angle instance Degreeable (InclToEcl a b) where toDeg (InclToEcl a b angle) = Degrees $ InclToEcl a b $ radToDeg angle instance Degreeable (MeanAno a b) where toDeg (MeanAno a b angle) = Degrees $ MeanAno a b $ radToDeg angle

With this one can write following to get Mercury’s longitude of ascending node on 29th of June, 2018 at noon:

*Planets> toDeg $ n Mercury Sun $ day 2018 6 29 12 Degrees (LongAscNode Mercury Sun 48.550575)

In order to find planet’s location in orbit, we need mean anomaly, eccentric anomaly and true anomaly. Mean anomaly we have in orbital elements, rest two needs to be calculated. And as always, we start by defining types to represent these:

data EccAnomaly a b = EccAnomaly a b Float deriving Show instance Degreeable (EccAnomaly a b) where toDeg (EccAnomaly a b angle) = Degrees $ EccAnomaly a b $ radToDeg angle data TrueAnomaly a b = TrueAnomaly a b Float deriving Show instance Degreeable (TrueAnomaly a b) where toDeg (TrueAnomaly a b angle) = Degrees $ TrueAnomaly a b $ radToDeg angle

Eccentric anomaly is calculated from mean anomaly and eccentricity, while true anomaly is calculated from eccentric anomaly and eccentricity (remember when I said that it’s tricky enough to keep various angles in order?):

eccAnomaly :: MeanAno a b -> Ecc a b -> EccAnomaly a b eccAnomaly (MeanAno a b m) (Ecc _ _ e) = EccAnomaly a b $ m + e * sin m * (1.0 + e * cos m) trueAnomaly :: EccAnomaly a b -> Ecc a b -> TrueAnomaly a b trueAnomaly (EccAnomaly a b ea) (Ecc _ _ e) = TrueAnomaly a b v where xv = cos ea - e yv = sqrt (1 - e ^ 2) * sin ea v = atan2 yv xv

When eccentric anomaly, eccentricity and semimajor axis are known, we can calculate distance of two objects:

data Distance a b = Distance a b Float deriving Show dist :: EccAnomaly a b -> Ecc a b -> SemiMajor a b -> Distance a b dist (EccAnomaly a b ea) (Ecc _ _ e) (SemiMajor _ _ semiMajorAxis) = Distance a b r where xv = semiMajorAxis * (cos ea - e) yv = semiMajorAxis * (sqrt (1 - e ^ 2) * sin ea) r = sqrt $ (xv ^ 2) + (yv ^ 2)

With true anomaly and distance we have solved location of the objects on the plane of orbit. Sometimes this is enough, but usually we want to translate these coordinates into ecliptic coordinate system. Here origin is set to location of Sun, x-axis runs towards vernal equinox. For more detailed information, Wikipedia has better explanation.

As always, we start by defining a new type to represent these coordinates. *EclCoord* has two parameters, *body* and *center*. This lets us encode what is orbiting what, which in turn makes it easier to translate between different coordinate systems (say, from Earth centered system to Mars centered system):

data EclCoord body center = EclCoord body center Float Float Float deriving Show toEclCoord :: TrueAnomaly a b -> Distance a b -> LongAscNode a b -> ArgPeri a b -> InclToEcl a b -> EclCoord a b toEclCoord (TrueAnomaly a b v) (Distance _ _ r) (LongAscNode _ _ n) (ArgPeri _ _ w) (InclToEcl _ _ i) = EclCoord a b x y z where x = r * (cos n * cos (v + w) - sin n * sin (v + w) * cos i) y = r * (sin n * cos (v + w) + cos n * sin (v + w) * cos i) z = r * (sin (v + w) * sin i) calcEclCoords :: Orbit a b => a -> b -> Day d -> EclCoord a b calcEclCoords body center day = toEclCoord ta di asc peri incl where ea = eccAnomaly (m body center day) ecc ta = trueAnomaly ea ecc ecc = e body center day di = dist ea ecc $ a body center day asc = n body center day peri = w body center day incl = i body center day

*calcEclCoords* is a convenience function that makes it easy to calculate ecliptic coordinates of given object.

Translating Sun centered coordinates to Earth centered system isn’t complicated, it’s matter of calculating relative location of Earth, Sun and the third object and then summing up the coordinates:

calcGeoCoord :: Orbit body Sun => body -> Day d -> EclCoord body Earth calcGeoCoord body day = EclCoord body Earth (x' + x) (y' + y) (z' + z) where (EclCoord _ _ x y z) = calcEclCoords body Sun day (EclCoord _ _ x' y' z') = calcEclCoords Sun Earth day

When viewing night sky, rectangular coordinates aren’t convenient. Instead, we rather translate them into right ascension and declination:

data EqCoord body = EqCoord body Float Float deriving Show instance Degreeable (EqCoord a) where toDeg (EqCoord a ra decl) = Degrees $ EqCoord a (radToDeg ra) $ radToDeg decl toEqCoordinates :: EclCoord body Earth -> Day Float -> EqCoord body toEqCoordinates (EclCoord a Earth x y z) (Day day) = EqCoord a ra decl where ra = atan2 ye xe decl = atan2 ze $ sqrt (xe^2 + ye^2) xe = x ye = y * cos ecl - z * sin ecl ze = y * sin ecl + z * cos ecl ecl = degToRad $ 23.4393 - 3.563E-7 * day calcEqCoords :: Orbit body Sun => body -> Day Float -> EqCoord body calcEqCoords body day = toEqCoordinates gCoords day where gCoords = calcGeoCoord body day

This concludes our calculations. There’s still (there always is) room for improvement. Calculations don’t take perturbations into account, which are needed in order to achiece better accuracy. Another nice addition would be calculating azimuth and altitude, which would make it even easier to locate the planet in the night sky.

I must say that for once there were not a single case where I accidentally had mixed coordinate systems or other values. Types kept track of them for me and I could concentrate on hunting typos in equations.

Whole sourcecode (for Sun, Mercury and Earth) is shown below.

{-# LANGUAGE MultiParamTypeClasses, FlexibleInstances, FlexibleContexts #-} module Planets where data Sun = Sun deriving Show data Earth = Earth deriving Show data Mercury = Mercury deriving Show data Day d = Day Float deriving Show data LongAscNode a b = LongAscNode a b Float deriving Show data InclToEcl a b = InclToEcl a b Float deriving Show data ArgPeri a b = ArgPeri a b Float deriving Show data SemiMajor a b = SemiMajor a b Float deriving Show data Ecc a b = Ecc a b Float deriving Show data MeanAno a b = MeanAno a b Float deriving Show class Orbit body center where n :: body -> center -> Day d -> LongAscNode body center i :: body -> center -> Day d -> InclToEcl body center w :: body -> center -> Day d -> ArgPeri body center a :: body -> center -> Day d -> SemiMajor body center e :: body -> center -> Day d -> Ecc body center m :: body -> center -> Day d -> MeanAno body center instance Orbit Sun Earth where n Sun Earth (Day _) = LongAscNode Sun Earth 0.0 i Sun Earth (Day _) = InclToEcl Sun Earth 0.0 w Sun Earth (Day d) = ArgPeri Sun Earth $ degToRad $ clamp $ 282.9404 + 4.70935E-5 * d a Sun Earth (Day _) = SemiMajor Sun Earth 1.0 e Sun Earth (Day d) = Ecc Sun Earth $ 0.016709 - 1.151E-9 * d m Sun Earth (Day d) = MeanAno Sun Earth $ degToRad $ clamp $ 356.0470 + 0.9856002585 * d instance Orbit Mercury Sun where n Mercury Sun (Day d) = LongAscNode Mercury Sun $ degToRad $ clamp $ 48.3313 + 3.24587E-5 * d i Mercury Sun (Day d) = InclToEcl Mercury Sun $ degToRad $ clamp $ 7.0047 + 5.00E-8 * d w Mercury Sun (Day d) = ArgPeri Mercury Sun $ degToRad $ clamp $ 29.1241 + 1.01444E-5 * d a Mercury Sun (Day _) = SemiMajor Mercury Sun 0.387098 e Mercury Sun (Day d) = Ecc Mercury Sun $ 0.205635 + 5.59E-10 * d m Mercury Sun (Day d) = MeanAno Mercury Sun $ degToRad $ clamp $ 168.6562 + 4.0923344368 * d day :: Int -> Int -> Int -> Float -> Day Float day year month day time = Day $ fromIntegral d' + t' where d' = (367 * year - 7 * (year + (month + 9) `quot` 12 ) `quot` 4 + 275 * month `quot` 9 + day - 730530) t' = time / 24.0 radToDeg :: Floating a => a -> a radToDeg angle = 180 / pi * angle degToRad :: Floating a => a -> a degToRad angle = angle * pi / 180 clamp :: Float -> Float clamp angle | angle 360 = clamp $ angle - 360 | otherwise = angle data Degrees a = Degrees a deriving Show class Degreeable element where toDeg :: element -> Degrees element instance Degreeable (ArgPeri a b) where toDeg (ArgPeri a b angle) = Degrees $ ArgPeri a b $ radToDeg angle instance Degreeable (LongAscNode a b) where toDeg (LongAscNode a b angle) = Degrees $ LongAscNode a b $ radToDeg angle instance Degreeable (InclToEcl a b) where toDeg (InclToEcl a b angle) = Degrees $ InclToEcl a b $ radToDeg angle instance Degreeable (MeanAno a b) where toDeg (MeanAno a b angle) = Degrees $ MeanAno a b $ radToDeg angle data EccAnomaly a b = EccAnomaly a b Float deriving Show eccAnomaly :: MeanAno a b -> Ecc a b -> EccAnomaly a b eccAnomaly (MeanAno a b m) (Ecc _ _ e) = EccAnomaly a b $ m + e * sin m * (1.0 + e * cos m) instance Degreeable (EccAnomaly a b) where toDeg (EccAnomaly a b angle) = Degrees $ EccAnomaly a b $ radToDeg angle data TrueAnomaly a b = TrueAnomaly a b Float deriving Show data Distance a b = Distance a b Float deriving Show dist :: EccAnomaly a b -> Ecc a b -> SemiMajor a b -> Distance a b dist (EccAnomaly a b ea) (Ecc _ _ e) (SemiMajor _ _ semiMajorAxis) = Distance a b r where xv = semiMajorAxis * (cos ea - e) yv = semiMajorAxis * (sqrt (1 - e ^ 2) * sin ea) r = sqrt $ (xv ^ 2) + (yv ^ 2) trueAnomaly :: EccAnomaly a b -> Ecc a b -> TrueAnomaly a b trueAnomaly (EccAnomaly a b ea) (Ecc _ _ e) = TrueAnomaly a b v where xv = cos ea - e yv = sqrt (1 - e ^ 2) * sin ea v = atan2 yv xv instance Degreeable (TrueAnomaly a b) where toDeg (TrueAnomaly a b angle) = Degrees $ TrueAnomaly a b $ radToDeg angle data EclCoord body center = EclCoord body center Float Float Float deriving Show toEclCoord :: TrueAnomaly a b -> Distance a b -> LongAscNode a b -> ArgPeri a b -> InclToEcl a b -> EclCoord a b toEclCoord (TrueAnomaly a b v) (Distance _ _ r) (LongAscNode _ _ n) (ArgPeri _ _ w) (InclToEcl _ _ i) = EclCoord a b x y z where x = r * (cos n * cos (v + w) - sin n * sin (v + w) * cos i) y = r * (sin n * cos (v + w) + cos n * sin (v + w) * cos i) z = r * (sin (v + w) * sin i) calcEclCoords :: Orbit a b => a -> b -> Day d -> EclCoord a b calcEclCoords body center day = toEclCoord ta di asc peri incl where ea = eccAnomaly (m body center day) ecc ta = trueAnomaly ea ecc ecc = e body center day di = dist ea ecc $ a body center day asc = n body center day peri = w body center day incl = i body center day calcGeoCoord :: Orbit body Sun => body -> Day d -> EclCoord body Earth calcGeoCoord body day = EclCoord body Earth (x' + x) (y' + y) (z' + z) where (EclCoord _ _ x y z) = calcEclCoords body Sun day (EclCoord _ _ x' y' z') = calcEclCoords Sun Earth day data EqCoord body = EqCoord body Float Float deriving Show instance Degreeable (EqCoord a) where toDeg (EqCoord a ra decl) = Degrees $ EqCoord a (radToDeg ra) $ radToDeg decl toEqCoordinates :: EclCoord body Earth -> Day Float -> EqCoord body toEqCoordinates (EclCoord a Earth x y z) (Day day) = EqCoord a ra decl where ra = atan2 ye xe decl = atan2 ze $ sqrt (xe^2 + ye^2) xe = x ye = y * cos ecl - z * sin ecl ze = y * sin ecl + z * cos ecl ecl = degToRad $ 23.4393 - 3.563E-7 * day calcEqCoords :: Orbit body Sun => body -> Day Float -> EqCoord body calcEqCoords body day = toEqCoordinates gCoords day where gCoords = calcGeoCoord body day