#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! #!! This script describes SMA data reduction using miriad to fill and #!! apply Tsys correction, and then writing uvfits. CASA is then used #!! to read uvfits, flag, calibrate, continuum subtract, and image the #!! data. Plotting for datasets with many spw such as SMA data !! #!! (plotxy) can be slow, where possibel the new plotter plotms has been #!! used. Please be patient, CASA is a work in progress. Although this #!! script can be run, it is probably more valuable to cut and paste and #!! read all the detailed comments. If you have difficulty with either #!! plotxy or plotms starting, type clearstat(). #!! This script requires Beta patch 3.0.0 or later. Last edited #!! 12/17/09. # # JO 2010-03-13 Release 3.0.1 changed uvcontsub2 and fallout # #!! The data themselves are a six-pointing mosaic of the infrared dark #!! cloud G19.3+0.07 with 12CO(2-1) in the USB from Brogan et al. in #!! prep. Please do not use these data for scientific purposes. # #!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ####################################################################### ####################################################################### # DATA LOADING, TSYS CORRECTION, AND UVFITS WRITING IN MIRIAD ####################################################################### ####################################################################### # # smalod in='../050819_04.23.28' out=g19_d2 sideband=2 rxif=0 # # smafix vis=g19_d2_rx0.usb out=g19_d2_rx0.usb.tsys device=/xs xaxis=time # \ yaxis=systemp nxy=3,3 options=tsyscor # # fits in=g19_d2_rx0.usb.tsys out=g19_d2usb.spw1 op=uvout select=win\(1\) # fits in=g19_d2_rx0.usb.tsys out=g19_d2usb.spw2 op=uvout select=win\(2\) # ... # fits in=g19_d2_rx0.usb.tsys out=g19_d2usb.spw24 op=uvout select=win\(24\) # #### Note that currently these really do have to be written out #### individually (of course you could write a little miriad script to #### do it). If not the spw are concatenated together into a single #### spw which is inappropriate before calibration. # # mkdir spw_fits_files # # mv g19_d2usb.spw* spw_fits_files/ # # the file spw_fits_files.tar is what is distributed as the "data" # ####################################################################### ####################################################################### #### INSIDE CASAPY ####################################################################### ####################################################################### os.system('rm -rf spw_fits.files') os.system('tar xvf spw_fits_files.tar') os.system('rm -rf g19_d2usb*') #### Read in each uvfits file (need to put upper range to one more #### than have, to satisfy ipython peculiarity). The myfiles makes a #### list of the filenames that can be used in a later call. os.system('rm -rf spw_ms_files') os.system('mkdir spw_ms_files') myfiles = [] for i in range(1,25): msfile = "spw_ms_files/g19_d2usb.spw"+str(i)+".ms" importuvfits(fitsfile="spw_fits_files/g19_d2usb.spw"+str(i), vis=msfile) myfiles.append(msfile) concat(vis=myfiles,concatvis='g19_d2usb.ms',timesort=True) #### Initialize the scratch coloumns clearcal(vis='g19_d2usb.ms') #### Get a listing of the properties of the concatenated files. listobs(vis='g19_d2usb.ms',verbose=False) # 0 jupiter 13:01:25.48 -05.18.53.00 J2000 # 1 1743-038 17:43:58.85 -03.50.04.62 J2000 # 2 g19.3a 18:25:58.70 -12.03.57.80 J2000 # 3 g19.3b 18:25:57.40 -12.04.18.50 J2000 # 4 g19.3c 18:25:56.09 -12.04.28.20 J2000 # 5 g19.3d 18:25:54.60 -12.04.23.10 J2000 # 6 1908-201 19:11:09.65 -20.06.55.11 J2000 # 7 g19.3e 18:25:54.80 -12.04.43.80 J2000 # 8 g19.3f 18:25:53.30 -12.04.51.00 J2000 # 9 3c454.3 22:53:57.74 +16.08.53.56 J2000 # 10 0359+509 03:59:29.75 +50.57.50.15 J2000 # 11 hip19167 04:06:35.04 +50.21.04.54 J2000 # 12 polaris 02:31:48.70 +89.15.50.72 J2000 # 13 uranus 22:44:40.51 -08.50.23.77 J2000 # # Spectral Windows: (24 unique spectral windows and 1 unique polarization setups) # SpwID #Chans Frame Ch1(MHz) ChanWid(kHz)TotBW(kHz) Ref(MHz) Corrs # 0 128 LSRK 229347.598 812.5 104000 229347.598 XX # 1 128 LSRK 229429.601 812.5 104000 229429.601 XX # 2 128 LSRK 229504.79 812.5 104000 229504.79 XX # 3 128 LSRK 229586.792 812.5 104000 229586.792 XX # 4 128 LSRK 229675.607 812.5 104000 229675.607 XX # 5 128 LSRK 229757.609 812.5 104000 229757.609 XX # 6 128 LSRK 229832.798 812.5 104000 229832.798 XX # 7 128 LSRK 229914.8 812.5 104000 229914.8 XX # 8 128 LSRK 230003.615 812.5 104000 230003.615 XX # 9 128 LSRK 230085.617 812.5 104000 230085.617 XX # 10 128 LSRK 230160.806 812.5 104000 230160.806 XX # 11 128 LSRK 230242.808 812.5 104000 230242.808 XX # 12 128 LSRK 230331.623 812.5 104000 230331.623 XX # 13 128 LSRK 230413.625 812.5 104000 230413.625 XX # 14 128 LSRK 230488.814 812.5 104000 230488.814 XX # 15 128 LSRK 230570.817 812.5 104000 230570.817 XX # 16 128 LSRK 230659.631 812.5 104000 230659.631 XX # 17 128 LSRK 230741.634 812.5 104000 230741.634 XX # 18 128 LSRK 230816.823 812.5 104000 230816.823 XX # 19 128 LSRK 230898.825 812.5 104000 230898.825 XX # 20 128 LSRK 230987.639 812.5 104000 230987.639 XX # 21 128 LSRK 231069.641 812.5 104000 231069.641 XX # 22 128 LSRK 231144.831 812.5 104000 231144.831 XX # 23 128 LSRK 231226.833 812.5 104000 231226.833 XX #### !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! #### Brief description of the purposes of the sources in this project; #### () gives source id number. #### Bandpass calibrator: 3c454.3 (9) # #### Gain calibrators: 1743-038 (1), 1908-201 (6) # #### Flux calibrators: 3c454.3 (9) and uranus (13) # #### Science target g19: (2,3,4,5,7,8) six-pointing mosaic # #### other sources were not part of this program #### !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ####################################################################### #### Initial inspection and flagging ####################################################################### # Here we average channels, and spw and clearstat() plotms(vis="g19_d2usb.ms",xaxis="time",yaxis="amp",field='1~9,13', avgchannel='3072') # Now you can use the locate features in plotms to find the bad # data points on field 8 and antenna ids 1-2 and 0-2 and either # flag in plotms or for script purposes write a flagdata command. # Pause script if you are running scriptmode # NOTE: if you are running interactively and cutting/pasting you do not # need to include these lines that come after plotting user_check=raw_input('Return to continue script\n') #### NOTE: locate gives antenna ids, but flag data looks first at the #### name column, since SMA antenna names are numbers this can be #### confusing -- the corresponding NAMES are antennas 2-3 and 1-3, #### which you can check with listobs. flagdata(vis='g19_d2usb.ms',mode='manualflag',field='8',spw='20~23', timerange='07:52:45~07:52:47',antenna='2&3;1&3') clearstat() plotms(vis="g19_d2usb.ms",xaxis="time",yaxis="amp",field='8', avgchannel='3072') # Pause script if you are running scriptmode user_check=raw_input('Return to continue script\n') # Note to confirm the data were flagged in plotms, you need to change # selection in plot (for example change field and then change back), # and then change back to original to clear the cache and force a # re-read of the data. #### Make an inital amplitude vs. channel plot of the bandpass #### calibrator 3c454.3 averaging all the data together in #### time. crossscan=T needed because there are two scans on #### 3c454.3. Note the rolloff of 24 individual spw. Try zooming to #### see better. #### NOTE: after you run plotxy, run clearstat() before trying to run #### plotms or you will get a table lock. clearstat() plotxy(vis="g19_d2usb.ms",xaxis="channel",yaxis="amp",datacolumn="data", iteration="baseline", antenna="",spw="",field="9", averagemode="vector",width="1",timebin="all", crossscans=T,crossbls=False,stackspw=F, subplot=331,plotcolor="darkcyan", overplot=False,showflags=False,interactive=True,figfile="") # Pause script if you are running scriptmode user_check=raw_input('Return to continue script\n') #### As you can see in this plot, SMA spw overlap by a fair amount, #### flag first 7 channels on either side of each spw to aid #### calibration. default('flagdata') flagdata(vis='g19_d2usb.ms', mode='manualflag',spw='0~23:0~6;121~127') #### On many baselines the first integration of every scan is #### low. These can be flagged with flagdata. The integration time of #### these data is 30.9s. Note, more recent SMA data rarely shows this #### problem. flagdata(vis="g19_d2usb.ms",mode="quack",field='', spw="0~23",quackinterval=40.0) #### The following can be used to plot the antenna positions, note #### that antenna 8 was in the barn for repair. clearstat() plotxy(vis="g19_d2usb.ms",xaxis="x") # Pause script if you are running scriptmode user_check=raw_input('Return to continue script\n') ####################################################################### #### BANDPASS CALIBRATION ####################################################################### #### Next plot phase vs. channel. Note the large descrepancy in phase #### for spw=4-7 on some baselines, this will be important later. clearstat() plotxy(vis="g19_d2usb.ms",xaxis="channel",yaxis="phase",datacolumn="data", iteration="baseline", antenna="",spw="",field="9", averagemode="vector",width="1",timebin="all", crossscans=T,crossbls=False,stackspw=F, subplot=331,plotcolor="darkcyan", overplot=False,showflags=False,interactive=True,figfile="") # Pause script if you are running in scriptmode user_check=raw_input('Return to continue script\n') #### Now plot amplitude and phase as a function of time. clearstat() plotms(vis="g19_d2usb.ms",xaxis="time",yaxis="amp",field='9', avgchannel='3072') # Pause script if you are running scriptmode user_check=raw_input('Return to continue script\n') #### The gap in 3C454.3 observation is due to a pause to determine #### pointing solutions. The flags below catch times just before and #### after when the amplitudes and phase are off. This flagging could #### also have been accomplished interactively. default('flagdata') flagdata(vis="g19_d2usb.ms",field="9",timerange='12:16:58~12:19:00') flagdata(vis="g19_d2usb.ms",field="9",timerange='12:08:00~12:09:30') #### Here because msplot cannot do multipanel iteration yet we use the #### slower plotxy clearstat() plotxy(vis="g19_d2usb.ms",xaxis="time",yaxis="phase",datacolumn="data", iteration="baseline", antenna="",spw="0~23",field="9", averagemode="vector",width="allspw",timebin="0", subplot=331,plotsymbol=".", multicolor="none",plotcolor="darkcyan", overplot=False,showflags=False,interactive=True,figfile="") # Pause script if you are running in scriptmode user_check=raw_input('Return to continue script\n') #### Start with an initial phase only solution for 3c454.3 to take out #### major phase variations with time. For maximum sensitivity we #### average all the spw together with the exception of spw 4-7 which #### show large phase offsets from the other spw in phase vs. channel. gaincal(vis='g19_d2usb.ms', caltable='g19_d2usb.ms.pcal', field='9', spw='0~3,8~23', gaintype='G', minsnr=2.0, refant='3', calmode='p',solint='int', combine='spw') clearstat() plotcal(caltable="g19_d2usb.ms.pcal",xaxis="time",yaxis="phase", field="9",antenna="",spw="",timerange="",subplot=331, overplot=False,clearpanel="Auto",iteration="antennas", plotsymbol="o",plotcolor="blue",showgui=True,figfile="") # Pause script if you are running in scriptmode user_check=raw_input('Return to continue script\n') #### A regular 'B' type solution can work well if you have lots of #### S/N on your bandpass calibrator. You may wish to make this table #### (uncomment commands below) and compare to 'BPOLY' solution #### further down. # bandpass(vis='g19_d2usb.ms', caltable='g19_d2usb.ms.bpcal', # field='9', spw='0~23', bandtype='B', solint='inf', combine='scan', # refant='3', gaincurve=False, opacity=0.0, # gaintable='g19_d2usb.ms.pcal',spwmap=[0]) # #plotcal(caltable="g19_d2usb.ms.bpcal",xaxis="freq",yaxis="amp", # field="9",antenna="",spw="",timerange="",subplot=311, # overplot=False,clearpanel="Auto",iteration="antenna", # plotsymbol="o",plotcolor="blue",showgui=True,figfile="") #### A polynomial solution is more sensitive than simply fitting a #### polynomial to the channel based solutions 'B' after the fact. For #### the mm/submm where the S/N is often low on the bandpass calibrator #### 'BPOLY' is often the best choice, but every dataset is different. #### You may especially need to play with degamp and degphase. The #### combine='scan' is needed to average together the two 3c454.3 scans. #### The spwmap=[0] is telling it to apply the phase solution that was #### found from the average of all spw to each spw before determining #### the bandpass solution. bandpass(vis='g19_d2usb.ms', caltable='g19_d2usb.ms.bpoly7', field='9', spw='0~23', bandtype='BPOLY', solint='inf', combine='scan', refant='3', degamp=7, degphase=7, gaintable='g19_d2usb.ms.pcal',spwmap=[0]) clearstat() plotcal(caltable="g19_d2usb.ms.bpoly7",xaxis="freq",yaxis="amp", field="9",antenna="",spw="",timerange="",subplot=311, overplot=False,clearpanel="Auto",iteration="antenna", plotsymbol="o",plotcolor="blue",showgui=True,figfile="") # Pause script if you are running in scriptmode user_check=raw_input('Return to continue script\n') #### Unfortunately, this kind of plot with many spw and colors #### takes a long time. Be patient for now its being worked on. You #### can safely ignore where the solution gets large at edges of #### spw -- its showing where we flagged channels -- this is #### being worked on. Zoom in if necessary to see that inner #### portions are nice and flat. #### The next step is to apply the bandpass solution while #### re-calculating the phase only solution for 3c454.3 on a short #### (integration time) time interval, but now including all spw #### (including 4~7 since antenna based solutions will have taken out #### phase vs.spw differences). gaincal(vis='g19_d2usb.ms', caltable='g19_d2usb.ms.pcal2', field='9', spw='0~23', gaintype='G', minsnr=2.0, refant='3', calmode='p',solint='int', combine='spw', gaintable=['g19_d2usb.ms.bpoly7']) clearstat() plotcal(caltable="g19_d2usb.ms.pcal2",xaxis="time",yaxis="phase", field="9",antenna="",spw="",timerange="",subplot=331, overplot=False,clearpanel="Auto",iteration="antennas", plotsymbol="o",plotcolor="blue",showgui=True,figfile="") # Pause script if you are running in scriptmode user_check=raw_input('Return to continue script\n') #### Next we apply the antenna based bandpass, and phase only solution #### to solve for amplitude solutions over a longer solint (to #### increase S/N).The spwmap=[[0],[]] tells it that the pcal solution #### is independent of spw while the bpoly7 solution is per spw, #### spwmap order must be in the same order as the tables in #### gaintable. gaincal(vis='g19_d2usb.ms', caltable='g19_d2usb.ms.apcal', field='9', spw='0~23', gaintype='G', minsnr=2.0, refant='3', calmode='a',solint='300s', combine='spw', gaintable=['g19_d2usb.ms.pcal2','g19_d2usb.ms.bpoly7'], spwmap=[[0],[]]) clearstat() plotcal(caltable="g19_d2usb.ms.apcal",xaxis="time",yaxis="amp", field="9",antenna="",spw="",timerange="",subplot=331, overplot=False,clearpanel="Auto",iteration="antennas", plotsymbol="o",plotcolor="blue",showgui=True,figfile="") # Pause script if you are running in scriptmode user_check=raw_input('Return to continue script\n') #### Although in previous steps solutions were applied "on-the-fly", #### the data were not actually changed. Running applycal will fill #### the corrected data column in the ms. Again be careful with #### specification of gaintable, and keeping the same order in spwmap #### and gainfield. NOTE this applycal is not necessary for later #### steps -- it is being done in order to allow plotting of data #### calibration to this point. applycal(vis='g19_d2usb.ms',spw='0~23', field='9', gaintable=['g19_d2usb.ms.pcal2','g19_d2usb.ms.apcal', 'g19_d2usb.ms.bpoly7'], spwmap=[[0],[0],[]],gainfield=['9','9','9']) clearstat() plotxy(vis="g19_d2usb.ms",xaxis="channel",yaxis="amp", datacolumn="corrected",iteration="baseline", antenna='',spw="0~23",field="9", averagemode="vector",width="1",timebin="all",crossscans=T, crossbls=False,stackspw=False,extendflag=F, extendchan="",extendspw="",extendant="",extendtime="", plotcolor='darkcyan',subplot=331, multicolor="none",figfile="") # Pause script if you are running in scriptmode user_check=raw_input('Return to continue script\n') #### As this plotxy shows, after antenna based bandpass, SMA data #### often show residual baseline based problems as a function of spw #### (i.e. a few abnormally low amplitudes for a few spw on a few #### baselines). To fix this problem we need to do a baseline based #### amplitude only solution per spw, but we'll need to remove an #### overall normalization factor first by combining all spw (this #### option will be incorporated into blcal one day so that only one #### run is necessary). blcal(vis='g19_d2usb.ms', caltable='g19_d2usb.ms.blcal1', field='9', spw='0~23',solint='inf', combine='spw,scan', gaintable=['g19_d2usb.ms.pcal2','g19_d2usb.ms.apcal', 'g19_d2usb.ms.bpoly7'], calmode='a',spwmap=[[0],[0],[]]) blcal(vis='g19_d2usb.ms', caltable='g19_d2usb.ms.blcal2', field='9', spw='0~23',solint='inf', combine='scan', gaintable=['g19_d2usb.ms.pcal2','g19_d2usb.ms.apcal', 'g19_d2usb.ms.blcal1','g19_d2usb.ms.bpoly7',], calmode='a',spwmap=[[0],[0],[0],[]]) clearstat() plotcal(caltable="g19_d2usb.ms.blcal1",xaxis="freq",yaxis="amp", field="9",antenna="",spw="",timerange="",subplot=511, overplot=False,clearpanel="Auto",iteration="antenna", plotsymbol="o",plotcolor="blue",showgui=True,figfile="") # Pause script if you are running in scriptmode user_check=raw_input('Return to continue script\n') #### Now we can do a new applycal including ONLY the second blcal2 #### solution. Note that applycal overwrites the corrected column so #### there is no need to clearcal. Like the antenna based bandpass #### the blcal solution is independent for each spw, so its spwmap #### entry is []. applycal(vis='g19_d2usb.ms',spw='0~23', field='9', gaintable=['g19_d2usb.ms.pcal2','g19_d2usb.ms.apcal', 'g19_d2usb.ms.bpoly7','g19_d2usb.ms.blcal2'], spwmap=[[0],[0],[],[]],gainfield=['9','9','9','9']) clearstat() plotxy(vis="g19_d2usb.ms",xaxis="channel",yaxis="amp", datacolumn="corrected",iteration="baseline", antenna='',spw="0~23",field="9", averagemode="vector",width="1",timebin="all",crossscans=T, crossbls=False,stackspw=False,extendflag=F, extendchan="",extendspw="",extendant="",extendtime="", plotcolor='darkcyan',subplot=331, multicolor="none",figfile="") # Pause script if you are running in scriptmode user_check=raw_input('Return to continue script\n') #### Note how flat the full bandpass is after all this calibration ####################################################################### #### INSPECTION OF GAIN CALIBRATORS AND SET FLUX SCALE ####################################################################### #### Apply antenna and baseline based bandpass solutions to the #### gain calibrators clearstat() applycal(vis='g19_d2usb.ms',spw='0~23', field='1,6', gaintable=['g19_d2usb.ms.bpoly7','g19_d2usb.ms.blcal2'], spwmap=[[],[]],gainfield=['9','9']) #### Have a look at both calibrators with channel vs. amp. clearstat() plotxy(vis="g19_d2usb.ms",xaxis="channel",yaxis="amp",datacolumn="corrected", iteration="baseline", antenna="",spw="",field="1", averagemode="vector",width="1",timebin="all", subplot=331,plotsymbol=".", multicolor="none",plotcolor="darkcyan", overplot=False,showflags=False,interactive=True,figfile="") # Pause script if you are running in scriptmode user_check=raw_input('Return to continue script\n') clearstat() plotxy(vis="g19_d2usb.ms",xaxis="channel",yaxis="amp",datacolumn="corrected", iteration="baseline", antenna="",spw="",field="6", averagemode="vector",width="1",timebin="all", subplot=331,plotsymbol=".", multicolor="none",plotcolor="darkcyan", overplot=False,showflags=False,interactive=True,figfile="") # Pause script if you are running in scriptmode user_check=raw_input('Return to continue script\n') #### Have a look at both calibrators in phase and amplitude #### vs. time. For the moment, unfortunately multicolor='field' #### doesn't work, so easiest to do seperately. clearstat() plotxy(vis="g19_d2usb.ms",xaxis="time",yaxis="amp",datacolumn="corrected", iteration="baseline", antenna="",spw="0~23",field="1", averagemode="vector",width="allspw",timebin="0", subplot=331,plotsymbol=".", multicolor="none",plotcolor="darkcyan", overplot=False,showflags=False,interactive=True,figfile="") # Pause script if you are running in scriptmode user_check=raw_input('Return to continue script\n') clearstat() plotxy(vis="g19_d2usb.ms",xaxis="time",yaxis="amp",datacolumn="corrected", iteration="baseline", antenna="",spw="0~23",field="6", averagemode="vector",width="allspw",timebin="0", subplot=331,plotsymbol=".", multicolor="none",plotcolor="red", overplot=F,showflags=False,interactive=True,figfile="") # Pause script if you are running in scriptmode user_check=raw_input('Return to continue script\n') clearstat() plotxy(vis="g19_d2usb.ms",xaxis="time",yaxis="phase",datacolumn="corrected", iteration="baseline", antenna="",spw="0~23",field="1", averagemode="vector",width="allspw",timebin="0", subplot=331,plotsymbol=".", multicolor="none",plotcolor="darkcyan", overplot=False,showflags=False,interactive=True,figfile="") # Pause script if you are running in scriptmode user_check=raw_input('Return to continue script\n') clearstat() plotxy(vis="g19_d2usb.ms",xaxis="time",yaxis="phase",datacolumn="corrected", iteration="baseline", antenna="",spw="0~23",field="6", averagemode="vector",width="allspw",timebin="0", subplot=331,plotsymbol=".", multicolor="none",plotcolor="darkcyan", overplot=False,showflags=False,interactive=True,figfile="") # Pause script if you are running in scriptmode user_check=raw_input('Return to continue script\n') #### At the time of these observations, the calibrator J1911 had the #### most recent and relatively stable monitor data (SMA website), so #### we will use it for absolute flux calibration. In future tools for #### using resolved planets will be implemented. #### 1911-201 #### 1mm 17 Jul 2005 11:38 SMA 274.5 1.88 +/- 0.13 mgurwell #### 1mm 05 Sep 2005 08:10 SMA 221.4 1.99 +/- 0.12 mgurwell setjy(vis='g19_d2usb.ms', field='6', spw='0~23', fluxdensity=[2.0,0.,0.,0.]) ####################################################################### #### FINAL GAIN AND FLUX CALIBRATION FOR ALL CALIBRATORS ####################################################################### #### To avoid decorrelation of the amplitude, the phase solutions fed #### to the amplitude calibration step should be at the shortest #### solint your S/N can support, typically the integration time. gaincal(vis='g19_d2usb.ms', caltable='g19_d2usb.ms.allpcal', field='1,6,9', spw='0~23', gaintype='G', minsnr=2.0, refant='3', calmode='p',solint='int', combine='spw', gaintable=['g19_d2usb.ms.bpoly7','g19_d2usb.ms.blcal2'], spwmap=[[],[]]) #### Unless you chose to do some more sophisticated smoothing, the #### phase solutions applied to your target might as well be on the #### scan time of your phase calibrator(s). Without combine='scan', #### the solution will not cross scan boundaries even though #### solint='inf'. This is an easy way to "manually" interpolate. gaincal(vis='g19_d2usb.ms', caltable='g19_d2usb.ms.allpcalscan', field='1,6,9', spw='0~23', gaintype='G', minsnr=2.0, refant='3', calmode='p',solint='inf', combine='spw', gaintable=['g19_d2usb.ms.bpoly7','g19_d2usb.ms.blcal2'], spwmap=[[],[]]) #### Plot the phase solutions of the phase calibrators. Can repeat for #### both short and long solution intervals. clearstat() plotcal(caltable="g19_d2usb.ms.allpcalscan",xaxis="time",yaxis="phase", field="1,6",antenna="",spw="",timerange="",subplot=331, overplot=False,clearpanel="Auto",iteration="antennas", plotsymbol="o",plotcolor="blue",showgui=True,figfile="") # Pause script if you are running in scriptmode user_check=raw_input('Return to continue script\n') #### Amplitude calibration with longer solint and applying #### short solint phase solution. gaincal(vis='g19_d2usb.ms', caltable='g19_d2usb.ms.allapcal', field='1,6,9', spw='0~23', gaintype='G', minsnr=2.0, refant='3', calmode='ap',solint='300s', combine='spw', gaintable=['g19_d2usb.ms.allpcal','g19_d2usb.ms.bpoly7', 'g19_d2usb.ms.blcal2'],spwmap=[[0],[],[]]) #### Plot the amplitude solutions clearstat() plotcal(caltable="g19_d2usb.ms.allapcal",xaxis="time",yaxis="amp", field="1,6",antenna="",spw="",timerange="",subplot=331, overplot=False,clearpanel="Auto",iteration="antennas", plotsymbol="o",plotcolor="blue",showgui=True,figfile="") # Pause script if you are running in scriptmode user_check=raw_input('Return to continue script\n') #### Derive absolute flux calibration (based on setjy above) fluxscale(vis='g19_d2usb.ms',caltable='g19_d2usb.ms.allapcal', fluxtable='g19_d2usb.ms.fluxcal',reference='6') ####################################################################### #### APPLY FINAL CALIBRATIONS ####################################################################### #### Final application of calibration for each source is done #### independently below. Flux densities from SMA flux monitoring are #### provided as checks of the absolute flux calibration. applycal(vis='g19_d2usb.ms',spw='0~23', field='9', gaintable=['g19_d2usb.ms.allpcal','g19_d2usb.ms.fluxcal', 'g19_d2usb.ms.bpoly7','g19_d2usb.ms.blcal2'],spwmap=[[0],[0],[],[]], gainfield=['9','9','9','9']) clearstat() plotms(vis="g19_d2usb.ms",xaxis="time",yaxis="amp",field='9', avgchannel='3072',ydatacolumn='corrected') # Pause script if you are running scriptmode user_check=raw_input('Return to continue script\n') #### 3C454.3 #### 1mm 19 Aug 2005 16:01 SMA 225.3 27.36 +/- 1.37 mgurwell applycal(vis='g19_d2usb.ms',spw='0~23', field='1', gaintable=['g19_d2usb.ms.allpcal','g19_d2usb.ms.fluxcal', 'g19_d2usb.ms.bpoly7','g19_d2usb.ms.blcal2'],spwmap=[[0],[0],[],[]], gainfield=['1','1','9','9']) clearstat() plotms(vis="g19_d2usb.ms",xaxis="time",yaxis="amp",field='1', avgchannel='3072',ydatacolumn='corrected') # Pause script if you are running scriptmode user_check=raw_input('Return to continue script\n') #### J1743-038 #### 1mm 26 Aug 2005 06:05 SMA 225.6 1.71 +/- 0.09 mgurwell applycal(vis='g19_d2usb.ms',spw='0~23', field='6', gaintable=['g19_d2usb.ms.allpcal','g19_d2usb.ms.fluxcal', 'g19_d2usb.ms.bpoly7','g19_d2usb.ms.blcal2'],spwmap=[[0],[0],[],[]], gainfield=['6','6','9','9']) clearstat() plotms(vis="g19_d2usb.ms",xaxis="time",yaxis="amp",field='6', avgchannel='3072',ydatacolumn='corrected') # Pause script if you are running in scriptmode user_check=raw_input('Return to continue script\n') #### J1911-201 #### 1mm 17 Jul 2005 11:38 SMA 274.5 1.88 +/- 0.13 mgurwell #### 1mm 05 Sep 2005 08:10 SMA 221.4 1.99 +/- 0.12 mgurwell #### Since Uranus is resolved we apply the solutions from 3C454.3 to #### achieve at lest some reduction in instrumental gain variations. applycal(vis='g19_d2usb.ms',spw='0~23', field='13', gaintable=['g19_d2usb.ms.allpcal','g19_d2usb.ms.fluxcal', 'g19_d2usb.ms.bpoly7','g19_d2usb.ms.blcal2'],spwmap=[[0],[0],[],[]], gainfield=['9','9','9','9']) clearstat() plotms(vis="g19_d2usb.ms",xaxis="uvdist",yaxis="amp",field='13', avgchannel='3072',ydatacolumn='corrected') # Pause script if you are running in scriptmode user_check=raw_input('Return to continue script\n') #### A UV-distance plot of the expected flux of Uranus can be obtained #### from the SMA webpage. In future we will implement tools to use #### planets for absolute flux calibration. applycal(vis='g19_d2usb.ms',spw='0~23', field='2,3,4,5,7,8', gaintable=['g19_d2usb.ms.allpcalscan','g19_d2usb.ms.fluxcal', 'g19_d2usb.ms.bpoly7','g19_d2usb.ms.blcal2'],spwmap=[[0],[0],[],[]], gainfield=['1,6','1,6','9','9']) ################################################################### #### SPLIT OFF TARGET DATA FOR FURTHER INSPECTION ################################################################### split(vis="g19_d2usb.ms",outputvis='g19_d2usb_targets.ms', field="2,3,4,5,7,8") listobs(vis="g19_d2usb_targets.ms",verbose=F) #### ID Code Name Right Ascension Declination Epoch #### 0 g19.3a 18:25:58.70 -12.03.57.80 J2000 #### 1 g19.3b 18:25:57.40 -12.04.18.50 J2000 #### 2 g19.3c 18:25:56.09 -12.04.28.20 J2000 #### 3 g19.3d 18:25:54.60 -12.04.23.10 J2000 #### 4 g19.3e 18:25:54.80 -12.04.43.80 J2000 #### 5 g19.3f 18:25:53.30 -12.04.51.00 J2000 #### Note that after splitting, the "corrected" data column in the #### original vis file will become the "data" and "corrected" columns #### and the field IDs are renumbered. clearstat() plotms(vis='g19_d2usb_targets.ms',xaxis="frequency",yaxis="amp",field='0', avgtime='1e7',avgscan=True,avgbaseline=True) # Pause script if you are running in scriptmode user_check=raw_input('Return to continue script\n') clearstat() plotms(vis='g19_d2usb_targets.ms',xaxis="time",yaxis="amp",field='0~5', avgchannel='3072',avgspw=True) # Pause script if you are running in scriptmode user_check=raw_input('Return to continue script\n') #### Data at times > 11:04:00 show lots of scatter: the elevation of #### these data are very low and are worth flagging. One can plotxy #### with xaxis="time",yaxis="elevation" to see this. flagdata(vis="g19_d2usb_targets.ms",field="",timerange='>11:04:00') #### For now to replot you must change selection. #### Note that because the Doppler tracked frequency was 13CO (+ 10GHz #### in the USB), the velocity diplayed in in plotxy will not be #### correct unless you set the restfreq to the line of #### interest. Plotms cannot do this yet. #### CO(2-1) restfreq=230.53797 GHz clearstat() plotxy(vis="g19_d2usb_targets.ms",xaxis="velocity",yaxis="amp", datacolumn="data",iteration="",selectdata=True, antenna='',spw="0~23",field="0",timerange="",averagemode="vector", restfreq='230.53797GHz', width="1",timebin="all",crossscans=T,crossbls=T, stackspw=False,subplot=111,interactive=T,figfile="") # Pause script if you are running in scriptmode user_check=raw_input('Return to continue script\n') #### The other weak line detected in this sideband is not visible in #### a vector average plot. #### A-CH3OH (8 -1 8 0 -- 7 +0 7 0) 229.75876 GHz ################################################################### #### CONTINUUM IMAGING AND SUBSTRACTION ################################################################### ################################################################### # Use field 0 to pick line free spw clearstat() plotxy(vis="g19_d2usb_targets.ms",xaxis="channel",yaxis="amp", datacolumn="data",iteration="",selectdata=True, antenna='',spw="0~23",field="0",timerange="",averagemode="vector", restfreq='230.53797GHz', width="1",timebin="all",crossscans=T,crossbls=T, stackspw=False,subplot=111,interactive=T,figfile="") # Pause script if you are running in scriptmode user_check=raw_input('Return to continue script\n') # The line is confined to spws 13~15. #### Now use new task uvcontsub2 to subtract continuum -- this task #### unlike uvcontsub will subtract the continuum estimate for #### spw that are **not** in the fitspw list. In future releases #### uvcontsub2 will like replace uvcontsub. uvcontsub2(vis='g19_d2usb_targets.ms',fitspw='0~12,16~23', solint='int',combine='spw',want_cont=True) # g19_d2usb_targets.ms.cont/ contains the estimated continuum # g19_d2usb_targets.ms.contsub/ contains the continuum subtracted # line data #### Note that for sparse mosaics like this one, the subparameter #### ftmachine='ft' will typically give better results than the #### default ftmachine='mosaic'. For Nyquist sampled mosaics, #### ftmachine='mosaic' is the best choice, and the only choice for #### heterogeneous arrays. Also note that in a mosaic, the rms noise #### changes as a function of position. Therefore, the rms noise that #### you calculate from the radiometer equation is only valid where #### the combined primary beam response of the mosaic is near #### one. Clean weights the decision of what to clean next by the #### primary beam response function (PBRF), so emission near PBRF=1 #### will be cleaned faster than elsewhere. Thus, if you have emission #### away from the "sweet spot(s)" you will need to clean to a #### threshold lower than you might expect from the anticipated rms #### noise. Development is on-going to provide other options which may #### work better for the case of strong emission near image #### boundaries. clean(vis='g19_d2usb_targets.ms.cont',imagename='g19_d2usb_cont', field='',spw='0~12,16~23', mode='mfs',niter=200,gain=0.1,threshold=0.0, psfmode='clark',imagermode='mosaic',scaletype='SAULT', ftmachine='ft', interactive=T,imsize=500,cell="0.5arcsec", phasecenter="J2000 18h25m56.09 -12d04m28.20", pbcor=F,minpb=0.2) viewer('g19_d2usb_cont.image') # Pause script if you are running in scriptmode user_check=raw_input('Return to continue script\n') # Here you can see one stronger continuum source to the NE (18:25:58.5 # -12:03:58.9) and one weak one to the SW (18:25:52, -12:04:53). The # weak one is only visible after one round of cleaning. Because the # mosaic is sparse you should make a clean mask to get the best # cleaned image. # Make a dirty image of the continuum subtracted data, "continuum # channels" clean(vis='g19_d2usb_targets.ms.contsub',imagename='g19_d2usb_contsub', field='',spw='0~12,16~23', mode='mfs',niter=0,gain=0.1,threshold=0.0, psfmode='clark',imagermode='mosaic',scaletype='SAULT', ftmachine='ft', interactive=F,imsize=500,cell="0.5arcsec", phasecenter="J2000 18h25m56.09 -12d04m28.20", pbcor=F,minpb=0.2) #### You can look at the image with the viewer; you should see #### nothing because the continuum has been subtracted. viewer('g19_d2usb_contsub.image') # Pause script if you are running in scriptmode user_check=raw_input('Return to continue script\n') ################################################################### #### Spectral line imaging ################################################################### #### This image selects a channel range suitable to image the #### 12CO(2-1) line. Notice use of restfreq parameter to get velocity #### scale correct. clean(vis='g19_d2usb_targets.ms.contsub',imagename='g19_d2usb_coline', field='',spw='', mode='channel',start=1730,nchan=40,width=3, niter=1000,gain=0.1,threshold=0.2, psfmode='clark',imagermode='mosaic',scaletype='SAULT', ftmachine='ft',restfreq='230.53797GHz', interactive=T,npercycle=400, imsize=500,cell="0.5arcsec", phasecenter="J2000 18h25m56.09 -12d04m28.20", pbcor=F,minpb=0.15) #### It is essential to run with interactive=T and make a clean #### mask. A clean mask is available in the original data directory #### spw_fits_files. If you would like to use it uncomment lines below. #clean(vis='g19_d2usb_targets.ms.contsub',imagename='g19_d2usb_coline', # field='',spw='', # mode='channel',start=1730,nchan=40,width=3, # niter=1000,gain=0.1,threshold=0.2, # psfmode='clark',imagermode='mosaic',scaletype='SAULT', # ftmachine='ft',restfreq='230.53797GHz', # interactive=F,npercycle=400, # imsize=500,cell="0.5arcsec", # mask='spw_fits_files/coline.mask', # phasecenter="J2000 18h25m56.09 -12d04m28.20", # pbcor=F,minpb=0.15) viewer('g19_d2usb_coline.image') # Pause script if you are running in scriptmode user_check=raw_input('Return to continue script\n') #### Now make 0th and 1st moment maps. immoments(imagename='g19_d2usb_coline.image',moments=[0],axis='spectral', chans='6~32',outfile='g19_d2usb_coline.mom0') viewer('g19_d2usb_coline.mom0') # Pause script if you are running in scriptmode user_check=raw_input('Return to continue script\n') immoments(imagename='g19_d2usb_coline.image',moments=[1],axis=3, chans='6~32',outfile='g19_d2usb_coline.mom1', includepix=[0,2]) immath(outfile='g19_d2usb_coline.mom1masked',mode='evalexpr', expr='"g19_d2usb_coline.mom1"["g19_d2usb_coline.mom0">8]') viewer('g19_d2usb_coline.mom1masked') # Pause script if you are running in scriptmode user_check=raw_input('Return to continue script\n') #### By clicking on "folder" icon in viewer window you can bring up #### the moment0 as a contour image. The strange appearence of the the #### mom1 map toward the NE source is probably due to more than one #### outflow being superposed. #### Note that one could and should correct the image for the #### primary beam response using immath to divide the image by the #### .flux image. One would then need to make a mask for the moment #### images in order to hide the ugly stuff at the edges.