# Determine spectral slope of science source

(Difference between revisions)

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.

## Overview

1. Determine spectral slope of bandpass calibrator and its flux at some reference frequency.

2. Determine spectral slope of phase calibrator and its flux at some reference frequency.

3. Determine spectral slope of science sources using the above values to prevent any slope due to calibration process.

## Determine Spectral Index of Passband Calibrator

```set flux=URANUS # Uranus was our flux calibrator
set pass=3C84 # 3C84 was our passband calibrator
set vis=uvdata # this is our visibility data
set bands = ( 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 ) # these are the 16 bands
```
```#Brightness temperature of URANUS at our observing frequencies to be used during flux calibration command "bootflux"
#Derived using function provided by A. Isella. He derived it from WMAP data (and references therein).
#We hope this is more accurate than the models provided in Miriad.
#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 picked as refant since it is at the center of this D-array and looks good
set refant=8
```
```#Prepare data
rm -Rf \$vis
```
```#Flagging of visibility data already done
```
```#LOOP OVER ALL 16 BANDS
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
gpplt vis=\$vis.\$band yaxis=phase device=1/xs nxy=4,4

#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?
#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
```