04-05-2016, 11:08 PM

This HP Prime program demonstrates how to interact with a routine called rkf78 which solves a system of ordinary differential equations. The computer program is based on the Runge-Kutta-Felhberg 7(8) numerical method. This algorithm is a variable step size numerical method that allows the user to specify a truncation error tolerance. This tolerance can be used to "calibrate" the algorithm.

This demo program solves the system of Keplerian, unperturbed orbital motion. For this example, a measure of the effectiveness of this algorithm can be determined by how well the orbital closes on itself after one complete orbital period.

Here's the source code for the equations of motion implemented in this example.

EXPORT keqm (t, y)

// first order equations of orbital motion

// Keplerian motion - no perturbations

// input

// t = simulation time (seconds)

// y = state vector

// output

// ydot = integration vector

/////////////////////////////

BEGIN

LOCAL rmag, rcubed;

LOCAL ydot := [0, 0, 0, 0, 0, 0];

rmag := sqrt(y(1) * y(1) + y(2) * y(2) + y(3) * y(3));

rcubed := rmag * rmag * rmag;

// integration vector

ydot(1) := y(4);

ydot(2) := y(5);

ydot(3) := y(6);

ydot(4) := -emu * y(1) / rcubed;

ydot(5) := -emu * y(2) / rcubed;

ydot(6) := -emu * y(3) / rcubed;

return ydot;

END;

The attached PDF document describes other forms of orbital equations of motion.

The following is the output created by this demo program. A typical elliptical orbit is propagated by one orbital period with a tolerance of 1.0e-10. The differences between the initial and final position and velocity vectors is displayed at the end of the printout.

demo_rkf78

initial classical orbital elements

semimajor axis = 12000 km

eccentricity = 0.15

inclination = 51.9434659014 degrees

argument of perigee = 347.329965384 degrees

raan = 212.885279071 degrees

true anomaly = 347.74343648 degrees

initial eci state vector

eci position vector (km)

r_x = −9233.9364482

r_y = −2805.2805302

r_z = −3395.02813242

eci velocity vector (km/sec)

v_x = −0.169315078972

v_y = −4.61622930329

v_z = 4.83421581734

final eci state vector

eci position vector (km)

r_x = −9233.93644828

r_y = −2805.28052594

r_z = −3395.02813705

eci velocity vector (km/sec)

v_x = −0.16931508232

v_y = −4.61622930428

v_z = 4.83421581607

final classical orbital elements

semimajor axis = 12000.0000006 km

eccentricity = 0.150000000001

inclination = 51.9434659015 degrees

argument of perigee = 347.329965374 degrees

raan = 212.885279071 degrees

true anomaly = 347.743436457 degrees

position vector differences (meters)

delta r_x = −0.00008

delta r_y = 0.00426

delta r_z = −0.00463

velocity vector differences (meters/second)

delta v_x = −0.000003348

delta v_y = −0.00000099

delta v_z = −0.00000127

This demo program solves the system of Keplerian, unperturbed orbital motion. For this example, a measure of the effectiveness of this algorithm can be determined by how well the orbital closes on itself after one complete orbital period.

Here's the source code for the equations of motion implemented in this example.

EXPORT keqm (t, y)

// first order equations of orbital motion

// Keplerian motion - no perturbations

// input

// t = simulation time (seconds)

// y = state vector

// output

// ydot = integration vector

/////////////////////////////

BEGIN

LOCAL rmag, rcubed;

LOCAL ydot := [0, 0, 0, 0, 0, 0];

rmag := sqrt(y(1) * y(1) + y(2) * y(2) + y(3) * y(3));

rcubed := rmag * rmag * rmag;

// integration vector

ydot(1) := y(4);

ydot(2) := y(5);

ydot(3) := y(6);

ydot(4) := -emu * y(1) / rcubed;

ydot(5) := -emu * y(2) / rcubed;

ydot(6) := -emu * y(3) / rcubed;

return ydot;

END;

The attached PDF document describes other forms of orbital equations of motion.

The following is the output created by this demo program. A typical elliptical orbit is propagated by one orbital period with a tolerance of 1.0e-10. The differences between the initial and final position and velocity vectors is displayed at the end of the printout.

demo_rkf78

initial classical orbital elements

semimajor axis = 12000 km

eccentricity = 0.15

inclination = 51.9434659014 degrees

argument of perigee = 347.329965384 degrees

raan = 212.885279071 degrees

true anomaly = 347.74343648 degrees

initial eci state vector

eci position vector (km)

r_x = −9233.9364482

r_y = −2805.2805302

r_z = −3395.02813242

eci velocity vector (km/sec)

v_x = −0.169315078972

v_y = −4.61622930329

v_z = 4.83421581734

final eci state vector

eci position vector (km)

r_x = −9233.93644828

r_y = −2805.28052594

r_z = −3395.02813705

eci velocity vector (km/sec)

v_x = −0.16931508232

v_y = −4.61622930428

v_z = 4.83421581607

final classical orbital elements

semimajor axis = 12000.0000006 km

eccentricity = 0.150000000001

inclination = 51.9434659015 degrees

argument of perigee = 347.329965374 degrees

raan = 212.885279071 degrees

true anomaly = 347.743436457 degrees

position vector differences (meters)

delta r_x = −0.00008

delta r_y = 0.00426

delta r_z = −0.00463

velocity vector differences (meters/second)

delta v_x = −0.000003348

delta v_y = −0.00000099

delta v_z = −0.00000127