Source code for gwcs.spectroscopy

# Licensed under a 3-clause BSD style license - see LICENSE.rst
"""
Spectroscopy related models.
"""

import numpy as np
from astropy.modeling.core import Model
from astropy.modeling.parameters import Parameter
import astropy.units as u


__all__ = ['WavelengthFromGratingEquation', 'AnglesFromGratingEquation3D',
           'Snell3D', 'SellmeierGlass', 'SellmeierZemax']


[docs]class WavelengthFromGratingEquation(Model): r""" Solve the Grating Dispersion Law for the wavelength. .. Note:: This form of the equation can be used for paraxial (small angle approximation) as well as oblique incident angles. With paraxial systems the inputs are ``sin`` of the angles and it transforms to :math:`(\sin(alpha_in) + \sin(alpha_out)) / (groove_density * spectral_order)`. With oblique angles the inputs are the direction cosines of the angles. Parameters ---------- groove_density : int Grating ruling density in units of 1/length. spectral_order : int Spectral order. Examples -------- >>> from astropy.modeling.models import math >>> model = WavelengthFromGratingEquation(groove_density=20000*1/u.m, spectral_order=-1) >>> alpha_in = (math.Deg2radUfunc() | math.SinUfunc())(.0001 * u.deg) >>> alpha_out = (math.Deg2radUfunc() | math.SinUfunc())(.0001 * u.deg) >>> lam = model(alpha_in, alpha_out) >>> print(lam) -1.7453292519934437e-10 m """ _separable = False linear = False n_inputs = 2 n_outputs = 1 groove_density = Parameter(default=1) """ Grating ruling density in units of 1/m.""" spectral_order = Parameter(default=1) """ Spectral order.""" def __init__(self, groove_density, spectral_order, **kwargs): super().__init__(groove_density=groove_density, spectral_order=spectral_order, **kwargs) self.inputs = ("alpha_in", "alpha_out") """ Sine function of the angles or the direction cosines.""" self.outputs = ("wavelength",) """ Wavelength."""
[docs] def evaluate(self, alpha_in, alpha_out, groove_density, spectral_order): return (alpha_in + alpha_out) / (groove_density * spectral_order)
@property def return_units(self): if self.groove_density.unit is None: return None return {'wavelength': u.Unit(1 / self.groove_density.unit)}
[docs]class AnglesFromGratingEquation3D(Model): """ Solve the 3D Grating Dispersion Law in Direction Cosine space for the refracted angle. Parameters ---------- groove_density : int Grating ruling density in units of 1/m. order : int Spectral order. Examples -------- >>> from astropy.modeling.models import math >>> model = AnglesFromGratingEquation3D(groove_density=20000*1/u.m, spectral_order=-1) >>> alpha_in = (math.Deg2radUfunc() | math.SinUfunc())(.0001 * u.deg) >>> beta_in = (math.Deg2radUfunc() | math.SinUfunc())(.0001 * u.deg) >>> lam = 2e-6 * u.m >>> alpha_out, beta_out, gamma_out = model(lam, alpha_in, beta_in) >>> print(alpha_out, beta_out, gamma_out) 0.04000174532925199 -1.7453292519934436e-06 0.9991996098716049 """ _separable = False linear = False n_inputs = 3 n_outputs = 3 groove_density = Parameter(default=1) """ Grating ruling density in units 1/ length.""" spectral_order = Parameter(default=1) """ Spectral order.""" def __init__(self, groove_density, spectral_order, **kwargs): super().__init__(groove_density=groove_density, spectral_order=spectral_order, **kwargs) self.inputs = ("wavelength", "alpha_in", "beta_in") """ Wavelength and 2 angle coordinates going into the grating.""" self.outputs = ("alpha_out", "beta_out", "gamma_out") """ Two angles coming out of the grating. """
[docs] def evaluate(self, wavelength, alpha_in, beta_in, groove_density, spectral_order): if alpha_in.shape != beta_in.shape: raise ValueError("Expected input arrays to have the same shape.") if isinstance(groove_density, u.Quantity): alpha_in = u.Quantity(alpha_in) beta_in = u.Quantity(beta_in) alpha_out = -groove_density * spectral_order * wavelength + alpha_in beta_out = - beta_in gamma_out = np.sqrt(1 - alpha_out ** 2 - beta_out ** 2) return alpha_out, beta_out, gamma_out
@property def input_units(self): if self.groove_density.unit is None: return None return { 'wavelength': 1 / self.groove_density.unit, 'alpha_in': u.Unit(1), 'beta_in': u.Unit(1) }
[docs]class Snell3D(Model): """ Snell model in 3D form. Inputs are index of refraction and direction cosines. Returns ------- alpha_out, beta_out, gamma_out : float Direction cosines. """ _separable = False linear = False n_inputs = 4 n_outputs = 3 def __init__(self, **kwargs): super().__init__(**kwargs) self.inputs = ('n', 'alpha_in', 'beta_in', 'gamma_in') self.outputs = ('alpha_out', 'beta_out', 'gamma_out')
[docs] @staticmethod def evaluate(n, alpha_in, beta_in, gamma_in): # Apply Snell's law through front surface, # eq 5.3.3 II in Nirspec docs alpha_out = alpha_in / n beta_out = beta_in / n gamma_out = np.sqrt(1.0 - alpha_out**2 - beta_out**2) return alpha_out, beta_out, gamma_out
[docs]class SellmeierGlass(Model): """ Sellmeier equation for glass. Parameters ---------- B_coef : ndarray Iterable of size 3 containing B coefficients. C_coef : ndarray Iterable of size 3 containing c coefficients in units of `u.um**2`. Returns ------- n : float Refractive index. Examples -------- >>> import astropy.units as u >>> b_coef = [0.58339748, 0.46085267, 3.8915394] >>> c_coef = [0.00252643, 0.010078333, 1200.556] * u.um**2 >>> model = SellmeierGlass(b_coef, c_coef) >>> model(2 * u.m) <Quantity 2.43634758> References ---------- .. [1] https://en.wikipedia.org/wiki/Sellmeier_equation Notes ----- Model formula: .. math:: n(\\lambda)^2 = 1 + \\frac{(B1 * \\lambda^2 )}{(\\lambda^2 - C1)} + \\frac{(B2 * \\lambda^2 )}{(\\lambda^2 - C2)} + \\frac{(B3 * \\lambda^2 )}{(\\lambda^2 - C3)} """ _separable = False standard_broadcasting = False linear = False n_inputs = 1 n_outputs = 1 B_coef = Parameter(default=np.array([1, 1, 1])) """ B1, B2, B3 coefficients. """ C_coef = Parameter(default=np.array([0, 0, 0])) """ C1, C2, C3 coefficients in units of um ** 2. """ def __init__(self, B_coef, C_coef, **kwargs): super().__init__(B_coef, C_coef) self.inputs = ('wavelength',) self.outputs = ('n',)
[docs] @staticmethod def evaluate(wavelength, B_coef, C_coef): B1, B2, B3 = B_coef[0] C1, C2, C3 = C_coef[0] n = np.sqrt(1. + B1 * wavelength ** 2 / (wavelength ** 2 - C1) + B2 * wavelength ** 2 / (wavelength ** 2 - C2) + B3 * wavelength ** 2 / (wavelength ** 2 - C3) ) return n
@property def input_units(self): if self.C_coef.unit is None: return None return {'wavelength': u.um}
[docs]class SellmeierZemax(Model): """ Sellmeier equation used by Zemax. Parameters ---------- temperature : float Temperature of the material in `u.Kelvin`. ref_temperature : float Reference emperature of the glass in `u.Kelvin`. ref_pressure : float Reference pressure in ATM. pressure : float Measured pressure in ATM. B_coef : ndarray Iterable of size 3 containing B coefficients. C_coef : ndarray Iterable of size 3 containing C coefficients in units of `u.um**2`. D_coef : ndarray Iterable of size 3 containing constants to describe the behavior of the material. E_coef : ndarray Iterable of size 3 containing constants to describe the behavior of the material. Returns ------- n : float Refractive index. """ _separable = False standard_broadcasting = False linear = False n_inputs = 1 n_outputs = 1 temperature = Parameter(default=0) ref_temperature = Parameter(default=0) ref_pressure = Parameter(default=0) pressure = Parameter(default=0) B_coef = Parameter(default=[1, 1, 1]) C_coef = Parameter(default=[0, 0, 0]) D_coef = Parameter(default=[0, 0, 0]) E_coef = Parameter(default=[1, 1, 1]) def __init__(self, temperature=temperature, ref_temperature=ref_temperature, ref_pressure=ref_pressure, pressure=pressure, B_coef=B_coef, C_coef=C_coef, D_coef=D_coef, E_coef=E_coef, **kwargs): super().__init__(temperature=temperature, ref_temperature=ref_temperature, ref_pressure=ref_pressure, pressure=pressure, B_coef=B_coef, C_coef=C_coef, D_coef=D_coef, E_coef=E_coef, **kwargs) self.inputs = ('wavelength',) self.outputs = ('n',)
[docs] def evaluate(self, wavelength, temp, ref_temp, ref_pressure, pressure, B_coef, C_coef, D_coef, E_coef): """ Input ``wavelength`` is in units of microns. """ if isinstance(temp, u.Quantity): temp = temp.to(u.Celsius) ref_temp = ref_temp.to(u.Celsius) else: KtoC = 273.15 # kelvin to celcius conversion temp -= KtoC ref_temp -= KtoC delt = temp - ref_temp D0, D1, D2 = D_coef[0] E0, E1, lam_tk = E_coef[0] nref = 1. + (6432.8 + 2949810. * wavelength ** 2 / (146.0 * wavelength ** 2 - 1.) + (5540.0 * wavelength ** 2) / (41.0 * wavelength ** 2 - 1.)) * 1e-8 # T should be in C, P should be in ATM nair_obs = 1.0 + ((nref - 1.0) * pressure) / (1.0 + (temp - 15.) * 3.4785e-3) nair_ref = 1.0 + ((nref - 1.0) * ref_pressure) / (1.0 + (ref_temp - 15) * 3.4785e-3) # Compute the relative index of the glass at Tref and Pref using Sellmeier equation I. lamrel = wavelength * nair_obs / nair_ref nrel = SellmeierGlass.evaluate(lamrel[0], B_coef, C_coef) # Convert the relative index of refraction at the reference temperature and pressure # to absolute. nabs_ref = nrel * nair_ref # Compute the absolute index of the glass delnabs = (0.5 * (nrel ** 2 - 1.) / nrel) * \ (D0 * delt + D1 * delt ** 2 + D2 * delt ** 3 + \ (E0 * delt + E1 * delt ** 2) / (lamrel ** 2 - lam_tk ** 2)) nabs_obs = nabs_ref + delnabs # Define the relative index at the system's operating T and P. n = nabs_obs / nair_obs return n