The SimIM Interface for Halo Catalogs

One of SimIM’s features is a uniform format for cosmological simulations. Once simulations are written into this format the simim.siminterface.SimHandler class can interface with any simulation in a homogeneous manner.

This tutorial demonstrates some basic features of the SimHandlers and SnapHandlers SimIM uses. Many of these concepts can also be applied to interacting with light cones built in SimIM.

In this example and throughout the SimIM documentation we use TNG100-3-Dark simulation. This is a simulation box with 100 Mpc sides, a dark matter particle mass of \(4\times10^8\) \(M_\odot\), and no Baryonic physics. We use this simulation because the total data volume is relatively small (~3 GB compared to ~1 TB for the full physics, full resolution TNG300-1 simulation). For most scientific applications it is advisable to utilize a higher resolution simulation. However the smaller, low-resolution TNG100-3 can be downloaded, formatted, and loaded much more quickly, making it useful for demonstration and testing purposes.

import numpy as np
import textwrap

import matplotlib.pyplot as plt
from matplotlib import rcParams
pltsty = {
    'font.size' : 18,

    # Axis appearance
    'xtick.direction' : 'in',
    'ytick.direction' : 'in',
    'xtick.top' : True,
    'ytick.right' : True,

    # Legends
    'legend.frameon' : False,

    # Lines
    'lines.dashed_pattern' : [5,3],
}
rcParams.update(pltsty)

# We've downloaded and formatted the TNG100-3-Dark simulation
# in a previous tutorial. Now we can simply initialize a
# SimHandler and hit the ground running
import simim.siminterface as sim
tng100 = sim.SimHandler('TNG100-3-Dark')

Simulation and Snapshot Metadata

A variety of metadata is stored in the SimHandler instance:

for key, value in tng100.metadata.items():
    print(f"{key}: {value}")

box_edge: 75.0
cosmo_h: 0.6774
cosmo_name: Planck2015
cosmo_omega_baryon: 0.0486
cosmo_omega_lambda: 0.6911
cosmo_omega_matter: 0.3089
groupcat_number_files: 4
name: TNG100-3-Dark
number_snaps: 100
snapshots: [( 0, 2.0046492e+01, 1.7087862e+01, 2.0046492e+01,  0.17781216,  0.22342461, 7250.3267  , 7434.8086  , 7250.3267  , 7434.8086  )
 ( 1, 1.4989173e+01, 1.3289812e+01, 1.7087862e+01,  0.22342461,  0.3186268 , 6933.908   , 7250.3267  , 6933.908   , 7250.3267  )
 ( 2, 1.1980213e+01, 1.1452663e+01, 1.3289812e+01,  0.3186268 ,  0.39193386, 6731.1284  , 6933.908   , 6731.1284  , 6933.908   )
 ( 3, 1.0975643e+01, 1.0460055e+01, 1.1452663e+01,  0.39193386,  0.44409138, 6601.795   , 6731.1284  , 6601.795   , 6731.1284  )
 ( 4, 9.9965906e+00, 9.6818905e+00, 1.0460055e+01,  0.44409138,  0.49361545, 6488.039   , 6601.795   , 6488.039   , 6601.795   )
 ( 5, 9.3887711e+00, 9.1909819e+00, 9.6818905e+00,  0.49361545,  0.5297872 , 6409.671   , 6488.039   , 6409.671   , 6488.039   )
 ( 6, 9.0023403e+00, 8.7160978e+00, 9.1909819e+00,  0.5297872 ,  0.56917566, 6328.2827  , 6409.671   , 6328.2827  , 6409.671   )
 ( 7, 8.4494762e+00, 8.2243586e+00, 8.7160978e+00,  0.56917566,  0.61536866, 6237.477   , 6328.2827  , 6237.477   , 6328.2827  )
 ( 8, 8.0121727e+00, 7.7974753e+00, 8.2243586e+00,  0.61536866,  0.66076237, 6152.563   , 6237.477   , 6152.563   , 6237.477   )
 ( 9, 7.5951071e+00, 7.4109192e+00, 7.7974753e+00,  0.66076237,  0.70688754, 6070.171   , 6152.563   , 6070.171   , 6152.563   )
 (10, 7.2362761e+00, 7.1187997e+00, 7.4109192e+00,  0.70688754,  0.74541545, 6004.053   , 6070.171   , 6004.053   , 6070.171   )
 (11, 7.0054169e+00, 6.7378883e+00, 7.1187997e+00,  0.74541545,  0.8011656 , 5912.296   , 6004.053   , 5912.296   , 6004.053   )
 (12, 6.4915977e+00, 6.2412467e+00, 6.7378883e+00,  0.8011656 ,  0.8850064 , 5781.999   , 5912.296   , 5781.999   , 5912.296   )
 (13, 6.0107574e+00, 5.9274750e+00, 6.2412467e+00,  0.8850064 ,  0.94579893, 5692.5874  , 5781.999   , 5692.5874  , 5781.999   )
 (14, 5.8466139e+00, 5.6835189e+00, 5.9274750e+00,  0.94579893,  0.9980251 , 5618.7886  , 5692.5874  , 5618.7886  , 5692.5874  )
 (15, 5.5297656e+00, 5.3742218e+00, 5.6835189e+00,  0.9980251 ,  1.0714717 , 5519.2417  , 5618.7886  , 5519.2417  , 5618.7886  )
 (16, 5.2275810e+00, 5.1090288e+00, 5.3742218e+00,  1.0714717 ,  1.1418936 , 5427.986   , 5519.2417  , 5427.986   , 5519.2417  )
 (17, 4.9959335e+00, 4.8243771e+00, 5.1090288e+00,  1.1418936 ,  1.2264551 , 5323.2437  , 5427.986   , 5323.2437  , 5427.986   )
 (18, 4.6645179e+00, 4.5431485e+00, 4.8243771e+00,  1.2264551 ,  1.3207208 , 5212.022   , 5323.2437  , 5212.022   , 5323.2437  )
 (19, 4.4280338e+00, 4.2987480e+00, 4.5431485e+00,  1.3207208 ,  1.4128404 , 5108.349   , 5212.022   , 5108.349   , 5212.022   )
 (20, 4.1768351e+00, 4.0906568e+00, 4.2987480e+00,  1.4128404 ,  1.5000024 , 5014.342   , 5108.349   , 5014.342   , 5108.349   )
 (21, 4.0079451e+00, 3.8526685e+00, 4.0906568e+00,  1.5000024 ,  1.6111543 , 4899.6235  , 5014.342   , 4899.6235  , 5014.342   )
 (22, 3.7087743e+00, 3.5966349e+00, 3.8526685e+00,  1.6111543 ,  1.7468175 , 4766.5815  , 4899.6235  , 4766.5815  , 4899.6235  )
 (23, 3.4908614e+00, 3.3839178e+00, 3.5966349e+00,  1.7468175 ,  1.8745713 , 4647.4937  , 4766.5815  , 4647.4937  , 4766.5815  )
 (24, 3.2830331e+00, 3.1399896e+00, 3.3839178e+00,  1.8745713 ,  2.041211  , 4500.088   , 4647.4937  , 4500.088   , 4647.4937  )
 (25, 3.0081310e+00, 2.9509807e+00, 3.1399896e+00,  2.041211  ,  2.1879046 , 4376.889   , 4500.088   , 4376.889   , 4500.088   )
 (26, 2.8957851e+00, 2.8123467e+00, 2.9509807e+00,  2.1879046 ,  2.3069239 , 4280.962   , 4376.889   , 4280.962   , 4376.889   )
 (27, 2.7331426e+00, 2.6531942e+00, 2.8123467e+00,  2.3069239 ,  2.4572814 , 4164.438   , 4280.962   , 4164.438   , 4280.962   )
 (28, 2.5772903e+00, 2.5092282e+00, 2.6531942e+00,  2.4572814 ,  2.6077237 , 4052.5771  , 4164.438   , 4052.5771  , 4164.438   )
 (29, 2.4442258e+00, 2.3787005e+00, 2.5092282e+00,  2.6077237 ,  2.7576826 , 3945.3455  , 4052.5771  , 3945.3455  , 4052.5771  )
 (30, 2.3161108e+00, 2.2609375e+00, 2.3787005e+00,  2.7576826 ,  2.905529  , 3843.4312  , 3945.3455  , 3843.4312  , 3945.3455  )
 (31, 2.2079256e+00, 2.1545560e+00, 2.2609375e+00,  2.905529  ,  3.050601  , 3746.8022  , 3843.4312  , 3746.8022  , 3843.4312  )
 (32, 2.1032696e+00, 2.0516453e+00, 2.1545560e+00,  3.050601  ,  3.202571  , 3648.8809  , 3746.8022  , 3648.8809  , 3746.8022  )
 (33, 2.0020282e+00, 1.9520924e+00, 2.0516453e+00,  3.202571  ,  3.361715  , 3549.682   , 3648.8809  , 3549.682   , 3648.8809  )
 (34, 1.9040896e+00, 1.8627039e+00, 1.9520924e+00,  3.361715  ,  3.5159051 , 3456.5928  , 3549.682   , 3456.5928  , 3549.682   )
 (35, 1.8226893e+00, 1.7824667e+00, 1.8627039e+00,  3.5159051 ,  3.6643636 , 3369.5762  , 3456.5928  , 3369.5762  , 3456.5928  )
 (36, 1.7435706e+00, 1.7044784e+00, 1.7824667e+00,  3.6643636 ,  3.8186705 , 3281.6665  , 3369.5762  , 3281.6665  , 3369.5762  )
 (37, 1.6666696e+00, 1.6350199e+00, 1.7044784e+00,  3.8186705 ,  3.9651728 , 3200.444   , 3281.6665  , 3200.444   , 3281.6665  )
 (38, 1.6042346e+00, 1.5671337e+00, 1.6350199e+00,  3.9651728 ,  4.1173496 , 3118.2454  , 3200.444   , 3118.2454  , 3200.444   )
 (39, 1.5312390e+00, 1.5132287e+00, 1.5671337e+00,  4.1173496 ,  4.2450213 , 3050.895   , 3118.2454  , 3050.895   , 3118.2454  )
 (40, 1.4955121e+00, 1.4540279e+00, 1.5132287e+00,  4.2450213 ,  4.3927426 , 2974.7043  , 3050.895   , 2974.7043  , 3050.895   )
 (41, 1.4140983e+00, 1.3854544e+00, 1.4540279e+00,  4.3927426 ,  4.574445  , 2883.402   , 2974.7043  , 2883.402   , 2974.7043  )
 (42, 1.3575766e+00, 1.3296057e+00, 1.3854544e+00,  4.574445  ,  4.7315087 , 2806.5056  , 2883.402   , 2806.5056  , 2883.402   )
 (43, 1.3023784e+00, 1.2750646e+00, 1.3296057e+00,  4.7315087 ,  4.8933887 , 2729.1067  , 2806.5056  , 2729.1067  , 2806.5056  )
 (44, 1.2484726e+00, 1.2271405e+00, 1.2750646e+00,  4.8933887 ,  5.04305   , 2659.1409  , 2729.1067  , 2659.1409  , 2729.1067  )
 (45, 1.2062581e+00, 1.1800886e+00, 1.2271405e+00,  5.04305   ,  5.1972003 , 2588.5964  , 2659.1409  , 2588.5964  , 2659.1409  )
 (46, 1.1546028e+00, 1.1341639e+00, 1.1800886e+00,  5.1972003 ,  5.3549986 , 2517.906   , 2588.5964  , 2517.906   , 2588.5964  )
 (47, 1.1141505e+00, 1.0940968e+00, 1.1341639e+00,  5.3549986 ,  5.4989595 , 2454.6987  , 2517.906   , 2454.6987  , 2517.906   )
 (48, 1.0744579e+00, 1.0547818e+00, 1.0940968e+00,  5.4989595 ,  5.6462502 , 2391.2437  , 2454.6987  , 2391.2437  , 2454.6987  )
 (49, 1.0355104e+00, 1.0162051e+00, 1.0547818e+00,  5.6462502 ,  5.7969084 , 2327.5562  , 2391.2437  , 2327.5562  , 2391.2437  )
 (50, 9.9729425e-01, 9.7361338e-01, 1.0162051e+00,  5.7969084 ,  5.9707108 , 2255.5513  , 2327.5562  , 2255.5513  , 2327.5562  )
 (51, 9.5053136e-01, 9.3666106e-01, 9.7361338e-01,  5.9707108 ,  6.1282086 , 2191.6008  , 2255.5513  , 2191.6008  , 2255.5513  )
 (52, 9.2300081e-01, 9.0476638e-01, 9.3666106e-01,  6.1282086 ,  6.2694335 , 2135.2666  , 2191.6008  , 2135.2666  , 2191.6008  )
 (53, 8.8689691e-01, 8.6900622e-01, 9.0476638e-01,  6.2694335 ,  6.433894  , 2070.82    , 2135.2666  , 2070.82    , 2135.2666  )
 (54, 8.5147089e-01, 8.3391738e-01, 8.6900622e-01,  6.433894  ,  6.601875  , 2006.23    , 2070.82    , 2006.23    , 2070.82    )
 (55, 8.1671000e-01, 8.0379409e-01, 8.3391738e-01,  6.601875  ,  6.7515683 , 1949.6843  , 2006.23    , 1949.6843  , 2006.23    )
 (56, 7.9106826e-01, 7.7408981e-01, 8.0379409e-01,  6.7515683 ,  6.9043756 , 1892.9116  , 1949.6843  , 1892.9116  , 1949.6843  )
 (57, 7.5744140e-01, 7.4494815e-01, 7.7408981e-01,  6.9043756 ,  7.0595274 , 1836.2161  , 1892.9116  , 1836.2161  , 1892.9116  )
 (58, 7.3263615e-01, 7.1621406e-01, 7.4494815e-01,  7.0595274 ,  7.2178154 , 1779.3257  , 1836.2161  , 1779.3257  , 1836.2161  )
 (59, 7.0010638e-01, 6.8802208e-01, 7.1621406e-01,  7.2178154 ,  7.378461  , 1722.5376  , 1779.3257  , 1722.5376  , 1779.3257  )
 (60, 6.7611039e-01, 6.6022646e-01, 6.8802208e-01,  7.378461  ,  7.5422544 , 1665.5889  , 1722.5376  , 1665.5889  , 1722.5376  )
 (61, 6.4464182e-01, 6.3295317e-01, 6.6022646e-01,  7.5422544 ,  7.708411  , 1608.7687  , 1665.5889  , 1608.7687  , 1665.5889  )
 (62, 6.2142873e-01, 6.0990566e-01, 6.3295317e-01,  7.708411  ,  7.8531947 , 1560.0133  , 1608.7687  , 1560.0133  , 1608.7687  )
 (63, 5.9854329e-01, 5.8718342e-01, 6.0990566e-01,  7.8531947 ,  8.000003  , 1511.2742  , 1560.0133  , 1511.2742  , 1560.0133  )
 (64, 5.7598084e-01, 5.6105018e-01, 5.8718342e-01,  8.000003  ,  8.174028  , 1454.3821  , 1511.2742  , 1454.3821  , 1511.2742  )
 (65, 5.4639220e-01, 5.3540426e-01, 5.6105018e-01,  8.174028  ,  8.350409  , 1397.6691  , 1454.3821  , 1397.6691  , 1454.3821  )
 (66, 5.2456582e-01, 5.1373357e-01, 5.3540426e-01,  8.350409  ,  8.503937  , 1349.0573  , 1397.6691  , 1349.0573  , 1397.6691  )
 (67, 5.0304753e-01, 4.9236873e-01, 5.1373357e-01,  8.503937  ,  8.65947   , 1300.5061  , 1349.0573  , 1300.5061  , 1349.0573  )
 (68, 4.8183295e-01, 4.7130546e-01, 4.9236873e-01,  8.65947   ,  8.816997  , 1252.026   , 1300.5061  , 1252.026   , 1300.5061  )
 (69, 4.6091780e-01, 4.5053947e-01, 4.7130546e-01,  8.816997  ,  8.976511  , 1203.6277  , 1252.026   , 1203.6277  , 1252.026   )
 (70, 4.4029784e-01, 4.3006656e-01, 4.5053947e-01,  8.976511  ,  9.138     , 1155.3215  , 1203.6277  , 1155.3215  , 1203.6277  )
 (71, 4.1996893e-01, 4.0988263e-01, 4.3006656e-01,  9.138     ,  9.301452  , 1107.1182  , 1155.3215  , 1107.1182  , 1155.3215  )
 (72, 3.9992696e-01, 3.8998356e-01, 4.0988263e-01,  9.301452  ,  9.466854  , 1059.0281  , 1107.1182  , 1059.0281  , 1107.1182  )
 (73, 3.8016787e-01, 3.7036535e-01, 3.8998356e-01,  9.466854  ,  9.634192  , 1011.0619  , 1059.0281  , 1011.0619  , 1059.0281  )
 (74, 3.6068764e-01, 3.5424355e-01, 3.7036535e-01,  9.634192  ,  9.774977  ,  971.22943 , 1011.0619  ,  971.22943 , 1011.0619  )
 (75, 3.4785384e-01, 3.3828175e-01, 3.5424355e-01,  9.774977  ,  9.917349  ,  931.4219  ,  971.22943 ,  931.4219  ,  971.22943 )
 (76, 3.2882974e-01, 3.1939328e-01, 3.3828175e-01,  9.917349  , 10.08977   ,  883.837   ,  931.4219  ,  883.837   ,  931.4219  )
 (77, 3.1007412e-01, 3.0387032e-01, 3.1939328e-01, 10.08977   , 10.234756  ,  844.3411  ,  883.837   ,  844.3411  ,  883.837   )
 (78, 2.9771769e-01, 2.8543562e-01, 3.0387032e-01, 10.234756  , 10.410904  ,  796.97845 ,  844.3411  ,  796.97845 ,  844.3411  )
 (79, 2.7335334e-01, 2.6732391e-01, 2.8543562e-01, 10.410904  , 10.588276  ,  749.9597  ,  796.97845 ,  749.9597  ,  796.97845 )
 (80, 2.6134327e-01, 2.5238788e-01, 2.6732391e-01, 10.588276  , 10.737853  ,  710.822   ,  749.9597  ,  710.822   ,  749.9597  )
 (81, 2.4354018e-01, 2.3471172e-01, 2.5238788e-01, 10.737853  , 10.9188385 ,  664.0792  ,  710.822   ,  664.0792  ,  710.822   )
 (82, 2.2598839e-01, 2.2018380e-01, 2.3471172e-01, 10.9188385 , 11.070893  ,  625.3165  ,  664.0792  ,  625.3165  ,  664.0792  )
 (83, 2.1442504e-01, 2.0580406e-01, 2.2018380e-01, 11.070893  , 11.224408  ,  586.64233 ,  625.3165  ,  586.64233 ,  625.3165  )
 (84, 1.9728418e-01, 1.8878537e-01, 2.0580406e-01, 11.224408  , 11.410067  ,  540.4763  ,  586.64233 ,  540.4763  ,  586.64233 )
 (85, 1.8038526e-01, 1.7479715e-01, 1.8878537e-01, 11.410067  , 11.565969  ,  502.21097 ,  540.4763  ,  502.21097 ,  540.4763  )
 (86, 1.6925204e-01, 1.6095297e-01, 1.7479715e-01, 11.565969  , 11.723276  ,  464.05576 ,  502.21097 ,  464.05576 ,  502.21097 )
 (87, 1.5274876e-01, 1.4729182e-01, 1.6095297e-01, 11.723276  , 11.881506  ,  426.12854 ,  464.05576 ,  426.12854 ,  464.05576 )
 (88, 1.4187621e-01, 1.3377218e-01, 1.4729182e-01, 11.881506  , 12.041101  ,  388.32434 ,  426.12854 ,  388.32434 ,  426.12854 )
 (89, 1.2575933e-01, 1.1777012e-01, 1.3377218e-01, 12.041101  , 12.233959  ,  343.23273 ,  388.32434 ,  343.23273 ,  388.32434 )
 (90, 1.0986994e-01, 1.0461648e-01, 1.1777012e-01, 12.233959  , 12.39578   ,  305.8876  ,  343.23273 ,  305.8876  ,  343.23273 )
 (91, 9.9401802e-02, 9.1600336e-02, 1.0461648e-01, 12.39578   , 12.558904  ,  268.68506 ,  305.8876  ,  268.68506 ,  305.8876  )
 (92, 8.3884433e-02, 7.8754269e-02, 9.1600336e-02, 12.558904  , 12.722881  ,  231.72787 ,  268.68506 ,  231.72787 ,  268.68506 )
 (93, 7.3661387e-02, 6.6043235e-02, 7.8754269e-02, 12.722881  , 12.888123  ,  194.92479 ,  231.72787 ,  194.92479 ,  231.72787 )
 (94, 5.8507323e-02, 5.3497560e-02, 6.6043235e-02, 12.888123  , 13.05419   ,  158.37294 ,  194.92479 ,  158.37294 ,  194.92479 )
 (95, 4.8523631e-02, 4.1084476e-02, 5.3497560e-02, 13.05419   , 13.221478  ,  121.98607 ,  158.37294 ,  121.98607 ,  158.37294 )
 (96, 3.3724371e-02, 2.8832177e-02, 4.1084476e-02, 13.221478  , 13.389566  ,   85.855835,  121.98607 ,   85.855835,  121.98607 )
 (97, 2.3974428e-02, 1.6710050e-02, 2.8832177e-02, 13.389566  , 13.558832  ,   49.900814,   85.855835,   49.900814,   85.855835)
 (98, 9.5216669e-03, 4.7442745e-03, 1.6710050e-02, 13.558832  , 13.728869  ,   14.207323,   49.900814,   14.207323,   49.900814)
 (99, 2.2204460e-16, 0.0000000e+00, 4.7442745e-03, 13.728869  , 13.797113  ,    0.      ,   14.207323,    0.      ,   14.207323)]

Additionally, we can review the fields in the halo catalogs:

fields = [f for f in tng100.extract_snap_keys()]
print(textwrap.fill(f"The halo catalog fields for this simulation are: " + ', '.join(fields),subsequent_indent='    '))

The halo catalog fields for this simulation are: cm_x, cm_y, cm_z,
    mass, mass_maxr, pos_x, pos_y, pos_z, r_hm, r_vmax, spin_x,
    spin_y, spin_z, v_x, v_y, v_z, vdisp, vmax

Simulation Snapshots and SnapHanlder

Each simulation is composed of a series of snapshots which contain the halo catalog frozen at a given point in time (or equivalently at a given redshift). Because the simulations are large, by default the SimHandler doesn’t load the data for the snapshots - this is done to reduce memory usage since many applications don’t require using data from all snapshots. To interact with data from an individual snapshot, we need to retrieve it.

First let’s determine the snapshot associated with our redshift of interest. We’ll search for the snapshot closest to z=1.0:

z = 1.0
snap = tng100.z_to_snap(z)
print(textwrap.fill(f"The snapshot closest to redshift {z :.1f} has index {snap}",subsequent_indent='    '))

The snapshot closest to redshift 1.0 has index 50

Now we can load this snapshot into memory. Each snapshot comes with its own SnapHandler object which is used to store both metadata and the halo catalog for that snapshot.

tng100_z1 = tng100.get_snap(snap)

### Equivalently we could have used
# tng100_z1 = tng100.get_snap_from_z(z)

SnapHandlers still don’t automatically load all of the data for a snapshot, again trying to save RAM because we’re working with large files. Instead they have a list of properties associated with the halo catalog, which they will load as needed.

We can see the properties associated with a snapshot:

fields = [f for f in tng100_z1.extract_keys()]
print(textwrap.fill("The halo catalog fields for this snapshot are: "+', '.join([f for f in fields]),subsequent_indent='    '))

The halo catalog fields for this snapshot are: cm_x, cm_y, cm_z, mass,
    mass_maxr, pos_x, pos_y, pos_z, r_hm, r_vmax, spin_x, spin_y,
    spin_z, v_x, v_y, v_z, vdisp, vmax

We can also specifically see what properties are loaded into memory, what properties are saved to the disk, and what additonal properties have been added since loading the snapshot, but aren’t saved to the disk:

fields = [f for f in tng100_z1.extract_keys(set='loaded')]
print(textwrap.fill("The loaded fields for this snapshot are: "+', '.join([f for f in fields]),subsequent_indent='    '))

fields = [f for f in tng100_z1.extract_keys(set='saved')]
print(textwrap.fill("The saved fields for this snapshot are: "+', '.join([f for f in fields]),subsequent_indent='    '))

fields = [f for f in tng100_z1.extract_keys(set='generated')]
print(textwrap.fill("The new fields for this snapshot are: "+', '.join([f for f in fields]),subsequent_indent='    '))

The loaded fields for this snapshot are:
The saved fields for this snapshot are: cm_x, cm_y, cm_z, mass,
    mass_maxr, pos_x, pos_y, pos_z, r_hm, r_vmax, spin_x, spin_y,
    spin_z, v_x, v_y, v_z, vdisp, vmax
The new fields for this snapshot are:

We can load fields into memory explicity:

tng100_z1.load_property('pos_x','pos_y','pos_z')

Now printing the loaded fields will list pos_x, pos_y, and pos_z. These are still on the disk and therefore still appear in the saved fields list as well:

fields = [f for f in tng100_z1.extract_keys(set='loaded')]
print(textwrap.fill("The loaded fields for this snapshot are: "+', '.join([f for f in fields]),subsequent_indent='    '))

fields = [f for f in tng100_z1.extract_keys(set='saved')]
print(textwrap.fill("The saved fields for this snapshot are: "+', '.join([f for f in fields]),subsequent_indent='    '))

fields = [f for f in tng100_z1.extract_keys(set='generated')]
print(textwrap.fill("The new fields for this snapshot are: "+', '.join([f for f in fields]),subsequent_indent='    '))

The loaded fields for this snapshot are: pos_x, pos_y, pos_z
The saved fields for this snapshot are: cm_x, cm_y, cm_z, mass,
    mass_maxr, pos_x, pos_y, pos_z, r_hm, r_vmax, spin_x, spin_y,
    spin_z, v_x, v_y, v_z, vdisp, vmax
The new fields for this snapshot are:

We can also make simple plots of these properties. For instance, we can plot the x and y positions of all halos:

tng100_z1.plot('pos_x','pos_y',
               axkws={'xlabel':'x-position [Mpc/h]','ylabel':'y-position [Mpc/h]','aspect':'equal'},
               plotkws={'marker':',','color':'k'})

This plots all halos x- and y-positions. Because the box is thick in the z-direction as well, it can result in a very dense looking plot. We can fix this by restricting the set of halos we use:

tng100_z1.set_property_range('pos_z',0,10) # Restrict catalog to only halos with pos_z<10 Mpc/h
tng100_z1.plot('pos_x','pos_y',
               axkws={'xlabel':'x-position [Mpc/h]','ylabel':'y-position [Mpc/h]','aspect':'equal'},
               plotkws={'marker':',','color':'k'})

To undo this restriction, simpy call set_property_range with no arguments. Note that by default set_property_range will overwrite the last set of restrictions set, so calling set_property_range('pos_x',0,10) followed by set_property_range('pos_y',0,10) would result in the halo catalog being restricted only along the y-direction and not the x direction:

tng100_z1.set_property_range('pos_x',0,10) # Restrict catalog to only halos with pos_x<10 Mpc/h
tng100_z1.plot('pos_x','pos_y',
                 axkws={'xlabel':'x-position [Mpc/h]','ylabel':'y-position [Mpc/h]','aspect':'equal','xlim':[0,75],'ylim':[0,75]},
                 plotkws={'marker':',','color':'k'})

tng100_z1.set_property_range('pos_y',0,10) # Restrict catalog to only halos with pos_y<10 Mpc/h
tng100_z1.plot('pos_x','pos_y',
                 axkws={'xlabel':'x-position [Mpc/h]','ylabel':'y-position [Mpc/h]','aspect':'equal','xlim':[0,75],'ylim':[0,75]},
                 plotkws={'marker':',','color':'k'})

To apply multiple cuts, set the reset parameter to False:

tng100_z1.set_property_range('pos_x',0,10) # Restrict catalog to only halos with pos_x<10 Mpc/h
tng100_z1.set_property_range('pos_y',0,10,reset=False) # Restrict catalog to only halos with pos_y<10 Mpc/h
tng100_z1.plot('pos_x','pos_y',
                 axkws={'xlabel':'x-position [Mpc/h]','ylabel':'y-position [Mpc/h]','aspect':'equal','xlim':[0,75],'ylim':[0,75]},
                 plotkws={'marker':',','color':'k'})

# Undo these cuts.
tng100_z1.set_property_range()

We can also make histograms. Here we’ll plot the distribution of halo masses in our catalog. We can see that below a mass of about \(10^{10}\) \(M_\odot\) the low mass resolution of the TNG100-3 simulation results in a sharp fall off of the number of halos.

tng100_z1.hist('mass',
                 axkws={'xlabel':'Halo Mass [M$_\odot$/h]','xscale':'log','yscale':'log'},
                 plotkws={'bins':np.logspace(9,15,61)})

Note that we never loaded the ‘mass’ property. The SnapHandler will load (and then unload) the property automatically when it is needed to execute a method call.

Direct Access to Halo/Galaxy Properties

The .plot and .hist methods are meant primarily for exploratory data visualization, and not detailed analysis. It’s always possible to extract the underlying data and make more sophisticated plots. The following code is an equivalent way to make the halo mass histogram shown above:

fig, ax = plt.subplots()
ax.set(xlabel='Halo Mass [M$_\odot$/h]',xscale='log',yscale='log')

masses = tng100_z1.return_property('mass')
ax.hist(masses, bins=np.logspace(9,15,61))

plt.show()

Any property associated with a SnapHanlder instance can be retrieved using the .return_property method.

Adding Properties to Simulations

Because TNG100-3-Dark is a dark matter only simulation, it doesn’t include important properties needed for modeling galaxy observations, such as star formation rates. We will need to add these ourselves. (Note: we could also use a different simulation which includes a model for star formation, such as the full-hydrodynamics versions of TNG)

SimIM handles adding new properties to snapshots (and other objects) using the galprops module. The key element of this module is the Prop class, which contains some function for generating a field in the halo catalog, given some existing halo catalog fields, along with basic instructions for how to execute the function.

The galprops module has a number of pre-defined properties, and we will work with a few of those first before seeing how to construct new ones.

Let’s start by assigning star formation rates to each halo. To do this we will use the prescription of Behroozi et al. 2013, which is built into the galprops module:

# Load the Behroozi et al. halo mass to SFR prescription
from simim.galprops import prop_behroozi_sfr

Adding this property to the halo catalog of our snapshot is very simple:

tng100_z1.make_property(prop_behroozi_sfr)

‘sfr’ will now appear in the list of properties associated with the snapshot. Note that it has been added to the list of loaded fields and the list of newly generated fields for the snapshot, but not to the list of fields saved to the disk. If we want to, we can save our new property using the command SnapHandler.write_property('sfr'). However, doing so would write the property for only a single snapshot, rather than the whole simulation, which is not generally desirable. We will revisit writing properties to whole simulations in a moment.

fields = [f for f in tng100_z1.extract_keys(set='loaded')]
print(textwrap.fill("The loaded fields for this snapshot are: "+', '.join([f for f in fields]),subsequent_indent='    '))

fields = [f for f in tng100_z1.extract_keys(set='saved')]
print(textwrap.fill("The saved fields for this snapshot are: "+', '.join([f for f in fields]),subsequent_indent='    '))

fields = [f for f in tng100_z1.extract_keys(set='generated')]
print(textwrap.fill("The new fields for this snapshot are: "+', '.join([f for f in fields]),subsequent_indent='    '))

The loaded fields for this snapshot are: pos_x, pos_y, pos_z, sfr
The saved fields for this snapshot are: cm_x, cm_y, cm_z, mass,
    mass_maxr, pos_x, pos_y, pos_z, r_hm, r_vmax, spin_x, spin_y,
    spin_z, v_x, v_y, v_z, vdisp, vmax
The new fields for this snapshot are: sfr

We can also add line luminosities. We can often model the luminosity of a given spectral line to a halo using a prescription of the form

\[\log L_\mathrm{line} = a\log\mathrm{SFR}+b+\epsilon(\sigma)\]

where \(a\) and \(b\) are an emperically determined coeficients of a power law fit and \(\epsilon\) is a random scatter drawn from an (assumed) normal distribution of width \(\sigma\).

Fits for the correlation between SFR and the 158µm [CII] line have been presented multiple places in the literature. For example DeLooze et al. 2014, provides fits for multiple far infrared lines and multiple galaxy populations. This prescription is also built in to galprops, and we can add it to our halos as follows:

from simim.galprops import prop_delooze_cii
tng100_z1.make_property(prop_delooze_cii,rename='LCII')

By default, prop_delooze_cii uses the fit for the whole galaxy sample presented in DeLooze et al. However, we can also use the fits from subsamples, for example the fit for only starburst galaxies. To do this, we need to specify a few parameters when calling the make_property method. The other_kws parameter is a dictionary of keyword arguments that are passed directly to the function call that generates the property - in this case, set='starbursts' will be passed to the function underlying prop_delooze_cii. Because we already assigned the property ‘LCII’ (the default property name for prop_delooze_cii), we will also specify a new property name using the rename parameter to avoid overwriting the previous result.

tng100_z1.make_property(prop_delooze_cii,other_kws={'set':'starbursts'},rename='LCII_sb')

fields = [f for f in tng100_z1.extract_keys(set='loaded')]
print(textwrap.fill("The loaded fields for this snapshot are: "+', '.join([f for f in fields]),subsequent_indent='    '))

fields = [f for f in tng100_z1.extract_keys(set='saved')]
print(textwrap.fill("The saved fields for this snapshot are: "+', '.join([f for f in fields]),subsequent_indent='    '))

fields = [f for f in tng100_z1.extract_keys(set='generated')]
print(textwrap.fill("The new fields for this snapshot are: "+', '.join([f for f in fields]),subsequent_indent='    '))

The loaded fields for this snapshot are: pos_x, pos_y, pos_z, sfr,
    LCII, LCII_sb
The saved fields for this snapshot are: cm_x, cm_y, cm_z, mass,
    mass_maxr, pos_x, pos_y, pos_z, r_hm, r_vmax, spin_x, spin_y,
    spin_z, v_x, v_y, v_z, vdisp, vmax
The new fields for this snapshot are: sfr, LCII, LCII_sb

Now we can compare the luminosities generated by the two prescriptions:

tng100_z1.plot('LCII','LCII_sb',
               axkws={'xlabel':'L$_\mathregular{CII}$ [L$_\odot$] - Prescription 1',
                      'ylabel':'L$_\mathregular{CII}$ [L$_\odot$] - Prescription 2',
                      'xscale':'log','yscale':'log'},
               plotkws={'marker':',','color':'k'})

Note that there is quite a bit of scatter between the two prescriptions. Both prescriptions add a random scatter to the luminosity of every halo, but as it stands we have not done anything to make sure that the same halos are scattered in the same way (i.e. fixed the seed for our random number generator). Let’s fix this.

rng = np.random.default_rng(85712)
tng100_z1.make_property(prop_delooze_cii,rename='LCII',other_kws={'rng':rng},overwrite=True)
rng = np.random.default_rng(85712) # Reseed the RNG to get the same numbers for the second application
tng100_z1.make_property(prop_delooze_cii,other_kws={'rng':rng,'set':'starbursts'},rename='LCII_sb',overwrite=True)

tng100_z1.plot('LCII','LCII_sb',
               axkws={'xlabel':'L$_\mathregular{CII}$ [L$_\odot$] - Prescription 1',
                      'ylabel':'L$_\mathregular{CII}$ [L$_\odot$] - Prescription 2',
                      'xscale':'log','yscale':'log'},
               plotkws={'marker':',','color':'k'})

/Users/keenan/Dropbox/4_research/2.1_simim_release/simim/_handlers.py:269: UserWarning: Property LCII already exists, overwriting
  warnings.warn("Property {} already exists, overwriting".format(name))
/Users/keenan/Dropbox/4_research/2.1_simim_release/simim/_handlers.py:269: UserWarning: Property LCII_sb already exists, overwriting
  warnings.warn("Property {} already exists, overwriting".format(name))

This clearly reduces the scatter between the two prescriptions. Some dispersion remains, because the amplitude of the random scatter added by the two prescriptions differs, but this should be a better starting point for making an apples-to-apples comparisons between results found using the different models.

Now let’s implement our own Prop instance for a different [CII] prescription. (Note: the prescription described here is already built into SimIM as simim.galprops.prop_schaerer_cii, we recreate it here for illustrative purposes)

We’ll use SFR-[CII] the scaling relation fit by Schaerer et al. 2020 using the z~4-6 galaxies from the ALPINE survey. We begin by defining a function to assign a CII luminosity from a given SFR using the Schaerer fit plus a lognormal scatter of 0.3 dex:

def lcii_schaerer(sfr,scatter_dex=0.3,rng=np.random.default_rng()):

    # Determine mean luminosity from Schaerer fit:
    a = 1.17
    b = 6.61
    lcii = 10**(a*np.log10(np.array(sfr))+b)

    # Add a lognormal scatter
    lcii = lcii / 10**(scatter_dex**2 * np.log(10)/2) # Adjustment required to preserve linear mean
    lcii = lcii * 10**(scatter_dex*rng.normal(loc=0,scale=1,size=lcii.shape))

    return lcii

Now we can wrap this function in a Prop class that supplies some additional information needed by the SnapHandler to evaluate it. In particular, we will give a default property name, a list of snapshot properties to be passed in the function call (in this case just ‘sfr’), the units of the resulting property (\(L_\odot\)), the dependence of the property on the Hubble parameter, and whether the arguments fed to the function call (SFRs) are expected to have units with little h (they are not) and the values returned by the function call have units with little h (again they do not). All of this little h accounting is because most simulations work with at least some properties in little h units, and keeping these can be useful for translating between different cosmologies. But in most cases similar to correlating SFRs and line luminosities we do not need the little h units, by default the prop class will assume no h dependence and that all arguments and returns are in units where the h dependence has already been evaluated.

from simim.galprops.galprops import Prop

prop_schaerer_cii = Prop(prop_name='LCII',
                         prop_function=lcii_schaerer,
                         kwargs=['sfr'],
                         units='Lsun',
                         h_dependence=-2,
                         give_args_in_h_units=False,
                         function_return_in_h_units=False)

From here we can use the new Prop instance in exactly the same manner as we did previously:

rng = np.random.default_rng(85712) # Reseed the RNG again
tng100_z1.make_property(prop_schaerer_cii,other_kws={'rng':rng},rename='LCII_z5')

tng100_z1.plot('LCII','LCII_z5',
               axkws={'xlabel':'L$_\mathregular{CII}$ [L$_\odot$] - Prescription 1',
                      'ylabel':'L$_\mathregular{CII}$ [L$_\odot$] - Prescription 3',
                      'xscale':'log','yscale':'log'},
               plotkws={'marker':',','color':'k'})

Interacting with a Simulation Across All Snapshots

Now that we’ve seen the basics of how adding properties to snapshots works, we can go about adding properties to the entire simulation. This is valuable for assessing the evolution of properties with time, and is also necessary if we want to add new properties to the simulations that can be inherited by light cones (survey volumes) we generate.

The SimHandler has methods for handling this - essentially just wrappers that iterate over each snapshot. We will add SFRs, and our three \(L_\mathrm{CII}\) prescriptions to the whole simulation. In the following, we will use our SimHandler instance instead of the SnapHandler instance for only z=1. Otherwise the method calls are quite similar to those used earlier, but the additional argument write=True will save the new properties we generate to the disk, making them useable again in the future.

rng1 = np.random.default_rng(58290)
tng100.make_property(prop_behroozi_sfr,other_kws={'rng':rng1},write=True,overwrite=True)

rng2 = np.random.default_rng(85712)
tng100.make_property(prop_delooze_cii,other_kws={'rng':rng2},rename='LCII',write=True,overwrite=True)
rng2 = np.random.default_rng(85712)
tng100.make_property(prop_delooze_cii,other_kws={'set':'starbursts','rng':rng2},rename='LCII_sb',write=True,overwrite=True)
rng2 = np.random.default_rng(85712)
tng100.make_property(prop_schaerer_cii,other_kws={'rng':rng2},rename='LCII_z5',write=True,overwrite=True)

Assigning props for Snapshot 99.  

Now that we’ve constructed our properties, let’s use assess their evolution over cosmic history. We’ll start with the cosmic star formation rate density:

# Function to compute the SFRD
def sfrdensity(sfr,box_edge,h):
    return np.sum(sfr)/(box_edge/h)**3

# 'box_edge', and 'h' are two special keywords that Handler objects recognize
# corresponding to the side length of the box and little h
sfrdensity_val, redshift = tng100.snap_stat(stat_function=sfrdensity,
                                            kwargs=['sfr','box_edge','h'])

Collecting sources from Snapshot 99.  

# We'll compare this result to the fit from Madau & Dickinson 2014
redshift_madau = np.linspace(0,10,1000)
sfrdensity_madau = 0.015 * (1+redshift_madau)**2.7 / (1+((1+redshift_madau)/2.9)**5.6)

# Plot the stuff
fig,ax = plt.subplots()
ax.set(xlabel='Redshift',ylabel=r'SFRD [M$_\odot$ yr$^{-1}$ Mpc$^{-3}$]',yscale='log',xlim=[0,5],ylim=[.01,.2])
ax.plot(redshift,sfrdensity_val,color='C0',label='Simulated')
ax.plot(redshift_madau,sfrdensity_madau,color='k',ls='-',lw=4,label='MD14')
ax.legend(fontsize='small')
plt.show()

We can see that our model’s SFRD deviates considerably from the fit to data from Madau and Dickinson 2014. These deviations are primarily at z>3 where uncertainties in the SFRD were large when Madau and Dickinson 2014 and Behroozi et al. 2012 were published. More recent data and models are better converged.

We can also implement more sophisticated statistical calculations, such as luminosity functions:

# Function to compute the CII luminosity function
bin_edges = np.logspace(4,11,46)
bin_centers = np.sqrt(bin_edges[1:]*bin_edges[:-1])
def cii_lf(LCII,box_edge,h):
    vol = (box_edge/h)**3
    lf_val, _ = np.histogram(LCII, bin_edges)
    lf_val = lf_val / vol / np.diff(np.log10(bin_edges))
    return lf_val

# 'box_edge', and 'h' are two special keywords that Handler objects recognize
# corresponding to the side length of the box and little h
lf_vals, redshift = tng100.snap_stat(stat_function=cii_lf,kwargs=['LCII','box_edge','h'])

Collecting sources from Snapshot 99.  

# Plot the redshift evolution of the CII luminosity function for our model

from matplotlib import colormaps

fig,ax = plt.subplots()
ax.set(
    xlabel='L$_\mathregular{CII}$ [L$_\odot$]',
    ylabel='$\Phi_\mathregular{CII}$ [Mpc$^{-3}$ dex$^{-1}$]',
    xscale='log',
    yscale='log',
    xlim=(1e6,1e11)
)

cmap = colormaps['Spectral']
for i,lf_val in enumerate(lf_vals):
    ax.plot(bin_centers,lf_val,color=cmap(i/len(lf_vals)))
plt.show()