How to make an algorithm that finds time to a true anomaly at epoch?

Question

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."

Answer 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!