casa6

###############################################

#                                             #

# Regression/Benchmarking Script for NGC 4826 #

#        Single Dish + Synthesis              #

###############################################

​

import time

import os

​

os.system('rm -rf n4826_t* mosaic.flux mosaic.model mosaic.image mosaic.residual n12m_gaussian.im')

​

pathname=os.environ.get('CASAPATH').split()[0]

datapath = pathname + '/data/regression/ATST3/NGC4826/'

​

print '--Copy data to local directory--'

mspath='cp -r '+datapath+'n4826_both.ms .'

os.system(mspath)

os.system('chmod -R a+wx n4826_both.ms')

​

startTime = time.time()

startProc = time.clock()

​

​

print '--Feather--'

# Starting from:

#    BIMA mosaic image (Moment 0) : n4826_mom0.im

#    NRAO 12m OTF image (Moment 0): n4826_12mmom0.im

default('feather')

feather('n4826_tfeather.im',datapath+'n4826_mom0.im',datapath+'n4826_12mmom0.im')

feathertime = time.time() 

#combo: Max:1.533498e+02        Flux:1.523515e+03 Jy    rms:1.187669e+01

#Pcombo:Max:1.511604e+02        Flux:2.016790e+03 Jy    rms:1.873132e+01

#BIMA:  Max:1.627327e+02        Flux:1.035314e+04 Jy    rms:1.587101e+01

#12m:   Max:9.939910e+02        Flux:2.600573e+03 Jy    rms:1.159703e+02

​

print '--Feather cube--'

# Starting from:

#    BIMA mosaic image (Moment 0) : n4826_bima.im

#    NRAO 12m OTF image (Moment 0): NGC4826.12motf.chan.fits

default('importfits')

importfits(datapath+'NGC4826.12motf.chan.fits','n4826_t12mchan.im')

default('imhead')

imhead('n4826_t12mchan.im',mode='put',hdkey='bunit',hdvalue='Jy/beam')

#imhead('n4826_t12mchan.im',mode='put',hdkey='beam',hdvalue='55arcsec, 55arcsec, 0deg')

###wow what an idea...more to type to change beam

imhead('n4826_t12mchan.im',mode='put',hdkey='beammajor',hdvalue='55arcsec')

imhead('n4826_t12mchan.im',mode='put',hdkey='beamminor',hdvalue='55arcsec')

imhead('n4826_t12mchan.im',mode='put',hdkey='beampa',hdvalue='0deg')

​

​

default('feather')

feather('n4826_tfeathercube.im',datapath+'n4826_bima.im','n4826_t12mchan.im')

feathercubetime = time.time()

# Maximum value   1.881626e+00 at [113, 115, 1, 25] (12:56:45.388, +21.40.51.100, I, 1.15172e+11Hz)

# Flux density  =   1.054628e+02 Jy

​

print '--Feather cube - create synth image--'

# Starting from:

#    BIMA calibrated visibilities: n4826_both.ms

#    NRAO 12m OTF cube: NGC4826.12motf.chan.fits

default('clean')

clean(vis='n4826_both.ms', imagename='mosaic',

      nchan=30, start=45, width=4, spw='0~2',

      field='0~6',

      cell=[1., 1.], imsize=[256, 256],

      stokes='I',

      mode='channel',

      interpolation='nearest',

      psfmode='clark',

      imagermode='mosaic',

      ftmachine='mosaic',

      niter=500,

      gain=0.1,                     # Specified for certainty, 6/10/2009

      scaletype='SAULT', minpb=0.1) # Changed from 0.01 6/10/2009

default('importfits')

importfits(datapath+'NGC4826.12motf.chan.fits', 'n4826_t12mchan2.im')

default('imhead')

imhead('n4826_t12mchan2.im', mode='put', hdkey='bunit', hdvalue='Jy/beam')

#imhead('n4826_t12mchan2.im',mode='put',hdkey='beam',hdvalue='55arcsec,55arcsec,0deg')

imhead('n4826_t12mchan2.im',mode='put',hdkey='beammajor',hdvalue='55arcsec')

imhead('n4826_t12mchan2.im',mode='put',hdkey='beamminor',hdvalue='55arcsec')

imhead('n4826_t12mchan2.im',mode='put',hdkey='beampa',hdvalue='0deg')

default('feather')

feather('n4826_tfeathersynth.im','mosaic.image','n4826_t12mchan2.im')

#Combo:   Max: 1.881626e+00   Flux: 1.054628e+02 Jy   Rms: 8.911050e-02

#Pcombo:  Max: 1.681789e+00   Flux: 1.454376e+02 Jy   Rms: 8.196454e-02

feathersynthtime = time.time()

​

print '--Single Dish as Model--'

# Starting from:

#    BIMA calibrated visibilities: n4826_both.ms