Help with RK5 program

Andy44 · Apr 26, 2008

I don't know if anyone can help with this, but I'm at wit's end. I've got Dan Boulet's book Methods of Orbit Determination for the Microcomputer (Willmann-Bell, 1991) which has BASIC programs for lots of cool stuff in the text. I'm trying to create a 2-body RK5 propagator in Excel from Boulet's program "RUNGE". Since Excel uses Visual BASIC for Applications, the code cannot be simply copied from the book, but it's very easy to translate. I get the inputs from cells and write the outputs to other cells, and the functions defined in the BASIC code I've written as seperate VBA functions. Also BASIC uses line numbers for GOTO statements, but you can set up a line label instead.

Anyway, I'm trying to propagate the orbit of Mars from JD 2446280.5 at 5-day intervals (using days, AUs, and solar masses for units), and it doesn't match the numerical example given in the book. I've gone over the code over and over again, and corrected all the typos. Same for the input data. I don't know if any of you has this book, let alone coded it up, but I'll post my code here anyway in case anyone can help. I tried calling the publisher to see if an errata sheet was available, but the line was disconnected and the website looks like it hasn't been updated in 2 years.

Code:

[FONT=Courier New]Public I As Integer, K As Integer, N As Integer
Public F, G, FP, GP, Rmag, Vmag, M, H, RKmag, VKmag
Public R(), V(), T(), FX(), F1(), F2(), F3(), F4(), F5(), F6()
Public RK(), VK(), GX(), G1(), G2(), G3(), G4(), G5(), G6()[/FONT]

[FONT=Courier New]Sub RUNGE()
'main routine[/FONT]
[FONT=Courier New]ReDim T(200), R(200, 3), V(200, 3)
ReDim FX(3), F1(3), F2(3), F3(3)
ReDim F4(3), F5(3), F6(3), RK(3)
ReDim GX(3), G1(3), G2(3), G3(3)
ReDim G4(3), G5(3), G6(3), VK(3)[/FONT]
[FONT=Courier New]'gather inputs
Name = Sheet1.Range("C8")
Equ = Sheet1.Range("C9")
Knum = Sheet1.Range("C10")
M = Sheet1.Range("C13")[/FONT]
[FONT=Courier New]T(0) = Sheet1.Range("C18")
R(0, 1) = Sheet1.Range("C19")
R(0, 2) = Sheet1.Range("C20")
R(0, 3) = Sheet1.Range("C21")
V(0, 1) = Sheet1.Range("C22")
V(0, 2) = Sheet1.Range("C23")
V(0, 3) = Sheet1.Range("C24")[/FONT]
[FONT=Courier New]DT = Sheet1.Range("C14")
NP = Sheet1.Range("C15")[/FONT]
[FONT=Courier New]Sheet2.Range("A2") = "Two-Body Motion by RK5 Integration"
Sheet2.Range("B3") = Name
Sheet2.Range("B4") = Equ[/FONT]
[FONT=Courier New]H = Knum * DT[/FONT]
[FONT=Courier New]For I = 0 To NP
    If I = 0 Then GoTo fourteenninety
    T(I) = T(I - 1) + DT
    Call RK5sub
    
    For K = 1 To 3
        R(I, K) = RK(K)
        V(I, K) = VK(K)
    Next K
    
fourteenninety:
    R(I, 0) = FNMG(R(I, 1), R(I, 2), R(I, 3))
    V(I, 0) = FNMG(V(I, 1), V(I, 2), V(I, 3))
    Sheet2.Cells(7 + I, 1) = T(I)
    Sheet2.Cells(7 + I, 2) = R(I, 1)
    Sheet2.Cells(7 + I, 3) = R(I, 2)
    Sheet2.Cells(7 + I, 4) = R(I, 3)
    Sheet2.Cells(7 + I, 8) = R(I, 0)
    Sheet2.Cells(7 + I, 5) = V(I, 1)
    Sheet2.Cells(7 + I, 6) = V(I, 2)
    Sheet2.Cells(7 + I, 7) = V(I, 3)
    Sheet2.Cells(7 + I, 9) = V(I, 0)[/FONT]
[FONT=Courier New]Next I[/FONT]
[FONT=Courier New]End Sub[/FONT]
[FONT=Courier New]Sub RK5sub()
'Runge-Kutta 5 Integration[/FONT]
[FONT=Courier New]'******STEP 1********[/FONT]
[FONT=Courier New]For K = 1 To 3
    RK(K) = R(I - 1, K)
    VK(K) = V(I - 1, K)
Next K[/FONT]
[FONT=Courier New]Call RK5subsub[/FONT]
[FONT=Courier New]For K = 1 To 3
    F1(K) = FX(K)
    G1(K) = GX(K)
Next K[/FONT]
[FONT=Courier New]'******STEP 2********[/FONT]
[FONT=Courier New]For K = 1 To 3
    RK(K) = R(I - 1, K) + (F1(K)) / 4
    VK(K) = V(I - 1, K) + (G1(K)) / 4
Next K[/FONT]
[FONT=Courier New]Call RK5subsub[/FONT]
[FONT=Courier New]For K = 1 To 3
    F2(K) = FX(K)
    G2(K) = GX(K)
Next K[/FONT]
[FONT=Courier New]'******STEP 3********[/FONT]
[FONT=Courier New]For K = 1 To 3
    RK(K) = R(I - 1, K) + (F1(K) + F2(K)) / 8
    VK(K) = V(I - 1, K) + (G1(K) + G2(K)) / 8
Next K[/FONT]
[FONT=Courier New]Call RK5subsub[/FONT]
[FONT=Courier New]For K = 1 To 3
    F3(K) = FX(K)
    G3(K) = GX(K)
Next K[/FONT]
[FONT=Courier New]'******STEP 4********[/FONT]
[FONT=Courier New]For K = 1 To 3
    RK(K) = R(I - 1, K) - (F2(K) - 2 * F3(K)) / 2
    VK(K) = V(I - 1, K) - (G2(K) - 2 * G3(K)) / 2
Next K[/FONT]
[FONT=Courier New]Call RK5subsub[/FONT]
[FONT=Courier New]For K = 1 To 3
    F4(K) = FX(K)
    G4(K) = GX(K)
Next K[/FONT]
[FONT=Courier New]'******STEP 5********[/FONT]
[FONT=Courier New]For K = 1 To 3
    RK(K) = R(I - 1, K) + (3 * F1(K) + 9 * F4(K)) / 16
    VK(K) = V(I - 1, K) + (3 * G1(K) + 9 * G4(K)) / 16
Next K[/FONT]
[FONT=Courier New]Call RK5subsub[/FONT]
[FONT=Courier New]For K = 1 To 3
    F5(K) = FX(K)
    G5(K) = GX(K)
Next K[/FONT]
[FONT=Courier New]'******STEP 6********[/FONT]
[FONT=Courier New]For K = 1 To 3
    RK(K) = R(I - 1, K) - (3 * F1(K) - 2 * F2(K) - 12 * F3(K) + 12 * F4(K) - 8 * F5(K)) / 7
    VK(K) = R(I - 1, K) - (3 * G1(K) - 2 * G2(K) - 12 * G3(K) + 12 * G4(K) - 8 * G5(K)) / 7
Next K[/FONT]
[FONT=Courier New]Call RK5subsub[/FONT]
[FONT=Courier New]For K = 1 To 3
    F6(K) = FX(K)
    G6(K) = GX(K)
Next K[/FONT]
[FONT=Courier New]'******RESULT********[/FONT]
[FONT=Courier New]For K = 1 To 3
    RK(K) = R(I - 1, K) + (7 * F1(K) + 32 * F3(K) + 12 * F4(K) + 32 * F5(K) _
    + 7 * F6(K)) / 90
    VK(K) = V(I - 1, K) + (7 * G1(K) + 32 * G3(K) + 12 * G4(K) + 32 * G5(K) _
    + 7 * G6(K)) / 90
Next K[/FONT]
[FONT=Courier New]End Sub[/FONT]
[FONT=Courier New]Sub RK5subsub()[/FONT]
[FONT=Courier New]RKmag = FNMG(RK(1), RK(2), RK(3))[/FONT]
[FONT=Courier New]For K = 1 To 3
    'FX(K) = H * FNF(VK(K))
    FX(K) = H * VK(K)
    GX(K) = H * FNG(M, RK(K), RKmag)
Next K[/FONT]
[FONT=Courier New]End Sub[/FONT]

And here are the function subroutines:

Code:

[FONT=Courier New]Function FNF(VK)
'f(v) = v[/FONT]
[FONT=Courier New]FNF = VK[/FONT]
[FONT=Courier New]End Function[/FONT]
[FONT=Courier New]Function FNMG(X, Y, Z)
'vector magnitude function (RSS)[/FONT]
[FONT=Courier New]FNMG = Sqr(X * X + Y * Y + Z * Z)[/FONT]
[FONT=Courier New]End Function[/FONT]
[FONT=Courier New]Function FNG(M, RK, R)[/FONT]
[FONT=Courier New]FNG = -M * RK / (R * R * R)[/FONT]
[FONT=Courier New]End Function[/FONT]
[FONT=Courier New]Function FNVS(X, Y, Z)
'vector squaring function[/FONT]
[FONT=Courier New]FNVS = (X ^ 2 + Y ^ 2 + Z ^ 2)[/FONT]
[FONT=Courier New]End Function[/FONT]

Thanks in advance for anyone who can help, and sorry if I'm wasting bandwidth.

Andy44 · Apr 26, 2008

Almost forgot: here is the input data:

Code:

[B][FONT=Arial][SIZE=2]Name:[/SIZE][/FONT][/B][FONT=Arial][SIZE=2]Mars[/SIZE][/FONT]
[B][FONT=Arial][SIZE=2]Equinox of Coordinates:[/SIZE][/FONT][/B][FONT=Arial][SIZE=2]J2000.0[/SIZE][/FONT]
[B][FONT=Arial][SIZE=2]Gravitational Constant k:[/SIZE][/FONT][/B][FONT=Arial][SIZE=2]0.01720209895[/SIZE][/FONT][SIZE=2][FONT=Arial]au3/2d-1m[/FONT][FONT=Wingdings 2]8[/FONT][FONT=Arial]-1/2[/FONT][/SIZE]
[B][SIZE=2][FONT=Arial]Central Body Mass m1:[/FONT][/SIZE][/B][FONT=Arial][SIZE=2]1.000000000[/SIZE][/FONT][SIZE=2][FONT=Arial]m[/FONT][FONT=Wingdings 2]8[/FONT][/SIZE]
[B][SIZE=2][FONT=Arial]Secondary Body Mass m2:[/FONT][/SIZE][/B][FONT=Arial][SIZE=2]0.000000323[/SIZE][/FONT][SIZE=2][FONT=Arial]m[/FONT][FONT=Wingdings 2]8[/FONT][/SIZE]
[B][FONT=Arial][SIZE=2]Combined Mass μ:[/SIZE][/FONT][/B][FONT=Arial][SIZE=2]1.000000323[/SIZE][/FONT][SIZE=2][FONT=Arial]m[/FONT][FONT=Wingdings 2]8[/FONT][/SIZE]
[B][FONT=Arial][SIZE=2]Step Interval dt (days):[/SIZE][/FONT][/B][FONT=Arial][SIZE=2]5.00[/SIZE][/FONT][FONT=Arial][SIZE=2]d[/SIZE][/FONT]
[B][FONT=Arial][SIZE=2]Number of steps[/SIZE][/FONT][/B][FONT=Arial][SIZE=2]144.00[/SIZE][/FONT]
[B][SIZE=2][FONT=Arial]Epoch t0:[/FONT][/SIZE][/B][FONT=Arial][SIZE=2]1985/08/03 00:00:00.000[/SIZE][/FONT][B][FONT=Arial][SIZE=2][COLOR=#808080]JD:[/COLOR][/SIZE][/FONT][/B][FONT=Arial][SIZE=2][COLOR=#808080]2446280.5[/COLOR][/SIZE][/FONT]
[B][SIZE=2][FONT=Arial]x0:[/FONT][/SIZE][/B][FONT=Arial][SIZE=2]-0.88884620[/SIZE][/FONT][FONT=Arial][SIZE=2]au[/SIZE][/FONT]
[B][SIZE=2][FONT=Arial]y0:[/FONT][/SIZE][/B][FONT=Arial][SIZE=2]1.24187820[/SIZE][/FONT][FONT=Arial][SIZE=2]au[/SIZE][/FONT]
[B][SIZE=2][FONT=Arial]z0:[/FONT][/SIZE][/B][FONT=Arial][SIZE=2]0.59365830[/SIZE][/FONT][FONT=Arial][SIZE=2]au[/SIZE][/FONT]
[B][SIZE=2][FONT=Arial]x-dot0:[/FONT][/SIZE][/B][FONT=Arial][SIZE=2]-0.65234240[/SIZE][/FONT][FONT=Arial][SIZE=2]au/d[/SIZE][/FONT]
[B][SIZE=2][FONT=Arial]y-dot0:[/FONT][/SIZE][/B][FONT=Arial][SIZE=2]-0.34498700[/SIZE][/FONT][FONT=Arial][SIZE=2]au/d[/SIZE][/FONT]
[B][SIZE=2][FONT=Arial]z-dot0:[/FONT][/SIZE][/B][FONT=Arial][SIZE=2]-0.14057410[/SIZE][/FONT][FONT=Arial][SIZE=2]au/d[/SIZE][/FONT]

Here is the expected first line (day 5) of output data:

Code:

[FONT=Courier New]JD            X          Y         Z         [/FONT]
[FONT=Courier New]2446285.50000 -0.9441931 1.2111721 0.5810728 [/FONT]

Here is the day 5 output I'm getting. Quite a difference:

Code:

[FONT=Courier New]2446285.50000 -0.94577527 1.2217878 0.5859845[/FONT]

Chode · Apr 26, 2008

First observations: your day 5 output is identical to your input. Also, your xyz-dots are on the order of 1 au/d, which doesn't seem right.

Regards

Andy44 · Apr 26, 2008

Chode said:
First observations: your day 5 output is identical to your input. Also, your xyz-dots are on the order of 1 au/d, which doesn't seem right.

Regards

Okay, I copied the wrong line from my spreadsheet, it's fixed now. I'll go look at the input again.

Andy44 · Apr 26, 2008

Okay, the velocity numbers match those given in the book exactly. I agree it seems pretty fast, but even if the state vector is wrong, it should produce the same output as the example in the book...

Help with RK5 program

Andy44

owner: Oil Creek Astronautix

Andy44

owner: Oil Creek Astronautix

Chode

Addon Developer

Andy44

owner: Oil Creek Astronautix

Andy44

owner: Oil Creek Astronautix

Similar threads