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

You know what? Screw it. Let's start again...

data State = State {velocity, position :: !Double}

data Derivative = Derivative {dvelocity, dposition :: !Double}

evaluate :: State -> Double -> Double -> Maybe Derivative -> Derivative
evaluate s t dt Nothing = Derivative {dvelocity = acceleration s t, 
dposition = velocity s}
evaluate s t dt (Just d) =
   let
     s' = State {velocity = velocity s + dvelocity d * dt, position = 
position s + dposition d * dt}
   in  Derivative {dvelocity = acceleration s' (t+dt), dposition = 
velocity s'}

integrate :: State -> Double -> Double -> State
integrate s t dt =
   let
     a = evaluate s t dt     Nothing
     b = evaluate s t (dt/2) (Just a)
     c = evaluate s t (dt/2) (Just b)
     d = evaluate s t dt     (Just c)
     dv = (dvelocity a + 2 * (dvelocity b + dvelocity c) + dvelocity d) / 6
     dp = (dposition a + 2 * (dposition b + dposition c) + dposition d) / 6
   in  State {velocity = velocity s + dv * dt, position = position s + 
dp * dt}

acceleration :: State -> Double -> Double
acceleration = error "not implemented yet"

Does that look correct?

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

Post a reply to this message

Invisible wrote:

> Does that look correct?

Well, I just ran it on a simple harmonic oscilator and got a sine wave, 
so that's always a good start...

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

Post a reply to this message

> I know I've asked this before, but what's the algorithm for this?

The wave equation.

http://en.wikipedia.org/wiki/Wave_equation

2nd time differential of "wave value" is proportional to the 2nd spatial 
differential of "wave value"

2nd time differential is just acceleration (use it in an integrator like RK4 
to get the new positions).

2nd spatial differential in 2D is just d2u/d2x + d2u/d2y

See attached POV code.

Post a reply to this message

Attachments:
Download 'us-ascii' (3 KB)

scott wrote:
>> I know I've asked this before, but what's the algorithm for this?
> 
> The wave equation.
> 
> http://en.wikipedia.org/wiki/Wave_equation
> 
> 2nd time differential of "wave value" is proportional to the 2nd spatial 
> differential of "wave value"
> 
> 2nd time differential is just acceleration (use it in an integrator like 
> RK4 to get the new positions).
> 
> 2nd spatial differential in 2D is just d2u/d2x + d2u/d2y
> 
> See attached POV code.

Right. So for each pixel, measure the difference against the four 
neighboring pixels, total that up, multiply by the wave propogation 
constant, and that's the total acceleration for this pixel. Seems simple 
enough.

What happens if the wave propogation medium propogates different 
frequencies at different velocities?

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

Post a reply to this message

> Right. So for each pixel, measure the difference against the four
> neighboring pixels, total that up, multiply by the wave propogation
> constant, and that's the total acceleration for this pixel. Seems simple
> enough.
> 
> What happens if the wave propogation medium propogates different
> frequencies at different velocities?

The wave propagation constant is not constant, but a function of 
frequency (dispersion relation).

Post a reply to this message

> What happens if the wave propogation medium propogates different 
> frequencies at different velocities?

Well then your c is a function of frequency, it's also probably a function 
of "wave value" and "wave steepness" too in real life. (higher and steeper 
waves might propogate faster than lower waves).  Good luck calculating the 
frequency at a point though :-)  I'm thinking you might be able to do 
something clever with fourier though...

Post a reply to this message

scott wrote:
>> What happens if the wave propogation medium propogates different 
>> frequencies at different velocities?
> 
> Well then your c is a function of frequency, it's also probably a 
> function of "wave value" and "wave steepness" too in real life. (higher 
> and steeper waves might propogate faster than lower waves).  Good luck 
> calculating the frequency at a point though :-)  I'm thinking you might 
> be able to do something clever with fourier though...

Well, Reaktor's "Steam Pipe" instrument simplates dispersion using a 
simple all-pass filter in the feedback loop. OTOH, it's only trying to 
simulate a 1D wave. At exactly 1 point in space. Hmm, maybe you could do 
some sort of interesting convulation operation?

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

Post a reply to this message

>> Well then your c is a function of frequency, it's also probably a 
>> function of "wave value" and "wave steepness" too in real life. (higher 
>> and steeper waves might propogate faster than lower waves).  Good luck 
>> calculating the frequency at a point though :-)  I'm thinking you might 
>> be able to do something clever with fourier though...
>
> Well, Reaktor's "Steam Pipe" instrument simplates dispersion using a 
> simple all-pass filter in the feedback loop. OTOH, it's only trying to 
> simulate a 1D wave. At exactly 1 point in space. Hmm, maybe you could do 
> some sort of interesting convulation operation?

I think if you are simulating a musical instrument, it would be better to 
get the resonant parts working first before you start trying to simulate 
different propagation speeds at different frequencies.  I would imagine that 
is a much smaller secondary effect.

Post a reply to this message

scott wrote:
>>> Well then your c is a function of frequency, it's also probably a 
>>> function of "wave value" and "wave steepness" too in real life. 
>>> (higher and steeper waves might propogate faster than lower waves).  
>>> Good luck calculating the frequency at a point though :-)  I'm 
>>> thinking you might be able to do something clever with fourier though...
>>
>> Well, Reaktor's "Steam Pipe" instrument simplates dispersion using a 
>> simple all-pass filter in the feedback loop. OTOH, it's only trying to 
>> simulate a 1D wave. At exactly 1 point in space. Hmm, maybe you could 
>> do some sort of interesting convulation operation?
> 
> I think if you are simulating a musical instrument, it would be better 
> to get the resonant parts working first before you start trying to 
> simulate different propagation speeds at different frequencies.  I would 
> imagine that is a much smaller secondary effect.

Oh yeah, sure. I'm just going off at a tangent now. ;-)

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

Post a reply to this message