Python code supplementing astronomy.stackexchange answer

asteroid_luminosity_calcs.py

# Python code supplementing an answer to a question on the Astronomy Stackexchange site

# https://astronomy.stackexchange.com/questions/44453/when-thermal-infrared-space-telescopes-spot-asteroids-are-they-seeing-the-body/44456#44456

# We want to compute the relative luminosities produced by an asteroid orbiting the Sun, comparing the monochromatic
# luminosity from reflected sunlight to that from the asteroid's own thermal emission.

# Requires: Numpy, Astropy

import math
import numpy as np
from astropy.constants import h, c, k_B, R_sun
import astropy.units as u


def BB_monochromatic( T, wavelength ):
	"""
	Compute monochromatic power/area emitted by the surface of a blackbody.
	
	Parameters
	----------
	T : float
		temperature in K
	
	wavelength: Quantity
		wavelength in length units
	
	Returns
	-------
	power/area : Quantity
		emitted power [J/s] per Hz per unit surface area
	"""
	
	nu = (c / wavelength.to(u.m)).to(u.hertz)
	T = T * u.K
	power_per_area = math.pi * ((2 * h * nu**3) / c**2) / (np.exp(h * nu / (k_B * T)) - 1.0)
	# reduce the units a bit
	return power_per_area.to( u.W / u.hertz / (u.m*u.m) )


def LumRatio( wavelength, R_ast=100*u.m, dist=1*u.au, T_ast=300.0, emissivity=0.9 ):
	"""
	Returns the monochromatic luminosity of thermal emission, monochromatic 
	luminosity from reflected sunlight, and the ratio of the former to the latter, 
	for the specified wavelength, for an asteroid with diameter R_ast located at dist 
	from the Sun, with blackbody temperature T_as.
	
	Parameters
	----------
	wavelength: Quantity
		wavelength in length units

	R_ast : Quantity, optional
		radius (in length units) for asteroid

	dist : Quantity, optional
		distance (in length units) of asteroid from Sun

	T_ast : Quantity, optional
		temperature in K for asteroid
  
	emissivity : float, optional
		emissivity at specified wavelength for asteroid (= 1 - albedo)
	
	Returns
	-------
	(L_therm, L_refl, ratio) : tuple of (Quantity, Quantity, Quantity)
		L_therm : thermal luminosity of asteroid (W Hz^-1)
		L_refl : reflected-sunlight luminosity of asteroid (W Hz^-1)
		L_therm : ratio of L_therm / L_refl
	"""
	dist_m = dist.to(u.m)
	albedo = 1 - emissivity
	
	L_sun = A_sun * BB_monochromatic(5800, wavelength)
	L_reflected_sun = albedo * math.pi * R_ast**2 * L_sun / (4 * math.pi * dist_m**2)
	
	L_ast = 4 * math.pi * R_ast**2 * BB_monochromatic(T_ast, wavelength) * emissivity
	
	return (L_ast, L_reflected_sun, L_ast/L_reflected_sun)


# Examples
# wavelength = 10 microns
# In [451]: LumRatio(10*u.micron)
# Out[451]: 
# (<Quantity 1.17616908e-06 W / Hz>,
#  <Quantity 3.01199752e-10 W / Hz>,
#  <Quantity 3904.94703845>)

# wavelength = 5 microns
# In [450]: LumRatio(5*u.micron)
# Out[450]: 
# (<Quantity 7.71157157e-08 W / Hz>,
#  <Quantity 1.0561264e-09 W / Hz>,
#  <Quantity 73.0175058>)