Created
July 26, 2018 00:15
-
-
Save abirkill/8cabfdf7474871d83d96038e0ce90e7d to your computer and use it in GitHub Desktop.
Calculate the visual magnitude of a satellite (e.g. the ISS) using PyEphem
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
#!/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),) |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment