Fit a Gaussian to Visibility data and plot it over the data

From CASA Guides
Revision as of 08:51, 28 July 2015 by Bmason (talk | contribs)
Jump to navigationJump to search
  • 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.
  2 Estimate starting values for the gaussian component fit
  3 do the fit.
  4 run ft() to put the results of the fit in the MS explicitly (as a MODEL column).
  5 Extract MODEL and DATA, make the plot

The optional first step is to shift the phase center of the visibility data to the peak of emission in the map. This will maximize the amplitude of the azimuthally averaged profile in UV space, but is not strictly necessary for the results to be correct. This can be done using FIXVIS() --

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

(note: we have not actually done this stip in this example)

Next we fit a Gaussian component to the visibilities using UVMODELFIT(). Note that UVMODELFIT is currently restricted to fitting a single component. It is important that the

default uvmodelfit

again, in this example we had not phase shifted to the peak; otherwise we would have put "0,0" for the initial component center given by the "0.1,-0.05" [arcsec] entries of the SOURCEPAR parameters.

run FT -- you must have usescratch=True

default ft

# did that work?
# 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 = ''

# set the uvrange in kilo lambda
#  NB: take care that the bin width and uvmin/max ranges ensure
#   you get no bins without data! the appropriate values may change
#   if you, e.g., change the SPW you are looking at
uvmin = 0
#uvmax = 1200
uvmax = 800
 # define steps in klambda

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 & build up
#  the binned data values (looking only at the Real part -- a pretty
#  good approach if the signal is approximately symmetric and you have
#  recentered the phase center on the peak of emission):
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'

# now build up the model values
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'

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


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)')

The result--