sincos.hs

import Data.List
import qualified Data.Vec as Vec

dt = pi/100
t0 = 0
tmax = 50

dydt :: (Double -> Vec.Vec2D) -> Double -> Vec.Vec2D
dydt y t =
  Vec.Vec2D
    (Vec.get Vec.n0 (y (t - 3*pi/2)))
    (Vec.get Vec.n0 (-(y (t))))

initialY :: [Vec.Vec2D]
initialY = map (\x -> Vec.Vec2D (sin x) (cos x)) [0, (-dt) .. (-2*pi) ]

--- Integrator
newval (oldYs, oldT) =
  let
    y t = oldYs !! round ((oldT - t)/dt)
    newY = (head oldYs) + (Vec.fromList $ repeat dt) * (dydt y oldT)
  in
    Just (newY, (newY:oldYs, oldT + dt))

integrated = unfoldr newval (initialY, t0)

--- Display results
main = sequence $ take (round ((tmax-t0)/dt)) $ map print integrated