{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Answer to https://space.stackexchange.com/a/30442/8609"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"{-# LANGUAGE TypeFamilies, FlexibleContexts #-}\n",
"import Math.LinearMap.Category\n",
"import Data.VectorSpace\n",
"import Linear.V3\n",
"import Data.AffineSpace\n",
"import Control.Arrow\n",
"import Data.Semigroup\n",
"import qualified Diagrams.Prelude as Dia"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Types for the physical quantities"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"type ℝ = Double\n",
"type Distance = ℝ -- in m\n",
"type Pos = V3 Distance\n",
"type Speed = ℝ -- in m/s\n",
"type Velo = V3 Speed\n",
"type PhaseSpace = (Pos, Velo)\n",
"type Mass = ℝ -- in kg"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Standard, 4.-Ordnung Runge-Kutta method"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"rk4 :: (AffineSpace y, RealSpace (Diff y), t ~ ℝ)\n",
" => (y -> Diff y) -> Diff t -> (t,y) -> [(t,y)]\n",
"rk4 f h = go\n",
" where go (t,y) = (t,y) : go\n",
" (t+h, y .+^ h/6 · (k₁ ^+^ 2·k₂ ^+^ 2·k₃ ^+^ k₄))\n",
" where k₁ = f y\n",
" k₂ = f $ y .+^ h/2 · k₁\n",
" k₃ = f $ y .+^ h/2 · k₂\n",
" k₄ = f $ y .+^ h · k₃"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"import Graphics.Dynamic.Plot.R2\n",
"import qualified Diagrams.Prelude as Dia\n",
"import Data.List\n",
"\n",
"takeEach :: Int -> [a] -> [a]\n",
"takeEach n = go n\n",
" where go i (x:xs) | i>=n = x : go 1 xs\n",
" | otherwise = go (i+1) $ xs\n",
"\n",
"trajectoryPlot :: Int -> [(String, Distance)] -> [[(ℝ,ℝ)]] -> DynamicPlottable\n",
"trajectoryPlot speed meta = plotLatest\n",
" . map ( transpose . take 80 >>>\n",
" \\chunkCompos -> {-plotDelay (1/5) $ -} plotMultiple\n",
" [ (if name/=\"\" then legendName name else id)\n",
" $ plot [ lineSegPlot chunkCompo\n",
" , shapePlot . Dia.moveTo (Dia.p2 $ last chunkCompo)\n",
" . Dia.opacity 0.6\n",
" $ Dia.circle radius ]\n",
" | ((name,radius), chunkCompo) <- zip meta chunkCompos ]\n",
" )\n",
" . iterate (drop 4) . takeEach speed"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"type ThreeBody = ((PhaseSpace, PhaseSpace), PhaseSpace)\n",
"\n",
"me, mm, ms :: Mass\n",
"me = 5.972e24\n",
"mm = 7.346e22\n",
"ms = 2e3\n",
"\n",
"moonDist, moonRadius, earthRadius :: Distance\n",
"moonDist = 404e6 -- in m\n",
"moonRadius = 1.737e6 -- in m\n",
"earthRadius = 6.371e6 -- in m\n",
"\n",
"moonSpeed :: Speed\n",
"moonSpeed = 1020 -- in m/s\n",
"\n",
"gravConst :: ℝ\n",
"gravConst = 6.674e-11 -- in N·m²/kg²"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Trajectory of a 3-body system\n",
"Using the gravitational acceleration\n",
"$$\n",
"\\vec F = G\\cdot M\\cdot m\\cdot \\frac{\\vec{\\Delta x}}{\\|\\vec{\\Delta x}\\|^3}\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"gravAcc :: Mass -> Diff Pos -> Diff Velo\n",
"gravAcc mt δx = (gravConst * mt / magnitude δx^3) *^ δx\n",
"\n",
"traject3Body :: ThreeBody -> [ThreeBody]\n",
"traject3Body xv₀ = snd <$>\n",
" rk4 (\\(((xe,ve), (xm,vm)), (xs,vs))\n",
" -> (((ve, gravAcc mm $ xm.-.xe)\n",
" , (vm, gravAcc me $ xe.-.xm))\n",
" , (vs, gravAcc me (xe.-.xs) ^+^ gravAcc mm (xm.-.xs)) ) )\n",
" 90\n",
" (0, xv₀)"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"trebuchetOrbit :: Velo -> [ThreeBody]\n",
"trebuchetOrbit dv = replicate 1000{-0-} startState ++ traject3Body startState\n",
" where startState\n",
" = ( ((zeroV,zeroV)\n",
" ,(V3 moonDist 0 0, V3 0 moonSpeed 0))\n",
" , ( V3 (moonDist+ny*moonRadius) (-nx*moonRadius) 0, V3 0 moonSpeed 0 ^-^ dv ) )\n",
" V3 nx ny nz = dv ^/ magnitude dv"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Geocentric view"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"GraphWindowSpecR2{lBound=-6.597292373896027e8, rBound=6.597292373896027e8, bBound=-4.94796928042202e8, tBound=4.94796928042202e8, xResolution=640, yResolution=480}"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"dv = V3 (-2275) 1200 0 :: Velo\n",
"plotWindow [ trajectoryPlot 32 [(\"Earth\",earthRadius), (\"Moon\",moonRadius), (\"Spacecraft\",1)]\n",
" [ [(xe,ye),(xm,ym),(xs,ys)]\n",
" | (((V3 xe ye _,_),(V3 xm ym _,_)),(V3 xs ys _,_))\n",
" <- trebuchetOrbit dv ]\n",
" , dynamicAxes, unitAspect, forceXRange (-2*moonDist, 2*moonDist)\n",
" , forceYRange (-moonDist, moonDist) ]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"[![Animation of an Earth-influenced lunar orbit][1]][1]\n",
" [1]: https://i.stack.imgur.com/3rEKg.gif"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Lunacentric view"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"GraphWindowSpecR2{lBound=-1.2581314850581405e8, rBound=7.83201368639614e7, bBound=-4.8371011854243636e8, tBound=-3.3061026550688004e8, xResolution=640, yResolution=480}"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"plotWindow [ trajectoryPlot 4 [(\"Earth\",earthRadius), (\"Moon\",moonRadius), (\"Spacecraft\",1)]\n",
" [ [(xe-xm,ye-ym),(0,0),(xs-xm,ys-ym)]\n",
" | (((V3 xe ye _,_),(V3 xm ym _,_)),(V3 xs ys _,_))\n",
" <- trebuchetOrbit dv ]\n",
" , dynamicAxes, unitAspect, forceXRange (-20*moonRadius, 20*moonRadius)\n",
" , forceYRange (-15*moonRadius, 15*moonRadius) ]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"[![The unstable orbit, seem from the moon][2]][2] \n",
" [2]: https://i.stack.imgur.com/ZGofa.gif"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"GraphWindowSpecR2{lBound=-4635285.003931116, rBound=4995677.898723568, bBound=-3460111.0884955074, tBound=3763111.0884955074, xResolution=640, yResolution=480}"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"import Text.Printf\n",
"\n",
"veloVis = 2e3\n",
"\n",
"launch = [ Dia.p2 (x₀, y₀)\n",
" , Dia.p2 (x₀ - v₀x*veloVis, y₀ - v₀y*veloVis) ]\n",
" where x₀ = ny*moonRadius\n",
" y₀ = -nx*moonRadius\n",
" V3 v₀x v₀y _ = dv\n",
" V3 nx ny nz = dv ^/ magnitude dv\n",
"\n",
"plotWindow [ legendName \"To Earth\"\n",
" . shapePlot . Dia.dashingO [3,7] 0\n",
" $ Dia.arrowBetween (Dia.p2(0,0)) (Dia.p2(-5e6,0))\n",
" , legendName \"Moon\"\n",
" (shapePlot $ Dia.arrowBetween Dia.origin (Dia.p2 (0, moonSpeed*veloVis)))\n",
" <> shapePlot (Dia.opacity 0.6 $ Dia.circle moonRadius)\n",
" , legendName (printf \"v₀ = %.0g m/s\" $ magnitude dv)\n",
" . shapePlot $ Dia.arrowBetween (head launch) (last launch)\n",
" , unitAspect ]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"[![Start configuration needed for crash-orbit][2]][2] \n",
" [2]: https://i.stack.imgur.com/pHJ1o.png"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Haskell",
"language": "haskell",
"name": "haskell"
},
"language_info": {
"codemirror_mode": "ihaskell",
"file_extension": ".hs",
"name": "haskell",
"version": "8.2.1"
}
},
"nbformat": 4,
"nbformat_minor": 2
}