Nodal period of an Earth satellite - numerical integration solution

[cdeaglejr]

This PPL program calculates the nodal period of a satellite using a numerical integration technique. The algorithm first uses the RKF78 method to propagate the first-order perturbed equations of satellite motion forward in time to bracket a "nodal objective function". Once a bracket is found, the algorithm then uses Brent's root-finding method to find the time at which the objective function is zero to within a user-defined tolerance. While looking for the root the software also uses the RKF78 method to propagate the satellite's orbit subject to the main oblateness gravity perturbation of the Earth.

The objective function is bracketed whenever the value of the function changes sign. It may be negative at the left bracket and positive at the right bracket, for example.

The objective function for this example is the instantaneous declination of the satellite as it moves along its orbit (r_z / |r|). The root-finding process also ensures that the z-component of the satellite's velocity vector is positive at the root. When the declination is zero and the z-component is positive, the satellite is at the ascending node. This can be verified by examining the z-component of the position vector, which should be very nearly zero, and the z-component of the velocity vector at the nodal crossing. This program provides this information to the user.

The user can provide the classical orbital elements starting at line 127 of the source code.

// ----------------------------------
// initial classical orbital elements
// ----------------------------------

// semimajor axis (kilometers)

oev(1) := 8000.0;

// orbital eccentricity (non-dimensional)

oev(2) := 0.015;

// orbital inclinations (degrees)

oev(3) := 28.5 * dtr;

// argument of perigee (degrees)

oev(4) := 270.0 * dtr;

Can you modify the following pfunc subroutine to calculate the anomalistic or sidereal orbital period? The attached PDF document describes how this can be done.

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

pfunc(x)

// declination objective function

// Orbital Mechanics with the HP Prime

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

BEGIN

LOCAL h, tetol, rmag, fx;

// step size guess (seconds)

h := 10.0;

// rkf78 truncation error tolerance

tetol := 1.0e-10;

// integrate from tiwrk to requested time = x

yfinal := rkf78(tiwrk, x, h, tetol, yi);

// compute current value of geocentric declination (radians)

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

fx := yfinal(3) / rmag;

return fx;

end;