A universe in 125 lines of haskell

Newtonian.hs

--A universe in 125 lines of haskell: http://i.imgur.com/9dJdaV1.png
--Note: I had made a better version but that computer died :/
--to build: cabal install random gloss && ghc -O3 -threaded Newtonian.hs && ./Newtonian
import Graphics.Gloss
import Graphics.Gloss.Data.Picture
import Graphics.Gloss.Data.Vector
import Graphics.Gloss.Interface.Pure.Simulate
import Data.List
import System.Random

--VECTOR OPERATIONS--
addV :: Vector -> Vector -> Vector
addV (x1, y1) (x2, y2) = (x1+x2, y1+y2)

subV :: Vector -> Vector -> Vector
subV (x1, y1) (x2, y2) = (x1-x2, y1-y2)

mulV :: Float -> Vector -> Vector
mulV s (x, y) = (x*s, y*s)

divV :: Vector -> Float -> Vector
divV (x, y) s = (x/s, y/s)


--DATA TYPES--
--a simple representation of a point mass, with a mass, location, and velocity
data Body = Body { mass :: !Float, loc :: !Point, vel :: !Vector } deriving (Eq, Show)

--a force (mass * acceleration)
type Force = Vector

--a Universe state (a list of bodies)
type Universe = [Body]

--a law of physics (given a body and a universe, it calculates the force acting on the body from the universe)
type Law = Body -> Universe -> Force


--SIMULATION FUNCTIONS--
--calculate the net force acting on a body given a list of forces
net :: [Force] -> Force
net = foldr addV (0,0)

--calculate the next state of a body, given the list of forces acting on it
stepBody :: Float -> [Force] -> Body -> Body
stepBody t vs m = Body (mass m)                                        --mass
                ((loc m) `addV` (t `mulV` (vel m)))                    --location + velocity*dt
                ((vel m) `addV` (t `mulV` ((net vs) `divV` (mass m)))) --velocity + dt*net_force/m

--calculate all of the forces acting on all of the bodies within a universe, given a set of laws
forces :: [Law] -> Universe -> [[Force]]
forces laws u = map forcesIndiv u
    where forcesIndiv m = [l m u | l <- laws]


--EXAMPLE FORCES
--simple downward acceleration (at 9.81 pixels / second^2)
fall :: Law
fall m u = (mass m) `mulV` (0, -9.81)

--gravity a la Newton
gravity :: Law
gravity m u = foldl addV (0, 0) $ map (attract m) u

attract :: Body -> Body -> Force
attract m1 m2
    | (loc m1) == (loc m2) = (0,0)
    | otherwise = g * ((mass m1) * (mass m2) / (magV diff)**2) `mulV` (normaliseV diff)
    where diff = subV (loc m2) (loc m1)
          g = 50 --gravitational constant

--keeping the bodies within view
cohese :: Law
cohese m u 
    | (loc m) == (0, 0) = (0,0)
    | otherwise         = (-1)*(mass m) `mulV` (loc m)

--physics is just a list of laws
physics :: [Law]
physics = [gravity] --only gravity, to start

--INITIAL UNIVERSE STATE--
--adjust parameters at will
count :: Int
count = 50 --number of masses

initial :: Universe
initial = zipWith5 (\a b c d e -> Body a (b, c) (d,e)) masses xls yls xvs yvs
    where   masses = randomListBounded (10, 20)      count $ mkStdGen 0        --upper and lower bound for mass
            xls    = randomListBounded ((-100), 100)     count $ mkStdGen 1    --upper and lower bound for initial x location
            yls    = randomListBounded ((-100), 100)     count $ mkStdGen 2    --y location
            xvs    = randomListBounded ((-10), 10)   count $ mkStdGen 3        --x velocity
            yvs    = randomListBounded ((-10), 10)   count $ mkStdGen 4        --y velocity
            randomListBounded :: (Float, Float) -> Int -> StdGen -> [Float]    --convenience
            randomListBounded (a,b) n = take n . unfoldr (Just . randomR (a, b)) 

--DRAWING FUNCTIONS--
renderForce :: Float -> Force -> Picture
renderForce r v = Translate rx ry $ Pictures [long, chev1, chev2]
    where   (rx, ry) = r `mulV` n
            long     = Line [(0,0), v]
            chev1    = Rotate 15 $ Line [(0,0), n]
            chev2    = Rotate (-15) $ Line [(0,0), n]
            n        = normaliseV v

renderBody :: Body -> [Force] -> Picture
renderBody m fs = Translate x y $ Pictures [op, velp, fsp]
    where   (x, y)	= loc m
            r		= sqrt $ mass m
            op		= Color red $ Circle r                                                --circle for body
            velp	= Color green $ Line [(0,0), r `mulV` (normaliseV $ vel m)]           --green line for velocity
            fsp		= Color white $ Pictures $ map (renderForce r) $ map (mulV (-1)) fs   --white arrows for forces

--FINAL FUNCTIONS--
--calculate the next state of the universe
step :: ViewPort -> Float -> Universe -> Universe
step v t u = zipWith (stepBody (t*1)) (forces physics u) u

--draw the universe with Gloss
render :: Universe -> Picture
render u = Pictures $ zipWith renderBody u $ forces physics u

--MAIN--
main :: IO ()
main = simulate (InWindow "Newtonian" (500,500) (200,200)) black 60 initial render step