One common task JCMT scientists carry out is estimating the dust mass or column density from the emission in a SCUBA-2 map. To do this, it is necessary to calculate the 850um dust flux of an object/field.
In general there are numerous scientific subtleties to which type of photometry is most appropriate for a given situation. Point sources are probably the most straight forward, followed by cases where you can integrate over an entire large area. The common (for galactic scientists) case of having extended structure throughout a large map, and hoping to break it down into separate structures lends itself to a clumpfinding approach, a filament finding approach or attempting to convert your map into hierarchical structures.
This guide will just quickly demonstrate some of the simpler cases and how to use Starlink software to calculate these. It is not a guide to which approach is appropriate for which scientific investigations, and it does not address the underlying uncertainties in these approaches.
First of all, we have to set up Starlink and KAPPA as normal. We will also create a working directory and change into it, and ensure that Starlink won't try and prompt us for any values.
(You should replace the /star
in the first line below with the location of your Starlink installation directory: e.g. ~/star-2018A
)
export STARLINK_DIR=/star
source $STARLINK_DIR/etc/profile
export ADAM_NOPROMPT=1
kappa
mkdir S2Fluxes
cd S2Fluxes/
SCUBA-2 maps are normally calibrated by our pipeline into units of either mJy/beam or mJy/arcsec$^2$. If you need more detail on the calibration process, please see Dempsey et al 2013, or our webpages at https://www.eaobservatory.org/jcmt/instrumentation/continuum/scuba-2/calibration.
If you want to find the flux of a point source, you would normally calibrate your map in units of Jy/beam. Each point source corresponds to one beam, and therefore the peak value in the source will tell you the total flux density of the source.
For this example, we are going to use a calibrator observation of the bright point source CRL618, in the file Data/jcmts20191022_00059_850_reduced001_nit_000.sdf
. You can download this file from CADC at this link here, or you can go to the JCMT Archive at CADC and download a reduced calibrator observation of your choice. If you download a file from the archive, you'll need to convert it from FITS format to NDF (see fits2ndf
in the CONVERT
package), and please ensure it is a point source.
The first step is to see what units your map is currently in. We can use the KAPPA command ndftrace
to do this. (Here we have also used the shell command head
to only show the first 10 lines).
ndftrace ../Data/jcmts20191022_00059_850_reduced001_nit_000.sdf|head -n 10
As you can see, this is already in mJy/beam, so we don't need to change the calibration. If we did, we would normally use the PICARD recipe CALIBRATE_SCUBA2_DATA
with the recipe parameters FCF_CALTYPE=BEAM,USEFCF=True
.
The nature of JCMT scan patterns mean that we normally have areas of higher noise around the map. There is also a PICARD recipe that can crop that out, CROP_SCUBA2_IMAGES
. Here we will call it asking for a circle of radius 250". (Please note this requires all the input observations to a map to have been towards the same source, or it will be hard to predict where the mask will be centered.)
picard -log sf -nodisplay CROP_SCUBA2_IMAGES --recpars='CROP_METHOD=CIRCLE,MAP_RADIUS=250' \
../Data/jcmts20191022_00059_850_reduced001_nit_000.sdf
The output file is jcmts20191022_00059_850_reduced001_nit_000_crop
. We can now take a quick look at our map to make sure it looks sensible.
lutwarm dev=PNG
palentry 0 white dev=PNG
palentry 1 black dev=PNG
display in=jcmts20191022_00059_850_reduced001_nit_000_crop dev=cropped.png/PNG mode=faint style=!
You can get a quick and crude estimate of the flux of the point source by just "reading off" the peak flux in the source, using e.g. KAPPA's stats
. The brightest pixel should roughly give you the flux of the source in mJy, assuming there is not a high background level.
stats jcmts20191022_00059_850_reduced001_nit_000_crop
The peak pixel has a value of 5556.168 mJy/beam, corresponding to a source flux of 5.6 Jy. This is quite a bit higher than our website value of 4.89 ± 0.24 Jy/beam for the peak brightness (see https://www.eaobservatory.org/jcmt/instrumentation/continuum/scuba-2/calibration/calibrators/, but within the correct range.
A better estimate would be to fit a 2-D Gaussian to the source and find the peak. We can do this with the KAPPA command beamfit
. We take the initial coordinates to be the position of the max pixel, as found from stats
above.
beamfit ndf=jcmts20191022_00059_850_reduced001_nit_000_crop mode=interface pos='"4:42:53.7,36:06:53"' beams=1
From the amplitude value, this gives a flux of 5.54+/-0.02 Jy for the source.
(You can see from the output that this has a background of 1.2mJy/beam, or or 0.001 Jy, so any residual background should not be a significant effect on this value).
If you want to calculate the flux integrated over resolved sources you have to consider what you mean by a source. One approach for a compact object would be to use annular photometry, including a correction for the background. This is what we do to calculate our arcsecond FCF (see Dempsey et al 2013). Another approach would be to calculate the flux within an entire map or object. Here we will demonstrate calculating the integrated flux in a 2-D Clump, such as that found when running the CUPID fellwalker on a map.
Here we will use a map of G34.3, an occassional SCUBA-2 calibrator source. The file is in Data/jcmts20191018_00009_850_reduced001_nit_000.sdf
First we need to ensure our map is in Jy/arcsec$^{2}$.
ndftrace ../Data/jcmts20191018_00009_850_reduced001_nit_000.sdf
You can see that it is in mJy/beam, from the 'Unit' component shown by ndftrace
. We will use the PICARD recipe CALIBRATE_SCUBA2_DATA
to switch the units to mJy/arced$^2$. (This recipe is very simple: it checks what the current calibration is (if any), undoes it if required by dividing by the FCF, then multiplies by the correct standard FCF factor for the requested units.)
picard -log sf -nodisplay CALIBRATE_SCUBA2_DATA --recpars='FCF_CALTYPE=ARCSEC,USEFCF=True' \
../Data/jcmts20191018_00009_850_reduced001_nit_000.sdf
The output file name is jcmts20191022_00059_850_reduced001_nit_000_uncal_cal.sdf
: this shows you that PICARD first removed the existing calibration _uncal
and then applied the new calibration _cal
.
We'll now run a very quick clumpfinding routine using the Fellwalker algorithm in CUPID. First we have to initialise the CUPID package:
cupid
Next we will run the findclumps
command with the fellwalker algorithm, but leaving most variables to their default value.
findclumps in=jcmts20191018_00009_850_reduced001_nit_000_uncal_cal method=FellWalker out=clumps outcat=! accept
You can see from the output that it found 8 clumps. We can take a quick look at the data values in the first clump by examining the CUPID.CLUMPS extension within the NDF like so:
display in=clumps.sdf.MORE.CUPID.CLUMPS\(1\) dev=clumps.png/PNG mode=faint style=!
stats clumps.sdf.MORE.CUPID.CLUMPS\(1\)
As you can see from when we ran ndftrace
above, we have 4" pixels. Therefore we need to integrate our flux over the whole region, to get a value in Jy instead of in Jy/arcsec$^{2}$. The flux in our clump is therefore the total sum of the emission, which is 8828 mJy/arcsec$^{2}$from stats, multipled by 16 square arcseconds = 141 Jy. This assumes that the background level is 0.