6

I have developed an orbit propagator, taking J2 perturbation into account according to the formulation as shown: acceleration equation with Runge-Kutta 4th order, timestep of 1 second as the integrator. Formulation as shown:

Formulation

With J2 = 0.0010826, Re = 6.378137E+6 and mu = 3.986004418000000e+14.

Subsequently, I tried to compared its orbit propagation accuracy with SGP4 propagator as well as the 2 Body propagator and I found out that the position error between "SGP4" and "Orbit Propagator with J2" is much larger compared to the position error between "SGP4" and "2 Body propagator".

Some of the details of the orbit propagation simulation are:

  1. TLE used for SGP4 propagator: enter image description here

  2. Propagation duration of 16 hours

  3. As the output of SGP4 is in TEME frame, there have been converted into J2000 frame when comparing the propagation error.

  4. The initial position and velocity for the "orbit propagator with J2" and the "2 Body propagator" is obtain from the initial position and velocity output of SGP4 converted to J2000 frame.

  5. SGP4 is a function from Matlab Aerospace toolbox

The position error in cartesian coordinates, with respect to J2000 is as shown: x y z

I have an impression that orbit propagation by taking J2 perturbation into account should be more accurate compared to 2 Body propagator and thus I am wondering if I have made a mistake somewhere? Or is there a possibility that introducing J2 perturbation will induce more error? Any help/advice/sharing based on your experience is much appreciated!

Chia Jiun Wei
  • 91
  • 6
  • 1
    *Cool!* Just a thought; what happens if you change the sign of your $J_2$? What is the sign of Earth's J2? You can also compare your acceleration equations to mine. The problem with gravitational potential expressions is that there can be multiplicative factor variations from one source to another that can throw monkey wrenches into calculation if one doesn't read very carefully and check intermediate steps. – uhoh Oct 21 '21 at 02:54
  • See this answer to Trouble deriving rectangular components of acceleration of satellite in orbit around Earth with J2 consideration and also see For the mathematical relationship between J2 (km^5/s^2) and dimensionless J2 - which one is derived from the other? If you can resolve your problem then please feel free to post an answer to your question. Thanks! – uhoh Oct 21 '21 at 02:55
  • 2
    You haven't shown your code, so there's no way to truly tell. One thing that you might be doing wrong is your use of the J2000 frame. The ECEF frame is becoming tilted with result to J2000 due to precession and nutation. While you do not need to model the Earth's rotation when modeling Earth gravity with J2, you do need to model the Earth's precession and nutation. – David Hammen Oct 21 '21 at 04:15
  • 1
    @uhoh Look at the OP's equations. $J_2$ obviously has to be unitless. $J_2$ is positive, unitless, and has a value of about 10$^{-3}$, 0.0010826359 to be precise. – David Hammen Oct 21 '21 at 06:24
  • 1
    Thanks for your suggestion. I have added my formulation and some extra information. As I am not an expert in this field, I will take a little more time to study the modelling of Earth's precession and nutation and see if I can reduce the error. – Chia Jiun Wei Oct 21 '21 at 07:12
  • 1
    There seems to be an error in your correction term for gravity with J2 (first Equation). The (a/r) is squared twice. – Ng Ph Oct 21 '21 at 09:14
  • 1
    @NgPh I'm fairly sure that the expression for gravitational acceleration is correct. It's not normally written that way, but it is correct. A more typical expression for the J2 perturbation (for example equation (3) in this masters degree thesis) is $$\vec p = \frac32 \frac{\mu J_2 a^2}{r^5} \left(\left(5\frac{z^2}{r^2}-1\right)(x\hat x + y \hat y) + \left(5\frac{z^2}{r^2}-3\right) z\hat z\right)$$ – David Hammen Oct 21 '21 at 13:25
  • 1
    @NgPh The error in the potential is not replicated in the gravitational acceleration equations. – David Hammen Oct 21 '21 at 13:26
  • 1
    You've shown us your code. Have you unit tested it? A nice unit test is to see if you can replicate gravitational acceleration on the surface of the Earth. You can get some of the values needed for such a unit test from this answer at Physics.SE. Things left out of that answer are the polar radius, which is easy given the equatorial radius and the flattening, and modifying $\mu$ to exclude the mass of the atmosphere, which is not so easy, but still doable. It is always good to unit test your code. – David Hammen Oct 21 '21 at 13:40
  • 1
    @David Hammen, I see. So, the 1st equation (for the potential $V_{J2}$) is wrong, but the equations for $g_i$ are correct (for which I agree). – Ng Ph Oct 21 '21 at 15:45
  • Just fyi normally we post a code block as text, not as a screenshot. This is for several reasons. Just copy your code, make sure it is all at least four spaces indented and paste it like I did here. – uhoh Oct 21 '21 at 22:18
  • In your plots of errors of 2body vs SGP4, I notice that the errors seems to swing around 0. This is strange to me since J2 would make your orbit plane to precess. I would expect that the X-axis error to grow linearly. Are you confident of your plotting routine? – Ng Ph Oct 22 '21 at 10:35

1 Answers1

3

Thanks everyone for your help and advice!

After some troubleshooting, I found out that the large position error of the "orbit propagator with J2 perturbation" is due to the bad initial position and velocity.

Apparently the initial position and velocity at TLE epoch time generated from the MATLAB Aerospace toolbox SGP4 is off by a few kilometres, hence the large propagation error when it is used in the "orbit propagator with J2 perturbation".

I have downloaded David Vallado's SGP4 code from here SGP4 reference code and use the initial PV generated from it for the "orbit propagator with J2 perturbation" as well as the "2Body Propagator". The position error comparison in all 3-axis is as shown:

x y z

Special thanks to Dr S.T.Goh from NUS STAR.

Chia Jiun Wei
  • 91
  • 6
  • Congrats, nice work! – uhoh Oct 22 '21 at 10:45
  • @uhoh, do we have any clue that the displayed numerical results are correct? – Ng Ph Oct 23 '21 at 20:17
  • 1
    @NgPh As an aside, no orbit propagation is "correct", but if 2 body + J2 agrees with SGP4 for 48 hours within a few km, that's pretty darn good! TLE+SGP4 is already known to only be good to circa 2 kilometers per day (https://space.stackexchange.com/a/29817/12102 and https://space.stackexchange.com/a/29783/12102) It contains coefficients up to (I think) 6th or 8th order in gravity, but does not numerically propagate at all. Instead it just evolves the mean orbital elements. So we don't know how much of the tiny disagreement shown is due to either one. – uhoh Oct 23 '21 at 22:04
  • @uhoh, then too bad. Only the OP understood the problem. Only the OP can understand the solution. Only the OP can produce the plots. This question and its answer has zero value to the Space SE community. Although, if somebody takes the OP displayed numerical results at face value, he may make conclusions (at his risk and peril). – Ng Ph Oct 23 '21 at 22:05
  • @NgPh well I wouldn't presume to know what other readers can glean from the OP's work, but all I can say is Welcome to Stack Exchange! If the OP wants someone to check their work they can post a new question. – uhoh Oct 23 '21 at 22:07
  • 1
    @uho, I just wanted to warn that the results are not verifiable, don't make hasty conclusions (such as J2-only matches SPG4 up to x days). – Ng Ph Oct 23 '21 at 22:11
  • 1
    @NgPh This post was meant to find out why does the "Orbit Propagator with J2 perturbation" has larger error compared with simple "2 body Propagator" using SGP4 as a baseline. From my understanding, taking J2 perturbation into account should yield better accuracy compared to 2 body. And the answer for it was that I used an initial position and velocity that is off by a few kilometres, hence the large propagation error. The numerical accuracy of the "orbit propagator with J2 perturbation" was not discussed at all and no conclusion about the performance is even mentioned. – Chia Jiun Wei Oct 24 '21 at 01:01
  • You observed unexplainable results and you asked SE why. You later found a discrepancy in the settings. I agree with you that it is enough to close the original question, independently of your new results. I agree that I should not draw any conclusion based on the new results (I don't). Nevertheless, they seem to display the trend that the 2body position comes close to the SGP4 position every orbital period. Also the amplitude of the "swings" of the 2body seems to grow faster in the latest set of results. Care to explain why? – Ng Ph Oct 24 '21 at 09:14
  • @uhoh, my concept of a place like SE is sharing knowledge. It would be sad that Stack Exchange is used as a place for debugging somebody's code. Techniques for checking/debugging, avoiding known pitfalls in certain programming/reasoning, ... on the other hand are knowledge. – Ng Ph Oct 24 '21 at 12:32
  • @NgPh when we post questions, it's often (but not always) because we simply don't know the answer. Perhaps 1% of the questions about orbit simulation here turn out to be something simple (I recall one patched conic question where the OP simply forgot to turn on Mars' gravity at patch point). We don't know ahead of time if a question is going to be of that 1% variety or not. The OP did exactly the right thing. They figured out what the problem was, it was not an issue of orbital mechanics, so they wrote a short answer to inform all what's going on. Issue closed. – uhoh Oct 24 '21 at 22:33
  • @NgPh "...It would be sad that Stack Exchange is used as a place for debugging somebody's code..." Yes, if 30% of all questions were about code debugging, that would have to be addressed. But there is nothing sad about a once-in-a-while phenomenon. There's no basis for you to keep asking them to provide you something more. If you're genuinely curious about SGP4 vs 2-body vs 2-body+J2, then post a new SE question about orbital mechanics, include a link back here for background, and let anybody answer. That's better than going after a brand new user and potentially making them feel unwelcome – uhoh Oct 24 '21 at 22:37
  • @ChiaJiunWei this comment chain is a little unusual and not at all the typical reception here in Space SE; I think you can safely ignore it. Welcome to Space! :-) – uhoh Oct 24 '21 at 22:42
  • @uhoh, I do not have the authority to qualify which comments are "unusual". But in case you were implying that I was disparaging what the OP did, then you are wrong. If you assume that I am disputing the correctness of the OP's new results, you are wrong. As you think I am requesting explanations, you are dead wrong. I don't like polemics, so I stop here, by concluding that you and I do not share the same opinion on how SE works, and that's Ok for me. – Ng Ph Oct 25 '21 at 09:53