Public Data Access Overview / Discussion Forum

I checked myself and just wanted to quickly confirm, is there a 90 percent radius of any sort within IllustrisTNG? 90 percent mass or 90 percent light would work.

That is correct, you would need to compute other radii like r10 or r90 yourself (from the particle data).

Hi Dylan, how would I calculate r90 from the particle data? Is there a formula for this? Also, is there a 50 percent light radius within IllustrisTNG?
Thanks,
Julian

Following the definition, you want the radius enclosing 90% of the "total" mass, for example. You should be careful to define "total".

You can compute the distance of every particle from the center of a galaxy, then compute the cumulative cum of mass as a function of enclosed distance, then find the radius as above.

Hi Dylan,
Thanks for the comment.
I have made some progress following the procedure you describe and task #6 in the API docs (https://www.tng-project.org/data/docs/api/#reference). My code is below. I am testing my code by seeing if it can reproduce r50 in the subhalo data: (https://www.tng-project.org/api/TNG100-1/snapshots/99/subhalos/447698/). However, there are several problems with it, which I list in no particular order:

• I can get it to give gas particles instead of stellar particles by changing “stars” with “gas” in the params (line 8) and then changing to “PartType0” in lines 20-23. I can’t get dark particles (changing to “dm” or “total” doesn’t work). But I don’t think it matters for this application as I would rather have r50 and r90 as light radii than mass radii. Is there a way to get these measures in IllustrisTNG?
• The total stellar/gas mass for the subhalo matches the one given in the link above but it slightly off. Why?
• The r50 value is extremely far off and I don’t know why. Do you have any ideas?

Ideally, I would do this with the API version as I don’t really want to download all the h5py files but will do so if necessary.

Thanks,
Julian

CODE:

import h5py
import numpy as np

id = 447698 #eventually I will loop over many subhalos

params = {'stars':'Coordinates,Masses'}


url = "http://www.tng-project.org/api/TNG100-1/snapshots/z=0" + "/subhalos/" + str(id)
sub = get(url) # get json response of subhalo properties
saved_filename = get(url + "/cutout.hdf5", params) # get and save HDF5 cutout file

rr100=0


with h5py.File(saved_filename) as g:
    dx100 = g['PartType4']['Coordinates'][:,0] - sub['pos_x']
    dy100 = g['PartType4']['Coordinates'][:,1] - sub['pos_y']
    dz100 = g['PartType4']['Coordinates'][:,2] - sub['pos_z']
    mass100 = g['PartType4']['Masses'][:]

    rr100 = np.sqrt(dx100**2 + dy100**2 + dz100**2)

cum_mass=0
radius=0

mass100_sorted_gas = np.sort(mass100_gas)
rr100_sorted_gas = np.sort(rr100_gas)

for i in range(len(rr100)):
            if cum_mass <= sub['massinhalfrad_stars']: #alternatively, this is np.sum(mass100): 

        radius = rr100_sorted_gas[i]
        cum_mass += mass100_sorted_gas[i]     
print("radius:", radius)
print("mass:", cum_mass)

You are sorting the masses and radii separately, which will give the wrong answer.

You should use np.argsort() instead, and also I would suggest np.cumsum().

In general you need to handle the periodic boundary conditions, if the subhalo position is near the box edges.

Hi Dylan,
Thanks for the suggestions. I had forgotten about those functions but they definitely help. I was able to reproduce many of the mass and radii values, but I still can't get r50 for gas. What is the problem? Also, I attempted to handle the periodic boundary conditions; did I do that correctly (I am using IllustrisTNG-100, which I think I read has a box size of 110.7 Mpc, which I set as the box size for the boundary conditions)?

More importantly, how do I get r50 and r90 as light radii instead of mass radii?

My code is below:

Thanks,
Julian

import h5py
import numpy as np

id = 447698
params = {'gas':'Coordinates,Masses', 'stars':'Coordinates,Masses' }




url = "http://www.tng-project.org/api/TNG100-1/snapshots/z=0" + "/subhalos/" + str(id)
sub = get(url) # get json response of subhalo properties
saved_filename = get(url + "/cutout.hdf5", params) # get and save HDF5 cutout file

rr100_gas=0
rr100_stars=0



#######  gas is PartType0
####### stars is PartType4



####### boxsize=110.7 Mpc ? 

####### scale_factor=1 at z=0, so I want boxsize* 1000/0.704 for Mpc -> kpc/h


boxsize = 110.7 \* 0.704 \* 1000  

x_pos=sub['pos_x']
y_pos=sub['pos_y']
z_pos=sub['pos_z']


with h5py.File(saved_filename) as g:

    dx100_gas = g['PartType0']['Coordinates'][:,0] - x_pos
    dy100_gas = g['PartType0']['Coordinates'][:,1] - y_pos
    dz100_gas = g['PartType0']['Coordinates'][:,2] - z_pos
    mass100_gas = g['PartType0']['Masses'][:]


    dx100_stars = g['PartType4']['Coordinates'][:,0] - x_pos
    dy100_stars = g['PartType4']['Coordinates'][:,1] - y_pos
    dz100_stars = g['PartType4']['Coordinates'][:,2] - z_pos
    mass100_stars = g['PartType4']['Masses'][:]


    for i in range(len(dx100_stars)):
        if (dx100_stars[i] <  -boxsize * 0.5): dx100_stars[i] = dx100_stars[i] + boxsize
        if (dx100_stars[i] >=  boxsize * 0.5): dx100_stars[i] = dx100_stars[i] - boxsize

        if (dy100_stars[i] <  -boxsize * 0.5): dy100_stars[i] = dy100_stars[i] + boxsize
        if (dy100_stars[i] >=  boxsize * 0.5): dy100_stars[i] = dy100_stars[i] - boxsize

        if (dz100_stars[i] <  -boxsize * 0.5): dz100_stars[i] = dz100_stars[i] + boxsize
        if (dz100_stars[i] >=  boxsize * 0.5): dz100_stars[i] = dz100_stars[i] - boxsize



    for i in range(len(dx100_gas)):
        if (dx100_gas[i] <  -boxsize * 0.5): dx100_gas[i] = dx100_gas[i] + boxsize
        if (dx100_gas[i] >=  boxsize * 0.5): dx100_gas[i] = dx100_gas[i] - boxsize

        if (dy100_gas[i] <  -boxsize * 0.5): dy100_gas[i] = dy100_gas[i] + boxsize
        if (dy100_gas[i] >=  boxsize * 0.5): dy100_gas[i] = dy100_gas[i] - boxsize

        if (dz100_gas[i] <  -boxsize * 0.5): dz100_gas[i] = dz100_gas[i] + boxsize
        if (dz100_gas[i] >=  boxsize * 0.5): dz100_gas[i] = dz100_gas[i] - boxsize



    rr100_gas = np.sqrt(dx100_gas**2 + dy100_gas**2 + dz100_gas**2)
    rr100_stars = np.sqrt(dx100_stars**2 + dy100_stars**2 + dz100_stars**2)


####### print("done!")

sorted_along_radius_stars=rr100_stars.argsort()

mass100_sorted_stars = mass100_stars[sorted_along_radius_stars]
rr100_sorted_stars = rr100_stars[sorted_along_radius_stars]

sorted_along_radius_gas=rr100_gas.argsort()

mass100_sorted_gas = mass100_gas[sorted_along_radius_gas]
rr100_sorted_gas = rr100_gas[sorted_along_radius_gas]

print("stellar half mass radius", rr100_sorted_stars[np.where(np.cumsum(mass100_sorted_stars) <= sub['massinhalfrad_stars']  )][-1])

print("gas half mass radius",   rr100_sorted_gas[   np.where(np.cumsum(  mass100_sorted_gas)  <= sub['massinhalfrad_gas'] )][-1])


print("***")

print("gas mass in stellar half mass radius:", np.cumsum(mass100_sorted_gas)[np.where(rr100_sorted_gas <= sub['halfmassrad_stars'])] [-1])


print("gas mass in twice half stellar mass radius:", np.cumsum(mass100_sorted_gas)[np.where(rr100_sorted_gas <= 2*sub['halfmassrad_stars'])] [-1])

print("***")
print("stellar mass in half stellar mass radius:", np.cumsum(mass100_sorted_stars)[np.where(rr100_sorted_stars <= sub['halfmassrad_stars'])] [-1])


print("stellar in twice half stellar mass radius:", np.cumsum(mass100_sorted_stars)[np.where(rr100_sorted_stars <= 2*sub['halfmassrad_stars'])] [-1])

print("***")
print("stellar mass in vmax radius:", np.cumsum(mass100_sorted_stars)[np.where(rr100_sorted_stars <= sub['vmaxrad'] )] [-1])
print("gas mass in vmax radius:", np.cumsum(mass100_sorted_gas)[np.where(rr100_sorted_gas <= sub['vmaxrad'] )] [-1])

The periodic BC seems ok, but you can use BoxSize = 75000.0 exactly (in "code units" of ckpc/h). Note you can use np.where() instead of any loops to do the periodic BC (faster).

I don't see any problem with the gas, perhaps there is an issue in how you treat wind particles. Wind particles are PartType4 but have negative ages (please see documentation). They are more correctly called gas, and not stars, and so the half mass radii in the catalogs likely

If you want half-light radii, you need to replace the mass arrays by some estimate of "light". For stars you could use one of the bands in GFM_StellarPhotometrics. For gas, you need to use some model to estimate its emission in whatever band/wavelength/line you are thinking of.

Hi Dylan,
I'm sorry for the very late reply but thanks so much; that helps a lot!

Julian