Planetary Orbits and Haskell

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
Advertisements

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out /  Change )

Google+ photo

You are commenting using your Google+ account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s