Calculate brightness magnitude of ISS (and whether Earth's shadow is falling on it) using pyephem

iss-brightness.py

# by cmdr2, anyone is free to use and modify this, without the need to give credit

import math
import ephem
import time
from datetime import datetime, timedelta

# consts
degrees_per_rad = 180.0 / math.pi
au = 149597892

home = ephem.Observer()
home.lat = '40.7128'
home.lon = '-74.0060'
home.elevation = 0

iss = ephem.readtle('ISS (ZARYA)', '1 25544U 98067A   20036.19908141  .00016717  00000-0  10270-3 0  9067', '2 25544  51.6420 286.4848 0005040 230.6842 129.3862 15.49145399 11399')

sun = ephem.Sun()

iss_std_mag = -1.3

sat = iss
std_mag = iss_std_mag

# all times in UTC
# compare with: https://heavens-above.com/passdetails.aspx?lat=40.7128&lng=-74.0060&loc=Unspecified&alt=0&tz=UCT&satid=25544&mjd=58883.9610767747&type=V
t = datetime(2020, 2, 4, 23, 3, 56)
home.date = t #datetime.utcnow()

sat.compute(home)
sun.compute(home)

actual_sat_range = (sat.range - ephem.earth_radius) / 1000

a = sun.earth_distance * au # - ephem.earth_radius/1000 - 695510 # Km
b = sat.range / 1000
angle_c = ephem.separation( (sat.az, sat.alt), (sun.az, sun.alt) )

c = math.sqrt(math.pow(a, 2) + math.pow(b, 2) - 2*a*b*math.cos(angle_c))

phase_angle = math.acos((math.pow(b, 2) + math.pow(c, 2) - math.pow(a, 2)) / (2 * b * c))

term1 = std_mag
term2 = 5.0 * math.log10(b / 1000)

arg = math.sin(phase_angle) + (math.pi - phase_angle) * math.cos(phase_angle)
term3 = -2.5 * math.log10(arg)

apparent_mag = term1 + term2 + term3

print('magnitude: ' + str(apparent_mag) + ', eclipsed: ' + str(sat.eclipsed))