# CASA sit-together 2 # analysis script for NGC 4258 = 12190+47182 # D. Petry # K-band spectral analysis, H2O Maser lines "1" and "3", i.e. 490 km/s and 1310 km/s # according to the proposal doall = False try: print "List of steps to be executed ... ", mysteps except: print "global variable mysteps not set." mysteps = [] if (mysteps==[]): print "mysteps empty. Executing all steps." doall = True # my globals basename = 'ah847_1-k-lines1+3' steptitle = { 0 : "data import", 1 : "data inspection", 2 : "data selection", 3 : "flagging", 4 : "set the flux scale", 5 : "bandpass calibration", 6 : "gain calibration", 7 : "apply calibration", 8 : "image the low-res spw", 9 : "image the hi-res spw" } # import mystep = 0 if(mystep in mysteps or doall): print "Step ", mystep, steptitle[mystep] os.system('rm -rf '+basename+'*') importvla( archivefiles = 'AH847_1', vis = basename+'.ms', bandname = 'K', project = '', starttime = '', stoptime = '', applytsys = True, autocorr = False, antnamescheme = 'new', keepblanks = False, evlabands = False, async = False ) # inspection mystep = 1 if(mystep in mysteps or doall): print "Step ", mystep, steptitle[mystep] listobs( vis = basename+'.ms' ) ## listobs Observer: unavailable Project: AH847 ## listobs Observation: VLA ## listobs Data records: 479700 Total integration time = 9370 seconds ## listobs Observed from 22-May-2004/01:04:35.0 to 22-May-2004/03:40:45.0 (TAI) ## listobs::ms::summary ## listobs ObservationID = 0 ArrayID = 0 ## listobs Date Timerange (TAI) Scan FldId FieldName nVis Int(s) SpwIds ## listobs 22-May-2004/01:04:35.0 - 01:05:45.0 1 0 1219+4KC 5200 10 [0, 1] ## listobs 01:06:05.0 - 01:09:15.0 2 1 12190+47182 13000 10 [0, 1] - target with continuum SPW ## listobs 01:09:25.0 - 01:10:15.0 3 2 12191+48299 3900 10 [0, 1] ## listobs 01:10:25.0 - 01:13:15.0 4 1 12190+47182 11700 10 [0, 1] ## listobs 01:13:25.0 - 01:14:15.0 5 2 12191+48299 3900 10 [0, 1] ## listobs 01:14:25.0 - 01:17:15.0 6 1 12190+47182 11700 10 [0, 1] ## listobs 01:17:25.0 - 01:18:15.0 7 2 12191+48299 3900 10 [0, 1] ## listobs 01:18:25.0 - 01:21:15.0 8 1 12190+47182 11700 10 [0, 1] ## listobs 01:21:25.0 - 01:22:15.0 9 2 12191+48299 3900 10 [0, 1] ## listobs 01:22:25.0 - 01:25:15.0 10 1 12190+47182 11700 10 [0, 1] ## listobs 01:25:25.0 - 01:26:15.0 11 2 12191+48299 3900 10 [0, 1] ## listobs 01:26:25.0 - 01:29:15.0 12 1 12190+47182 11700 10 [0, 1] ## listobs 01:29:25.0 - 01:30:15.0 13 2 12191+48299 3900 10 [0, 1] ## listobs 01:30:25.0 - 01:33:15.0 14 1 12190+47182 11700 10 [0, 1] ## listobs 01:33:25.0 - 01:34:15.0 15 2 12191+48299 3900 10 [0, 1] ## listobs 01:34:25.0 - 01:37:15.0 16 1 12190+47182 11700 10 [0, 1] ## listobs 01:37:25.0 - 01:38:15.0 17 2 12191+48299 3900 10 [0, 1] ## listobs 01:38:25.0 - 01:41:15.0 18 1 12190+47182 11700 10 [0, 1] ## listobs 01:41:25.0 - 01:41:55.0 19 2 12191+48299 2600 10 [0, 1] ## listobs 01:42:22.5 - 01:43:37.5 20 3 1219+48 3900 15 [2, 3] ## listobs 01:44:22.5 - 01:46:52.5 21 1 12190+47182 7150 15 [2, 3] - target with lower and higher res SPWs ## listobs 01:47:07.5 - 01:47:52.5 22 2 12191+48299 2600 15 [2, 3] - phase cal. - " - ## listobs 01:48:07.5 - 01:52:22.5 23 1 12190+47182 11700 15 [2, 3] ## listobs 01:52:37.5 - 01:53:22.5 24 2 12191+48299 2600 15 [2, 3] ## listobs 01:53:37.5 - 01:55:22.5 25 1 12190+47182 5200 15 [2, 3] ## listobs 01:55:37.5 - 01:56:22.5 26 2 12191+48299 2600 15 [2, 3] ## listobs 01:56:37.5 - 02:00:52.5 27 1 12190+47182 11700 15 [2, 3] ## listobs 02:01:07.5 - 02:01:52.5 28 2 12191+48299 2600 15 [2, 3] ## listobs 02:02:07.5 - 02:03:52.5 29 1 12190+47182 5200 15 [2, 3] ## listobs 02:04:07.5 - 02:04:52.5 30 2 12191+48299 2600 15 [2, 3] ## listobs 02:05:07.5 - 02:09:22.5 31 1 12190+47182 11700 15 [2, 3] ## listobs 02:09:37.5 - 02:10:22.5 32 2 12191+48299 2600 15 [2, 3] ## listobs 02:10:37.5 - 02:12:22.5 33 1 12190+47182 5200 15 [2, 3] ## listobs 02:12:37.5 - 02:13:22.5 34 2 12191+48299 2600 15 [2, 3] ## listobs 02:13:37.5 - 02:14:07.5 35 1 12190+47182 1950 15 [2, 3] ## listobs 02:14:52.5 - 02:16:07.5 36 3 1219+48 3900 15 [4, 3] ## listobs 02:16:37.5 - 02:18:52.5 37 1 12190+47182 6500 15 [4, 3] - target with two high res mode spws ## listobs 02:19:07.5 - 02:19:52.5 38 2 12191+48299 2600 15 [4, 3] - phase cal ## listobs 02:20:07.5 - 02:21:52.5 39 1 12190+47182 5200 15 [4, 3] ## listobs 02:22:07.5 - 02:22:37.5 40 2 12191+48299 1950 15 [4, 3] ## listobs 02:22:52.5 - 02:24:37.5 41 1 12190+47182 5200 15 [4, 3] ## listobs 02:24:52.5 - 02:25:22.5 42 2 12191+48299 1950 15 [4, 3] ## listobs 02:25:37.5 - 02:27:22.5 43 1 12190+47182 5200 15 [4, 3] ## listobs 02:27:37.5 - 02:28:22.5 44 2 12191+48299 2600 15 [4, 3] ## listobs 02:28:37.5 - 02:30:22.5 45 1 12190+47182 5200 15 [4, 3] ## listobs 02:30:37.5 - 02:31:07.5 46 2 12191+48299 1950 15 [4, 3] ## listobs 02:31:22.5 - 02:33:07.5 47 1 12190+47182 5200 15 [4, 3] ## listobs 02:33:22.5 - 02:33:52.5 48 2 12191+48299 1950 15 [4, 3] ## listobs 02:34:07.5 - 02:35:52.5 49 1 12190+47182 5200 15 [4, 3] ## listobs 02:36:07.5 - 02:36:52.5 50 2 12191+48299 2600 15 [4, 3] ## listobs 02:37:07.5 - 02:38:52.5 51 1 12190+47182 5200 15 [4, 3] ## listobs 02:39:07.5 - 02:39:37.5 52 2 12191+48299 1950 15 [4, 3] ## listobs 02:39:52.5 - 02:41:37.5 53 1 12190+47182 5200 15 [4, 3] ## listobs 02:41:52.5 - 02:42:22.5 54 2 12191+48299 1950 15 [4, 3] ## listobs 02:42:37.5 - 02:44:22.5 55 1 12190+47182 5200 15 [4, 3] ## listobs 02:44:37.5 - 02:45:22.5 56 2 12191+48299 2600 15 [4, 3] ## listobs 02:45:37.5 - 02:46:07.5 57 1 12190+47182 1950 15 [4, 3] ## listobs 02:46:45.0 - 02:47:45.0 58 0 1219+4KC 4550 10 [0, 1] ## listobs 02:48:05.0 - 02:51:15.0 59 1 12190+47182 13000 10 [0, 1] ## listobs 02:51:25.0 - 02:52:15.0 60 2 12191+48299 3900 10 [0, 1] ## listobs 02:52:25.0 - 02:55:15.0 61 1 12190+47182 11700 10 [0, 1] ## listobs 02:55:25.0 - 02:56:15.0 62 2 12191+48299 3900 10 [0, 1] ## listobs 02:56:25.0 - 02:59:15.0 63 1 12190+47182 11700 10 [0, 1] ## listobs 02:59:25.0 - 03:00:15.0 64 2 12191+48299 3900 10 [0, 1] ## listobs 03:00:25.0 - 03:03:15.0 65 1 12190+47182 11700 10 [0, 1] ## listobs 03:03:25.0 - 03:04:15.0 66 2 12191+48299 3900 10 [0, 1] ## listobs 03:04:25.0 - 03:07:15.0 67 1 12190+47182 11700 10 [0, 1] ## listobs 03:07:25.0 - 03:08:15.0 68 2 12191+48299 3900 10 [0, 1] ## listobs 03:08:25.0 - 03:11:15.0 69 1 12190+47182 11700 10 [0, 1] ## listobs 03:11:25.0 - 03:12:15.0 70 2 12191+48299 3900 10 [0, 1] ## listobs 03:12:25.0 - 03:15:15.0 71 1 12190+47182 11700 10 [0, 1] ## listobs 03:15:25.0 - 03:16:15.0 72 2 12191+48299 3900 10 [0, 1] ## listobs 03:16:25.0 - 03:19:15.0 73 1 12190+47182 11700 10 [0, 1] ## listobs 03:19:25.0 - 03:20:15.0 74 2 12191+48299 3900 10 [0, 1] ## listobs 03:20:25.0 - 03:23:15.0 75 1 12190+47182 11700 10 [0, 1] ## listobs 03:23:25.0 - 03:24:05.0 76 2 12191+48299 3250 10 [0, 1] ## listobs 03:30:25.0 - 03:32:25.0 77 4 1331+305 8450 10 [0, 1] - flux cal ## listobs 03:36:15.0 - 03:38:15.0 78 5 1256-05 8450 10 [2, 3] - bandpass calibrator ## listobs 03:38:45.0 - 03:40:45.0 79 5 1256-05 8450 10 [4, 3] ## listobs (nVis = Total number of time/baseline visibilities per scan) ## listobs Fields: 6 ## listobs ID Code Name RA Decl Epoch SrcId nVis ## listobs 0 A 1219+4KC 12:19:06.4147 +48.29.56.1640 J2000 0 9750 ## listobs 1 12190+47182 12:18:57.5046 +47.18.14.3030 J2000 1 328250 - target ## listobs 2 A 12191+48299 12:19:06.4147 +48.29.56.1634 J2000 2 108550 - phasecal ## listobs 3 A 1219+48 12:19:06.4147 +48.29.56.1640 J2000 3 7800 ## listobs 4 A 1331+305 13:31:08.2879 +30.30.32.9580 J2000 4 8450 - flux cal ## listobs 5 B 1256-05 12:56:11.1665 -05.47.21.5240 J2000 5 16900 - bandpass cal ## listobs (nVis = Total number of time/baseline visibilities per field) ## listobs Spectral Windows: (5 unique spectral windows and 2 unique polarization setups) ## listobs SpwID #Chans Frame Ch1(MHz) ChanWid(kHz)TotBW(kHz) Ref(MHz) Corrs ## listobs 0 1 TOPO 22485.1 50000 50000 22485.1 RR RL LR LL ## listobs 1 1 TOPO 22435.1 50000 50000 22435.1 RR RL LR LL ## listobs 2 63 LSRK 22197.224 48.8301013 3076.29638 22198.7376 RR LL ## listobs 3 127 LSRK 22141.6008 12.2075253 1550.35572 22142.3697 RR LL ## listobs 4 127 LSRK 22137.5216 12.2075267 1550.3559 22138.2904 RR LL ## listobs Sources: 6 ## listobs ID Name SpwId RestFreq(MHz) SysVel(km/s) ## listobs 0 1219+4KC any 0 100 ## listobs 1 12190+47182 any 0 100 ## listobs 2 12191+48299 any 0 100 ## listobs 3 1219+48 any 22235.08 490 ## listobs 4 1331+305 any 0 100 ## listobs 5 1256-05 any 22235.08 484.6 ## listobs Antennas: 26: ## listobs ID Name Station Diam. Long. Lat. ## listobs 0 VA01 VLA:_N18 25.0 m -107.37.12.0 +33.54.58.3 ## listobs 1 VA02 VLA:_E1 25.0 m -107.37.05.7 +33.53.59.2 ## listobs 2 VA03 VLA:_E7 25.0 m -107.36.52.4 +33.53.56.5 ## listobs 3 VA04 VLA:_W2 25.0 m -107.37.07.4 +33.54.00.9 ## listobs 4 VA05 VLA:_N6 25.0 m -107.37.06.9 +33.54.10.3 ## listobs 5 VA06 VLA:_W8 25.0 m -107.37.21.6 +33.53.53.0 ## listobs 6 VA07 VLA:_N10 25.0 m -107.37.08.2 +33.54.22.4 ## listobs 7 VA08 VLA:_W1 25.0 m -107.37.05.9 +33.54.00.5 ## listobs 8 VA09 VLA:_N8 25.0 m -107.37.07.5 +33.54.15.8 ## listobs 9 VA10 VLA:_E8 25.0 m -107.36.48.9 +33.53.55.1 ## listobs 10 VA11 VLA:_W4 25.0 m -107.37.10.8 +33.53.59.1 ## listobs 11 VA12 VLA:_E2 25.0 m -107.37.04.4 +33.54.01.1 ## listobs 12 VA15 VLA:_E4 25.0 m -107.37.00.8 +33.53.59.7 ## listobs 13 VA16 VLA:_N14 25.0 m -107.37.09.9 +33.54.38.5 ## listobs 14 VA17 VLA:_E3 25.0 m -107.37.02.8 +33.54.00.5 ## listobs 15 VA18 VLA:_N4 25.0 m -107.37.06.5 +33.54.06.1 ## listobs 16 VA19 VLA:_E5 25.0 m -107.36.58.4 +33.53.58.8 ## listobs 17 VA20 VLA:_W7 25.0 m -107.37.18.4 +33.53.54.8 ## listobs 18 VA21 VLA:_W6 25.0 m -107.37.15.6 +33.53.56.4 ## listobs 19 VA22 VLA:_W9 25.0 m -107.37.25.1 +33.53.51.0 ## listobs 20 VA23 VLA:_W3 25.0 m -107.37.08.9 +33.54.00.1 ## listobs 21 VA24 VLA:_E6 25.0 m -107.36.55.6 +33.53.57.7 ## listobs 22 VA25 VLA:_N2 25.0 m -107.37.06.2 +33.54.03.5 ## listobs 23 VA26 VLA:_N16 25.0 m -107.37.10.9 +33.54.48.0 ## listobs 24 VA27 VLA:_N12 25.0 m -107.37.09.0 +33.54.30.0 ## listobs 25 VA28 VLA:_E9 25.0 m -107.36.45.1 +33.53.53.6 # selection mystep = 2 if(mystep in mysteps or doall): print "Step ", mystep, steptitle[mystep] os.system('rm -rf '+basename+'-spectral.ms*') split( vis = basename+'.ms', outputvis = basename+'-spectral.ms', datacolumn = 'data', field = '1,2,5', # target, phase cal, bandpass cal spw = '2:2~61,4:5~115' # remove edge channels right here ) # record flags flagmanager( vis = basename+'-spectral.ms', mode = 'save', versionname = 'original_data', merge = 'replace' ) listobs( vis = basename+'-spectral.ms' ) mystep = 3 if(mystep in mysteps or doall): print "Step ", mystep, steptitle[mystep] # restore original flags flagmanager( vis = basename+'-spectral.ms', mode = 'restore', versionname = 'original_data', merge = 'replace' ) flagdata( vis = basename+'-spectral.ms', mode = 'quack', unflag = False, quackinterval = 15.0, # == the integration time quackmode = 'beg', quackincrement = False ) # flagdata( # vis = basename+'-spectral.ms', # mode = 'manualflag', # selectdata = True, # correlation = 'RL,LR' # ) # plotms: field 1: timestamps 2:01:37.5 and 2:03:36.5, and 2:10:07.5 look low # => flag scans 28 and 29 and scans 32 and 33 # also flag scans 34,35,37, 57 because they are not flanked by phase calib observations # viewer: nothing obvious flagdata( vis = basename+'-spectral.ms', mode = 'manualflag', selectdata = True, timerange = '02:01:07.5~02:03:52.5, 02:09:37.5~02:18:52.5, 02:45:37.5~02:46:07.5' ) # see in viewer that calibrated data for all baselines involving ID 15 is high flagdata( vis = basename+'-spectral.ms', mode = 'manualflag', selectdata = True, antenna = '0015' ) plotms( vis = basename+'-spectral.ms', xaxis = 'time', yaxis = 'amp', selectdata = True, field = '1' ) user_check = raw_input('Flag if you like, then press return to continue ... or ctrl-c + return to stop') # calibration # a) set flux scale mystep = 4 if(mystep in mysteps or doall): print "Step ", mystep, steptitle[mystep] clearcal( vis = basename+'-spectral.ms' ) setjy( vis = basename+'-spectral.ms', field = '1', # phase cal spw = '', # all spws modimage = '', # assume point source, models are in os.environ["CASAPATH"].split()[0]+'/data/nrao/VLA/CalModels/3C286_K.im', fluxdensity = [0.5968,0.,0.,0.], # take from continuum analysis of phase calib (from the K-band continuum data in this set) standard = 'Perley-Taylor 99' ) # b) bandpass mystep = 5 if(mystep in mysteps or doall): print "Step ", mystep, steptitle[mystep] os.system('rm -rf '+basename+'-spectral.bcal') # first try B solution, i.e. independent solution for each channel usebpoly = True if(not usebpoly): bandpass( vis = basename+'-spectral.ms', caltable = basename+'-spectral.bcal', field = '2', # bandpass cal spw = '', # all available, i.e. 0~1 opacity = 0.0, solint = 'inf', combine = '', # i.e. only combine scans, the default refant = 'VA04', minblperant = 4, # the default solnorm = True, # should only be True when bandpass is first step of cal bandtype = 'B', gaintable = '', # no previous gain table gainfield = '', # if there was a previous gain table, use all fields gaincurve = True ) # bandpass Found good B Jones solutions in 4 slots. else: # use bpoly bandpass( vis = basename+'-spectral.ms', caltable = basename+'-spectral.bcal', field = '2', # bandpass cal spw = '0', opacity = 0.0, solint = 'inf', combine = '', # i.e. only combine scans, the default refant = 'VA04', minblperant = 4, # the default solnorm = True, # should only be True when bandpass is first step of cal bandtype = 'BPOLY', degamp = 7, # degree of the polynomial for amp solution, default = 3 degphase = 7, # degree of the polynomial for phase solution, default = 3 maskedge = 5, # ignore edge channels, default = 5 (%) gaintable = '', # no previous gain table gainfield = '', # if there was a previous gain table, use all fields gaincurve = True ) bandpass( vis = basename+'-spectral.ms', caltable = basename+'-spectral.bcal', field = '2', # bandpass cal spw = '1', opacity = 0.0, solint = 'inf', combine = '', # i.e. only combine scans, the default refant = 'VA04', minblperant = 4, # the default solnorm = True, # should only be True when bandpass is first step of cal bandtype = 'BPOLY', degamp = 7, # degree of the polynomial for amp solution, default = 3 degphase = 7, # degree of the polynomial for phase solution, default = 3 maskedge = 5, # ignore edge channels, default = 5 (%) gaintable = '', # no previous gain table gainfield = '', # if there was a previous gain table, use all fields gaincurve = True, append = True ) # end if plotcal( caltable = basename+'-spectral.bcal', xaxis = 'chan', yaxis = 'amp', iteration = 'antenna, spw' ) user_check = raw_input('Quit plotcal, then press return.') plotcal( caltable = basename+'-spectral.bcal', xaxis = 'chan', yaxis = 'phase', iteration = 'antenna, spw' ) user_check = raw_input('Quit plotcal, then press return or ctrl-c + return to stop.') # gain calibration mystep = 6 if(mystep in mysteps or doall): print "Step ", mystep, steptitle[mystep] os.system('rm -rf '+basename+'.gcal') gaincal( vis = basename+'-spectral.ms', caltable = basename+'-spectral.gcal', field = '1', # phase calibrator spw = '', calmode = 'ap', # amplitude and phase solution (default) gaintype = 'G', solint = 'inf', combine = '', refant = 'VA04', # name, not index minsnr = 0.5, # minimum signal to noise (reduced to 0.3 to accomodate some time ranges in spw 1) gaincurve = True, # internal VLA gain curve correction gaintable = basename+'-spectral.bcal', # preapply the bandpass table from above opacity = 0.0 ) plotcal( caltable = basename+'-spectral.gcal', xaxis = 'time', yaxis = 'amp', iteration = 'antenna' ) user_check = raw_input('Quit plotcal, then press return ...') plotcal( caltable = basename+'-spectral.gcal', xaxis = 'time', yaxis = 'phase', iteration = 'antenna' ) user_check = raw_input('Quit plotcal, then press return ... or ctrl-c + return to stop') mystep = 7 if(mystep in mysteps or doall): print "Step ", mystep, steptitle[mystep] applycal( vis = basename+'-spectral.ms', gaintable = [basename+'-spectral.gcal',basename+'-spectral.bcal'], field = '0', # the target gaincurve = True, opacity = 0.0, spwmap = [], # one-to-one, the default ) # image the low-res spw mystep = 8 if(mystep in mysteps or doall): print "Step ", mystep, steptitle[mystep] os.system('rm -rf '+basename+'-spectral-target*') clean( vis=basename+'-spectral.ms', imagename=basename+'-spectral-target-lores', field='0', spw='0', cell=[0.5,0.5],imsize=[128,128], stokes='I', mode='frequency', nchan=-1,start='',width='', # the edge channels were already removed outframe='', psfmode='clark', imagermode='csclean', scaletype='SAULT', niter=100, gain = 0.1, restoringbeam = ['2.1arcsec'], threshold='0.5mJy', restfreq='22235.0MHz', # the rest freq of the H2O maser line phasecenter='', weighting='natural', interactive=F, minpb=0.1, pbcor=F ) print 'Created image ', basename+'-spectral-target-lores.image' # image the hires spw mystep = 9 if(mystep in mysteps or doall): print "Step ", mystep, steptitle[mystep] os.system('rm -rf '+basename+'-spectral-target-hires*') clean( vis=basename+'-spectral.ms', imagename=basename+'-spectral-target-hires', field='0', spw='1', cell=[0.5,0.5],imsize=[128,128], stokes='I', mode='frequency', nchan=-1,start='',width='', outframe='', psfmode='clark',imagermode='csclean', scaletype='SAULT', niter=500, gain = 0.1, restoringbeam = ['2.1arcsec'], threshold='0.5mJy', restfreq='22235.0MHz', # the rest freq of the H2O maser line phasecenter='', weighting='natural', interactive=F, minpb=0.1, pbcor=F ) print 'Created image ', basename+'-spectral-target-hires.image'