Right. Well maybe it would be helpful to fix a few typos...

>   evaluate (v,p) t dt (dv,dp) =
>     let
>       v' = v + dt*dv
>       p' = p + dt*dp
>     in  (v', a (t+dt) v' p')

Err... that's incorrect. The acceleration function ("a") needs to be a 
parameter somewhere. Or it needs to be a top-level function. Let's make 
it top-level, like in the C++ version:

   acceleration t v p = ???

   evaluate (v,p) t dt Nothing = (v, acceleration t)
   evaluate (v,p) t dt (Just (dv,dp)) =
     let
       v' = v + dt*dv
       p' = p + dt*dp
     in  (v', acceleration (t+dt) v' p')

[I made the derivatives optional too. Did I get it right?]

>   integrate t dt a (v,p) =
>     let
>       a = evaluate (v,p) t dt -- er... where's the rest?
>       b = evaluate (v,p) t (dt/2) a
>       c = evaluate (v,p) t (dt/2) b
>       d = evaluate (v,p) t dt     c
>       dv = (fst a + 2 * (fst b + fst c) + fst d) / 6
>       dp = (snd a + 2 * (snd b + snd c) + snd d) / 6
>     in (v + dt*dv, p + dt*dp)

...and if we're using a top-level acceleration function, you don't need 
that "a" parameter there (which incidentally causes a name clash).

   integrate t dt (v,p) =
     let
       a = evaluate (v,p) t dt     Nothing
       b = evaluate (v,p) t (dt/2) (Just a)
       c = evaluate (v,p) t (dt/2) (Just b)
       d = evaluate (v,p) t dt     (Just c)
       dv = (fst a + 2 * (fst b + fst c) + fst d) / 6
       dp = (snd a + 2 * (snd b + snd c) + snd d) / 6
     in  (v + dt*dv, p + dt*dp)

I'll have to try this out later and see if it actually works...

-- 
http://blog.orphi.me.uk/
http://www.zazzle.com/MathematicalOrchid*

Post a reply to this message