NGC3256Band3 for CASA 3.3: Difference between revisions
| Line 820: | Line 820: | ||
</source> | </source> | ||
Make images of the CN high emission | |||
<source lang="python"> | <source lang="python"> | ||
# In CASA | # In CASA | ||
Revision as of 15:43, 20 May 2011
Overview
[To be written by Eric]
Retrieving the Data
The data were taken in six different datasets over two consecutive nights: April 16-17, 2011. There are three datasets for April 16th and three for April 17th. Here we provide you with "starter" datasets, where we have taken the raw data in ALMA Science Data Model (ASDM) format and converted them to CASA Measurement Sets (MS). We did this using the importasdm task in CASA.
[What else are we going to do to the data we provide?]
Along with the Measurement Sets, we also provide some tables that you will need for the calibration. These include the System Temperature (Tsys) tables, which contain corrections for atmospheric opacity, and Water Vapor Radiometer (WVR) tables, which contain the atmospheric phase corrections determined by the water vapor radiometers on each antenna.
You can download the data here: [Provide link to the raw .ms files in tar'd, gzip'd format]
Once the download has finished, unpack the file:
# In a terminal outside CASA
tar -xvf ngc3256band3.tgz
[Also provide links to the calibrated data (but maybe not here...maybe better at end of calibration page?)]
Initial Inspection and A priori Flagging
We will eventually concatenate the six datasets used here into one large dataset. However, we will keep them separate for now, as some of the steps to follow require individual datasets (specifically, the application of the Tsys and WVR tables). We therefore start by defining an array 'basename' that includes the names of the six files. This will simplify the following steps by allowing us to loop through the files using a simple for-loop in python.
# In CASA
basename=['uid___A002_X1d54a1_X5','uid___A002_X1d54a1_X174','uid___A002_X1d54a1_X2e3',
'uid___A002_X1d5a20_X5','uid___A002_X1d5a20_X174','uid___A002_X1d5a20_X330']
The usual first step is then to get some basic information about the data. We do this using the task listobs, which will output a detailed summary of each dataset supplied.
# In CASA
for name in basename:
listobs(vis=name+'.ms')
The output will be sent to the CASA logger. You will have to scroll up to see the individual output for each of the six datasets. Here is an example of the most relevant output for the first file in the list.
Fields: 3 ID Code Name RA Decl Epoch SrcId nVis 0 none 1037-295 10:37:16.0790 -29.34.02.8130 J2000 0 38759 1 none Titan 00:00:00.0000 +00.00.00.0000 J2000 1 16016 2 none NGC3256 10:27:51.6000 -43.54.18.0000 J2000 2 151249 (nVis = Total number of time/baseline visibilities per field) Spectral Windows: (9 unique spectral windows and 2 unique polarization setups) SpwID #Chans Frame Ch1(MHz) ChanWid(kHz)TotBW(kHz) Ref(MHz) Corrs 0 4 TOPO 184550 1500000 7500000 183300 I 1 128 TOPO 113211.988 15625 2000000 113204.175 XX YY 2 1 TOPO 114188.55 1796875 1796875 113204.175 XX YY 3 128 TOPO 111450.813 15625 2000000 111443 XX YY 4 1 TOPO 112427.375 1796875 1796875 111443 XX YY 5 128 TOPO 101506.187 15625 2000000 101514 XX YY 6 1 TOPO 100498.375 1796875 1796875 101514 XX YY 7 128 TOPO 103050.863 15625 2000000 103058.675 XX YY 8 1 TOPO 102043.05 1796875 1796875 103058.675 XX YY Sources: 48 ID Name SpwId RestFreq(MHz) SysVel(km/s) 0 1037-295 0 - - 0 1037-295 9 - - 0 1037-295 10 - - 0 1037-295 11 - - 0 1037-295 12 - - 0 1037-295 13 - - 0 1037-295 14 - - 0 1037-295 15 - - 0 1037-295 1 - - 0 1037-295 2 - - 0 1037-295 3 - - 0 1037-295 4 - - 0 1037-295 5 - - 0 1037-295 6 - - 0 1037-295 7 - - 0 1037-295 8 - - 1 Titan 0 - - 1 Titan 9 - - 1 Titan 10 - - 1 Titan 11 - - 1 Titan 12 - - 1 Titan 13 - - 1 Titan 14 - - 1 Titan 15 - - 1 Titan 1 - - 1 Titan 2 - - 1 Titan 3 - - 1 Titan 4 - - 1 Titan 5 - - 1 Titan 6 - - 1 Titan 7 - - 1 Titan 8 - - 2 NGC3256 0 - - 2 NGC3256 9 - - 2 NGC3256 10 - - 2 NGC3256 11 - - 2 NGC3256 12 - - 2 NGC3256 13 - - 2 NGC3256 14 - - 2 NGC3256 15 - - 2 NGC3256 1 - - 2 NGC3256 2 - - 2 NGC3256 3 - - 2 NGC3256 4 - - 2 NGC3256 5 - - 2 NGC3256 6 - - 2 NGC3256 7 - - 2 NGC3256 8 - - Antennas: 7: ID Name Station Diam. Long. Lat. 0 DV04 J505 12.0 m -067.45.18.0 -22.53.22.8 1 DV06 T704 12.0 m -067.45.16.2 -22.53.22.1 2 DV07 J510 12.0 m -067.45.17.8 -22.53.23.5 3 DV08 T703 12.0 m -067.45.16.2 -22.53.23.9 4 DV09 N602 12.0 m -067.45.17.4 -22.53.22.3 5 PM02 T701 12.0 m -067.45.18.8 -22.53.22.2 6 PM03 J504 12.0 m -067.45.17.0 -22.53.23.0
This output shows that three fields were observed: 1037-295, Titan, and NGC3256. Field 0 (1037-295) will serve as the gain calibrator and bandpass calibrator; field 1 (Titan) will serve as the flux calibrator; and field 2 (NGC3256) is, of course, the science target.
Note that there are more than four SpwIDs even though the observations were set up to have four spectral windows. The spectral line data themselves are found in spectral windows 1,3,5,7, which have 128 channels each. The first one (spw 1) is centered on the CO(1-0) emission line in the galaxy NGC 3256 and is our highest frequency spectral window. There is one additional spectral window (spw 3) in the Upper Side Band (USB), and there are two spectral windows (spw 5 and 7) in the Lower Side Band (LSB). These additional spectral windows are used to measure the continuum emission in the galaxy, and may contain other emission lines as well.
Spectral windows 2,4,6,8 contain channel averages of the data in spectral windows 1,3,5,7, respectively. These are not useful for the offline data reduction. Spectral window 0 contains the WVR data. You may notice that there are additional SpwIDs listed in the "Sources" section which are not listed in the "Spectral Windows" section. These spectral windows are reserved for the WVRs of each antenna (seven in our case). At the moment, all WVRs point to spw 0, which contains nominal frequencies. The additional spectral windows (spw 9-15) are therefore not used and can be ignored.
Another important thing to note is that the position of Titan is listed as 00:00:00.0000 +00.00.00.0000. This is due to the fact that for ephemeris objects, the positions are currently not stored in the asdm. This will be handled correctly in the near future, but at present, we have to fix this offline. We will correct the coordinates below by running the procedure fixplanet.
There were seven antennas in the array for these observations. Note that numbering in python always begins with "0", so the antennas have IDs 0-6. To see what the antenna configuration looked like at the time of the observations, we will use the task plotants. Since the configuration did not change during the course of the observations, we will simply look at the first dataset from the list. [ACTUALLY IT DID. THE SECOND DAY HAS 8 ANTENNAS. WE CAN JUST MENTION THAT HERE]
# In CASA
plotants(vis=basename[0]+'.ms', figfile='ngc3256band3_plotants.png')
This will plot the antenna configuration on your screen as well as save it under the specified filename for future reference. This will be important later on when we need to choose a reference antenna, since the reference antenna should be close to the center of the array (as well as stable and present for the entire observation).
The first editing we will do is some a priori flagging. We will start by flagging the shadowed data and the autocorrelation data:
# In CASA
for name in basename:
flagdata(vis=name+'.ms', flagbackup = F, mode = 'shadow')
flagautocorr(vis=name+'.ms')
There are a number of scans in the data that were used by the online system for pointing calibration. These scans are no longer needed, and we can flag them easily by selecting on 'intent':
# In CASA
for name in basename:
flagdata(vis=name+'.ms', mode='manualflag', flagbackup = F, intent='*POINTING*')
Similarly, we can flag the scans corresponding to atmospheric calibration:
# In CASA
for name in basename:
flagdata(vis=name+'.ms', mode='manualflag', flagbackup = F, intent='*ATMOSPHERE*')
We will then store the current flagging state for each dataset using the flagmanager:
# In CASA
for name in basename:
flagmanager(vis = name+'.ms', mode = 'save', versionname = 'Apriori')
We will continue with some initial flagging/corrections specific to these datasets. For uid___A002_X1d54a1_X174.ms there is a outlying feature in spw=7, antenna DV04. This corresponds to scans 5 and 9, so we flag those data:
# In CASA
flagdata(vis='uid___A002_X1d54a1_X174.ms', mode='manualflag',
antenna='DV04', flagbackup = F, scan='5,9', spw='7')
Antenna DV07 shows large delays for the first three datasets. We correct this by calculating a K-type delay calibration table with gencal. The parameters are the delays measured in nanoseconds, first cycling over polarization product, and then over spectral window (thus giving eight numbers in total). Before creating these tables, make sure to delete any existing versions.
# In CASA
for i in range(3): # loop over the first three ms's
name=basename[i]
os.system('rm -rf '+name+'_del.K')
gencal(vis=name+'.ms', caltable=name+'_del.K',
caltype='sbd', antenna='DV07', pol='X,Y', spw='1,3,5,7',
parameter=[0.99, 1.10, -3.0, -3.0, -3.05, -3.05, -3.05, -3.05])
The delay values were simply estimated by eye. The purpose here is to get the phases approximately flat as a function of frequency. Any additional phase variations will be corrected for later when we do the bandpass calibration.
[MARTIN, HOW DID YOU GET THESE PARAMETERS? I ASSUME YOU FIT THE WRAPPING SOMEHOW AND PERHAPS WE SHOULD MENTION IT - OK, I simply guestimated these, see above. Maybe we can phrase that better.]
[ALSO, WHEN I RUN THIS, CASA SAYS "No scr col generation!!!" DO YOU KNOW WHAT THIS MEANS? I SEE YOU NOTE IT LATER IN YOUR SCRIPT AS WELL - Yes, this means that the field association is not made in the table, that is, you can apply te tables t any of the fields in your data. This is not really important, so I think we can a) ignore it, or b) mention it here and tell people to ignore it. What do you think?]
We will apply these K tables to the data in the next section.
WVR Correction and Tsys Calibration
First, we apply the delay correction table and the WVR calibration tables to the data. We do this in two steps, first cycling over the three datasets from the first day of observations because we have to correct the delay error for DV07 for those data. For the last three datasets, taken during the second day, we do not need to correct the delays, so we just apply the WVR tables. For both types of tables, we will use interpolation "nearest"... [BECAUSE WHY? - I think for the WVR data you don't want any linear interpolation, but just apply the nearest value. For the delay it does not really matter what you use because there is no time dependence ]
# In CASA
for i in range(3): # loop over the first three data sets
name=basename[i]
applycal(vis=name+'.ms', flagbackup=F, spw='1,3,5,7',
interp=['nearest','nearest'], gaintable=[name+'_del.K',name+'.W'])
for i in range(3,6): # loop over the last three data sets
name=basename[i]
applycal(vis=name+'.ms', flagbackup=F, spw='1,3,5,7',
interp='nearest', gaintable=name+'.W')
Now we split out the datasets with delays and WVR tables applied. These datasets are given the extention "_K_WVR" to indicate that the delay tables and WVR tables have been applied. Again, we are careful to remove any previous versions of the split ms's before running the split command.
# In CASA
for name in basename:
os.system('rm -rf '+name+'_K_WVR.ms*')
split(vis=name+'.ms', outputvis=name+'_K_WVR.ms',
datacolumn='corrected', spw='0~8')
[WHEN I RUN THIS, I GET A MESSAGE ABOUT SPW 9-15 MISSING, AND IT'S GOING TO IGNORE THEM. THIS IS FINE, OF COURSE, BUT DO WE NEED TO SAY ANYTHING ABOUT IT? _ Yes, this is fine. I think we should just mention that, or we can say spw='0~8' in the split (as above)]
Next we do the Tsys calibration. Tsys measurements correct for the atmospheric opacity (to first-order) and allow the calibration sources to be measured at elevations that differ from the science target. The Tsys tables for these datasets were provided with the downloadable data. We will start by inspecting them:
# In CASA
for name in basename:
plotcal(caltable='tsys_'+name+'.cal', xaxis='freq', yaxis='amp',
spw='1,3,5,7', timerange='<2020', subplot=221, overplot=False,
iteration='spw', plotrange=[0, 0, 40, 180], plotsymbol='.',
figfile='tsys_per_spw'+name+'.png')
Note that we only plot the spectral windows with the spectral line data, and we set timerange="<2020" because the SYSCAL table currently includes values from March 2081. [IS THERE ANY WAY TO GET AROUND EXPLAINING THIS EMBARRASSING ERROR HERE? -- I think this is fixed now with the new aU tools to fix the Tsys tables. So, we can take out the 'timerange<2020'] In addition to plotting on your screen, the above command will also produce a plot file (png) for each of the six datasets.
The plots look acceptable upon examination, so we will apply the Tsys tables with applycal: We do this for each field separately so that the appropriate calibration data are applied to the right fields. The "field" parameter specifies the field to which we will apply the calibration, and the "gainfield" parameter specifies the field from which you wish to take the calibration solutions from the gaintable.
# In CASA
for name in basename:
for field in ['Titan','1037*','NGC*']:
applycal(vis=name+'_K_WVR.ms', spw='1,3,5,7', flagbackup=F, field=field, gainfield=field,
interp='nearest', gaintable=['tsys_'+name+'.cal'])
[**I COULDNT GET THIS TO WORK: KEPT GETTING SEVERE ERROR: SEVERE applycal::Calibrater::setapply(type, applypar) Error loading 'tsys_uid___A002_X1d5a20_X330.cal' with parsed selection: (FIELD_ID IN [2]). 2011-05-20 01:11:47 SEVERE applycal::Calibrater::setapply(type, applypar)+ (/var/rpmbuild/BUILD/casapy-stable/casapy-stable-32.0.15107/code/include/calibration/CalTables/CalSet.tcc : 320) Failed AlwaysAssert numberAnt==nElem_ Check inputs and try again. 2011-05-20 01:11:47 SEVERE applycal:::: Exception Reported: Error in Calibrater::setapply.
IT SEEMS TO BE COMPLAINING ABOUT THE FIELD SPECIFICATION IN THE TSYS SCRIPT. I GENERATED THE TSYS TABLES THE WAY MARTIN DID IN THE SCRIPT...IS THAT NOT CORRECT? No, you need to fix them with the new tools that Stuartt wrote. See my email.]
We then split out spectral windows 1,3,5,7. This will get rid of the extraneous spectral windows, including the channel averaged spectral windows and spw 0, which is the one for the WVR data.
# In CASA
for name in basename:
os.system('rm -rf '+name+'_line.ms*')
split(vis=name+'_K_WVR.ms', outputvis=name+'_line.ms',
datacolumn='corrected', spw='1,3,5,7')
The WVR and Tsys tables are now applied in the DATA column of the resultant measurement sets. The new data sets have the extension "_line" to indicate that they only contain the line data, and no longer the "channel average" spectral windows. These measurement sets therefore have four spectral windows.
Now that we have applied the Tsys calibration and WVR corrections, we can concatenate the six individual data sets into one big measurement set. We define an array "comvis" that contains the names of the measurement sets we wish to concatenate, and then we run the task concat.
# In CASA
comvis=[]
for name in basename:
comvis.append(name+'_line.ms')
os.system('rm -rf ngc3256_line.ms*')
concat(vis=comvis, concatvis='ngc3256_line.ms')
Additional Data Inspection
Now that the data are concatenated into one dataset, we will do some additional inspection with plotms. First we will plot amplitude versus channel, averaging over time and baselines in order to speed up the plotting process.
# In CASA
plotms(vis='ngc3256_line.ms', xaxis='channel', yaxis='amp',
averagedata=T, avgbaseline=T, avgtime='1e8', avgscan=T)
From this plot we see that the edge channels have abnormally high amplitudes. We will use flagdata to remove some channels from both sides of the bandpass:
# In CASA
flagdata(vis = 'ngc3256_line.ms', flagbackup = F, spw = ['*:0~10','*:125~127'])
Next, we will look at amplitude versus time, averaging over all channels and colorizing by field. Since the observations take place over two days, you will need to zoom in to examine the data in more detail. In particular, note the difference in Titan's amplitude between the two days and the change in amplitude during the second day.
# In CASA
plotms(vis='ngc3256_line.ms', xaxis='time', yaxis='amp',
averagedata=T, avgchannel='128', coloraxis='field')
Titan is our primary flux calibrator. However, for the second day of observations, Titan had moved to close to Saturn, and Saturn's rings moved into the primary beam. Another way to see this is to plot amplitude versus uv-distance and colorize by scan: [MARTIN, IS THIS WHAT YOU HAD IN MIND FOR THIS PLOT? - yes!]
# In CASA
plotms(vis='ngc3256_line.ms', xaxis='uvdist', yaxis='amp', field='1',
averagedata=T, avgchannel='128', avgtime='1e8', coloraxis='scan')
If you use the "Mark a region" and "Locate" tools in plotms, you will see the signature of a bright, resolved source during the latter three (last day's) scans. We will therefore need to flag the Titan scans for the second day:
# In CASA
flagdata(vis = 'ngc3256_line.ms', flagbackup = F,
timerange='>2011/04/16/12:00:00', field='Titan')
Next, we will fix the position of Titan in the combined dataset. Recall that the position of the Titan field is currently set to 00:00:00.0000 +00.00.00.0000. The following procedure will replace this with the actual position observed by the telescopes and, at the same time, it will recalculate the uvw coordinates:
# In CASA
execfile(os.getenv("CASAPATH").split(' ')[0]+"/lib/python2.6/recipes/fixplanets.py")
fixplanets('ngc3256_line.ms', 'Titan', True)
Continue to inspect the data with plotms, plotting different axes and colorizing by the different parameters. Don't forget to average the data if possible to speed the plotting process. You will find the following:
- Baselines with DV07 have very high amplitudes in spw 3, correlation YY
- Baselines with DV08 have very low amplitudes in spw 3, correlation YY, but only for the last observation
- Baselines with PM03 have low amplitudes at 2011/04/17/02:15:00 for spw 0
- Baselines with PM03 have low amplitudes at 2011/04/16/04:15:15 for spw 2 and 3
We flag the bad data with the following commands:
# In CASA
flagdata(vis='ngc3256_line.ms', flagbackup=F, spw='3',
correlation='YY', mode='manualflag', selectdata=T,
antenna='DV07', timerange='')
flagdata(vis='ngc3256_line.ms', flagbackup=F, spw='3',
correlation='YY', mode='manualflag', selectdata=T,
antenna='DV08', timerange='>2011/04/17/03:00:00')
flagdata(vis='ngc3256_line.ms', flagbackup=F, spw='0',
mode='manualflag', selectdata=T, antenna='PM03',
timerange='2011/04/17/02:15:00~02:15:50')
flagdata(vis='ngc3256_line.ms', flagbackup=F, spw='2,3',
mode='manualflag', selectdata=T, antenna='PM03',
timerange='2011/04/16/04:13:50~04:18:00')
Bandpass Calibration
We are now ready to begin the bandpass calibration. The first thing we will do is run the task gaincal on the bandpass calibrator to determine phase-only gain solutions. We will use solint='int' for the solution interval, which means that one gain solution will be determined for every integration time. This will correct for any phase variations in the bandpass calibrator as a function of time, a step which will prevent decorrelation of the vector-averaged bandpass solutions. We will then apply these solutions on-the-fly when we run bandpass.
We will use the average of channels 40 to 80 to increase our signal-to-noise in the determination of the antenna-based phase solutions. Averaging over a subset of channels near the center of the bandpass is acceptable when the phase variation as a function of channel is small, which it is here. For our reference antenna, we choose PM03. We call the output calibration table "ngc3256.G1".
# In CASA
gaincal(vis = 'ngc3256_line.ms', caltable = 'ngc3256.G1', spw = '*:40~80', field = '1037*',
selectdata=T, solint= 'int', refant = 'PM03', calmode = 'p')
We check the time variations of the phase solutions with plotcal. We will plot the XX and YY polarization products separately and make different subplots for each of the spectral windows. This is done by setting the 'iteration' parameter to 'spw' and specifying subplot=221. By setting the parameter "figfile" to a non-blank value, it will also generate png files of the plots.
# In CASA
plotcal(caltable = 'ngc3256.G1', xaxis = 'time', yaxis = 'phase',
poln='X', plotsymbol='o', plotrange = [0,0,-180,180], iteration = 'spw',
figfile='phase_vs_time_XX.G1.png', subplot = 221)
plotcal(caltable = 'ngc3256.G1', xaxis = 'time', yaxis = 'phase',
poln='Y', plotsymbol='o', plotrange = [0,0,-180,180], iteration = 'spw',
figfile='phase_vs_time_YY.G1.png', subplot = 221)
Now that we have a first measurement of the phase variations as a function of time, we can determine the bandpass solutions with bandpass, applying the phase calibration table on-the-fly.
First, we will plot the phase and amplitude as a function of frequency for the bandpass calibrator. We use avgscan=T and avgtime='1E6' to average in time over all scans, and we specify coloraxis='baseline' to colorize by baseline. [WHY ARE WE DOING THIS? JUST TO SEE WHAT IT LOOKS LIKE BEFORE THE BANDPASS CALIBRATION? Yes, just a basic inspection. Happy to skip this if you think it's not useful.]
# In CASA
plotms(vis='ngc3256_line.ms', xaxis='freq', yaxis='phase', selectdata=True,
field='1037*', avgtime='1E6', avgscan=T, coloraxis='baseline', iteraxis='antenna')
plotms(vis='ngc3256_line.ms', xaxis='freq', yaxis='amp', selectdata=True, spw='*:10~120',
field='1037*', avgtime='1E6', avgscan=T, coloraxis='baseline', iteraxis='antenna')
Now we do the bandpass calibration, applying the phase calibration solutions on-the-fly with the parameter "gaintable". We will use solint='inf' to determine a solution for each scan separately. [...Need to discuss parameters in more detail]
# In CASA
bandpass(vis = 'ngc3256_line.ms', caltable = 'ngc3256.B1', gaintable = 'ngc3256.G1',
field = '1037*', minblperant=3, minsnr=1, solint='inf',
bandtype='B', fillgaps=1, refant = 'PM03', solnorm = F)
We then plot the bandpass solutions with the following commands:
# In CASA
plotcal(caltable = 'ngc3256.B1', xaxis='freq',yaxis='phase', spw='',
subplot=212, overplot=False, plotrange = [0,0,-70,70],
plotsymbol='.', timerange='')
plotcal(caltable = 'ngc3256.B1', xaxis='freq',yaxis='amp', spw='',
subplot=211, overplot=False,
figfile='bandpass.B1.png', plotsymbol='.', timerange='')
Gain Calibration
First, get the flux density for Titan using the Butler-JPL-Horizons 2010 model. The flux density of Titan is 296 mJy
# In CASA
setjy(vis='ngc3256_line.ms', field='Titan', spw='', modimage='',
scalebychan=False, fluxdensity=-1,
standard='Butler-JPL-Horizons 2010')
Now we will do a new gain calibration, this time applying the bandpass calibration solutions on-the-fly:
# In CASA
gaincal(vis = 'ngc3256_line.ms', caltable = 'ngc3256.G2', spw =
'*:30~90', field = '1037*,Titan',
solint= 'inf', selectdata=T, solnorm=False, refant = 'PM03',
gaintable = ['ngc3256.B1'], calmode = 'ap')
[ YOU ARE JUST SOLVING FOR A&P AT THE SAME TIME? YOU DON'T SOLVE FOR P FIRST USING A SHORTER TIME INTERVAL, THEN USE THAT TO DERIVE THE AMPLITUDE SOLUTIONS PER SCAN? MAYBE IT DOESN'T MAKE A BIG DIFFERENCE HERE? - thought it would not make much of a difference, but please try if you have time. Otherwise just leave as is.]
Now we will examine the amplitude and phase solutions as a function of time, iterating on spectral window and plotting the XX and YY correlations separately for clarity:
# In CASA
plotcal(caltable = 'ngc3256.G2', xaxis = 'time', yaxis = 'phase',
poln='X', plotsymbol='o', plotrange = [0,0,-180,180], iteration
= 'spw', figfile='phase_vs_time_XX.G2.png', subplot = 221)
plotcal(caltable = 'ngc3256.G2', xaxis = 'time', yaxis = 'phase',
poln='Y', plotsymbol='o', plotrange = [0,0,-180,180], iteration
= 'spw', figfile='phase_vs_time_YY.G2.png', subplot = 221)
plotcal(caltable = 'ngc3256.G2', xaxis = 'time', yaxis = 'amp',
poln='X', plotsymbol='o', plotrange = [], iteration = 'spw',
figfile='amp_vs_time_XX.G2.png', subplot = 221)
plotcal(caltable = 'ngc3256.G2', xaxis = 'time', yaxis = 'amp',
poln='Y', plotsymbol='o', plotrange = [], iteration = 'spw',
figfile='amp_vs_time_YY.G2.png', subplot = 221)
Finally, we will bootstrap the fluxes of the secondary calibrator and science target to that of Titan using the task fluxscale:
# In CASA
fluxscale( vis='ngc3256_line.ms', caltable='ngc3256.G2',
fluxtable='ngc3256.G2.flux', reference='Titan',
transfer='1037*', append=False)
[WE WILL HAVE TO SAY SOMETHING ABOUT THE FACT THAT WE APPLY THE TITAN CAL TO BOTH DAYS...]
We find that the flux density of 1037-295 is 1.8 Jy, which agrees very well with the 3mm measurements from the ATCA and OVRO of a few years ago.
Applying the calibrations
apply the bandpass and gaincal tables. Apply the calibrations to the phase cal and science target
# In CASA
applycal( vis='ngc3256_line.ms', flagbackup=F, field='NGC*,1037*',
interp='nearest', gainfield = ['1037*', '1037*', '1037*'],
gaintable=['ngc3256.G2.flux','ngc3256.B1'])
Plot amplitude as function of freq for spw 0 and 1
# In CASA
plotms(vis="ngc3256_line.ms", xaxis="freq", yaxis="amp",
ydatacolumn="corrected", selectdata=True, field="NGC*",
spw="0,1", averagedata=True, avgchannel="", avgtime="10000s",
avgscan=True, avgbaseline=T)
There are three emission lines. There is no emission in spw 2 and 3
There are some additional bad data points
# In CASA
flagdata( vis='ngc3256_line.ms', mode='manualflag',
timerange='2011/04/16/04:13:35~04:13:45', flagbackup = F)
tget(flagdata)
timerange='2011/04/16/05:21:13~05:21:19'
flagdata()
timerange='2011/04/16/04:16:40~04:16:49'
flagdata()
Image the phase calibrator
# In CASA
os.system('rm -rf phasecal_cont*')
clean(vis="ngc3256_line.ms", imagename="phasecal_cont",
field="1037*", spw='*:20~120',
selectdata=T, mode="mfs", niter=1000, gain=0.1,
psfmode="clark", interactive=False, mask=[126, 126, 130,130],
imsize=256, cell="1arcsec", weighting="briggs", robust=0.0)
Get some statistics on the phasecal image
# In CASA
calstat=imstat(imagename="phasecal_cont.image",region="",box="160,50,210,210")
rms=(calstat['rms'][0])
print ">> rms in phase calibrator image: "+str(rms)
calstat=imstat(imagename="phasecal_cont.image",region="")
peak=(calstat['max'][0])
print ">> Peak in phase calibrator image: "+str(peak)
print ">> Dynamic range in phase calibrator image: "+str(peak/rms)
The image dynamic range in calibrator image ~2500
Make a nice image of the phase cal:phasecal_map.png
# In CASA
imview(raster={'file': 'phasecal_cont.image', 'colorwedge':T,
'range':[-0.004, 0.250], 'scaling':-2.5, 'colormap':
'Rainbow 2'},
out='phasecal_map.png', zoom=1)
Applycal on Titan. Use the phase solution on Titan this time.
# In CASA
applycal(vis='ngc3256_line.ms', flagbackup=F, field='Titan',
interp='nearest', gainfield = ['Titan', '1037*'],
gaintable=['ngc3256.G2.flux','ngc3256.B1'])
Image the amplitude calibrator (Titan)
# In CASA
os.system('rm -rf ampcal_cont*')
clean( vis="ngc3256_line.ms", imagename="ampcal_cont", field="Titan",
spw="0:20~120,1:20~120,2:20~120,3:20~120", selectdata=T,
mode="mfs", niter=500, gain=0.1, psfmode="clark",
interactive=False, mask=[126, 126, 130,130], imsize=256,
cell="1arcsec", weighting="briggs", robust=0.0)
Get some statistics on the ampcal image
# In CASA
calstat=imstat(imagename="ampcal_cont.image",region="",box="160,50,210,210")
rms=(calstat['rms'][0])
print ">> rms in amp calibrator image: "+str(rms)
calstat=imstat(imagename="ampcal_cont.image",region="")
peak=(calstat['max'][0])
print ">> Peak in amp calibrator image: "+str(peak)
print ">> Dynamic range in amp calibrator image: "+str(peak/rms)
Image dynamic range in calibrator image ~64
# In CASA
imview(raster={'file': 'ampcal_cont.image', 'colorwedge':T,
'range':[-0.020, 0.250], 'scaling':-2.5, 'colormap':
'Rainbow 2'},
out='ampcal_map.png', zoom=1)
Continuum image of the galaxy
Split out the science target
# In CASA
os.system('rm -rf ngc3256_line_target.ms*')
split(vis='ngc3256_line.ms', outputvis='ngc3256_line_target.ms',
field='NGC*')
Clean the continuum image
# In CASA
os.system('rm -rf ngc3256.cont*')
clean( vis="ngc3256_line_target.ms", imagename="ngc3256_cont",
field='NGC*', spw="0:20~53;71~120,1:70~120,2:20~120,3:20~120",
selectdata=T, mode="mfs", niter=500, gain=0.1,
psfmode="clark", mask=[118,119,147,138], interactive=False,
imsize=256, cell="1arcsec", weighting="briggs", robust=0.0)
# In CASA
calstat=imstat(imagename="ngc3256_cont.image",region="",box="160,50,210,210")
rms=(calstat['rms'][0])
print ">> rms in continuum image: "+str(rms)
calstat=imstat(imagename="ngc3256_cont.image",region="")
peak=(calstat['max'][0])
print ">> Peak in continuum image: "+str(peak)
print ">> Dynamic range in continuum image: "+str(peak/rms)
Image dynamic range is around 34. Peak flux is ~6.7 mJy
# In CASA
imview(raster={'file': 'ngc3256_cont.image', 'colorwedge':T,
'range':[-0.0003, 0.005], 'scaling':-1, 'colormap':
'Rainbow 2'},
out='ngc3256_cont.png', zoom=2)
- Do some phase-only self-cal on the continuum image.
- First use a solution interval of 15 mins
# In CASA
gaincal(vis="ngc3256_line_target.ms", field='NGC*',
caltable="ngc3256_cont_900.Gp",
spw="0:20~53;71~120,1:70~120,2:20~120,3:20~120",
selectdata=True, solint='900s', refant='PM03', calmode='p',
minblperant=3)
# In CASA
plotcal(caltable = 'ngc3256_cont_900.Gp', xaxis = 'time', yaxis =
'phase', poln='X', plotsymbol='o', plotrange = [0,0,-180,180],
iteration = 'spw', figfile='phase_vs_time_XX_900_Gp.png',
subplot = 221)
plotcal(caltable = 'ngc3256_cont_900.Gp', xaxis = 'time', yaxis =
'phase', poln='Y', plotsymbol='o', plotrange = [0,0,-180,180],
iteration = 'spw', figfile='phase_vs_time_YY_900_Gp.png',
subplot = 221)
Phase-only self-cal solutions look good.
# In CASA
applycal(vis='ngc3256_line_target.ms', interp='nearest',
gaintable=['ngc3256_cont_900.Gp'])
Make a new continuum map
# In CASA
os.system('rm -rf ngc3256.cont_sc1*')
clean( vis="ngc3256_line_target.ms", imagename="ngc3256_cont_sc1",
spw="0:20~53;71~120,1:70~120,2:20~120,3:20~120", selectdata=T,
mode="mfs", niter=500, gain=0.1, psfmode="clark",
mask=[118,119,147,138], interactive=False, imsize=256,
cell="1arcsec", weighting="briggs", robust=0.0)
# In CASA
calstat=imstat(imagename="ngc3256_cont_sc1.image",region="",box="160,50,210,210")
rms=(calstat['rms'][0])
print ">> rms in continuum image: "+str(rms)
calstat=imstat(imagename="ngc3256_cont_sc1.image",region="")
peak=(calstat['max'][0])
print ">> Peak in continuum image: "+str(peak)
print ">> Dynamic range in continuum image: "+str(peak/rms)
Continuum image has improved a lot. Dynamic range is ~110 Peak flux is ~9 mJy
# In CASA
imview(raster={'file': 'ngc3256_cont_sc1.image', 'colorwedge':T,
'range':[-0.0003, 0.007], 'scaling':-1, 'colormap':
'Rainbow 2'},
out='ngc3256_cont_sc1.png', zoom=2)
Line images of the galaxy
We split to apply the latest self-cal tables
# In CASA
os.system('rm -rf ngc3256_line_target_sc.ms*')
split(vis='ngc3256_line_target.ms', outputvis='ngc3256_line_target_sc.ms')
Subtract the continuum
# In CASA
uvcontsub(vis = 'ngc3256_line_target_sc.ms',
fitspw="0:20~53;71~120,1:70~120,2:20~120,3:20~120", solint ='inf',
fitorder = 1, fitmode = 'subtract')
Clean the CO line emission
# In CASA
os.system('rm -rf ngc3256_line_CO.*')
clean(vis="ngc3256_line_target_sc.ms", imagename="ngc3256_line_CO",
outframe='LSRK', spw="0:48~78", start='', nchan=31, width='',
restfreq='115.271201800GHz', selectdata=T, mode="channel",
niter=500, gain=0.1, psfmode="clark", mask=[53,50,87,83],
interactive=True, imsize=128, cell="1arcsec",
weighting="briggs", robust=0.0)
Make the moment maps of the CO (1-0) emission
# In CASA
os.system('rm -rf ngc3256_CO1-0.mom.*')
immoments(imagename="ngc3256_line_CO.image", moments=[0,1],
chans="5~24", box="38,38,90,90", axis="spectral",
includepix=[0.045, 10000], outfile="ngc3256_CO1-0.mom")
Clean the CN high line emission
# In CASA
os.system('rm -rf ngc3256_line_CNhi.*')
clean(vis="ngc3256_line_target_sc.ms", imagename="ngc3256_line_CNhi",
outframe='LSRK', spw="1:50~76", start='', nchan=27, width='',
restfreq='113.49097GHz', selectdata=T, mode="channel",
niter=500, gain=0.1, psfmode="clark", mask=[53,50,87,83],
interactive=True, imsize=128, cell="1arcsec",
weighting="briggs", robust=0.0)
Make the moment maps of the CN (high) emission
# In CASA
os.system('rm -rf ngc3256_CNhi.mom.*')
immoments( imagename="ngc3256_line_CNhi.image", moments=[0,1],
chans="5~18", axis="spectral",box="38,38,90,90",
includepix=[0.005, 10000], outfile="ngc3256_CNhi.mom")
Clean the CN low emission
# In CASA
os.system('rm -rf ngc3256_line_CNlo.*')
clean( vis="ngc3256_line_target_sc.ms", imagename="ngc3256_line_CNlo",
outframe='LSRK', spw="1:29~54", start='', nchan=26, width='',
restfreq='113.18666GHz', selectdata=T, mode="channel",
niter=300, gain=0.1, psfmode="clark", mask=[53,50,87,83],
interactive=True, imsize=128, cell="1arcsec",
weighting="briggs", robust=0.0)
Make the moment maps of the CN (low) emission
# In CASA
os.system('rm -rf ngc3256_CNlo.mom.*')
immoments( imagename="ngc3256_line_CNlo.image", moments=[0,1],
chans="8~18", axis="spectral", box="38,38,90,90",
includepix=[0.003, 10000], outfile="ngc3256_CNlo.mom")
Make images of the CO(1-0) emission
# In CASA
imview(contour={'file': 'ngc3256_CO1-0.mom.integrated','range': []},
raster={'file': 'ngc3256_CO1-0.mom.weighted_coord','range': [2630,2920],'colorwedge':T,
'colormap': 'Rainbow 2'},
out='CO_velfield.png')
# In CASA
imview(contour={'file': 'ngc3256_CO1-0.mom.weighted_coord','levels': [2650,2700,2750,2800,2850,2900],'base':0,'unit':1},
raster={'file': 'ngc3256_CO1-0.mom.integrated','colorwedge':T,
'colormap': 'Rainbow 2'},
out='CO_map.png')
Make images of the CN high emission
# In CASA
imview(contour={'file': 'ngc3256_CNhi.mom.integrated','range': []},
raster={'file': 'ngc3256_CNhi.mom.weighted_coord','range': [2630,2920],'colorwedge':T,
'colormap': 'Rainbow 2'},
out='CNhi_velfield.png')
# In CASA
imview(contour={'file': 'ngc3256_CNhi.mom.weighted_coord','levels': [2650,2700,2750,2800,2850,2900],'base':0,'unit':1},
raster={'file': 'ngc3256_CNhi.mom.integrated','colorwedge':T,
'colormap': 'Rainbow 2'},
out='CNhi_map.png')
Make images of the CN low emission
# In CASA
imview(contour={'file': 'ngc3256_CNlo.mom.integrated','range': []},
raster={'file': 'ngc3256_CNlo.mom.weighted_coord','range': [2630,2920],'colorwedge':T,
'colormap': 'Rainbow 2'},
out='CNlo_velfield.png')
# In CASA
imview(contour={'file': 'ngc3256_CNlo.mom.weighted_coord','levels': [2650,2700,2750,2800,2850,2900],'base':0,'unit':1},
raster={'file': 'ngc3256_CNlo.mom.integrated','colorwedge':T,
'colormap': 'Rainbow 2'},
out='CNlo_map.png')
