In [1]:
{-# LANGUAGE TypeFamilies, FlexibleContexts #-}
import Math.LinearMap.Category
import Data.VectorSpace
import Linear.V3
import Data.AffineSpace
import Control.Arrow
import Data.Semigroup
import qualified Diagrams.Prelude as Dia


### Types for the physical quantities¶

In [2]:
type ℝ = Double
type Distance = ℝ  -- in m
type Pos = V3 Distance
type Speed = ℝ -- in m/s
type Velo = V3 Speed
type PhaseSpace = (Pos, Velo)
type Mass = ℝ   -- in kg


### Standard, 4.-Ordnung Runge-Kutta method¶

In [3]:
rk4 :: (AffineSpace y, RealSpace (Diff y), t ~ ℝ)
=> (y -> Diff y) -> Diff t -> (t,y) -> [(t,y)]
rk4 f h = go
where go (t,y) = (t,y) : go
(t+h, y .+^ h/6 · (k₁ ^+^ 2·k₂ ^+^ 2·k₃ ^+^ k₄))
where k₁ = f y
k₂ = f $y .+^ h/2 · k₁ k₃ = f$ y .+^ h/2 · k₂
k₄ = f $y .+^ h · k₃  In [4]: import Graphics.Dynamic.Plot.R2 import qualified Diagrams.Prelude as Dia import Data.List takeEach :: Int -> [a] -> [a] takeEach n = go n where go i (x:xs) | i>=n = x : go 1 xs | otherwise = go (i+1)$ xs

trajectoryPlot :: Int -> [(String, Distance)] -> [[(ℝ,ℝ)]] -> DynamicPlottable
trajectoryPlot speed meta = plotLatest
. map ( transpose . take 80 >>>
\chunkCompos -> {-plotDelay (1/5) $-} plotMultiple [ (if name/="" then legendName name else id)$ plot [ lineSegPlot chunkCompo
, shapePlot . Dia.moveTo (Dia.p2 $last chunkCompo) . Dia.opacity 0.6$ Dia.circle radius ]
| ((name,radius), chunkCompo) <- zip meta chunkCompos ]
)
. iterate (drop 4) . takeEach speed

In [5]:
type ThreeBody = ((PhaseSpace, PhaseSpace), PhaseSpace)

me, mm, ms :: Mass
me = 5.972e24
mm = 7.346e22
ms = 2e3

moonDist = 404e6 -- in m
moonRadius = 1.737e6 -- in m
earthRadius = 6.371e6 -- in m

moonSpeed :: Speed
moonSpeed = 1020 -- in m/s

gravConst :: ℝ
gravConst = 6.674e-11  -- in N·m²/kg²


### Trajectory of a 3-body system¶

Using the gravitational acceleration $$\vec F = G\cdot M\cdot m\cdot \frac{\vec{\Delta x}}{\|\vec{\Delta x}\|^3}$$

In [6]:
gravAcc :: Mass -> Diff Pos -> Diff Velo
gravAcc mt δx = (gravConst * mt / magnitude δx^3) *^ δx

traject3Body :: ThreeBody -> [ThreeBody]
traject3Body xv₀ = snd <$> rk4 (\(((xe,ve), (xm,vm)), (xs,vs)) -> (((ve, gravAcc mm$ xm.-.xe)
, (vm, gravAcc me $xe.-.xm)) , (vs, gravAcc me (xe.-.xs) ^+^ gravAcc mm (xm.-.xs)) ) ) 90 (0, xv₀)  In [7]: trebuchetOrbit :: Velo -> [ThreeBody] trebuchetOrbit dv = replicate 1000{-0-} startState ++ traject3Body startState where startState = ( ((zeroV,zeroV) ,(V3 moonDist 0 0, V3 0 moonSpeed 0)) , ( V3 (moonDist+ny*moonRadius) (-nx*moonRadius) 0, V3 0 moonSpeed 0 ^-^ dv ) ) V3 nx ny nz = dv ^/ magnitude dv  ### Geocentric view¶ In [8]: dv = V3 (-2275) 1200 0 :: Velo plotWindow [ trajectoryPlot 32 [("Earth",earthRadius), ("Moon",moonRadius), ("Spacecraft",1)] [ [(xe,ye),(xm,ym),(xs,ys)] | (((V3 xe ye _,_),(V3 xm ym _,_)),(V3 xs ys _,_)) <- trebuchetOrbit dv ] , dynamicAxes, unitAspect, forceXRange (-2*moonDist, 2*moonDist) , forceYRange (-moonDist, moonDist) ]  GraphWindowSpecR2{lBound=-6.597292373896027e8, rBound=6.597292373896027e8, bBound=-4.94796928042202e8, tBound=4.94796928042202e8, xResolution=640, yResolution=480} ### Lunacentric view¶ In [9]: plotWindow [ trajectoryPlot 4 [("Earth",earthRadius), ("Moon",moonRadius), ("Spacecraft",1)] [ [(xe-xm,ye-ym),(0,0),(xs-xm,ys-ym)] | (((V3 xe ye _,_),(V3 xm ym _,_)),(V3 xs ys _,_)) <- trebuchetOrbit dv ] , dynamicAxes, unitAspect, forceXRange (-20*moonRadius, 20*moonRadius) , forceYRange (-15*moonRadius, 15*moonRadius) ]  GraphWindowSpecR2{lBound=-1.2581314850581405e8, rBound=7.83201368639614e7, bBound=-4.8371011854243636e8, tBound=-3.3061026550688004e8, xResolution=640, yResolution=480} In [10]: import Text.Printf veloVis = 2e3 launch = [ Dia.p2 (x₀, y₀) , Dia.p2 (x₀ - v₀x*veloVis, y₀ - v₀y*veloVis) ] where x₀ = ny*moonRadius y₀ = -nx*moonRadius V3 v₀x v₀y _ = dv V3 nx ny nz = dv ^/ magnitude dv plotWindow [ legendName "To Earth" . shapePlot . Dia.dashingO [3,7] 0$ Dia.arrowBetween (Dia.p2(0,0)) (Dia.p2(-5e6,0))
, legendName "Moon"
(shapePlot $Dia.arrowBetween Dia.origin (Dia.p2 (0, moonSpeed*veloVis))) <> shapePlot (Dia.opacity 0.6$ Dia.circle moonRadius)
, legendName (printf "v₀ = %.0g m/s" $magnitude dv) . shapePlot$ Dia.arrowBetween (head launch) (last launch)
, unitAspect ]

GraphWindowSpecR2{lBound=-4635285.003931116, rBound=4995677.898723568, bBound=-3460111.0884955074, tBound=3763111.0884955074, xResolution=640, yResolution=480}

In [ ]: