How to calculate cone angle between two satellites given their look angles?

Question

I have the azimuth, elevation and distance of 2 satellites relative to a ground station, defined by latitude and longitude, as I've discussed in my previous question. I'm using the satellites' TLEs and the Python package Skyfield's .altaz() method to obtain their alt/az/el.

How can I calculate the cone angle between 2 satellites relative to the mentioned ground station?

UPDATE

As Brandon Rhodes answered here, it's not necessary to use look angles. Skyfield is able to calculate the separation angle from the positions.

Answer 1

note: This answer addresses the question directly:

How to calculate cone angle between two satellites given their look angles?

If you need to use the look angles, this is a good way to do it. This better answer explains to the OP that if you are using Skyfield, that you should not use the look angles but instead use the coordinates in their original form.

The cosine of the angle between two vectors is given by the dot product of their norms.

$$\cos(\theta_{12}) = \mathbf{\hat{v}_1} \cdot \mathbf{\hat{v}_2} = \frac{\mathbf{v_1} \cdot \mathbf{v_2}}{|\mathbf{v_1}| \ |\mathbf{v_2}|} = \frac{\mathbf{v_1} \cdot \mathbf{v_2}}{v_1 \ v_2}$$

$$\theta_{12} = \cos^{-1}\left( \frac{\mathbf{v_1} \cdot \mathbf{v_2}}{v_1 \ v_2} \right)$$

But if you drop the distance and just use $az, el$ you get a normalized vectors automatically:

$$\mathbf{\hat{v}_i} = \cos(el_i)\left(\cos(az_i)\mathbf{\hat{x}}+\sin(az_i)\mathbf{\hat{y}}\right) + \sin(el_i)\mathbf{\hat{z}}$$

I think this is no different than what @AdamTrhon has already described in this answer! You can use spherical trig and possibly the law of cosines, but sometimes those expressions lead to computational errors due to such things as subtraction of nearly-equal numbers and divides by nearly zero, whereas this way - working in cartesian as much as possible, in the way shown here at least, - there are no chances of that happening.

In Python that would be something like:

def nvec(elaz):
    (cel, sel), (caz, saz) = [[f(q) for f in (np.cos, np.sin)] for q in elaz]   # parentheses for Py3
    return np.array([cel*caz, cel*saz, sel])

def angle(elaz1, elaz2):
    v1, v2 = [nvec(elaz) for elaz in (elaz1, elaz2)]  # parentheses for Py3
    return np.arccos((v1*v2).sum(axis=0))

So if you download TLEs for three TDRS satellites plus the ISS and calculate the cone angle between TDRS pairs, and between the ISS and each TDRS, you'd get something like what's shown below.

The objects also have their geocentric positions stored, so you can make a 3D map like is shown in this example.

---

EDIT: I have used WhiteSands.at(times).observe(sat.ICRF).apparent().altaz() and this would be the recommended way to get the apparent optical position. The .apparent() method includes a variety of effects, including the light-time delay, astronomical aberration, and even... wait for it... gravitational effects of massive bodies that might deflect the path, as well as atmospheric refraction. You can read more about this in the documentation at http://rhodesmill.org/skyfield/api-position.html#skyfield.positionlib.Astrometric.apparent

---

TLEs = """TDRS 5                  
1 21639U 91054B   18086.36437858  .00000071  00000-0  00000-0 0  9995
2 21639  14.5306  18.4626 0026343 345.4651 144.6288  1.00281508 97593
TDRS 10                 
1 27566U 02055A   18087.12756861  .00000056  00000-0  00000+0 0  9998
2 27566   5.5204  57.1630 0011308 250.2541 109.7547  1.00278469 56111
TDRS 11                 
1 39070U 13004A   18086.87347718  .00000063  00000-0  00000-0 0  9994
2 39070   5.0128 328.6219 0008993 321.0893  38.7221  1.00272889 16583
ISS (ZARYA)             
1 25544U 98067A   18088.22902370  .00003740  00000-0  63642-4 0  9999
2 25544  51.6415  57.3234 0001506 271.5382 195.6957 15.54152785106088"""

lines           = TLEs.splitlines()
names, L1s, L2s = [[x.strip() for x in lines[i::3]] for i in range(3)]
triplets        = zip(names, L1s, L2s)

class Sat(object):
    def __init__(self, name):
        self.name = name

def nvec(elaz):
    (cel, sel), (caz, saz) = [[f(q) for f in (np.cos, np.sin)] for q in elaz]   # parentheses for Py3
    return np.array([cel*caz, cel*saz, sel])

def angle(elaz1, elaz2):
    v1, v2 = [nvec(elaz) for elaz in (elaz1, elaz2)]  # parentheses for Py3
    return np.arccos((v1*v2).sum(axis=0))

import numpy as np
import matplotlib.pyplot as plt
from skyfield.api import Loader, Topos, EarthSatellite
import itertools

halfpi, pi, twopi = [f*np.pi for f in (0.5, 1, 2)]  # parentheses for Py3
degs, rads = 180/pi, pi/180

load = Loader('~/Documents/fishing/SkyData')  # avoids multiple copies of large files
ts   = load.timescale()

data    = load('de421.bsp')
earth   = data['earth']
ts      = load.timescale()

WhiteSands  = earth + Topos(latitude_degrees   =   32.4,
                             longitude_degrees = -106.5,
                             elevation_m       = 1300.0  )

minutes = np.arange(0, 1441, 1)
times   = ts.utc(2018, 3, 29, 0, minutes)

sats = []
for name, L1, L2 in triplets:
    sat = Sat(name)
    sats.append(sat)
    sat.Geo  = EarthSatellite(L1, L2)
    sat.ICRF = earth + EarthSatellite(L1, L2)
    sat.obs  = WhiteSands.at(times).observe(sat.ICRF)
    sat.elaz = [x.radians for x in sat.obs.apparent().altaz()[:2]]
    sat.below = sat.elaz[0] <= 0.
    sat.pos   = sat.Geo.at(times).position.km

ISS   = [sat for sat in sats if 'ISS' in sat.name][0]
TDRSs = [sat for sat in sats if 'TDRS' in sat.name]
TDRSpairs = list(itertools.combinations(TDRSs, 2))

intra_TDRS_cones = []
for pair in TDRSpairs:
    name = ''.join([x.name + ' ' for x in pair])[:-1]
    elaz1, elaz2 = [s.elaz for s in pair]
    cone = angle(elaz1, elaz2)
    intra_TDRS_cones.append((name, cone))

ISS_TDRS_cones = []
for TDRS in TDRSs:
    name = ISS.name + ' ' + TDRS.name
    cone = angle(ISS.elaz, TDRS.elaz)
    ISS_TDRS_cones.append((name, cone))

if True:
    plt.figure()
    plt.subplot(2, 1, 1)
    for name, cone in intra_TDRS_cones:
        plt.plot(minutes/60., degs*cone)
    plt.title("intra-TDRS cone angles", fontsize=16)
    plt.xlim(0, 24)
    plt.subplot(2, 1, 2)
    for name, cone in ISS_TDRS_cones:
        plt.plot(minutes/60., degs*cone)
    plt.title("ISS-TDRS cone angles", fontsize=16)
    plt.xlim(0, 24)
    plt.show()

Answer 2

Here is manual approach:

1. Setup orthogonal coordinate system:
   1. Unit is 1Km (but it does not matter much).
   2. Origin is in ground station.
   3. x-axis points to 0° azimuth, y-axis to 90°, z-axis straight up.
   4. Fact: x-y-plane is tangent plane of Earth.
2. Compute unit vectors pointing of each satellite. Beware of radians/degrees.
   1. Unit vector's x-coordinate is cosine of azimuth.
   2. Unit vector's y-coordinate is sine of azimuth.
   3. Fact: The vector now points to the satellite, but only in the x-y-plane.
   4. Vector's z-coordinate is tangent of elevation.
   5. Fact: Now the vector points to the satellite in 3D space, but it is no longer a unit vector.
   6. Normalize it: compute it's length and divide each coordinate by this length.
3. Compute the cone angle:
   1. Compute unit vectors' dot product
   2. Arccos of the dot product is cone angle.

Answer 3

Computing the angle between two vectors is going to be difficult if you first transform their x, y, and z coordinates into angles, because you'll then need to dive in the formulae of spherical trigonometry. Skyfield natively considers all positions to be x,y,z vectors, and it's often easier to compute if you leave them as "position" objects until you're ready to display results.

If you compute two satellite positions relative to an Earth observer, you can then get the angle between the satellites using Skyfield's .separation_from() method that position objects carry:

from skyfield import api

# Time.
ts = api.load.timescale()
t = ts.utc(2018, 3, 30, 23, 8)

# Satellites.
sats = api.load.tle('https://celestrak.com/NORAD/elements/stations.txt')
s1 = sats['ISS']
s2 = sats['ASTERIA']

# Observe the satellites from a position on the Earth's surface.
usno = api.Topos('38.9215 N', '77.0669 W', elevation_m=92.0)
pos1 = (s1 - usno).at(t)
pos2 = (s2 - usno).at(t)

# How far apart are the satellites in the sky?
print(pos1.separation_from(pos2))