Calculate the visual magnitude of a satellite (e.g. the ISS) using PyEphem

visual_magnitude.py

#!/usr/bin/env python
import ephem
import math

'''
  Takes two PyEphem objects, one for the sun and one for the satellite, and the
  satellite's standard magnitude value (-1.3 for the International Space Station),
  and calculates the visual magnitude of the satellite.
'''
def calculate_visual_magnitude(sun, satellite, standard_magnitude):
  # Calculate distance from observer to the sun, and to the satellite (length a and b of our triangle).
  sun_distance = (sun.earth_distance * ephem.meters_per_au) - ephem.earth_radius
  satellite_distance = satellite.range

  # Calculate angle separating sun and satellite (angle C of our triangle)
  separation_angle = ephem.separation(satellite, sun)

  # Calculate distance between sun and satellite (length c of our triangle)
  # c = sqrt(a*a + b*b - 2 * a * b * cos(C))
  sun_satellite_distance = math.sqrt((sun_distance*sun_distance) + (satellite_distance*satellite_distance) - (2 * sun_distance * satellite_distance * math.cos(separation_angle)))

  # Now we know the length of all sides of the triangle, calculate the phase angle (angle A of our triangle)
  # A = acos((b*b + c*c - a*a) / (2 * b * c))
  phase_angle = math.acos(((satellite_distance*satellite_distance) + (sun_satellite_distance*sun_satellite_distance) - (sun_distance*sun_distance)) / (2 * satellite_distance * sun_satellite_distance))

  # Calculate visual magnitude (source: Robert Matson, <http://www.satobs.org/seesat/Apr-2001/0313.html>)
  magnitude = standard_magnitude - 15 + (5 * math.log10(satellite.range / 1000)) - 2.5 * math.log10(math.sin(phase_angle) + ((math.pi-phase_angle)*math.cos(phase_angle)))

  return round(magnitude, 1)

''' Example usage '''
sun = ephem.Sun()
satellite = ephem.readtle("ISS (ZARYA)",
  "1 25544U 98067A   18205.83140740  .00001117  00000-0  24221-4 0  9998",
  "2 25544  51.6402 191.0806 0004072 354.9923  73.3091 15.54007570124377")

obs = ephem.Observer()
obs.lon = '-123.138330'
obs.lat = '49.315697'
obs.elevation = 111

sun.compute(obs)
satellite.compute(obs)

print 'Current satellite magnitude is %s' % (calculate_visual_magnitude(sun, satellite, -1.3),)