09-16-2022, 04:49 PM
This HP Prime program demonstrates how to interact with an rkf78 subroutine which solves first-order systems of ordinary differential equations. This is an HPPL implementation of the Runge-Kutta-Fehlberg method of order 7 with an 8th order error estimate. This estimate is used to determine the variable step size required to satisfy a user-defined "convergence" criterion.
Here's the source code for a typical system of first-order differential equations included in this demo program. This code implements Keplerian (unperturbed) orbital motion.
keqm (t, y)
// first order equations of orbital motion
// Keplerian motion - no perturbations
// input
// t = simulation time (seconds)
// y = state vector (kilometers and kilometers/second)
// output
// ydot = integration vector (kilometers/second/second)
/////////////////////////////
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 main program defines typical orbital elements of an Earth satellite and propagates the orbit for 10 Keplerian periods. The errors between the initial and final conditions indicate how well the algorithm performs.
Here's the source code for a typical system of first-order differential equations included in this demo program. This code implements Keplerian (unperturbed) orbital motion.
keqm (t, y)
// first order equations of orbital motion
// Keplerian motion - no perturbations
// input
// t = simulation time (seconds)
// y = state vector (kilometers and kilometers/second)
// output
// ydot = integration vector (kilometers/second/second)
/////////////////////////////
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 main program defines typical orbital elements of an Earth satellite and propagates the orbit for 10 Keplerian periods. The errors between the initial and final conditions indicate how well the algorithm performs.