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

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.

• So as it turns out this is an XY problem. The question states "... given their look angles?" but the accepted and best answer does not use look angles at all. I'd recommend you remove "look angles" from your title, and instead add "using Skyfield", since your accepted answer only works with Skyfield and does not in any way help if you are starting with look angles. – uhoh Apr 1 '18 at 1:46
• @uhoh You're right. I would remain the question title and accept your answer. However, edited the question and mentioned Brandon's answer – Leeloo Apr 1 '18 at 8:33
• It's your choice. The best guidance is to leave a question and its answers in such a state that a person searching for your question will find the answers helpful and useful, and so all is well; all three answers are here. If you remove "look angles" and add "Skyfield" to the title (which best reflects what you are doing and wanted to know how to do), this question will be more specific. That's what I'd recommend. But it's not a big deal either way. The title could be something like "Cone angle between to satellites with Skyfield; do I need the look angles, or is there a better way?" – uhoh Apr 1 '18 at 9:12

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

earth   = data['earth']

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]

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 = []
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()

• To find the elevation angle relative to ground station, I just used sat_i = EarthSatellite(L1, L2) ground_i = Topos(g_i[0],g_i[1]) diff_i = sat_i - ground_i. And then elevation=diff_i.at(time).altaz()[0].degrees. – Leeloo Mar 31 '18 at 13:28
• In rhodesmill.org/skyfield/earth-satellites.html it's said that observer() is useless and too expensive for Earth satellites. – Leeloo Mar 31 '18 at 13:32
• And why did you use elevation_m = 1300.0 for ground station? – Leeloo Mar 31 '18 at 13:36
• @Leeloo If you have a different, better, simpler, way to do this, then it would be great to post it as an additional answer. There is nothing wrong with answering your own questions from time to time - in fact it is encouraged if you have the solution that solves your problem/question the best. You can even accept it as well. It looks like you've read the documentation more than I have, consider posting your was as an additional answer. – uhoh Mar 31 '18 at 16:52
• "And why did you use elevation_m = 1300.0 for ground station?" Because the elevation of a ground station in White Sands would be about 1300 meters. – uhoh Mar 31 '18 at 16:55

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.
• Thanks! Could you add some test data from an independent source, so I could check? – Leeloo Mar 29 '18 at 7:08
• This certainly looks good to me! – uhoh Mar 29 '18 at 10:19
• ...except for items 2.2 and 2.3. Those need to also be multiplied by cosine of elevation in addition to what is already shown, as shown here. – uhoh Mar 29 '18 at 11:54
• @uhoh Thank you for review, indeed there was a mistake. Now it should be fine. – Adam Trhon Mar 30 '18 at 7:27

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.
t = ts.utc(2018, 3, 30, 23, 8)

# Satellites.
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))

• Great! What about the angle between 2 satellites from the center of the Earth? – Leeloo Mar 31 '18 at 23:19
• What about the angle between s1 and usno from the Earth center? – Leeloo Mar 31 '18 at 23:43
• This is of course a better way to do it using Skyfield. I've answered the OP's question as asked "... given their look angles?" without using spherical trig here. Would this be ok under the (unnecessary) constraints of the question? – uhoh Apr 1 '18 at 1:43
• @Leeloo You would do pos1 = s1.at(t) and pos2 = s2.at(t) for the vectors to the satellites from the Earth's center. The value usno.at(t) is, similarly, the vector from the Earth's center to the USNO. – Brandon Rhodes Apr 1 '18 at 12:14
• @BrandonRhodes Thanks! And pos1=s1.at(t); pos2=(s1-usno).at(t); pos1.separation_from(pos2); gives the angle between ground station-satellite vector and satellite-Earth center vector? – Leeloo Apr 1 '18 at 22:41