File Exchange

image thumbnail

Kepler State Transition Matrix (MEX)

version 1.3 (133 KB) by

Compute Kepler state transition matrix for an arbitrary number of time steps. M and C++ version.



 The Kepler State transition matrix provides a way to progress any given state vector for a given time step, without having to perform a lengthy triple-coordinate conversion (from Cartesian coordinates to Kepler elements, progressing, and back). This procedure has been generalized and greatly optimized for computational efficiency by Shepperd (1985)[1], who
 used continued fractions to calculate the corresponding solution to Kepler's equation. The resulting algorithm is an order of magnitude faster that mentioned triple coordinate conversion, which makes it very well suited as a basis for a number of other algorithms that require frequent progressing of state vectors.
 This procedure implements a robust version of this algorithm. A small correction to the original version was made: instead of Newton's method to update the Kepler-loop, a Halley-step is used. This change makes the algorithm much more robust while increasing the rate of convergence even further. If the algorithm fails however (which SHOULD never occur), the slower triple-coordinate conversion is automatically started.
Both a M-version and a C++-version have been included in this submission. The C++ version can be compiled by issuing
mex progress_orbit.c

at the MATLAB command prompt (when in the correct directory). After compiling (let me know what goes wrong) you can run a small demo (test_STM.m) which will produce something similar to the screenshot on top here; a plot of the factor by which the MEX version outperforms the M-version. On my computer this was 5.5 times faster on average, sometimes peaking to no less than 260 times faster.

[1] S.W. Shepperd, "Universal Keplerian State Transition Matrix". Celestial Mechanics 35(1985) pp. 129--144, DOI: 0008-8714/85.15

If you find this work useful, please consider a donation:

Comments and Ratings (4)

Dear Oldenhuis, the code is nice.

I would like to point out the following (which is confusing for me) and read your opinion about it.

In control theory the state-transition matrix is a matrix whose product with the state vector x at an initial time t_0 gives x at a later time t. (i.e. spacecraft/satellite orbit propagation in this case )

However, in Shepperd (1985), it appears to me that the state-transition matrix is not the one which propagate the orbit, but is the one containing the first-order sensitivities of the spacecraft's position and velocity vectors at some point t with respect to variations in the initial position/velocity vector at some previous time t0. (i.e. is the matrix which propagate the state deviations in position and velocity)
Eventually, the calculation of this matrix incidentally also provides a robust method for the Keplerian propagation of spacecraft/satellite.
In fact, the combination of successive State Transition Matrix (as defined by Shepperd) is used in software like GALLOP or EMTG to calculate the derivatives required to solve the two-point boundary value problem formed from the necessary conditions for optimality.

Is it a notation error? I Hope you can help me understand.

Guangjun Liu

Rody Oldenhuis

Rody Oldenhuis (view profile)

Dear Necat,

I've read your ideas and proofs on ScienceRay ("The Orbit of The Planets", "Are Kepler's Laws Right?" and "Is Kepler Right"). I'm sorry to disappoint you, but there are a few mistakes in the proofs you outline:

- The formula you use (h=-1/2*g*t^2+Vh0*t+h0) only applies in a homogeneous gravitational field (e.g., g=constant), but this is only a good approximation when you're near the surface of a planet; the further away you get, the weaker a planet's gravitational pull (g=g(r), r=distance to planet). Moreover, the formula (and the parabola it describes) only results when you assume a flat Earth; but we live on a spherical Earth. The way in which you apply the formula does not take either of these things into account, hence, you can not use it in this way. The *tricky* math you describe is actually the correct way of doing it.

- You transform coordinates from Cartesian to Polar, and continue by interpreting the polar plot as if it describes Cartesian coordinates; you simply confuse what both plots represent. For example, you can make the formula y=x look like a spiral, but in that polar representation it just means y=x would *look* like a spiral if you (the observer) were rotating; which usually doesn't make anything more clear.

Moreover, if what you propose is indeed true, isn't it then a bit odd that the Sun always looks the same from the Earth? Your theory would imply that the Sun would constantly change its apparent size, since the Earth would continually be moving closer and farther from the Sun as time progresses. This would of course also be true for the Moon. Clearly, this is not the case. How does your theory explain seasons, for example? Or the fact that anyone with binoculars and patience can see Jupiter's moons move in a near-circular trajectory around Jupiter? Or the fact that the GPS system is able to determine your position within a few millimeters if you want, while the GPS system is based on Newton's laws (and then some)?

The overwhelming amount of observational evidence clearly is in Kepler's/Newton's favor. Your theory simply makes predictions that do not agree with observation, hence, the theory is incorrect (or flawed at least).

Will you please remember that Kepler's laws are wrong: the orbits are not elliptical,the Sun is not at a focus of such ellipses,no eccentricity,no aphelion,no perihelion,no equality of swept out areas in equal time,...etc.Then what are you doing before considering such arguments? Please See a recapitulation on Kepler's laws right?/dated 2010.01.09



[linked to Github]


Updated contact info (this time, WITH the file!)


Updated contact info

MATLAB Release
MATLAB 7.9 (R2009b)

Download apps, toolboxes, and other File Exchange content using Add-On Explorer in MATLAB.

» Watch video