I'm wondering if there is a way to implement a simulation of a short duration targeted but high resolution observation and get the appropriate noise and sky models. For instance, if we did a short (hours) targeted scan of a small patch (like I would do if I were scanning a planet) it would have the appropriate noise and map depth etc. from the short scan and include the appropriate sky components (CMB, dust, etc). I'm particularly interested in producing several realizations. I can do this in Toast in TDS, but it is computationally intensive to generate and map multiple realizations over all optics tubes/frequencies. I see that MBS can take in hitmaps from e.g. toast (I have hit maps) but the "on the fly" tutorial seemed to indicate that it used the hit maps only for relative weighting, not for scaling the actual noise (which I think is set by years of observation etc.). Maybe MBS isn't the right way to do this and I should just keep using TDS?

There are broken links to documentation at the bottom of this page: https://github.com/simonsobs/map_based_simulations/blob/master/201911_lensed_cmb/README.md

See documentation at:

https://github.com/simonsobs/map_based_simulations/blob/master/mbs-s0012-20230321/README.md

Still working on verification. The README has links to Notebooks with full sky power spectra comparison.
The issue I noticed in those plots is that there are features at high ell in the synchrotron maps.

• @zonca copy masks from scratch to: /project/projectdirs/sobs/v4_sims/mbs/201904_highres_foregrounds_equatorial/{output_nside}/masks
• @NicolettaK add documentation about masks (how they were generated) in https://github.com/simonsobs/map_based_simulations/blob/master/201904_highres_foregrounds_equatorial/README.md#available-galactic-masks (we can then link to here from other releases)

I plan to switch to nside 1024 for SAT and 8192 for LAT in the next future, as that agrees with the data interface document.

Just keeping track of this for now, let's discuss in August/September

@NicolettaK should I run tophat bandpass simulations also for the SO_?1 models?

some files, even if in Equatorial, say reference frame is G in the FITS file.
Double-check for next release

In #5 we agreed to output all maps in Equatorial.
WebSky Sims instead were not rotated and assumed in Equatorial.
This creates an issue, see my note here: de26612

Same applies for the TOD simulations.
I could hack the TOD simulations script to get a consistent result but looks error prone.
I'd rather run the TOD simulations rotating everything to Equatorial, then create another extragalactic release rotated to Equatorial so that we have consistent sky in TOD Sims and mbs.

Anyone has suggestions or better solution to this?

• Update all the foregrounds models to 8192
• Use realistic bandpasses from the instrument model instead of tophat

We should have a standard CAR geometry to use for both the map-based simulations and the time-domain simulations.

We can be guided by two points:

1. The analysis interface document (see section 2.1.1) specifies that we use a CAR geometry that spans roughly -70 <= DEC <= 30 with X, Y resolution in deg/pix of CDELT1 = −1/120 , CDELT2 = 1/120 and reference pixel index CRPIX1 = 21601.0, CRPIX2 = 7561.0 at RA, DEC = (0, 0), i.e. CRVAL1 = 0.0, CRVAL2 = 0.0. Note the negative sign of CDELT1 indicating that RA decreases with increasing pixel index.
2. The quadrature rules (integration weights) for the spherical harmonic transforms. Not all versions of CAR have a corresponding quadrature rule. Without a quadrature rule spherical harmonic analysis ("map2alm") becomes inaccurate. The geometry specified in point 1 will conform to the Clenshaw–Curtis quadrature rule, but will not correspond to a known rule when the maps are downgraded using enmap.downgrade. This is a very common step in many analysis pipelines. A more suitable geometry can be obtained by shifting the y pixel centers by half a pixel. This puts the pixel centers half a pixel width away from the north and south pole, which corresponds to the Fejer1 quadrature rule. This geometry stays a valid Fejer1 geometry when downgraded. Transforming a Clenshaw–Curtis CAR map to a Fejer1 map corresponds to: CRPIX2 -= 0.5, NAXIS2 -= 1. Note that strictly speaking the second step is only necessary when the input map includes the north pole. Although I cannot find a written statement, I am quite sure we have decided that for SO we want to use the Fejer1 CAR variant. The ACT DR6 maps will likely also be also be released in the Fejer1 CAR variant.

Zach Atkins pointed out to me that the current scan-s0003 TOD sims use a different CAR geometry. They use a positive CDELT1 value and CRVAL1 = 180.0. This goes against the above specification where where RA=0 is in the middle of the map instead of at the edge. It also goes against the convention used within ACT and thus some of the plotting tools have some trouble with the positive CDELT1 value (thanks to Zach for pointing this out). Reijo has made the geometry for these maps using this script which is called using this script

I propose we update the scripts used by Reijo such that it produces output that conforms to point 1, except for the half pixel shift that makes it conform to the Fejer1 rule.

The new script would be:

import argparse
import numpy as np
from pixell import utils, enmap, curvedsky

parser = argparse.ArgumentParser()
parser.add_argument("box", type=str)
parser.add_argument("ofile")
parser.add_argument("-r", "--res",  type=float, default=0.5)
parser.add_argument("-p", "--proj", type=str, default="car")

args = parser.parse_args()
box = utils.parse_box(args.box) * utils.degree
shape, wcs = enmap.geometry(box, res=args.res * utils.arcmin, proj=args.proj)

if wcs.wcs.cdelt[0] > 0:
    raise ValueError(f'By convention the X CDELT value should be negative, got : {wcs.wcs.cdelt[0]}')
if wcs.wcs.cdelt[1] < 0:
    raise ValueError(f'By convention the Y CDELT value should be positive, got : {wcs.wcs.cdelt[1]}')

# Check if geometry conforms to Clenshaw–Curtis. This will raise error if no quad rule can be found.
# This should never happen because the above always creates CC CAR maps.
curvedsky.get_minfo(shape, wcs, quad=True)

# Convert Clenshaw–Curtis CAR map to a Feyer1 CAR map.
wcs.wcs.crpix[1] -= 0.5

# Check if input map included the north pole, if so, we need to subtract one row.
# We do this by checking if the y pixel center of the top corner of the map has been shifted above +90 deg.
if np.degrees(enmap.corners(shape, wcs, corner=False))[1,0] > 90:
    shape = (shape[0] - 1, shape[1])

# Check if geometry conforms to Feyer1. 
curvedsky.get_minfo(shape, wcs, quad=True)

enmap.write_map(args.ofile, enmap.zeros(shape, wcs, np.uint8))

The above script can be called using an updated version of the calling script:

#!/bin/bash

# The coordinates are given in format Dec_min:Dec_max,RA_min:RA_max

python3 get_area_file.py --res 15  --proj car " -63:25,180:-180" area_file_15.0.fits
python3 get_area_file.py --res 5   --proj car " -63:25,180:-180" area_file_05.0.fits
python3 get_area_file.py --res 1.5 --proj car " -63:25,180:-180" area_file_01.5.fits
python3 get_area_file.py --res 0.5 --proj car " -63:25,180:-180" area_file_00.5.fits

The last line of the script: python3 get_area_file.py --res 0.5 --proj car " -63:25,180:-180" area_file_00.5.fits will now produce a geometry with:

CDELT1 = -0.00833333
CDELT2 = 0.00833333
CRPIX1 = 21601.0
CRPIX2 = 7560.5
CRVAL1 = 0.0
CRVAL2 = 0.0
NAXIS1 = 43200
NAXIS2 = 10560

Which matches the above requirements.

The current scan-s0003 maps can be modified in place I think. We just update the WCS. They are noise maps, so a half-pixel shift won't matter. No need to rerun those.

Should we reproject or create simulations in CAR natively?

The last run of map-based sims had an inconsistency where the Galactic foregrounds were in Galactic coordinates, but the noise sims were in Equatorial coordinates: simonsobs/so_pysm_models#35 (The coordinates of the other components do not matter because the extragalactic universe is isotropic)

Let's use this issue to decide what coordinates our map-based simulations will be in (and by extension, possibly also decide for all SO maps from data and from time-domain sims). The decision here should be documented in the analysis_interface document in https://github.com/simonsobs/tod2maps_docs .

Some things to note:

1. ACT currently produces maps in Equatorial coordinates
2. Planck maps and PySM templates have been produced in Galactic coordinates
3. Reprojection from one coordinate system to the other is possible, but does not preserve every property of the map (i.e. additional characterization may be needed if a map containing data is reprojected, depending on the level of precision and scales used)
4. Instrument/atmospheric noise, ground pickup and other systematics for a ground-based experiment arguably live naturally in Equatorial coordinates.

Due to point 4, my vote is for Equatorial coordinates since operations like local Fourier space filtering are much more natural in this coordinate system. Point 3 does not affect Planck maps and PySM templates because the scales affected by reprojection are scales that are either low S/N (in the data) or are uncertain (in the simulations).

The CO component is available in so_pysm_models, see the documentation, it has been implemented by @giuspugl.

We need to check more realistic bandpasses and see if we really need to add CO to our simulations. @giuspugl would you volunteer to do this?

@simonsobs/mbs @jodunkley

Plan is to have a short summary by end of this week, so we can submit to publication committee, I'll post updates here.

Data release:

• The main product is 201906_highres_foregrounds_extragalactic_tophat
• The second is noise, unfortunately 201906_noise_no_lowell has both old hitmaps and outdated noise figures.
  I would like to have the new noise sims completed (see #21), also I would like to produce them directly in the new pixelization (variable nside, see #18). At this point I would also rerun 201906_highres_foregrounds_extragalactic_tophat with the new pixelization, it is pretty quick to do.
• Should we release the simulations for Planck (201909_highres_foregrounds_extragalactic_planck_deltabandpass) as well? (I would release them)

We should allow the mapsims API to accept a number of splits argument and provide independent noise realizations after scaling the input noise curves by the number of splits.

Some things to decide:

1. default number of splits per tube (proposal, nsplits=4)
2. whether we should save split noise sims to disk

For (2), it depends on how PS (@thibautlouis @xzackli @stevekchoi ), component separation (@jcolinhill @dpole ) and lensing (@ajvanengelen and me) plan to use the SO noise sims. For lensing, we plan to directly call the mapsims code and generate them on the fly, so we don't need them saved to disk.

For noise, can someone volunteer to make a pull request to https://github.com/simonsobs/mapsims/blob/master/mapsims/noise.py to generate noise in CAR? I'll then run the sims. You need to install pysm 3 from github.com/healpy/pysm (clone master and pip install . it)

Originally posted by @zonca in #13 (comment)

In the recent effort to inform LAT scan strategies through TOD simulations, @zatkins2 @AdriJD and others have produced models for generating map-based simulations. These are much more realistic and should replace the simplistic model. It would be good to build an interface to these models, if this has not already been done.

cc: @keskitalo @smsimon

For the next TOD simulation, we would like to use top-hat bandpasses.
All the map based simulations we have are delta-bandpass, I would like to run:

• high resolution foregrounds with fixed spectral index
• high resolution foregrounds with variable spectral index
• extragalactic

with top-hat bandpasses.

2 questions:

• what figures should I use for bandwidth?
• how do we define flat bandpasses? There is a discussion on the CMB-S4 wiki: https://cmb-s4.org/wiki/index.php/Bandpass_Convention_-_What_does_flat_mean , they actually recommend nu^-2

@AdriJD @aiolasimone I have the output of the last MBS run on Popeye.
It is 2.1T

125G    ame/
125G    ame_high/
125G    cib/
125G    cmb/
125G    cmb_unlensed/
125G    co/
125G    co_low/
125G    dust/
125G    dust_high/
125G    dust_low/
125G    freefree/
125G    ksz/
125G    radio/
125G    synchrotron/
125G    synchrotron_high/
125G    synchrotron_low/
125G    tsz/

I tried to copy to NERSC under sobs, but the quota is exceeded.
Doing some space by backing up the older simulations to tape.

Maybe we can keep all the components on Popeye and only copy combination maps (which I still need to create) to NERSC?

• 3 galactic emission maps (low, medium, high complexity)
• extragalactic
• CMB lensed
• CMB unlensed

I think it would be useful to save the 10 Gaussian realizations produced for the LAT instrument also at the SAT resolution, in order to have maps with the same foreground realization for the two instruments.
This imply the application of a simple degrade function (with hp.ud_grade to go from Nside=4096 to Nside=512) to both the dust and synchrotron maps.
@zonca could you do that?

For the LAT likelihood, we would like to test co-addition with Planck data, so we would need Planck-like simulations. Maybe we should start by generating foregrounds simulations corresponding to Planck frequency. Having Planck noise simulations would also be useful.

See the 201909_highres_foregrounds_extragalactic_planck_deltabandpass folder in the repository.

See a notebook with loglog power spectra from NSIDE 512 maps.

Please check the maps and provide feedback here or opening a dedicated issue in the repository

Just ran 10 high resolution and 100 low resolution simulations with Gaussian foregrounds, lensed CMB and noise, see https://github.com/simonsobs/map_based_simulations/tree/master/201901_gaussian_fg_lensed_cmb_realistic_noise.

They are not very expensive to run (see the benchmarks, should we run more? If so please specify how many and which components you would like to simulate.

• Implement new noise model from: https://github.com/simonsobs/V3_calc/blob/master/so_noise_lat_v3_1_CAND.py
• New hitmaps, from TOD simulations?
• Frequency channel correlations?

Hi, just wanted to point out that there seems to be some inconsistencies in the maps provided between the Planck delta-bandpass simulations and the SO tophat runs.

• Planck 545 and 857 GHz maps also have polarization
• Non-polarized foregrounds in one case have polarization equal to zero, in the other polarization is missing

Nothing crucial but maybe it is useful to keep it in mind for future runs as it requires some error-prone exception handling on the user side.

We have new noise spectra from https://github.com/simonsobs/so_noise_models

New noise spectra plots

You can run the code in the demos/ folder to reproduce the plots and play with the curves.

The main difference is a fix in the low ell contribution from the atmosphere that should be now higher at high frequencies and lower at low frequencies.

For the last noise simulations run, @thibautlouis @dpole reported we had too much low-ell noise, see #14 and I re-ran a new version of the sims with zeroed noise for ell < 30, https://github.com/simonsobs/map_based_simulations/tree/master/201906_noise_no_lowell

Is it now clear what was the cause of that?
What should we do for the next noise run which in now in preparation (#21) ?

@thibautlouis @dpole @jcolinhill @msyriac @keskitalo @mhasself

We found (Me and Davide Poletti) that the noise simulations are a bit difficult to analyse due to the very large scale modes (we are currently checking our end of the pipeline, but the noise power spectra are very very red). A good tool to debug would be to have noise simulation with no power for l<30. Could you produce that ? (the real maps will probably not have power on these scales anyway due to filtering or lack of convergence of maximum likelihood maps). Two simulations at nside=4096 per frequency channel would be enough for me. I will send more details on this later this week (but the sims would be useful in any case).

We should provide either an inverse variance map or a hit count map for each frequency as an output for the map-based simulations (both through the API and as a product saved to disk).

The mapsims origin currently says:
"It also includes simulated relative hitmaps simulated in time domain with TOAST."

Could we add more information somewhere on these time domain simulations? Is there a version number and person associated with the run? (I was specifically interested in knowing the duration of the scan simulated).

@marcelo-alvarez here are the currently available CIB frequencies:

https://portal.nersc.gov/project/cmb/so_pysm_models_data/websky/0.3/

I'm now running with bandpasses (#12), and the lower channel gets down to 21.7GHz. Can you please extend the WebSky CIB until 20 or 21 GHz?

Also, it would be nice to add a few more frequencies also between channels, so that the interpolation is not across a wide frequency range.