0

I am having trouble creating a function that finds the time for an inputed true anomaly when I know my current eccentric and mean anomaly (my equations are in vector form).

Note that I only intend to use this function when the orbit is an ellipse (so when eccentricity is less than 1 ) Only standard SI units are used.

I am trying to express all of my angles (that includes true, eccentric and mean anomalies in the ranges of 0 to 2π radians)

The problem is that for the true, mean and eccentric anomalies you need to do logic checks to convert the angles into the proper ranges (that being from 0 to 2π radians).

Before reading further note that :

Mu = 6.67E-11 * Mass_of_Earth 
Mean_orbital_motion = sqrt(pow(abs(Semi_major_axis),3)/Mu)

$\mu = G M_{earth} = 6.67480 \times 10^{-11} \cdot 5.972 \times 10^{24}$

$n =\sqrt\frac{\mid a^{3}\mid}{\mu} $

My True anomaly is defined as follows :

Current_True_anomaly = if dot(v:Position,v:Velocity)>=0 then acos(dot(v:Eccentricity_vector,v:Position)/(magnitude(v:Eccentricity_vector)*magnitude(v:Position))) else (2*(pi)) - acos(dot(v:Eccentricity_vector,v:Position)/(magnitude(v:Eccentricity_vector)*magnitude(v:Position)))
# The "if statement" checks to see if we have past the periapsis 
# If we have passed it then the dot(v:Position,Velocity) would be less than 0.
# Or else if dot(v:Position,Velocity) would be greater than 0 if we haven't passed it yet.

My current eccentric anomaly :


Current_Eccentric_anomaly = if dot(v:Position,v:Velocity)>=0 then  2*atan(sqrt(1-Eccentricity/1+Eccentricity)*tan(Current_true_anomaly/2)) else (2*(pi) + 2*atan(sqrt(1-Eccentricity/1+Eccentricity)*tan(Current_true_anomaly/2))
# The "if statement" checks to see if we have passed peripasis.
# It will do the correct adjustments to fix the angle.

Here is the pseudocode for the function:

def Time_to(true_anomaly):
    Eccentric_anomaly_input =  if true_anomaly>0 then 2*atan(sqrt(1-Eccentricity/1+Eccentricity)*tan(true_anomaly/2)) else (2*pi) + 2*atan(sqrt(1-Eccentricity/1+Eccentricity)*tan(true_anomaly/2))
    #The if statement here (if true_anomaly>0) is to see if we have passed the periapsis.

Mean_anomaly_input = Eccentric_anomaly_input - Eccentricity*(sin(Eccentric_anomaly_input)

Current_mean_anomaly_time = Current_mean_anomaly * Mean_orbital_motion
Mean_anomaly_time_input
# Converting the mean anomalies into the form of time since periapsis
Difference = Mean_aonamly_time_input - Current_mean_anomaly_time

time_to_true_anomaly = if Difference<0 then Difference + Orbital_period else Difference
#logic check to see if our inputted true anomaly is past the periapsis.

The problem is that this function does not work at all. I'm not too sure why it breaks down. If anyone could give a solution that keeps the anomalies in the ranges of 0 to 2π radians that would be appreciated.

{edit} The questioner left a comment on Jan 22, 2022, reading "Thankfully I have found a solution. I will update my post soon."

terry-s
  • 1,102
  • 8
  • 15
John
  • 113
  • 4
  • You may find some things in this answer to What is the analytical closed-form solution of the two-body problem to verify its numerical integration results? potentially helpful. I'm not good at debugging software so I can't tell what's wrong here. Are you sure your math is correct? If you work the problem by hand, say with a calculator, do you get the right answer? That would suggest a typo. If no, then write your equations out as math rather than code and then it will be easier for others to read them and help look for the problem. – uhoh Jan 17 '22 at 22:26
  • 1
    I will also need to double check my results since I have made a few modifications since my last execution .MathJax is basically like Latex? If that is the case then it should not be too hard to turn my code into understandable equations. Though I will take some time so I will get it done latter during the day. – John Jan 18 '22 at 04:37
  • yes MathJax is a lot like LaTeX so if you are familiar with that it should be easy. Thanks! – uhoh Jan 18 '22 at 04:43
  • 1
    Yep that Is fine. I am just a bit worried that because I left it quite late my question might not be answered at all. – John Jan 18 '22 at 04:49
  • Your value for Earth's standard gravitational parameter $\mu = GM$ of about 3.986E+14 m^3/s^2 seems right, is your semimajor axis in meters as well, or might you have used kilometers? In meters it should be about 1.5E+11. – uhoh Jan 18 '22 at 09:55
  • 3
    When you update, please also provide a set of input values that return a result you find incorrect. It will be much easier to diagnose what's going wrong with your calculation if we know the values of the position vector, velocity vector, eccentricity vector, and eccentricity that you are feeding your calculations. – notovny Jan 18 '22 at 09:58
  • 4
    As a minor note, do not use $\mu = GM$. The masses of the Sun and the planets are determined by $M_b = \mu_b/G$, where the subscript $_b$ denotes the body in question. It's much better to use the appropriate standard gravitational parameter value on the linked wikipedia page. The values aren't up to date, but they're a lot better than using \mu=GM. – David Hammen Jan 18 '22 at 15:00
  • 3
    A much more severe problem is that the very long lines make your code hard to read, so hard that I cannot see the bug. There's a reason for the Python rule of 79 characters, max, per line. It's a human factors thing. In a code review (one of the things I do for a living) I would have sent this back as unreadable; try again. – David Hammen Jan 18 '22 at 15:12
  • In addition to the post uhoh linked to, I think this post : https://space.stackexchange.com/questions/54448/calculating-the-time-of-flight-between-two-anomalies-gives-a-negative-result may be useful. It deals with calculating time of flight between two true anomalies, using eccentric anomalies; which (correct me if I'm wrong) looks close to what you are dealing with. – Krafpy Jan 18 '22 at 20:33
  • Sorry for the late reply. My week has been relatively busy. Thankfully I have found a solution. I will update my post soon. – John Jan 22 '22 at 08:04

1 Answers1

1

sounds like studying these two libraries could be helpful for you:

https://github.com/dfm/kepler.py

{edit:} ('kepler.py' links to a small and fairly easily accessible piece of software billed as a "Fast and stable solver for Kepler's equation extracted from 'exoplanet'.")

https://esa.github.io/pykep/

The latter is from the European Space Agency.

{edit:} 'pykep' is a very large software collection, and to get to a place where the orbital equation solvers actually 'live', one has to visit the code installation page at (https://github.com/esa/pykep.git). There's over 30 MB of software in there, and two or three of the routines are useful for the elliptical orbital calculations.

Enjoy!

terry-s
  • 1,102
  • 8
  • 15
Matteo Manzi
  • 11
  • 1
  • Your answer could be improved with additional supporting information. Please [edit] to add further details, such as citations or documentation, so that others can confirm that your answer is correct. You can find more information on how to write good answers in the help center. – Community Jul 10 '22 at 12:30
  • @Matteo Manzi : I'm sorry you got downvoted by the 'evil' community-bot, especially as you're a new user who did actually put up a couple of relevant links. It is true, though, that a little more detail in the information could have been even more helpful. Those modern packages on github and elsewhere have masses of software housekeeping, they are a big forest of trees in which it is difficult to see the relevant part. It would be useful to pinpoint the particular routines of relevance here. I did look myself at the two links, the results were surprising, I will post them in an answer. – terry-s Jul 12 '22 at 21:40
  • @Matteo Manzi : I also took the liberty of offering some extra information about the software to which you gave links. I hope that is helpful. – terry-s Jul 12 '22 at 22:03