Skip to content

Instantly share code, notes, and snippets.

@jkseppan
Created May 17, 2016 05:31
  • Star 0 You must be signed in to star a gist
  • Fork 1 You must be signed in to fork a gist
Star You must be signed in to star a gist
Save jkseppan/c12a3877d5e5ca6c358e2a1bf5a01afd to your computer and use it in GitHub Desktop.
from collections import defaultdict, namedtuple
import numpy as np
def visible_cap_central_angle(altitude):
earth_radius = 6371.0
return np.arctan2(np.sqrt(altitude * (altitude + 2 * earth_radius)),
earth_radius)
def hav(x):
return np.sin(x/2)**2
def inv_hav(x):
return np.arcsin(np.sqrt(x))*2
def angle_between(location0, location1):
lat0, lon0 = location0
lat1, lon1 = location1
x = hav(lat1 - lat0) + np.cos(lat0) * np.cos(lat1) * hav(lon1 - lon0)
return inv_hav(x)
Satellite = namedtuple('Satellite', ['lat', 'lon', 'name', 'central_angle'])
def sat(name, lat, lon, alt):
return Satellite(name=name, lat=lat, lon=lon,
central_angle=visible_cap_central_angle(alt))
def visible(x, y):
angle = angle_between((x.lat, x.lon), (y.lat, y.lon))
return angle < x.central_angle + y.central_angle
def read_data(filename):
nodes = set()
with open(filename) as data:
for line in data:
if line.startswith('#'):
continue
fields = line.split(',')
if fields[0].startswith('SAT'):
lat, lon, alt = map(float, fields[1:])
nodes.add(sat(fields[0], np.radians(lat), np.radians(lon), alt))
elif fields[0] == 'ROUTE':
lat0, lon0, lat1, lon1 = map(np.radians, map(float, fields[1:]))
nodes.add(sat('source', lat0, lon0, 0))
nodes.add(sat('destination', lat1, lon1, 0))
else:
assert False, line
adjacent = defaultdict(list)
for x in nodes:
for y in nodes - set(x):
if visible(x, y):
adjacent[x.name].append(y.name)
return adjacent
def dijkstra(adjacent, source, destination):
unvisited = set(adjacent.keys())
distance = { k: np.Infinity for k in adjacent.keys() }
path = { k: None for k in adjacent.keys() }
distance[source] = 0
path[source] = []
while unvisited:
min_dist = min(distance[n] for n in unvisited)
node = [ n for n in unvisited if distance[n] == min_dist ][0]
unvisited.remove(node)
for neighbor in adjacent[node]:
new_dist = distance[node] + 1
if new_dist < distance[neighbor]:
distance[neighbor] = new_dist
path[neighbor] = path[node] + [node]
return distance[destination], path[destination]
if __name__ == '__main__':
adjacent = read_data('data.csv')
dist, path = dijkstra(adjacent, 'source', 'destination')
print(dist, ','.join(path[1:]))
@jkseppan
Copy link
Author

How does this work?

First let's solve which points on the Earth are visible from a satellite. Imagine projecting a light from the satellite towards the Earth: you get a spherical cap delimited by tangents drawn from the satellite to the Earth. Choose some plane through the center of the Earth and the satellite, and draw a right triangle whose vertices are the center of the Earth, the satellite, and either point of tangency in this plane. Since we know the radius of the Earth and the altitude of the satellite, we can solve this triangle, and we're going to be interested in the central angle, which is computed by thevisible_cap_central_angle function. (I chose to use arctan because arccos has numerical difficulties with small angles.)

To decide if a point on the surface is within the spherical cap, look at both the satellite and the point from the center of the Earth, and compute the angle between these directions. This can be done with the haversine formula (which again is numerically better than the spherical law of cosines for small angles, but there are more complicated ways to do this if you want to be precise in a larger range of angles).

The next question is when two satellites see each other. Shine a light from both satellites to get two spherical caps on the Earth. These caps tell you where the Earth blocks the light and casts a shadow. If the caps intersect, the Earth is not between the satellites and they can see each other. If they don't intersect, the Earth blocks the line of sight between the satellites. Thus this question reduces to the previous one.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment