Public Data Access Overview / Discussion Forum

Hello, I am writing in with a quick question regarding a subhalo's face-on orientation. I have been looking through the literature for the Data Specs as well as the web-based API and have not been able to find a way to output the rotation values required to have the subhalo at a face-on orientation. I have noticed that this is included in the subhalo's visualization when I request 'galaxy_stellar_light_faceon' which indicates the rotation needed to output the galaxy face-on. I would like to know if there is a way I can retrieve this information as an output from the galaxy metadata. I am asking this to be able to project the axis vectors to any line of sight and calculate apparent major and minor axes. Thank you for your help.

For reference: I am looking at TNG-100, z = 0 snapshot, subhalo id: 9.

Hi Jean,

The actual values used for these visualizations are not otherwise available, you're right. For this particular subhalo, the rotation matrix is:

matrix([[-0.91728508, -0.35546651, -0.1795316 ],
        [ 0.33745363, -0.93319982,  0.12354423],
        [ 0.21145469, -0.0527417 , -0.97596377]])

Do you need the code in generality? It is based on constructing the moment of inertia tensor for all star-forming gas cells within twice the stellar half mass radius, if there are enough, otherwise the stars within the stellar half mass radius. Diagonalizing this tensor then gives a 3x3 rotation matrix (above) which can bring a flattened system into an edge-on or face-on type of view.

If you would be able to provide me with the existing code for this, I would greatly appreciate it.

Hi Jean,

I think you can get the idea from this snippet:

def momentOfInertiaTensor(sP, subhaloID=None):
    """ Calculate the moment of inertia tensor (3x3 matrix) for a subhalo-scope particle set. """

    # load required particle data for this subhalo
    subhalo = sP.groupCatSingle(subhaloID=subhaloID)
    rHalf = subhalo['SubhaloHalfmassRadType'][sP.ptNum('stars')]
    shPos = subhalo['SubhaloPos']

    gas = sP.snapshotSubset('gas', fields=['mass','pos','sfr'], subhaloID=subhaloID)

    assert gas['count'] > 1

    # distances
    rad_gas = sP.periodicDists(shPos, gas['Coordinates'])

    wGas = np.where( (rad_gas <= 2.0*rHalf) & (gas['StarFormationRate'] > 0.0) )[0]

    masses = gas['Masses'][wGas]
    xyz = gas['Coordinates'][wGas,:]

    # shift
    xyz = np.squeeze(xyz)

    if xyz.ndim == 1:
        xyz = np.reshape( xyz, (1,3) )

    for i in range(3):
        xyz[:,i] -= shPos[i]

    # if coordinates wrapped box boundary before shift:
    sP.correctPeriodicDistVecs(xyz)

    # construct moment of inertia
    I = np.zeros( (3,3), dtype='float32' )

    I[0,0] = np.sum( masses * (xyz[:,1]*xyz[:,1] + xyz[:,2]*xyz[:,2]) )
    I[1,1] = np.sum( masses * (xyz[:,0]*xyz[:,0] + xyz[:,2]*xyz[:,2]) )
    I[2,2] = np.sum( masses * (xyz[:,0]*xyz[:,0] + xyz[:,1]*xyz[:,1]) )
    I[0,1] = -1 * np.sum( masses * (xyz[:,0]*xyz[:,1]) )
    I[0,2] = -1 * np.sum( masses * (xyz[:,0]*xyz[:,2]) )
    I[1,2] = -1 * np.sum( masses * (xyz[:,1]*xyz[:,2]) )
    I[1,0] = I[0,1]
    I[2,0] = I[0,2]
    I[2,1] = I[1,2]

    return I


def rotationMatricesFromInertiaTensor(I):
    """ Calculate 3x3 rotation matrix by a diagonalization of the moment of inertia tensor.
    Note the resultant rotation matrices are hard-coded for projection with axes=[0,1] e.g. along z. """

    # get eigen values and normalized right eigenvectors
    eigen_values, rotation_matrix = np.linalg.eig(I)

    # sort ascending the eigen values
    sort_inds = np.argsort(eigen_values)
    eigen_values = eigen_values[sort_inds]

    # permute the eigenvectors into this order, which is the rotation matrix which orients the
    # principal axes to the cartesian x,y,z axes, such that if axes=[0,1] we have face-on
    new_matrix = np.matrix( (rotation_matrix[:,sort_inds[0]],
                             rotation_matrix[:,sort_inds[1]],
                             rotation_matrix[:,sort_inds[2]]) )

    # make a random edge on view
    phi = np.random.uniform(0, 2*np.pi)
    theta = np.pi / 2
    psi = 0

    A_00 =  np.cos(psi)*np.cos(phi) - np.cos(theta)*np.sin(phi)*np.sin(psi)
    A_01 =  np.cos(psi)*np.sin(phi) + np.cos(theta)*np.cos(phi)*np.sin(psi)
    A_02 =  np.sin(psi)*np.sin(theta)
    A_10 = -np.sin(psi)*np.cos(phi) - np.cos(theta)*np.sin(phi)*np.cos(psi)
    A_11 = -np.sin(psi)*np.sin(phi) + np.cos(theta)*np.cos(phi)*np.cos(psi)
    A_12 =  np.cos(psi)*np.sin(theta)
    A_20 =  np.sin(theta)*np.sin(phi)
    A_21 = -np.sin(theta)*np.cos(phi)
    A_22 =  np.cos(theta)

    random_edgeon_matrix = np.matrix( ((A_00, A_01, A_02), (A_10, A_11, A_12), (A_20, A_21, A_22)) )

    # prepare return with a few other useful versions of this rotation matrix
    r = {}
    r['face-on'] = new_matrix
    r['edge-on'] = np.matrix( ((1,0,0),(0,0,1),(0,-1,0)) ) * r['face-on'] # disk along x-hat
    r['edge-on-smallest'] = np.matrix( ((0,1,0),(0,0,1),(1,0,0)) ) * r['face-on']
    r['edge-on-y'] = np.matrix( ((0,0,1),(1,0,0),(0,-1,0)) ) * r['face-on'] # disk along y-hat
    r['edge-on-random'] = random_edgeon_matrix * r['face-on']
    r['phi'] = phi
    r['identity'] = np.matrix( np.identity(3) )

    return r

This is exactly what I need, thank you for your help Dylan.

Hello Dylan,
I recently revisited the code after a hiatus and saw that the above code snippet seems to be deprecated. Would you be able to provide me with the updated version of the above code to return the rotation matrices from the mass of inertia tensor? Thank you for your help.

Dear Dylan,

When I took the determinant of that actual numerical rotation matrix you posted above, I got minus 1. From what I read, e.g., https://en.wikipedia.org/wiki/Rotation_matrix, that quality defines an "improper rotation," combining "proper rotations with reflections (which invert orientation)." I would think, naively, that I only need to rotate the figure to see it face on or edge on. Is this particular matrix a special case, or am I not thinking about the mechanics correctly?

Thanks.

-- Denis

Perhaps this is just a consequence of how I construct rotation matrices with the above code. You could certainly do it in another way, the code is just meant as an example, which works well enough in practice.

Thanks for the info and the rapid response. Yes, that’s what I figured, but I wanted to be sure (my software was confused). Thanks again. — dfc