[yt-svn] commit/yt: xarthisius: Merged in cosmosquark/yt (pull request #1585)
commits-noreply at bitbucket.org
commits-noreply at bitbucket.org
Wed Jun 3 09:39:52 PDT 2015
1 new commit in yt:
https://bitbucket.org/yt_analysis/yt/commits/8d3663ac217d/
Changeset: 8d3663ac217d
Branch: yt
User: xarthisius
Date: 2015-06-03 16:39:43+00:00
Summary: Merged in cosmosquark/yt (pull request #1585)
[bugfix][improvement] Specific Angular Momentum [xyz] computed relative to a normal vector and some bug fixes
Affected #: 3 files
diff -r f000bf5609433df69390a45396dc3fe60c87b0b5 -r 8d3663ac217ddab3405fc79be7e7b299aa09bbd0 yt/fields/particle_fields.py
--- a/yt/fields/particle_fields.py
+++ b/yt/fields/particle_fields.py
@@ -23,7 +23,8 @@
ValidateSpatial
from yt.units.yt_array import \
- uconcatenate
+ uconcatenate, \
+ ucross
from yt.utilities.math_utils import \
get_sph_r_component, \
@@ -118,10 +119,8 @@
units = "g")
def particle_density(field, data):
- pos = data[ptype, coord_name]
- mass = data[ptype, mass_name]
- pos.convert_to_units("code_length")
- mass.convert_to_units("code_mass")
+ pos = data[ptype, coord_name].convert_to_units("code_length")
+ mass = data[ptype, mass_name].convert_to_units("code_mass")
d = data.deposit(pos, [data[ptype, mass_name]], method = "sum")
d = data.ds.arr(d, "code_mass")
d /= data["index", "cell_volume"]
@@ -244,17 +243,26 @@
units = "cm / s",
particle_type=True)
+def get_angular_momentum_components(ptype, data, spos, svel):
+ if data.has_field_parameter("normal"):
+ normal = data.get_field_parameter("normal")
+ else:
+ normal = data.ds.arr([0.0,0.0,1.0],"code_length") # default to simulation axis
+ bv = data.get_field_parameter("bulk_velocity")
+ pos = data.ds.arr([data[ptype, spos % ax] for ax in "xyz"]).T
+ vel = data.ds.arr([data[ptype, svel % ax] - bv[iax] for iax, ax in enumerate("xyz")]).T
+ return pos, vel, normal, bv
+
def standard_particle_fields(registry, ptype,
spos = "particle_position_%s",
svel = "particle_velocity_%s"):
# This function will set things up based on the scalar fields and standard
# yt field names.
+ # data.get_field_parameter("bulk_velocity") defaults to YTArray([0,0,0] cm/s)
def _particle_velocity_magnitude(field, data):
""" M{|v|} """
bulk_velocity = data.get_field_parameter("bulk_velocity")
- if bulk_velocity is None:
- bulk_velocity = np.zeros(3)
return np.sqrt((data[ptype, svel % 'x'] - bulk_velocity[0])**2
+ (data[ptype, svel % 'y'] - bulk_velocity[1])**2
+ (data[ptype, svel % 'z'] - bulk_velocity[2])**2 )
@@ -270,19 +278,13 @@
Calculate the angular of a particle velocity. Returns a vector for each
particle.
"""
- if data.has_field_parameter("bulk_velocity"):
- bv = data.get_field_parameter("bulk_velocity")
- else: bv = np.zeros(3, dtype=np.float64)
- xv = data[ptype, svel % 'x'] - bv[0]
- yv = data[ptype, svel % 'y'] - bv[1]
- zv = data[ptype, svel % 'z'] - bv[2]
center = data.get_field_parameter('center')
- coords = self.ds.arr([data[ptype, spos % d] for d in 'xyz'],
- dtype=np.float64, units='cm')
- new_shape = tuple([3] + [1]*(len(coords.shape)-1))
- r_vec = coords - np.reshape(center,new_shape)
- v_vec = self.ds.arr([xv,yv,zv], dtype=np.float64, units='cm/s')
- return np.cross(r_vec, v_vec, axis=0).T
+ pos, vel, normal, bv = get_angular_momentum_components(ptype, data, spos, svel)
+ L, r_vec, v_vec = modify_reference_frame(center, normal, P=pos, V=vel)
+ # adding in the unit registry allows us to have a reference to the dataset
+ # and thus we will always get the correct units after applying the cross product.
+ return -ucross(r_vec, v_vec, registry=data.ds.unit_registry)
+
registry.add_field((ptype, "particle_specific_angular_momentum"),
function=_particle_specific_angular_momentum,
@@ -290,87 +292,29 @@
units="cm**2/s",
validators=[ValidateParameter("center")])
- def _particle_specific_angular_momentum_x(field, data):
- if data.has_field_parameter("bulk_velocity"):
- bv = data.get_field_parameter("bulk_velocity")
- else: bv = np.zeros(3, dtype=np.float64)
- center = data.get_field_parameter('center')
- y = data[ptype, spos % "y"] - center[1]
- z = data[ptype, spos % "z"] - center[2]
- yv = data[ptype, svel % "y"] - bv[1]
- zv = data[ptype, svel % "z"] - bv[2]
- return yv*z - zv*y
-
- registry.add_field((ptype, "particle_specific_angular_momentum_x"),
- function=_particle_specific_angular_momentum_x,
- particle_type=True,
- units="cm**2/s",
- validators=[ValidateParameter("center")])
-
- def _particle_specific_angular_momentum_y(field, data):
- if data.has_field_parameter("bulk_velocity"):
- bv = data.get_field_parameter("bulk_velocity")
- else: bv = np.zeros(3, dtype=np.float64)
- center = data.get_field_parameter('center')
- x = data[ptype, spos % "x"] - center[0]
- z = data[ptype, spos % "z"] - center[2]
- xv = data[ptype, svel % "x"] - bv[0]
- zv = data[ptype, svel % "z"] - bv[2]
- return -(xv*z - zv*x)
-
- registry.add_field((ptype, "particle_specific_angular_momentum_y"),
- function=_particle_specific_angular_momentum_y,
- particle_type=True,
- units="cm**2/s",
- validators=[ValidateParameter("center")])
-
- def _particle_specific_angular_momentum_z(field, data):
- if data.has_field_parameter("bulk_velocity"):
- bv = data.get_field_parameter("bulk_velocity")
- else: bv = np.zeros(3, dtype=np.float64)
- center = data.get_field_parameter('center')
- x = data[ptype, spos % "x"] - center[0]
- y = data[ptype, spos % "y"] - center[1]
- xv = data[ptype, svel % "x"] - bv[0]
- yv = data[ptype, svel % "y"] - bv[1]
- return xv*y - yv*x
-
- registry.add_field((ptype, "particle_specific_angular_momentum_z"),
- function=_particle_specific_angular_momentum_z,
- particle_type=True,
- units="cm**2/s",
- validators=[ValidateParameter("center")])
-
- create_magnitude_field(registry, "particle_specific_angular_momentum",
- "cm**2/s", ftype=ptype, particle_type=True)
-
- def _particle_angular_momentum_x(field, data):
- return data[ptype, "particle_mass"] * \
- data[ptype, "particle_specific_angular_momentum_x"]
- registry.add_field((ptype, "particle_angular_momentum_x"),
- function=_particle_angular_momentum_x,
- units="g*cm**2/s", particle_type=True,
- validators=[ValidateParameter('center')])
-
- def _particle_angular_momentum_y(field, data):
- return data[ptype, "particle_mass"] * \
- data[ptype, "particle_specific_angular_momentum_y"]
- registry.add_field((ptype, "particle_angular_momentum_y"),
- function=_particle_angular_momentum_y,
- units="g*cm**2/s", particle_type=True,
- validators=[ValidateParameter('center')])
-
- def _particle_angular_momentum_z(field, data):
- return data[ptype, "particle_mass"] * \
- data[ptype, "particle_specific_angular_momentum_z"]
- registry.add_field((ptype, "particle_angular_momentum_z"),
- function=_particle_angular_momentum_z,
- units="g*cm**2/s", particle_type=True,
- validators=[ValidateParameter('center')])
+ def _get_spec_ang_mom_comp(axi, ax, _ptype):
+ def _particle_specific_angular_momentum_component(field, data):
+ return data[_ptype, "particle_specific_angular_momentum"][:, axi]
+ def _particle_angular_momentum_component(field, data):
+ return data[_ptype, "particle_mass"] * \
+ data[ptype, "particle_specific_angular_momentum_%s" % ax]
+ return _particle_specific_angular_momentum_component, \
+ _particle_angular_momentum_component
+ for axi, ax in enumerate("xyz"):
+ f, v = _get_spec_ang_mom_comp(axi, ax, ptype)
+ registry.add_field(
+ (ptype, "particle_specific_angular_momentum_%s" % ax),
+ particle_type = True, function=f, units="cm**2/s",
+ validators=[ValidateParameter("center")]
+ )
+ registry.add_field((ptype, "particle_angular_momentum_%s" % ax),
+ function=v, units="g*cm**2/s", particle_type=True,
+ validators=[ValidateParameter('center')])
def _particle_angular_momentum(field, data):
- return (data[ptype, "particle_mass"] *
- data[ptype, "particle_specific_angular_momentum"].T).T
+ am = data[ptype, "particle_mass"] * data[ptype, "particle_specific_angular_momentum"].T
+ return am.T
+
registry.add_field((ptype, "particle_angular_momentum"),
function=_particle_angular_momentum,
particle_type=True,
@@ -405,9 +349,7 @@
"""
normal = data.get_field_parameter('normal')
center = data.get_field_parameter('center')
- pos = spos
- pos = YTArray([data[ptype, pos % ax] for ax in "xyz"])
- pos = pos.T
+ pos = data.ds.arr([data[ptype, spos % ax] for ax in "xyz"]).T
L, pos = modify_reference_frame(center, normal, P=pos)
return pos
@@ -418,81 +360,6 @@
units="cm",
validators=[ValidateParameter("normal"), ValidateParameter("center")])
- def _particle_position_relative_x(field, data):
- """The x component of the particle positions in a rotated reference
- frame
-
- Relative to the coordinate system defined by the *normal* vector and
- *center* field parameters.
-
- Note that the orientation of the x and y axes are arbitrary.
- """
- normal = data.get_field_parameter('normal')
- center = data.get_field_parameter('center')
- pos = spos
- pos = YTArray([data[ptype, pos % ax] for ax in "xyz"])
- pos = pos.T
- L, pos, = modify_reference_frame(center, normal, P=pos)
- pos = pos.T
- return pos[0]
-
- registry.add_field(
- (ptype, "particle_position_relative_x"),
- function=_particle_position_relative_x,
- particle_type=True,
- units="cm",
- validators=[ValidateParameter("normal"), ValidateParameter("center")])
-
- def _particle_position_relative_y(field, data):
- """The y component of the particle positions in a rotated reference
- frame
-
- Relative to the coordinate system defined by the *normal* vector and
- *center* field parameters.
-
- Note that the orientation of the x and y axes are arbitrary.
- """
- normal = data.get_field_parameter('normal')
- center = data.get_field_parameter('center')
- pos = spos
- pos = YTArray([data[ptype, pos % ax] for ax in "xyz"])
- pos = pos.T
- L, pos = modify_reference_frame(center, normal, P=pos)
- pos = pos.T
- return pos[1]
-
- registry.add_field((ptype, "particle_position_relative_y"),
- function=_particle_position_relative_y,
- particle_type=True, units="cm",
- validators=[ValidateParameter("normal"),
- ValidateParameter("center")])
-
-
- def _particle_position_relative_z(field, data):
- """The z component of the particle positions in a rotated reference
- frame
-
- Relative to the coordinate system defined by the *normal* vector and
- *center* field parameters.
-
- Note that the orientation of the x and y axes are arbitrary.
- """
- normal = data.get_field_parameter('normal')
- center = data.get_field_parameter('center')
-
- pos = spos
- pos = YTArray([data[ptype, pos % ax] for ax in "xyz"])
- pos = pos.T
- L, pos = modify_reference_frame(center, normal, P=pos)
- pos = pos.T
- return pos[2]
-
- registry.add_field((ptype, "particle_position_relative_z"),
- function=_particle_position_relative_z,
- particle_type=True, units="cm",
- validators=[ValidateParameter("normal"),
- ValidateParameter("center")])
-
def _particle_velocity_relative(field, data):
"""The vector particle velocities in an arbitrary coordinate system
@@ -504,10 +371,7 @@
normal = data.get_field_parameter('normal')
center = data.get_field_parameter('center')
bv = data.get_field_parameter("bulk_velocity")
- vel = svel
- vel = YTArray([data[ptype, vel % ax] for ax in "xyz"])
- vel = vel - np.reshape(bv, (3, 1))
- vel = vel.T
+ vel = data.ds.arr([data[ptype, svel % ax] - bv[iax] for iax, ax in enumerate("xyz")]).T
L, vel = modify_reference_frame(center, normal, V=vel)
return vel
@@ -517,83 +381,22 @@
validators=[ValidateParameter("normal"),
ValidateParameter("center")])
- def _particle_velocity_relative_x(field, data):
- """The x component of the particle velocities in an arbitrary coordinate
- system
- Relative to the coordinate system defined by the *normal* vector,
- *bulk_velocity* vector and *center* field parameters.
+ def _get_coord_funcs_relative(axi, _ptype):
+ def _particle_pos_rel(field, data):
+ return data[_ptype, "particle_position_relative"][:, axi]
+ def _particle_vel_rel(field, data):
+ return data[_ptype, "particle_velocity_relative"][:, axi]
+ return _particle_vel_rel, _particle_pos_rel
+ for axi, ax in enumerate("xyz"):
+ v, p = _get_coord_funcs_relative(axi, ptype)
+ registry.add_field((ptype, "particle_velocity_relative_%s" % ax),
+ particle_type = True, function = v,
+ units = "code_velocity")
+ registry.add_field((ptype, "particle_position_relative_%s" % ax),
+ particle_type = True, function = p,
+ units = "code_length")
- Note that the orientation of the x and y axes are arbitrary.
- """
- normal = data.get_field_parameter('normal')
- center = data.get_field_parameter('center')
- bv = data.get_field_parameter("bulk_velocity")
- vel = svel
- vel = YTArray([data[ptype, vel % ax] for ax in "xyz"])
- vel = vel - np.reshape(bv, (3, 1))
- vel = vel.T
- L, vel = modify_reference_frame(center, normal, V=vel)
- vel = vel.T
- return vel[0]
-
- registry.add_field((ptype, "particle_velocity_relative_x"),
- function=_particle_velocity_relative_x,
- particle_type=True, units="cm/s",
- validators=[ValidateParameter("normal"),
- ValidateParameter("center")])
-
- def _particle_velocity_relative_y(field, data):
- """The y component of the particle velocities in an arbitrary coordinate
- system
-
- Relative to the coordinate system defined by the *normal* vector,
- *bulk_velocity* vector and *center* field parameters.
-
- Note that the orientation of the x and y axes are arbitrary.
- """
- normal = data.get_field_parameter('normal')
- center = data.get_field_parameter('center')
- bv = data.get_field_parameter('bulk_velocity')
- vel = svel
- vel = YTArray([data[ptype, vel % ax] for ax in "xyz"])
- vel = vel - np.reshape(bv, (3, 1))
- vel = vel.T
- L, vel = modify_reference_frame(center, normal, V=vel)
- vel = vel.T
- return vel[1]
-
- registry.add_field((ptype, "particle_velocity_relative_y"),
- function=_particle_velocity_relative_y,
- particle_type=True, units="cm/s",
- validators=[ValidateParameter("normal"),
- ValidateParameter("center")])
-
- def _particle_velocity_relative_z(field, data):
- """The z component of the particle velocities in an arbitrary coordinate
- system
-
- Relative to the coordinate system defined by the *normal* vector,
- *bulk_velocity* vector and *center* field parameters.
-
- Note that the orientation of the x and y axes are arbitrary.
- """
- normal = data.get_field_parameter('normal')
- center = data.get_field_parameter('center')
- bv = data.get_field_parameter("bulk_velocity")
- vel = svel
- vel = YTArray([data[ptype, vel % ax] for ax in "xyz"])
- bv = vel - np.reshape(bv, (3, 1))
- vel = vel.T
- L, vel = modify_reference_frame(center, normal, V=vel)
- vel = vel.T
- return vel[2]
-
- registry.add_field((ptype, "particle_velocity_relative_z"),
- function=_particle_velocity_relative_z,
- particle_type=True, units="cm/s",
- validators=[ValidateParameter("normal"),
- ValidateParameter("center")])
# this is just particle radius but we add it with an alias for the sake of
# consistent naming
@@ -621,8 +424,7 @@
"""
normal = data.get_field_parameter("normal")
center = data.get_field_parameter("center")
- pos = spos
- pos = YTArray([data[ptype, pos % ax] for ax in "xyz"])
+ pos = data.ds.arr([data[ptype, spos % ax] for ax in "xyz"])
pos = pos - np.reshape(center, (3, 1))
return data.ds.arr(get_sph_theta(pos, normal), "")
@@ -651,8 +453,7 @@
"""
normal = data.get_field_parameter("normal")
center = data.get_field_parameter("center")
- pos = spos
- pos = YTArray([data[ptype, pos % ax] for ax in "xyz"])
+ pos = data.ds.arr([data[ptype, spos % ax] for ax in "xyz"])
pos = pos - np.reshape(center, (3, 1))
return data.ds.arr(get_sph_phi(pos, normal), "")
@@ -683,10 +484,8 @@
normal = data.get_field_parameter('normal')
center = data.get_field_parameter('center')
bv = data.get_field_parameter("bulk_velocity")
- pos = spos
- pos = YTArray([data[ptype, pos % ax] for ax in "xyz"])
- vel = svel
- vel = YTArray([data[ptype, vel % ax] for ax in "xyz"])
+ pos = data.ds.arr([data[ptype, spos % ax] for ax in "xyz"])
+ vel = data.ds.arr([data[ptype, svel % ax] for ax in "xyz"])
theta = get_sph_theta(pos, center)
phi = get_sph_phi(pos, center)
pos = pos - np.reshape(center, (3, 1))
@@ -728,10 +527,8 @@
normal = data.get_field_parameter('normal')
center = data.get_field_parameter('center')
bv = data.get_field_parameter("bulk_velocity")
- pos = spos
- pos = YTArray([data[ptype, pos % ax] for ax in "xyz"])
- vel = svel
- vel = YTArray([data[ptype, vel % ax] for ax in "xyz"])
+ pos = data.ds.arr([data[ptype, spos % ax] for ax in "xyz"])
+ vel = data.ds.arr([data[ptype, svel % ax] for ax in "xyz"])
theta = get_sph_theta(pos, center)
phi = get_sph_phi(pos, center)
pos = pos - np.reshape(center, (3, 1))
@@ -765,8 +562,8 @@
normal = data.get_field_parameter('normal')
center = data.get_field_parameter('center')
bv = data.get_field_parameter("bulk_velocity")
- pos = YTArray([data[ptype, spos % ax] for ax in "xyz"])
- vel = YTArray([data[ptype, svel % ax] for ax in "xyz"])
+ pos = data.ds.arr([data[ptype, spos % ax] for ax in "xyz"])
+ vel = data.ds.arr([data[ptype, svel % ax] for ax in "xyz"])
phi = get_sph_phi(pos, center)
pos = pos - np.reshape(center, (3, 1))
vel = vel - np.reshape(bv, (3, 1))
@@ -798,7 +595,7 @@
"""
normal = data.get_field_parameter("normal")
center = data.get_field_parameter('center')
- pos = YTArray([data[ptype, spos % ax] for ax in "xyz"])
+ pos = data.ds.arr([data[ptype, spos % ax] for ax in "xyz"])
pos = pos - np.reshape(center, (3, 1))
return data.ds.arr(get_cyl_r(pos, normal),
'code_length')
@@ -818,7 +615,7 @@
"""
normal = data.get_field_parameter("normal")
center = data.get_field_parameter('center')
- pos = YTArray([data[ptype, spos % ax] for ax in "xyz"])
+ pos = data.ds.arr([data[ptype, spos % ax] for ax in "xyz"])
pos = pos - np.reshape(center, (3, 1))
return data.ds.arr(get_cyl_theta(pos, normal), "")
@@ -837,7 +634,7 @@
"""
normal = data.get_field_parameter("normal")
center = data.get_field_parameter('center')
- pos = YTArray([data[ptype, spos % ax] for ax in "xyz"])
+ pos = data.ds.arr([data[ptype, spos % ax] for ax in "xyz"])
pos = pos - np.reshape(center, (3, 1))
return data.ds.arr(get_cyl_z(pos, normal),
'code_length')
@@ -858,10 +655,8 @@
normal = data.get_field_parameter('normal')
center = data.get_field_parameter('center')
bv = data.get_field_parameter("bulk_velocity")
- pos = spos
- pos = YTArray([data[ptype, pos % ax] for ax in "xyz"])
- vel = svel
- vel = YTArray([data[ptype, vel % ax] for ax in "xyz"])
+ pos = data.ds.arr([data[ptype, spos % ax] for ax in "xyz"])
+ vel = data.ds.arr([data[ptype, svel % ax] for ax in "xyz"])
theta = get_cyl_theta(pos, center)
pos = pos - np.reshape(center, (3, 1))
vel = vel - np.reshape(bv, (3, 1))
@@ -884,10 +679,8 @@
normal = data.get_field_parameter('normal')
center = data.get_field_parameter('center')
bv = data.get_field_parameter("bulk_velocity")
- pos = spos
- pos = YTArray([data[ptype, pos % ax] for ax in "xyz"])
- vel = svel
- vel = YTArray([data[ptype, vel % ax] for ax in "xyz"])
+ pos = data.ds.arr([data[ptype, spos % ax] for ax in "xyz"])
+ vel = data.ds.arr([data[ptype, svel % ax] for ax in "xyz"])
theta = get_cyl_theta(pos, center)
pos = pos - np.reshape(center, (3, 1))
vel = vel - np.reshape(bv, (3, 1))
@@ -920,10 +713,8 @@
normal = data.get_field_parameter('normal')
center = data.get_field_parameter('center')
bv = data.get_field_parameter("bulk_velocity")
- pos = spos
- pos = YTArray([data[ptype, pos % ax] for ax in "xyz"])
- vel = svel
- vel = YTArray([data[ptype, vel % ax] for ax in "xyz"])
+ pos = data.ds.arr([data[ptype, spos % ax] for ax in "xyz"])
+ vel = data.ds.arr([data[ptype, svel % ax] for ax in "xyz"])
pos = pos - np.reshape(center, (3, 1))
vel = vel - np.reshape(bv, (3, 1))
cylz = get_cyl_z_component(vel, normal)
@@ -1044,3 +835,4 @@
units = "g/cm**3")
return [field_name]
+
diff -r f000bf5609433df69390a45396dc3fe60c87b0b5 -r 8d3663ac217ddab3405fc79be7e7b299aa09bbd0 yt/units/yt_array.py
--- a/yt/units/yt_array.py
+++ b/yt/units/yt_array.py
@@ -1224,6 +1224,19 @@
v = validate_numpy_wrapper_units(v, arrs)
return v
+def ucross(arr1,arr2, registry=None):
+ """Applies the cross product to two YT arrays.
+
+ This wrapper around numpy.cross preserves units.
+ See the documentation of numpy.cross for full
+ details.
+ """
+
+ v = np.cross(arr1,arr2)
+ units = arr1.units * arr2.units
+ arr = YTArray(v,units, registry=registry)
+ return arr
+
def uintersect1d(arr1, arr2, assume_unique=False):
"""Find the sorted unique elements of the two input arrays.
diff -r f000bf5609433df69390a45396dc3fe60c87b0b5 -r 8d3663ac217ddab3405fc79be7e7b299aa09bbd0 yt/utilities/math_utils.py
--- a/yt/utilities/math_utils.py
+++ b/yt/utilities/math_utils.py
@@ -301,22 +301,37 @@
[ 0.00000000e+00, 0.00000000e+00, 0.00000000e+00]])
"""
- if (L == np.array([0, 0, 1.])).all():
- # Whew! Nothing to do!
- if V is None:
+ # First translate the positions to center of mass reference frame.
+ if P is not None:
+ P = P - CoM
+
+ # is the L vector pointing in the Z direction?
+ if np.inner(L[:2], L[:2]) == 0.0:
+ # the reason why we need to explicitly check for the above
+ # is that formula is used in denominator
+ # this just checks if we need to flip the z axis or not
+ if L[2] < 0.0:
+ # this is just a simple flip in direction of the z axis
+ if P is not None:
+ P = -P
+ if V is not None:
+ V = -V
+
+ # return the values
+ if V is None and P is not None:
return L, P
- if P is None:
+ elif P is None and V is not None:
return L, V
else:
return L, P, V
- # First translate the positions to center of mass reference frame.
- if P is not None:
- P = P - CoM
+
+ # Normal vector is not aligned with simulation Z axis
+ # Therefore we are going to have to apply a rotation
# Now find the angle between modified L and the x-axis.
LL = L.copy()
- LL[2] = 0.
- theta = np.arccos(np.inner(LL, [1.,0,0])/np.inner(LL,LL)**.5)
- if L[1] < 0:
+ LL[2] = 0.0
+ theta = np.arccos(np.inner(LL, [1.0, 0.0, 0.0]) / np.inner(LL, LL) ** 0.5)
+ if L[1] < 0.0:
theta = -theta
# Now rotate all the position, velocity, and L vectors by this much around
# the z axis.
@@ -326,16 +341,18 @@
V = rotate_vector_3D(V, 2, theta)
L = rotate_vector_3D(L, 2, theta)
# Now find the angle between L and the z-axis.
- theta = np.arccos(np.inner(L, [0,0,1])/np.inner(L,L)**.5)
+ theta = np.arccos(np.inner(L, [0.0, 0.0, 1.0]) / np.inner(L, L) ** 0.5)
# This time we rotate around the y axis.
if P is not None:
P = rotate_vector_3D(P, 1, theta)
if V is not None:
V = rotate_vector_3D(V, 1, theta)
L = rotate_vector_3D(L, 1, theta)
- if V is None:
+
+ # return the values
+ if V is None and P is not None:
return L, P
- if P is None:
+ elif P is None and V is not None:
return L, V
else:
return L, P, V
Repository URL: https://bitbucket.org/yt_analysis/yt/
--
This is a commit notification from bitbucket.org. You are receiving
this because you have the service enabled, addressing the recipient of
this email.
More information about the yt-svn
mailing list