[yt-users] spin parameter quantity

[Andrew Davis]

Hi folks,

I'm looking at the spin of some objects and have a few questions about the
quantities BaryonSpinParameter, and ParticleSpinParameter.  I've only measured
the spin of DM particles in the past, which is straightfoward compared to the
baryons.  In looking at the code (posted below from derived_quantities.py)
however, I don't think either of the two functions actually calculate the spin
parameter as usually defined: spin = | J | sqrt( | Etot | ) / (G M^5/2).

In the both the Baryon and Particle Spin Parameter functions, the energy
calculated is only the kinetic energy (and it's missing a factor of 0.5),
whereas the definition of spin includes the total energy -- for DM this is just
potential and kinetic energies.  I'm assuming for baryons, it is potential +
kinetic + internal energy.  Neither functions include the potential energy, or
the internal energy.

Also, it looks like the second function (ParticleSpinParameter) really is just
calculating the spin paramter of the DM particles (or star particles I guess),
so it may make more sense just to not include the gas mass in the first line,
and then re-word the short description so that it says that it really is just
calculating the particle spin, ignoring any gas.

If I'm wrong, please let me know -- otherwise, I'd be happy to submit a modified
version of the two functions.

Cheers,

Andrew


def _BaryonSpinParameter(data):
    """
    This function returns the spin parameter for the baryons, but it uses
    the particles in calculating enclosed mass.
    """
    m_enc = data["CellMassMsun"].sum() + data["ParticleMassMsun"].sum()
    am = data["SpecificAngularMomentum"]*data["CellMassMsun"]
    j_mag = am.sum(axis=1)
    e_term_pre = na.sum(data["CellMassMsun"]*data["VelocityMagnitude"]**2.0)
    weight=data["CellMassMsun"].sum()
    return j_mag, m_enc, e_term_pre, weight
def _combBaryonSpinParameter(data, j_mag, m_enc, e_term_pre, weight):
    # Because it's a vector field, we have to ensure we have enough dimensions
    if len(j_mag.shape) < 2: j_mag = na.expand_dims(j_mag, 0)
    W = weight.sum()
    M = m_enc.sum()
    J = na.sqrt(((j_mag.sum(axis=0))**2.0).sum())/W
    E = na.sqrt(e_term_pre.sum()/W)
    G = 6.67e-8 # cm^3 g^-1 s^-2
    spin = J * E / (M*1.989e33*G)
    return spin
add_quantity("BaryonSpinParameter", function=_BaryonSpinParameter,
             combine_function=_combBaryonSpinParameter, n_ret=4)

def _ParticleSpinParameter(data):
    """
    This function returns the spin parameter for the baryons, but it uses
    the particles in calculating enclosed mass.
    """
    m_enc = data["CellMassMsun"].sum() + data["ParticleMassMsun"].sum()
    am = data["ParticleSpecificAngularMomentum"]*data["ParticleMassMsun"]
    if am.size == 0: return (na.zeros((3,), dtype='float64'), m_enc, 0, 0)
    j_mag = am.sum(axis=1)
    e_term_pre = na.sum(data["ParticleMassMsun"]
                       *data["ParticleVelocityMagnitude"]**2.0)
    weight=data["ParticleMassMsun"].sum()
    return j_mag, m_enc, e_term_pre, weight
add_quantity("ParticleSpinParameter", function=_ParticleSpinParameter,
             combine_function=_combBaryonSpinParameter, n_ret=4)

*************************************************************
**  Andrew J. Davis             andrew.davis at yale.edu      **
**  Dept. of Astronomy          203-432-5119               **
**  Yale University             www.astro.yale.edu/adavis  **
*************************************************************