Lagrange-L4.ipynb

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