Difference between revisions of "Fit a Gaussian to Visibility data and plot it over the data"

From CASA Guides
Jump to navigationJump to search
(Created page with " *DRAFT* B.Mason, K.Ward-Duong, & J.Patience This CASA guide illustrates how to run UVMODELFIT to fit a single Gaussian component to UV data (i.e. the raw visibilities); and...")
 
Line 3: Line 3:
  
 
B.Mason, K.Ward-Duong, & J.Patience
 
B.Mason, K.Ward-Duong, & J.Patience
 +
July 2015
 +
developed in CASA 4.2.2
  
This CASA guide illustrates how to run UVMODELFIT to fit a single Gaussian component to UV data (i.e. the raw visibilities); and subsequently how to make a plot of the binned, "azimuthally" averaged (in UV space) UV data together with the fitted model.
+
This CASA guide illustrates how to run UVMODELFIT to fit a single Gaussian component to UV data (i.e. the raw visibilities); and subsequently how to make a plot of the binned, "azimuthally" averaged (in UV space) UV data together with the fitted model. In broad outlines, the steps required are
 +
  1 Estimate the peak location of emission. shift the phase center here if desired.
 +
  1 Estimate starting values for the gaussian component fit
 +
  1 do the fit.
 +
  1 run ft() to put the results of the fit in the MS explicitly (as a MODEL column).
 +
  1 Extract MODEL and DATA, make the plot
 +
 
 +
PHASE SHIFT:
 +
<pre>
 +
fixvis(vis="J04292165_cont.ms",outputvis="J04292165_cont_recenter.ms",field="",refcode="",reuse=True,phasecenter="J2000 4h29m21.66s +27d01m25.72s",datacolumn="all")
 +
</pre>
 +
(note: not done in this example)
 +
 
 +
do the fit -- this is for all spws together:
 +
<pre>
 +
default uvmodelfit
 +
vis=myvis
 +
field='0'
 +
comptype='G'
 +
sourcepar=[0.006,0.1,-0.05,0.4,0.3,0.0]
 +
varypar=[]
 +
outfile='J04292165_cont_allspw.uv.cl'
 +
inp
 +
go
 +
</pre>
 +
in this example we had not phase shifted.
 +
 
 +
run FT -- you must have usescratch=True
 +
<pre>
 +
default ft
 +
vis=myvis
 +
complist='J04292165_cont_allspw.uv.cl'
 +
usescratch=True
 +
inp
 +
go
 +
 
 +
# did that work?
 +
plotms(vis=myvis,xaxis='uvwave',yaxis='real',ydatacolumn='model')
 +
# yes it did.
 +
</pre>
 +
 
 +
Finally, use VISSTAT() to extract model and data values in bins of uv-radius:
 +
<pre>
 +
import numpy as np
 +
import matplotlib.pyplot as plt
 +
import glob
 +
 
 +
testms = 'J04292165_cont.ms'
 +
 
 +
# set the uvrange in kilo lambda
 +
uvmin = 0
 +
#uvmax = 1200
 +
uvmax = 800
 +
# define steps in klambda
 +
duv=40
 +
 
 +
uvsteps = np.arange(uvmin, uvmax, duv)
 +
 
 +
avg_amps = [] # define an empty list in which to put the averaged amplitudes
 +
 
 +
stddev_amps = [] # define an empty list in which to put the amp stddev
 +
 
 +
numpoints = []
 +
 
 +
uvpoints = uvsteps[:-1] + np.diff(uvsteps)/2 # get array of uvmidpoints over which avg taken
 +
 
 +
model_amps = [] # define list to get model amplitudes after fourier transforming into the model data column
 +
 
 +
# iterate over the uvrange:
 +
for ii in range(len(uvsteps[:-1])):
 +
        tmp = visstat(vis = testms,axis = 'real',
 +
                      uvrange = str(uvsteps[ii]) + '~' + str(uvsteps[ii+1]) + 'klambda',
 +
                      datacolumn = 'data')
 +
        print str(uvsteps[ii]) + '~' + str(uvsteps[ii+1]) + 'klambda'
 +
        avg_amps.append(tmp['DATA']['mean'])
 +
        stddev_amps.append(tmp['DATA']['stddev'])
 +
        numpoints.append(tmp['DATA']['npts'])
 +
 
 +
 
 +
for ii in range(len(uvsteps[:-1])):
 +
        tmp = visstat(vis = testms,
 +
                        axis = 'real', # you may want to change this to real...?
 +
                        uvrange = str(uvsteps[ii]) + '~' + str(uvsteps[ii+1]) + 'klambda',
 +
                        datacolumn = 'model')
 +
        print str(uvsteps[ii]) + '~' + str(uvsteps[ii+1]) + 'klambda, from the model column'
 +
        #avg_amps.append(tmp['DATA']['mean'])
 +
        #stddev_amps.append(tmp['DATA']['stddev'])
 +
        #numpoints.append(tmp['DATA']['npts'])
 +
        model_amps.append(tmp['MODEL']['mean'])
 +
 
 +
error_amps = stddev_amps/(np.sqrt(numpoints) - 1)
 +
 
 +
plt.clf()
 +
plt.cla()
 +
 
 +
avg_amps_mjy = [x * 1000 for x in avg_amps]
 +
model_amps_mjy = [x * 1000 for x in model_amps]
 +
error_amps_mjy = [x * 1000 for x in error_amps]
 +
 
 +
plt.errorbar(uvpoints, avg_amps_mjy, yerr = error_amps_mjy, mfc='k', fmt='o', label='data')
 +
plt.plot(uvpoints, model_amps_mjy, 'r-', label=r'model from $uvmodelfit$')
 +
 
 +
plt.legend(loc = 3, numpoints = 1)
 +
 
 +
plt.xlabel(r'UV Distance (k$\lambda$)')
 +
plt.ylabel('Real Visibility (mJy)')
 +
 
 +
plt.show()
 +
</pre>

Revision as of 13:22, 24 July 2015

  • DRAFT*

B.Mason, K.Ward-Duong, & J.Patience July 2015 developed in CASA 4.2.2

This CASA guide illustrates how to run UVMODELFIT to fit a single Gaussian component to UV data (i.e. the raw visibilities); and subsequently how to make a plot of the binned, "azimuthally" averaged (in UV space) UV data together with the fitted model. In broad outlines, the steps required are

  1 Estimate the peak location of emission. shift the phase center here if desired.
  1 Estimate starting values for the gaussian component fit
  1 do the fit.
  1 run ft() to put the results of the fit in the MS explicitly (as a MODEL column).
  1 Extract MODEL and DATA, make the plot

PHASE SHIFT:

fixvis(vis="J04292165_cont.ms",outputvis="J04292165_cont_recenter.ms",field="",refcode="",reuse=True,phasecenter="J2000 4h29m21.66s +27d01m25.72s",datacolumn="all")

(note: not done in this example)

do the fit -- this is for all spws together:

default uvmodelfit
vis=myvis
field='0'
comptype='G'
sourcepar=[0.006,0.1,-0.05,0.4,0.3,0.0]
varypar=[]
outfile='J04292165_cont_allspw.uv.cl'
inp
go

in this example we had not phase shifted.

run FT -- you must have usescratch=True

default ft
vis=myvis
complist='J04292165_cont_allspw.uv.cl'
usescratch=True
inp
go

# did that work?
plotms(vis=myvis,xaxis='uvwave',yaxis='real',ydatacolumn='model')
# yes it did.

Finally, use VISSTAT() to extract model and data values in bins of uv-radius:

import numpy as np
import matplotlib.pyplot as plt
import glob

testms = 'J04292165_cont.ms'

# set the uvrange in kilo lambda
uvmin = 0
#uvmax = 1200
uvmax = 800
 # define steps in klambda
duv=40

uvsteps = np.arange(uvmin, uvmax, duv)

avg_amps = [] # define an empty list in which to put the averaged amplitudes

stddev_amps = [] # define an empty list in which to put the amp stddev

numpoints = []

uvpoints = uvsteps[:-1] + np.diff(uvsteps)/2 # get array of uvmidpoints over which avg taken

model_amps = [] # define list to get model amplitudes after fourier transforming into the model data column

# iterate over the uvrange:
for ii in range(len(uvsteps[:-1])):
        tmp = visstat(vis = testms,axis = 'real',
                      uvrange = str(uvsteps[ii]) + '~' + str(uvsteps[ii+1]) + 'klambda',
                      datacolumn = 'data')
        print str(uvsteps[ii]) + '~' + str(uvsteps[ii+1]) + 'klambda'
        avg_amps.append(tmp['DATA']['mean'])
        stddev_amps.append(tmp['DATA']['stddev'])
        numpoints.append(tmp['DATA']['npts'])


for ii in range(len(uvsteps[:-1])):
        tmp = visstat(vis = testms,
                        axis = 'real', # you may want to change this to real...?
                        uvrange = str(uvsteps[ii]) + '~' + str(uvsteps[ii+1]) + 'klambda',
                        datacolumn = 'model')
        print str(uvsteps[ii]) + '~' + str(uvsteps[ii+1]) + 'klambda, from the model column'
        #avg_amps.append(tmp['DATA']['mean'])
        #stddev_amps.append(tmp['DATA']['stddev'])
        #numpoints.append(tmp['DATA']['npts'])
        model_amps.append(tmp['MODEL']['mean'])

error_amps = stddev_amps/(np.sqrt(numpoints) - 1)

plt.clf()
plt.cla()

avg_amps_mjy = [x * 1000 for x in avg_amps]
model_amps_mjy = [x * 1000 for x in model_amps]
error_amps_mjy = [x * 1000 for x in error_amps]

plt.errorbar(uvpoints, avg_amps_mjy, yerr = error_amps_mjy, mfc='k', fmt='o', label='data')
plt.plot(uvpoints, model_amps_mjy, 'r-', label=r'model from $uvmodelfit$')

plt.legend(loc = 3, numpoints = 1)

plt.xlabel(r'UV Distance (k$\lambda$)')
plt.ylabel('Real Visibility (mJy)')

plt.show()