casa学习代码记录

标题

casa是专门用于处理VLA望远镜的集成软件，相应的教程对于干涉阵处理来说可谓是非常详细，甚至对于其他干涉阵的pipline来说，一些参数的意义甚至可以参考回来casa教程的的说明：以下记录一下相关的代码和主要的流程：

ps：很多具体到实际的物理意义我也说不清楚太多东西，以后可能理解深了，会尝试添加一下说明。

https://casaguides.nrao.edu/index.php/Karl_G._Jansky_VLA_Tutorials#Self-calibration_of_VLA_Data
这个网址是教程大全，还是要参考着官方的教程来看

https://data.nrao.edu/portal/#/
这里可以下载到开放的数据，选择你需要的项目的数据
下面开始全部流程

查看观测日志

https://www.vla.nrao.edu/cgi-bin/oplogs.cgi
先确定你数据的日期(可以先看（2）确定日期也可以在你下载数据那里点击scan也有一些信息)，然后在这个网站上查询当天本项目望远镜的情况，例如下面写着data lost的那就记得要先flag掉（也有可能在你下载数据之前望远镜已经flag了）

例如这个就写着天线X在XX时间lost，你在后面的操作中就要flag掉这些数据

确定数据的各种信息

先对数据进行解压tar -xvf 文件名，获得.ms文件

msfile = 'spw6.ms'
listobs(vis=msfile)

一共有这些源，前面几个是校准源（建议记住这些源名字，虽然你可以使用filed0，1，2这些编号去索引，但是容易搞混）

G16.7+0.1(targe), J1331+3030, J1407+2827,J1822-0938(phase calibration),0137+331=3C48(flux calibration)

一共15个频率窗口，前面和后面是重复的

选择这个窗口spw6

阵列天线分布图，找到一个最中间作为参考天线,例如e25


plotants(vis=msfile,figfile='plotants_antenna_layout.png')

ea25

本教程只处理一个窗口的数据（当你数据量不大的时候，也可以全部窗口一些处理）

简单看看spw6的10~55通道数据的情况（因为每个spw前后的一些通道总是不太好）

plotms(vis =msfile, xaxis='time',yaxis='amp', correlation='RR,LL',coloraxis='field', avgchannel='64')

plotms(vis =msfile, xaxis='channel',yaxis='amp', correlation='RR,LL',coloraxis='field', avgtime='1e8')

感觉这个窗口的数据还可以，先把他分出来成新的文件吧，减少文件的大小

split(vis=msfile,outputvis='XXX_spw6.ms',spw='6:10~55',datacolumn='all')
#这个spw的格式是6窗口的10到55通道
msfile = 'XXX_spw6.ms' #后面就用这个新的.ms文件了

Initi Flag

对数据进行一个基础的flag

flagdata(vis=msfile, mode='quack', quackinterval=5.0, quackmode='beg') 
flagdata(vis=msfile, mode='manual', antenna='ea13')
#flagdata(vis=msfile, mode='manual', scan='1')

#分别是对每段时间的起始几秒的数据flag；对观测日志里面不好的数据进行flag；
#对第一个scan数据(从listobs中发现观测时间只有几秒）flag

#对校准源的数据做一个更加细致的flag，flag掉那些一眼看上去的坏的数据点，以及标记看起来不好的数据点。通过plotms打开校准源。

```python
plotms(vis=msfile,field='',
       xaxis='chan',yaxis='amp',correlation='RR,LL',
       avgtime='1e8',spw='', coloraxis='baseline', iteraxis='spw')

plotms(vis=msfile,
       xaxis='time',yaxis='amp',correlation='RR,LL',
       avgchannel='64', coloraxis='field', iteraxis='spw')

利用下面的工具栏确定异常数据是来源是哪里（eg：天线？scan？时间？）

ps：flag之后的数据以修改的数据的格式存放在一个flag文件里面，每次做了flag或者校准命令都会保存，但是不显示具体实现了什么操作，你可以通过以下两个例子来给你当前flag状态的保存或者回溯到上一个flag，当然casa也会自动保存flag的版本，但是只有编号没有flag的信息。

flagmanager(vis = msfile,
            mode = 'save',
            versionname = 'initialflags',
            comment = 'initial flags, including quack',
            merge='replace')
            
flagmanager(vis = msfile,
            mode = 'restore',
            versionname = 'flags.flagdata_2',
            merge='replace')

# versionname就是flag的版本，第一个是保存，第二个是回溯。

Initial Flux Density Scaling

gencal(vis=msfile, caltable='spw6.antpos',caltype='antpos') 
# 天线的位置校准，有可能没有

setjy(vis=msfile, listmodels=True)
setjy(vis=msfile, field='0137+331=3C48',standard='Perley-Butler 2017',model='3C48_C.im',usescratch=True,scalebychan=True,spw='')
#先看有那些模型，然后确定流量校准源的模型

refantenna='ea25'

Initial Phase Calibration

gaincal(vis='spw6.ms', caltable='spw6.G0all',field='0137+331=3C48, J1822-0938', refant=refantenna, spw='0:27~36',gaintype='G',calmode='p', solint='int',minsnr=5,gaintable=['spw6.antpos'])
#实际上并没有用到这个输出的增益表，只是简单看看增益的情况

plotms(vis='spw6.G0all', xaxis='time',yaxis='phase',
coloraxis='corr',iteraxis='antenna',plotrange=[-1,-1,-180,180])
#每个gain后面都值得plotms一下表格确定一下输出情况

gaincal(vis='spw6.ms', caltable='spw6.G0',field='0137+331=3C48', refant=refantenna, spw='0:27~36', calmode='p', solint='int',minsnr=5, timerange='', gaintable=['spw6.antpos'])
#相位校准，3C138是流量校准源，这个时间的选择不一定需要，如果这个源有两端的观测，你可以像我一样用时间去限制只要前面一段的观测。

plotms(vis='spw6.G0',xaxis='time',yaxis='phase',coloraxis='corr',iteraxis='antenna',plotrange=[-1,-1,-180,180])

Delay calibration

gaincal(vis=msfile, caltable='spw6.K0',field='0137+331=3C48',refant=refantenna,spw='0:10~55',gaintype='K',solint='inf',combine='scan',minsnr=5,gaintable=['spw6.antpos','spw6.G0'])

plotms(vis='spw6.K0',xaxis='antenna1',yaxis='delay',coloraxis='baseline')

Bandpass calibration

bandpass(vis=msfile,caltable='spw6.B0',field='0137+331=3C48',spw='0:10~55',refant=refantenna,combine='scan',solint='inf',bandtype='B',gaintable=['spw6.antpos','spw6.G0','spw6.K0'])

plotms(vis='spw6.B0',xaxis='chan',yaxis='amp',coloraxis='corr',gridrows=2,gridcols=2, spw='0:10~55',iteraxis='antenna')

plotms(vis='spw6.B0',xaxis='chan',yaxis='phase',coloraxis='corr',gridrows=2,gridcols=2, spw='0:10~55',iteraxis='antenna',plotrange=[-1,-1,-180,180])

flagdata(vis=msfile, mode='manual', antenna='ea04&ea25')

flag ea04&ea25？

Gain Calibration

gaincal(vis=msfile,caltable='spw6.G1',field='0137+331=3C48',spw='0:10~55',
solint='inf',refant=refantenna,gaintype='G',calmode='ap',solnorm=False,
gaintable=['spw6.antpos','spw6.K0','spw6.B0'],interp=['','','nearest'])

gaincal(vis=msfile,caltable='spw6.G1',field='J1822-0938',spw='0:10~55',
solint='inf',refant=refantenna,gaintype='G',calmode='ap',gaintable=['spw6.antpos','spw6.K0','spw6.B0'],interp=['','','nearest'],solnorm=False,append=True)

#照常看一下输出的情况
plotms(vis='spw6.G1',xaxis='time',yaxis='amp', correlation='/',coloraxis='baseline')
plotms(vis='spw6.G1',xaxis='time',yaxis='phase',correlation='/',coloraxis='baseline', plotrange=[-1,-1,-180,180])

Scaling the Amplitude Gains

myscale = fluxscale(vis='spw6.ms',caltable='spw6.G1',fluxtable='spw6.fluxscale1',reference='0137+331=3C48',transfer=['J1822-0938'],incremental=False)

plotms(vis='spw6.fluxscale1',xaxis='time',yaxis='amp',
correlation='R',coloraxis='baseline')

https://science.nrao.edu/facilities/vla/docs/manuals/observing/callist

VLA公布的校准源流量密度等信息

Applying the Calibration

applycal(vis=msfile,field='0137+331=3C48',gaintable=['spw6.antpos','spw6.fluxscale1','spw6.K0','spw6.B0'],gainfield=['','0137+331=3C48','',''],interp=['','nearest','',''],calwt=False)

applycal(vis=msfile,field='J1822-0938',gaintable=['spw6.antpos','spw6.fluxscale1','spw6.K0','spw6.B0'],gainfield=['','J1822-0938','',''],interp=['','nearest','',''],calwt=False)

#对target
applycal(vis=msfile,field='G16.7+0.1',gaintable=['spw6.antpos','spw6.fluxscale1','spw6.K0','spw6.B0'],gainfield=['','J1822-0938','',''],interp=['','nearest','',''],calwt=False)

checking

plotms(vis='spw6.ms',field='*3C48',correlation='RR,LL',avgtime='60',xaxis='channel',yaxis='amp',ydatacolumn='corrected',coloraxis='corr', spw='0:10~55')

plotms(vis='spw6.ms',field='*3C48',correlation='RR,LL',xaxis='time',yaxis='amp',ydatacolumn='corrected',coloraxis='corr', spw='0:10~55',avgchannel='64')

#first flag !ea17&ea25

#second flag scan134 ? ea25&ea28,ea25&ea26,ea24&ea25,ea25? 
#or change refantenna ea25 from target

baseline,antenna,scan,channel,time

flagmanager(vis = msfile,mode = 'save',versionname = 'initialflag_1',comment='first flag',merge='replace')

flagdata(vis=msfile, mode='manual', antenna='ea17&ea25')

mapping


split(vis=msfile,outputvis='spw6_cal.ms',datacolumn='corrected',field='0137+331=3C48,J1822-0938', spw='0:10~55')

split(vis='spw6.ms',outputvis='spw6_target.ms',datacolumn='corrected',field='G16.7+0.1', spw='0:10~55')

statwt(vis='spw6_target.ms',datacolumn='data')

cleansource='M31_spw6_cal.ms'

#计算cell可以使用
plotms(vis='spw6_target.ms',xaxis='uvwave',yaxis='amp',ydatacolumn='data', field='G16.7+0.1',avgtime='30',correlation='RR')

plotms(vis='spw6_target.ms',xaxis='real',yaxis='imag',ydatacolumn='data', field='G16.7+0.1',avgtime='30',correlation='RR')

plotms(vis='spw6_target.ms',xaxis='uwave',yaxis='vwave',ydatacolumn='data', field='G16.7+0.1',avgtime='30',correlation='RR')

np.degrees(1/17500)*3600 ≈ 11.7
11.7/5 ≈ 3

tclean(vis='spw6_cal.ms',field='0137+331=3C48',imagename='cal_0137+331=3C48',specmode='mfs', datacolumn='data',imsize=100,cell=['3.0arcsec'],weighting='briggs',robust=0,niter=1000,interactive=True)

tclean(vis='spw6_cal.ms',field='J1822-0938',imagename='cal_J1822-0938',specmode='mfs', datacolumn='data',imsize=100,cell=['3.0arcsec'],weighting='briggs',robust=0,niter=1000,interactive=True)

tclean(vis='spw6_target.ms', imagename='SNR',specmode='mfs',niter=20000,gain=0.1, threshold='0.0mJy',deconvolver='multiscale',scales=[0, 5, 15, 45], smallscalebias=0.9,interactive=True,imsize=[200,200], cell=['3.0arcsec','3.0arcsec'],stokes='I',weighting='briggs',robust=0.5,pbcor=False,savemodel='modelcolumn')