MoonEscape.hs

{-# 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 ()

MoonTrebuchet.ipynb

{
 "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
}

To run the .hs file:

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

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