Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
{-# 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
import Graphics.Dynamic.Plot.R2
import qualified Diagrams.Prelude as Dia
import Data.List
import Text.Printf
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
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₃
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
type ThreeBody = ((PhaseSpace, PhaseSpace), PhaseSpace)
me, mm, ms :: Mass
me = 5.972e24
mm = 7.346e22
ms = 2e3
moonDist, moonRadius, earthRadius :: Distance
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²
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₀)
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
dv = V3 (-2275) 1200 0 :: Velo
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
main = do
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) ]
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) ]
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 ]
return ()
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
@leftaroundabout

This comment has been minimized.

Copy link
Owner Author

leftaroundabout commented Sep 4, 2018

To run the .hs file:

  1. Install the Haskell Platform
  2. cabal install dynamic-plot
  3. runhaskell MoonEscape.hs
@carusocomputing

This comment has been minimized.

Copy link

carusocomputing commented Sep 4, 2018

Thanks! I'll check this out, I haven't used Haskell for anything other than code-golf! (Octopus Urn here)

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.