Z value of astrometric.velocity.m_per_s

[Captjackz]

To Whom, I have written a python simulation of Starlink launch, insertion, final orbit adjust Crew_Dragon launch and ISS rendezvous Artemis launch-moon orbit Perseverance Earth-Mars launch and mars rendezvous THIS FAILS

In perfecting the Earth-Mars system I needed a needed a precise location of Sun-Earth-Mars on the Perseverance launch date and other dates. I turned to SKYFIELD.

COORDINATE SYSTEM: Sun Centered - Earth Ecliptic plane. In this plane, the earth should be rotating with an sinusoidal X-Y position and velocity. There should be no or little Z-value. I get huge z values for earth, mars, venus and mercury. So I wrote a small test program. (below) That shows a particular bad value for Vz of 11581.3, to point you to the problem. Switching from EMB to earth does not change the result. In my simplified model (for this request), the earth should have NO (little) Z displacement or velocity. I have tried dozens of options and tried to buy numerous books understand my coordinate system problem. If I set the Z values for mercury, Venus, earth, mars to zero my simulation comes close, but now a viable solution. Can you please HELP. If you could recommend a text to clarify my confusion of the reference frame I would appreciate it.

Jack Finley

from skyfield.api import load
from skyfield.framelib import ecliptic_frame
from skyfield import almanac
from skyfield.api import N, W, wgs84, load
import math
planets = load('de421.bsp')
sun, earth,EMB = planets['sun'],planets['EARTH BARYCENTER'],planets['EARTH BARYCENTER']

AU  = 149597870700
ts = load.timescale()
t = ts.utc(2021, 6, 30, 15, 19)
t = ts.now()        #a test vector to see if repeated requests changed results over time.

astrometric = sun.at(t).observe(EMB)
ra, dec, distance = astrometric.radec()

position = sun.at(t).observe(earth)
ra, dec, distance = position.radec()
d = sun.at(t).observe(EMB).apparent().distance()
earthau = d.au
earthang = 360*(ra.hours)/24
print("Earth          %6.2f         %s         %s"%( earthang, dec, distance))

earthpx, earthpy, earthpz = astrometric.position.m
earthvx, earthvy, earthvz = astrometric.velocity.m_per_s

earth.dst = math.sqrt(earthpx**2 + earthpy**2 + earthpz**2)/AU
earth.vel = math.sqrt(earthvx**2 + earthvy**2 + earthvz**2)

print('Earth   Pos %16.2f  %16.2f  %16.2f  %8.2f'%(earthpx/AU, earthpy/AU, earthpz/AU, earth.dst))                 
print('Earth   Vel %16.1f  %16.1f  %16.1f  %8.2f'%(earthvx, earthvy, earthvz, earth.vel))

#barycentric = sun.at(t)
#astrometric = barycentric.observe(EMB)
#astrometric = sun.at(t).observe(EMB)
#apparent = astrometric.apparent()

RESULT:

Earth          349.78         -04deg 23' 52.0"         1.00659 au
Earth   Pos             0.99             -0.18             -0.08      1.01
Earth   Vel           5258.7           26716.1           11581.3  29589.38

[Captjackz]

sky2.zip

[Captjackz]

Still learning this system.

[brandon-rhodes]

(I've edited your question to add triple-backticks to format your Python as code, but otherwise I think I left your question unaltered.)

The x,y,z coordinates of a Skyfield position are J2000 = ICRS coordinates, as explained at:

https://rhodesmill.org/skyfield/positions.html#the-icrs-reference-system-and-j2000

If you want x,y,z coordinates in another reference frame, you have to do a conversion. Here's code I just tested that I think does what you want:

from skyfield.api import load
from skyfield.framelib import ecliptic_frame
import math

planets = load('de421.bsp')
sun, earth, EMB = (
    planets['sun'],
    planets['earth'],
    planets['earth barycenter'],
)

ts = load.timescale()
t = ts.utc(2021, 6, 30, 15, 19)
t = ts.now()

astrometric = sun.at(t).observe(EMB)

r, v = astrometric.frame_xyz_and_velocity(ecliptic_frame)
earthvx, earthvy, earthvz = v.m_per_s
earth.vel = math.sqrt(earthvx**2 + earthvy**2 + earthvz**2)

print('Earth   Vel %16.1f  %16.1f  %16.1f  %8.2f' % (
    earthvx, earthvy, earthvz, earth.vel
))

The z-axis velocity should print as zero.

How could the documentation be clearer to help future users understand how to produce ecliptic x,y,z coordinates? Could you point to where you looked in the documentation so that I can add it there? Thanks!

[brandon-rhodes]

I'm going to close this in the optimistic hope that my answer solved your problem, but feel free to comment further if you have a follow-up question or want to report that the answer worked for you.