a .astronomy hack from sydney!

spaceprobes.py

from astropy.utils.data import download_file
import matplotlib.pyplot as plt
import json
import numpy as np
import astropy.units as u
from astropy.time import Time
from astropy.coordinates import (SkyCoord, CartesianRepresentation,
                                 SphericalRepresentation, AltAz, Angle,
                                 EarthLocation, CylindricalRepresentation)
import ephem

# Download the DSN data from spaceprob.es
dsn_url = 'http://murmuring-anchorage-8062.herokuapp.com/dsn/probes.json'
f = download_file(dsn_url, cache=False)
dsn = json.loads(open(f).read())['dsn_by_probe']

# Parse out the data that we care about
spacecraft = []
distances = []
altitudes = []
azimuths = []
stations = []
times = []
for probe in dsn:
    if 'downlegRange' in dsn[probe] and 'azimuthAngle' in dsn[probe]:
        distance = dsn[probe]['downlegRange']

        if distance != '-1.0':
            spacecraft.append(probe)
            distances.append(float(distance))
            stations.append(dsn[probe]['station'])
            altitude = float(dsn[probe]['elevationAngle'])
            if altitude*u.deg > 90*u.deg:
                    altitude = 90
            altitudes.append(altitude)

            azimuths.append(float(dsn[probe]['azimuthAngle']))
            times.append(Time(dsn[probe]['last_conact'][:-1], format='isot'))

# Make the alt/az/dist astropy quantities
distances = u.Quantity(distances, unit=u.km)
altitudes = u.Quantity(altitudes, unit=u.deg)
azimuths = u.Quantity(azimuths, unit=u.deg)

# Define the positions of the DSN dishes
canberra = EarthLocation.from_geodetic(-149.1244*u.deg, -35.3075*u.deg, 0*u.m)
goldstone = EarthLocation.from_geodetic(Angle('-116d53m24s'),
                                        Angle('35d25m36s'), 0*u.m)
madrid = EarthLocation.from_geodetic(Angle('-4d14m57s'),
                                     Angle('40d25m45s'), 0*u.m)
observatory_dict = dict(Canberra=canberra,
                        Goldstone=goldstone,
                        Madrid=madrid)

# Make SkyCoord objects for each spacecraft
coords = []
for i in range(len(spacecraft)):
    c = SkyCoord(alt=altitudes[i], az=azimuths[i], distance=distances[i],
                 frame=AltAz(obstime=times[i],
                             location=observatory_dict[stations[i]]))
    coords.append(c.icrs)

coords = SkyCoord(ra=[coord.ra for coord in coords],
                  dec=[coord.dec for coord in coords],
                  distance=[coord.distance for coord in coords],
                  frame='icrs')
# Convert to ecliptic cartesian coordinates
cartesian_coords = coords.barycentrictrueecliptic.cartesian

# Make the plot!
fig, ax = plt.subplots(1, 1, figsize=(10, 10))
ax.plot(cartesian_coords.x.to(u.AU).value,
         cartesian_coords.y.to(u.AU).value, 'w.')

print(spacecraft, [(i.x.to(u.AU), i.y.to(u.AU)) for i in cartesian_coords])

object_positions = dict()
for name, coord in zip(spacecraft, cartesian_coords.xyz):
    object_positions[name] = {}
    for axis, c in zip(['x', 'y', 'z'], coord):
        object_positions[name][axis] = c.to(u.AU).value

# Plot the orbits of the planets (assuming circular)
semimajor_axes = u.Quantity([0.38, 0.72, 1, 1.5, 5.2, 9.5, 19.1, 30.0],
                            unit=u.AU)
n_points = 100
theta = np.linspace(0, 2*np.pi, n_points)*u.radian
r = np.ones(n_points)[np.newaxis,:]*np.atleast_2d(semimajor_axes).T
z = np.zeros(n_points)*u.m
orbits_cyl_coords = SkyCoord(rho=r, phi=theta, z=z,
                             representation='cylindrical')
orbits_cartesian = orbits_cyl_coords.icrs.cartesian

plt.plot(orbits_cartesian.x.T, orbits_cartesian.y.T)

# Label each spacecraft
for i, name, dist in zip(range(len(spacecraft)), spacecraft, distances):
    ax.text(cartesian_coords.x[i].to(u.AU).value,
             cartesian_coords.y[i].to(u.AU).value, name,
             color='w')

# plot Earth
earth_coord = SkyCoord(ra=0*u.deg, dec=0*u.deg, distance=0*u.m, frame='gcrs',
                       obstime=Time.now())
cartesian_earth = earth_coord.icrs.barycentrictrueecliptic.cartesian
ax.plot(cartesian_earth.x.to(u.AU).value,
         cartesian_earth.y.to(u.AU).value, 'o', markersize=10, color='cyan')

# plot other planets
pyephem_planets = [ephem.Mercury(), ephem.Venus(), ephem.Mars(),
                   ephem.Jupiter(), ephem.Saturn(), ephem.Uranus(),
                   ephem.Neptune()]
observer = ephem.Observer()
observer.date = Time.now().datetime

planet_coords = []
for planet in pyephem_planets:
    planet.compute(observer)
    c = SkyCoord(ra=float(planet.ra)*u.rad,
                 dec=float(planet.dec)*u.rad,
                 distance=planet.earth_distance*u.AU, frame='gcrs',
                 obstime=Time.now())
    planet_cartesian = c.barycentrictrueecliptic.cartesian
    ax.plot(planet_cartesian.x.to(u.AU).value,
            planet_cartesian.y.to(u.AU).value, 'o', markersize=10)
    object_positions[planet.name] = {}
    for axis, coord in zip(['x', 'y', 'z'], planet_cartesian.xyz):
        object_positions[planet.name][axis] = coord.to(u.AU).value

ax.plot(0, 0, 'o', markersize=10, color='yellow')

ax.set_axis_bgcolor('k')
ax.set(xlabel='AU', ylabel='AU')
ax.set_aspect('equal')
plt.show()