build a general position object from: r_vec, v_vec, time epoch, and gm

[beluej123]

Is there a better workflow to build a general position object to get the osculating elements of a transfer orbit.
-- I answered my original questions --
My task is to design a series of 1st order solar system transfer orbits (Keplerian transfer for initial design).

From Spacecraft Mission Design, C.D. Brown, AIAA, 1992, pp.106-113.
Example, Earth to Venus transfer, 1988-04-08->1988-07-26 (109 days).

My workflow process:

1. Skyfield, get Earth ephemeris; 1988-04-08 (t1), need r1_vec.
2. Skyfield, get Venus ephemeris; 1988-07-26 (t2), need r2_vec.
3. Use a Lambert solver to compute transfer velocity vectors (v1_vec and v2_vec), given:
   tof (time-of-flight), r1_vec, r2_vec, .
   Several solvers out there; I generally use my forked version of
   https://github.com/jorgepiloto/lamberthub
4. Convert r1_vec and v1_vec to units expected by skyfield.positionlib ICRF()
5. Skyfield, build position object with ICRF() from r1_vec, v1_vec, t1, and barycenter gm.
   a) calculate osculating elements

Thank you for your help.

The code below works fine at this point. I get the orbit transfer elements, then propagate the transfer orbit using Kepler elements (using true anomaly, aka; nu or ta), then rotate elements to 3D cartesian space for plotting. If anyone is interested I can provide the complete set of code - including plot code.

import numpy as np # np only used for print, in this code
from skyfield.positionlib import ICRF # convert vectors to position object
# Osculating elements, https://rhodesmill.org/skyfield/elements.html
from skyfield.elementslib import osculating_elements_of
from skyfield.api import Loader # to set my folder for ephemeris files
# 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
# sun gravitational parameter; D.A. Vallado [4], 5th ed., 2022
mu_s = 1.32712428e11  # [km^3/s^2], sun, Vallado [4] p.1059, tbl.D-5
# units conversion factor
au_ = 149598023.0  # [km], Vallado [4] p.1059, tbl.D-5

ts = load.timescale()
t1 = ts.utc(1988, 4, 8, 0, 0)  # launch time
t2 = ts.utc(1988, 7, 26, 0, 0)  # arrive time
tof = (t2 - t1) * 24 * 3600  # Time of flight in seconds

# get planets ephemeris
earth, venus = planets["earth"], planets["venus"]
# calculate earth/venus positions at launch and arrival; fool with units as needed
posE_t1 = earth.at(t1)  # earth at t1, launch
posE_launch_km = earth.at(t1).position.km  # earth ICRF position, launch
posV_arrive_km = venus.at(t2).position.km  # venus ICRF position, arrive

# ************* Lambert Solutions ************
# get transfer orbit from Lambert Problem/Solution
from lamberthub import gooding1990


r1_vec = posE_launch_km  # initial position vector
r2_vec = posV_arrive_km  # final position vector

v1_vec, v2_vec = gooding1990(
    mu_s,
    r1_vec,
    r2_vec,
    tof,
    M=0,
    prograde=True,
    low_path=True,  # Type of solution
    maxiter=35,
    atol=1e-5,
    rtol=1e-7,
    full_output=False,  # Iteration config
)
# print departure/arrival positions for sanity checks
np.set_printoptions(precision=4, floatmode='fixed')
print(f"Earth departure position: {t1.utc_strftime()}")
print(f"    r1_vec: {r1_vec} [km]")
print(f"    v1_vec: {v1_vec} [km/s]")
print(f"    v1_mag: {np.linalg.norm(v1_vec):.6g} [km/s]")
print(f"Venus arrival position: {t2.utc_strftime()}")
print(f"    r2_vec: {r2_vec} [km]")
print(f"    v2_vec: {v2_vec} [km/s]")
print(f"    v2_mag : {np.linalg.norm(v2_vec):.6g} [km/s]")

# convert r1_vec and v1_vec to units expected by ICRF()
r1_vec_au=r1_vec/au_
v1_vec_au_day=(v1_vec/au_)*(60*60*24)

pos1=ICRF(position_au=r1_vec_au, velocity_au_per_d=v1_vec_au_day, t=t1, center=None)
elements = osculating_elements_of(position=pos1, reference_frame=None, gm_km3_s2=mu_s)
# list attributes/values in elements
print(f"\nTransfer orbital elements:")
print(f"   Calculated periapsis time: {elements.periapsis_time.utc_strftime()}")
for attr in dir(elements):
    if not attr.startswith("_"):
        print(f"   {attr}, {getattr(elements, attr)}")

The following graphic helps provide context. Easy to visualize the departure/arrival ellipse (ecc ~= 0.17, sma ~= 1.3e8 [km], calculated manually), but the point is to have available the "full set" of skyfield tools for a host of other calculations. Thank you.

[brandon-rhodes]

Is there a better workflow to build a general position object to get the osculating elements of a transfer orbit.

I don't see any immediate improvements to what you're doing, the code is easy to read and I understand both the numbers that are pulled from Skyfield for use by the other library, and then the numbers that get plugged back in. This will be a fun example for other folks to look at who might be interested in the same kind of problem solving. One tweak I might make is to say center=0 on the position you build so that Skyfield knows it's centered at the SSB — or, even better, you could use Barycentric instead of ICRF and then the SSB is the default center without you needing to specify it.