Determine spectral slope of science source

Work in progress! (12/07/2011, Shaye)

These notes are from a project to determine the spectral slope of DG Tau, which is a protostellar object.

We assume your data has been properly flagged.

There are 16 continuum windows (500 MHz BW) in this example.

Determine Spectral Index of Passband Calibrator

set flux=URANUS

set pass=3C84

set vis=uvdata

set bands = ( 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 )

#Brightness temperature of URANUS at our observing frequencies for flux calibration

#Derived using function provided by A. Isella. He derived it from WMAP data (and references therein).

   #Tb=b1 + b2*log(x) + b3*log(x)^2 + b4*log(x)^3
   #b1=154.813 b2=106.598 b3=-91.0474 b4=15.0796

set uranustb = ( 125.02683 125.20088 125.37588 125.55175 125.72856 125.90630 126.08495 126.26458 121.98721 121.82910 121.67174 121.51508 121.35952 121.20437 121.04997 120.89625 )

#C8 at center of this D-array

set refant=8

#Prepare data

rm -Rf $vis uvcat vis=c0800.4D_102HH211.2.miriad out=$vis select="-auto" options=nopass,nocal,nopol

#Flag visibility data

foreach band ( $bands )

 #Select individual window
 /bin/rm -Rf $vis.$band
 uvcat vis=$vis out=$vis.$band select="'win($band)'" options=nocal,nopol,nopass

 #Run mfcal on this window on very short timescale to minimize decorrelation
 mfcal vis=$vis.$band select="source("$pass")" interval=0.1,10000. refant=$refant edge=2 line=chan,39,1,1,1

 #Check mfcal solution
 #gpplt vis=$vis.$band yaxis=amp device=1/xs nxy=4,4  
 #set text = "Press ENTER .."
 #if ($#argv > 0) set text = "$* >>" 
 #printf "\n$text" 
 #set junk = ($<)

 #Apply mfcal solution to correct for variation across the bandpass
 /bin/rm -Rf $vis.$band.pb 
 uvcat vis=$vis.$band out=$vis.$band.pb options=nocal,nopol

 #STARTING STEPS FOR BOOTFLUX 

   #BANDPASS SOURCE PHASE CALIBRATION
   #Run selfcal on bandpass source do derive gain solutions for the bandpass source
   #Interval can be short since just applying it to the same source (?)
   mselfcal vis=$vis.$band.pb select="source("$pass")" interval=1.01 options=apriori,noscale,phase refant=$refant line=chan,1,1,39,39
   
   #Check selfcal solution for bandpass source
   
   #Apply selfcal phase solution from the bandpass calibrator to the bandpass calibrator
   /bin/rm -Rf $vis.$band.$pass.scp
   uvcat vis=$vis.$band.pb out=$vis.$band.$pass.scp select="source("$pass")" options=nopass,nopol
   
   #FLUX SOURCE PHASE CALIBRATION
   #Run selfcal on flux source to derive gain solutions for the flux source
   #Interval can be short since just applying it to the same source (?)
   mselfcal vis=$vis.$band.pb select="source("$flux")" interval=1.01 options=apriori,noscale,phase refant=$refant line=chan,1,1,39,39

   #Check selfcal solution for flux source

   #Apply selfcal phase solution from the flux calibrator to the flux calibrator
   /bin/rm -Rf $vis.$band.$flux.scp
   uvcat vis=$vis.$band.pb out=$vis.$band.$flux.scp select="source("$flux")" options=nopass,nopol

   #UVCAT bandpass and phase corrected visiblity data for the BANDPASS and FLUX calibrator together
   #This is the visibility file to use in bootflux
   /bin/rm -Rf $vis.$band.btrdy
   uvcat vis=$vis.$band.$pass.scp,$vis.$band.$flux.scp out=$vis.$band.btrdy options=nopass,nopol,nocal

   #BOOTFLUX STEP
   #print out Tb of URANUS at this window
   echo $uranustb[$band]
   #run bootflux and dump output into log file
   bootflux vis=$vis.$band.btrdy select="source("$pass,$flux")" line=chan,1,1,39,39 taver=5.01 primary=$flux,$uranustb[$band] badres=70 log=bootfluxJ.$pass.log2
   #grep that log file and extract the bootstrapped Average Flux and Flux Error of 3C84 at this window
   grep Average bootfluxJ.$pass.log2 | grep Flux | cut -c 17-31 > script1_bootflux.$pass.$band
   #print the result to the screen
   set tempfluxpass = `awk 'NR==1' script1_bootflux.$pass.$band`
   echo $tempfluxpass
   
   #clean up
   /bin/rm -Rf $vis.$band $vis.$band.pb $vis.$band.$pass.scp $vis.$band.$flux.scp $vis.$band.btrdy bootflux.$pass.log2

end

Determine Spectral Index of Phase Calibrator

set onsource=IRAS4

set flux=URANUS

set pass1=3C84

set pass2=3C84

set cal=3C84

set vis=uvdata

set bands = ( 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 )

set uranustb = ( 124.89169 125.05010 125.20929 125.36919 125.52988 125.69135 125.85355 126.01654 121.49650 121.35417 121.21234 121.07124 120.93065 120.79076 120.65140 120.51263 )

set refant=8

mfcal vis=$vis select="source("$pass1")" interval=0.1,10000. refant=$refant edge=2 line=chan,624,1,1,1 flux=1.868,103.0,-0.695 /bin/rm -Rf $vis.fullpb uvcat vis=$vis out=$vis.fullpb options=nopol,nocal

foreach band ($bands)

 #Select individual window
 /bin/rm -Rf $vis.$band
 uvcat vis=$vis.fullpb out=$vis.$band select="'win($band)'" options=nocal,nopol,nopass

 mselfcal vis=$vis.$band select="source("$cal")" interval=1.01 options=noscale,phase,apriori refant=$refant line=chan,1,1,39,39
 #smagpplt vis=$vis.fullpb yaxis=phase device=3/xs nxy=4,4
 /bin/rm -Rf $vis.$cal.$band.scp
 uvcat vis=$vis.$band out=$vis.$cal.$band.scp select="source("$cal")" options=nopass,nopol 

 mselfcal vis=$vis.$band select="source("$flux")" interval=1.01 options=noscale,phase,apriori refant=$refant line=chan,1,1,39,39
 #smagpplt vis=$vis.fullpb yaxis=phase device=3/xs nxy=4,4
 /bin/rm -Rf $vis.$flux.$band.scp
 uvcat vis=$vis.$band out=$vis.$flux.$band.scp select="source("$flux")" options=nopass,nopol 

 /bin/rm -Rf $vis.$band.btrdy
 uvcat vis=$vis.$cal.$band.scp,$vis.$flux.$band.scp out=$vis.$band.btrdy options=nopass,nopol,nocal

 bootflux vis=$vis.$band.btrdy select="source("$cal,$flux")" line=chan,1,1,39,39 taver=5.01 badres=70 primary=$flux,$uranustb[$band] log=bootfluxJ.$cal.log2
 grep Average bootfluxJ.$cal.log2 | grep Flux | cut -c 17-31 > script2_bootflux.$cal.$band
 set tempfluxcal = `awk 'NR==1' script2_bootflux.$cal.$band`
 echo $tempfluxcal
   
 /bin/rm -Rf $vis.$band $vis.$cal.$band.scp $vis.$flux.$band.scp $vis.$band.btrdy
 /bin/rm -Rf bootfluxJ.$cal.log2

end

/bin/rm -Rf $vis.fullpb

Determine Spectral Index of Science Source

set onsource=IRAS03 set flux=URANUS set cal=3C84 set vis=uvdata

set bands = ( 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 )

set uranustb = ( 125.027 125.201 125.376 125.552 125.729 125.906 126.085 126.265 121.987 121.829 121.672 121.515 121.360 121.204 121.050 120.896 )

set sources = ( IRAS03 B1c HH211 IRAS2A )

set refant=8

set fg=11.868

foreach onsource ( $sources )

   mfcal vis=$vis select="source("$cal")" interval=0.1,10000. refant=$refant edge=2 line=chan,624,1,1,1 flux=11.868,103.0,-0.695
   /bin/rm -Rf $vis.fullpb
   uvcat vis=$vis out=$vis.fullpb options=nopol,nocal

   mselfcal vis=$vis.fullpb select="source("$cal")" interval=18.01 options=noscale,phase,apriori refant=$refant line=chan,1,1,624,624
   /bin/rm -Rf $vis.scp
   uvcat vis=$vis.fullpb select="source("$cal,$onsource")" out=$vis.scp options=nopass,nopol

   mselfcal vis=$vis.scp select="source("$cal")" interval=40.01 options=noscale,amplitude refant=$refant line=chan,1,1,624,624 flux=$fg
   /bin/rm -Rf $vis.$onsource.scpa
   uvcat vis=$vis.scp select="source("$cal,$onsource")" out=$vis.$onsource.scpa options=nopass,nopol

   #Remove DGTAUB with uvmodel
   #/bin/rm -Rf $vis.scpam $vis.scpam2 $vis.calibrators $vis.scpamf
   #uvmodel vis=$vis.scpa out=$vis.scpam select="source("$onsource")" options=subtract flux=0.015 offset=-28.592,-45.934
   #uvmodel vis=$vis.scpam out=$vis.scpam2 select="source("$onsource")" options=subtract flux=0.0010 offset=-28.887,-43.999
   #uvcat vis=$vis.scpa out=$vis.calibrators select="source("$cal,$pass1,$pass2")" options=nopass,nocal,nopol
   #uvcat vis=$vis.calibrators,$vis.scpam2 out=$vis.scpamf options=nopass,nocal,nopol
   
   #Make image of IRAS2A to see peak flux is at expected phase center
   /bin/rm -Rf $vis.fullmap $vis.fullbm $vis.fullcln $vis.$onsource.scpamres
   invert vis=$vis.$onsource.scpa select="source("$onsource")" map=$vis.fullmap beam=$vis.fullbm imsize=257,257 cell=1,1 options=systemp robust=0 line=chan,1,1,624,624
   clean map=$vis.fullmap beam=$vis.fullbm out=$vis.fullcln gain=0.1 niters=2000
   restor model=$vis.fullcln beam=$vis.fullbm map=$vis.fullmap out=$vis.$onsource.scpamres

   #Is phase center of DGTau where expected?
   #imhead in=$vis.scpamres > script3_imhead
   #imfit in=$vis.scpamres object=gaussian region="rel,box(-20,-20,21,21)" > script3_imfit.tbl
   
   #extract ra and dec position of center position from imfit
   #grep Declination script3_imfit.tbl | cut -c 35-36 > script3_deccenter1
   #grep Declination script3_imfit.tbl | cut -c 38-39 > script3_deccenter2
   #grep Declination script3_imfit.tbl | cut -c 41-46 > script3_deccenter3
   #set deccenterval1 = `awk 'NR==1' script3_deccenter1` 
   #set deccenterval2 = `awk 'NR==1' script3_deccenter2`
   #set deccenterval3 = `awk 'NR==1' script3_deccenter3`
   #grep Ascension script3_imfit.tbl | cut -c 30-35 > script3_racenter1
   #grep Ascension script3_imfit.tbl | cut -c 37-38 > script3_racenter2
   #grep Ascension script3_imfit.tbl | cut -c 40-45 > script3_racenter3
   #set racenterval1 = `awk 'NR==1' script3_racenter1`
   #set racenterval2 = `awk 'NR==1' script3_racenter2`
   #set racenterval3 = `awk 'NR==1' script3_racenter3`
   #Correct phase center
   #/bin/rm -Rf temp.scpamf
   #uvedit vis=$vis.scpamf out=temp.scpamf source=$onsource ra=$racenterval1,$racenterval2,$racenterval3 dec=$deccenterval1,$deccenterval2,$deccenterval3 options=nouv
   #/bin/rm -Rf $vis.scpamf
   #uvcat vis=temp.scpamf out=$vis.scpamf options=nopass,nocal,nopol
   
   foreach band ($bands)

uvamp vis=$vis.$onsource.scpa select="source("$onsource"),win("$band")" bin=10,5,klam log=script3_uvamp.$onsource.$band uvamp vis=$vis.$onsource.scpa select="source("$cal"),win("$band")" bin=1,50,klam log=script3_uvamp.$cal.$band

   end
   
   /bin/rm -Rf $vis.fullpb $vis.scp 
   /bin/rm -Rf $vis.scpam $vis.scpam2 $vis.calibrators 
   /bin/rm -Rf $vis.fullmap $vis.fullbm $vis.fullcln $vis.scpamres
   /bin/rm -Rf temp.scpamf

end