[Yt-svn] yt: Adding angular momentum component fields for particles. spi...

[hg at spacepope.org]

hg Repository: yt
details:   yt/rev/e7a157d48788
changeset: 3750:e7a157d48788
user:      John Wise <jwise at astro.princeton.edu>
date:
Sat Feb 19 12:37:21 2011 -0500
description:
Adding angular momentum component fields for particles.  spin derived
quantities now use the AM component fields, so they work in parallel.

diffstat:

 yt/data_objects/derived_quantities.py |  20 ++++++++----
 yt/data_objects/universal_fields.py   |  54 +++++++++++++++++++++++++++++++++++
 2 files changed, 67 insertions(+), 7 deletions(-)

diffs (115 lines):

diff -r 8d6cef145d37 -r e7a157d48788 yt/data_objects/derived_quantities.py
--- a/yt/data_objects/derived_quantities.py	Fri Feb 18 14:41:21 2011 -0700
+++ b/yt/data_objects/derived_quantities.py	Sat Feb 19 12:37:21 2011 -0500
@@ -215,8 +215,10 @@
     """
     This function returns the mass-weighted average angular momentum vector.
     """
-    am = data["SpecificAngularMomentum"]*data["CellMassMsun"]
-    j_mag = am.sum(axis=1)
+    amx = data["SpecificAngularMomentumX"]*data["CellMassMsun"]
+    amy = data["SpecificAngularMomentumY"]*data["CellMassMsun"]
+    amz = data["SpecificAngularMomentumZ"]*data["CellMassMsun"]
+    j_mag = [amx.sum(), amy.sum(), amz.sum()]
     return [j_mag]
 def _combAngularMomentumVector(data, j_mag):
     if len(j_mag.shape) < 2: j_mag = na.expand_dims(j_mag, 0)
@@ -232,8 +234,10 @@
     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)
+    amx = data["SpecificAngularMomentumX"]*data["CellMassMsun"]
+    amy = data["SpecificAngularMomentumY"]*data["CellMassMsun"]
+    amz = data["SpecificAngularMomentumZ"]*data["CellMassMsun"]
+    j_mag = na.array([amx.sum(), amy.sum(), amz.sum()])
     e_term_pre = na.sum(data["CellMassMsun"]*data["VelocityMagnitude"]**2.0)
     weight=data["CellMassMsun"].sum()
     return j_mag, m_enc, e_term_pre, weight
@@ -256,9 +260,11 @@
     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)
+    amx = data["ParticleSpecificAngularMomentumX"]*data["ParticleMassMsun"]
+    if amx.size == 0: return (na.zeros((3,), dtype='float64'), m_enc, 0, 0)
+    amy = data["ParticleSpecificAngularMomentumY"]*data["ParticleMassMsun"]
+    amz = data["ParticleSpecificAngularMomentumZ"]*data["ParticleMassMsun"]
+    j_mag = na.array([amx.sum(), amy.sum(), amz.sum()])
     e_term_pre = na.sum(data["ParticleMassMsun"]
                        *data["ParticleVelocityMagnitude"]**2.0)
     weight=data["ParticleMassMsun"].sum()
diff -r 8d6cef145d37 -r e7a157d48788 yt/data_objects/universal_fields.py
--- a/yt/data_objects/universal_fields.py	Fri Feb 18 14:41:21 2011 -0700
+++ b/yt/data_objects/universal_fields.py	Sat Feb 19 12:37:21 2011 -0500
@@ -608,6 +608,43 @@
           function=_ParticleSpecificAngularMomentum, particle_type=True,
           convert_function=_convertSpecificAngularMomentumKMSMPC, vector_field=True,
           units=r"\rm{km}\rm{Mpc}/\rm{s}", validators=[ValidateParameter('center')])
+
+def _ParticleSpecificAngularMomentumX(field, data):
+    if data.has_field_parameter("bulk_velocity"):
+        bv = data.get_field_parameter("bulk_velocity")
+    else: bv = na.zeros(3, dtype='float64')
+    center = data.get_field_parameter('center')
+    y = data["particle_position_y"] - center[1]
+    z = data["particle_position_z"] - center[2]
+    yv = data["particle_velocity_y"] - bv[1]
+    zv = data["particle_velocity_z"] - bv[2]
+    return yv*z - zv*y
+def _ParticleSpecificAngularMomentumY(field, data):
+    if data.has_field_parameter("bulk_velocity"):
+        bv = data.get_field_parameter("bulk_velocity")
+    else: bv = na.zeros(3, dtype='float64')
+    center = data.get_field_parameter('center')
+    x = data["particle_position_x"] - center[0]
+    z = data["particle_position_z"] - center[2]
+    xv = data["particle_velocity_x"] - bv[0]
+    zv = data["particle_velocity_z"] - bv[2]
+    return -(xv*z - zv*x)
+def _ParticleSpecificAngularMomentumZ(field, data):
+    if data.has_field_parameter("bulk_velocity"):
+        bv = data.get_field_parameter("bulk_velocity")
+    else: bv = na.zeros(3, dtype='float64')
+    center = data.get_field_parameter('center')
+    x = data["particle_position_x"] - center[0]
+    y = data["particle_position_y"] - center[1]
+    xv = data["particle_velocity_x"] - bv[0]
+    yv = data["particle_velocity_y"] - bv[1]
+    return xv*y - yv*x
+for ax in 'XYZ':
+    n = "ParticleSpecificAngularMomentum%s" % ax
+    add_field(n, function=eval("_%s" % n), particle_type=True,
+              convert_function=_convertSpecificAngularMomentum,
+              units=r"\rm{cm}^2/\rm{s}", validators=[ValidateParameter("center")])
+
 def _ParticleAngularMomentum(field, data):
     return data["ParticleMass"] * data["ParticleSpecificAngularMomentum"]
 add_field("ParticleAngularMomentum",
@@ -620,6 +657,23 @@
           units=r"M_{\odot}\rm{km}\rm{Mpc}/\rm{s}",
           particle_type=True, validators=[ValidateParameter('center')])
 
+def _ParticleAngularMomentumX(field, data):
+    return data["CellMass"] * data["ParticleSpecificAngularMomentumX"]
+add_field("ParticleAngularMomentumX", function=_ParticleAngularMomentumX,
+         units=r"\rm{g}\/\rm{cm}^2/\rm{s}", particle_type=True,
+         validators=[ValidateParameter('center')])
+def _ParticleAngularMomentumY(field, data):
+    return data["CellMass"] * data["ParticleSpecificAngularMomentumY"]
+add_field("ParticleAngularMomentumY", function=_ParticleAngularMomentumY,
+         units=r"\rm{g}\/\rm{cm}^2/\rm{s}", particle_type=True,
+         validators=[ValidateParameter('center')])
+def _ParticleAngularMomentumZ(field, data):
+    return data["CellMass"] * data["ParticleSpecificAngularMomentumZ"]
+add_field("ParticleAngularMomentumZ", function=_ParticleAngularMomentumZ,
+         units=r"\rm{g}\/\rm{cm}^2/\rm{s}", particle_type=True,
+         validators=[ValidateParameter('center')])
+
+
 def _ParticleRadius(field, data):
     center = data.get_field_parameter("center")
     DW = data.pf.domain_right_edge - data.pf.domain_left_edge