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
Engineer's Journey 