casa天文软件全代码记录

最近的工作使用casa的频率不少，但是我经常写了代码或者弄好的脚本想要复用找起来很麻烦。我想记录一下成图所有使用过的东西，到时候我复用起来比较方便。

【自动clean每个field】

msfile= '20A-346.60247-spw22.ms'
msfile_h =  msfile[:-3]+'-hanning.ms'
refantenna='ea15'
amp_cals = '3C48'
phase_cals = 'J2355+4950'
all_cals = '3C138,J0319+4130,3C48,J2355+4950'
cals_1 = '3C138'
cals_2 = 'J0319+4130'
targe_n = '*LARGE*'cal_out_name = 'M31_spw22'
fig_dir = 'fig/'
cleansource=cal_out_name+'_hanning_cal.ms'
targetsource = cal_out_name+'_hanning_target.ms'for i in range(0,49):outname = 'single-field-tclean-10cell-1500niter/M31_spw22_field-%s' %(i)field_name = '*LARGE_%s*' %(i)tclean(vis=targetsource, imagename=outname, specmode='mfs',niter=1500, gain=0.1, threshold='0mJy',gridder='standard',deconvolver='multiscale',scales=[0, 5, 15, 45], smallscalebias=0.9,interactive=False, imsize=400, cell=['10.0arcsec','10.0arcsec'], stokes='I',weighting='briggs',robust=2, pbcor=False,savemodel='modelcolumn', field=field_name, pblimit=0)

【合并多个field成为一张图】

例如我有好多个filed的图像，想要合并在一起。

#!/usr/bin/env pythonimport globfrom casatools import linearmosaic as LM
lm = LM()"""
ss_image = 'M31_spw6_field-XX.image'
ss_pb = 'M31_spw6_field-XX.pb'
ptList = [25, 37, 38, 39, 40, 41]
images = []
weightimages = []
for pts in ptList:images.append(ss_image.replace('XX',str(pts)))weightimages.append(ss_pb.replace('XX',str(pts)))
print(images, weightimages)
"""images = glob.glob('./*.image')
weightimages =  []for term in images:weightimages.append(term.replace('image','pb'))
#需要两个文件单个filed的image文件和对应的pb文件（主波束模型）
#lm.defineoutputimage(nx=4000,ny=4000,cellx='2arcsec', imagecenter='00h41m48.00 +40d28m21.00',outputimage='test.linmos')lm.defineoutputimage(nx=7200,ny=7200,cellx='2arcsec', imagecenter='00h42m44.00 +41d16m08.00',outputimage='m31.all.linmos')
#这个像素点数量和单个间隔，图像中心，输出文件名词都需要更改一下lm.makemosaic(images=images,weightimages=weightimages)

【自校准】

1计算N值

#生成脏图(这里我只对38field做了自校准)
tclean(vis='M31_spw22_hanning_target.ms',imagename='obj.dirty', datacolumn='data', imsize=400, cell='10arcsec', pblimit=0, gridder='standard', field='*38*')#起始模型，robust取0比较好
tclean(vis='M31_spw22_hanning_target.ms', imagename='obj.prelim_clean', datacolumn='data', imsize=400, cell='10arcsec', pblimit=0, gridder='standard', deconvolver='mtmfs', nterms=2, niter=1000, interactive=True, weighting='briggs', robust=0, savemodel='modelcolumn', field='*38*')#打开carta看看图像的峰值和rms的比值，这个叫做动态范围。#验证模型是否正确保存进去
plotms(vis='M31_spw22_hanning_target.ms', xaxis='UVwave', yaxis='amp', ydatacolumn='model', avgchannel='64',field='*38*')#初始校准表
gaincal(vis='M31_spw22_hanning_target.ms', caltable='selfcal_initial.tb', solint='int', refant='ea15', calmode='p', gaintype='G', minsnr=0, field='*38*')gaincal(vis='M31_spw22_hanning_target.ms',caltable='selfcal_combine_pol_solint_3.tb',solint='int',refant='ea15',calmode='p',gaintype='T', minsnr=0, field='*38*')
gaincal(vis='M31_spw22_hanning_target.ms',caltable='selfcal_combine_pol_solint_6.tb',solint='6s',refant='ea15',calmode='p',gaintype='T', minsnr=0, field='*38*')
gaincal(vis='M31_spw22_hanning_target.ms',caltable='selfcal_combine_pol_solint_12.tb',solint='12s',refant='ea15',calmode='p',gaintype='T', minsnr=0, field='*38*')
gaincal(vis='M31_spw22_hanning_target.ms',caltable='selfcal_combine_pol_solint_24.tb',solint='24s',refant='ea15',calmode='p',gaintype='T', minsnr=0, field='*38*')
gaincal(vis='M31_spw22_hanning_target.ms',caltable='selfcal_combine_pol_solint_48.tb',solint='48s',refant='ea15',calmode='p',gaintype='T', minsnr=0, field='*38*')
gaincal(vis='M31_spw22_hanning_target.ms',caltable='selfcal_combine_pol_solint_96.tb',solint='96s',refant='ea15',calmode='p',gaintype='T', minsnr=0, field='*38*')#求解解的间隔选取
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
tb.open( 'selfcal_combine_pol_solint_6.tb' )
snr_6s = tb.getcol( 'SNR' ).ravel()
tb.close()
tb.open( 'selfcal_combine_pol_solint_12.tb' )
snr_12s = tb.getcol( 'SNR' ).ravel()
tb.close()
tb.open( 'selfcal_combine_pol_solint_24.tb' )
snr_24s = tb.getcol( 'SNR' ).ravel()
tb.close()
tb.open( 'selfcal_combine_pol_solint_48.tb' )
snr_48s = tb.getcol( 'SNR' ).ravel()
tb.close()
plt.hist( snr_6s, bins=50, density=True, histtype='step', label='6 seconds' )
plt.hist( snr_12s, bins=50, density=True, histtype='step', label='12 seconds' )
plt.hist( snr_24s, bins=50, density=True, histtype='step', label='24 seconds' )
plt.hist( snr_48s, bins=50, density=True, histtype='step', label='48 seconds' )
plt.legend( loc='upper right' )
plt.xlabel( 'SNR' )
print( 'P(<=6) = {0}  ({1})'.format( stats.percentileofscore( snr_6s, 6 ), '6s' ) )
print( 'P(<=6) = {0}  ({1})'.format( stats.percentileofscore( snr_12s, 6 ), '12s' ) )
print( 'P(<=6) = {0}  ({1})'.format( stats.percentileofscore( snr_24s, 6 ), '24s' ) )
print( 'P(<=6) = {0}  ({1})'.format( stats.percentileofscore( snr_48s, 6 ), '48s' ) )#应用这个求解间隔
gaincal(vis='M31_spw22_hanning_target.ms', caltable='selfcal_combine_pol_solint_12_minsnr_6.tb', solint='12s', refant='ea15', calmode='p', gaintype='T', minsnr=6,field='*38*')#应用
applycal(vis='M31_spw22_hanning_target.ms', gaintable='selfcal_combine_pol_solint_12_minsnr_6.tb',field='*38*')#第一次自校准的检查
tclean(vis='M31_spw22_hanning_target.ms', imagename='obj.selfcal1_clean.3arcmin', datacolumn='corrected', imsize=400, cell='10arcsec', pblimit=0, gridder='standard', deconvolver='mtmfs', nterms=2, niter=1000, interactive=True, weighting='briggs', robust=0, savemodel='none',field='*38*')

【linearmosaic】

把下面的写出.py脚本，然后在casa里面使用execfile(‘script.py’)运行脚本就可以了。注意这个输入不能fits文件

#!/usr/bin/env pythonimport globfrom casatools import linearmosaic as LM
lm = LM()images = glob.glob('./m31.spw22.f*.wsclean-MFS-image')
weightimages =  glob.glob('./M31_spw22_field-*.pb')print(images)
print(weightimages)
#lm.defineoutputimage(nx=4000,ny=4000,cellx='2arcsec', imagecenter='00h41m48.00 +40d28m21.00',outputimage='test.linmos')lm.defineoutputimage(nx=3000,ny=3000,cellx='5arcsec', imagecenter='00h42m44.00 +41d16m08.00',outputimage='m31.all.20A.5cell.wsclean.linmos.image')lm.makemosaic(images=images,weightimages=weightimages)