Use Newtonian gravity in model

Main.hs

@@ -8,32 +8,49 @@ import Graphics.Gloss.Geometry.Angle
 
 data PolarCoord = PolarCoord { radius :: Float, angle :: Float }
 
+data State = State { position :: PolarCoord
+                   , radialMomentum :: Float
+                   , angularMomentum :: Float }
+
 majorMass :: Float
-majorMass = 10
+majorMass = 200
 
 minorMass :: Float
-minorMass = 4
+minorMass = 20
+
+initialState :: State
+initialState = State (PolarCoord 300 0) 0 0.001
 
 framesPerSecond :: Int
 framesPerSecond = 120
 
-angularVelocity :: Float
-angularVelocity = pi
+render :: State -> Picture
+render state = pictures [ circleSolid (sqrt majorMass)
+                        , translatePolar (position state) $ circleSolid (sqrt minorMass) ]
 
-render :: PolarCoord -> Picture
-render state = pictures [ circleSolid majorMass
-                        , translatePolar state $ circleSolid minorMass ]
+step :: a -> b -> State -> State
+step _ _ state = state { position = nextPosition state
+                       , radialMomentum = nextRadialMomentum state }
 
-step :: a -> b -> PolarCoord -> PolarCoord
-step _ _ state = let radiansPerFrame = angularVelocity / fromIntegral framesPerSecond
-                     nextAngle = normalizeAngle $ angle state + radiansPerFrame
-                 in state { angle = nextAngle }
+nextPosition :: State -> PolarCoord
+nextPosition state = let r = radius . position $ state
+                         a = angle . position $ state
+                         rDot = radialMomentum state / minorMass
+                         aDot = angularMomentum state / (minorMass * r ^ 2)
+                     in PolarCoord (r + rDot) (a + aDot)
+
+nextRadialMomentum :: State -> Float
+nextRadialMomentum state = let r = radius . position $ state
+                               p = radialMomentum state
+                               l = angularMomentum state
+                               pDot = (l ^ 2 / (2 * minorMass * r) - minorMass * majorMass) / r ^ 2
+                           in p + pDot
 
 translatePolar :: PolarCoord -> Picture -> Picture
 translatePolar q = rotate (radToDeg $ angle q) . translate (radius q) 0
 
 window :: Display
-window = InWindow "Foo" (200, 200) (100, 100)
+window = InWindow "Foo" (800, 800) (100, 100)
 
 main :: IO ()
-main = simulate window white framesPerSecond (PolarCoord 80 0) render step
+main = simulate window white framesPerSecond initialState render step