Local observer coordinate system at satellite panel (LVLH coordinate system)

Question

In skyfield I want to calculate alt-az distance from one reference satellite to next (for antenna pointing simulation). My current workaround is to create the reference observer at sat1 height and calculate the alt-az angles:

sat_observer = sat_ref.subpoint()
ref_PoV = Topos(sat_observer.latitude, sat_observer.longitude, elevation_m=sat_observer.elevation.m)
satellite = Sat(name)
orbit = (satellite - ref_PoV).at(time[0])
el, az, distance = orbit.altaz()

Code is able to calculate relative pointing angles at one time instance. This method is not working for absolute values as the alt-az reference always points to Earth north (the local observer coordinate system at satellite panel will rotate around alt axis in the next time instant).

What would be the most efficient way to define local coordinate system (LVLH frame) at satellite in skyfield to obtain the alt-az angles respective to this LVLH frame?

This is an interesting question! Does "It works..." mean that it didn't fail, or that it returns exactly the values you are expecting? If you like you can also consider going to Skyfield's github page and post an issue there, and include a link to your question here. There is some overlap with users here and in Astronomy SE and in Skyfield's github. (see this question for links) — uhoh, Nov 27 '20 at 00:38
Also, there is more than one possible local observer coordinate system. Can you explain in more detail what you are looking for? See for example Are LVLH and RSW coordinate systems the same thing? — uhoh, Nov 27 '20 at 00:41
@uhoh Thanks for the reference to different coordinate systems, I updated my question and will link it to the skyfield github. — matbru, Nov 27 '20 at 15:11

score 4 · Answer 1 · answered Dec 02 '20 at 15:04

So after help from @uhoh, digging into this post and the discussion here, I managed to produce this minimal working example. Comments appreciated.

from skyfield.api import Loader, EarthSatellite
from skyfield.api import Topos, load
from skyfield.timelib import Time
import skyfield.functions as sf
from sklearn import preprocessing
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
%matplotlib inline
halfpi, pi, twopi = [f*np.pi for f in [0.5, 1, 2]]
degs, rads = 180/pi, pi/180
ts = load.timescale()
line1 = '1 25544U 98067A   14020.93268519  .00009878  00000-0  18200-3 0  5082'
line2 = '2 25544  51.6498 109.4756 0003572  55.9686 274.8005 15.49815350868473'
satellite = EarthSatellite(line1, line2, 'ISS (ZARYA)', ts)
print(satellite)
line1 = '1 43205U 18017A   18038.05572532 +.00020608 -51169-6 +11058-3 0  9993'
line2 = '2 43205 029.0165 287.1006 3403068 180.4827 179.1544 08.75117793000017'
satellite2 = EarthSatellite(line1, line2, 'Roadster', ts)
print(satellite2)
time = ts.utc(2020, 24, 11, np.arange(0, 1, 0.01))
#calculate LVLH reference frame for the reference sat
#Z = - R / ||R||
#Y = Z X V / ||Z X V||
#X = Y X Z
R = satellite.at(time).position.km.T
V = satellite.at(time).velocity.km_per_s.T
Z = -preprocessing.normalize(R, norm='l2')
Y = preprocessing.normalize(np.cross(Z, V), norm='l2')
X = np.cross(Y, Z)
Rpos = satellite2.at(time).position.km.T
#check: LVLH coordinate frame at n events
fig = plt.figure(figsize=[10, 8])  # [12, 10]
ax  = fig.add_subplot(1, 1, 1, projection='3d')
axis_length=20
for i in range(0,5):
    x, y, z = R[i,:]
    u, v, w = X[i,:]
    ax.quiver(x, y, z, u, v, w, length=axis_length, color='red')
    u, v, w = Y[i,:]
    ax.quiver(x, y, z, u, v, w, length=axis_length, color='blue')
    u, v, w = Z[i,:]
    ax.quiver(x, y, z, u, v, w, length=axis_length, color='green')

#construct the rotation matrix at time 0
RM = np.array([X[0,:],Y[0,:],Z[0,:]]).T
#view vector PoV = R_sat - R_ref
PoV = Rpos[:] - R
#rotate PoV to LVLH coordinate system
PoV_LHLV = RM.dot(PoV[0,:])
#go to spherical CS
r1, el1, az1 = sf.to_spherical(PoV_LHLV)
Plot the view angles in polar plot.
plt.figure()
ax = plt.subplot(111, projection='polar')
ax.set_rlim([-90, 90])
ax.set_theta_zero_location('N')
ax.set_theta_direction(1)
ax.set_title('Visibility of satellite2 form satellite PoV', y=1.1)
ax.plot(az1, el1*degs, 'r+')

Local observer coordinate system at satellite panel (LVLH coordinate system)

1 Answers1

Plot the view angles in polar plot.