Last active
September 4, 2018 22:59
-
-
Save leftaroundabout/3955d27877e19be39d0f61fdafce069e to your computer and use it in GitHub Desktop.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
{-# 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 () |
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
To run the
.hs
file:cabal install dynamic-plot
runhaskell MoonEscape.hs