Calculate interstellar reddening from the Na D1 line strength.

gistfile1.py

from numpy import arange
from math import fractorial

def calculate_interstellar_extinction(Na_D_equivalent_width, num_terms=100, precision=0.01):
	"""Calculates interstellar extinction from the strength of the Na D1 line
	at 5889.95 \AA based on the relationship found by Munari & Zwitter (2006)."""


	alpha, err_alpha = 0.354, 0.01
	beta, err_beta = 11.0, 1.0

	def calculate_expected_ew(Ebmv, num_terms, coefficients):
		
		alpha, beta = coefficients

		equivalent_width = 0
		for n in xrange(1, num_terms + 1):
			equivalent_width += alpha * pow(-1, n - 1) * (pow(beta * Ebmv, n)  / (factorial(n) * pow(n, 0.5)))

		return equivalent_width


	def calculate_ebmv(Na_D_equivalent_width, ebmv_range, num_terms, coefficients):

		Ebmv_values = arange(*ebmv_range)
		for Ebmv in Ebmv_values:
			if calculate_expected_ew(Ebmv, num_terms=num_terms, coefficients=coefficients) > Na_D_equivalent_width:
				break

		if Ebmv == Ebmv_values[-1]:
			raise ValueError("Could not find reddenning - expected value exceeds %1.2f" % (Ebmv, ))

		return Ebmv


	# Get the value
	Ebmv = calculate_ebmv(Na_D_equivalent_width, (0.0, 2.0, precision), num_terms, (alpha, beta))

	# Uncertainties

	Ebmv_alpha = calculate_ebmv(Na_D_equivalent_width, (0.0, 2.0, precision), num_terms, (alpha + err_alpha, beta + err_beta))
	Ebmv_beta = calculate_ebmv(Na_D_equivalent_width, (0.0, 2.0, precision), num_terms, (alpha - err_alpha, beta - err_beta))

	return Ebmv, Ebmv_alpha - Ebmv, Ebmv_beta - Ebmv