5. Correcting Observations

Examples are provided on how to correct astronomical observations for atmospheric effects. For demonstration purposes, the following examples exclusively consider correcting observations of a black body.

5.1. Spectrographic Observations

In order to correct spectrographic observations for atmospheric effects, observed spectra are divided by the atmospheric transmission function. As an example, consider the spectral energy distribution (SED) of a black body under the effects of atmospheric absorption. For a temperature of 8,000 K, and a PWV level of 15 mm, this can be found as:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
>>> from pwv_kpno import blackbody_with_atm as bb_atm
>>> import numpy as np
>>>
>>> temp = 8000
>>> wavelengths = np.arange(7000, 10000, 100)
>>> pwv = 15
>>>
>>> sed = bb_atm.sed(temp, wavelengths, pwv)
>>> print(sed)

  [5406750.98492275, 5646734.81222486, ...,
   2460033.92176797  2383125.54323284]

Since we know the PWV level used to generate this spectra, we can determine the corresponding transmission function using the trans_for_pwv function (although in a real world setting the trans_for_date function may be more appropriate).

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
>>> from pwv_kpno import pwv_atm
>>>
>>> modeled_trans = pwv_atm.trans_for_pwv(15)
>>> print(modeled_trans)

          wavelength        transmission
           Angstrom
      ------------------ ------------------
                  7000.0 0.9198165020431311
       7001.000333444482 0.8590254136327501
       7002.000666888963 0.9967736288238565
                     ...                ...

In order to divide these two results, the SED at and the transmission function must be known for the same wavelength values. Since the SED is a well behaved function, we interpolate the SED to match the wavelength sampling of the transmission function.

Using the numpy package, we interpolate as follows:

1
2
3
>>> sampled_sed = np.interp(modeled_trans["wavelength"],
>>>                         wavelengths,
>>>                         sed)

The corrected spectrum can then be found by dividing the observed flux by the transmission on a wavelength by wavelength basis.

1
>>> corrected_spec = np.divide(sampled_sed, modeled_trans["transmission"])

5.2. Photometric Observations

For an atmospheric transmission \(T(\lambda)\), the photometric correction for an object with an intrinsic spectral energy distribution \(S(\lambda)\) is given by

\[C = \frac{\int S(\lambda) \cdot T(\lambda) \, d\lambda} {\int S(\lambda) \, d\lambda}\]

where the integration bounds are defined by the wavelength range of the photometric bandpass. In practice an SED of the desired object may not be available. In such a case spectral templates should be used instead.

Note that product in the numerator \(S(\lambda) \cdot T(\lambda)\) represents the SED under the influence of atmospheric effects. For a photometric observation taken in the i band (7,000 to 8,500 \(\unicode{x212B}\)), the integration arguments can be found as:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
>>> from pwv_kpno import blackbody_with_atm as bb_atm
>>> import numpy as np
>>>
>>> pwv = 15
>>> temp = 8000
>>>
>>> i_band = (7000, 8500)
>>> sample_rate = 100
>>> wavelengths = np.arange(i_band[0], i_band[1], sample_rate)
>>>
>>> sed_with_atm = bb_atm.sed(temp, wavelengths, pwv)  # S(lambda) T(lambda)
>>> intrinsic_sed = bb_atm.sed(temp, wavelengths, 0)  # S(lambda)

Trapezoidal integration of array like objects in Python can be performed using the Numpy package. Using results from the spectrographic example we have:

1
2
3
4
>>> numerator = np.trapz(sed_with_atm, i_band, sample_rate)
>>> denominator = np.trapz(intrinsic_sed, i_band, sample_rate)
>>>
>>> photo_corr = np.divide(numerator, denominator)