Produce maps from a catalog of sources

[zonca]

We should implement a generic component that can take a list of coordinates and flux from an input file and generate a map already smoothed with the instrument beam.
This can for example be used for radio galaxies.

Looking for an existing implementation that we can port into PySM.

[seclark]

This is a great idea!

[giuspugl]

I have a tool handy that does that , can post a notebook and we can discuss about it.

[zonca]

@giuspugl yes, please, that would be very helpful, no hurry

[zonca]

@giuspugl @xzackli could you please share the code that you have to generate maps from a catalog of point sources?
@keskitalo do you have a pointer to the documentation of the new TOAST operator that handles point sources from a catalog?

I would like to gather all the code we already have, then come out with a plan on what we would like to port into PySM.

[giuspugl]

@zonca thanks for refreshing this issue, can surely point you to my code that does that, i can surely help with porting it into pysm3. The code assumes that all sources are unresolved and projected as point-sources with a gaussian profile.

• reads fluxes from a catalog
• converts fluxes into brightness temperature
• estimates polarization fluxes from the latest measurements of average pol. fraction
• projects sources as circles with size proportional to beam_solid_angle
• applies a gaussian profile to the source
• makes QU maps assuming uniformly distributed pol.angle over 0,pi .
  Note that the sources are projected with the beam already convolved, so we need to be careful to handle this
  see
  https://gitlab.com/giuse.puglisi/PS4C/-/blob/master/ps4c/catalogforecaster.py#L51

[keskitalo]

The Simons Observatory operator for simulating TOD from a catalog of sources is tersely documented here: https://sotodlib.readthedocs.io/en/latest/toast.html#reference. There is no anchor so you need to search for SimCatalog. The sources can be static, variable or even transient but their position must be held fixed to the Celestial sphere.

[zonca]

@xzackli 's code to generate maps from a catalog is https://github.com/WebSky-CITA/XGPaint.jl/blob/main/src/radio.jl#L254-L329, this puts all the flux for each source into a HEALPix pixel.

[zonca]

see also : https://github.com/simonsobs/pixell/blob/fa8e884d90127282209f6446c2233c319fb73c40/pixell/pointsrcs.py#L35
and
https://github.com/amaurea/tenki/blob/master/sim_websky.py

[zonca]

I think we can use pixell for this, see https://www.zonca.dev/posts/2024-05-02-pysm-point-source-pixell for a demo.
The low level functionality is implemented in C, it can simulate 10000 sources in half a second on Jupyter@NERSC.

[zonca]

next we need to decide the catalog format.

TOAST uses a TOML file which also includes frequency scaling for sources:

https://github.com/simonsobs/sotodlib/blob/071f2efc2dcd6cf33898a09e235fb004ac4e97b8/sotodlib/toast/ops/sim_catalog.py#L68-L82

@giuspugl I see websky has catalogs as a function of frequency: https://portal.nersc.gov/project/sobs/users/Radio_WebSky/matched_catalogs_2/, do you also have or it is possible to produce a single catalog that assumes a model (or some fit) for frequency dependence?

currently the catalogs have the columns: ['flux', 'phi', 'polarized flux', 'theta']

I think PySM should include a single Websky catalog for Radio Galaxies, already with a pre-defined flux cut to reduce size and with parameters for PySM to estimate on the fly the flux at each frequency. Then we can use pixell to produce beam-smoothed maps.

[xzackli]

Yes, it's possible to generate that kind of catalog for the bright sources. I think the easiest thing to do is actually to fit the existing sources with a log polynomial, since they have some moderately complicated spectral behavior, and we didn't save the original parameters which generated those spectra.

[zonca]

thanks @xzackli, that would be appropriate, would you be able to create that? For testing purposes even having a catalog with just 10 sources would be enough.

[giuspugl]

the plan sounds fine to me.
i think that 3rd or 4th order log-polynomial should encompass most of the fluxes we've sampled with websky.

[zonca]

I think I can give it a try and you can review it,
should I use the catalogs in https://portal.nersc.gov/project/sobs/users/Radio_WebSky/matched_catalogs_2/?

[zonca]

@msyriac as you have experience with other simulations, do you think "coordinates + coefficients of 4th order polynomial in log(freq) for flux" is ok as a standard representation of a catalog? Or do you have a better suggestion?

[msyriac]

What exactly do you mean by standard here? Existing radio catalogs from e.g. the Sehgal+ and Agora simulations are likely not in that format, but could be post-processed to be.

[zonca]

@msyriac we need to make a decision on the format of this catalog. do you agree with "coordinates + coefficients of 4th order polynomial in log(freq) for flux" or do you have something better in mind? I see that both Sehgal and Agora just give fluxes at some reference frequencies.