Primary beam correction

Is your feature request related to a problem or opportunity? Please describe. Synthesized images are not yet corrected for primary beam sensitivity.

In a CLEAN framework, the primary beam is "corrected" for after the deconvolution process is complete, by division with a normalized a primary beam profile.

In an RML framework, I think we want to treat the primary beam as part of the forward modeling process.

If this is how visibilities are generated from a sky brightness distribution,

Then the primary beam profile just needs to be multiplied against the surface brightness distribution before taking the FFT.

For smaller-extent images (most protoplanetary disks), the pbcor isn't mission critical and can be ignored if one is just looking at morphologies. To get correct fluxes, though, the pbcor is important.

Update: This ALMA knowledge-base article has helpful leads for characterizing the ALMA primary beam. There are varying degrees of approximation that we can do here. The simplest solution, which we should implement first, is to divide by the Airy disk calculated for a 12m illumination pattern. Section 3.2 of the ALMA Technical Handbook describes the Single dish response (Fig 3.3). Strangely enough, I don't see the actual functional form in the Technical Handbook. Wikipedia has notes on the Airy pattern, which we should try to check. We can also investigate the actual functional form used by CASA by using tclean to make both a pbcor'd image and a non-pbcor'd image, and then we can divide the two.

Mathematical description

@kdesoto-psu, could you and your group members reply to this issue, updating with your current thinking re: the mathematical form of the primary beam profile for an ALMA 12 meter dish?

Implementation

The primary beam should be implemented as a new PyTorch "layer," probably as part of the Images module: https://mpol-dev.github.io/MPoL/api.html#module-mpol.images

When initialized, the primary beam layer is to calculate the normalized primary beam pattern (i.e. 1 at the center) over the extent of the image. To do this, the necessary information I think you would need as arguments includes

some coords object (to get image size)
diameter of the antenna, if using a uniformly illuminated circular aperture, should default to ALMA 12m antennas
two diameters if using a secondary blockage illumination pattern, should default to ALMA 12m antennas
frequency of the observations in Hz. This argument can either be a float or an array of floats. You might want to use np.at_least_1d to handle these arguments, since you'll want to store both internally as a length (nchan) array.

The init routine of this layer should then

use the length of the frequency array to calculate nchan
use the other quantities to fill up an array of size (nchan, npix, npix) with the primary beam transmission values, and store this to the layer as a variable using self, e.g., self.primary_beam_profile

Then, the forward method of the layer should take in an ImageCube, multiply it self.primary_beam_profile and then instantiate and return a new ImageCube that has the primary beam applied.

If you really wanted to be thorough, you could subclass ImageCube to create a PrimaryBeamAppliedImageCube (or something like that), instantiate that with your result, and then return that.

Finally, once the Primary Beam layer is complete, it can be added as part of a new "Precomposed Module" like SimpleNet. The primary beam correction would come after the ImageCube and before the FourierCube: https://github.com/MPoL-dev/MPoL/blob/main/src/mpol/precomposed.py#L50

Tests

We like test driven development at MPoL, so please outline some of the tests that make sense for a primary beam correction layer.

@iancze Here's the current status:

Implementation We have begun creating a cube in images.py called PBCorrectedCube (if you prefer another name let me know), which goes between ImageCube and FourierCube in the SimpleNet. It has the same functions as ImageCube to convert to a sky cube and FITS files, and a forward method for easy integration in the SimpleNet. The current arguments are coords, nchan, and telescope, which I'm currently setting to ALMA_12 for the ALMA 12 in. telescope. The init() function calculates the primary beam correction, so it should take in frequency channel information instead of nchan, and then just set self. nchan=len(freq_channels) within the init function. The forward() function takes the cube as an argument, applies the PB correction with a simple multiplication, and then returns the corrected cube.

The SimpleNet would then have to take in the channel frequencies instead of nchan as an argument, to initialize its PBCorrectedCube attribute correctly. To test this, however, it would be nice to know how to import the channel frequencies from CASA data with the rest of the data import from the tutorial, which I currently do not know how to do. Even though test cases can be created without this import, as the data import itself is outside of MPoL, I want to know the format and expected value ranges of these frequency imports for successful implementation.

Update 10/28: The frequencies are available as part of the examine.get_processed_visibilities() dictionary that is used to import the real visibility data examples in mpoldatasets. We will generate a separate Zenodo with .npz files with this freq info included to test our updated pipeline on.

The current primary beam correction I'm implementing is an Airy disk with 10.7m beam diameter, with a 0.75m blockage. The form of the Airy disk depends on both this beam radius and the wavelength lambda of the channel. All image values within the blockage will be set to 0.

Tests A simple test that adjusts an image of all 1s to a PB-corrected image (in this case a partially blocked Airy disk) would be good to include. For completion, one that ensures an array of all 0s converts to still all 0s should be included. If channel frequency data is also included, we would need to add tests that check for proper Airy disk generation for different channel frequencies.

Changes can be tracked here: https://github.com/kdesoto-psu/MPoL

Re: tests: I like all of these ideas! You might want to add a stub for them in the tests/ folder. It probably makes sense to create an entirely new test file like primary_beam_test.py that contains these ideas. You will probably want to take advantage of some of the pytest fixtures we created, too.

Re: frequencies: The frequencies aspect of MPoL is something that we are treating as an additional bit of information, rather than a core bit of information that would be an attribute of a particular object. I think the main reason for this is that there exist a large set of useful RML applications (like non-PB corrected continuum observations) that do not require frequency information, and so it's just easier to be conceptually minimalist (and abstract) about it for the sake of modular simplicity.

Moreover, even some of our image-cube oriented routines don't necessarily require the channel frequency information. If/when we eventually develop regularizers for image cubes, we'd probably be more interested in the velocities of each channel rather than their individual frequencies, and therefore just ask the user to provide velocities directly.

For the PB project and applications, you should feel free to make up your own frequency array. An example for a Band 6 ALMA continuum observation might be something like 8 channels, with chan0 = 220.00 GHz with 0.5 GHz channel spacing. You might be able to detect subtle differences for the primary beam within each channel.

But you are right, for a real data application we will need to be clear about how to propagate the frequency information from the CASA measurement set to a raw numpy array and then into MPoL.

Re: telescope type. Rather than getting specific about ALMA, I think it makes sense to have an abstracted Primary Beam type. You could have a "Uniform dish" type which takes as an argument the dish diameter in meters, and then a CentralyBlockedDisk (or better name) that takes in two arguments, the dish diameter and the diameter of the central blockage.

Thanks for the nice work! I'll provide some more comments on your fork.

MPoL-dev / MPoL