Ecliptic & Equatorial (ICRF) and position attributes/methods for r_vec & v_vec

[beluej123]

My overall question is, what is your suggestion to convert between Ecliptic and Equatorial xyz frames (xyz & velocity); I read thru the documentation, but it is not clear, to me the most straight forward techniques. BTW, thank you for https://rhodesmill.org/skyfield/design.html , and of course https://rhodesmill.org/skyfield/coordinates.html , but it is not clear to me the "smartest work-flow" to convert: Ecliptic->Equatorial AND Equatorial->Ecliptic (xyz for position and velocity).
FYI: I work mostly in Ecliptic xyz for r_vectors and v_vectors (mission planning and the like).

1. I find a skyfield position object has methods ecliptic_position(), ecliptic_velocity(), but they do NOT give me the same results as position.frame_xyz_and_velocity(ecliptic_frame). The differences were big enough to make we wonder where ecliptic_position() came from - I could not figure it out from either the code or the documentation. When I look at the position attributes/methods with the python command dir(position) I see a bunch of methods/attributes that are not clear to me... However, I would like your suggestion on a smart work-flow to go back an forth between Ecliptic/Equatorial.
   Thank you.
   In addition to my code, below that I show my output. I hope it makes sense to you.

# Calculate earth departure & venus arrival; prepare for plots.
#   Explore conversions between equatorial (ICRF) and ecliptic planes.
import numpy as np
from skyfield.api import Loader  # to set my folder for ephemeris files
from skyfield.framelib import ecliptic_frame
# prevent multiple copies of large data files, get ephemerous data from common folder
#   https://rhodesmill.org/skyfield/api-iokit.html
ephem_dir = r"~\Documents\!Data\OrbMech\!Orb_data\spice_data"
load = Loader(ephem_dir)
planets = load("de421.bsp")
# JPL ephemeris files: https://rhodesmill.org/skyfield/planets.html
#   DE421 valid dates 1899-2053, small file, 17MB, circa 2008
#   DE422 valid dates -3000 to 3000, medium file, 623MB, circa 2008
#   DE441 valid dates -13200 to 17191, huge file 3.1GB, circa 2020
#   DE440s valid dates 1849 to 2150, smallish file 32MB, circa 2020
earth = planets["earth"]
venus = planets["venus"]
ts = load.timescale()
t1 = ts.utc(1988, 4, 8, 0, 0)  # launch time
t2 = ts.utc(1988, 7, 26, 0, 0)  # arrive time

# calculate planet positions at t1, t2
posE_launch = earth.at(t1)  # pos earth at launch, ICRF equatorial
posE_arrive = earth.at(t2)  # pos earth at arrive, ICRF equatorial
posV_launch = venus.at(t1)  # pos venus at launch, ICRF equatorial
posV_arrive = venus.at(t2)  # pos venus at arrive, ICRF equatorial
# print(f"dir(posE_launch):\n{dir(posE_launch)}")

# earth ICRF equatorial (default)
np.set_printoptions(precision=7, floatmode="fixed")
print(f"t1: {t1.utc_strftime()}")
print(f"  Earth launch, ICRF position (equatorial frame):")
print(f"    posE r_vec {posE_launch.xyz.au} [au]")
print(f"    posE r_vec {posE_launch.xyz.km} [km]")
print(f"    earth r_mag: {np.linalg.norm(posE_launch.xyz.au):.6g} [au]")
print(f"    posE v_vec {posE_launch.velocity.km_per_s} [km/s]")
print(f"    posE v_vec {posE_launch.velocity.au_per_d} [au/day]")
# earth ICRF ecliptic (need to convert to ecliptic)
print(f"  Earth launch, ICRF position (ecliptic frame):")
print(f"    ?1 posE r_vec {posE_launch.frame_xyz(ecliptic_frame).au} [au]")
print(f"    ?2 posE r_vec {posE_launch.frame_xyz_and_velocity(ecliptic_frame)[0].au} [au]")
print(f"    ?3 posE r_vec {posE_launch.ecliptic_position().au} [au]")
print(f"") # this line for output readability
print(f"    ?1 posE v_vec {posE_launch.frame_xyz_and_velocity(ecliptic_frame)[1].au_per_d} [au/day]")
print(f"    ?2 posE v_vec {posE_launch.ecliptic_velocity().au_per_d} [au/day]")

# venus ICRF equatorial (default)
print(f"\nt2: {t2.utc_strftime()}")
print(f"  Venus arrive, ICRF position (equatorial frame):")
print(f"    posV r_vec {posV_arrive.xyz.au} [au]")
print(f"    posV r_vec {posV_arrive.xyz.km} [km]")
print(f"    posV r_mag: {np.linalg.norm(posV_arrive.xyz.au):.6g} [au]")
print(f"  Venus arrive, ICRF position (ecliptic frame):")
print(f"    ?1 posV r_vec {posV_arrive.frame_xyz(ecliptic_frame).au} [au]")
print(f"    ?2 posV r_vec {posV_arrive.frame_xyz_and_velocity(ecliptic_frame)[0].au} [au]")
print(f"    ?3 posV r_vec {posV_arrive.ecliptic_position().au} [au]")

# angle between earth and venus at t1 (equatorial)
posE_r_mag_t1 = np.linalg.norm(posE_launch.xyz.au) # assume equatorial
posV_r_mag_t1 = np.linalg.norm(posV_launch.xyz.au)
ang_E_V_t1= np.arccos(np.dot(posE_launch.xyz.au, posV_launch.xyz.au)
                   / (posE_r_mag_t1 * posV_r_mag_t1))
print(f"\nAngle earth->venus (t1) (equatorial)")
print(f"  {np.degrees(ang_E_V_t1)} [deg]")

# angle between earth and venus at t1 (ecliptic)
posE_r_mag_t1 = np.linalg.norm(posE_launch.frame_xyz(ecliptic_frame).au)
posV_r_mag_t1 = np.linalg.norm(posV_launch.frame_xyz(ecliptic_frame).au)
ang_E_V_t1= np.arccos(np.dot(posE_launch.xyz.au, posV_launch.xyz.au)
                   / (posE_r_mag_t1 * posV_r_mag_t1))
print(f"Angle earth->venus (t1) (ecliptic)")
print(f"  {np.degrees(ang_E_V_t1)} [deg]")

Code output:

t1: 1988-04-08 00:00:00 UTC
  Earth launch, ICRF position (equatorial frame):
    posE r_vec [-0.9530037 -0.2902019 -0.1258114] [au]
    posE r_vec [-1.4256733e+08 -4.3413585e+07 -1.8821120e+07] [km]
    earth r_mag: 1.00412 [au]
    posE v_vec [  8.9820136 -26.0225478 -11.2830480] [km/s]
    posE v_vec [ 0.0051875 -0.0150293 -0.0065165] [au/day]
  Earth launch, ICRF position (ecliptic frame):
    ?1 posE r_vec [-9.5390212e-01 -3.1358012e-01  1.6693759e-05] [au]
    ?2 posE r_vec [-9.5390212e-01 -3.1358012e-01  1.6693759e-05] [au]
    ?3 posE r_vec [-9.5300372e-01 -3.1629994e-01  5.9630137e-06] [au]

    ?1 posE v_vec [ 5.1407969e-03 -1.6395942e-02 -5.2450942e-08] [au/day]
    ?2 posE v_vec [ 5.1875470e-03 -1.6381211e-02 -4.7305387e-07] [au/day]

t2: 1988-07-26 00:00:00 UTC
  Venus arrive, ICRF position (equatorial frame):
    posV r_vec [ 0.6286743 -0.3114309 -0.1801112] [au]
    posV r_vec [ 94048335.2777013 -46589395.7963608 -26944256.9570653] [km]
    posV r_mag: 0.724335 [au]
  Venus arrive, ICRF position (ecliptic frame):
    ?1 posV r_vec [ 0.6276841 -0.3591135 -0.0413610] [au]
    ?2 posV r_vec [ 0.6276841 -0.3591135 -0.0413610] [au]
    ?3 posV r_vec [ 0.6286743 -0.3573764 -0.0413687] [au]

Angle earth->venus (t1) (equatorial)
  42.15031656468812 [deg]
Angle earth->venus (t1) (ecliptic)
  42.150316564688055 [deg]

[brandon-rhodes]

Your calls to .frame_xyz(ecliptic_frame) and .frame_xyz_and_velocity(ecliptic_frame) are fine, and indeed match the documentation at:

https://rhodesmill.org/skyfield/coordinates.html#ecliptic-coordinates-1

The old ecliptic_position() method was deprecated and removed from the documentation about a half-decade ago because of its confusing behavior: it defaults to the ecliptic of J2000, unless given a time argument to select a different epoch. I suggest ignoring that method.

It looks like I should take this opportunity to expand the documentation at the bottom of this section of the docs:

https://rhodesmill.org/skyfield/coordinates.html#turning-coordinates-into-a-position

The call you want to convert the other way is ICRF.from_time_and_frame_vectors(t, ecliptic_frame, r, v). See if it works for you!

[beluej123]

To answer your remark to see how the conversions work for me.
It will take me a day or two to work thru a couple sanity checks. I will probably choose a couple of [4] Vallado examples; maybe chap 2.5 and 2.6, Kepler; or chap 3.4, coordinates; or chap 12 (example 12-8), interplanetary. Or [3] Curtis maybe chap 4 orbits in 3D (ex. 4.11), or chap 8 interplanetary trajectories (starting with ex. 8.4 continuing thru ex. 8.8).

FYI:
This effort began as I wanted sanity checks against various orbital calculations, including patched conic estimations in [1] Bate, Mueller, White and [3] Curtis and [4] Vallado and [8] Brown (you may note my references list below). In planetary orbital plots I noticed the orbital plane was tilted, which reminded me that my position coordinates in the equatorial plane...

If you are curious, I keep a references list shown below, of course there are others not listed that were not directly used for this skyfield effort. My go-to reference is usually [4] Vallado or [1] BMWS (Horizons for a "gold standard", for now).

References collection that were used somewhere in ...
References:
----------
    [1] BMWS; Bate, R. R., Mueller, D. D., White, J. E., & Saylor, W. W. (2020, 2nd ed.),
        Fundamentals of Astrodynamics. Dover Publications Inc.
    [2] Vallado, David A., (2013, 4th ed.),
        Fundamentals of Astrodynamics and Applications, Microcosm Press.
    [3] Curtis, H.W. (2009 2nd ed.), yes later versions exist,
        Orbital Mechanics for Engineering Students. Elsevier Ltd.
    [4] Vallado, David A., (2022, 5th ed., 1st printing),
        Fundamentals of Astrodynamics and Applications, Microcosm Press.
            Focus on chap.12 interplanetary mission analysis, but support from
            the whole book; be careful there are errors
    [5] Wertz, James R., (2001, 5th ed.),
        Mission Geometry; Orbit and Constellation Design and Management,
        Space Technology Library, Microcosm Press et.al.
        Julian dates, chap.4, section time, p.185,
            review https://en.wikipedia.org/wiki/Julian_day
    [6] Meeus, Jean, (1998 2nd. ed.) Astronomical Algorithms,
        Willmann-Bell.  Julian dates, pp.61-62 in particular
    [7] Duffett-Smith et.al. (2017 4th ed.),
        Practical Astronomy with your Calculator or Spreadsheet,
            section 4, Julian dates (CDJD), pp.8-31
    [8] Brown, C.D., (1992 2nd printing), AIAA Education Series,
        Spacecraft Mission Design,
            section 6, Interplanetary Trajectories, pp.95-131
Notes:
----------
    Julian date conversions (BCE, pre-Gregorian/Gregorian) is not as
        straight forward as I expected.
        [6] Meeus outlines reasonable conversion boundry tests, but it took
            [7] Duffett-Smith to make the calculation implementation clear.
            I understand the conversion concepts, but the algorythm explanation
            is still not clear to me.  Note some modern conversions may not pass
            the boundary tests.

"My" fixed ecliptic frame graphic. The left graphic is intended to show the planets in nearly the same plane. The right graphic is intended to show that to move between planets still requires important inclination changes.

Code for the graphics.

# magic command has format problem with vscode black formatter,
#   but the interactive plot works ok; though is it commented out for now
# %matplotlib widget
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D # helps with subtitles :--)
from skyfield.api import Loader
from skyfield.framelib import ecliptic_frame  # conversion to ecliptic

# Load planets ephemeris
# prevent multiple copies of large data files, get ephemerous data from common folder
#   https://rhodesmill.org/skyfield/api-iokit.html
ephem_dir = r"~\Documents\!Data\OrbMech\!Orb_data\spice_data"
load = Loader(ephem_dir)
planets = load("de421.bsp")
# JPL ephemeris files: https://rhodesmill.org/skyfield/planets.html
#   DE421 valid dates 1899-2053, small file, 17MB, circa 2008
#   DE422 valid dates -3000 to 3000, medium file, 623MB, circa 2008
#   DE441 valid dates -13200 to 17191, huge file 3.1GB, circa 2020
#   DE440s valid dates 1849 to 2150, smallish file 32MB, circa 2020
venus = planets["venus"]
earth = planets["earth"]
mars = planets["mars"]
jupiter = planets["JUPITER BARYCENTER"]
# Generate time steps
ts = load.timescale()
times = ts.utc(2023, 1, range(1, 3660))  # start 2023, for xx days/years
# Calculate positions across time range
venus_positions = [venus.at(time).frame_xyz(ecliptic_frame).km for time in times]
earth_positions = [earth.at(time).frame_xyz(ecliptic_frame).km for time in times]
mars_positions = [mars.at(time).frame_xyz(ecliptic_frame).km for time in times]
jupiter_positions = [jupiter.at(time).frame_xyz(ecliptic_frame).km for time in times]

def build_3D_plot(ax):
    # Plot planet position, at t[0]
    ax.scatter(
        venus_positions[0][0],
        venus_positions[0][1],
        venus_positions[0][2],
        c="r",
        marker="o",
    )
    ax.scatter(
        earth_positions[0][0],
        earth_positions[0][1],
        earth_positions[0][2],
        c="r",
        marker="o",
    )
    ax.scatter(
        mars_positions[0][0],
        mars_positions[0][1],
        mars_positions[0][2],
        c="r",
        marker="o",
    )
    ax.scatter(
        jupiter_positions[0][0],
        jupiter_positions[0][1],
        jupiter_positions[0][2],
        c="r",
        marker="o",
    )
    # Plot orbits
    ax.plot(
        [pos[0] for pos in venus_positions],
        [pos[1] for pos in venus_positions],
        [pos[2] for pos in venus_positions],
        label="Venus",
    )
    ax.plot(
        [pos[0] for pos in earth_positions],
        [pos[1] for pos in earth_positions],
        [pos[2] for pos in earth_positions],
        label="Earth",
    )
    ax.plot(
        [pos[0] for pos in mars_positions],
        [pos[1] for pos in mars_positions],
        [pos[2] for pos in mars_positions],
        label="Mars",
    )
    ax.plot(
        [pos[0] for pos in jupiter_positions],
        [pos[1] for pos in jupiter_positions],
        [pos[2] for pos in jupiter_positions],
        label="Jupiter",
    )
    # Set labels and title
    ax.set_xlabel("X-axis [km]")
    ax.set_ylabel("Y-axis [km]")
    ax.set_zlabel("Z-axis [km]")

    ax.legend()  # puts in plot label's
    return


fig = plt.figure(figsize=(10, 6), constrained_layout=True)
# ax1 = 1row, 2columns, 1st column
ax1 = fig.add_subplot(121, projection="3d")
build_3D_plot(ax1)
ax1.set_aspect("equal")
ax1.view_init(elev=30, azim=-45, roll=0)
# set ax1 figure title and subtitle
ax1.set_title("XYZ-axis equal ranges", fontsize=12)
ax1.text2D(0.5, 0.98, "(ecliptic frame; flat orbital inclination)", transform=ax1.transAxes,
            horizontalalignment='center', fontsize=10, color='black')

# ax2 = 1row, 2columns, 2nd column
ax2 = fig.add_subplot(122, projection="3d")
build_3D_plot(ax2)
ax2.set_title("Z-axis autoscale: exaggerated inclination", fontsize=12)
# ax2.text2D(0.5, 0.98, "(ecliptic frame; note orbital plane tilt)", transform=ax1.transAxes,
#             horizontalalignment='center', fontsize=10, color='black')
ax2.set_aspect("auto")
ax2.view_init(elev=30, azim=-45, roll=0)
fig.suptitle("Planetary Orbits: inclination comparison: 2023-2033 (ecliptic)", fontsize=14)
plt.show()