OK, this isn’t actually the ultimate resolved source cheatsheet … yet. I’m just trying to centralize my notes on dealing with extended sources. I’ve been doing this lately, and I had to rederive a bunch of things that I know I’ve figured out before. Hopefully, next time I’ll remember that I wrote all this stuff down here.
Units
Radio astronomical images typically come in two unit systems: Jy/pixel and Jy/beam. For various questions one usually wants to convert to Jy/arcsec² or Jy/sr. Jy/pixel is not trivially equivalent to Jy/arcsec², because the pixel size can change as a function of position in many map projections.
There must be something profound about why maps come in Jy/bm and models come in Jy/pixel, but I don’t see it right now.
Source fluxes are best measured as total fluxes which come in units of plain Jy. This should be obvious! Flux is conserved! In the general case peak fluxes are always dependent on the image properties.
The important exception, however, is unresolved sources which by their nature have a peak flux equal to their total flux no matter what unit system you’re using. Total flux of 1 Jy, peak flux of 1 Jy/px spread over 1 pixel, or peak flux of 1 Jy/bm spread out over 1 beam.
Resolved sources will have a peak flux that depends on the particular image units and projection, because the total flux will be spread out over a number of pixels that depends on the image particulars. In a Jy/px image, you need to figure out how many pixels the source subtends, and mutatis mutandis for Jy/bm. Generically, pkflux = totflux * arcsec2perpix / arcsec2insource. If you have a 2D Gaussian source where the axes are FWHMs in arcsec, arcsec2insource = 2 pi major minor / (8 ln 2).
I think the best way to work with Jy/bm is to go through Jy/px: the beam volume in pixels can be derived from the beam parameters and the pixel volume: pixperbm = arcsec2perbm / arcsec2perpix. For various operation I feel like this quantity will cancel out, though. Once again, due to the position-dependent pixel volume in most maps, this parameter is not constant across the map.
2D Gaussian Parametrizations
There are a few ways to parametrize a 2D Gaussian. I’ll omit centering and overall scaling in the following equations for simplicity.
Tersest. The simplest expression in terms of symbols is: 
This parametrization is computationally efficient but not intuitive. One particular danger to note is that the parameters behave a bit oddly, and there are domain limitations that make it awkward for numerical optimizers (which are likely to tweak the parameters outside of their domains): A and C must be negative, and, if I’m doing my algebra right, you must have B² < 4AC.
Mathematical Major/Minor/PA. We often want to express the shape of the 2D Gaussian as a major axis σ1, minor axis σ2, and position angle θ. Here σ1 ≥ σ2 > 0 and θ is the rotation between the major axis and +x (with +θ being towards +y). Because of the symmetry of the shape, there’s a 180° degeneracy in θ.
We find ![A = -\frac{1}{2}\left[\frac{\cos^2 \theta}{\sigma_1^2} + \frac{\sin^2 \theta}{\sigma_2^2}\right]](http://s.wordpress.com/latex.php?latex=A%20%3D%20-%5Cfrac%7B1%7D%7B2%7D%5Cleft%5B%5Cfrac%7B%5Ccos%5E2%20%5Ctheta%7D%7B%5Csigma_1%5E2%7D%20%2B%20%5Cfrac%7B%5Csin%5E2%20%5Ctheta%7D%7B%5Csigma_2%5E2%7D%5Cright%5D&bg=ffffff&fg=000000&s=0 "A = -\frac{1}{2}\left[\frac{\cos^2 \theta}{\sigma_1^2} + \frac{\sin^2 \theta}{\sigma_2^2}\right]")
, 
, and ![C = -\frac{1}{2}\left[\frac{\sin^2 \theta}{\sigma_1^2} + \frac{\cos^2 \theta}{\sigma_2^2}\right]](http://s.wordpress.com/latex.php?latex=C%20%3D%20-%5Cfrac%7B1%7D%7B2%7D%5Cleft%5B%5Cfrac%7B%5Csin%5E2%20%5Ctheta%7D%7B%5Csigma_1%5E2%7D%20%2B%20%5Cfrac%7B%5Ccos%5E2%20%5Ctheta%7D%7B%5Csigma_2%5E2%7D%5Cright%5D&bg=ffffff&fg=000000&s=0 "C = -\frac{1}{2}\left[\frac{\sin^2 \theta}{\sigma_1^2} + \frac{\cos^2 \theta}{\sigma_2^2}\right]")
.
Inverting, we find 
, 
, and 
. During a roundtrip, θ may flip by 180° depending on the numerics, but as mentioned above this doesn’t change anything.
Astronomical Major/Minor/PA. This is the same as the above, except that by the definition of PA, the x axis is declination and the y axis is right ascension. Be careful to write 
.
Numerically-friendly. If you’re fitting for Gaussian parameters, it’s very helpful to numerical stability to have nice continuous behavior over all of your parameters. The B parameters above seem to be problematic, so here’s a parametrization that is pretty cheap to compute and behaves nicely: 
, 
, and 
. It shouldn’t be hard to see that 
and 
.
Here, the optimizer may end up with either of the σ values being the larger. If σ2 > σ1, you can exchange them and add (or subtract) π/2 to θ.
Wikipedia. The Wikipedia page on Gaussian Functions uses: ![y = \exp\left(-\left[ax^2 + 2bxy + cy^2\right]\right)](http://s.wordpress.com/latex.php?latex=y%20%3D%20%5Cexp%5Cleft%28-%5Cleft%5Bax%5E2%20%2B%202bxy%20%2B%20cy%5E2%5Cright%5D%5Cright%29&bg=ffffff&fg=000000&s=0 "y = \exp\left(-\left[ax^2 + 2bxy + cy^2\right]\right)")
Here the condition on the parameters is that the matrix [a b \ b c] be positive definite, or that 
(or, equivalently I think, that a and c be positive and b² < ac).
Wikipedia gives a conversion formula to this system from major/minor/PA, but there are two caveats: the Wikipedia system names σ{1,2} as σ{x,y} somewhat misleadingly, and inverts θ to be a clockwise rotation from +x to -y. Avoid those equations since the ones I give above are clearer.
Gaussian Convolutions
Resolved sources are often modeled as 2D Gaussians. These are convolved with the 2D Gaussian of your synthesized beam, so one often wants to compute this convolution or deconvolution analytically. The file src/subs/gaupar.for in MIRIAD has routines for going both directions, and they are not trivial. The deconvolution case needs special handling for sources that are close to (or smaller than!) point sources. The existing heuristics seem only so-so. See mmm/pwilliams/ms/mssfextract for a Python deconvolution implementation.
To manage flux through such a deconvolution, work in total flux units. Nothing changes! The peak flux is scaled in two ways: once, to convert from Jy/bm to Jy/px (assuming this is what you’re doing), which involves the beam volume and the pixel volume; and once to convert the source area, which involves the original volume and the deconvolved volume.
Gaussian Coordinate Transformations
When you have a Gaussian model for a source, this usally needs to be rendered into actual pixels or something.
Doing this purely geometrically is difficult. For instance, in a given map projection, the RA/Dec axes may not be precisely aligned with the X/Y pixel axes. If you want to draw a Gaussian source as an ellipse, you need to rotate it a little bit. If the source is very large, these effects change over the source extent and things get scary. It’s better to calculate the RA/Dec of each pixel and brute-force your way through many of these issues.
To do this brute forcing, you need to find the value of the parameter that goes in the exponential. You can think of this as computing a distance in from the origin in the transformed space where the Gaussian is a circle. You un-translate, then un-rotate, then un-scale:
def gaussdist (x0, y0, maj, min, pa, x, y):
# x0, y0: Gaussian center
# maj, min: major/minor axes
# pa: position angle, from +x toward +y
dx, dy = x - x0, y - y0
c, s = N.cos (pa), N.sin (pa)
a = c * dx + s * dy
b = -s * dx + c * dy
d2 = (a / maj)**2 + (b / min)**2
# d2: "distance" from Gaussian
return d2
The astronomical “East from North” PA convention means that you can replace X with declination and Y with right ascension. This formulation can be turned into the one in the Wikipedia article on Gaussian Fuction but is a lot simpler to think about.
A: Converting Jy/pixel to Jy/arcsec²
To do this you need to calculate the size of a pixel at the position of the source you’re considering. Code to do this is in mmm/pwilliams/ms/msimgen. It looks something like this:
def pixelvolume (pixel2world, pixelcoords):
delta = 1e-5
w1 = pixel2world (pixelcoords)
pixelcoords[X] += delta
pixelcoords[Y] += delta
w2 = pixel2world (pixelcoords)
dra = w2[RA] - w1[RA]
ddec = w2[DEC] - w2[DEC]
return (dra**2 + ddec**2) / (2 * delta**2)
B: Error Ellipses From Covariant Variables
This is almost entirely unrelated to everything else here, except that it involves 2D Gaussians. Given two correlated variables, what are the parameters of the ellipse that describes a confidence region?
def ellpar (x, y, sx, sy, cxy, cl):
# x: expectation value of x variable
# y: expectation value of y variable
# sx: std dev (not variance) of x var
# sy: std dev (not variance) of y var
# cxy: covariance (not corr. coeff.) of x and y
# cl: desired confidence limit (0 < cl < 1)
from numpy import arctan2, log, sqrt
th = arctan2 (2 * cxy, sx**2 - sy**2)
rcl2 = -2 * log (1 - cl)
tmp1 = sqrt ((sx**2 - sy**2)**2 + 4*cxy**2))
tmp2 = sx2**2 + sy**2
maj = sqrt (0.5 * rcl2 * (tmp2 + tmp1)
min = sqrt (0.5 * rcl2 * (tmp2 - tmp1)
# x: same as input
# y: same as input
# maj, min: ellipse major minor axes as
# sigma values, not FWHMs
# th: position angle, rotating from +x to +y
return x, y, maj, min, th
This is derived by transforming X and Y to have a mean of 0 and stddev of 1, then considering the probability distribution of r = x²+y², which follows P(r) = r exp (-r²/2), and has CDF(r) = 1 – exp (-r²/2), allowing the r for a given confidence limit to be calculated.
![\frac{F}{2\pi\sigma_x\sigma_y} \exp\left(-\frac{1}{2}\left[\frac{x^2}{\sigma_x^2}+\frac{y^2}{\sigma_y^2}\right]\right)](http://s.wordpress.com/latex.php?latex=%5Cfrac%7BF%7D%7B2%5Cpi%5Csigma_x%5Csigma_y%7D%20%5Cexp%5Cleft%28-%5Cfrac%7B1%7D%7B2%7D%5Cleft%5B%5Cfrac%7Bx%5E2%7D%7B%5Csigma_x%5E2%7D%2B%5Cfrac%7By%5E2%7D%7B%5Csigma_y%5E2%7D%5Cright%5D%5Cright%29&bg=ffffff&fg=000000&s=0 "\frac{F}{2\pi\sigma_x\sigma_y} \exp\left(-\frac{1}{2}\left[\frac{x^2}{\sigma_x^2}+\frac{y^2}{\sigma_y^2}\right]\right)")
. Assuming circularity and simplifying, we get ![A\exp\left(B\left[x^2+y^2\right]\right)](http://s.wordpress.com/latex.php?latex=A%5Cexp%5Cleft%28B%5Cleft%5Bx%5E2%2By%5E2%5Cright%5D%5Cright%29&bg=ffffff&fg=000000&s=0 "A\exp\left(B\left[x^2+y^2\right]\right)")
, where 
and 
.
