I’m attempting to use astropy to calculate the solar & lunar illumination over a sequence of regularly spaced intervals. My code seems to work, but it seems pretty slow - I’m calculating both values every 15 minutes over 30 days and the lunar calculations take around 10 seconds, while the solar calculations take around 24s.
Is there anything I can do to speed things up? I should note, I’m just a software guy - I have no experience in astronomy and/or astropy so feel free to explain the astronomy side to me like I’m a child! Full code below:
import numpy as np
from astropy.coordinates import get_body, get_sun, EarthLocation, AltAz
from astropy.time import Time, TimeDelta
from astropy.coordinates import solar_system_ephemeris
import astropy.units as u
from datetime import datetime, timezone
import time as tm
def calculate_solar_millilux(time, location) -> float:
"""
Calculates solar illuminance in millilux.
"""
# 1. Define Location and Time
loc = location
# 2. Get Sun position (Alt/Az)
sun_altaz = get_sun(time).transform_to(AltAz(obstime=time, location=loc))
elevation = sun_altaz.alt.deg
# 3. Model: Illuminance based on elevation
# This is a simplified model (e.g., Kasten & Young or empirical approximation)
if elevation <= 0:
return 0.0, elevation # Sun is below horizon
# Approximate solar illuminance (lux) as a function of elevation (degrees)
# At zenith (90 deg), direct sun is ~100,000 - 130,000 lux
# Formula approximation for bright day:
lux = 130000 * np.sin(np.radians(elevation)) * (1 - 0.1 * np.exp(-elevation/10))
# Convert to millilux
millilux = lux * 1000
return millilux, elevation
def calculate_lunar_illumination(time, location):
# Get moon position
with solar_system_ephemeris.set('builtin'):
moon = get_body('moon', time, location)
altaz = AltAz(obstime=time, location=location)
moon_altaz = moon.transform_to(altaz)
# Calculate phase (0.0 to 1.0)
# 1.0 is full moon, 0.0 is new moon
sun = get_sun(time)
elongation = moon.separation(sun, origin_mismatch="ignore")
phase = (1 - np.cos(elongation)) / 2.0
# Illumination calculation (approximate formula)
# Full moon zenith illumination ≈ 300-400 millilux (mlx)
max_lux = 350 # millilux at zenith
altitude_deg = moon_altaz.alt.deg
# If moon is below horizon, illumination is 0
if altitude_deg < 0:
millilux = 0.0
else:
# Basic correction for altitude (sin(altitude))
# Note: Atmospheric extinction not fully modeled here
millilux = max_lux * phase * np.sin(np.radians(altitude_deg))
return millilux, altitude_deg, phase
if __name__ == "__main__":
# Set location (observer)
location = EarthLocation(lat=42.28, lon=-83.74, height=250) # Example: Ann Arbor, MI
today = datetime.combine(datetime.today(), datetime.min.time(), tzinfo=timezone.utc)
time = Time(today)
interval_minutes = 15
numDays = 30
times = [time + TimeDelta(interval_minutes * i * u.minute) for i in range(0, numDays*24*60//interval_minutes)]
print(f"Calculating illuminance for {location} starting at {today.isoformat()} (UTC)")
print(f" {len(times)} points ({numDays} days at {interval_minutes} minute intervals)")
print(f" --------------------------------------------------------------------------")
solar_illuminations = []
lunar_illuminations = []
start_time = tm.perf_counter()
for t in times:
sol_illum, sol_elev = calculate_solar_millilux(t, location)
solar_illuminations.append((t, sol_illum))
end_time = tm.perf_counter()
elapsed_time = end_time - start_time
print(f" Calculated solar illuminance for {len(times)} time points in {elapsed_time:.6f} seconds.")
start_time = tm.perf_counter()
for t in times:
lunar_illum, lunar_alt, lunar_phase = calculate_lunar_illumination(t, location)
lunar_illuminations.append((t, lunar_illum))
end_time = tm.perf_counter()
elapsed_time = end_time - start_time
print(f" Calculated illuminance for {len(times)} time points in {elapsed_time:.6f} seconds.")