Calculating Excitation Temperatures from 12CO cube¶
Ignoring self absorption and assuming infinite optical depth and Local Thermodynamic Equilibrium, it is possible to estimate the excitation temperature of a line from:
$$T_{B} =\frac{h\nu}{k_B} \left( \frac{1}{\exp[h\nu/k_{B}T_{ex}]} - \frac{1}{\exp[h\nu/k_{B}T_{CMB}]}\right)$$
Where $T_{B}$ is the RJ brightness temperature of the line, $\nu$ is the frequency of the line, $T_{CMB}$ is the CMB temperature, and $T_{ex}$ is the excitation temperature of the line.
We can use the peak temperature measured in our cube as the RJ brightness temperature $T_{B}$. For $^{12}$CO J=3$\rightarrow$2 as observed with HARP, we can then simplify and rearrange this formula to:
$$T_{ex} = \frac{16.59\,K}{ln\left(1 + \frac{16.59\,K}{T_{peak} + 0.036\,K}\right)}$$
To calculate this, we need a 12CO cube, which we will then get the peak temperature of each spectra in it. We then need to apply the above formula to this. This tutorial will take you through the steps of doing this within Starlink.
First of all, we must set up Starlink, change into our working directory, and initialise the KAPPA package of commands. We will also tell Starlink not to prompt us for user input, as we are working inside a jupyter notebook. (You should replace the /star with the correct location of your Starlink installation directory.)
export STARLINK_DIR=/star
source $STARLINK_DIR/etc/profile
export ADAM_NOPROMPT=1
mkdir CubeAnalysis
cd CubeAnalysis
kappa
We can then collapse our 12CO cube along the velocity axis. We want to end up with a 2-D cube where the pixel values are the max pixel value found in that spectra. We can use the command collapse with the estimator=max option.
collapse in=../Data/g34-3-12co_small_trim.sdf axis=VRAD estimator=max out=g34-3-12co_peakcoll.sdf WLIM=0.0
Lets take a quick look at what this map looks like, using KAPPA's display command, and setting some defaults appropriate for PNG output.
lutwarm dev=peak.png/PNG
palentry 0 White device=PNG
palentry 1 Black device=PNG
echo 'border=1' > style.dat
echo 'color(border) = black' >>style.dat
echo 'color(numlab) = black' >>style.dat
echo 'color(ticks) = grey' >>style.dat
display in=g34-3-12co_peakcoll.sdf style=^style.dat dev=peak.png/PNG mode=faint
Before we calculate our excitation temperature, lets check what temperature scale our data is in. We can check the label attribute of the NDF, as shown by ndftrace.
ndftrace g34-3-12co_peakcoll.sdf
You can see it is $T^*_A$ temperature scale, in units of $K$.
We now need to apply a standard calibration to go from the $T^*_A$ (corrected antenna temperature) scale to a telescope-independent scale. We will use $\eta_{MB}$, which allows us to convert to main beam temperature, which assumes a source that fills the beam. For HARP, we use the standard value of $\eta_{MB}=0.67$. We need to divide by the efficiency (cdiv) and set the new temperature scale into the label attribute of the NDF.
cdiv in=g34-3-12co_peakcoll.sdf scalar=0.67 out=12co_peakcoll_tmb.sdf
setlabel ndf=12co_peakcoll_tmb.sdf label='"T_MB Main Beam temperature"'
You can check that this has updated the label by rerunning ndftrace:
ndftrace 12co_peakcoll_tmb.sdf |grep Label
In practice, we probably don't want to include points with a poor detection. We can use errclip to set as BAD all pixels without a SNR of 5.
errclip in=12co_peakcoll_tmb.sdf out=12co_peakcoll_tmb_clip limit=5.0 mode=SNR
We can now calculate our excitation temperature in each pixel. The maths command lets us apply an arbitrary mathematical formula to our map. To apply the formula above, we use:
maths exp="pa/(log(1+pa/(ia+pb)))" pa=16.59 pb=0.036 ia=12co_peakcoll_tmb_clip.sdf out=12co_excitation_temp.sdf
The exp (expression) argument defines the mathemtaical function (using fortran style syntax). The parameters pa and pb refer to constants, and ia refers to an input NDF. These are then given as separate argument to the command. This is a very powerful command, that can take up to 26 separate constant parameters and up to 26 (ia, ib ... iz) input NDFs.
We can now display the excitation temperature map we have made:
palentry 0 White dev=PNG
palentry 1 Black dev=PNG
lutheat dev=PNG
display in=12co_excitation_temp.sdf style=^style.dat device=Tex.png/PNG mode=faint
We can then calculate the histogram of these data points using the histogram command, with 100 bins and using the same style.dat file as we used before. This command will create a plot of the histogram, and also print information to screen about the bins and the number of pixels in each one.
histogram 12co_excitation_temp 100 style=^style.dat accept dev=Tex_histogram.png/PNG
