1
$\begingroup$

I am following Rene Schwarz' formulas and I've verified my calculations work using the converter by Tony Dunn.

Everything is fine when I use a very small number Epsilon for the z-axis values of the state vectors. However, once I switch to fully 2D and set the z value for the state vectors to 0 I can no longer calculate the position of a body at t (seconds). I believe this has something to do with both lan and w becoming undefined/NaN - when this happens I set them to 0 instead. I think this is causing the tranformation step (#6 in "Orbital Elements to State Vectors" link) to fail to properly transform the position/velocity vectors. Unfortunately I have no mathamatical background and am lucky to have gotten this far, so while I believe I know where the problem is I have no idea how to resolve it.

Am I stuck using Epsilon as the state vector z-axis values or is there a better way to do this for 2D orbits? If there is, what are the changes to the formulas at the links provided or what formulas do I need to use instead?

$\endgroup$
1
  • 1
    $\begingroup$ The general solution to making n-vectors work using algorithms designed for (n+1) dimensions is to set the position of the unused elements to 1, not 0. E.g. in your last link, set Rz to 1, and Vz to 0. $\endgroup$ Apr 25 at 1:53

1 Answer 1

2
$\begingroup$

The full solution is to choose a different set of orbital elements, but in this case, you can get away with just tweaking how you calculate the standard ones.

In Step 1a, the specific angular momentum, $\mathbf{\vec{h}}$, is calculated normally. When $\mathbf{\vec{r}}$ and $\mathbf{\dot{\vec{r}}}$ have no $z$ components, then $\mathbf{\vec{h}}$ will have only a $z$ component.

In Step 1b, the eccentricity vector, $\mathbf{\vec{e}}$, is calculated normally, but be sure you pick the right gravitational parameter $\mu$ for the central body you want. The one in the example only works for the sun.

In Step 1c, when $\mathbf{\vec{h}}$ is along the $z$ axis, the cross product with the $z$ axis direction, $\mathbf{\vec{n}}$, is necessarily zero. Equation 3 for true anomaly, $\nu$, is unchanged --- unless you also have a perfectly circular orbit ($e=0$), which is trivial to solve.

In Step 2, the inclination, $i$, equals arccos(1), which is 0, as expected for an orbit in the equatorial plane. That we can ignore.

Step 3 for eccentricity and eccentric anomaly are unchanged.

Step 4 is the only thing we need to change, because when $\mathbf{\vec{n}}=0$, as it does here, these equations (#6) divide by that zero, so you need to do something else. The $\nu$ you found in Step 1c tells you how far around the ellipse the orbiting body is from periapsis.

You can set either $\Omega$ or $\omega$ to zero, but you can't set them both to zero at the same time. Something (which you may choose to call $\Omega$, $\omega$, or any other letter) needs to equal the angular distance from periapsis to whatever reference direction you choose for your coordinate system, like zero degrees longitude and latitude for an ECF frame or the first point of Aries for an ECI frame.

If you are interested in the more complicated way to define different "equinoctial" orbital elements that work better in this instance, I recommend you start with "A Nonsingular Set of Orbit Elements" (Cohen and Hubbard 1962). If that doesn't make it clear, you might have a look at online documentation for OreKit, FreeFlyer, or STK software.

$\endgroup$
4
  • $\begingroup$ There shouldn't be any problem with both $\Omega$ and $\omega$ being zero at the same time; it just means that your periapsis is the same point on your orbit as your Ascending node., and that your line of apsides is parallel to the reference direction. $\endgroup$
    – notovny
    Apr 25 at 0:36
  • $\begingroup$ @notovny It can happen that both are zero, if you choose your line of apsides to lie in the reference direction. I mean that in general, for arbitrary input vectors, there will be some difference between your preferred reference direction and whatever line of apsides the vectors tell you that particular ellipse has, so you need some parameter to keep track of that offset. $\endgroup$
    – Ryan C
    Apr 25 at 0:56
  • 1
    $\begingroup$ @RyanC - While I don't fully understand your second last paragraph I tested out setting Ω to 90° (1.57 rad), that seems to have done the trick. I am plotting points using the orbital elements -> state vectors conversion looking at the system from the top down. Prior to this fix the plot was definitely rotated incorrectly but this has fixed it and the real simulated path of the object matches the plotted one. Can you confirm I have applied the proper modification given that the reference direction is from the top down looking at the orbit from Inf Z to -Inf Z? Sorry I am out of my depth here. $\endgroup$
    – Wafer
    Apr 25 at 2:31
  • $\begingroup$ @RyanC - Oops I didn't describe my reference frame correctly. I'll try to clarify it better. It is looking at the entire flat surface of the x,y plane from the z axis such that y decreases to the north, x decreases to the west, and z increases the closer the reference "viewer" gets towards the x,y plane. $\endgroup$
    – Wafer
    Apr 25 at 2:51

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.