##########################################################################
# This casapy script was prepared for the EU ARC Cycle 2 PI CASA tutorial.
# It was first made to work correctly.
# Then, a number of problems were introduced.
# It is your task to fix them.
# You will find a list of the problems further below.
# Look for the string 'PROBLEM' to find all problems.
#
# The script uses some infrastructure to make it easier for you
# to run and repeat the script step by step:
# Further below you find the "step_title" Python dictionary which
# shows you what is happening in each analysis step.
# You can start this script by starting casapy and entering
# execfile('reduce-m100-casa4.2.2-with-problems.py')
# By default it will run all steps.
# (Of course it will only run correctly when you have fixed all problems.)
# If you would like to run individual steps, you can set the
# variable mysteps before starting the script, e.g.:
#
# mysteps = [4,5,6]
# execfile('reduce-m100-casa4.2.2-with-problems.py')
#
# This example would only execute steps 4, 5, and 6.
# Setting mysteps to an empty list [] will execute all steps.
#
#
# List of problems in the individual analysis steps of this script:
# Step 0: Write a flagdata command to flag the channels
# 239, 447/448, 720/721 and 2847/2848 in all SPWs.
# Step 4: Flag the CO line of Titan and noisy data at uvdist < 40 klambda.
# Step 5: Determine the solint and calmode parameters in gaincal.
# Step 6: Determine the field, solint, and calmode parameters in gaincal.
# Step 7: Determine the gaintable parameter in gaincal.
# Step 8: Determine the right values for reference and transfer in fluxscale.
# Step 9: Decide on the value for calwt and complete the last applycal command.
# Step 13: Determine the missing parameters in split.
# Step 15: Determine the solint parameter in gaincal.
# Step 17: Determine the mode, imagermode, and mask parameters in clean.
# Step 18: Determine the fitspw parameter in uvcontsub.
# Step 19: The mask test-M100line-orig.mask is missing. Generate it using
# interactive clean.
# Steps 19+20: Determine the field, outframe, and weighting parameters in clean
# Step 21: Determine the axis and includepix parameters in immoments
#
##########################################################################
# EU ALMA Regional Centre CASA Tutorial for ALMA Cycle 2
# Sample ALMA data PI analysis script
# "M100"
step_title = { 0: 'Flagging \n'\
' (Certain channels are bad, certain antennas have problems during a few short periods)',
1: 'Rebin to a reduced resolution of approx. 10 km/s \n'\
' (in order to save processing time in the following steps)',
2: 'Fast phase-only gaincal for bandpass \n'\
' (This is essentially a one-step self-cal on the bandpass calibrator to take out fast phase fluctuations)',
3: 'Bandpass \n'\
' (Frequency-dependence calibration: applying the phase cal table from the previous step on-the-fly, \n'\
' we generate the bandpass cal table)',
4: 'Setjy \n'\
' (Titan is not a point source. After flagging a strong CO line in channels 200~479 of Titan, we set \n'\
' the correct model for Titan in all SPWs)',
5: 'Fast phase-only gaincal \n'\
' (Similar to step 2, we do a phase-only selfcal on the bandpass, the primary phase, and the flux \n'\
' calibrator to improve S/N)',
6: 'Slow phase-only gaincal \n'\
' (The phase-only cal to be used for the secondary phase calibrator and the target. Long solint in order \n'\
' to improve stability of the interpol.solution for the target)',
7: 'Slow amp and phase gaincal \n'\
' (Time-dependence calibration: using the bandpass solution from step 3 and the fast phase-only solution \n'\
' from step 5 on-the-fly)',
8: 'Fluxscale \n'\
' (using the amp and phase solution from step 7 on-the-fly, the flux scale is transferred from Titan to \n'\
' the phase and the bandpass calibrator)',
9: 'Applycal \n'\
' (apply the bandpass solution from step 3, the phase solutions from steps 5 or 6 respectively, and the \n'\
' flux solution from step 8)',
10: 'Test image of the secondary phase cal',
11: 'Test image of the primary phase cal',
12: 'Test image of Titan',
13: 'Split out corrected data and time average \n'\
' (Now that the time dependence is calibrated, we can reduce the time resolution to 60s and save time\n'\
' in the following steps)',
14: 'Concatenate data \n'\
' (To prepare imaging, the two datasets X220 and X54 are merged)',
15: 'Adjust fluxscale \n'\
' (Since we calibrated the flux scale independently for X220 and X54, we need to align the two flux \n'\
' calibrations assuming that the flux of the phasecal did not change within the short time between X220 and X54)',
16: 'Split out the corrected M100 data \n'\
' (in order to further shrink the input MS for imaging)',
17: 'Continuum image of M100 \n'\
' (excluding the channels of the CO(1-0) line determined with plotms)',
18: 'Determine and subtract continuum \n'\
' (only for spw 0 where the CO line is)',
19: 'Test image of central field \n'\
' (averaging over the CO line channels)',
20: 'Clean CO(1-0) line cube mosaic \n'\
' (for the CO line channels)',
21: 'Make moment maps \n'\
' (for the integrated sum and the gradient - the velocity field - of the cube;\n'\
' also demonstrate how to make overlay plots of the maps using imview)'
}
# global defs
basename=['X54','X220']
makeplots=True
therefant = 'DV01'
############################################################################
# Some infrastructure to make repeating individual parts
# of this workflow more convenient.
totaltime = 0
inittime = time.time()
ttime = inittime
steptime = []
thesteps = []
def timing():
global totaltime
global inittime
global ttime
global steptime
global step_title
global mystep
global thesteps
thetime = time.time()
dtime = thetime - ttime
steptime.append(dtime)
totaltime += dtime
ttime = thetime
casalog.origin('TIMING')
casalog.post( 'Step '+str(mystep)+': '+step_title[mystep].splitlines()[0], 'WARN')
casalog.post( 'Time now: '+str(ttime), 'WARN')
casalog.post( 'Time used this step: '+str(dtime), 'WARN')
casalog.post( 'Total time used so far: ' + str(totaltime), 'WARN')
casalog.post( 'Step Time used (s) Fraction of total time (percent) [description]', 'WARN')
for i in range(0, len(steptime)):
casalog.post( ' '+str(thesteps[i])+' '+str(steptime[i])+' '+str(steptime[i]/totaltime*100.)
+' ['+step_title[thesteps[i]].splitlines()[0]+']', 'WARN')
try:
print 'List of steps to be executed ...', mysteps
thesteps = mysteps
except:
print 'global variable mysteps not set.'
if(thesteps==[]):
thesteps = range(0,len(step_title))
print 'mysteps empty. Executing all steps: ', thesteps
# The Python variable 'mysteps' will control which steps
# are executed when you start the script using
# execfile('reduce-m100.py')
# e.g. setting
# mysteps = [2,3,4]
# before starting the script will make the script execute
# only steps 2, 3, and 4
# Setting mysteps = [] will make it execute all steps.
##########################################################################
# Begin analysis
# flagging
mystep = 0
if(mystep in thesteps):
print 'Step ', mystep, step_title[mystep]
os.system('ln -sf uid___A002_X2a5c2f_X220.ms.split X220-line.ms')
os.system('cp -R uid___A002_X2a5c2f_X220.ms.split.flagversions X220-line.ms.flagversions')
os.system('ln -sf uid___A002_X2a5c2f_X54.ms.split X54-line.ms')
os.system('cp -R uid___A002_X2a5c2f_X54.ms.split.flagversions X54-line.ms.flagversions')
for name in basename:
flagmanager(vis=name+'-line.ms', mode='restore', versionname='apriori')
# Edge channels
flagdata(vis=name+'-line.ms',
field='', mode='manual',
spw='0~3:0~10;3800~3839',
flagbackup = F)
# Channels 239, 447/448, 720/721 and 2847/2848 are off in all SPWs in both MSs
# PROBLEM: write a flagdata command to flag these channels
# some ints are off
flagdata(vis='X220-line.ms',
mode='manual',
timerange='19:52:55~19:53:04',
flagbackup = F)
flagdata(vis='X54-line.ms',
antenna='PM01',
timerange='19:03:35~19:03:42',
mode='manual',
flagbackup = F)
flagdata(vis='X54-line.ms',
antenna='DV04',
timerange='19:38:45~19:38:55',
mode='manual',
flagbackup = F)
timing()
# Bin it up to lower spectral resolution to about 10 km/s
mystep = 1
if(mystep in thesteps):
print 'Step ', mystep, step_title[mystep]
for name in basename:
os.system('rm -rf '+name+'-line-vs.ms')
split(vis=name+'-line.ms',
outputvis=name+'-line-vs.ms',
datacolumn='data',
width='8')
flagmanager(vis=name+'-line-vs.ms', mode='save', versionname='step1')
timing()
# Fast phase-only gaincal for bandpass
mystep = 2
if(mystep in thesteps):
print 'Step ', mystep, step_title[mystep]
for name in basename:
os.system('rm -rf cal-'+name+'-BPint.Gp')
gaincal(vis=name+'-line-vs.ms',
caltable='cal-'+name+'-BPint.Gp',
spw='*:190~310',
field='*Bandpass*',
selectdata=F, solint='int', refant=therefant, calmode='p')
if(makeplots):
plotcal(caltable='cal-'+name+'-BPint.Gp',
xaxis = 'time', yaxis = 'phase',
poln='X', plotsymbol='o', plotrange = [0,0,-180,180],
iteration = 'spw',
figfile='cal-'+name+'-phase_vs_time_XX.BPint.Gp.png', subplot = 221)
plotcal(caltable='cal-'+name+'-BPint.Gp',
xaxis = 'time', yaxis = 'phase',
poln='Y', plotsymbol='o', plotrange = [0,0,-180,180],
iteration = 'spw',
figfile='cal-'+name+'-phase_vs_time_YY.BPint.Gp.png', subplot = 221)
timing()
# Bandpass
mystep = 3
if(mystep in thesteps):
print 'Step ', mystep, step_title[mystep]
for name in basename:
os.system('rm -rf cal-'+name+'.B1')
bandpass(vis=name+'-line-vs.ms',
caltable='cal-'+name+'.B1',
field='*Bandpass*',
bandtype='B', fillgaps=10, solnorm = T, combine='',
selectdata=F,
solint='inf',
refant=therefant,
gaintable='cal-'+name+'-BPint.Gp')
if(makeplots):
plotbandpass(caltable = 'cal-'+name+'.B1',
xaxis='freq', yaxis='phase',
showatm=True, interactive=False,
plotrange = [0,0,-70,70],
figfile='cal-'+name+'-phase.B1')
plotbandpass(caltable = 'cal-'+name+'.B1',
xaxis='freq', yaxis='amp',
showatm=True, interactive=False,
figfile='cal-'+name+'-amplitude.B1')
timing()
# Setjy
# Noisy for uvdistances less than 40 klambda
# Strong line for Titan is obvious in spw 0
mystep = 4
if(mystep in thesteps):
print 'Step ', mystep, step_title[mystep]
for name in basename:
flagmanager(vis=name+'-line-vs.ms', mode='restore', versionname='step1')
flagdata(vis=name+'-line-vs.ms',
field=PROBLEM,
mode=PROBLEM,
uvrange=PROBLEM,
spw=PROBLEM,
flagbackup = F)
flagdata(vis=name+'-line-vs.ms',
field=PROBLEM,
mode=PROBLEM,
uvrange=PROBLEM,
spw=PROBLEM,
flagbackup = F)
setjy(vis=name+'-line-vs.ms',
field='Titan',
standard='Butler-JPL-Horizons 2012',
spw='0,1,2,3')
timing()
# Fast phase-only gaincal
mystep = 5
if(mystep in thesteps):
print 'Step ', mystep, step_title[mystep]
for name in basename:
os.system('rm -rf cal-'+name+'-int.Gp')
gaincal(vis=name+'-line-vs.ms',
caltable='cal-'+name+'-int.Gp',
spw='*:25~455',
field='*Phase*,*Band*,Titan',
gaintable='cal-'+name+'.B1',
selectdata=F, solint=PROBLEM,
refant=therefant, calmode=PROBLEM)
if(makeplots):
plotcal(caltable='cal-'+name+'-int.Gp',
xaxis = 'time', yaxis = 'phase',
poln='X', plotsymbol='o', plotrange = [0,0,-180,180],
iteration = 'spw',
figfile='cal-'+name+'-phase_vs_time_XX.int.Gp.png', subplot = 221)
plotcal(caltable='cal-'+name+'-int.Gp',
xaxis = 'time', yaxis = 'phase',
poln='Y', plotsymbol='o', plotrange = [0,0,-180,180],
iteration = 'spw',
figfile='cal-'+name+'-phase_vs_time_YY.int.Gp.png', subplot = 221)
timing()
# Slow phase-only gaincal
mystep = 6
if(mystep in thesteps):
print 'Step ', mystep, step_title[mystep]
for name in basename:
os.system('rm -rf cal-'+name+'-scan.Gp')
gaincal(vis=name+'-line-vs.ms',
caltable='cal-'+name+'-scan.Gp',
spw='*:25~455',
field=PROBLEM,
gaintable='cal-'+name+'.B1',
selectdata=F, solint='PROBLEM',
refant=therefant, calmode='PROBLEM')
if(makeplots):
plotcal(caltable='cal-'+name+'-scan.Gp',
xaxis = 'time', yaxis = 'phase',
poln='X', plotsymbol='o', plotrange = [0,0,-180,180],
iteration = 'spw',
figfile='cal-'+name+'-phase_vs_time_XX.scan.Gp.png', subplot = 221)
plotcal(caltable='cal-'+name+'-scan.Gp',
xaxis = 'time', yaxis = 'phase',
poln='Y', plotsymbol='o', plotrange = [0,0,-180,180],
iteration = 'spw',
figfile='cal-'+name+'-phase_vs_time_YY.scan.Gp.png', subplot = 221)
timing()
# Slow amplitude and phase gaincal
mystep = 7
if(mystep in thesteps):
print 'Step ', mystep, step_title[mystep]
for name in basename:
os.system('rm -rf cal-'+name+'-scan.Gap')
gaincal(vis=name+'-line-vs.ms',
caltable='cal-'+name+'-scan.Gap',
spw='*:25~455',
field='*Phase*,*Band*,Titan',
gaintable=['PROBLEM'],
selectdata=F, solint='inf',
refant=therefant, calmode='ap')
if(makeplots):
plotcal(caltable='cal-'+name+'-scan.Gap',
xaxis = 'time', yaxis = 'amp',
poln='X', plotsymbol='o',
iteration = 'spw',
figfile='cal-'+name+'-amp_vs_time_XX.scan.Gap.png', subplot = 221)
plotcal(caltable='cal-'+name+'-scan.Gap',
xaxis = 'time', yaxis = 'amp',
poln='Y', plotsymbol='o',
iteration = 'spw',
figfile='cal-'+name+'-amp_vs_time_YY.scan.Gap.png', subplot = 221)
timing()
# Fluxscale
mystep = 8
if(mystep in thesteps):
print 'Step ', mystep, step_title[mystep]
for name in basename:
fluxscale(vis=name+'-line-vs.ms',caltable='cal-'+name+'-scan.Gap',
fluxtable='cal-'+name+'.flux',reference=PROBLEM,transfer=PROBLEM,
incremental=False)
timing()
# Applycal
mystep = 9
if(mystep in thesteps):
print 'Step ', mystep, step_title[mystep]
for name in basename:
# to the bandpass cal
applycal(vis=name+'-line-vs.ms',field='*Band*',
gaintable=['cal-'+name+'.B1','cal-'+name+'-int.Gp','cal-'+name+'.flux'],
interp='linear,linear',
gainfield=['*Band*','*Band*','*Band*'],
calwt=PROBLEM,
flagbackup=T)
# to the secondary phase cal
applycal(vis=name+'-line-vs.ms',field='3c273 - Phase',
gaintable=['cal-'+name+'.B1','cal-'+name+'-scan.Gp','cal-'+name+'.flux'],
interp='linear,linear',
gainfield=['*Band*','1224*','1224*'],
calwt=PROBLEM,
flagbackup=T)
# to the primary phase cal
applycal(vis=name+'-line-vs.ms',field='1224*',
gaintable=['cal-'+name+'.B1','cal-'+name+'-int.Gp','cal-'+name+'.flux'],
interp='linear,linear',
gainfield=['*Band*','1224*','1224*'],
calwt=PROBLEM,
flagbackup=T)
# to Titan
applycal(vis=name+'-line-vs.ms',field='Titan',
gaintable=['cal-'+name+'.B1','cal-'+name+'-int.Gp','cal-'+name+'.flux'],
interp='linear,linear',
gainfield=['*Band*','Titan','Titan'],
calwt=PROBLEM,
flagbackup=T)
# to M100
applycal(PROBLEM)
timing()
# Test image of the secondary phase cal
mystep = 10
if(mystep in thesteps):
print 'Step ', mystep, step_title[mystep]
# Note: the WARNING about the beam is benign
for name in basename:
os.system('rm -rf test-'+name+'-sec_phasecal*')
clean(vis=name+'-line-vs.ms',
imagename='test-'+name+'-sec_phasecal',
field='3c*Ph*',spw='0~3',
nterms=2,
mode='mfs',niter=100,
interactive=False,
mask=[96, 96, 104, 104],imsize=200,cell='0.5arcsec')
if makeplots:
for name in basename:
imview(raster={'file': 'test-'+name+'-sec_phasecal.image.tt0', 'colorwedge':T,
'range':[-0.02, 8.0], 'scaling':-1.5, 'colormap':'Rainbow 2'},
out='test-'+name+'-sec_phasecal.png', zoom=1)
timing()
# Test image of the primary phase cal
mystep = 11
if(mystep in thesteps):
print 'Step ', mystep, step_title[mystep]
for name in basename:
os.system('rm -rf test-'+name+'-prim_phasecal*')
clean(vis=name+'-line-vs.ms',
imagename='test-'+name+'-prim_phasecal',
field='1224*',spw='0~3',
nterms=2,
mode='mfs',niter=100,
interactive=False,
mask=[96, 96, 104, 104],imsize=200,cell='0.5arcsec')
if makeplots:
for name in basename:
imview(raster={'file': 'test-'+name+'-prim_phasecal.image.tt0', 'colorwedge':T,
'range':[-0.005, 0.9], 'scaling':-2.5, 'colormap':'Rainbow 2'},
out='test-'+name+'-prim_phasecal.png', zoom=1)
calstat=imstat(imagename='test-'+name+'-prim_phasecal.image.tt0', region='', box='30,30,170,80')
rms=(calstat['rms'][0])
print '>> rms in phase calibrator image: '+str(rms)
calstat=imstat(imagename='test-'+name+'-prim_phasecal.image.tt0', region='')
peak=(calstat['max'][0])
print '>> Peak in phase calibrator image: '+str(peak)
print '>> Dynamic range in phase calibrator image: '+str(peak/rms)
timing()
# Test image of Titan
mystep = 12
if(mystep in thesteps):
print 'Step ', mystep, step_title[mystep]
for name in basename:
os.system('rm -rf test-'+name+'-Titan*')
clean(vis=name+'-line-vs.ms',
imagename='test-'+name+'-Titan',
field='Titan',spw='0~3',
mode='mfs',niter=100,
interactive=False,
mask=[96, 96, 104, 104],imsize=200,cell='0.5arcsec')
timing()
# Split out corrected data and time average at the same time
mystep = 13
if(mystep in thesteps):
print 'Step ', mystep, step_title[mystep]
for name in basename:
os.system('rm -rf '+name+'-corrected.ms*')
split(vis=name+'-line-vs.ms', outputvis=name+'-corrected.ms') # PROBLEM
timing()
# Concat
mystep = 14
if(mystep in thesteps):
print 'Step ', mystep, step_title[mystep]
os.system('rm -rf all.ms*')
concat(vis=['X54-corrected.ms', 'X220-corrected.ms'],
concatvis='all.ms', copypointing=False)
timing()
# Adjust fluxscale
mystep = 15
if(mystep in thesteps):
print 'Step ', mystep, step_title[mystep]
#X54
#INFO fluxscale:::: Flux density for 1224+213 Phase in SpW=0 is: 1.12613 +/- 0.0111589 (SNR = 100.918, N = 26)
#INFO fluxscale:::: Flux density for 1224+213 Phase in SpW=1 is: 1.13117 +/- 0.0113115 (SNR = 100.002, N = 26)
#INFO fluxscale:::: Flux density for 1224+213 Phase in SpW=2 is: 1.18905 +/- 0.00829696 (SNR = 143.311, N = 26)
#INFO fluxscale:::: Flux density for 1224+213 Phase in SpW=3 is: 1.19364 +/- 0.00717581 (SNR = 166.342, N = 26)
#X220
#INFO fluxscale:::: Flux density for 1224+213 Phase in SpW=0 is: 1.08145 +/- 0.0104661 (SNR = 103.329, N = 26)
#INFO fluxscale:::: Flux density for 1224+213 Phase in SpW=1 is: 1.10837 +/- 0.0083772 (SNR = 132.308, N = 26)
#INFO fluxscale:::: Flux density for 1224+213 Phase in SpW=2 is: 1.17593 +/- 0.00666044 (SNR = 176.554, N = 26)
#INFO fluxscale:::: Flux density for 1224+213 Phase in SpW=3 is: 1.1746 +/- 0.00597671 (SNR = 196.529, N = 26)
# INFO listobs::ms::summary+ SpwID #Chans Frame Ch1(MHz) ChanWid(kHz) TotBW(kHz) Corrs
# INFO listobs::ms::summary+ 0 480 TOPO 113728.128 3906.25 1875000 XX YY
# INFO listobs::ms::summary+ 1 480 TOPO 111853.128 3906.25 1875000 XX YY
# INFO listobs::ms::summary+ 2 480 TOPO 103661.722 -3906.25 1875000 XX YY
# INFO listobs::ms::summary+ 3 480 TOPO 101849.222 -3906.25 1875000 XX YY
# calculate spectral index of primary phase cal
# S = fluxdensity * (freq/reffreq)**spix
# (alternatively run setjy separately on each SPW)
freq0 = 101849.222 - 480./2. * 3906.25E-3
freq1 = 113728.128 + 480./2. * 3906.25E-3
flux0 = 1.19364
flux1 = 1.12613
myspix1 = log(flux0/flux1)/log(freq0/freq1)
flux0 = 1.1746
flux1 = 1.08145
myspix2 = log(flux0/flux1)/log(freq0/freq1)
myspix = (myspix1 + myspix2)/2.
print "Spectral index of primary phase cal: ", myspix
# find myspix approx. -0.55, SMA database says -0.51 +- 0.05 at 3 mm, -0.67 +- 0.04 at 1 mm; we are at ca. 2.9 mm
setjy(vis='all.ms',
field='1224*',
standard = 'manual',
fluxdensity=flux0,
reffreq=str(freq0)+'MHz',
spix=myspix,
spw='0,1,2,3')
os.system('rm -rf cal-fluxadjust.Ga')
gaincal(vis='all.ms',
caltable='cal-fluxadjust.Ga',
spw='*:25~455',
field='1224*',
gaintable=[],
selectdata=F, solint=PROBLEM,
refant=therefant, calmode='a')
if(makeplots):
plotcal(caltable='cal-fluxadjust.Ga',
xaxis = 'time', yaxis = 'amp',
poln='X', plotsymbol='o',
iteration = 'spw',
figfile='cal-all-amp_vs_time_XXfluxadjust.Ga.png', subplot = 221)
applycal(vis='all.ms',field='M100',
gaintable=['cal-fluxadjust.Ga'],
interp=['linear'],
gainfield=['1224*'],
calwt=T,
flagbackup=T)
timing()
# split out the corrected M100 data
mystep = 16
if(mystep in thesteps):
print 'Step ', mystep, step_title[mystep]
os.system('rm -rf M100all_lores.ms*')
split(vis='all.ms', outputvis='M100all_lores.ms',
field = 'M100',
datacolumn='corrected')
timing()
# Continuum image
mystep = 17
if(mystep in thesteps):
print 'Step ', mystep, step_title[mystep]
os.system('rm -rf M100cont.*')
clean(vis = 'M100all_lores.ms',
imagename = 'M100cont',
field='2~47',
spw='0:10~210;256~440,1~3:10~460', # exclude CO(1-0) line emission
mode = 'PROBLEM',
niter = 1000,
mask=[0,0,0,0], #PROBLEM
imagermode = 'PROBLEM',
interactive = F,
imsize = 200,
cell = '0.5arcsec',
phasecenter='J2000 12h22m54.9 +15d49m15')
# Continuum peak is 0.5 mJy. Too weak for self-cal...
timing()
# uvcontsub
mystep = 18
if(mystep in thesteps):
print 'Step ', mystep, step_title[mystep]
os.system('rm -rf M100all_lores.ms.c*')
uvcontsub(vis='M100all_lores.ms',field='',fitspw='PROBLEM',
combine='',solint='inf',fitorder=1,spw='0',want_cont=False)
timing()
# Test image of central field
mystep = 19
if(mystep in thesteps):
print 'Step ', mystep, step_title[mystep]
os.system('rm -rf test-M100line.*')
clean(
vis='M100all_lores.ms.contsub',
imagename='test-M100line',
field='PROBLEM',
spw='0:231~248',
mode='mfs',
niter=500,gain=0.1,threshold='0.0mJy',
imagermode='csclean',
interactive=False,
mask='test-M100line-orig.mask', # PROBLEM: mask missing
outframe='PROBLEM',
veltype='radio',
imsize=200,cell='0.5arcsec',
phasecenter='',
stokes='I',
weighting='PROBLEM',robust=0.5,
npercycle=100,cyclefactor=1.5,cyclespeedup=-1)
timing()
# Clean line cube mosaic
mystep = 20
if(mystep in thesteps):
print 'Step ', mystep, step_title[mystep]
os.system('rm -rf M100line.*')
clean(vis='M100all_lores.ms.contsub',imagename='M100line',
field='PROBLEM',
spw='0:220~259',
mode='channel',
niter=1000,gain=0.1,threshold='5mJy',psfmode='clark',
imagermode='mosaic',ftmachine='mosaic',
mosweight=False,
scaletype='SAULT',
interactive=False, mask='M100line300x300.mask',
nchan=40,start=220,
width=1,
outframe='PROBLEM',
veltype='radio',
imsize=300,cell='1arcsec', # could be set to 600, 0.5 arcsec but we need to save time
phasecenter='J2000 12h22m54.9 +15d49m10',
restfreq='115.271201800GHz',
stokes='I',
weighting='PROBLEM',
robust=0.5,
pbcor=False,
minpb=0.2,
npercycle=100,cyclefactor=1.5,cyclespeedup=-1)
timing()
# Moments
mystep = 21
if(mystep in thesteps):
print 'Step ', mystep, step_title[mystep]
os.system('rm -rf M100-CO.mom?')
immoments(imagename='M100line.image',
moments=[0],
axis='PROBLEM',
region='',box='50,55,255,250',
chans='7~35',
mask='',
outfile='M100-CO.mom0',
includepix=[PROBLEM])
immoments(imagename='M100line.image',
moments=[1],
axis='PROBLEM',
region='',box='50,55,255,250',
chans='7~35',
mask='',
outfile='M100-CO.mom1',
includepix=[PROBLEM])
if makeplots:
os.system('rm -rf M100-CO_velfield.png')
imview(contour={'file': 'M100-CO.mom0','levels':
[1,2,5,10,20,40,80,160],'base':0,'unit':1},
raster={'file': 'M100-CO.mom1','range': [1440,1700],
'colorwedge':T, 'colormap': 'Rainbow 2'}, out='M100-CO_velfield.png')
os.system('rm -rf M100-CO_map.png')
imview(contour={'file': 'M100-CO.mom1','levels':
[1430,1460,1490,1520,1550,1580,1610,1640,1670,1700],'base':0,'unit':1},
raster={'file': 'M100-CO.mom0', 'colorwedge':T,
'colormap': 'Rainbow 2','scaling':-1.8,'range': [0.5,20]},
out='M100-CO_map.png')
os.system('rm -rf M100-CO_contmap.png')
imview(contour={'file': 'M100cont.image','levels':
[0.00025,0.0004],'base':0,'unit':1},
zoom=3,
raster={'file': 'M100-CO.mom0', 'colorwedge':T,
'colormap': 'Rainbow 2','scaling':0,'range': [0.8,40]},
out='M100-CO_contmap.png')
timing()
print 'Done.'