How to manually extract the flux (in Jy) from a datacube of flux density (in Jy/beam)

I'm working with interferometric data from NOEMA. I have a datacube with intensity values that are in the units of "Jy/beam" (janskys per beam). My goal was to extract the flux (in regular janskys, or Jy) corresponding to the central source in these observations. This involved two hang-ups:

1. Selecting the pixels corresponding to the clean beam shape and size
2. Correctly calculating the beam area to match convention and get the proper number of janskys out.

I'm pretty sure that tools like CASA and gildas/mapping have already figured this stuff out, but for my purposes I wanted to extract the flux using my own tooling so that I could have more control (and have it be more reproducible); additionally, I'm not familiar enough with CASA or gildas to know where to start.

For issue #1, I needed to figure out how to select a rotated ellipse with major/minor axes corresponding to my beam. I found this math-stackexchange post super helpful in getting the Cartesian coordinates defining such an ellipse:

http://math.stackexchange.com/questions/426150/what-is-the-general-equation-of-the-ellipse-that-is-not-in-the-origin-and-rotate

In the following code snippet, you can see how I pulled this off:

https://github.com/tomr-stargazer/scratch-code/blob/master/ellipse_experiment.py

For Issue #2, I needed to know how exactly the beam area was defined, which wasn't obvious. I was told (by gildas/mapping and by the FITS header of my datacube) that the beam major and minor FWHM were 1.9 and 1.0 arcseconds, respectively. But what ellipse, exactly, defines the beam area? Is it the ellipse with major and minor axes equal to the FWHM listed previously, or was it an ellipse that used Gaussian sigmas as its major and minor axes? A FWHM and a Gaussian sigma (or "c") differ by the following factor:

(see https://en.wikipedia.org/wiki/Gaussian_function#Properties)

so the impact on the area calculation could be fairly large.

Guided by the following astropy documentation page:
http://docs.astropy.org/en/stable/api/astropy.units.equivalencies.brightness_temperature.html
I concluded that it was the ellipse with this area:

>>> bmaj = 15*u.arcsec
>>> bmin = 15*u.arcsec
>>> fwhm_to_sigma = 1./(8*np.log(2))**0.5
>>> beam_area = 2.*np.pi*(bmaj*bmin/fwhm_to_sigma**2)

(where I assume 'bmaj' and 'bmin' are the FWHM of the beam).
In the above example, there is a typo in the last line (already noted here: https://github.com/astropy/astropy/pull/4766 and ready to be fixed in the next astropy release) where the fwhm_to_sigma factor should be multiplied, not divided.

This page is probably also super helpful:
https://www.cv.nrao.edu/course/astr534/Interferometers2.html
especially this bit:

The beam solid angle of a Gaussian beam with HPBW 0 is

A=024ln2 

Ultimately I chose to extract flux corresponding to a "three-sigma" ellipse around the source, in order to capture the whole beam's flux.

I confirmed that my conversion was correct by looking at the integrated flux of a bright line in the "go view" mapping tool and comparing it to the integrated flux that I computed.

Blue: integrated flux (Jy) over the whole source size. Green: flux density (Jy/beam) at the central pixel.

gildas/mapping, showing basically the same as above, except that here the integrated flux is from a rectangular selection (the zoom window) rather than a more careful ellipse.