Subhalo face-on vector values

Jean Donet
  • 9 Sep '19

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.

Dylan Nelson
  • 10 Sep '19

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.

Jean Donet
  • 10 Sep '19

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

Dylan Nelson
  • 14 Sep '19

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:

    # 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[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
Jean Donet
  • 16 Sep '19

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

Jean Donet
  • 14 Jan

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.

  • Page 1 of 1