{-# LANGUAGE FlexibleContexts, GADTs #-}
import Graphics.Dynamic.Plot.R2
import Math.LinearMap.Category
import Linear (V2(..))
import Linear.Affine (Point(..))
import Data.VectorSpace
import qualified Diagrams.Prelude as Dia
type Time = Double
type Δ x = x
rk4 :: (VectorSpace v, Scalar v ~ Time)
=> Δ Time -> (Time -> v -> Δ v) -> (Time,v) -> [(Time,v)]
rk4 h f (t,y) = (t,y) : rk4 h f (t+h, y ^+^ (h/6)*^(k₁ ^+^ 2*^k₂ ^+^ 2*^k₃ ^+^ k₄))
where k₁ = f t y
k₂ = f (t+h/2) (y ^+^ (h/2)*^k₁)
k₃ = f (t+h/2) (y ^+^ (h/2)*^k₂)
k₄ = f (t+h) (y ^+^ h*^k₃)
type Space = V2 Double
type Velo = Δ Space
type Mass = Double
type Accel = Δ Velo
gravity :: Mass -> Space -> Accel
gravity m r = (-m / magnitude r^3) *^ r
m₁, m₂ :: Mass
m₁ = 0.03
m₂ = 0.003
bodies :: [(Time, ((Space, Velo), (Space, Velo), (Space, Velo)))]
bodies = takeEach 20 $ rk4 0.005
(\_ ((rSun, vSun), (rEarth, vEarth), (rSat, vSat))
-> ( (vSun, gravity m₁ (rSun^-^rEarth) ^+^ gravity m₂ (rSun^-^rSat))
, (vEarth, gravity 1 (rEarth^-^rSun) ^+^ gravity m₂ (rEarth^-^rSat))
, (vSat, gravity 1 (rSat^-^rSun) ^+^ gravity m₁ (rSat^-^rEarth)) ) )
( 0, ( (V2 (-0.02) 0, V2 0 (-0.03))
, (V2 1 0, V2 0 1)
, (V2 0.499 0.885, V2 (-0.885) 0.5)) )
takeEach n l = h : takeEach n r
where (h:_,r) = splitAt n l
plotWindow [
plotLatest [ {-legendName (show t) . -} plotDelay 0.05
$ plotMultiple [ shapePlot . Dia.moveTo (P r) $ Dia.circle d
| (d,r) <- [(0.04,rEarth), (0.1,rSun), (0.02,rMoon)] ]
| (t, ((rSun,_), (rEarth,_), (rMoon,_))) <- bodies ]
, dynamicAxes
, xInterval (-1.4,1.4)
, yInterval (-1,1)
]
GraphWindowSpecR2{lBound=-1.8666666666666663, rBound=1.8666666666666663, bBound=-1.3333333333333335, tBound=1.3333333333333335, xResolution=640, yResolution=480}