Satellite position computation using Runge-Kutta 4

Question

my issue is related to Runge-Kutta 4 (RK4) method and the correct iteration steps required for the state vector of an orbiting satellite. The below code (in Python) describes the motion based on the description as per this link (http://www.navipedia.net/index.php/GLONASS_Satellite_Coordinates_Computation):

    if total_step_number != 0:   
        for i in range(1, total_step_number+1):                             
            #Calculate k1                
            k1[0] = (-cs.GM_GLONASS * XYZ[0] / radius**3) \
             + ((3/2) * cs.C_20 * cs.GM_GLONASS * cs.SEMI_MAJOR_AXIS_GLONASS**2 * XYZ[0] * (1 - (5*(XYZ[2]**2) / (radius**2))) / radius**5) \
             + XYZDDot[0] + (cs.OMEGAE_DOT**2 * XYZ[0]) + (2 * cs.OMEGAE_DOT * XYZDot[1])
            k1[1] = (-cs.GM_GLONASS * XYZ[1] / radius**3) \
             + ((3/2) * cs.C_20 * cs.GM_GLONASS * cs.SEMI_MAJOR_AXIS_GLONASS**2 * XYZ[1] * (1 - (5*(XYZ[2]**2) / (radius**2))) / radius**5) \
             + XYZDDot[1] + (cs.OMEGAE_DOT**2 * XYZ[1]) - (2 * cs.OMEGAE_DOT * XYZDot[0])
            k1[2] = (-cs.GM_GLONASS * XYZ[2] / radius**3) \
             + ((3/2) * cs.C_20 * cs.GM_GLONASS * cs.SEMI_MAJOR_AXIS_GLONASS**2 * XYZ[2] * (3 - (5*(XYZ[2]**2) / (radius**2))) / radius**5) \
             + XYZDDot[2]

            #Intermediate step to bridge k1 to k2
            XYZ2[0] = XYZ[0] + (XYZDot[0] * h / 2) + (k1[0] * h**2 / 8)
            XYZDot2[0] = XYZDot[0] + (k1[0] * h / 2)
            XYZ2[1] = XYZ[1] + (XYZDot[1] * h / 2) + (k1[1] * h**2 / 8)
            XYZDot2[1] = XYZDot[1] + (k1[1] * h / 2)
            XYZ2[2] = XYZ[2] + (XYZDot[2] * h / 2) + (k1[2] * h**2 / 8)
            XYZDot2[2] = XYZDot[2] + (k1[2] * h / 2)
            radius = np.sqrt((XYZ2[0]**2)+(XYZ2[1]**2)+(XYZ2[2]**2))

             ....

There is more code however I want to limit what I show for now since it's the intermediate steps I'm most interested in resolving. Basically, for those familiar with state vectors and using RK4, you can see that the position and velocity is updated in the intermediate step, but not the acceleration. My question is related to the calculation required in order to update too the acceleration. It would begin:

XYZDDot[0] = ...
XYZDDot[1] = ...
XYZDDot[2] = ...

...but what exactly comes after is not very clear. Any advice welcome.

Below is the full code:

        for j in h_step_values:
            h = j    
            if h > 0:
                one_way_iteration_steps = one_way_iteration_steps -1
            elif h < 0:
                one_way_iteration_steps = one_way_iteration_steps +1
                XYZ = initial_XYZ
            #if total_step_number != 0:   
            for i in range(0, one_way_iteration_steps):                             
                #Calculate k1                
                k1[0] = (-cs.GM_GLONASS * XYZ[0] / radius**3) \
                 + ((3/2) * cs.C_20 * cs.GM_GLONASS * cs.SEMI_MAJOR_AXIS_GLONASS**2 * XYZ[0] * (1 - (5*(XYZ[2]**2) / (radius**2))) / radius**5) \
                 + XYZDDot[0] + (cs.OMEGAE_DOT**2 * XYZ[0]) + (2 * cs.OMEGAE_DOT * XYZDot[1])
                k1[1] = (-cs.GM_GLONASS * XYZ[1] / radius**3) \
                 + ((3/2) * cs.C_20 * cs.GM_GLONASS * cs.SEMI_MAJOR_AXIS_GLONASS**2 * XYZ[1] * (1 - (5*(XYZ[2]**2) / (radius**2))) / radius**5) \
                 + XYZDDot[1] + (cs.OMEGAE_DOT**2 * XYZ[1]) - (2 * cs.OMEGAE_DOT * XYZDot[0])
                k1[2] = (-cs.GM_GLONASS * XYZ[2] / radius**3) \
                 + ((3/2) * cs.C_20 * cs.GM_GLONASS * cs.SEMI_MAJOR_AXIS_GLONASS**2 * XYZ[2] * (3 - (5*(XYZ[2]**2) / (radius**2))) / radius**5) \
                 + XYZDDot[2]

                #Intermediate step to bridge k1 to k2
                XYZ2[0] = XYZ[0] + (XYZDot[0] * h / 2) + (k1[0] * h**2 / 8)
                XYZDot2[0] = XYZDot[0] + (k1[0] * h / 2)
                XYZDDot2[0] = XYZDDot[0] + (k1[0] * h / 2)
                XYZ2[1] = XYZ[1] + (XYZDot[1] * h / 2) + (k1[1] * h**2 / 8)
                XYZDot2[1] = XYZDot[1] + (k1[1] * h / 2)
                XYZ2[2] = XYZ[2] + (XYZDot[2] * h / 2) + (k1[2] * h**2 / 8)
                XYZDot2[2] = XYZDot[2] + (k1[2] * h / 2)
                radius = np.sqrt((XYZ2[0]**2)+(XYZ2[1]**2)+(XYZ2[2]**2))

                #Calculate k2  
                k2[0] = (-cs.GM_GLONASS * XYZ2[0] / radius**3) \
                 + ((3/2) * cs.C_20 * cs.GM_GLONASS * cs.SEMI_MAJOR_AXIS_GLONASS**2 * XYZ2[0] * (1 - (5*(XYZ2[2]**2) / (radius**2))) / radius**5) \
                 + XYZDDot[0] + (cs.OMEGAE_DOT**2 * XYZ2[0]) + (2 * cs.OMEGAE_DOT * XYZDot2[1])
                k2[1] = (-cs.GM_GLONASS * XYZ2[1] / radius**3) \
                 + ((3/2) * cs.C_20 * cs.GM_GLONASS * cs.SEMI_MAJOR_AXIS_GLONASS**2 * XYZ2[1] * (1 - (5*(XYZ2[2]**2) / (radius**2))) / radius**5) \
                 + XYZDDot[1] + (cs.OMEGAE_DOT**2 * XYZ2[1]) - (2 * cs.OMEGAE_DOT * XYZDot2[0])
                k2[2] = (-cs.GM_GLONASS * XYZ2[2] / radius**3) \
                 + ((3/2) * cs.C_20 * cs.GM_GLONASS * cs.SEMI_MAJOR_AXIS_GLONASS**2 * XYZ2[2] * (3 - (5*(XYZ2[2]**2) / (radius**2))) / radius**5) \
                 + XYZDDot[2]

                #Intermediate step to bridge k2 to k3
                XYZ2[0] = XYZ[0] + (XYZDot[0] * h / 2) + (k2[0] * h**2 / 8)
                XYZDot2[0] = XYZDot[0] + (k2[0] * h / 2)
                XYZ2[1] = XYZ[1] + (XYZDot[1] * h / 2) + (k2[1] * h**2 / 8)
                XYZDot2[1] = XYZDot[1] + (k2[1] * h / 2)
                XYZ2[2] = XYZ[2] + (XYZDot[2] * h / 2) + (k2[2] * h**2 / 8)
                XYZDot2[2] = XYZDot[2] + (k2[2] * h / 2)
                radius = np.sqrt((XYZ2[0]**2)+(XYZ2[1]**2)+(XYZ2[2]**2))

                #Calculate k3  
                k3[0] = (-cs.GM_GLONASS * XYZ2[0] / radius**3) \
                 + ((3/2) * cs.C_20 * cs.GM_GLONASS * cs.SEMI_MAJOR_AXIS_GLONASS**2 * XYZ2[0] * (1 - (5*(XYZ2[2]**2) / (radius**2))) / radius**5) \
                 + XYZDDot[0] + (cs.OMEGAE_DOT**2 * XYZ2[0]) + (2 * cs.OMEGAE_DOT * XYZDot2[1]) 
                k3[1] = (-cs.GM_GLONASS * XYZ2[1] / radius**3) \
                 + ((3/2) * cs.C_20 * cs.GM_GLONASS * cs.SEMI_MAJOR_AXIS_GLONASS**2 * XYZ2[1] * (1 - (5*(XYZ2[2]**2) / (radius**2))) / radius**5) \
                 + XYZDDot[1] + (cs.OMEGAE_DOT**2 * XYZ2[1]) - (2 * cs.OMEGAE_DOT * XYZDot2[0])
                k3[2] = (-cs.GM_GLONASS * XYZ2[2] / radius**3) \
                 + ((3/2) * cs.C_20 * cs.GM_GLONASS * cs.SEMI_MAJOR_AXIS_GLONASS**2 * XYZ2[2] * (3 - (5*(XYZ2[2]**2) / (radius**2))) / radius**5) \
                 + XYZDDot[2]

                #Intermediate step to bridge k3 to k4
                XYZ2[0] = XYZ[0] + (XYZDot[0] * h) + (k3[0] * h**2 / 2)
                XYZDot2[0] = XYZDot[0] + (k3[0] * h)
                XYZ2[1] = XYZ[1] + (XYZDot[1] * h) + (k3[1] * h**2 / 2)
                XYZDot2[1] = XYZDot[1] + (k3[1] * h)
                XYZ2[2] = XYZ[2] + (XYZDot[2] * h) + (k3[2] * h**2 / 2)
                XYZDot2[2] = XYZDot[2] + (k3[2] * h)
                radius = np.sqrt((XYZ2[0]**2)+(XYZ2[1]**2)+(XYZ2[2]**2))

                #Calculate k4 
                k4[0] = (-cs.GM_GLONASS * XYZ2[0] / radius**3) \
                 + ((3/2) * cs.C_20 * cs.GM_GLONASS * cs.SEMI_MAJOR_AXIS_GLONASS**2 * XYZ2[0] * (1 - (5*(XYZ2[2]**2) / (radius**2))) / radius**5) \
                 + XYZDDot[0] + (cs.OMEGAE_DOT**2 * XYZ2[0]) + (2 * cs.OMEGAE_DOT * XYZDot2[1])
                k4[1] = (-cs.GM_GLONASS * XYZ2[1] / radius**3) \
                 + ((3/2) * cs.C_20 * cs.GM_GLONASS * cs.SEMI_MAJOR_AXIS_GLONASS**2 * XYZ2[1] * (1 - (5*(XYZ2[2]**2) / (radius**2))) / radius**5) \
                 + XYZDDot[1] + (cs.OMEGAE_DOT**2 * XYZ2[1]) - (2 * cs.OMEGAE_DOT * XYZDot2[0]) 
                k4[2] = (-cs.GM_GLONASS * XYZ2[2] / radius**3) \
                 + ((3/2) * cs.C_20 * cs.GM_GLONASS * cs.SEMI_MAJOR_AXIS_GLONASS**2 * XYZ2[2] * (3 - (5*(XYZ2[2]**2) / (radius**2))) / radius**5) \
                 + XYZDDot[2]


                for p in range(3):
                    XYZ[p] = XYZ[p] + XYZDot[p] * h + h**2 * ((k1[p] + 2*k2[p] + 2*k3[p] + k4[p]) / 12)
                    XYZDot[p] = XYZDot[p] + (h * (k1[p] + 2*k2[p] + 2*k3[p] + k4[p]) / 6)

                radius = np.sqrt((XYZ[0])**2 + (XYZ[0])**2 + (XYZ[0])**2)

Answer 1

The equation you are solving is of the type

ddot x = a(x)

where a(x) is the acceleration which is computed in your k1 computation. Indeed, the first order system would be

dot v = a(x)
dot x = v

The RK4 implementation thus starts with

k1 = a(x)
l1 = v

k2 = a(x+l1*h/2) = a(x+v*h/2)
l2 = v+k1*h/2

etc. The use of the l1,l2,... seems implicit in the code, inserting these linear combinations directly where they occur.

---

In short, you are not missing the acceleration computation, it is the main part of the code fragment.

---

Update: (8/22) To come closer to the intention of the intermediate bridge steps, the abstract code should read ( with (* .. *) denoting comments or unnecessary computations)

k1 = a(x)                    (* l1 = v *)

x2 = x + v*h/2               (* v2 = v + k1*h/2 *)

k2 = a(x2)                   (* l2 = v2 *)

x3 (* = x + l2*h/2 *) 
   = x + v*h/2 + k1*h^2/4    (* v3 = v + k2*h/2 *)

k3 = a(x3)                   (* l3 = v3 *)

x4 (* = x + l3*h *)
   = x + v*h + k2*h^2/2      (* v4 = v + k3*h *)

k4 = a(x4)                   (* l4 = v4 *)


delta_v = ( k1+2*(k2+k3)+k4 ) * h/6

delta_x (* = ( l1+2*(l2+l3)+l4 ) * h/6 *)
        = v*h + (k1+k2+k3) * h^2/6