Front Office Football Central - Who's good with math and vba?

Apparently I'm not anymore.

I've got a series of equations that I'm trying to evaluate at work using vba and it's just not working out for me. Most of the evaluation I can get to work out, it's just the velocity term that is really getting away from me.

The velocity (dv/dt) term is shown below:

The biggest contributions to the velocity term is the Av term and the Can term.

K is a function of dv/dt where K = Ko*sin(theta) and is shown below.

and x can be found through dx/dt = v

I am evaluating the velocity term using a 4th order runge-kutta function with a step of 0.01 as v = f(t,v).

There is a term to check the growth of the error which I am evaluating using Simpson's 1/3 rule, which is shown below:

I have 1 series of data so I've been able to check my Can and Caer terms and they are giving the right numbers. The velocity I obtain from the dv/dt is much larger than it should be. The velocity I obtain from the check of the velocity is slightly higher than it should be. I've also included my code. I'm not really sure if I'm doing the RK method right as I've never evaluated an equation with the d/dt terms on the RHS before and am unsure how to approach them.

I'm sorry if this message is a big jumbled, let me know if there's anything you need to clarify it. I'm getting pretty frustrated as I been working on this off and on for a few days now and am really not sure what I'm doing, let alone what I am doing wrong.

Thanks for any help you can provide!

Code:

Sub VelocityEval()



Dim x As Variant

Dim Can As Variant

Dim Pmx As Variant

Dim Lam As Variant

Dim psi As Variant

Dim rn As Variant

Dim kn As Variant

With Sheets("Sheet1")

Can = Array(0, 0, 0) ' Anaerobic Energy

kn = Array(0, 0, 0, 0) ' k for R-K

Pmx = Array(.Cells(3, 18).Value, .Cells(3, 19).Value, .Cells(3, 20).Value) ' Maximum power constants

Lam = Array(.Cells(3, 21).Value, .Cells(3, 22).Value, .Cells(3, 23).Value) ' Lamda term

psi = Array(.Cells(3, 24).Value, .Cells(3, 25).Value, .Cells(3, 26).Value) ' Psi term

rn = Array(.Cells(3, 27).Value, .Cells(3, 28).Value, .Cells(3, 29).Value) ' rn = psi/lam

x = 0 ' horizontal distance

dT = 0.01 ' Time step

'Constants

R = 25

ho = 0.65

hcm = 1

grav = 9.81

Ko = 0.0029

A = 3.96

T = 0.01 ' Starting time

j = 0

out = 10.17 'First output at 10m + 0.17 starting gap

cout = 7

Vw = .Cells(3, 32).Value 'Wind term

intT = 0

hvo = 0

Cano = 0



Do While x < 100.17

' Calculate Aerobic Term

Caer = R * T - R * (1 - Exp(-1 * Lam(2) * T)) / Lam(2)

'Calculate Anerobic Term

For i = 0 To 2 ' Loop through 3 Pmx, psi, lam and rn's

    Temp1 = (Pmx(i) * rn(i) ^ (1 / rn(i))) * ((1 + rn(i)) / rn(i)) ^ ((1 + rn(i)) / rn(i))

    Temp2 = (rn(i) / (Lam(i) + psi(i)))

    Temp3 = Exp(-1 * Lam(i) * T) * (psi(i) + Lam(i) * (1 - Exp(-1 * psi(i) * T))) / (Lam(i) * (Lam(i) + psi(i)))

    Can(i) = Temp1 * (Temp2 - Temp3)

  

Next i



' Sum of the Anaerobic term

CanT = Can(0) + Can(1) + Can(2)



If CanT < 0 Then CanT = 0 'Anaerobic term must always be +ve



'For the first evaluation v = the following

If j = 0 Then

    vel = (2 * (CanT + Caer)) ^ 0.5

    dx = vel * dT ' Change is horizontal distance

    K = Ko

    

Else

    If x < 5 Then hv = hcm - (hcm - ho) * (5 - x) / 5 ' Change in hv

    delT = 0 ' change in T for RK method where v = fn(T,v)

    kmult = 0 ' term to be added to velocity term in RK method

    K = Ko * Sin(Atn(grav / (delV + K * (vi - Vw) ^ 2))) ' Calc K based on change in theta

    

    'RK evaluation for dv/dt

    For d = 0 To 3

        

        kn(d) = (-1 / (vi + kmult)) * (A * (vi + kmult) + K * (vi + kmult) * ((vi + kmult) - Vw) ^ 2 + grav * _

        (hv - ho) - CanT - R * (1 - Exp(-Lam(2) * (T + delT))))

        'Set the variation for the next RK loop for time

        If d = 2 Then delT = dT Else delT = dT / 2

        'and for velocity

        kmult = kn(d) * delT

    Next d

    'Calculate the new velocity from RK

    vel = vi + (1 / 6) * (kn(0) + 2 * (kn(1) + kn(2)) + kn(3)) * dT

    

    'Use simpson's rule to evaluate the error in the velocity and update the velocity term

    int1 = K * vel * ((vel - Vw) ^ 2) * T - grav * (hv - ho)

    int2 = K * vel * ((vel - Vw) ^ 2) * (T + dT / 2) - grav * (hv - ho)

    int3 = K * vel * ((vel - Vw) ^ 2) * (T + dT) - grav * (hv - ho)

    intT = dT * (int1 + 4 * int2 + int3) / 6

    vel = (2 * (CanT + Caer - A * x - intT)) ^ 0.5

    'Change in horizontal distance

    dx = vel * dT

End If



x = x + dx ' New horizontal distance

delV = vel - vi ' Change in velocity for evaluation of new K term

vi = vel ' Velocity for next iteration of RK

 





j = j + 1

'Output every 10m

If x > out Then

    .Cells(3, cout).Value = T

    .Cells(4, cout) = vel

    

    

    For m = 0 To 2

        .Cells(6 + m, cout).Value = Can(m)

    Next m

    .Cells(9, cout).Value = Caer

    cout = cout + 1

    out = out + 10

End If

T = T + dT





Loop



End With





End Sub