ALMA simulation脚本

#define working direction path:
path = '/home/mengyao/test/ALMA_simulation/W51n/uv-range_sim/'
#define project name: 
project_name='gauss_point'
#define parameters for tclean:
imsize=[1024,1024]
cell='0.005arcsec'
threshold='0.15mJy'
uvrange='> 910klambda'
#mask='circle[[512pix,512pix],400pix]'
mask = '/home/mengyao/test/ALMA_simulation/W51n/uv-range_sim/mask.region'


#生成一个component列表进行观测模拟
def make_component():
    os.system('rm -rf '+ project_name +'_components.cl')
    cl.done()
    cl.addcomponent(dir='J2000 19h23m40.05s +14d31m05.50s',flux=5000, fluxunit='mJy',freq='226.379GHz',shape='Gaussian',majoraxis='0.6arcsec',minoraxis='0.3arcsec',positionangle='-47.0deg')
    cl.addcomponent(dir='J2000 19h23m40.041s +14d31m05.647s',flux=5000, fluxunit='mJy',freq='226.379GHz',shape='Gaussian',majoraxis='0.6arcsec',minoraxis='0.3arcsec',positionangle='-47.0deg')
    cl.addcomponent(dir='J2000 19h23m40.056s +14d31m05.379s',flux=5000, fluxunit='mJy',freq='226.379GHz',shape='Gaussian',majoraxis='0.6arcsec',minoraxis='0.3arcsec',positionangle='-45.0deg')
    cl.addcomponent(dir='J2000 19h23m40.026s +14d31m05.744s',flux=5000, fluxunit='mJy',freq='226.379GHz',shape='Gaussian',majoraxis='0.6arcsec',minoraxis='0.3arcsec',positionangle='-47.0deg')
    cl.addcomponent(dir='J2000 19h23m40.067s +14d31m05.229s',flux=5000, fluxunit='mJy',freq='226.379GHz',shape='Gaussian',majoraxis='0.6arcsec',minoraxis='0.3arcsec',positionangle='-47.0deg')
    cl.addcomponent(dir='J2000 19h23m40.05s +14d31m05.50s',flux=1, fluxunit='mJy',freq='226.379GHz',shape='point')
    cl.addcomponent(dir='J2000 19h23m40.041s +14d31m05.647s',flux=1, fluxunit='mJy',freq='226.379GHz',shape='point')
    cl.addcomponent(dir='J2000 19h23m40.056s +14d31m05.379s',flux=1, fluxunit='mJy',freq='226.379GHz',shape='point')
    cl.rename(project_name +'_components.cl')
    cl.done()

#开始模拟，模拟可以使用component列表也可以使用模型fits文件，
def simulation():
    default('simobserve') 
    simobserve(project=project_name + 'project',
    complist=project_name +'_components.cl',
    compwidth='1GHz', #设置带宽
    direction='J2000 19h23m40.05s +14d31m05.50s', #设置观测中心点
    obsmode='int',  #观测模式，int表示干涉仪观测
    setpointings=False, #False表示单点观测，True表示马赛克观测
    ptgfile = 'w51n.ptg.txt',  #自己制作的单点观测文件，包含了观测中心点的坐标，如果是马赛克观测模式此处可以不设定。
    integration='10s',          #单次积分时间，
    refdate='2015/10/31',   #观测日期
    antennalist='test.cfg',   #望远镜配置文件，可以在/casa-release-5.4.0-70.el6/data/alma/simmos路径下找到，如果没有你需要的望远镜配置文件，可以参考其他的自己制作一个，望远镜的数据可以用listobs来查看。
    totaltime='8756s',    #总积分时间
    mapsize='5.0arcsec',   #观测尺寸
    thermalnoise='',
    overwrite=True)


def do_tclean():
    filename=path + project_name + 'project/'+ project_name + 'project.test.ms'
    os.system('rm -rf '+ project_name+'.*')
    tclean(field='0',
       spw='0',
       vis=filename,
       imagename=project_name,
       imsize=imsize,
       cell=cell,
       threshold=threshold,
       outframe='LSRK',
       gridder='standard',
       niter=1000000,
       datacolumn='data',
       savemodel='none',
       deconvolver='hogbom',
       weighting='briggs',
       robust=0.5,
       mask=mask,
       restoringbeam=['0.0286arcsec','0.0218arcsec','-46.688deg'],
       specmode='mfs',
       interactive=False)


def do_tclean_910k():
    filename=path + project_name + 'project/'+ project_name + 'project.test.ms'
    os.system('rm -rf '+ project_name+'_910k.*')
    tclean(field='0',
       spw='0',
       vis=filename,
       imagename=project_name+'_910k',
       uvrange=uvrange,
       imsize=imsize,
       cell=cell,
       threshold=threshold,
       outframe='LSRK',
       gridder='standard',
       niter=1000000,
       datacolumn='data',
       savemodel='none',
       deconvolver='hogbom',
       weighting='briggs',
       robust=0.5,
       mask=mask,
       restoringbeam=['0.0286arcsec','0.0218arcsec','-46.688deg'],
       specmode='mfs',
       interactive=False)


def export_fits():
    os.system('rm -rf sim2_*.fits')
    exportfits(imagename=project_name+'.image',fitsimage='sim2_full.fits')
    exportfits(imagename=project_name+'_910k.image',fitsimage='sim2_910k.fits')
    os.system('rm -rf '+project_name+'.* ' +project_name+'_910k.*')


make_component()
simulation()
do_tclean()
do_tclean_910k()
export_fits()