Radiative transfer analysis

There are numerous available radiative transfer models, yet these packages typically require complex installation and compilation procedures for them to operate, and they are particularly restrictive in operational scope (e.g., planet, type of atmosphere/surface) and wavelength (e.g., spectral database, file formats). For instance at the internet encyclopedia (Wikipedia), there is an entry for "Atmospheric radiative transfer codes", listing dozens of packages and their capabilities (wavelength range, geometry, scattering, polarization, accessibility/licensing, etc.). There have been several attempts to quantify the differences between different packages (e.g., Alvarado et al. 2013), and a commercial internet facility (www.spectralcalc.com) implements a small subset relevant to Earth/Mars science.
With this online tool, we pursue to provide to the community an accurate and comprehensive (UV/Vis/IR/mm) planetary tool, allowing to synthesize realistic spectra for a broad range of objects (rocky planets, gas giants, icy bodies, TNOs, comets, exoplanets, etc.), and assisting with the planning and execution of current and future NASA missions. This is achieved by integrating two core radiative-transfer models: the accurate and versatile PUMAS atmospheric/scattering model, and the cometary emission model (CEM).
Radiative transfer models
Components of the PSG radiative transfer: the diagram shows the different components considered by the radiative transfer modules. By performing a layer-by-layer analysis, PUMAS intrinsically integrates and calculates the different flux contributions across the wavelength grid. For comets, the molecular calculation is performed separately by CEM from the surface fluxes, and later added to compute integrated fluxes.

Atmospheric radiative transfer models

In order to compute atmospheric transmittances and radiances, the molecular modules ingest the parameters defined in the "atmosphere" section and compute observable fluxes:
Planetary and Universal Model of Atmospheric Scattering
PUMAS integrates the latest radiative-transfer methods and spectroscopic parameterizations, in order to compute high resolution spectra via line-by-line calculations, and utilizes the efficient correlated-k method at moderate resolutions. The scattering analysis is based on a Martian scattering model (Smith et al. 2009), while the line-by-line calculations have been validated and benchmarked with the accurate GENLN2 model (Edwards 1992). Villanueva, G. L. et al., Science, Volume 348, Issue 6231, pp. 218-221 (2015).
Smith, M. D. et al., JGR, Vol. 114, E00D03, (2009).
Edwards, D. P., "GENLN2: A general line-by-line atmospheric transmittance and radiance model, Version 3.0 description and users guide", NCAR/TN-367-STR, National Center for Atmospheric Research, Boulder, Co. (1992).
Cometary Emission Model
CEM incorporates excitation processes via non-LTE line-by-line fluorescence model at short wavelengths (employing GSFC databases), and ingests HITRAN, JPL and CDMS spectral databases to compute line-by-line LTE fluxes. It operates with expanding coma atmospheres and temperatures lower than 300K, and accurately computes photodissociation processes for parent and daugther species released in the coma. Villanueva, G. L. et al., The molecular composition of Comet C/2007 W1 (Boattini): Evidence of a peculiar outgassing and a rich chemistry. Icarus, Volume 216, Issue 1, p. 227-240 (2011)
Villanueva, G. L., The High Resolution Spectrometer for SOFIA-GREAT: Instrumentation, Atmospheric Modeling and Observations. PhD Thesis, Albert-Ludwigs-Universitaet zu Freiburg, ISBN 3-936586-34-9, Copernicus GmbH Verlag (2004)

Surface models

The surface defines one of the "boundary" conditions of the radiative-transfer analysis. For expanding atmospheres, these parameters relate to those of the nucleus and to the dust particles. Currently, PSG describes the surfaces as Lambertian (with the corresponding albedo) and emitting as a black-body (with the corresponding emissivity). We have developed a versatile surface module that combines a realistic Hapke scattering model (Protopapa et al. 2017) and the capability to ingest a broad range of optical constants, permitting PSG to accurately compute surface reflectances and emissitivities. The surface model considers areal mixing and a one-term Henyey-Greenstein phase funtion.

Instrument parameters

In this section, the user defines the desired characteristics of the synthetic spectra (wavelength range, resolution, desired radiance flux, noise performance, etc.). When observing with ground-based observatories, PSG allows to affect the synthetic spectra by telluric absorption. The tool has access to a database of telluric transmittances pre-computed for 5 altitudes and 4 columns of water for each case (20 cases in total). The tool can also perform noise (and signal-to-noise ratio) calculation by providing details about the detector and the telescope performance. The FOV at diffraction limit is computed for the center wavelength.
PSG allows to define three type of telescope/instrument modes: a) single monolithic telescope, b) interferometric array, and c) a coronograph instrument/telescope. In all cases, the integration of the fluxes is done over bounded and finite field-of-views and spectral ranges, with no convolutions applied to the fields. The model for the coronagraph is extremely simplified, and it is mainly intended for identifying regimes of operation - it assumes that the throughput is minimum (1/contrast) within half the inner-working-angle (IWA), it reaches 50% at the IWA, and the throughput is maximum (100%) at 1.5 times the IWA.
We have compiled a database of instrument models for a diverse range of telescope/instrument combinations, assisting the user when defining the basic parameters of the instrument. Additional instrument models are being developed.
Instrument simulator
Instrument model: We have identified 14 key parametes that are sufficient to describe the overall capabilities and performance of a particular telescope/instrument combination, and we are now in the process of compiling a database of instrument models, that will assist the user when defining the basic parameters of the instrument.

We have developed an echelle simulator for iSHELL (Immersion Grating Echelle Spectrograph) at the NASA Infrared Telescope Facility (IRTF). It employs a set of ab-initio grating equations and simplified parameters that are accurate enough to describe the overall behavior of the cross-dispersors and the gratings.


Radiance and wavelength units

PSG allows to compute synthetic fluxes in a broad range of possible units. The constants used for the conversion are: λ is the wavelength in microns (μm), c is the speed of light (299792458 m/s), ASR is arcseconds2 per steradian (4.2545166E+10), h is Planck's constant (6.6260693E-34 W s2), k is Boltzmann's constant (1.380658E-23 J/K), ATele is the total collecting area of the observatory (m2, nTele⋅π⋅[DTele/2]2), Ω is the field-of-view of the observations (steradian).
Spectral radiance
W / sr / m2 / μmL - This is the intrinsic unit of the modulesWsrm2um
W / sr / m2 / cm-1L' = L ⋅ λ2 / 1E4Wsrm2cm
W / sr / m2 / HzL' = L ⋅ λ2 / (1E6 ⋅ c)Wsrm2Hz
Jy / arcsec2L' = L ⋅ λ2 / (ASR ⋅ c) ⋅ 1E20Jyarc
K (brightness)L' = PT / ln(1 + (PF/L))
PT = 1E6 ⋅ hc/kλ
PF = 2E24 ⋅ hc2 / λ5
K (Rayleigh-Jeans)L' = L ⋅ 1E-18 ⋅ λ4 / 2kcKRJ
W / sr / m2L' = L ⋅ dλWsrm2
Spectral intensity
W / sr / μmI = L ⋅ ATeleWsrum
W / sr / cm-1I = L ⋅ ATele ⋅ λ2 / 1E4Wsrcm
Radiant intensity
W / srI = L ⋅ ATele ⋅ dλWsr
Spectral flux
W / μmF = L ⋅ ATele ⋅ ΩWum
W / cm-1F = L ⋅ ATele ⋅ λ2 / 1E4 ⋅ ΩWcm
Radiant flux
WF = L ⋅ ATele ⋅ Ω ⋅ dλW
Photons / sF = L ⋅ ATele ⋅ Ω ⋅ dλ ⋅ λ ⋅ 1E-6 / hcph
Integrated flux
PhotonsF = L ⋅ texp ⋅ nexp ⋅ ATele ⋅ Ω ⋅ dλ ⋅ λ ⋅ 1E-6 / hcpt
Photons measuredF = L ⋅ neff ⋅ texp ⋅ nexp ⋅ ATele ⋅ Ω ⋅ dλ ⋅ λ ⋅ 1E-6 / hcpm
W / m2E = L ⋅ dλ ⋅ ΩWm2
erg/s / cm2E = L ⋅ dλ ⋅ Ω ⋅ 1E3erg
Spectral irradiance
W / m2 / μmE = L ⋅ ΩWm2um
W / m2 / cm-1E = L ⋅ Ω ⋅ λ2 / 1E4 ⋅Wm2cm
JyE = L ⋅ Ω ⋅ λ2 / c ⋅ 1E20Jy
mJyE = L ⋅ Ω ⋅ λ2 / c ⋅ 1E23mJy
ContrastT = L / Lstellarrel
Transmittance (I/F)T = L / Lcontrif
Magnitudem = -2.5 ⋅ log(L ⋅ Ω / LVega)V


Online unit conversion tool

Radiation value:
From radiation unit:

To radiation unit:



Resolving power:


Telescope diameter [m, effective]:


Frequency / wavelength value:
From unit:       To unit:


Telluric absorption

When observing with ground-based observatories, PSG allows to affect the synthetic spectra (or simply show) by telluric absorption. The tool has access to a database of telluric transmittances pre-computed for 5 altitudes and 4 columns of water for each case (20 cases in total). The altitudes include that of Mauna-Kea/Hawaii (4200 m), Paranal/Chile (2600 m), SOFIA (14,000 m) and balloon observatories (35,000 m), while the water vapor column was established by scaling the tropical water profile by a factor of 0.1, 0.3 and 0.7 and 1.
Opacities at 225 GHz, a typical metric to quantify water at radio wavelengths, can be estimated from the reported water column as τ225PSG = 0.0642 x PWV, where PWV is the ammount of water in precipitable millimeters.

Noise simulator

PSG currently includes a noise calculator for quantum and thermal detectors, with the primary goal of providing users with quick look simulations for planning observations, and to assist with the development of new instrument/telescope concepts. The user can provide a constant value across all wavelengths (RMS), or a constant value with background noise (BKG), or can choose from several detector noise simulators.
At short wavelengths (e.g., optical or near IR), the background photon counts follow a Poissson distribution, and the fluctuations are given by √N where N is the mean number of photons received (see review in Zmuidzinas et al. 2003). This Poisson distribution holds only in the case that the mean photon mode occupation number is small, n<<1. For a thermal background, the occupation number is given by the Bose-Einstein formula, nth(v,T) = [exp(hv/kT)-1]-1, so the opposite classical limit n>>1 is the usual situation at longer wavelengths for which hv<<kT. When n>>1, the photons do not arrive independently according to a Poisson process but instead are strongly bunched, and the fluctuations are of order N, instead of √N. This is why the Dicke equation is used to calculate sensitivities for the receiver temperature mode (TRX), which states that the noise is proportional to the background power rather than its square root. The formalism employed for the TRX module is based on the ALMA sensitivity calculator (Yatagai et al. 2011).
PSG assumes that the instrument has a defined spatial resolution (beam [FWHM]), defined by the user for the center wavelength. This spatial resolution will change across wavelength (proportional to λ), and therefore the solid-angle (Ω) will be proportional to λ2. Since Ω is changing with wavelength, and the instrument has a fixed spatial resolution (defined by the optical design and detector properites), the number of pixels encompassing Ω will be also proportional to λ2.
Type of noiseParametersDetector specific noise formulas
Receiver temperature (radio)
TRX [K]: receiver temperature
g: sideband factor (0:SSB, 1:DSB)
npol: 1 (number of polarizations)
fN = 1 (number of baselines)
For interferometric systems (e.g., ALMA):
npol: 2 (assuming dual / full configuration)
fN = ntele(ntele-1)
LRJ = 1E-18 ⋅ λ4 / 2kc
Tsource = L ⋅ LRJ
Tback = LbackLRJ + Tground(1 - trnground)
ksys = (1+g)/(ηTotal trnground)
Tsys = ksys [TRX + εopticsToptics + Tsource + Tback]   [K]
fΩ = (ΩTele / Ω)      (Diffraction FOV / observations FOV) correction
dv = 1E6 ⋅ c ⋅ dλ / λ2
Ntotal = Tsys ⋅ fΩ / √(fN ⋅ npol ⋅ dv ⋅ nexp ⋅ texp)   [K]
Noise Equivalent Power
NEP [W / √Hz]: sensitivity
ND = npixels ⋅ nexp ⋅ texp ⋅ (NEP ⋅ λ ⋅ 1E-6 / hc)2   [e-2]
D* - Detectivity D* [cm √Hz / W]: detectivity
S [μm]: pixel size
ND = npixels ⋅ nexp ⋅ texp ⋅ ((S ⋅ 1E-4 / D*) ⋅ λ ⋅ 1E-6 / hc)2
Imager - Charge image sensor (e.g., CCD, CMOS, EMCCD, ICCD / MCP) Read-noise [e- / pixel]: read noise
Dark [e- / s / pixel]: dark current
ND = npixels ⋅ nexp ⋅ [ Nread2 + (Dark ⋅ texp) ]   [e-2]

The noise components with Poisson statistics (i.e., UV, optical, IR) are calculated as:
Le- = Ω ⋅ ATele ⋅ ηeff ⋅ dλ ⋅ λ ⋅ texp ⋅ nexp ⋅ 1E-6 / hc        Radiance to detector electrons conversion factor
Nsource = L ⋅ Le-        Noise introduced by the source itself [e-2]
Nback = (Lback + nezo⋅Lzodi) ⋅ Le-        Noise introduced by background sky sources [e-2]
Noptics = εoptics ⋅ Le- ⋅ (2E24 ⋅ h ⋅ c2 / λ5) / (exp(1E6 ⋅ h ⋅ c / (k ⋅ Toptics ⋅ λ)) - 1)        Noise introduced by the telescope [e-2]
Nground = Le- ⋅ (1 - trnground) ⋅ (2E24 ⋅ h ⋅ c2 / λ5) / (exp(1E6 ⋅ h ⋅ c / (k ⋅ Tground ⋅ λ)) - 1)    Noise for ground observations [e-2]
NTotal = √(ND + Nsource + Nback + Noptics + Nground)        Total noise [e-]
Assuming these parameters and units:
    L [W / sr / m2 / μm]: spectral radiance of the source
    Lback [W / sr / m2 / μm]: spectral radiance of the background sources
    texp [s]: time per exposure
    nexp: total number of exposures
    npixels: total number of pixels for Ω and dλ.
    nezo: Exozodiacal dust scaler relative to Solar System zodiacal dust
    Toptics [K]: temperature of the optics
    εoptics: emissivity of the optics
    ηeff: total throughput of the system (including quantum efficiencies)
    Ω [steradian]: is the solid angle of the observations. It is wavelength dependent.
    ATele [m2]: is the total collecting area of the observatory (nTele⋅π⋅[DTele/2]2)
    λ [μm]: is the wavelength in microns
    trnground: terrestrial transmittance
    Tground [K]: temperature of the terrestrial atmosphere - 280
    h [W s2]: is Planck's constant - 6.6260693E-34
    c [m / s]: is the speed of light - 299792458
    k [J / K]: is Boltzmann's constant - 1.380658E-23
Background noise sources
Background noise sources: When observing faint astronomical sources, the sensitivity is affected by the shot noise introduced by background and diffuse sources (Leiner et al. 1998, A&A, v.127, p.1-99). From space, the background is dominated by the faint and diffuse emission (thermal and scattered sunlight) from zodiacal dust, while airglow (a mixture of photoionization emissions, chemiluminescence and scattered sunlight) dominates the background for ground-based observations. Zodiacal dust fluxes depend greatly on the ecliptic longitude/latitude - in PSG, the noise simulator considers a scaling of 2 with respect to the minimum ecliptic pole values. PSG also employs a rudimentary (as shown), yet relatively effective, approximation for atmospheric airglow.