-- NBODY SIMULATOR IN HASKELL, mid MMXVI
-- written by Philip K., [http://pk.3n.cc]
-- for licensing, refer to the README

module Main where

import System.Environment
import Text.Printf
import Linear.V3
import Linear.Vector
import Linear.Metric

-- CONSTANTS & UTILITIES
gc :: Double
gc = 6.67408e-11

-- DATA TYPES & VECTOR OPERATIONS
type Body = (Double, Double,  V3 Double, V3 Double)
type State = [Body]

-- LINEAR CHANGE CALCULATOR
step :: Double -> State -> State
step dt prev = map (\b@(m,r,v,d) -> let
                       v' = v ^+^ dt *^ (sumV $ map (acc b) prev)
                       d' = d ^+^ dt *^ v ^+^ (dt ** 2 / 2) *^ (sumV $ map (acc b) prev)
                     in (m,r,v',d')) $ collide prev
  where
    acc (_,_,_,d1) (m2,_,_,d2) = if d1 == d2 then zero 
      else ((gc * m2) / ((sqrt . foldr1 (+) . fmap (** 2)) $ d1 ^-^ d2) ** 3) *^ (d2 ^-^ d1)
    collision b1@(m1,r1,v1,d1) (m2,r2,v2,d2) =
      if d1 `distance` d2 > r1 + r2 || d1 == d2 then b1
      else (let nr = d2 ^-^ d1
                vn v = nr ^* (nr `dot` v) ^/ (quadrance nr)
                v1' = v1 ^-^ vn v1 ^+^ (m1 *^ vn v1 ^+^ (m2 *^ (2 *^ vn v2 ^-^ vn v1))) ^/ (m1 + m2)
             in (m1, r1, v1', d1))
    collide = map (\b -> foldl collision b prev)

-- SIMULATION CONTROLLER
ctrl :: Double -> Int -> State -> IO ()
ctrl dt sk = mapM_ output . skip . iterate (step dt)
  where output = putStrLn . unwords . map (\(_,_,_,(V3 d1 d2 d3)) -> printf "%.8f %.8f %.8f" d1 d2 d3)
        skip [] = []
        skip (x:xs) = x : skip (drop sk xs)

-- SIMULATION CONFIGURATION
setup :: [String] -> State
setup = map (\s -> (read $ procArg s !! 0,
                    read $ procArg s !! 1,
                    procVec $ procArg s !! 2,
                    procVec $ procArg s !! 3))
  where procVec vec = let p = map read $ splitOn (',' ==) vec
                      in V3 (p !! 0) (p !! 1) (if length p >= 3 then p !! 2 else 0)
        procArg = splitOn (';' ==)
        splitOn _ [] = []
        splitOn f s = w : splitOn f (if length r == 0 then [] else tail r)
          where (w,r) = break f s

main :: IO ()
main = getArgs >>= \(dt:skip:bodies) -> ctrl (read dt) (read skip) (setup bodies)