Python implementations of findE.m, ymdhms2jd.m, findEarth.m, findMars.m

earth_mars_utils.py

"""Python implementations of findE.m, ymdhms2jd.m, findEarth.m, findMars.m

Based on Dr. Hartzell's files. Written by Nate Wunderly.
"""

from numpy import abs, floor, pi, sin, cos, arccos


"""findE.m

function E=findE(M,e)
%%%Find the eccentric anomaly from the mean motion and the eccentricity, M
%%%should be in radians
%%%% created by Dr. Christine Hartzell (2010)
%%%% Distributed to ENAE 404 for use in project assignment

if M>-pi&&M<0 ||M>pi
    E=M-e;
else
    E=M+e;
end

step=(M-E+e*sin(E))/(1-e*cos(E));
while abs(step)>1E-13
    E=E+step;
    step=(M-E+e*sin(E))/(1-e*cos(E));
end
"""


def findE(M, e):
    """Find the eccentric anomaly from the mean motion and the eccentricity, M."""
    if -pi < M < 0 or M > pi:
        E = M-e
    else:
        E = M+e

    step=(M-E+e*sin(E))/(1-e*cos(E))
    while abs(step) > 1e-13:
        E += step
        step = (M-E+e*sin(E)) / (1-e*cos(E))

    return E


"""ymdhms2jd.m

function jd = ymdhms2jd(year, month, day, hour, minute, second)
%
% - converts year, month, day, hour, minute, second to julian date
% - valid for all positive julian dates
%
% jd = ymdhms2jd(year, month, day, hour, minute, second)
%
% 2000 - Rodney L. Anderson
%
day = day + ((second/60 + minute)/60 + hour)/24;

if(month < 3)
  year = year - 1;
  month = month + 12;
end

A = floor(year/100);
B = 2 - A + floor(A/4);

if(year < 1583)
  B = 0;
  if(year == 1582 & month > 9)
    if(day > 14)
      B = 2 - A + floor(A/4);
    end
  end
end
jd = floor(365.25*(year+4716)) + floor(30.6001*(month+1)) + day + B - 1524.5;
"""


def ymdhms2jd(year, month, day, hour, minute, second):
    """converts year, month, day, hour, minute, second to julian date."""
    day = day + ((second / 60 + minute) / 60 + hour) / 24

    if month < 3:
        year = year - 1
        month = month + 12

    A = floor(year / 100)
    B = 2 - A + floor(A / 4)

    if year < 1583:
        B = 0
        if year == 1582 and month > 9:
            if day > 14:
                B = 2 - A + floor(A / 4)

    return floor(365.25 * (year + 4716)) + floor(30.6001 * (month + 1)) + day + B - 1524.5


"""findEarth.m

function[r,v]=findEarth(JD)
%%%%gives the inertial position and velocity of Earth at a given Julian day
%%%% created by Dr. Christine Hartzell (2010)
%%%% Distributed to ENAE 404 for use in project assignment
%%%% update with your own oe2cart code

TDB=(JD-2451545)/36525;
AU=1.4959787E8; %km
mu=132712440018; %km^3/s^2 (this is for the sun)
a=1.000001018*AU;
e=0.01670862- 0.000042037*TDB- 0.0000001236*TDB^2+0.00000000004*TDB^3;
i=0.0130546*TDB- 0.00000931*TDB^2 - 0.000000034*TDB^3;%deg
Omega=174.873174-0.2410908*TDB+0.00004067*TDB^2-0.000001327*TDB^3;%deg
wtilde=102.937348+0.3225557*TDB+0.00015026*TDB^2+0.000000478*TDB^3; %deg
L=100.466449+35999.3728519*TDB- 0.00000568*TDB^2; %deg

M=L-wtilde;
M=rem(M,360);%deg
w=wtilde-Omega;

E=findE(M*pi/180,e);%radians
E_check=rem(E,2*pi);
nu=acos((cos(E)-e)/(1-e*cos(E)));%rad
%%%Check E and nu should be in the same half plane
if (E_check>pi &&nu<pi) |(E_check<pi && nu>pi)
    nu=2*pi-nu;
end

[r,v]=oe2cart(a,e,i,Omega,w,nu*180/pi);
"""


def findEarth(JD, oe2cart):
    """Gives the inertial position and velocity of Earth at a given Julian day.

    Orbital elements are passed into oe2cart in degrees.
    oe2cart should return a tuple of (r, v)

    Returns r and v returned by oe2cart.
    """
    TDB = (JD-2451545)/36525
    AU = 1.4959787E8  # km
    mu = 132712440018  # km^3/s^2 (this is for the sun)

    a = 1.000001018 * AU
    e = 0.01670862 - 0.000042037*TDB - 0.0000001236*TDB**2 + 0.00000000004*TDB**3
    i = 0.0130546*TDB - 0.00000931*TDB**2 - 0.000000034*TDB**3  # deg
    Omega = 174.873174-0.2410908*TDB + 0.00004067*TDB**2 - 0.000001327*TDB**3  # deg
    wtilde = 102.937348 + 0.3225557*TDB + 0.00015026*TDB**2 + 0.000000478*TDB**3  # deg
    L = 100.466449 + 35999.3728519*TDB - 0.00000568*TDB**2  # deg

    M = L - wtilde
    M = M % 360
    w = wtilde - Omega

    E = findE(M*pi/180, e)  # rad
    E_check = E % 2*pi
    nu = arccos((cos(E)-e)/(1-e*cos(E)))  # rad

    # Check E and nu should be in the same half plane
    if (E_check > pi > nu) or (E_check < pi < nu):
        nu = 2*pi-nu

    r, v = oe2cart(a, e, i, Omega, w, nu*180/pi)
    return r, v


"""findMars.m

function[r,v]=findMars(JD)
%%%%gives the inertial position and velocity of Earth at a given Julian day
%%%% created by Dr. Christine Hartzell (2010)
%%%% Distributed to ENAE 404 for use in project assignment
%%%% update with your own oe2cart code
TDB=(JD-2451545)/36525;
AU=1.4959787E8; %km
mu=132712440018; %km^3/s^2 (this is for the sun)

a=1.523679342*AU;
e=0.09340062+0.000090483*TDB- 0.0000000806*TDB^2 - 0.00000000035*TDB^3;
i= 	 1.849726 - 0.0081479*TDB - 0.00002255*TDB^2- 0.000000027*TDB^3;%deg
Omega=49.558093- 0.2949846*TDB- 0.00063993*TDB^2- 0.000002143*TDB^3;%deg
wtilde=336.060234+0.4438898*TDB - 0.00017321*TDB^2+0.000000300*TDB^3; %deg
L=355.433275+19140.2993313*TDB+0.00000261*TDB^2- 0.000000003*TDB^3; %deg

M=L-wtilde;
M=rem(M,360);%deg
w=wtilde-Omega;

E=findE(M*pi/180,e);%radians
E_check=rem(E,2*pi);
nu=acos((cos(E)-e)/(1-e*cos(E)));%rad
%%%Check E and nu should be in the same half plane
if (E_check>pi &&nu<pi) |(E_check<pi && nu>pi)
    nu=2*pi-nu;
end

[r,v]=oe2cart(a,e,i,Omega,w,nu*180/pi);
"""


def findMars(JD, oe2cart):
    """Gives the inertial position and velocity of Earth at a given Julian day.

    Orbital elements are passed into oe2cart in degrees.
    oe2cart should return a tuple of (r, v)

    Returns r and v returned by oe2cart.
    """
    TDB = (JD-2451545) / 36525
    AU = 1.4959787E8  # km
    mu = 132712440018  # km^3/s^2 (this is for the sun)

    a = 1.523679342*AU
    e = 0.09340062 + 0.000090483*TDB - 0.0000000806*TDB**2 - 0.00000000035*TDB**3
    i = 1.849726 - 0.0081479*TDB - 0.00002255*TDB**2 - 0.000000027*TDB**3  # deg
    Omega = 49.558093 - 0.2949846*TDB - 0.00063993*TDB**2 - 0.000002143*TDB**3  # deg
    wtilde = 336.060234 + 0.4438898*TDB - 0.00017321*TDB**2 + 0.000000300*TDB**3  # deg
    L = 355.433275 + 19140.2993313*TDB + 0.00000261*TDB**2 - 0.000000003*TDB**3  # deg

    M = L - wtilde
    M = M % 360  # deg
    w = wtilde - Omega

    E = findE(M*pi/180, e)  # radians
    E_check = E % 2*pi
    nu = arccos((cos(E)-e)/(1-e*cos(E)))  # rad

    # Check E and nu should be in the same half plane
    if (E_check > pi > nu) or (E_check < pi < nu):
        nu=2*pi-nu

    r, v = oe2cart(a, e, i, Omega, w, nu*180/pi)
    return r, v