```\section{Runge-Kutta}

Haskell implementation of RK4.

\begin{code}
module RSAGL.RK4
(rk4,
integrateRK4,
rk4',
integrateRK4')
where

import RSAGL.AbstractVector
import RSAGL.Time
\end{code}

\begin{code}
genericRK4 :: (AbstractVector v) => (Time -> p -> v -> p) -> (Time -> p -> Rate v) -> p -> Time -> Time -> p
genericRK4 addPV diffF p0 t0 t1 = addPV h p0 \$ (scalarMultiply (recip 6) (abstractSum [k1,k2,k2,k3,k3,k4]) `over` h)
where k1 = diffF t0 p0
k2 = diffF tmid (addPV h2 p0 \$ k1 `over` h2)
k3 = diffF tmid (addPV h2 p0 \$ k2 `over` h2)
k4 = diffF t1 (addPV h p0 \$ k3 `over` h)
h = t1 `sub` t0
h2 = scalarMultiply (recip 2) h
tmid = t0 `add` h2

rk4 :: (AbstractVector v) => (p -> v -> p) -> (Time -> p -> Rate v) -> p -> Time -> Time -> p

rk4' :: (AbstractVector v) => (p -> v -> p) -> (Time -> p -> Rate v -> Acceleration v) -> (p,Rate v) -> Time -> Time -> (p,Rate v)
rk4' addPV diffF = genericRK4
(\t (p,old_v) delta_v -> let new_v = old_v `add` delta_v in (addPV p \$ (scalarMultiply (recip 2) \$ old_v `add` new_v) `over` t,new_v))
(\t (p,v) -> diffF t p v)

genericIntegrate :: (p -> Time -> Time -> p) -> p -> Time -> Time -> Integer -> p
genericIntegrate _ pn _ _ 0 = pn
genericIntegrate f p0 t0 tn n = genericIntegrate f p1 t1 tn (n-1)
where t1 = t0 `add` (scalarMultiply (recip \$ realToFrac n) \$ tn `sub` t0)
p1 = f p0 t0 t1

integrateRK4 :: (AbstractVector v) => (p -> v -> p) -> (Time -> p -> Rate v) -> p -> Time -> Time -> Integer -> p
integrateRK4 addPV diffF = genericIntegrate \$ rk4 addPV diffF

integrateRK4' :: (AbstractVector v) => (p -> v -> p) -> (Time -> p -> Rate v -> Acceleration v) -> (p,Rate v) -> Time -> Time -> Integer -> (p,Rate v)
integrateRK4' addPV diffF = genericIntegrate \$ rk4' addPV diffF
\end{code}
```