This script goes through a data reduction in MIRIAD of CARMA data. It can be modified in order to fit the specifics of various observations - depending on what needs to be flagged, which calibrator should be the passband and flux calibrator, etc.
# ------------------------------------------------------------- # [B] Overview of Millimeter Wavelength radio data reduction # ------------------------------------------------------------- # There are only a couple basic steps that must be done # (0) baseline solutions - apply if online ones were wrong/not applied # Needed if slope seen in baseline-baseline # pairs. # (1) bandpass/passband calibration # -accomplish this with task mfcal # -bootflux to scale the calibrator's amplitude # or a more involved FLUX calibration here # Next two steps create GAINS tables that can be applied to the source # (2) phase calibration # (3) amplitude calibration # -self-cal on the calibrator # # (4) And Magic - the Fourier transform! # - task "invert" to create the map # # The rest of the commands are simply to apply calibration # to the desired source/calibrator, to look at the data, # etc. # # (5) Flux calibration (see step 1)
Set variables for the script
set vis = ct012.arp193.1.miriad set refant = 9 set source_name = Arp193 set bandpass_name = 3C273 set phasecal_name = 3C273 set fluxcal_name = 3C279 set outfile = ct012_arp193_feb5 set antpos = antpos.070115
# The following commands selects only the non-auto-correlation data # and re-writes to the file data.mir uvcat vis=$vis select=-auto out=data_auto.mir
Baseline Correction applied. In this example, data from June 9, 2007 were used and needd a baseline correction.
# See also: http://cedarflat.mmarray.org/observing/baseline/antpos.070509 if (0) then uvedit vis=data_auto.mir apfile=antpos.070115 out=baseline_data.mir else cp -r data_auto.mir baseline_data.mir endif
#----------------------------------------------------------- # Rest Frequency Correction #---------------------------------------------------------- # Do a uvlist vis=xxx.mir options=spc to see the rest frequency # and starting frequency of each channel. To put in the # proper rest frequency, do the following: # # uvputhd in=xxx.mir hdvar=restfreq varval=myrestfreq out=yyy.mir uvputhd vis=baseline_data.mir hdvar=restfreq varval=225.282 out=data.mir
# ---------------------------------------------------------- # [E] Preliminary Examination of Data # ---------------------------------------------------------- # Here we look only at the calibrators to just check on weather, # system temp if ($4 == <2) then uvindex vis=data.mir log=uvindex.log # scans uvdata file, keywords: vis, interval, refant, log, options listobs vis=data.mir log=listobs.log uvlist vis=data.mir options=spectra log=uvlist.log # What other variables are possible to plot? smavarplt vis=data.mir device=systemp.ps/cps nxy=2,3 yaxis=systemp options=compress # removing compress prints all baselines smavarplt vis=data.mir device=/xs nxy=2,3 yaxis=systemp options=compress yrange=0,1000 uvlist vis=data.mir options=spec uvflag vis=data.mir options=noapply flagval=flag #how many flags applied? #uvplt axis variables: time, dtime, amplitude, real, imag, phase, uu, vv, uvdistance, uvangle, hangle, dhangle, parang smauvplt vis=data.mir device=/xs axis=time,phase select='-source(Arp193),-source(noise)' #line=wide,1,1 smauvplt vis=data.mir device=/xs axis=time,amp select='-source(Arp193),-source(noise)' line=wide,1,1 # closure vis=$vis line=wide,1,1,1 # to look at closure solutions # options=nowrap # to make plots not wrap - such as phase # Putting the correct rest frequency into the header #uvputhd vis=data.mir hdvar=restfreq type=d varval=226.434 out=data2.mir endif
# ------------------------------------------------------------
# [F] Flagging noticeably Bad Data
# ------------------------------------------------------------
# options=noapply # does not apply to data
# Determine what needs to be flagged by carefully examining the
# data. Look at Tsys vs. time, Amplitudes and phases versus
# time, and the pointing.
# OVROs and C12 are out
uvflag vis=data.mir flagval=flag "select=ant(1,2,3,4,5,6,12)"
# This step does not seem necessary for this data set, flagging based
# on pointing, because examining data with smavarplt, yaxis=axisrms,
# I get that the absolute magnitude of variations never exceeds 1.
# BAZ - 7/13/2007
uvflag vis=data.mir flagval=flag "select=-pointing(0,3),-source(NOISE)"
# To flag antennas during a certain time range
# Flagging all during these three minutes of bandpass observation
# because there is a peak in amplitudes - triangular.
#uvflag vis=data.mir flagval=flag "select=time(07:15:00,07:18:00)"
# If I don't get rid of all OVROs, I have to do something special with the
# differing beam sizes if and only if I am doing a mosaic
# Unflag following to auto-flag bad antennas
#set badants = "" # bad antennas to flag
#uvflag vis=$vis select=anten'('$badants')' flagval=flag
# noticing in uvplt that baseline 12-14 and last scan is bad
#uvflag vis=data.mir flagval=flag "select=ant(12)(14)"
#uvflag vis=data.mir flagval=flag "select=time(11:00:00,11:15:00)"
#uvflag vis=data.mir flagval=flag "select=amp(30)"
#-----------------------------------------------------------------
# [G] Data Handling
# ----------------------------------------------------------------
# Split up the data. One can also leave it all together and use
# the appropriate select commands while executing various tasks.
# Or one can use the mega split program, ProjectExplode to separate
# by window and source.
uvcat vis=data.mir "select=source("$source_name")" out=source.mir
uvcat vis=data.mir "select=source("$bandpass_name")" out=bandpass.mir
uvcat vis=data.mir "select=source("$phasecal_name")" out=phasecal.mir
uvcat vis=data.mir "select=source("$fluxcal_name")" out=fluxcal.mir
# fluxcal testing testing
set fluxvis=bpcalib_fluxcal
set phasevis=fluxstep_phasecal
uvcat vis=data.mir "select=source("$fluxcal_name")" out=$fluxvis.mir
uvcat vis=data.mir "select=source("$phasecal_name")" out=$phasevis.mir
# ------------------------------------------------------------- # [H] STEP 1 -> Bandpass/Passband calibration # ------------------------------------------------------------- # super-wideband mfcal passband - 1 minute interval # I choose a 1 minute interval here because I have observed my bandpass # for a total of 10 minutes. And why? mfcal vis=bandpass.mir interval=1 refant=$refant # look at the bandpass data if ($1 < 2) then gpplt vis=bandpass.mir options=bandpass yaxis=phase nxy=3,2 yrange=-360,360 device=/xs gpplt vis=bandpass.mir options=bandpass yaxis=phase nxy=3,2 yrange=-360,360 device=bandpassSolution.ps/ps uvspec vis=bandpass.mir interval=15 options=nopass device=/xs nxy=3,3 axis=channel,amplitude uvspec vis=bandpass.mir interval=15 device=/xs nxy=3,3 axis=channel,amplitude uvspec vis=bandpass.mir interval=1000 device=/xs nxy=3,3 axis=channel,phase yrange=-180,180 endif # copy and apply bandpass calibrator solution to phase calibrator # -no gain calibration gpcopy vis=bandpass.mir out=phasecal.mir options=nocal uvcat vis=phasecal.mir options=nocal out=phasecal_bp.mir # copy and apply bandpass calibrator solution to the source, no cal gpcopy vis=bandpass.mir out=source.mir options=nocal uvcat vis=source.mir out=source_bp.mir # copy and apply bandpass calibrator to the flux calibrator gpcopy vis=bandpass.mir out=fluxcal.mir options=nocal uvcat vis=fluxcal.mir out=fluxcal_bp.mir # making second copy to have for the super calibrator sandwich try uvcat vis=source.mir out=source2_bp.mir
#------------------------------------------------------------------
# [H (continued)] STEP 1b -> Flux calibration
#------------------------------------------------------------------
# Cleaning up before starting this section again
rm -rf *.gain *.gains *.sw *.wide *.flux *.medians
# Perhaps set these ahead of time
#set flux = 8.4
set flux = 12.03
set calint = 0.2
set vcalint = 18
set fcalint = 1
set superwidewin = "1,2,3,4,5,6"
set superwidechan = "1,1,90"
set lsbfluxchan = "1,1,45,45"
set usbfluxchan = "1,46,45,45"
# Following all done using test files. No work done on the main files, yet.
mfcal vis=$fluxvis.mir interval=$calint refant=$refant
# Examining the data vs. freq and time
gpplt vis=$fluxvis.mir options=bandpass yaxis=phase nxy=3,2 yrange=-360,360 device=/xs
gpplt vis=$fluxvis.mir yaxis=phase nxy=3,2 yrange=-360,360 device=/xs
# Copying this derived bp solution from the flux calibrator to our phase calibrator
gpcopy vis=$fluxvis.mir out=$phasevis.mir options=nocal,nopol
uvcat vis=$phasevis.mir out=$phasevis.wide options=nocal
uvcat vis=$fluxvis.mir out=$fluxvis.gain options=nocal "select=win($superwidewin)"
selfcal vis=$fluxvis.gain refant=$refant interval=$fcalint "select=source($fluxcal_name)" options=noscale,amplitude,apriori flux=$flux
gplist vis=$fluxvis.gain options=zeropha,amp > $fluxvis.gains
# Unix pipe to get median values
grep Medians $fluxvis.gains | tr -d Medians: > $fluxvis.medians
cat $fluxvis.medians
# straightening out the phase of the phase calibrator - use a phase only selfcal with
# a fairly long integration time
uvcat vis=$phasevis.wide out=$phasevis.sw "select=win($superwidewin)"
selfcal vis=$phasevis.sw line=channel,$superwidechan interval=$vcalint options=phase refant=$refant
# Now the amplitude gains derived from the flux calibrator can be applied to the phase calibrator
# by replacing the amplitudes and keeping the phases determined from the selfcal solution
gplist vis=$phasevis.sw options=replace jyperk=@$fluxvis.medians
#The following special program, uvflux, can gather statistics on this phase calibrator
uvflux vis=$phasevis.sw options=nopol line=chan,$lsbfluxchan
uvflux vis=$phasevis.sw options=nopol line=chan,$usbfluxchan
uvflux vis=$phasevis.sw options=nopol > $phasevis.flux
# This value obtained by averaging the results form the LSB and USB above.
set visflux=11.26
echo "come back to this part"
# Finally, checking the time variance of the phase calibrator
echo "Checking the time variance of the phase calibrator"
uvcat vis=$phasevis.wide out=$phasevis.wide.gain
selfcal vis=$phasevis.wide.gain refant=$refant interval=$vcalint "select=source($phasecal_name)" \
options=noscale,amplitude flux=$visflux
gplist vis=$phasevis.wide.gain options=zeropha,amp > $phasevis.gains
# ----------------------------------------------------------------- # [I] STEPS 2 & 3 -> Amplitude & Phase Calibration # ----------------------------------------------------------------- # Create the GAINS tables # selfcal # You want the interval to be about equal to source-calibrator cycle # Use gplist to look at the time intervals (vis=file options=amp or phase) selfcal vis=phasecal_bp.mir interval=13 refant=$refant #selfcal vis=mastercalib_bp.mir interval=13 refant=13 gpplt vis=phasecal_bp.mir device=/xs yaxis=phase yrange=-720,720 gpplt vis=phasecal_bp.mir device=gainSolutions.ps/ps yrange=-180,180 # Copy self cal solution to # SOURCE # copy selfcal solution to source gpcopy vis=phasecal_bp.mir out=source_bp.mir options=nopass mode=copy echo "************************" echo "*** Apply flux gains to phase cal***" grep Medians $phasevis.gains | tr -d Medians: > $phasevis.medians cat $phasevis.medians gplist vis=phasecal_bp.mir options=replace jyperk=@$phasevis.medians echo "************************" echo "********Applying phase gains to source*********" grep Medians $phasevis.gains | tr -d Medians: > $phasevis.medians cat $phasevis.medians gplist vis=source_bp.mir options=replace jyperk=@$phasevis.medians # play with calibrator # cleaning up rm -rf phasecal.mp phasecal.bm phasecal.sl phasecal.cm phasecal.fits #linetype,nchan,start,width,step # step=width if you don't specify. # Try the defaults # options=systemp (will weight by the inverse of the system temp.) invert vis=phasecal_bp.mir map=phasecal.mp beam=phasecal.bm cell=.33 imsize=512 line=channel,1,1,90 #invert vis=weakcal_bp.mir map=weakcal.mp beam=weakcal.bm cell=.33 imsize=512 line=channel,1,1,90 clean map=phasecal.mp beam=phasecal.bm out=phasecal.sl niters=1000 restor map=phasecal.mp beam=phasecal.bm model=phasecal.sl out=phasecal.cm # look at some stats imstat in=phasecal.cm region=quarter imstat in=phasecal.cm region=box'(180,180,200,200)' ellint in=phasecal.cm # to look at the uv coverage fits in=phasecal.cm op=xyout out=phasecal.fits
# ----------------------------------------------------------------- # [J] FINAL step - invert # ----------------------------------------------------------------- # (1) # Check calibration on the weak calibrator # copy bandpass solution to the weak calibrator #gpcopy vis=bandpass.mir out=weakcal.mir options=nocal #uvcat vis=weakcal.mir out=weakcal_bp.mir # copy selfcal solution to weak calibrator #gpcopy vis=phasecal_bp.mir out=weakcal_bp.mir options=nopass mode=copy #rm -rf weakcal.mp weakcal.bm weakcal.sl weakcal.cm weakcal.fits #invert vis=weakcal_bp.mir map=weakcal.mp beam=weakcal.bm cell=0.2 imsize=512 #clean map=weakcal.mp beam=weakcal.bm out=weakcal.sl niter=1000 #restor map=weakcal.mp beam=weakcal.bm model=weakcal.sl out=weakcal.cm # (2) # play with the source rm -rf $source_name.mp $source_name.bm $source_name.sl $source_name.cm $source_name.fits $outfile.cm #invert vis=source_bp.mir map=arp220.mp beam=arp220.bm cell=0.2 imsize=512 line=channel,1,1,90 invert vis=source_bp.mir map=Arp193.mp beam=Arp193.bm cell=0.2 imsize=512 line=channel,1,1,90 sup=0 #ine=velocity,30,-4.617,43.450 clean map=Arp193.mp beam=Arp193.bm out=Arp193.sl niters=10000 restor map=Arp193.mp beam=Arp193.bm model=Arp193.sl out=$outfile.cm # Creating the cube rm -rf $source_name.cube.mp $source_name.cube.bm $source_name.cube.sl $source_name.cube.cm $source_name.cube.fits $outfile.cube.cm #invert vis=source_bp.mir map=arp220.mp beam=arp220.bm cell=0.2 imsize=512 line=channel,1,1,90 invert vis=source_bp.mir map=Arp193.cube.mp beam=Arp193.cube.bm cell=0.2 imsize=512 line=channel,90,1,1 sup=0 #ine=velocity,30,-4.617,43.450 clean map=Arp193.cube.mp beam=Arp193.cube.bm out=Arp193.cube.sl niters=10000 restor map=Arp193.cube.mp beam=Arp193.cube.bm model=Arp193.cube.sl out=$outfile.cube.cm #fits in=M80.cm op=xyout out=M80.fits #rm -rf arp193_2.mp arp193_2.bm arp193_2.sl arp193_2.cm arp193_2.fits #invert vis=source2_bp.mir map=arp193_2.mp beam=arp193_2.bm cell=0.2 imsize=512 #clean map=arp193_2.mp beam=arp193_2.bm out=arp193_2.sl niter=1000 #restor map=arp193_2.mp beam=arp193_2.bm model=arp193_2.sl out=arp193_2.cm #fits in=arp193_2.cm op=xyout out=arp193_2.fits # Then look at the final image in ds9 or some other FITS format # viwer. Statistics can be examined, etc. # ------------------------------------------------------------------ # [K] DS9 Viewing Notes # ------------------------------------------------------------------ # You can look at the *.*m maps in ds9 by doing a # mirds9 <filename> # once ds9 is open. # FRAME -> TILE to plot more than 1 # FRAME -> BLINK to blink back and forth.
# ------------------------------------------------------------------ # [L] Misc. Notes for LATER # ------------------------------------------------------------------ # # For CARMA array: # note if you have 3 beam sizes # mospsf # imfit # so restor does not use first beam size and apply for all # (this is not an issue for Arp 220 on 25 Apr 2007 because I only # had BIMA dishes)
7.-0