<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 On-the-fly mapping exposure calculator (Radio astronomy)
Step1:  Introduction
Step2:  Atmospheric temperature is approximately given by ambient temperature at ground.
Step3:  Calculate telescope sensitivity (aka Kelvins per Jansky).
Step4:  The conversion between antenna temperatures and main-beam brightness temperatures is given by
Step5:  Define spectrometer properties
Step6:  If a certain velocity resolution is desired, we first have to infer the desired spectral resolution.
Step7:  This means, we can bin the original spectrum by a factor of
Step8:  which will further decrease the noise.
Step9:  Case 1
Step10:  Note, Tsys is higher for low elevation (more air mass).
Step11:  Calculate raw $T_\mathrm{A}$ noise
Step12:  For dual polarization observations, we can divide by $\sqrt{2}$.
Step13:  We also have to account for the position switch (division by noisy reference spectrum).
Step14:  Finally, we convert to main-beam brightness temperature, $\Delta T_\mathrm{B}$, and flux-density, $\Delta S$, noise
Step15:  The astronomical signal is furthermore attenuated by the atmosphere. There are two ways to handle this. (a) If the signal strength is known (or can be expected to have a certain value), just apply the attenuation factor and compare to the noise levels. (b) calculate an effective noise level, by increasing Tb and flux-density noise accordingly. Here, we follow the second approach, as the true signal level is unknown. The effective RMS values are not true RMS estimates, but merely serve to indicate the impact of Earth's atmosphere on the sensitivity. Close to zenith, the effect is a few percent only, but for very low elevations, almost half of the signal is lost!
Step16:  Case 2
Step17:  Again, divide by two for dual polarization
Step18:  and account for PSW (division by noisy reference spectrum)
Step19:  On-the-fly mapping
For OTF mapping, it is not straight-forward to calculate the effective noise in the map after gridding. Therefore, we simulate the process, use cygrid to produce a map, and will measure the noise accordingly.
Map setup
In Effelsberg, we do so-called zig-zag scanning. For position switching, one observes a reference position every $n$ scanlines. In contrast to the pointed observations, it is possible to spend a different amount of time on the reference position (this will have an impact on the degree of correlated spatial noise in the final map!).
Step20:  The spacing between the scans should not be larger than half the beamwidth, a third is even better.
Step21:  Likewise, the scan velocity must not be too high. This is determined by the sampling rate and the beam size (one wants at least three independent samples, per beam, otherwise the resulting map will have some smearing along the scan lines).
Let's calculate the maximal scan speed
Step22:  With this, we can calculate the minimal duration of one scan line
Step23:  However, each telescope will have a duty cycle between two scan lines, which is defined by the time the telescope needs for re-positioning. At Effelsberg, this is about 15 s, which means that the duration per scanline should at least be one or two minutes, otherwise the observing efficiency would become really poor.
Choose a meaningful scanline duration
Step24:  We now have to choose the time spent on the reference position. Furthermore, it is possible to do the reference scan only after every $n$ scanlines (to save time).
Let's do a relatively long refpos integration, but only after every 2 scanlines (avoids some duty cycles).
Step25:  Such a map would need the following total observing time.
Step26:  Now, we create the raw data. Note, the absolute value of the RMS in counts (for the P quantities) is not important, but the RMS ratio between On and Ref spectra is (the absolute number is unimportant, because it gets calibrated away in the (On-Ref/Ref) equation).
Compared to the On scan, a Ref position is integrated over 'refpos_duration', which means that noise in the Ref spectrum is $\sqrt{\texttt{refpos_duration} / \texttt{sampling_interval}}$ smaller.
To increase the noise estimate accuracy, we will produce several noise maps at once.
Step27:  Here, a dummy $T_\mathrm{sys}$ value is used. Later, the map (and thus the RMS) can simply be scaled to match the true $T_\mathrm{sys}$.
Step28:  We will now create spectral noise (in arbitrary units). This must be much smaller than $T_\mathrm{sys}$ to avoid numerical problems. It is also important, that each of the raw-data spectra has an offset (the $T_\mathrm{sys}$ level in counts, if you will). If one would neglect this, one had a division by very small numbers in the line where the reduced spectra are calculated. The magnitude of the noise is furthermore tied to the $T_\mathrm{sys}$ level - it is based on the "radiometer" factor $\left(1/\sqrt{\tau\Delta f}\right)$
Step29:  We can test this, by just plotting the "raw" data for one channel.
Step30:  Now do the real gridding. First, prepare a WCS header
Step31:  We already define a WCS object for later use in our plots
Step32:  Setup the gridder and define kernel sizes (half the beamsize is always a good choice).
Step33:  The gridder needs the coordinates as flat arrays. The data to be gridded must be a 2D array (first dimension has to match the number of coordinate samples, second dimension is the number of channels/maps in the desired data cube).
Step34:  Again, as a sanity check, we plot one of the channel maps.
Step35:  Last but not least, we can measure the noise. There are two possibilities to do this.
First, we can calculate the RMS over the full data cube.
Step36:  Second, we can calculate the RMS per plane (and the average over all planes)
Step37:  The ultimate questions is, what is the noise level in such a map, if we account for real $T_\mathrm{sys}$ and atmospheric effects, etc.
Step38:  $T_\mathrm{A}$ noise is now the measured noise from the gridded maps multiplied with the $T_\mathrm{sys}$
Step39:  Note, we don't need to account for the position switch again (division by noisy reference spectrum was already performed in the reduced-spectra array creation)!
<ASSISTANT_TASK:>
Python Code:
%load_ext autoreload
%autoreload 2
%matplotlib inline
%config InlineBackend.figure_format = 'retina'
import numpy as np
import matplotlib.pyplot as pl
from astropy.wcs import WCS
from scipy import constants
import cygrid
np.set_printoptions(precision=1)
def gaincurve(elev, a0, a1, a2):
    '''
    Radio telescope sensitivity is usually a function of elevation (parametrized as parabula).
    '''
    return a0 + a1 * elev + a2 * elev * elev
def dv_to_df(restfreq, velo_kms):
    '''
    Convert velocity resolution to frequency resolution.
    '''
    return restfreq * velo_kms * 1.e3 / constants.c
def setup_header(mapcenter, mapsize, beamsize_fwhm):
    '''
    Produce a FITS header that contains the target field.
    '''
    
    # define target grid (via fits header according to WCS convention)
    # a good pixel size is a third of the FWHM of the PSF (avoids aliasing)
    pixsize = beamsize_fwhm / 3.
    dnaxis1 = int(mapsize[0] / pixsize)
    dnaxis2 = int(mapsize[1] / pixsize)
    header = {
        'NAXIS': 2,
        'NAXIS1': dnaxis1,
        'NAXIS2': dnaxis2,
        'CTYPE1': 'RA---SIN',
        'CTYPE2': 'DEC--SIN',
        'CUNIT1': 'deg',
        'CUNIT2': 'deg',
        'CDELT1': -pixsize,
        'CDELT2': pixsize,
        'CRPIX1': (dnaxis1 + 1) / 2.,
        'CRPIX2': (dnaxis2 + 1) / 2.,
        'CRVAL1': mapcenter[0],
        'CRVAL2': mapcenter[1],
        }
    
    return header
Explanation: On-the-fly mapping exposure calculator (Radio astronomy)
End of explanation
dual_pol = True
restfreq = 23.7e9  # Hz
opacity = 0.07  # assume reasonably good weather
Tsys_zenith = 60.
Explanation: Introduction: Position switching
Let's start with an illustrative example, that shows how exposure calculation is done for a pointed observation (position switch). Usually, the backend provides the measured (spectral) intensities, $P$, in an uncalibrated fashion, e.g., in units of counts. Therefore one needs to derive a conversion factor. Furthermore, one has to remove the influence of the system bandpass (frequency-dependent gain). Both problems can be solved with the so-called position-switching technique:
\begin{equation}
T = T_\mathrm{sys}\frac{P_\mathrm{on} - P_\mathrm{ref}}{P_\mathrm{ref}}\,.
\end{equation}
The system temperature, $T_\mathrm{sys}$, is a property of the receiver, which determines the base noise level. It depends on the receiver noise itself, but also has conributions from ground and atmosphere, as well as astronomical background (Galactic continuum, CMB). The noise level, $\Delta T$, of the derived spectral intensity, $T$, will decrease the longer one integrates (radiometer equation):
\begin{equation}
\Delta T = \frac{T_\mathrm{sys}}{\sqrt{\tau\Delta f}}\,.
\end{equation}
However, for position switching one divides the On and Off spectra, which increases the noise by a factor of $\sqrt{2}$. This factor is absorbed by the fact, that one often has two polarization channels available, which can be averaged.
Example: observation of NH3 lines at 23.7 GHz (at Effelsberg)
For this one would use a K-band receiver, which has about 60 K system temperature (zenith). For non-zenith observations, the airmass will increase the effective system temperature, depending on the atmospheric opacity, $\tau$.
End of explanation
T_amb = 290.  # K
T_atm = T_amb - 17.
Explanation: Atmospheric temperature is approximately given by ambient temperature at ground.
End of explanation
Gamma = 1.12  # K/Jy
eta_MB = 0.79  # main beam efficiency
Gamma_MB = Gamma / eta_MB
Explanation: Calculate telescope sensitivity (aka Kelvins per Jansky).
End of explanation
Ta_to_Tb = 1. / eta_MB
Explanation: The conversion between antenna temperatures and main-beam brightness temperatures is given by
End of explanation
nchan = 2 ** 16  # 64 k
bandwidth = 5e8  # 500 MHz
spec_reso = bandwidth / nchan * 1.16  # true spectral resolution 16% worse than channel width
print('spec_reso = {:.1f} kHz'.format(spec_reso * 1.e-3))
Explanation: Define spectrometer properties
End of explanation
desired_vel_resolution = 1.  # km/s
desired_freq_resolution = dv_to_df(restfreq, desired_vel_resolution)
print('desired_freq_resolution = {:.1f} kHz'.format(desired_freq_resolution * 1.e-3))
Explanation: If a certain velocity resolution is desired, we first have to infer the desired spectral resolution.
End of explanation
smooth_nbin = int(desired_freq_resolution / spec_reso + 0.5)
Explanation: This means, we can bin the original spectrum by a factor of
End of explanation
print('smooth_nbin', smooth_nbin)
Explanation: which will further decrease the noise.
End of explanation
exposure = 60.  # seconds
elevations = np.array([10, 20, 30, 40, 50, 60, 90])
AM = 1. / np.sin(np.radians(elevations))
gain_correction = gaincurve(elevations, 0.954, 3.19E-3, -5.42E-5)
Explanation: Case 1: Calculate the noise after a certain integration time. The result will depend on the elevation of the source (because different airmasses modify $T_\mathrm{sys}$).
End of explanation
Tsys_corr = Tsys_zenith + T_atm * (np.exp(opacity * AM) - np.exp(opacity * 1))
print('Tsys_corr', Tsys_corr)
# Tsys_corr = Tsys_zenith + opacity * T_atm * (AM - 1)  # approximate formula, for small opacity * AM
# print(Tsys_corr)
Explanation: Note, Tsys is higher for low elevation (more air mass).
End of explanation
Ta_rms = Tsys_corr / np.sqrt(spec_reso * smooth_nbin * exposure)
Explanation: Calculate raw $T_\mathrm{A}$ noise:
End of explanation
if dual_pol:
    Ta_rms /= np.sqrt(2.)
Explanation: For dual polarization observations, we can divide by $\sqrt{2}$.
End of explanation
Ta_rms *= np.sqrt(2.)
Explanation: We also have to account for the position switch (division by noisy reference spectrum).
End of explanation
Tb_rms = Ta_to_Tb * Ta_rms / gain_correction
S_rms = Tb_rms / Gamma_MB
Explanation: Finally, we convert to main-beam brightness temperature, $\Delta T_\mathrm{B}$, and flux-density, $\Delta S$, noise:
End of explanation
atm_atten = np.exp(-opacity * AM)
print('{0:>8s} {1:>8s} {2:>10s} {3:>10s} {4:>10s} {5:>10s} {6:>10s} {7:>10s} {8:>10s}'.format(
    'Elev', 'Airmass', 'Tsys', 'Ta RMS', 'Tb RMS', 'S RMS', 'AtmAtten', 'Tb_eff RMS', 'S_eff RMS'
    ))
print('{0:>8s} {1:>8s} {2:>10s} {3:>10s} {4:>10s} {5:>10s} {6:>10s} {7:>10s} {8:>10s}'.format(
    '[d]', '', '[K]', '[K]', '[K]', '[Jy]', '', '[K]', '[Jy]'
    ))
for idx in range(len(elevations)):
    print(
        '{0:>8.2f} {1:>8.2f} {2:>10.4f} {3:>10.4f} {4:>10.4f} '
        '{5:>10.4f} {6:>10.4f} {7:>10.4f} {8:>10.4f}'.format(
            elevations[idx], AM[idx], Tsys_corr[idx],
            Ta_rms[idx], Tb_rms[idx], S_rms[idx],
            atm_atten[idx],
            Tb_rms[idx] / atm_atten[idx], S_rms[idx] / atm_atten[idx],
            ))
print('Ta RMS = Antenna temp. noise')
print('Tb RMS = Brightness temp. noise')
print('S RMS = Flux density noise')
Explanation: The astronomical signal is furthermore attenuated by the atmosphere. There are two ways to handle this. (a) If the signal strength is known (or can be expected to have a certain value), just apply the attenuation factor and compare to the noise levels. (b) calculate an effective noise level, by increasing Tb and flux-density noise accordingly. Here, we follow the second approach, as the true signal level is unknown. The effective RMS values are not true RMS estimates, but merely serve to indicate the impact of Earth's atmosphere on the sensitivity. Close to zenith, the effect is a few percent only, but for very low elevations, almost half of the signal is lost!
End of explanation
Tb_eff_rms_desired = 0.01  # 10 mK
Ta_eff_rms_desired = Tb_eff_rms_desired * gain_correction / Ta_to_Tb
Explanation: Case 2: Calculate, how long we have to integrate, to reach a certain (effective) noise level.
Start with a given noise level
End of explanation
exposure = (Tsys_corr / Ta_eff_rms_desired) ** 2 / (spec_reso * smooth_nbin)
if dual_pol:
    exposure /= 2.
Explanation: Again, divide by two for dual polarization
End of explanation
exposure *= 2.
print('Exposure time needed to reach an effective MB brightness temperature noise level of {:.1f} mK'.format(
    Tb_eff_rms_desired * 1.e3))
print('{0:>8s} {1:>10s}'.format('Elev [d]', 'Time [min]'))
for idx in range(len(elevations)):
    print('{0:>8.2f} {1:>10.1f}'.format(
        elevations[idx], (exposure / 60.)[idx]
        ))
Explanation: and account for PSW (division by noisy reference spectrum)
End of explanation
map_width, map_height = 100., 100.  # arcsec
beamsize_fwhm = 38.  # arcsec; at the frequency given in our example
Explanation: On-the-fly mapping
For OTF mapping, it is not straight-forward to calculate the effective noise in the map after gridding. Therefore, we simulate the process, use cygrid to produce a map, and will measure the noise accordingly.
Map setup
In Effelsberg, we do so-called zig-zag scanning. For position switching, one observes a reference position every $n$ scanlines. In contrast to the pointed observations, it is possible to spend a different amount of time on the reference position (this will have an impact on the degree of correlated spatial noise in the final map!).
End of explanation
num_scan_lines = int(3 * map_height / beamsize_fwhm + 0.5)
print('num_scan_lines', num_scan_lines)
Explanation: The spacing between the scans should not be larger than half the beamwidth, a third is even better.
End of explanation
sampling_interval = 1.  # s (== 4 x 250 ms at Effelsberg)
max_speed = beamsize_fwhm / 3 / sampling_interval
print('max_speed = {:.2f} arcsec per s'.format(max_speed))
Explanation: Likewise, the scan velocity must not be too high. This is determined by the sampling rate and the beam size (one wants at least three independent samples, per beam, otherwise the resulting map will have some smearing along the scan lines).
Let's calculate the maximal scan speed
End of explanation
min_duration = map_width / max_speed
print('min_duration = {:.1f} s'.format(min_duration))
Explanation: With this, we can calculate the minimal duration of one scan line
End of explanation
scanline_duration = 90.  # seconds
samples_per_scanline = int(scanline_duration / sampling_interval + 0.5)
print('samples_per_scanline', samples_per_scanline)
Explanation: However, each telescope will have a duty cycle between two scan lines, which is defined by the time the telescope needs for re-positioning. At Effelsberg, this is about 15 s, which means that the duration per scanline should at least be one or two minutes, otherwise the observing efficiency would become really poor.
Choose a meaningful scanline duration
End of explanation
refpos_duration = 90.  # seconds
refpos_interval = 2 
Explanation: We now have to choose the time spent on the reference position. Furthermore, it is possible to do the reference scan only after every $n$ scanlines (to save time).
Let's do a relatively long refpos integration, but only after every 2 scanlines (avoids some duty cycles).
End of explanation
duty_cycle = 15.  # s
total_on_time = (scanline_duration + duty_cycle) * num_scan_lines
total_ref_time = (refpos_duration + duty_cycle) * (num_scan_lines // refpos_interval)
total_time = total_on_time + total_ref_time
print('Total time necessary for map: {:.1f} min'.format(total_time / 60.))
Explanation: Such a map would need the following total observing time.
End of explanation
num_maps = 128
Explanation: Now, we create the raw data. Note, the absolute value of the RMS in counts (for the P quantities) is not important, but the RMS ratio between On and Ref spectra is (the absolute number is unimportant, because it gets calibrated away in the (On-Ref/Ref) equation).
Compared to the On scan, a Ref position is integrated over 'refpos_duration', which means that noise in the Ref spectrum is $\sqrt{\texttt{refpos_duration} / \texttt{sampling_interval}}$ smaller.
To increase the noise estimate accuracy, we will produce several noise maps at once.
End of explanation
dummy_tsys = 1.  
Explanation: Here, a dummy $T_\mathrm{sys}$ value is used. Later, the map (and thus the RMS) can simply be scaled to match the true $T_\mathrm{sys}$.
End of explanation
on_noise = dummy_tsys / np.sqrt(spec_reso * smooth_nbin * sampling_interval)
ref_noise = dummy_tsys / np.sqrt(spec_reso * smooth_nbin * refpos_duration)
print('on_noise = {:.2e}, ref_noise = {:.2e}'.format(on_noise, ref_noise))
reduced_specs = np.empty((num_scan_lines, samples_per_scanline, num_maps))
xcoords, ycoords = np.empty((2, num_scan_lines, samples_per_scanline))
lons = np.linspace(-map_width / 2, map_width / 2, samples_per_scanline)
lats = np.linspace(-map_height / 2, map_height / 2, num_scan_lines)
for scan_line in range(num_scan_lines):
    
    if scan_line % refpos_interval == 0:
        ref_spec = np.random.normal(0., ref_noise, num_maps) + dummy_tsys
    
    on_specs = np.random.normal(0., on_noise, (samples_per_scanline, num_maps)) + dummy_tsys    
    reduced_specs[scan_line] = dummy_tsys * (on_specs - ref_spec) / ref_spec
    
    xcoords[scan_line] = lons
    ycoords[scan_line] = np.repeat(lats[scan_line], lons.size)
   
Explanation: We will now create spectral noise (in arbitrary units). This must be much smaller than $T_\mathrm{sys}$ to avoid numerical problems. It is also important, that each of the raw-data spectra has an offset (the $T_\mathrm{sys}$ level in counts, if you will). If one would neglect this, one had a division by very small numbers in the line where the reduced spectra are calculated. The magnitude of the noise is furthermore tied to the $T_\mathrm{sys}$ level - it is based on the "radiometer" factor $\left(1/\sqrt{\tau\Delta f}\right)$
End of explanation
tmp_map = reduced_specs[..., 0]
cabs = np.max(np.abs(tmp_map))
fig = pl.figure(figsize=(10, 5))
ax = fig.add_axes((0.1, 0.1, 0.8, 0.8))
_ = ax.scatter(
    xcoords, ycoords, c=tmp_map, 
    cmap='bwr', edgecolor='none',
    vmin=-cabs, vmax=cabs,
    )
Explanation: We can test this, by just plotting the "raw" data for one channel.
End of explanation
target_header = setup_header((0, 0), (map_width / 3600., map_height / 3600.), beamsize_fwhm / 3600.)
Explanation: Now do the real gridding. First, prepare a WCS header
End of explanation
target_wcs = WCS(target_header)
# print(target_header)
Explanation: We already define a WCS object for later use in our plots:
End of explanation
gridder = cygrid.WcsGrid(target_header, naxis3=num_maps)
kernelsize_fwhm = beamsize_fwhm / 2  
kernelsize_fwhm /= 3600.  # need to convert to degree
kernelsize_sigma = kernelsize_fwhm / np.sqrt(8 * np.log(2))
support_radius = 4. * kernelsize_sigma
healpix_reso = kernelsize_sigma / 2.
gridder.set_kernel(
    'gauss1d',
    (kernelsize_sigma,),
    support_radius,
    healpix_reso,
    )
Explanation: Setup the gridder and define kernel sizes (half the beamsize is always a good choice).
End of explanation
gridder.grid(xcoords.flatten() / 3600, ycoords.flatten() / 3600, reduced_specs.reshape((-1, num_maps)))
cygrid_cube = gridder.get_datacube()
Explanation: The gridder needs the coordinates as flat arrays. The data to be gridded must be a 2D array (first dimension has to match the number of coordinate samples, second dimension is the number of channels/maps in the desired data cube).
End of explanation
tmp_map = cygrid_cube[0]
cabs = np.max(np.abs(tmp_map))
fig = pl.figure(figsize=(10, 10))
ax = fig.add_subplot(111, projection=target_wcs.celestial)
ax.imshow(tmp_map, cmap='bwr', interpolation='nearest', origin='lower', vmin=-cabs, vmax=cabs)
lon, lat = ax.coords
lon.set_axislabel('R.A. [deg]')
lat.set_axislabel('Dec [deg]')
pl.show()
Explanation: Again, as a sanity check, we plot one of the channel maps.
End of explanation
rms_cube = np.std(cygrid_cube, ddof=1)
Explanation: Last but not least, we can measure the noise. There are two possibilities to do this.
First, we can calculate the RMS over the full data cube.
End of explanation
rms_plane = np.mean(np.std(cygrid_cube, ddof=1, axis=(1, 2)))
print('rms_cube', rms_cube, 'rms_plane', rms_plane)
Explanation: Second, we can calculate the RMS per plane (and the average over all planes)
End of explanation
elevations = np.array([10, 20, 30, 40, 50, 60, 90])
AM = 1. / np.sin(np.radians(elevations))
gain_correction = gaincurve(elevations, 0.954, 3.19E-3, -5.42E-5)
Tsys_corr = Tsys_zenith + T_atm * (np.exp(opacity * AM) - np.exp(opacity * 1))
print('Tsys_corr', Tsys_corr)
Explanation: The ultimate questions is, what is the noise level in such a map, if we account for real $T_\mathrm{sys}$ and atmospheric effects, etc.
End of explanation
Ta_rms = rms_cube * Tsys_corr
if dual_pol:
    Ta_rms /= np.sqrt(2.)
Explanation: $T_\mathrm{A}$ noise is now the measured noise from the gridded maps multiplied with the $T_\mathrm{sys}$
End of explanation
Tb_rms = Ta_to_Tb * Ta_rms / gain_correction
S_rms = Tb_rms / Gamma_MB
atm_atten = np.exp(-opacity * AM)
print('-' * 95)
print('RMS per map')
print('-' * 95)
print('{0:>8s} {1:>8s} {2:>10s} {3:>10s} {4:>10s} {5:>10s} {6:>10s} {7:>10s} {8:>10s}'.format(
    'Elev', 'Airmass', 'Tsys', 'Ta RMS', 'Tb RMS', 'S RMS', 'AtmAtten', 'Tb_eff RMS', 'S_eff RMS'
    ))
print('{0:>8s} {1:>8s} {2:>10s} {3:>10s} {4:>10s} {5:>10s} {6:>10s} {7:>10s} {8:>10s}'.format(
    '[d]', '', '[K]', '[K]', '[K]', '[Jy]', '', '[K]', '[Jy]'
    ))
for idx in range(len(elevations)):
    print(
        '{0:>8.2f} {1:>8.2f} {2:>10.4f} {3:>10.4f} {4:>10.4f} '
        '{5:>10.4f} {6:>10.4f} {7:>10.4f} {8:>10.4f}'.format(
            elevations[idx], AM[idx], Tsys_corr[idx],
            Ta_rms[idx], Tb_rms[idx], S_rms[idx],
            atm_atten[idx],
            Tb_rms[idx] / atm_atten[idx], S_rms[idx] / atm_atten[idx],
            ))
print('Ta RMS = Antenna temp. noise')
print('Tb RMS = Brightness temp. noise')
print('S RMS = Flux density noise')
Explanation: Note, we don't need to account for the position switch again (division by noisy reference spectrum was already performed in the reduced-spectra array creation)!
End of explanation
<END_TASK>