# From dshepher@aoc.nrao.edu Wed Jan 17 15:14:25 2001
# Date: Tue, 31 Oct 2000 11:57:15 -0700 (MST)
# Subject: final aips++ script for Orion mosaic
#========================================================================
# Orion X-band Mosaic script for AIPS++
# 
# Created by Debra Shepherd based on a script developed by Tim Cornwell
# 
# This script fills, edits, calibrates, and images the 
# Orion X-band continuum mosaic using AIPS++.  
# 
#========================================================================
# initialize the aips++ path with a command like:
source /aips++/version/aipsinit.csh
# start aips++:
aips++

The following aips++ commands should be run from the 
glish command line interface (CLI) in the AIPS++ xterm directly. 
Alternatively, you can run each task from the GUI interface.  

dowait := T;
#
# Fill data from a remote tape drive:
#
#you must have the line:
   source /aips++/version/aipsinit.csh
#   (or what ever path you have setup for AIPS++)
#   in the VERY TOP of your .cshrc (or .tcshrc) file.  
#   "version" refers to the version of AIPS++ you are using. 
#   This is needed to set up the remote machine's tape drive. 
#   Ignore this if you are filling data from your tape drive
#
include 'vlafiller.g';
myvlafiller:=vlafiller(host="remote_machine" , forcenewserver=F);
ok:=myvlafiller.tapeinput(device="/dev/tape_drive_name" , files=1);
ok:=myvlafiller.output(msname="/fullpathname/orion/orion.ms" ,
overwrite=F);
ok:=myvlafiller.selectproject(project="DSTST" );
ok:=myvlafiller.fill(verbose=T, async=F);
#
# (full path must be specified for remote machine). 
  # Data should be loaded into the measurement set orion.ms 
  # RESULTS should be:
  # The measurement set contains 1094418 rows.
  # The antenna sub-table contains 26 entries
  # The field sub-table contains 12 entries
  # The polarization sub-table contains 1 entries
  # The spectral window sub-table contains 2 entries
myvlafiller.done();

# In an xterm window: copy ms to a save file, just in case:
cp -r orion.ms orion.save.ms

#
# Do the obvious flagging: quack, autocorrelations, known bad antennas 
# quack 1st ints (5 sec ints, 3min/field so set scaninterval=5s)
#

include 'flagger.g';
myflagger:=flagger(msfile="orion.ms" );

# The logger window will show a summary of the data: 
#
==============================================================================
#            MeasurementSet Name:  orion.ms      MS Version 2
#
==============================================================================
#    Observer: unavailable      Project: DSTST   
# Observation: VLA(26 antennas)
# Fields: 12  Name: 0518+165  RA: 05:21:09.89      Dec: +16.38.22.04
(J2000   )
#             Name: 0539-057  RA: 05:41:38.09      Dec: -05.41.49.43
(J2000   )
#             Name: ORION1    RA: 05:35:07.42      Dec: -05.25.36.07
(J2000   )
#             Name: ORION2    RA: 05:35:17.42      Dec: -05.25.36.79
(J2000   )
#             Name: ORION3    RA: 05:35:27.42      Dec: -05.25.37.52
(J2000   )
#             Name: ORION4    RA: 05:35:27.47      Dec: -05.23.07.52
(J2000   )
#             Name: ORION5    RA: 05:35:17.47      Dec: -05.23.06.79
(J2000   )
#             Name: ORION6    RA: 05:35:07.47      Dec: -05.23.06.07
(J2000   )
#             Name: ORION7    RA: 05:35:07.52      Dec: -05.20.36.07
(J2000   )
#             Name: ORION8    RA: 05:35:17.52      Dec: -05.20.36.80
(J2000   )
#             Name: ORION9    RA: 05:35:27.52      Dec: -05.20.37.52
(J2000   )
#             Name: ORION10   RA: 05:35:32.61      Dec: -05.16.07.88
(J2000   )
# Data records: 1094418       Total integration time = 8550 seconds
#    Observed from   21-Sep-2000/11:15:48   to   21-Sep-2000/13:38:18
# Spectral Windows: 2
#    Ref.Freq    #Chans  Resolution  TotalBW     
#    8435.1 MHz  1       60500  kHz  50000  kHz  
#    8485.1 MHz  1       60500  kHz  50000  kHz  
# Polarization setups: 1
#    Correlations    
#    RR  LL  RL  LR  
#
------------------------------------------------------------------------------
# 
# Tables(rows):   (-1 = table absent)
#
MAIN(1094418)   ANTENNA(26)   DATA_DESCRIPTION(2)   DOPPLER(-1)   FEED(26)   FIELD(12)
#
FLAG_CMD(0)   FREQ_OFFSET(-1)   HISTORY(0)   OBSERVATION(1)   POINTING(0)   POLARIZATION(1)   PROCESSOR(0)   SOURCE(-1)   SPECTRAL_WINDOW(2)   STATE(0)   SYSCAL(-1)   WEATHER(-1)
#
==============================================================================

# Quack first scan interval:
ok:=myflagger.quack(scaninterval="5.5s" , delta=[value=5.0, unit="s" ],
trial=F);
  # takes a while...
  # RESULTS: 
  # Flagging 32292 rows -- out of 1094418 rows. about 0.3% -- OK
# flag autocorrelations for easier plotting (can't use them here). 
ok:=myflagger.auto(trial=F);
  # RESULTS: 
  # Flagging 81068 rows -- OK
# Flag outrageously large points. 
# No data should be over 1000 Jy, for sure, or equal to zero:
ok:=myflagger.filter(column="DATA" , operation="range" ,
comparison="Amplitude" ,
	range='1e-6Jy 1e3Jy', trial=F);
  # Filtering (in range mode) has flagged 33630 points. about 0.3% -- OK
myflagger.done();

# Look at the data with visplot:
include 'visplot.g';
myvisplot:=visplot(msfile="orion.ms" , nrows=16384, edit=T, 
  flagfile=[i_am_unset="i_am_unset" ], displaywidth=600,
displayheight=350);
  # Interactively flag bad data
  # Be sure to flag the 400-1000 Jy points in RR scattered in the data
set!
myvisplot.done();

# ANT 21 looks bad, flag data with flagger:
include 'flagger.g';
myflagger:=flagger(msfile="orion.ms" );
# select what antenna to flag:
ok:=myflagger.setantennas(ants=21);
# flag the bad antenna at all time ranges:
ok:=myflagger.timerange(starttime="21-SEP-2000/11:15:48" ,
                        endtime="21-SEP-2000/13:38:18" , trial=F);
myflagger.done();

#
# Set the fluxes for the flux and phase calibraters
#
include 'imagerwizard.g';
myimager:=imager(filename="orion.ms" );
ok:=myimager.setjy(fieldid=1, spwid=-1, fluxdensity=-1.0);
  # 3c138 (expect 2.52 Jy at X band):
  # RESULTS: 
  # 0518+165  spwid=  1  [I=2.517, Q=0, U=0, V=0] Jy, (Perley-Taylor 95)
  # 0518+165  spwid=  2  [I=2.505, Q=0, U=0, V=0] Jy, (Perley-Taylor 95)
ok:=myimager.setjy(fieldid=2, spwid=-1, fluxdensity=-1.0);
  # phase cal (expect 1 Jy at X band): 
  # RESULTS: 
  # 0539-057  spwid=  1  [I=1, Q=0, U=0, V=0] Jy, (default)
  # 0539-057  spwid=  2  [I=1, Q=0, U=0, V=0] Jy, (default)
myimager.done();

#
# Solve for the gains and apply the calibration:
#
include 'calibrater.g';
mycalibrater:=calibrater(filename="orion.ms" );
ok:=mycalibrater.setdata(mode="none" , nchan=1, start=1, step=1, 
	mstart=[value=0.0, unit="km/s" ], mstep=[value=0.0, unit="km/s" ], 
	uvrange=0, msselect='FIELD_ID in [1,2]');
ok:=mycalibrater.setsolve(type="G" , t=60, preavg=0.0, phaseonly=F, 
	table="orion.gcal" , append=F);
ok :=mycalibrater.solve();
# bootstrap flux density scale from 3c138 to the gain cal:
ok:=mycalibrater.fluxscale(tablein="orion.gcal" ,
tableout="orion.ref.gcal" ,
	reference='0518+165 ', transfer='0539-057 ');
  # RESULTS: 
  # Set the flux density scale using reference calibrators
  # Flux density for 0539-057 (spw=1) is:     1.077 +/-  0.003 Jy
  # Flux density for 0539-057 (spw=2) is:     1.075 +/-  0.004 Jy
  # looks OK.
# select all source uv data to be corrected:
ok:=mycalibrater.setdata(mode="none" , nchan=1, start=1, step=1, 
	mstart=[value=0.0, unit="km/s" ], mstep=[value=0.0, unit="km/s" ], 
	uvrange=0, msselect='');
# select the flux-scaled gain cal solutions to be used in interpolation 
# to the source uv data:
ok:=mycalibrater.setapply(type="G" , t=0.0, table="orion.ref.gcal" , 
	select="FIELD_NAME=='0539-057'" );
ok:=mycalibrater.correct();
mycalibrater.done();

# Look at the data with visplot to isolate more bad data based on
residuals:
include 'visplot.g';
myvisplot:=visplot(msfile="orion.ms" , nrows=16384, edit=T, 
  flagfile=[i_am_unset="i_am_unset" ], displaywidth=600,
displayheight=350);
  # plot x = uvdist and time, y = residual amplitudes 
  # delete points with obvious bad residuals. 
myvisplot.done(); applied flags. 
  # This gets rid of most of the spurious stuff, 
  # You may want to flag antenna 18 also since the residuals are high. 
  # do this with flagger.setantennas and flagger.timerange as above. 
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # 
#
# Make a MEM mosaic image
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # 
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # 
# # # # 
# # # # Some details about the mosaic: 
# # # #   primary beam at X-band = 5.4', 
# # # #   mosaic field spacing = 2.5'
# # # #   total mosaic size = approx. 9.5' = 570"
# # # #   synthesized beam size = 8.4" in D config at 3.6 cm, 8.3 GHz
# # # #   cell size = 2" and nx,ny = 300 (600" field size)
# # # #   phase center = center field: J2000 05:35:17.470, -005.23.06.790
# # # #   NOTE: field 10 is outside of the 9 point primary mosaic (sitting 
# # # #    on M43 -- but the flux is resolved out so there is no use to 
# # # #    add it to the mosaic.  The script below leaves it out.  
# # # # 
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # 
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # 

include 'imagerwizard.g';
myimager:=imager(filename="orion.ms" );
# select which fields and IFs to image in the mosaic: 
ok:=myimager.setdata(mode="none" , nchan=1, start=1, step=1, 
	mstart=[value=0.0, unit="km/s" ], 
	mstep=[value=0.0, unit="km/s" ], spwid=[1, 2] , 
	fieldid=[3, 4, 5, 6, 7, 8, 9, 10, 11] , msselect='');
# set cell and image size, reference fieldid = 7, mosaic center:
ok:=myimager.setimage(nx=300, ny=300, cellx=[value=2.0, unit="arcsec" ], 
	celly=[value=2.0, unit="arcsec" ], stokes="I" , doshift=T, 
	phasecenter=[m0=[value=83.8227917, unit="deg" ], 
	m1=[value=-5.38521944, unit="deg" ], type="direction" , 
	refer="J2000" ], shiftx=[value=0.0, unit="arcsec" ], 
	shifty=[value=0.0, unit="arcsec" ], mode="mfs" , nchan=1, 
	start=1, step=1, mstart=[value=0.0, unit="km/s" ], 
	mstep=[value=0.0, unit="km/s" ], spwid=1, fieldid=7, facets=1);
  # 
  # NOTE, the current AIPS++ algorithm can only correctly mosaic one IF 
  #       at a time.  Thus, at the present time, you must create 2 images
  #       (one for each spwid) and then average them together in the image
  #       plane.  This script only images a single IF. 
  #
# set the voltage pattern for the VLA antennas:
ok:=myimager.setvp(dovp=T, usedefaultvp=T, vptable='', 
	dosquint=F, parangleinc=[value=360.0, unit="deg" ]);
# use robust weighting: 
ok:=myimager.weight(type="briggs" , rmode="norm" , 
	noise=[value=0.0, unit="Jy" ], robust=-1, 
	fieldofview=[value=0.0, unit="rad" ], npixels=0);
  # Weighting MS: IMAGING_WEIGHT column will be changed
  # Briggs weighting: sidelobes will be suppressed over full image
  #                  : using 300 pixels in the uv plane
# NOTE: use defaults on cyclefactor and cyclespeedup so no need to 
#  run the setup function myimager.setmfcontrol
#
# Run MEM with full image allowed in clean box:
ok:=myimager.mem(algorithm="mfentropy" , niter=30, 
	sigma=[value=0.04, unit="mJy" ], 
	targetflux=[value=1.0, unit="Jy" ], constrainflux=F, 
	displayprogress=T, model="orion.mem" , fixed=F, 
	complist='', prior='', mask='' , 
	image="orion.mem.restored" , residual="orion.mem.residual" );
# Now, create a mask for MEM clean region:  
# How to make a mask:
# 	- Set up imager as you need for mosaicing
# 	- Use imager.regionmask 
# 	- At the region input: 
#         from the wrench select "From image": you'll be asked for an
# 	  image to display. Select that as usual.
# 	- On the resulting viewer, make a polygon region and doubleclick.
# 	- It will then show up in the inputs for regionmask as "myregion".
# 	- Set a name for the mask (orion.mask) and press GO. When it
finishes,
# 	  you can use the wrench-View to verify that the mask is 
# 	  as you expect.
# 	- Run imager.mem again with mask set to orion.mask.
# note: region below not correct syntax, scripter won't save correctly.
ok:=myimager.regionmask(mask="orion.mask4" , region=myregion , value=1 );
ok:=myimager.mem(algorithm="mfentropy" , niter=100, 
	sigma=[value=0.04, unit="mJy" ], 
	targetflux=[value=1.0, unit="Jy" ], constrainflux=F, 
	displayprogress=T, model="orion.mem2" , fixed=F, 
	complist='', prior='', mask='orion.mask' , 
	image="orion.mem2.restored" , residual="orion.mem2.residual" );
  # watch the display progress plot and logger window results as MEM runs.  

# view image with the wrench-View on orion.mem2.restored: 
# Get image statistics: 
  #	flux = 179.5 Jy
  #	max  = 0.8245 Jy/beam
  #	min  = -0.044 Jy/beam

# view the residual image with the wrench-View on orion.mem2.residual:
# Get image statistics:
  #	rms  = 0.0081 Jy/beam = 8.1 mJy/beam
  #	mean = -2.0e-4 Jy/beam
  #	min  = -4.4e-2
  #	max  = 4.0e-2
  # This is a complicated emission region.  If you find patches of
emission 
  # outside of the region you selected, then remake the regionmask
  # using the residual image, orion.mem2.residual.  Call the mask 
  # orion.mask2 and rerun MEM.  

ok:=myimager.mem(algorithm="mfentropy" , niter=100, 
	sigma=[value=0.04, unit="mJy" ], targetflux=[value=1.0,
unit="Jy" ], 
	constrainflux=F, displayprogress=T, model="orion.mem3" , fixed=F, 
	complist='', prior='', mask="orion.mask2" , 
	image="orion.mem3.restored" , residual="orion.mem3.residual" );

# view residual image to make sure new mask covers all emission.  
  #	statistics:
  #	rms	= 6.1 mJy/beam
  #	mean	= -2.7e-4 Jy/beam
# view restored image: 
#  note: Use: Adjust resampling mode = bilinear to make a smooth image.  
  #	statistics:
  #	flux	= 189.2 Jy
  #	max	= 0.834 Jy/beam
  #	min	= -2.775e-2 Jy/beam
  # save this image to a ps file (using print button on viewer).  
  # You're done.