[yt-svn] commit/yt: 14 new changesets
Bitbucket
commits-noreply at bitbucket.org
Mon Oct 15 16:10:40 PDT 2012
14 new commits in yt:
https://bitbucket.org/yt_analysis/yt/changeset/ae7c2028ce27/
changeset: ae7c2028ce27
branch: yt
user: ngoldbaum
date: 2012-10-04 22:37:21
summary: Rearranging the way coordinate transformations are handled in universal_fields. Adding generate coordinate transformation code in math_utils. Adding radial and tangential velocity fields for cylindrical coordinates.
affected #: 2 files
diff -r ac4b312b5d6fac59499a2b8bc5aa6cc3faabd0f2 -r ae7c2028ce277873ac88ae0f7f4c327e6d728988 yt/data_objects/universal_fields.py
--- a/yt/data_objects/universal_fields.py
+++ b/yt/data_objects/universal_fields.py
@@ -194,13 +194,6 @@
function=_TangentialOverVelocityMagnitude,
take_log=False)
-def _TangentialVelocity(field, data):
- return np.sqrt(data["VelocityMagnitude"]**2.0
- - data["RadialVelocity"]**2.0)
-add_field("TangentialVelocity",
- function=_TangentialVelocity,
- take_log=False, units=r"\rm{cm}/\rm{s}")
-
def _Pressure(field, data):
"""M{(Gamma-1.0)*rho*E}"""
return (data.pf["Gamma"] - 1.0) * \
@@ -227,10 +220,7 @@
data['y'] - center[1],
data['z'] - center[2]]).transpose()
- ## The spherical coordinates radius is simply the magnitude of the
- ## coords vector.
-
- return np.sqrt(np.sum(coords**2,axis=-1))
+ return get_sph_r_component(vectors, center)
def _Convert_sph_r_CGS(data):
return data.convert("cm")
@@ -249,16 +239,7 @@
data['y'] - center[1],
data['z'] - center[2]]).transpose()
- ## The angle (theta) with respect to the normal (J), is the arccos
- ## of the dot product of the normal with the normalized coords
- ## vector.
-
- tile_shape = list(coords.shape)[:-1] + [1]
- J = np.tile(normal,tile_shape)
-
- JdotCoords = np.sum(J*coords,axis=-1)
-
- return np.arccos( JdotCoords / np.sqrt(np.sum(coords**2,axis=-1)) )
+ return get_sph_theta_component(vectors, center, normal)
add_field("sph_theta", function=_sph_theta,
validators=[ValidateParameter("center"),ValidateParameter("normal")])
@@ -273,27 +254,7 @@
data['y'] - center[1],
data['z'] - center[2]]).transpose()
- ## We have freedom with respect to what axis (xprime) to define
- ## the disk angle. Here I've chosen to use the axis that is
- ## perpendicular to the normal and the y-axis. When normal ==
- ## y-hat, then set xprime = z-hat. With this definition, when
- ## normal == z-hat (as is typical), then xprime == x-hat.
- ##
- ## The angle is then given by the arctan of the ratio of the
- ## yprime-component and the xprime-component of the coords vector.
-
- xprime = np.cross([0.0,1.0,0.0],normal)
- if np.sum(xprime) == 0: xprime = np.array([0.0, 0.0, 1.0])
- yprime = np.cross(normal,xprime)
-
- tile_shape = list(coords.shape)[:-1] + [1]
- Jx = np.tile(xprime,tile_shape)
- Jy = np.tile(yprime,tile_shape)
-
- Px = np.sum(Jx*coords,axis=-1)
- Py = np.sum(Jy*coords,axis=-1)
-
- return np.arctan2(Py,Px)
+ return get_sph_phi_component(vectors, center, normal)
add_field("sph_phi", function=_sph_phi,
validators=[ValidateParameter("center"),ValidateParameter("normal")])
@@ -309,14 +270,7 @@
data['y'] - center[1],
data['z'] - center[2]]).transpose()
- ## The cross product of the normal (J) with the coords vector
- ## gives a vector of magnitude equal to the cylindrical radius.
-
- tile_shape = list(coords.shape)[:-1] + [1]
- J = np.tile(normal,tile_shape)
-
- JcrossCoords = np.cross(J,coords)
- return np.sqrt(np.sum(JcrossCoords**2,axis=-1))
+ return get_cyl_r_component(coords, center, normal)
def _Convert_cyl_R_CGS(data):
return data.convert("cm")
@@ -324,6 +278,9 @@
add_field("cyl_R", function=_cyl_R,
validators=[ValidateParameter("center"),ValidateParameter("normal")],
convert_function = _Convert_cyl_R_CGS, units=r"\rm{cm}")
+add_field("cyl_RCode", function=_cyl_R,
+ validators=[ValidateParameter("center"),ValidateParameter("normal")],
+ units=r"Radius (code)")
### cylindrical coordinates: z (height above the cylinder's plane)
@@ -335,13 +292,7 @@
data['y'] - center[1],
data['z'] - center[2]]).transpose()
- ## The dot product of the normal (J) with the coords vector gives
- ## the cylindrical height.
-
- tile_shape = list(coords.shape)[:-1] + [1]
- J = np.tile(normal,tile_shape)
-
- return np.sum(J*coords,axis=-1)
+ return get_cyl_z_component(coords,center,normal)
def _Convert_cyl_z_CGS(data):
return data.convert("cm")
@@ -352,14 +303,19 @@
### cylindrical coordinates: theta (angle in the cylinder's plane)
-### [This is identical to the spherical coordinate's 'phi' angle.]
def _cyl_theta(field, data):
- return data['sph_phi']
+ center = data.get_field_parameter("center")
+ normal = data.get_field_parameter("normal")
+
+ coords = np.array([data['x'] - center[0],
+ data['y'] - center[1],
+ data['z'] - center[2]]).transpose()
+
+ return get_cyl_theta_component(coords, center, normal)
add_field("cyl_theta", function=_cyl_theta,
validators=[ValidateParameter("center"),ValidateParameter("normal")])
-
### The old field DiskAngle is the same as the spherical coordinates'
### 'theta' angle. I'm keeping DiskAngle for backwards compatibility.
def _DiskAngle(field, data):
@@ -392,6 +348,65 @@
ValidateParameter("normal")],
units=r"AU", display_field=False)
+def _cyl_RadialVelocity(data):
+ center = data.get_field_parameter("center")
+ normal = data.get_field_parameter("normal")
+ bulk_velocity = data.get_field_parameter("bulk_velocity")
+ if bulk_velocity == None:
+ bulk_velocity = np.zeros(3)
+
+ velocities = np.array([data['x-velocity'] - bulk_velocity[0],
+ data['y-velocity'] - bulk_velocity[1],
+ data['z-velocity'] - bulk_velocity[2]]).transpose()
+
+ return get_cyl_r_component(velocities, center, normal)
+
+def _cyl_RadialVelocityABS(field, data):
+ return np.abs(_cyl_RadialVelocity(field, data))
+def _Convert_cyl_RadialVelocityKMS(data):
+ return 1e-5
+add_field("cyl_RadialVelocity", function=_cyl_RadialVelocity,
+ units=r"\rm{cm}/\rm{s}",
+ validators=[ValidateParameter("center"),ValidateParameter("normal")])
+add_field("cyl_RadialVelocityABS", function=_cyl_RadialVelocityABS,
+ units=r"\rm{cm}/\rm{s}",
+ validators=[ValidateParameter("center"),ValidateParameter("normal")])
+add_field("cyl_RadialVelocityKMS", function=_cyl_RadialVelocity,
+ convert_function=_cyl_ConvertRadialVelocityKMS, units=r"\rm{km}/\rm{s}",
+ validators=[ValidateParameter("center"),ValidateParameter("normal")])
+add_field("cyl_RadialVelocityKMSABS", function=_cyl_RadialVelocityABS,
+ convert_function=_cyl_ConvertRadialVelocityKMS, units=r"\rm{km}/\rm{s}",
+ validators=[ValidateParameter("center"),ValidateParameter("normal")])
+
+def _cyl_TangentialVelocity(data):
+ center = data.get_field_parameter("center")
+ normal = data.get_field_parameter("normal")
+ bulk_velocity = data.get_field_parameter("bulk_velocity")
+ if bulk_velocity == None:
+ bulk_velocity = np.zeros(3)
+
+ velocities = np.array([data['x-velocity'] - bulk_velocity[0],
+ data['y-velocity'] - bulk_velocity[1],
+ data['z-velocity'] - bulk_velocity[2]]).transpose()
+
+ return get_cyl_theta_component(velocities, center, normal)
+
+def _cyl_TangentialVelocityABS(field, data):
+ return np.abs(_cyl_TangentialVelocity(field, data))
+def _Convert_cyl_TangentialVelocityKMS(data):
+ return 1e-5
+add_field("cyl_TangentialVelocity", function=_cyl_TangentialVelocity,
+ units=r"\rm{cm}/\rm{s}",
+ validators=[ValidateParameter("center"),ValidateParameter("normal")])
+add_field("cyl_TangentialVelocityABS", function=_cyl_TangentialVelocityABS,
+ units=r"\rm{cm}/\rm{s}",
+ validators=[ValidateParameter("center"),ValidateParameter("normal")])
+add_field("cyl_TangentialVelocityKMS", function=_cyl_TangentialVelocity,
+ convert_function=_cyl_ConvertTangentialVelocityKMS, units=r"\rm{km}/\rm{s}",
+ validators=[ValidateParameter("center"),ValidateParameter("normal")])
+add_field("cyl_TangentialVelocityKMSABS", function=_cyl_TangentialVelocityABS,
+ convert_function=_cyl_ConvertTangentialVelocityKMS, units=r"\rm{km}/\rm{s}",
+ validators=[ValidateParameter("center"),ValidateParameter("normal")])
def _DynamicalTime(field, data):
"""
@@ -887,13 +902,13 @@
bulk_velocity = data.get_field_parameter("bulk_velocity")
if bulk_velocity == None:
bulk_velocity = np.zeros(3)
- new_field = ( (data['x']-center[0])*(data["x-velocity"]-bulk_velocity[0])
- + (data['y']-center[1])*(data["y-velocity"]-bulk_velocity[1])
- + (data['z']-center[2])*(data["z-velocity"]-bulk_velocity[2])
- )/data["RadiusCode"]
- if np.any(np.isnan(new_field)): # to fix center = point
- new_field[np.isnan(new_field)] = 0.0
- return new_field
+
+ velocities = np.array([data['x-velocity'] - bulk_velocity[0],
+ data['y-velocity'] - bulk_velocity[1],
+ data['z-velocity'] - bulk_velocity[2]]).transpose()
+
+ return get_sph_r_component(velocities, center)
+
def _RadialVelocityABS(field, data):
return np.abs(_RadialVelocity(field, data))
def _ConvertRadialVelocityKMS(data):
@@ -911,6 +926,13 @@
convert_function=_ConvertRadialVelocityKMS, units=r"\rm{km}/\rm{s}",
validators=[ValidateParameter("center")])
+def _TangentialVelocity(field, data):
+ return np.sqrt(data["VelocityMagnitude"]**2.0
+ - data["RadialVelocity"]**2.0)
+add_field("TangentialVelocity",
+ function=_TangentialVelocity,
+ take_log=False, units=r"\rm{cm}/\rm{s}")
+
def _CuttingPlaneVelocityX(field, data):
x_vec, y_vec, z_vec = [data.get_field_parameter("cp_%s_vec" % (ax))
for ax in 'xyz']
diff -r ac4b312b5d6fac59499a2b8bc5aa6cc3faabd0f2 -r ae7c2028ce277873ac88ae0f7f4c327e6d728988 yt/utilities/math_utils.py
--- a/yt/utilities/math_utils.py
+++ b/yt/utilities/math_utils.py
@@ -674,3 +674,69 @@
[uz*ux*(1-cost)-uy*sint, uz*uy*(1-cost)+ux*sint, cost+uz**2*(1-cost)]])
return R
+
+def get_sph_r_component(vectors, center):
+ # The spherical coordinates radius is simply the magnitude of the
+ # vector.
+
+ return np.sqrt(np.sum(vectors**2, axis=-1))
+
+
+def get_sph_theta_component(vectors, center, normal):
+ # The angle (theta) with respect to the normal (J), is the arccos
+ # of the dot product of the normal with the normalized
+ # vector.
+
+ tile_shape = list(vectors.shape)[:-1] + [1]
+ J = np.tile(normal,tile_shape)
+
+ JdotVectors = np.sum(J*vectors,axis=-1)
+
+ return np.arccos( JdotVectors / np.sqrt(np.sum(Vectors**2,axis=-1)) )
+
+def get_sph_phi_component(vectors, center, normal):
+ # We have freedom with respect to what axis (xprime) to define
+ # the disk angle. Here I've chosen to use the axis that is
+ # perpendicular to the normal and the y-axis. When normal ==
+ # y-hat, then set xprime = z-hat. With this definition, when
+ # normal == z-hat (as is typical), then xprime == x-hat.
+ #
+ # The angle is then given by the arctan of the ratio of the
+ # yprime-component and the xprime-component of the vector.
+
+ xprime = np.cross([0.0,1.0,0.0],normal)
+ if np.sum(xprime) == 0: xprime = np.array([0.0, 0.0, 1.0])
+ yprime = np.cross(normal,xprime)
+
+ tile_shape = list(vectors.shape)[:-1] + [1]
+ Jx = np.tile(xprime,tile_shape)
+ Jy = np.tile(yprime,tile_shape)
+
+ Px = np.sum(Jx*vectors,axis=-1)
+ Py = np.sum(Jy*vectors,axis=-1)
+
+ return np.arctan2(Py,Px)
+
+def get_cyl_r_component(vectors, center, normal):
+ # The cross product of the normal (J) with a vector
+ # gives a vector of magnitude equal to the cylindrical radius.
+
+ tile_shape = list(vector.shape)[:-1] + [1]
+ J = np.tile(normal, tile_shape)
+
+ JcrossVectors = np.cross(J, vectors)
+ return np.sqrt(np.sum(JcrossVectors**2, axis=-1))
+
+def get_cyl_z_component(vectors, center, normal):
+ # The dot product of the normal (J) with the vector gives
+ # the cylindrical height.
+
+ tile_shape = list(vectors.shape)[:-1] + [1]
+ J = np.tile(normal, tile_shape)
+
+ return np.sum(J*vectors, axis=-1)
+
+def get_cyl_theta_component(vectors, center, normal):
+ # This is identical to the spherical phi component
+
+ return get_sph_phi_component(vectors, center, normal):
https://bitbucket.org/yt_analysis/yt/changeset/8a9bc98c1f24/
changeset: 8a9bc98c1f24
branch: yt
user: ngoldbaum
date: 2012-10-05 02:05:10
summary: Bugfixes for field coordinate conversion code.
affected #: 2 files
diff -r ae7c2028ce277873ac88ae0f7f4c327e6d728988 -r 8a9bc98c1f248a6f4080cfa35e79521a99b823e7 yt/data_objects/universal_fields.py
--- a/yt/data_objects/universal_fields.py
+++ b/yt/data_objects/universal_fields.py
@@ -55,6 +55,14 @@
G, \
rho_crit_now, \
speed_of_light_cgs
+
+from yt.utilities.math_utils import \
+ get_sph_r_component, \
+ get_sph_theta_component, \
+ get_sph_phi_component, \
+ get_cyl_r_component, \
+ get_cyl_z_component, \
+ get_cyl_theta_component
# Note that, despite my newfound efforts to comply with PEP-8,
# I violate it here in order to keep the name/func_name relationship
@@ -220,7 +228,7 @@
data['y'] - center[1],
data['z'] - center[2]]).transpose()
- return get_sph_r_component(vectors, center)
+ return get_sph_r(vectors, center)
def _Convert_sph_r_CGS(data):
return data.convert("cm")
@@ -239,7 +247,7 @@
data['y'] - center[1],
data['z'] - center[2]]).transpose()
- return get_sph_theta_component(vectors, center, normal)
+ return get_sph_theta(coords, center, normal)
add_field("sph_theta", function=_sph_theta,
validators=[ValidateParameter("center"),ValidateParameter("normal")])
@@ -254,7 +262,7 @@
data['y'] - center[1],
data['z'] - center[2]]).transpose()
- return get_sph_phi_component(vectors, center, normal)
+ return get_sph_phi(coords, center, normal)
add_field("sph_phi", function=_sph_phi,
validators=[ValidateParameter("center"),ValidateParameter("normal")])
@@ -270,7 +278,7 @@
data['y'] - center[1],
data['z'] - center[2]]).transpose()
- return get_cyl_r_component(coords, center, normal)
+ return get_cyl_r(coords, center, normal)
def _Convert_cyl_R_CGS(data):
return data.convert("cm")
@@ -292,7 +300,7 @@
data['y'] - center[1],
data['z'] - center[2]]).transpose()
- return get_cyl_z_component(coords,center,normal)
+ return get_cyl_z(coords, center, normal)
def _Convert_cyl_z_CGS(data):
return data.convert("cm")
@@ -311,7 +319,7 @@
data['y'] - center[1],
data['z'] - center[2]]).transpose()
- return get_cyl_theta_component(coords, center, normal)
+ return get_cyl_theta(coords, center, normal)
add_field("cyl_theta", function=_cyl_theta,
validators=[ValidateParameter("center"),ValidateParameter("normal")])
@@ -348,7 +356,7 @@
ValidateParameter("normal")],
units=r"AU", display_field=False)
-def _cyl_RadialVelocity(data):
+def _cyl_RadialVelocity(field, data):
center = data.get_field_parameter("center")
normal = data.get_field_parameter("normal")
bulk_velocity = data.get_field_parameter("bulk_velocity")
@@ -359,7 +367,11 @@
data['y-velocity'] - bulk_velocity[1],
data['z-velocity'] - bulk_velocity[2]]).transpose()
- return get_cyl_r_component(velocities, center, normal)
+ coords = np.array([data['x'] - center[0],
+ data['y'] - center[1],
+ data['z'] - center[2]]).transpose()
+
+ return get_cyl_r_component(velocities, coords, center, normal)
def _cyl_RadialVelocityABS(field, data):
return np.abs(_cyl_RadialVelocity(field, data))
@@ -372,13 +384,13 @@
units=r"\rm{cm}/\rm{s}",
validators=[ValidateParameter("center"),ValidateParameter("normal")])
add_field("cyl_RadialVelocityKMS", function=_cyl_RadialVelocity,
- convert_function=_cyl_ConvertRadialVelocityKMS, units=r"\rm{km}/\rm{s}",
+ convert_function=_Convert_cyl_RadialVelocityKMS, units=r"\rm{km}/\rm{s}",
validators=[ValidateParameter("center"),ValidateParameter("normal")])
add_field("cyl_RadialVelocityKMSABS", function=_cyl_RadialVelocityABS,
- convert_function=_cyl_ConvertRadialVelocityKMS, units=r"\rm{km}/\rm{s}",
+ convert_function=_Convert_cyl_RadialVelocityKMS, units=r"\rm{km}/\rm{s}",
validators=[ValidateParameter("center"),ValidateParameter("normal")])
-def _cyl_TangentialVelocity(data):
+def _cyl_TangentialVelocity(field, data):
center = data.get_field_parameter("center")
normal = data.get_field_parameter("normal")
bulk_velocity = data.get_field_parameter("bulk_velocity")
@@ -389,7 +401,11 @@
data['y-velocity'] - bulk_velocity[1],
data['z-velocity'] - bulk_velocity[2]]).transpose()
- return get_cyl_theta_component(velocities, center, normal)
+ coords = np.array([data['x'] - center[0],
+ data['y'] - center[1],
+ data['z'] - center[2]]).transpose()
+
+ return get_cyl_theta_component(velocities, coords, center, normal)
def _cyl_TangentialVelocityABS(field, data):
return np.abs(_cyl_TangentialVelocity(field, data))
@@ -402,10 +418,10 @@
units=r"\rm{cm}/\rm{s}",
validators=[ValidateParameter("center"),ValidateParameter("normal")])
add_field("cyl_TangentialVelocityKMS", function=_cyl_TangentialVelocity,
- convert_function=_cyl_ConvertTangentialVelocityKMS, units=r"\rm{km}/\rm{s}",
+ convert_function=_Convert_cyl_TangentialVelocityKMS, units=r"\rm{km}/\rm{s}",
validators=[ValidateParameter("center"),ValidateParameter("normal")])
add_field("cyl_TangentialVelocityKMSABS", function=_cyl_TangentialVelocityABS,
- convert_function=_cyl_ConvertTangentialVelocityKMS, units=r"\rm{km}/\rm{s}",
+ convert_function=_Convert_cyl_TangentialVelocityKMS, units=r"\rm{km}/\rm{s}",
validators=[ValidateParameter("center"),ValidateParameter("normal")])
def _DynamicalTime(field, data):
@@ -899,6 +915,9 @@
def _RadialVelocity(field, data):
center = data.get_field_parameter("center")
+ normal = data.get_field_parameter("normal")
+ if normal == None:
+ normal = [0,0,1]
bulk_velocity = data.get_field_parameter("bulk_velocity")
if bulk_velocity == None:
bulk_velocity = np.zeros(3)
@@ -906,8 +925,12 @@
velocities = np.array([data['x-velocity'] - bulk_velocity[0],
data['y-velocity'] - bulk_velocity[1],
data['z-velocity'] - bulk_velocity[2]]).transpose()
+
+ coords = np.array([data['x'] - center[0],
+ data['y'] - center[1],
+ data['z'] - center[2]]).transpose()
- return get_sph_r_component(velocities, center)
+ return get_sph_r_component(velocities, coords, center, normal)
def _RadialVelocityABS(field, data):
return np.abs(_RadialVelocity(field, data))
diff -r ae7c2028ce277873ac88ae0f7f4c327e6d728988 -r 8a9bc98c1f248a6f4080cfa35e79521a99b823e7 yt/utilities/math_utils.py
--- a/yt/utilities/math_utils.py
+++ b/yt/utilities/math_utils.py
@@ -675,26 +675,33 @@
return R
-def get_sph_r_component(vectors, center):
+def get_ortho_basis(normal):
+ xprime = np.cross([0.0,1.0,0.0],normal)
+ if np.sum(xprime) == 0: xprime = np.array([0.0, 0.0, 1.0])
+ yprime = np.cross(normal,xprime)
+ zprime = normal
+ return (xprime, yprime, zprime)
+
+def get_sph_r(coords, center):
# The spherical coordinates radius is simply the magnitude of the
- # vector.
+ # coordinate vector.
- return np.sqrt(np.sum(vectors**2, axis=-1))
+ return np.sqrt(np.sum(coords**2, axis=-1))
-def get_sph_theta_component(vectors, center, normal):
+def get_sph_theta(coords, center, normal):
# The angle (theta) with respect to the normal (J), is the arccos
- # of the dot product of the normal with the normalized
+ # of the dot product of the normal with the normalized coordinate
# vector.
- tile_shape = list(vectors.shape)[:-1] + [1]
+ tile_shape = list(coords.shape)[:-1] + [1]
J = np.tile(normal,tile_shape)
- JdotVectors = np.sum(J*vectors,axis=-1)
+ JdotCoords = np.sum(J*coords,axis=-1)
- return np.arccos( JdotVectors / np.sqrt(np.sum(Vectors**2,axis=-1)) )
+ return np.arccos( JdotCoords / np.sqrt(np.sum(coords**2,axis=-1)) )
-def get_sph_phi_component(vectors, center, normal):
+def get_sph_phi(coords, center, normal):
# We have freedom with respect to what axis (xprime) to define
# the disk angle. Here I've chosen to use the axis that is
# perpendicular to the normal and the y-axis. When normal ==
@@ -702,41 +709,134 @@
# normal == z-hat (as is typical), then xprime == x-hat.
#
# The angle is then given by the arctan of the ratio of the
- # yprime-component and the xprime-component of the vector.
+ # yprime-component and the xprime-component of the coordinate
+ # vector.
- xprime = np.cross([0.0,1.0,0.0],normal)
- if np.sum(xprime) == 0: xprime = np.array([0.0, 0.0, 1.0])
- yprime = np.cross(normal,xprime)
+ (xprime, yprime, zprime) = get_ortho_basis(normal)
+ tile_shape = list(coords.shape)[:-1] + [1]
+ Jx = np.tile(xprime,tile_shape)
+ Jy = np.tile(yprime,tile_shape)
+
+ Px = np.sum(Jx*coords,axis=-1)
+ Py = np.sum(Jy*coords,axis=-1)
+
+ return np.arctan2(Py,Px)
+
+def get_cyl_r(coords, center, normal):
+ # The cross product of the normal (J) with a coordinate vector
+ # gives a vector of magnitude equal to the cylindrical radius.
+
+ tile_shape = list(coords.shape)[:-1] + [1]
+ J = np.tile(normal, tile_shape)
+
+ JcrossCoords = np.cross(J, coords)
+ return np.sqrt(np.sum(JcrossCoords**2, axis=-1))
+
+def get_cyl_z(coords, center, normal):
+ # The dot product of the normal (J) with the coordinate vector
+ # gives the cylindrical height.
+
+ tile_shape = list(coords.shape)[:-1] + [1]
+ J = np.tile(normal, tile_shape)
+
+ return np.sum(J*coords, axis=-1)
+
+def get_cyl_theta(coords, center, normal):
+ # This is identical to the spherical phi component
+
+ return get_sph_phi(coords, center, normal)
+
+
+def get_cyl_r_component(vectors, coords, center, normal):
+ # The r of a vector is the vector dotted with rhat
+
+ theta = np.tile(get_cyl_theta(coords, center, normal), (3,1)).transpose()
+
+ (xprime, yprime, zprime) = get_ortho_basis(normal)
+
tile_shape = list(vectors.shape)[:-1] + [1]
Jx = np.tile(xprime,tile_shape)
Jy = np.tile(yprime,tile_shape)
+
+ rhat = Jx*np.cos(theta) + Jy*np.sin(theta)
+
+ return np.sum(vectors*rhat,axis=-1)
+
+def get_cyl_theta_component(vectors, coords, center, normal):
+ # The theta component of a vector is the vector dotted with thetahat
- Px = np.sum(Jx*vectors,axis=-1)
- Py = np.sum(Jy*vectors,axis=-1)
-
- return np.arctan2(Py,Px)
+ theta = np.tile(get_cyl_theta(coords, center, normal), (3,1)).transpose()
-def get_cyl_r_component(vectors, center, normal):
- # The cross product of the normal (J) with a vector
- # gives a vector of magnitude equal to the cylindrical radius.
+ (xprime, yprime, zprime) = get_ortho_basis(normal)
- tile_shape = list(vector.shape)[:-1] + [1]
- J = np.tile(normal, tile_shape)
-
- JcrossVectors = np.cross(J, vectors)
- return np.sqrt(np.sum(JcrossVectors**2, axis=-1))
+ tile_shape = list(vectors.shape)[:-1] + [1]
+ Jx = np.tile(xprime,tile_shape)
+ Jy = np.tile(yprime,tile_shape)
+
+ thetahat = -Jx*np.sin(theta) + Jy*np.cos(theta)
+
+ return np.sum(vectors*thetahat, axis=-1)
def get_cyl_z_component(vectors, center, normal):
- # The dot product of the normal (J) with the vector gives
- # the cylindrical height.
+ # The z component of a vector is the vector dotted with zhat
+
+ (xprime, yprime, zprime) = get_ortho_basis(normal)
+
+ tile_shape = list(vectors.shape)[:-1] + [1]
+ zhat = np.tile(zprime, tile_shape)
+
+ return np.sum(vectors*zhat, axis=-1)
+
+def get_sph_r_component(vectors, coords, center, normal):
+ # The r component of a vector is the vector dotted with rhat
+ theta = get_sph_theta(coords, center, normal)
+ phi = get_sph_phi(coords, center, normal)
+
+ (xprime, yprime, zprime) = get_ortho_basis(normal)
+
tile_shape = list(vectors.shape)[:-1] + [1]
- J = np.tile(normal, tile_shape)
+ Jx = np.tile(xprime,tile_shape)
+ Jy = np.tile(yprime,tile_shape)
+ Jz = np.tile(zprime,tile_shape)
- return np.sum(J*vectors, axis=-1)
+ rhat = Jx*np.sin(theta)*np.cos(phi) + \
+ Jy*np.sin(theta)*np.sin(phi) + \
+ Jz*np.cos(theta)
-def get_cyl_theta_component(vectors, center, normal):
- # This is identical to the spherical phi component
+ return np.sum(vectors*rhat, axis=-1)
- return get_sph_phi_component(vectors, center, normal):
+def get_sph_phi_component(vectors, coords, center, normal):
+ # The phi component of a vector is the vector dotted with phihat
+
+ phi = get_sph_phi(coords, center, normal)
+
+ (xprime, yprime, zprime) = get_ortho_basis(normal)
+
+ tile_shape = list(vectors.shape)[:-1] + [1]
+ Jx = np.tile(xprime,tile_shape)
+ Jy = np.tile(yprime,tile_shape)
+
+ phihat = -Jx*np.sin(phi) + Jy*np.cos(phi)
+
+ return np.sum(vectors*phihat, axis=-1)
+
+def get_sph_theta_component(vectors, coords, center, normal):
+ # The theta component of a vector is the vector dotted with thetahat
+
+ theta = get_sph_theta(coords, center, normal)
+ phi = get_sph_phi(coords, center, normal)
+
+ (xprime, yprime, zprime) = get_ortho_basis(normal)
+
+ tile_shape = list(vectors.shape)[:-1] + [1]
+ Jx = np.tile(xprime,tile_shape)
+ Jy = np.tile(yprime,tile_shape)
+ Jz = np.tile(zprime,tile_shape)
+
+ thetahat = Jx*np.cos(theta)*np.cos(phi) + \
+ Jy*np.cos(theta)*np.sin(theta) - \
+ Jz*np.sin(theta)
+
+ return np.sum(vectors*thetahat, axis=-1)
https://bitbucket.org/yt_analysis/yt/changeset/5d49ffb0208e/
changeset: 5d49ffb0208e
branch: yt
user: ngoldbaum
date: 2012-10-05 02:24:24
summary: Adding fields that grab the magnetic field in spherical coordinates.
affected #: 1 file
diff -r 8a9bc98c1f248a6f4080cfa35e79521a99b823e7 -r 5d49ffb0208ea4c4eafd593b38880d0dadd542d4 yt/data_objects/universal_fields.py
--- a/yt/data_objects/universal_fields.py
+++ b/yt/data_objects/universal_fields.py
@@ -267,8 +267,6 @@
add_field("sph_phi", function=_sph_phi,
validators=[ValidateParameter("center"),ValidateParameter("normal")])
-
-
### cylindrical coordinates: R (radius in the cylinder's plane)
def _cyl_R(field, data):
center = data.get_field_parameter("center")
@@ -1071,6 +1069,51 @@
display_name=r"\rm{Magnetic}\/\rm{Energy}",
units="\rm{ergs}\/\rm{cm}^{-3}")
+def _BPoloidal(field,data):
+ center = data.get_field_parameter("center")
+ normal = data.get_field_parameter("normal")
+
+ coords = np.array([data['x'] - center[0],
+ data['y'] - center[1],
+ data['z'] - center[2]]).transpose()
+
+ Bfields = np.array([data['Bx'], data['By'], data['Bz']])
+
+ return get_sph_theta_component(Bfields, coords, center, normal)
+add_field("BPoloidal", function=_BPoloidal,
+ units=r"\rm{Gauss}",
+ validators=[ValidateParameter("center"),ValidateParameter("normal")])
+
+def _BToroidal(field,data):
+ center = data.get_field_parameter("center")
+ normal = data.get_field_parameter("normal")
+
+ coords = np.array([data['x'] - center[0],
+ data['y'] - center[1],
+ data['z'] - center[2]]).transpose()
+
+ Bfields = np.array([data['Bx'], data['By'], data['Bz']])
+
+ return get_sph_phi_component(Bfields, coords, center, normal)
+add_field("BToroidal", function=_BToroidal,
+ units=r"\rm{Gauss}",
+ validators=[ValidateParameter("center"),ValidateParameter("normal")])
+
+def _BRadial(field,data):
+ center = data.get_field_parameter("center")
+ normal = data.get_field_parameter("normal")
+
+ coords = np.array([data['x'] - center[0],
+ data['y'] - center[1],
+ data['z'] - center[2]]).transpose()
+
+ Bfields = np.array([data['Bx'], data['By'], data['Bz']])
+
+ return get_sph_r_component(Bfields, coords, center, normal)
+add_field("BRadial", function=_BPoloidal,
+ units=r"\rm{Gauss}",
+ validators=[ValidateParameter("center"),ValidateParameter("normal")])
+
def _VorticitySquared(field, data):
mylog.debug("Generating vorticity on %s", data)
# We need to set up stencils
https://bitbucket.org/yt_analysis/yt/changeset/3bd4fe744fa7/
changeset: 3bd4fe744fa7
branch: yt
user: ngoldbaum
date: 2012-10-06 11:32:22
summary: Merging
affected #: 2 files
diff -r b8175a113d8ed98c0e8d514579ab607b01494565 -r 3bd4fe744fa7b1e36f6cbe86289405028ab2c149 yt/data_objects/universal_fields.py
--- a/yt/data_objects/universal_fields.py
+++ b/yt/data_objects/universal_fields.py
@@ -55,6 +55,14 @@
G, \
rho_crit_now, \
speed_of_light_cgs
+
+from yt.utilities.math_utils import \
+ get_sph_r_component, \
+ get_sph_theta_component, \
+ get_sph_phi_component, \
+ get_cyl_r_component, \
+ get_cyl_z_component, \
+ get_cyl_theta_component
# Note that, despite my newfound efforts to comply with PEP-8,
# I violate it here in order to keep the name/func_name relationship
@@ -194,13 +202,6 @@
function=_TangentialOverVelocityMagnitude,
take_log=False)
-def _TangentialVelocity(field, data):
- return np.sqrt(data["VelocityMagnitude"]**2.0
- - data["RadialVelocity"]**2.0)
-add_field("TangentialVelocity",
- function=_TangentialVelocity,
- take_log=False, units=r"\rm{cm}/\rm{s}")
-
def _Pressure(field, data):
"""M{(Gamma-1.0)*rho*E}"""
return (data.pf["Gamma"] - 1.0) * \
@@ -227,10 +228,7 @@
data['y'] - center[1],
data['z'] - center[2]]).transpose()
- ## The spherical coordinates radius is simply the magnitude of the
- ## coords vector.
-
- return np.sqrt(np.sum(coords**2,axis=-1))
+ return get_sph_r(vectors, center)
def _Convert_sph_r_CGS(data):
return data.convert("cm")
@@ -249,16 +247,7 @@
data['y'] - center[1],
data['z'] - center[2]]).transpose()
- ## The angle (theta) with respect to the normal (J), is the arccos
- ## of the dot product of the normal with the normalized coords
- ## vector.
-
- tile_shape = list(coords.shape)[:-1] + [1]
- J = np.tile(normal,tile_shape)
-
- JdotCoords = np.sum(J*coords,axis=-1)
-
- return np.arccos( JdotCoords / np.sqrt(np.sum(coords**2,axis=-1)) )
+ return get_sph_theta(coords, center, normal)
add_field("sph_theta", function=_sph_theta,
validators=[ValidateParameter("center"),ValidateParameter("normal")])
@@ -273,33 +262,11 @@
data['y'] - center[1],
data['z'] - center[2]]).transpose()
- ## We have freedom with respect to what axis (xprime) to define
- ## the disk angle. Here I've chosen to use the axis that is
- ## perpendicular to the normal and the y-axis. When normal ==
- ## y-hat, then set xprime = z-hat. With this definition, when
- ## normal == z-hat (as is typical), then xprime == x-hat.
- ##
- ## The angle is then given by the arctan of the ratio of the
- ## yprime-component and the xprime-component of the coords vector.
-
- xprime = np.cross([0.0,1.0,0.0],normal)
- if np.sum(xprime) == 0: xprime = np.array([0.0, 0.0, 1.0])
- yprime = np.cross(normal,xprime)
-
- tile_shape = list(coords.shape)[:-1] + [1]
- Jx = np.tile(xprime,tile_shape)
- Jy = np.tile(yprime,tile_shape)
-
- Px = np.sum(Jx*coords,axis=-1)
- Py = np.sum(Jy*coords,axis=-1)
-
- return np.arctan2(Py,Px)
+ return get_sph_phi(coords, center, normal)
add_field("sph_phi", function=_sph_phi,
validators=[ValidateParameter("center"),ValidateParameter("normal")])
-
-
### cylindrical coordinates: R (radius in the cylinder's plane)
def _cyl_R(field, data):
center = data.get_field_parameter("center")
@@ -309,14 +276,7 @@
data['y'] - center[1],
data['z'] - center[2]]).transpose()
- ## The cross product of the normal (J) with the coords vector
- ## gives a vector of magnitude equal to the cylindrical radius.
-
- tile_shape = list(coords.shape)[:-1] + [1]
- J = np.tile(normal,tile_shape)
-
- JcrossCoords = np.cross(J,coords)
- return np.sqrt(np.sum(JcrossCoords**2,axis=-1))
+ return get_cyl_r(coords, center, normal)
def _Convert_cyl_R_CGS(data):
return data.convert("cm")
@@ -324,6 +284,9 @@
add_field("cyl_R", function=_cyl_R,
validators=[ValidateParameter("center"),ValidateParameter("normal")],
convert_function = _Convert_cyl_R_CGS, units=r"\rm{cm}")
+add_field("cyl_RCode", function=_cyl_R,
+ validators=[ValidateParameter("center"),ValidateParameter("normal")],
+ units=r"Radius (code)")
### cylindrical coordinates: z (height above the cylinder's plane)
@@ -335,13 +298,7 @@
data['y'] - center[1],
data['z'] - center[2]]).transpose()
- ## The dot product of the normal (J) with the coords vector gives
- ## the cylindrical height.
-
- tile_shape = list(coords.shape)[:-1] + [1]
- J = np.tile(normal,tile_shape)
-
- return np.sum(J*coords,axis=-1)
+ return get_cyl_z(coords, center, normal)
def _Convert_cyl_z_CGS(data):
return data.convert("cm")
@@ -352,14 +309,19 @@
### cylindrical coordinates: theta (angle in the cylinder's plane)
-### [This is identical to the spherical coordinate's 'phi' angle.]
def _cyl_theta(field, data):
- return data['sph_phi']
+ center = data.get_field_parameter("center")
+ normal = data.get_field_parameter("normal")
+
+ coords = np.array([data['x'] - center[0],
+ data['y'] - center[1],
+ data['z'] - center[2]]).transpose()
+
+ return get_cyl_theta(coords, center, normal)
add_field("cyl_theta", function=_cyl_theta,
validators=[ValidateParameter("center"),ValidateParameter("normal")])
-
### The old field DiskAngle is the same as the spherical coordinates'
### 'theta' angle. I'm keeping DiskAngle for backwards compatibility.
def _DiskAngle(field, data):
@@ -392,6 +354,73 @@
ValidateParameter("normal")],
units=r"AU", display_field=False)
+def _cyl_RadialVelocity(field, data):
+ center = data.get_field_parameter("center")
+ normal = data.get_field_parameter("normal")
+ bulk_velocity = data.get_field_parameter("bulk_velocity")
+ if bulk_velocity == None:
+ bulk_velocity = np.zeros(3)
+
+ velocities = np.array([data['x-velocity'] - bulk_velocity[0],
+ data['y-velocity'] - bulk_velocity[1],
+ data['z-velocity'] - bulk_velocity[2]]).transpose()
+
+ coords = np.array([data['x'] - center[0],
+ data['y'] - center[1],
+ data['z'] - center[2]]).transpose()
+
+ return get_cyl_r_component(velocities, coords, center, normal)
+
+def _cyl_RadialVelocityABS(field, data):
+ return np.abs(_cyl_RadialVelocity(field, data))
+def _Convert_cyl_RadialVelocityKMS(data):
+ return 1e-5
+add_field("cyl_RadialVelocity", function=_cyl_RadialVelocity,
+ units=r"\rm{cm}/\rm{s}",
+ validators=[ValidateParameter("center"),ValidateParameter("normal")])
+add_field("cyl_RadialVelocityABS", function=_cyl_RadialVelocityABS,
+ units=r"\rm{cm}/\rm{s}",
+ validators=[ValidateParameter("center"),ValidateParameter("normal")])
+add_field("cyl_RadialVelocityKMS", function=_cyl_RadialVelocity,
+ convert_function=_Convert_cyl_RadialVelocityKMS, units=r"\rm{km}/\rm{s}",
+ validators=[ValidateParameter("center"),ValidateParameter("normal")])
+add_field("cyl_RadialVelocityKMSABS", function=_cyl_RadialVelocityABS,
+ convert_function=_Convert_cyl_RadialVelocityKMS, units=r"\rm{km}/\rm{s}",
+ validators=[ValidateParameter("center"),ValidateParameter("normal")])
+
+def _cyl_TangentialVelocity(field, data):
+ center = data.get_field_parameter("center")
+ normal = data.get_field_parameter("normal")
+ bulk_velocity = data.get_field_parameter("bulk_velocity")
+ if bulk_velocity == None:
+ bulk_velocity = np.zeros(3)
+
+ velocities = np.array([data['x-velocity'] - bulk_velocity[0],
+ data['y-velocity'] - bulk_velocity[1],
+ data['z-velocity'] - bulk_velocity[2]]).transpose()
+
+ coords = np.array([data['x'] - center[0],
+ data['y'] - center[1],
+ data['z'] - center[2]]).transpose()
+
+ return get_cyl_theta_component(velocities, coords, center, normal)
+
+def _cyl_TangentialVelocityABS(field, data):
+ return np.abs(_cyl_TangentialVelocity(field, data))
+def _Convert_cyl_TangentialVelocityKMS(data):
+ return 1e-5
+add_field("cyl_TangentialVelocity", function=_cyl_TangentialVelocity,
+ units=r"\rm{cm}/\rm{s}",
+ validators=[ValidateParameter("center"),ValidateParameter("normal")])
+add_field("cyl_TangentialVelocityABS", function=_cyl_TangentialVelocityABS,
+ units=r"\rm{cm}/\rm{s}",
+ validators=[ValidateParameter("center"),ValidateParameter("normal")])
+add_field("cyl_TangentialVelocityKMS", function=_cyl_TangentialVelocity,
+ convert_function=_Convert_cyl_TangentialVelocityKMS, units=r"\rm{km}/\rm{s}",
+ validators=[ValidateParameter("center"),ValidateParameter("normal")])
+add_field("cyl_TangentialVelocityKMSABS", function=_cyl_TangentialVelocityABS,
+ convert_function=_Convert_cyl_TangentialVelocityKMS, units=r"\rm{km}/\rm{s}",
+ validators=[ValidateParameter("center"),ValidateParameter("normal")])
def _DynamicalTime(field, data):
"""
@@ -884,16 +913,23 @@
def _RadialVelocity(field, data):
center = data.get_field_parameter("center")
+ normal = data.get_field_parameter("normal")
+ if normal == None:
+ normal = [0,0,1]
bulk_velocity = data.get_field_parameter("bulk_velocity")
if bulk_velocity == None:
bulk_velocity = np.zeros(3)
- new_field = ( (data['x']-center[0])*(data["x-velocity"]-bulk_velocity[0])
- + (data['y']-center[1])*(data["y-velocity"]-bulk_velocity[1])
- + (data['z']-center[2])*(data["z-velocity"]-bulk_velocity[2])
- )/data["RadiusCode"]
- if np.any(np.isnan(new_field)): # to fix center = point
- new_field[np.isnan(new_field)] = 0.0
- return new_field
+
+ velocities = np.array([data['x-velocity'] - bulk_velocity[0],
+ data['y-velocity'] - bulk_velocity[1],
+ data['z-velocity'] - bulk_velocity[2]]).transpose()
+
+ coords = np.array([data['x'] - center[0],
+ data['y'] - center[1],
+ data['z'] - center[2]]).transpose()
+
+ return get_sph_r_component(velocities, coords, center, normal)
+
def _RadialVelocityABS(field, data):
return np.abs(_RadialVelocity(field, data))
def _ConvertRadialVelocityKMS(data):
@@ -911,6 +947,13 @@
convert_function=_ConvertRadialVelocityKMS, units=r"\rm{km}/\rm{s}",
validators=[ValidateParameter("center")])
+def _TangentialVelocity(field, data):
+ return np.sqrt(data["VelocityMagnitude"]**2.0
+ - data["RadialVelocity"]**2.0)
+add_field("TangentialVelocity",
+ function=_TangentialVelocity,
+ take_log=False, units=r"\rm{cm}/\rm{s}")
+
def _CuttingPlaneVelocityX(field, data):
x_vec, y_vec, z_vec = [data.get_field_parameter("cp_%s_vec" % (ax))
for ax in 'xyz']
@@ -1026,6 +1069,51 @@
display_name=r"\rm{Magnetic}\/\rm{Energy}",
units="\rm{ergs}\/\rm{cm}^{-3}")
+def _BPoloidal(field,data):
+ center = data.get_field_parameter("center")
+ normal = data.get_field_parameter("normal")
+
+ coords = np.array([data['x'] - center[0],
+ data['y'] - center[1],
+ data['z'] - center[2]]).transpose()
+
+ Bfields = np.array([data['Bx'], data['By'], data['Bz']])
+
+ return get_sph_theta_component(Bfields, coords, center, normal)
+add_field("BPoloidal", function=_BPoloidal,
+ units=r"\rm{Gauss}",
+ validators=[ValidateParameter("center"),ValidateParameter("normal")])
+
+def _BToroidal(field,data):
+ center = data.get_field_parameter("center")
+ normal = data.get_field_parameter("normal")
+
+ coords = np.array([data['x'] - center[0],
+ data['y'] - center[1],
+ data['z'] - center[2]]).transpose()
+
+ Bfields = np.array([data['Bx'], data['By'], data['Bz']])
+
+ return get_sph_phi_component(Bfields, coords, center, normal)
+add_field("BToroidal", function=_BToroidal,
+ units=r"\rm{Gauss}",
+ validators=[ValidateParameter("center"),ValidateParameter("normal")])
+
+def _BRadial(field,data):
+ center = data.get_field_parameter("center")
+ normal = data.get_field_parameter("normal")
+
+ coords = np.array([data['x'] - center[0],
+ data['y'] - center[1],
+ data['z'] - center[2]]).transpose()
+
+ Bfields = np.array([data['Bx'], data['By'], data['Bz']])
+
+ return get_sph_r_component(Bfields, coords, center, normal)
+add_field("BRadial", function=_BPoloidal,
+ units=r"\rm{Gauss}",
+ validators=[ValidateParameter("center"),ValidateParameter("normal")])
+
def _VorticitySquared(field, data):
mylog.debug("Generating vorticity on %s", data)
# We need to set up stencils
diff -r b8175a113d8ed98c0e8d514579ab607b01494565 -r 3bd4fe744fa7b1e36f6cbe86289405028ab2c149 yt/utilities/math_utils.py
--- a/yt/utilities/math_utils.py
+++ b/yt/utilities/math_utils.py
@@ -674,3 +674,169 @@
[uz*ux*(1-cost)-uy*sint, uz*uy*(1-cost)+ux*sint, cost+uz**2*(1-cost)]])
return R
+
+def get_ortho_basis(normal):
+ xprime = np.cross([0.0,1.0,0.0],normal)
+ if np.sum(xprime) == 0: xprime = np.array([0.0, 0.0, 1.0])
+ yprime = np.cross(normal,xprime)
+ zprime = normal
+ return (xprime, yprime, zprime)
+
+def get_sph_r(coords, center):
+ # The spherical coordinates radius is simply the magnitude of the
+ # coordinate vector.
+
+ return np.sqrt(np.sum(coords**2, axis=-1))
+
+
+def get_sph_theta(coords, center, normal):
+ # The angle (theta) with respect to the normal (J), is the arccos
+ # of the dot product of the normal with the normalized coordinate
+ # vector.
+
+ tile_shape = list(coords.shape)[:-1] + [1]
+ J = np.tile(normal,tile_shape)
+
+ JdotCoords = np.sum(J*coords,axis=-1)
+
+ return np.arccos( JdotCoords / np.sqrt(np.sum(coords**2,axis=-1)) )
+
+def get_sph_phi(coords, center, normal):
+ # We have freedom with respect to what axis (xprime) to define
+ # the disk angle. Here I've chosen to use the axis that is
+ # perpendicular to the normal and the y-axis. When normal ==
+ # y-hat, then set xprime = z-hat. With this definition, when
+ # normal == z-hat (as is typical), then xprime == x-hat.
+ #
+ # The angle is then given by the arctan of the ratio of the
+ # yprime-component and the xprime-component of the coordinate
+ # vector.
+
+ (xprime, yprime, zprime) = get_ortho_basis(normal)
+
+ tile_shape = list(coords.shape)[:-1] + [1]
+ Jx = np.tile(xprime,tile_shape)
+ Jy = np.tile(yprime,tile_shape)
+
+ Px = np.sum(Jx*coords,axis=-1)
+ Py = np.sum(Jy*coords,axis=-1)
+
+ return np.arctan2(Py,Px)
+
+def get_cyl_r(coords, center, normal):
+ # The cross product of the normal (J) with a coordinate vector
+ # gives a vector of magnitude equal to the cylindrical radius.
+
+ tile_shape = list(coords.shape)[:-1] + [1]
+ J = np.tile(normal, tile_shape)
+
+ JcrossCoords = np.cross(J, coords)
+ return np.sqrt(np.sum(JcrossCoords**2, axis=-1))
+
+def get_cyl_z(coords, center, normal):
+ # The dot product of the normal (J) with the coordinate vector
+ # gives the cylindrical height.
+
+ tile_shape = list(coords.shape)[:-1] + [1]
+ J = np.tile(normal, tile_shape)
+
+ return np.sum(J*coords, axis=-1)
+
+def get_cyl_theta(coords, center, normal):
+ # This is identical to the spherical phi component
+
+ return get_sph_phi(coords, center, normal)
+
+
+def get_cyl_r_component(vectors, coords, center, normal):
+ # The r of a vector is the vector dotted with rhat
+
+ theta = np.tile(get_cyl_theta(coords, center, normal), (3,1)).transpose()
+
+ (xprime, yprime, zprime) = get_ortho_basis(normal)
+
+ tile_shape = list(vectors.shape)[:-1] + [1]
+ Jx = np.tile(xprime,tile_shape)
+ Jy = np.tile(yprime,tile_shape)
+
+ rhat = Jx*np.cos(theta) + Jy*np.sin(theta)
+
+ return np.sum(vectors*rhat,axis=-1)
+
+def get_cyl_theta_component(vectors, coords, center, normal):
+ # The theta component of a vector is the vector dotted with thetahat
+
+ theta = np.tile(get_cyl_theta(coords, center, normal), (3,1)).transpose()
+
+ (xprime, yprime, zprime) = get_ortho_basis(normal)
+
+ tile_shape = list(vectors.shape)[:-1] + [1]
+ Jx = np.tile(xprime,tile_shape)
+ Jy = np.tile(yprime,tile_shape)
+
+ thetahat = -Jx*np.sin(theta) + Jy*np.cos(theta)
+
+ return np.sum(vectors*thetahat, axis=-1)
+
+def get_cyl_z_component(vectors, center, normal):
+ # The z component of a vector is the vector dotted with zhat
+
+ (xprime, yprime, zprime) = get_ortho_basis(normal)
+
+ tile_shape = list(vectors.shape)[:-1] + [1]
+ zhat = np.tile(zprime, tile_shape)
+
+ return np.sum(vectors*zhat, axis=-1)
+
+def get_sph_r_component(vectors, coords, center, normal):
+ # The r component of a vector is the vector dotted with rhat
+
+ theta = get_sph_theta(coords, center, normal)
+ phi = get_sph_phi(coords, center, normal)
+
+ (xprime, yprime, zprime) = get_ortho_basis(normal)
+
+ tile_shape = list(vectors.shape)[:-1] + [1]
+ Jx = np.tile(xprime,tile_shape)
+ Jy = np.tile(yprime,tile_shape)
+ Jz = np.tile(zprime,tile_shape)
+
+ rhat = Jx*np.sin(theta)*np.cos(phi) + \
+ Jy*np.sin(theta)*np.sin(phi) + \
+ Jz*np.cos(theta)
+
+ return np.sum(vectors*rhat, axis=-1)
+
+def get_sph_phi_component(vectors, coords, center, normal):
+ # The phi component of a vector is the vector dotted with phihat
+
+ phi = get_sph_phi(coords, center, normal)
+
+ (xprime, yprime, zprime) = get_ortho_basis(normal)
+
+ tile_shape = list(vectors.shape)[:-1] + [1]
+ Jx = np.tile(xprime,tile_shape)
+ Jy = np.tile(yprime,tile_shape)
+
+ phihat = -Jx*np.sin(phi) + Jy*np.cos(phi)
+
+ return np.sum(vectors*phihat, axis=-1)
+
+def get_sph_theta_component(vectors, coords, center, normal):
+ # The theta component of a vector is the vector dotted with thetahat
+
+ theta = get_sph_theta(coords, center, normal)
+ phi = get_sph_phi(coords, center, normal)
+
+ (xprime, yprime, zprime) = get_ortho_basis(normal)
+
+ tile_shape = list(vectors.shape)[:-1] + [1]
+ Jx = np.tile(xprime,tile_shape)
+ Jy = np.tile(yprime,tile_shape)
+ Jz = np.tile(zprime,tile_shape)
+
+ thetahat = Jx*np.cos(theta)*np.cos(phi) + \
+ Jy*np.cos(theta)*np.sin(theta) - \
+ Jz*np.sin(theta)
+
+ return np.sum(vectors*thetahat, axis=-1)
https://bitbucket.org/yt_analysis/yt/changeset/809bff333e8b/
changeset: 809bff333e8b
branch: yt
user: ngoldbaum
date: 2012-10-06 12:20:53
summary: Rearranging things. Fixing the imports so that the cylindrical and spherical coordinate fields don't break. Fixing some typos.
This still reads x, y, z, vx, vy, and vz multiple times for my test script.
affected #: 2 files
diff -r 3bd4fe744fa7b1e36f6cbe86289405028ab2c149 -r 809bff333e8b5022ed54799209c3e69a668a78f1 yt/data_objects/universal_fields.py
--- a/yt/data_objects/universal_fields.py
+++ b/yt/data_objects/universal_fields.py
@@ -62,7 +62,10 @@
get_sph_phi_component, \
get_cyl_r_component, \
get_cyl_z_component, \
- get_cyl_theta_component
+ get_cyl_theta_component, \
+ get_cyl_r, get_cyl_theta, \
+ get_cyl_z, get_sph_r, \
+ get_sph_theta, get_sph_phi
# Note that, despite my newfound efforts to comply with PEP-8,
# I violate it here in order to keep the name/func_name relationship
@@ -247,7 +250,7 @@
data['y'] - center[1],
data['z'] - center[2]]).transpose()
- return get_sph_theta(coords, center, normal)
+ return get_sph_theta(coords, normal)
add_field("sph_theta", function=_sph_theta,
validators=[ValidateParameter("center"),ValidateParameter("normal")])
@@ -262,7 +265,7 @@
data['y'] - center[1],
data['z'] - center[2]]).transpose()
- return get_sph_phi(coords, center, normal)
+ return get_sph_phi(coords, normal)
add_field("sph_phi", function=_sph_phi,
validators=[ValidateParameter("center"),ValidateParameter("normal")])
@@ -276,7 +279,7 @@
data['y'] - center[1],
data['z'] - center[2]]).transpose()
- return get_cyl_r(coords, center, normal)
+ return get_cyl_r(coords, normal)
def _Convert_cyl_R_CGS(data):
return data.convert("cm")
@@ -298,7 +301,7 @@
data['y'] - center[1],
data['z'] - center[2]]).transpose()
- return get_cyl_z(coords, center, normal)
+ return get_cyl_z(coords, normal)
def _Convert_cyl_z_CGS(data):
return data.convert("cm")
@@ -317,7 +320,7 @@
data['y'] - center[1],
data['z'] - center[2]]).transpose()
- return get_cyl_theta(coords, center, normal)
+ return get_cyl_theta(coords, normal)
add_field("cyl_theta", function=_cyl_theta,
validators=[ValidateParameter("center"),ValidateParameter("normal")])
@@ -355,7 +358,6 @@
units=r"AU", display_field=False)
def _cyl_RadialVelocity(field, data):
- center = data.get_field_parameter("center")
normal = data.get_field_parameter("normal")
bulk_velocity = data.get_field_parameter("bulk_velocity")
if bulk_velocity == None:
@@ -365,11 +367,9 @@
data['y-velocity'] - bulk_velocity[1],
data['z-velocity'] - bulk_velocity[2]]).transpose()
- coords = np.array([data['x'] - center[0],
- data['y'] - center[1],
- data['z'] - center[2]]).transpose()
+ theta = np.tile(data['cyl_theta'], (3, 1)).transpose()
- return get_cyl_r_component(velocities, coords, center, normal)
+ return get_cyl_r_component(velocities, theta, normal)
def _cyl_RadialVelocityABS(field, data):
return np.abs(_cyl_RadialVelocity(field, data))
@@ -377,19 +377,18 @@
return 1e-5
add_field("cyl_RadialVelocity", function=_cyl_RadialVelocity,
units=r"\rm{cm}/\rm{s}",
- validators=[ValidateParameter("center"),ValidateParameter("normal")])
+ validators=[ValidateParameter("normal")])
add_field("cyl_RadialVelocityABS", function=_cyl_RadialVelocityABS,
units=r"\rm{cm}/\rm{s}",
- validators=[ValidateParameter("center"),ValidateParameter("normal")])
+ validators=[ValidateParameter("normal")])
add_field("cyl_RadialVelocityKMS", function=_cyl_RadialVelocity,
convert_function=_Convert_cyl_RadialVelocityKMS, units=r"\rm{km}/\rm{s}",
- validators=[ValidateParameter("center"),ValidateParameter("normal")])
+ validators=[ValidateParameter("normal")])
add_field("cyl_RadialVelocityKMSABS", function=_cyl_RadialVelocityABS,
convert_function=_Convert_cyl_RadialVelocityKMS, units=r"\rm{km}/\rm{s}",
- validators=[ValidateParameter("center"),ValidateParameter("normal")])
+ validators=[ValidateParameter("normal")])
def _cyl_TangentialVelocity(field, data):
- center = data.get_field_parameter("center")
normal = data.get_field_parameter("normal")
bulk_velocity = data.get_field_parameter("bulk_velocity")
if bulk_velocity == None:
@@ -399,11 +398,9 @@
data['y-velocity'] - bulk_velocity[1],
data['z-velocity'] - bulk_velocity[2]]).transpose()
- coords = np.array([data['x'] - center[0],
- data['y'] - center[1],
- data['z'] - center[2]]).transpose()
-
- return get_cyl_theta_component(velocities, coords, center, normal)
+ theta = np.tile(data['cyl_theta'], (3, 1)).transpose()
+
+ return get_cyl_theta_component(velocities, theta, normal)
def _cyl_TangentialVelocityABS(field, data):
return np.abs(_cyl_TangentialVelocity(field, data))
@@ -411,16 +408,16 @@
return 1e-5
add_field("cyl_TangentialVelocity", function=_cyl_TangentialVelocity,
units=r"\rm{cm}/\rm{s}",
- validators=[ValidateParameter("center"),ValidateParameter("normal")])
+ validators=[ValidateParameter("normal")])
add_field("cyl_TangentialVelocityABS", function=_cyl_TangentialVelocityABS,
units=r"\rm{cm}/\rm{s}",
- validators=[ValidateParameter("center"),ValidateParameter("normal")])
+ validators=[ValidateParameter("normal")])
add_field("cyl_TangentialVelocityKMS", function=_cyl_TangentialVelocity,
convert_function=_Convert_cyl_TangentialVelocityKMS, units=r"\rm{km}/\rm{s}",
- validators=[ValidateParameter("center"),ValidateParameter("normal")])
+ validators=[ValidateParameter("normal")])
add_field("cyl_TangentialVelocityKMSABS", function=_cyl_TangentialVelocityABS,
convert_function=_Convert_cyl_TangentialVelocityKMS, units=r"\rm{km}/\rm{s}",
- validators=[ValidateParameter("center"),ValidateParameter("normal")])
+ validators=[ValidateParameter("normal")])
def _DynamicalTime(field, data):
"""
@@ -912,7 +909,6 @@
display_name = "Radius (code)")
def _RadialVelocity(field, data):
- center = data.get_field_parameter("center")
normal = data.get_field_parameter("normal")
if normal == None:
normal = [0,0,1]
@@ -924,28 +920,23 @@
data['y-velocity'] - bulk_velocity[1],
data['z-velocity'] - bulk_velocity[2]]).transpose()
- coords = np.array([data['x'] - center[0],
- data['y'] - center[1],
- data['z'] - center[2]]).transpose()
+ theta = np.tile(data['sph_theta'], (3, 1)).transpose()
+ phi = np.tile(data['sph_phi'], (3, 1)).transpose()
- return get_sph_r_component(velocities, coords, center, normal)
+ return get_sph_r_component(velocities, theta, phi, normal)
def _RadialVelocityABS(field, data):
return np.abs(_RadialVelocity(field, data))
def _ConvertRadialVelocityKMS(data):
return 1e-5
add_field("RadialVelocity", function=_RadialVelocity,
- units=r"\rm{cm}/\rm{s}",
- validators=[ValidateParameter("center")])
+ units=r"\rm{cm}/\rm{s}")
add_field("RadialVelocityABS", function=_RadialVelocityABS,
- units=r"\rm{cm}/\rm{s}",
- validators=[ValidateParameter("center")])
+ units=r"\rm{cm}/\rm{s}")
add_field("RadialVelocityKMS", function=_RadialVelocity,
- convert_function=_ConvertRadialVelocityKMS, units=r"\rm{km}/\rm{s}",
- validators=[ValidateParameter("center")])
+ convert_function=_ConvertRadialVelocityKMS, units=r"\rm{km}/\rm{s}")
add_field("RadialVelocityKMSABS", function=_RadialVelocityABS,
- convert_function=_ConvertRadialVelocityKMS, units=r"\rm{km}/\rm{s}",
- validators=[ValidateParameter("center")])
+ convert_function=_ConvertRadialVelocityKMS, units=r"\rm{km}/\rm{s}")
def _TangentialVelocity(field, data):
return np.sqrt(data["VelocityMagnitude"]**2.0
@@ -1070,49 +1061,45 @@
units="\rm{ergs}\/\rm{cm}^{-3}")
def _BPoloidal(field,data):
- center = data.get_field_parameter("center")
normal = data.get_field_parameter("normal")
- coords = np.array([data['x'] - center[0],
- data['y'] - center[1],
- data['z'] - center[2]]).transpose()
-
Bfields = np.array([data['Bx'], data['By'], data['Bz']])
- return get_sph_theta_component(Bfields, coords, center, normal)
+ theta = np.tile(data['sph_theta'], (3, 1)).transpose()
+ phi = np.tile(data['sph_phi'], (3, 1)).transpose()
+
+ return get_sph_theta_component(Bfields, theta, phi, normal)
+
add_field("BPoloidal", function=_BPoloidal,
units=r"\rm{Gauss}",
- validators=[ValidateParameter("center"),ValidateParameter("normal")])
+ validators=[ValidateParameter("normal")])
def _BToroidal(field,data):
- center = data.get_field_parameter("center")
normal = data.get_field_parameter("normal")
- coords = np.array([data['x'] - center[0],
- data['y'] - center[1],
- data['z'] - center[2]]).transpose()
-
Bfields = np.array([data['Bx'], data['By'], data['Bz']])
- return get_sph_phi_component(Bfields, coords, center, normal)
+ phi = np.tile(data['sph_phi'], (3, 1)).transpose()
+
+ return get_sph_phi_component(Bfields, phi, normal)
+
add_field("BToroidal", function=_BToroidal,
units=r"\rm{Gauss}",
- validators=[ValidateParameter("center"),ValidateParameter("normal")])
+ validators=[ValidateParameter("normal")])
def _BRadial(field,data):
- center = data.get_field_parameter("center")
normal = data.get_field_parameter("normal")
- coords = np.array([data['x'] - center[0],
- data['y'] - center[1],
- data['z'] - center[2]]).transpose()
-
Bfields = np.array([data['Bx'], data['By'], data['Bz']])
- return get_sph_r_component(Bfields, coords, center, normal)
+ theta = np.tile(data['sph_theta'], (3, 1)).transpose()
+ phi = np.tile(data['sph_phi'], (3, 1)).transpose()
+
+ return get_sph_r_component(Bfields, theta, phi, normal)
+
add_field("BRadial", function=_BPoloidal,
units=r"\rm{Gauss}",
- validators=[ValidateParameter("center"),ValidateParameter("normal")])
+ validators=[ValidateParameter("normal")])
def _VorticitySquared(field, data):
mylog.debug("Generating vorticity on %s", data)
diff -r 3bd4fe744fa7b1e36f6cbe86289405028ab2c149 -r 809bff333e8b5022ed54799209c3e69a668a78f1 yt/utilities/math_utils.py
--- a/yt/utilities/math_utils.py
+++ b/yt/utilities/math_utils.py
@@ -682,14 +682,14 @@
zprime = normal
return (xprime, yprime, zprime)
-def get_sph_r(coords, center):
+def get_sph_r(coords):
# The spherical coordinates radius is simply the magnitude of the
# coordinate vector.
return np.sqrt(np.sum(coords**2, axis=-1))
-def get_sph_theta(coords, center, normal):
+def get_sph_theta(coords, normal):
# The angle (theta) with respect to the normal (J), is the arccos
# of the dot product of the normal with the normalized coordinate
# vector.
@@ -701,7 +701,7 @@
return np.arccos( JdotCoords / np.sqrt(np.sum(coords**2,axis=-1)) )
-def get_sph_phi(coords, center, normal):
+def get_sph_phi(coords, normal):
# We have freedom with respect to what axis (xprime) to define
# the disk angle. Here I've chosen to use the axis that is
# perpendicular to the normal and the y-axis. When normal ==
@@ -723,7 +723,7 @@
return np.arctan2(Py,Px)
-def get_cyl_r(coords, center, normal):
+def get_cyl_r(coords, normal):
# The cross product of the normal (J) with a coordinate vector
# gives a vector of magnitude equal to the cylindrical radius.
@@ -733,7 +733,7 @@
JcrossCoords = np.cross(J, coords)
return np.sqrt(np.sum(JcrossCoords**2, axis=-1))
-def get_cyl_z(coords, center, normal):
+def get_cyl_z(coords, normal):
# The dot product of the normal (J) with the coordinate vector
# gives the cylindrical height.
@@ -742,17 +742,15 @@
return np.sum(J*coords, axis=-1)
-def get_cyl_theta(coords, center, normal):
+def get_cyl_theta(coords, normal):
# This is identical to the spherical phi component
- return get_sph_phi(coords, center, normal)
+ return get_sph_phi(coords, normal)
-def get_cyl_r_component(vectors, coords, center, normal):
+def get_cyl_r_component(vectors, theta, normal):
# The r of a vector is the vector dotted with rhat
- theta = np.tile(get_cyl_theta(coords, center, normal), (3,1)).transpose()
-
(xprime, yprime, zprime) = get_ortho_basis(normal)
tile_shape = list(vectors.shape)[:-1] + [1]
@@ -763,11 +761,9 @@
return np.sum(vectors*rhat,axis=-1)
-def get_cyl_theta_component(vectors, coords, center, normal):
+def get_cyl_theta_component(vectors, theta, normal):
# The theta component of a vector is the vector dotted with thetahat
- theta = np.tile(get_cyl_theta(coords, center, normal), (3,1)).transpose()
-
(xprime, yprime, zprime) = get_ortho_basis(normal)
tile_shape = list(vectors.shape)[:-1] + [1]
@@ -778,9 +774,8 @@
return np.sum(vectors*thetahat, axis=-1)
-def get_cyl_z_component(vectors, center, normal):
+def get_cyl_z_component(vectors, normal):
# The z component of a vector is the vector dotted with zhat
-
(xprime, yprime, zprime) = get_ortho_basis(normal)
tile_shape = list(vectors.shape)[:-1] + [1]
@@ -788,12 +783,9 @@
return np.sum(vectors*zhat, axis=-1)
-def get_sph_r_component(vectors, coords, center, normal):
+def get_sph_r_component(vectors, theta, phi, normal):
# The r component of a vector is the vector dotted with rhat
- theta = get_sph_theta(coords, center, normal)
- phi = get_sph_phi(coords, center, normal)
-
(xprime, yprime, zprime) = get_ortho_basis(normal)
tile_shape = list(vectors.shape)[:-1] + [1]
@@ -807,10 +799,10 @@
return np.sum(vectors*rhat, axis=-1)
-def get_sph_phi_component(vectors, coords, center, normal):
+def get_sph_phi_component(vectors, phi, normal):
# The phi component of a vector is the vector dotted with phihat
- phi = get_sph_phi(coords, center, normal)
+ phi = get_sph_phi(coords, normal)
(xprime, yprime, zprime) = get_ortho_basis(normal)
@@ -822,11 +814,11 @@
return np.sum(vectors*phihat, axis=-1)
-def get_sph_theta_component(vectors, coords, center, normal):
+def get_sph_theta_component(vectors, theta, phi, normal):
# The theta component of a vector is the vector dotted with thetahat
- theta = get_sph_theta(coords, center, normal)
- phi = get_sph_phi(coords, center, normal)
+ theta = get_sph_theta(coords, normal)
+ phi = get_sph_phi(coords, normal)
(xprime, yprime, zprime) = get_ortho_basis(normal)
https://bitbucket.org/yt_analysis/yt/changeset/28ee8407e8ac/
changeset: 28ee8407e8ac
branch: yt
user: ngoldbaum
date: 2012-10-06 21:42:52
summary: Moving the code that obtains relative velocities to cython. Refactoring the coordinate transformation code. Using the appropriate unit converstion constants from physical_constants.py. Thanks for the suggestions Kacper!
affected #: 3 files
diff -r 809bff333e8b5022ed54799209c3e69a668a78f1 -r 28ee8407e8ac563031468f9a1f535e52b3f4eaae yt/data_objects/universal_fields.py
--- a/yt/data_objects/universal_fields.py
+++ b/yt/data_objects/universal_fields.py
@@ -32,7 +32,7 @@
from yt.funcs import *
-from yt.utilities.lib import CICDeposit_3, obtain_rvec
+from yt.utilities.lib import CICDeposit_3, obtain_rvec, obtain_rv_vec
from yt.utilities.cosmology import Cosmology
from field_info_container import \
add_field, \
@@ -54,7 +54,8 @@
kboltz, \
G, \
rho_crit_now, \
- speed_of_light_cgs
+ speed_of_light_cgs, \
+ km_per_cm
from yt.utilities.math_utils import \
get_sph_r_component, \
@@ -190,12 +191,8 @@
def _VelocityMagnitude(field, data):
"""M{|v|}"""
- bulk_velocity = data.get_field_parameter("bulk_velocity")
- if bulk_velocity == None:
- bulk_velocity = np.zeros(3)
- return ( (data["x-velocity"]-bulk_velocity[0])**2.0 + \
- (data["y-velocity"]-bulk_velocity[1])**2.0 + \
- (data["z-velocity"]-bulk_velocity[2])**2.0 )**(1.0/2.0)
+ velocities = obtain_rv_vec(data)
+ return np.sqrt(np.sum(velocities**2,axis=-1))
add_field("VelocityMagnitude", function=_VelocityMagnitude,
take_log=False, units=r"\rm{cm}/\rm{s}")
@@ -227,9 +224,7 @@
def _sph_r(field, data):
center = data.get_field_parameter("center")
- coords = np.array([data['x'] - center[0],
- data['y'] - center[1],
- data['z'] - center[2]]).transpose()
+ coords = obtain_rvec(data).transpose()
return get_sph_r(vectors, center)
@@ -246,9 +241,7 @@
center = data.get_field_parameter("center")
normal = data.get_field_parameter("normal")
- coords = np.array([data['x'] - center[0],
- data['y'] - center[1],
- data['z'] - center[2]]).transpose()
+ coords = obtain_rvec(data).transepose()
return get_sph_theta(coords, normal)
@@ -261,10 +254,8 @@
center = data.get_field_parameter("center")
normal = data.get_field_parameter("normal")
- coords = np.array([data['x'] - center[0],
- data['y'] - center[1],
- data['z'] - center[2]]).transpose()
-
+ coords = obtain_rvec(data).transpose()
+
return get_sph_phi(coords, normal)
add_field("sph_phi", function=_sph_phi,
@@ -275,9 +266,7 @@
center = data.get_field_parameter("center")
normal = data.get_field_parameter("normal")
- coords = np.array([data['x'] - center[0],
- data['y'] - center[1],
- data['z'] - center[2]]).transpose()
+ coords = obtain_rvec(data).transpose()
return get_cyl_r(coords, normal)
@@ -297,9 +286,7 @@
center = data.get_field_parameter("center")
normal = data.get_field_parameter("normal")
- coords = np.array([data['x'] - center[0],
- data['y'] - center[1],
- data['z'] - center[2]]).transpose()
+ coords = obtain_rvec(data).transpose()
return get_cyl_z(coords, normal)
@@ -316,9 +303,7 @@
center = data.get_field_parameter("center")
normal = data.get_field_parameter("normal")
- coords = np.array([data['x'] - center[0],
- data['y'] - center[1],
- data['z'] - center[2]]).transpose()
+ coords = obtain_rvec(data).transpose()
return get_cyl_theta(coords, normal)
@@ -359,13 +344,7 @@
def _cyl_RadialVelocity(field, data):
normal = data.get_field_parameter("normal")
- bulk_velocity = data.get_field_parameter("bulk_velocity")
- if bulk_velocity == None:
- bulk_velocity = np.zeros(3)
-
- velocities = np.array([data['x-velocity'] - bulk_velocity[0],
- data['y-velocity'] - bulk_velocity[1],
- data['z-velocity'] - bulk_velocity[2]]).transpose()
+ velocities = obtain_rv_vec(data).transpose()
theta = np.tile(data['cyl_theta'], (3, 1)).transpose()
@@ -374,7 +353,7 @@
def _cyl_RadialVelocityABS(field, data):
return np.abs(_cyl_RadialVelocity(field, data))
def _Convert_cyl_RadialVelocityKMS(data):
- return 1e-5
+ return km_per_cm
add_field("cyl_RadialVelocity", function=_cyl_RadialVelocity,
units=r"\rm{cm}/\rm{s}",
validators=[ValidateParameter("normal")])
@@ -390,14 +369,7 @@
def _cyl_TangentialVelocity(field, data):
normal = data.get_field_parameter("normal")
- bulk_velocity = data.get_field_parameter("bulk_velocity")
- if bulk_velocity == None:
- bulk_velocity = np.zeros(3)
-
- velocities = np.array([data['x-velocity'] - bulk_velocity[0],
- data['y-velocity'] - bulk_velocity[1],
- data['z-velocity'] - bulk_velocity[2]]).transpose()
-
+ velocities = obtain_rv_vec(data).transpose()
theta = np.tile(data['cyl_theta'], (3, 1)).transpose()
return get_cyl_theta_component(velocities, theta, normal)
@@ -405,7 +377,7 @@
def _cyl_TangentialVelocityABS(field, data):
return np.abs(_cyl_TangentialVelocity(field, data))
def _Convert_cyl_TangentialVelocityKMS(data):
- return 1e-5
+ return km_per_cm
add_field("cyl_TangentialVelocity", function=_cyl_TangentialVelocity,
units=r"\rm{cm}/\rm{s}",
validators=[ValidateParameter("normal")])
@@ -666,13 +638,7 @@
take_log=False, display_field=False)
def obtain_velocities(data):
- if data.has_field_parameter("bulk_velocity"):
- bv = data.get_field_parameter("bulk_velocity")
- else: bv = np.zeros(3, dtype='float64')
- xv = data["x-velocity"] - bv[0]
- yv = data["y-velocity"] - bv[1]
- zv = data["z-velocity"] - bv[2]
- return xv, yv, zv
+ return obtain_rv_vec(data)
def _convertSpecificAngularMomentum(data):
return data.convert("cm")
@@ -737,7 +703,7 @@
# convert_function=_convertSpecificAngularMomentum, vector_field=True,
# units=r"\rm{cm}^2/\rm{s}", validators=[ValidateParameter('center')])
def _convertSpecificAngularMomentumKMSMPC(data):
- return data.convert("mpc")/1e5
+ return km_per_cm*data.convert("mpc")
#add_field("ParticleSpecificAngularMomentumKMSMPC",
# function=_ParticleSpecificAngularMomentum, particle_type=True,
# convert_function=_convertSpecificAngularMomentumKMSMPC, vector_field=True,
@@ -912,14 +878,7 @@
normal = data.get_field_parameter("normal")
if normal == None:
normal = [0,0,1]
- bulk_velocity = data.get_field_parameter("bulk_velocity")
- if bulk_velocity == None:
- bulk_velocity = np.zeros(3)
-
- velocities = np.array([data['x-velocity'] - bulk_velocity[0],
- data['y-velocity'] - bulk_velocity[1],
- data['z-velocity'] - bulk_velocity[2]]).transpose()
-
+ velocities = obtain_rv_vec(data).transpose()
theta = np.tile(data['sph_theta'], (3, 1)).transpose()
phi = np.tile(data['sph_phi'], (3, 1)).transpose()
@@ -928,7 +887,7 @@
def _RadialVelocityABS(field, data):
return np.abs(_RadialVelocity(field, data))
def _ConvertRadialVelocityKMS(data):
- return 1e-5
+ return km_per_cm
add_field("RadialVelocity", function=_RadialVelocity,
units=r"\rm{cm}/\rm{s}")
add_field("RadialVelocityABS", function=_RadialVelocityABS,
diff -r 809bff333e8b5022ed54799209c3e69a668a78f1 -r 28ee8407e8ac563031468f9a1f535e52b3f4eaae yt/utilities/lib/geometry_utils.pyx
--- a/yt/utilities/lib/geometry_utils.pyx
+++ b/yt/utilities/lib/geometry_utils.pyx
@@ -338,3 +338,45 @@
rg[2,i,j,k] = zg[i,j,k] - c[2]
return rg
+ at cython.boundscheck(False)
+ at cython.wraparound(False)
+ at cython.cdivision(True)
+def obtain_rv_vec(data):
+ # This is just to let the pointers exist and whatnot. We can't cdef them
+ # inside conditionals.
+ cdef np.ndarray[np.float64_t, ndim=1] vxf
+ cdef np.ndarray[np.float64_t, ndim=1] vyf
+ cdef np.ndarray[np.float64_t, ndim=1] vzf
+ cdef np.ndarray[np.float64_t, ndim=2] rvf
+ cdef np.ndarray[np.float64_t, ndim=3] vxg
+ cdef np.ndarray[np.float64_t, ndim=3] vyg
+ cdef np.ndarray[np.float64_t, ndim=3] vzg
+ cdef np.ndarray[np.float64_t, ndim=4] rvg
+ cdef np.float64_t bv[3]
+ cdef int i, j, k
+ bulk_velocity = data.get_field_parameter("bulk_velocity")
+ bv[0] = bulk_velocity[0]; bv[1] = bulk_velocity[1]; bv[2] = bulk_velocity[2]
+ if len(data['x-velocity'].shape) == 1:
+ # One dimensional data
+ vxf = data['x-velocity']
+ vyf = data['y-velocity']
+ vzf = data['z-velocity']
+ rvf = np.empty((3, vxf.shape[0]), 'float64')
+ for i in range(vxf.shape[0]):
+ rvf[0, i] = vxf[i] - bv[0]
+ rvf[1, i] = vyf[i] - bv[1]
+ rvf[2, i] = vzf[i] - bv[2]
+ return rvf
+ else:
+ # Three dimensional data
+ vxg = data['x']
+ vyg = data['y']
+ vzg = data['z']
+ rvg = np.empty((3, vxg.shape[0], vxg.shape[1], vxg.shape[2]), 'float64')
+ for i in range(vxg.shape[0]):
+ for j in range(vxg.shape[1]):
+ for k in range(vxg.shape[2]):
+ rvg[0,i,j,k] = vxg[i,j,k] - bv[0]
+ rvg[1,i,j,k] = vyg[i,j,k] - bv[1]
+ rvg[2,i,j,k] = vzg[i,j,k] - bv[2]
+ return rvg
diff -r 809bff333e8b5022ed54799209c3e69a668a78f1 -r 28ee8407e8ac563031468f9a1f535e52b3f4eaae yt/utilities/lib/misc_utilities.pyx
--- a/yt/utilities/lib/misc_utilities.pyx
+++ b/yt/utilities/lib/misc_utilities.pyx
@@ -233,49 +233,6 @@
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
-def obtain_rvec(data):
- # This is just to let the pointers exist and whatnot. We can't cdef them
- # inside conditionals.
- cdef np.ndarray[np.float64_t, ndim=1] xf
- cdef np.ndarray[np.float64_t, ndim=1] yf
- cdef np.ndarray[np.float64_t, ndim=1] zf
- cdef np.ndarray[np.float64_t, ndim=2] rf
- cdef np.ndarray[np.float64_t, ndim=3] xg
- cdef np.ndarray[np.float64_t, ndim=3] yg
- cdef np.ndarray[np.float64_t, ndim=3] zg
- cdef np.ndarray[np.float64_t, ndim=4] rg
- cdef np.float64_t c[3]
- cdef int i, j, k
- center = data.get_field_parameter("center")
- c[0] = center[0]; c[1] = center[1]; c[2] = center[2]
- if len(data['x'].shape) == 1:
- # One dimensional data
- xf = data['x']
- yf = data['y']
- zf = data['z']
- rf = np.empty((3, xf.shape[0]), 'float64')
- for i in range(xf.shape[0]):
- rf[0, i] = xf[i] - c[0]
- rf[1, i] = yf[i] - c[1]
- rf[2, i] = zf[i] - c[2]
- return rf
- else:
- # Three dimensional data
- xg = data['x']
- yg = data['y']
- zg = data['z']
- rg = np.empty((3, xg.shape[0], xg.shape[1], xg.shape[2]), 'float64')
- for i in range(xg.shape[0]):
- for j in range(xg.shape[1]):
- for k in range(xg.shape[2]):
- rg[0,i,j,k] = xg[i,j,k] - c[0]
- rg[1,i,j,k] = yg[i,j,k] - c[1]
- rg[2,i,j,k] = zg[i,j,k] - c[2]
- return rg
-
- at cython.boundscheck(False)
- at cython.wraparound(False)
- at cython.cdivision(True)
def kdtree_get_choices(np.ndarray[np.float64_t, ndim=3] data,
np.ndarray[np.float64_t, ndim=1] l_corner,
np.ndarray[np.float64_t, ndim=1] r_corner):
https://bitbucket.org/yt_analysis/yt/changeset/c9edc60694ed/
changeset: c9edc60694ed
branch: yt
user: ngoldbaum
date: 2012-10-06 21:47:32
summary: Forgot a stray tranpose()
affected #: 1 file
diff -r 28ee8407e8ac563031468f9a1f535e52b3f4eaae -r c9edc60694edef326a39c43080831dd472de5583 yt/data_objects/universal_fields.py
--- a/yt/data_objects/universal_fields.py
+++ b/yt/data_objects/universal_fields.py
@@ -191,7 +191,7 @@
def _VelocityMagnitude(field, data):
"""M{|v|}"""
- velocities = obtain_rv_vec(data)
+ velocities = obtain_rv_vec(data).transpose()
return np.sqrt(np.sum(velocities**2,axis=-1))
add_field("VelocityMagnitude", function=_VelocityMagnitude,
take_log=False, units=r"\rm{cm}/\rm{s}")
https://bitbucket.org/yt_analysis/yt/changeset/3b427b92d9da/
changeset: 3b427b92d9da
branch: yt
user: ngoldbaum
date: 2012-10-09 03:53:45
summary: Merging.
affected #: 4 files
diff -r b9838b15dbde418c6f7a8862b1499e67a5d1e6db -r 3b427b92d9dab845eff0d6bd6a3105c90735c014 yt/data_objects/universal_fields.py
--- a/yt/data_objects/universal_fields.py
+++ b/yt/data_objects/universal_fields.py
@@ -32,7 +32,7 @@
from yt.funcs import *
-from yt.utilities.lib import CICDeposit_3, obtain_rvec
+from yt.utilities.lib import CICDeposit_3, obtain_rvec, obtain_rv_vec
from yt.utilities.cosmology import Cosmology
from field_info_container import \
add_field, \
@@ -54,7 +54,19 @@
kboltz, \
G, \
rho_crit_now, \
- speed_of_light_cgs
+ speed_of_light_cgs, \
+ km_per_cm
+
+from yt.utilities.math_utils import \
+ get_sph_r_component, \
+ get_sph_theta_component, \
+ get_sph_phi_component, \
+ get_cyl_r_component, \
+ get_cyl_z_component, \
+ get_cyl_theta_component, \
+ get_cyl_r, get_cyl_theta, \
+ get_cyl_z, get_sph_r, \
+ get_sph_theta, get_sph_phi
# Note that, despite my newfound efforts to comply with PEP-8,
# I violate it here in order to keep the name/func_name relationship
@@ -179,12 +191,8 @@
def _VelocityMagnitude(field, data):
"""M{|v|}"""
- bulk_velocity = data.get_field_parameter("bulk_velocity")
- if bulk_velocity == None:
- bulk_velocity = np.zeros(3)
- return ( (data["x-velocity"]-bulk_velocity[0])**2.0 + \
- (data["y-velocity"]-bulk_velocity[1])**2.0 + \
- (data["z-velocity"]-bulk_velocity[2])**2.0 )**(1.0/2.0)
+ velocities = obtain_rv_vec(data).transpose()
+ return np.sqrt(np.sum(velocities**2,axis=-1))
add_field("VelocityMagnitude", function=_VelocityMagnitude,
take_log=False, units=r"\rm{cm}/\rm{s}")
@@ -194,13 +202,6 @@
function=_TangentialOverVelocityMagnitude,
take_log=False)
-def _TangentialVelocity(field, data):
- return np.sqrt(data["VelocityMagnitude"]**2.0
- - data["RadialVelocity"]**2.0)
-add_field("TangentialVelocity",
- function=_TangentialVelocity,
- take_log=False, units=r"\rm{cm}/\rm{s}")
-
def _Pressure(field, data):
"""M{(Gamma-1.0)*rho*E}"""
return (data.pf["Gamma"] - 1.0) * \
@@ -223,14 +224,9 @@
def _sph_r(field, data):
center = data.get_field_parameter("center")
- coords = np.array([data['x'] - center[0],
- data['y'] - center[1],
- data['z'] - center[2]]).transpose()
+ coords = obtain_rvec(data).transpose()
- ## The spherical coordinates radius is simply the magnitude of the
- ## coords vector.
-
- return np.sqrt(np.sum(coords**2,axis=-1))
+ return get_sph_r(vectors, center)
def _Convert_sph_r_CGS(data):
return data.convert("cm")
@@ -245,20 +241,9 @@
center = data.get_field_parameter("center")
normal = data.get_field_parameter("normal")
- coords = np.array([data['x'] - center[0],
- data['y'] - center[1],
- data['z'] - center[2]]).transpose()
+ coords = obtain_rvec(data).transepose()
- ## The angle (theta) with respect to the normal (J), is the arccos
- ## of the dot product of the normal with the normalized coords
- ## vector.
-
- tile_shape = list(coords.shape)[:-1] + [1]
- J = np.tile(normal,tile_shape)
-
- JdotCoords = np.sum(J*coords,axis=-1)
-
- return np.arccos( JdotCoords / np.sqrt(np.sum(coords**2,axis=-1)) )
+ return get_sph_theta(coords, normal)
add_field("sph_theta", function=_sph_theta,
validators=[ValidateParameter("center"),ValidateParameter("normal")])
@@ -269,54 +254,21 @@
center = data.get_field_parameter("center")
normal = data.get_field_parameter("normal")
- coords = np.array([data['x'] - center[0],
- data['y'] - center[1],
- data['z'] - center[2]]).transpose()
-
- ## We have freedom with respect to what axis (xprime) to define
- ## the disk angle. Here I've chosen to use the axis that is
- ## perpendicular to the normal and the y-axis. When normal ==
- ## y-hat, then set xprime = z-hat. With this definition, when
- ## normal == z-hat (as is typical), then xprime == x-hat.
- ##
- ## The angle is then given by the arctan of the ratio of the
- ## yprime-component and the xprime-component of the coords vector.
+ coords = obtain_rvec(data).transpose()
- xprime = np.cross([0.0,1.0,0.0],normal)
- if np.sum(xprime) == 0: xprime = np.array([0.0, 0.0, 1.0])
- yprime = np.cross(normal,xprime)
-
- tile_shape = list(coords.shape)[:-1] + [1]
- Jx = np.tile(xprime,tile_shape)
- Jy = np.tile(yprime,tile_shape)
-
- Px = np.sum(Jx*coords,axis=-1)
- Py = np.sum(Jy*coords,axis=-1)
-
- return np.arctan2(Py,Px)
+ return get_sph_phi(coords, normal)
add_field("sph_phi", function=_sph_phi,
validators=[ValidateParameter("center"),ValidateParameter("normal")])
-
-
### cylindrical coordinates: R (radius in the cylinder's plane)
def _cyl_R(field, data):
center = data.get_field_parameter("center")
normal = data.get_field_parameter("normal")
- coords = np.array([data['x'] - center[0],
- data['y'] - center[1],
- data['z'] - center[2]]).transpose()
+ coords = obtain_rvec(data).transpose()
- ## The cross product of the normal (J) with the coords vector
- ## gives a vector of magnitude equal to the cylindrical radius.
-
- tile_shape = list(coords.shape)[:-1] + [1]
- J = np.tile(normal,tile_shape)
-
- JcrossCoords = np.cross(J,coords)
- return np.sqrt(np.sum(JcrossCoords**2,axis=-1))
+ return get_cyl_r(coords, normal)
def _Convert_cyl_R_CGS(data):
return data.convert("cm")
@@ -324,6 +276,9 @@
add_field("cyl_R", function=_cyl_R,
validators=[ValidateParameter("center"),ValidateParameter("normal")],
convert_function = _Convert_cyl_R_CGS, units=r"\rm{cm}")
+add_field("cyl_RCode", function=_cyl_R,
+ validators=[ValidateParameter("center"),ValidateParameter("normal")],
+ units=r"Radius (code)")
### cylindrical coordinates: z (height above the cylinder's plane)
@@ -331,17 +286,9 @@
center = data.get_field_parameter("center")
normal = data.get_field_parameter("normal")
- coords = np.array([data['x'] - center[0],
- data['y'] - center[1],
- data['z'] - center[2]]).transpose()
+ coords = obtain_rvec(data).transpose()
- ## The dot product of the normal (J) with the coords vector gives
- ## the cylindrical height.
-
- tile_shape = list(coords.shape)[:-1] + [1]
- J = np.tile(normal,tile_shape)
-
- return np.sum(J*coords,axis=-1)
+ return get_cyl_z(coords, normal)
def _Convert_cyl_z_CGS(data):
return data.convert("cm")
@@ -352,14 +299,17 @@
### cylindrical coordinates: theta (angle in the cylinder's plane)
-### [This is identical to the spherical coordinate's 'phi' angle.]
def _cyl_theta(field, data):
- return data['sph_phi']
+ center = data.get_field_parameter("center")
+ normal = data.get_field_parameter("normal")
+
+ coords = obtain_rvec(data).transpose()
+
+ return get_cyl_theta(coords, normal)
add_field("cyl_theta", function=_cyl_theta,
validators=[ValidateParameter("center"),ValidateParameter("normal")])
-
### The old field DiskAngle is the same as the spherical coordinates'
### 'theta' angle. I'm keeping DiskAngle for backwards compatibility.
def _DiskAngle(field, data):
@@ -392,6 +342,54 @@
ValidateParameter("normal")],
units=r"AU", display_field=False)
+def _cyl_RadialVelocity(field, data):
+ normal = data.get_field_parameter("normal")
+ velocities = obtain_rv_vec(data).transpose()
+
+ theta = np.tile(data['cyl_theta'], (3, 1)).transpose()
+
+ return get_cyl_r_component(velocities, theta, normal)
+
+def _cyl_RadialVelocityABS(field, data):
+ return np.abs(_cyl_RadialVelocity(field, data))
+def _Convert_cyl_RadialVelocityKMS(data):
+ return km_per_cm
+add_field("cyl_RadialVelocity", function=_cyl_RadialVelocity,
+ units=r"\rm{cm}/\rm{s}",
+ validators=[ValidateParameter("normal")])
+add_field("cyl_RadialVelocityABS", function=_cyl_RadialVelocityABS,
+ units=r"\rm{cm}/\rm{s}",
+ validators=[ValidateParameter("normal")])
+add_field("cyl_RadialVelocityKMS", function=_cyl_RadialVelocity,
+ convert_function=_Convert_cyl_RadialVelocityKMS, units=r"\rm{km}/\rm{s}",
+ validators=[ValidateParameter("normal")])
+add_field("cyl_RadialVelocityKMSABS", function=_cyl_RadialVelocityABS,
+ convert_function=_Convert_cyl_RadialVelocityKMS, units=r"\rm{km}/\rm{s}",
+ validators=[ValidateParameter("normal")])
+
+def _cyl_TangentialVelocity(field, data):
+ normal = data.get_field_parameter("normal")
+ velocities = obtain_rv_vec(data).transpose()
+ theta = np.tile(data['cyl_theta'], (3, 1)).transpose()
+
+ return get_cyl_theta_component(velocities, theta, normal)
+
+def _cyl_TangentialVelocityABS(field, data):
+ return np.abs(_cyl_TangentialVelocity(field, data))
+def _Convert_cyl_TangentialVelocityKMS(data):
+ return km_per_cm
+add_field("cyl_TangentialVelocity", function=_cyl_TangentialVelocity,
+ units=r"\rm{cm}/\rm{s}",
+ validators=[ValidateParameter("normal")])
+add_field("cyl_TangentialVelocityABS", function=_cyl_TangentialVelocityABS,
+ units=r"\rm{cm}/\rm{s}",
+ validators=[ValidateParameter("normal")])
+add_field("cyl_TangentialVelocityKMS", function=_cyl_TangentialVelocity,
+ convert_function=_Convert_cyl_TangentialVelocityKMS, units=r"\rm{km}/\rm{s}",
+ validators=[ValidateParameter("normal")])
+add_field("cyl_TangentialVelocityKMSABS", function=_cyl_TangentialVelocityABS,
+ convert_function=_Convert_cyl_TangentialVelocityKMS, units=r"\rm{km}/\rm{s}",
+ validators=[ValidateParameter("normal")])
def _DynamicalTime(field, data):
"""
@@ -640,13 +638,7 @@
take_log=False, display_field=False)
def obtain_velocities(data):
- if data.has_field_parameter("bulk_velocity"):
- bv = data.get_field_parameter("bulk_velocity")
- else: bv = np.zeros(3, dtype='float64')
- xv = data["x-velocity"] - bv[0]
- yv = data["y-velocity"] - bv[1]
- zv = data["z-velocity"] - bv[2]
- return xv, yv, zv
+ return obtain_rv_vec(data)
def _convertSpecificAngularMomentum(data):
return data.convert("cm")
@@ -711,7 +703,7 @@
# convert_function=_convertSpecificAngularMomentum, vector_field=True,
# units=r"\rm{cm}^2/\rm{s}", validators=[ValidateParameter('center')])
def _convertSpecificAngularMomentumKMSMPC(data):
- return data.convert("mpc")/1e5
+ return km_per_cm*data.convert("mpc")
#add_field("ParticleSpecificAngularMomentumKMSMPC",
# function=_ParticleSpecificAngularMomentum, particle_type=True,
# convert_function=_convertSpecificAngularMomentumKMSMPC, vector_field=True,
@@ -883,33 +875,34 @@
display_name = "Radius (code)")
def _RadialVelocity(field, data):
- center = data.get_field_parameter("center")
- bulk_velocity = data.get_field_parameter("bulk_velocity")
- if bulk_velocity == None:
- bulk_velocity = np.zeros(3)
- new_field = ( (data['x']-center[0])*(data["x-velocity"]-bulk_velocity[0])
- + (data['y']-center[1])*(data["y-velocity"]-bulk_velocity[1])
- + (data['z']-center[2])*(data["z-velocity"]-bulk_velocity[2])
- )/data["RadiusCode"]
- if np.any(np.isnan(new_field)): # to fix center = point
- new_field[np.isnan(new_field)] = 0.0
- return new_field
+ normal = data.get_field_parameter("normal")
+ if normal == None:
+ normal = [0,0,1]
+ velocities = obtain_rv_vec(data).transpose()
+ theta = np.tile(data['sph_theta'], (3, 1)).transpose()
+ phi = np.tile(data['sph_phi'], (3, 1)).transpose()
+
+ return get_sph_r_component(velocities, theta, phi, normal)
+
def _RadialVelocityABS(field, data):
return np.abs(_RadialVelocity(field, data))
def _ConvertRadialVelocityKMS(data):
- return 1e-5
+ return km_per_cm
add_field("RadialVelocity", function=_RadialVelocity,
- units=r"\rm{cm}/\rm{s}",
- validators=[ValidateParameter("center")])
+ units=r"\rm{cm}/\rm{s}")
add_field("RadialVelocityABS", function=_RadialVelocityABS,
- units=r"\rm{cm}/\rm{s}",
- validators=[ValidateParameter("center")])
+ units=r"\rm{cm}/\rm{s}")
add_field("RadialVelocityKMS", function=_RadialVelocity,
- convert_function=_ConvertRadialVelocityKMS, units=r"\rm{km}/\rm{s}",
- validators=[ValidateParameter("center")])
+ convert_function=_ConvertRadialVelocityKMS, units=r"\rm{km}/\rm{s}")
add_field("RadialVelocityKMSABS", function=_RadialVelocityABS,
- convert_function=_ConvertRadialVelocityKMS, units=r"\rm{km}/\rm{s}",
- validators=[ValidateParameter("center")])
+ convert_function=_ConvertRadialVelocityKMS, units=r"\rm{km}/\rm{s}")
+
+def _TangentialVelocity(field, data):
+ return np.sqrt(data["VelocityMagnitude"]**2.0
+ - data["RadialVelocity"]**2.0)
+add_field("TangentialVelocity",
+ function=_TangentialVelocity,
+ take_log=False, units=r"\rm{cm}/\rm{s}")
def _CuttingPlaneVelocityX(field, data):
x_vec, y_vec, z_vec = [data.get_field_parameter("cp_%s_vec" % (ax))
@@ -1026,6 +1019,47 @@
display_name=r"\rm{Magnetic}\/\rm{Energy}",
units="\rm{ergs}\/\rm{cm}^{-3}")
+def _BPoloidal(field,data):
+ normal = data.get_field_parameter("normal")
+
+ Bfields = np.array([data['Bx'], data['By'], data['Bz']])
+
+ theta = np.tile(data['sph_theta'], (3, 1)).transpose()
+ phi = np.tile(data['sph_phi'], (3, 1)).transpose()
+
+ return get_sph_theta_component(Bfields, theta, phi, normal)
+
+add_field("BPoloidal", function=_BPoloidal,
+ units=r"\rm{Gauss}",
+ validators=[ValidateParameter("normal")])
+
+def _BToroidal(field,data):
+ normal = data.get_field_parameter("normal")
+
+ Bfields = np.array([data['Bx'], data['By'], data['Bz']])
+
+ phi = np.tile(data['sph_phi'], (3, 1)).transpose()
+
+ return get_sph_phi_component(Bfields, phi, normal)
+
+add_field("BToroidal", function=_BToroidal,
+ units=r"\rm{Gauss}",
+ validators=[ValidateParameter("normal")])
+
+def _BRadial(field,data):
+ normal = data.get_field_parameter("normal")
+
+ Bfields = np.array([data['Bx'], data['By'], data['Bz']])
+
+ theta = np.tile(data['sph_theta'], (3, 1)).transpose()
+ phi = np.tile(data['sph_phi'], (3, 1)).transpose()
+
+ return get_sph_r_component(Bfields, theta, phi, normal)
+
+add_field("BRadial", function=_BPoloidal,
+ units=r"\rm{Gauss}",
+ validators=[ValidateParameter("normal")])
+
def _VorticitySquared(field, data):
mylog.debug("Generating vorticity on %s", data)
# We need to set up stencils
diff -r b9838b15dbde418c6f7a8862b1499e67a5d1e6db -r 3b427b92d9dab845eff0d6bd6a3105c90735c014 yt/utilities/lib/geometry_utils.pyx
--- a/yt/utilities/lib/geometry_utils.pyx
+++ b/yt/utilities/lib/geometry_utils.pyx
@@ -338,3 +338,45 @@
rg[2,i,j,k] = zg[i,j,k] - c[2]
return rg
+ at cython.boundscheck(False)
+ at cython.wraparound(False)
+ at cython.cdivision(True)
+def obtain_rv_vec(data):
+ # This is just to let the pointers exist and whatnot. We can't cdef them
+ # inside conditionals.
+ cdef np.ndarray[np.float64_t, ndim=1] vxf
+ cdef np.ndarray[np.float64_t, ndim=1] vyf
+ cdef np.ndarray[np.float64_t, ndim=1] vzf
+ cdef np.ndarray[np.float64_t, ndim=2] rvf
+ cdef np.ndarray[np.float64_t, ndim=3] vxg
+ cdef np.ndarray[np.float64_t, ndim=3] vyg
+ cdef np.ndarray[np.float64_t, ndim=3] vzg
+ cdef np.ndarray[np.float64_t, ndim=4] rvg
+ cdef np.float64_t bv[3]
+ cdef int i, j, k
+ bulk_velocity = data.get_field_parameter("bulk_velocity")
+ bv[0] = bulk_velocity[0]; bv[1] = bulk_velocity[1]; bv[2] = bulk_velocity[2]
+ if len(data['x-velocity'].shape) == 1:
+ # One dimensional data
+ vxf = data['x-velocity']
+ vyf = data['y-velocity']
+ vzf = data['z-velocity']
+ rvf = np.empty((3, vxf.shape[0]), 'float64')
+ for i in range(vxf.shape[0]):
+ rvf[0, i] = vxf[i] - bv[0]
+ rvf[1, i] = vyf[i] - bv[1]
+ rvf[2, i] = vzf[i] - bv[2]
+ return rvf
+ else:
+ # Three dimensional data
+ vxg = data['x']
+ vyg = data['y']
+ vzg = data['z']
+ rvg = np.empty((3, vxg.shape[0], vxg.shape[1], vxg.shape[2]), 'float64')
+ for i in range(vxg.shape[0]):
+ for j in range(vxg.shape[1]):
+ for k in range(vxg.shape[2]):
+ rvg[0,i,j,k] = vxg[i,j,k] - bv[0]
+ rvg[1,i,j,k] = vyg[i,j,k] - bv[1]
+ rvg[2,i,j,k] = vzg[i,j,k] - bv[2]
+ return rvg
diff -r b9838b15dbde418c6f7a8862b1499e67a5d1e6db -r 3b427b92d9dab845eff0d6bd6a3105c90735c014 yt/utilities/lib/misc_utilities.pyx
--- a/yt/utilities/lib/misc_utilities.pyx
+++ b/yt/utilities/lib/misc_utilities.pyx
@@ -233,49 +233,6 @@
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
-def obtain_rvec(data):
- # This is just to let the pointers exist and whatnot. We can't cdef them
- # inside conditionals.
- cdef np.ndarray[np.float64_t, ndim=1] xf
- cdef np.ndarray[np.float64_t, ndim=1] yf
- cdef np.ndarray[np.float64_t, ndim=1] zf
- cdef np.ndarray[np.float64_t, ndim=2] rf
- cdef np.ndarray[np.float64_t, ndim=3] xg
- cdef np.ndarray[np.float64_t, ndim=3] yg
- cdef np.ndarray[np.float64_t, ndim=3] zg
- cdef np.ndarray[np.float64_t, ndim=4] rg
- cdef np.float64_t c[3]
- cdef int i, j, k
- center = data.get_field_parameter("center")
- c[0] = center[0]; c[1] = center[1]; c[2] = center[2]
- if len(data['x'].shape) == 1:
- # One dimensional data
- xf = data['x']
- yf = data['y']
- zf = data['z']
- rf = np.empty((3, xf.shape[0]), 'float64')
- for i in range(xf.shape[0]):
- rf[0, i] = xf[i] - c[0]
- rf[1, i] = yf[i] - c[1]
- rf[2, i] = zf[i] - c[2]
- return rf
- else:
- # Three dimensional data
- xg = data['x']
- yg = data['y']
- zg = data['z']
- rg = np.empty((3, xg.shape[0], xg.shape[1], xg.shape[2]), 'float64')
- for i in range(xg.shape[0]):
- for j in range(xg.shape[1]):
- for k in range(xg.shape[2]):
- rg[0,i,j,k] = xg[i,j,k] - c[0]
- rg[1,i,j,k] = yg[i,j,k] - c[1]
- rg[2,i,j,k] = zg[i,j,k] - c[2]
- return rg
-
- at cython.boundscheck(False)
- at cython.wraparound(False)
- at cython.cdivision(True)
def kdtree_get_choices(np.ndarray[np.float64_t, ndim=3] data,
np.ndarray[np.float64_t, ndim=1] l_corner,
np.ndarray[np.float64_t, ndim=1] r_corner):
diff -r b9838b15dbde418c6f7a8862b1499e67a5d1e6db -r 3b427b92d9dab845eff0d6bd6a3105c90735c014 yt/utilities/math_utils.py
--- a/yt/utilities/math_utils.py
+++ b/yt/utilities/math_utils.py
@@ -674,3 +674,161 @@
[uz*ux*(1-cost)-uy*sint, uz*uy*(1-cost)+ux*sint, cost+uz**2*(1-cost)]])
return R
+
+def get_ortho_basis(normal):
+ xprime = np.cross([0.0,1.0,0.0],normal)
+ if np.sum(xprime) == 0: xprime = np.array([0.0, 0.0, 1.0])
+ yprime = np.cross(normal,xprime)
+ zprime = normal
+ return (xprime, yprime, zprime)
+
+def get_sph_r(coords):
+ # The spherical coordinates radius is simply the magnitude of the
+ # coordinate vector.
+
+ return np.sqrt(np.sum(coords**2, axis=-1))
+
+
+def get_sph_theta(coords, normal):
+ # The angle (theta) with respect to the normal (J), is the arccos
+ # of the dot product of the normal with the normalized coordinate
+ # vector.
+
+ tile_shape = list(coords.shape)[:-1] + [1]
+ J = np.tile(normal,tile_shape)
+
+ JdotCoords = np.sum(J*coords,axis=-1)
+
+ return np.arccos( JdotCoords / np.sqrt(np.sum(coords**2,axis=-1)) )
+
+def get_sph_phi(coords, normal):
+ # We have freedom with respect to what axis (xprime) to define
+ # the disk angle. Here I've chosen to use the axis that is
+ # perpendicular to the normal and the y-axis. When normal ==
+ # y-hat, then set xprime = z-hat. With this definition, when
+ # normal == z-hat (as is typical), then xprime == x-hat.
+ #
+ # The angle is then given by the arctan of the ratio of the
+ # yprime-component and the xprime-component of the coordinate
+ # vector.
+
+ (xprime, yprime, zprime) = get_ortho_basis(normal)
+
+ tile_shape = list(coords.shape)[:-1] + [1]
+ Jx = np.tile(xprime,tile_shape)
+ Jy = np.tile(yprime,tile_shape)
+
+ Px = np.sum(Jx*coords,axis=-1)
+ Py = np.sum(Jy*coords,axis=-1)
+
+ return np.arctan2(Py,Px)
+
+def get_cyl_r(coords, normal):
+ # The cross product of the normal (J) with a coordinate vector
+ # gives a vector of magnitude equal to the cylindrical radius.
+
+ tile_shape = list(coords.shape)[:-1] + [1]
+ J = np.tile(normal, tile_shape)
+
+ JcrossCoords = np.cross(J, coords)
+ return np.sqrt(np.sum(JcrossCoords**2, axis=-1))
+
+def get_cyl_z(coords, normal):
+ # The dot product of the normal (J) with the coordinate vector
+ # gives the cylindrical height.
+
+ tile_shape = list(coords.shape)[:-1] + [1]
+ J = np.tile(normal, tile_shape)
+
+ return np.sum(J*coords, axis=-1)
+
+def get_cyl_theta(coords, normal):
+ # This is identical to the spherical phi component
+
+ return get_sph_phi(coords, normal)
+
+
+def get_cyl_r_component(vectors, theta, normal):
+ # The r of a vector is the vector dotted with rhat
+
+ (xprime, yprime, zprime) = get_ortho_basis(normal)
+
+ tile_shape = list(vectors.shape)[:-1] + [1]
+ Jx = np.tile(xprime,tile_shape)
+ Jy = np.tile(yprime,tile_shape)
+
+ rhat = Jx*np.cos(theta) + Jy*np.sin(theta)
+
+ return np.sum(vectors*rhat,axis=-1)
+
+def get_cyl_theta_component(vectors, theta, normal):
+ # The theta component of a vector is the vector dotted with thetahat
+
+ (xprime, yprime, zprime) = get_ortho_basis(normal)
+
+ tile_shape = list(vectors.shape)[:-1] + [1]
+ Jx = np.tile(xprime,tile_shape)
+ Jy = np.tile(yprime,tile_shape)
+
+ thetahat = -Jx*np.sin(theta) + Jy*np.cos(theta)
+
+ return np.sum(vectors*thetahat, axis=-1)
+
+def get_cyl_z_component(vectors, normal):
+ # The z component of a vector is the vector dotted with zhat
+ (xprime, yprime, zprime) = get_ortho_basis(normal)
+
+ tile_shape = list(vectors.shape)[:-1] + [1]
+ zhat = np.tile(zprime, tile_shape)
+
+ return np.sum(vectors*zhat, axis=-1)
+
+def get_sph_r_component(vectors, theta, phi, normal):
+ # The r component of a vector is the vector dotted with rhat
+
+ (xprime, yprime, zprime) = get_ortho_basis(normal)
+
+ tile_shape = list(vectors.shape)[:-1] + [1]
+ Jx = np.tile(xprime,tile_shape)
+ Jy = np.tile(yprime,tile_shape)
+ Jz = np.tile(zprime,tile_shape)
+
+ rhat = Jx*np.sin(theta)*np.cos(phi) + \
+ Jy*np.sin(theta)*np.sin(phi) + \
+ Jz*np.cos(theta)
+
+ return np.sum(vectors*rhat, axis=-1)
+
+def get_sph_phi_component(vectors, phi, normal):
+ # The phi component of a vector is the vector dotted with phihat
+
+ phi = get_sph_phi(coords, normal)
+
+ (xprime, yprime, zprime) = get_ortho_basis(normal)
+
+ tile_shape = list(vectors.shape)[:-1] + [1]
+ Jx = np.tile(xprime,tile_shape)
+ Jy = np.tile(yprime,tile_shape)
+
+ phihat = -Jx*np.sin(phi) + Jy*np.cos(phi)
+
+ return np.sum(vectors*phihat, axis=-1)
+
+def get_sph_theta_component(vectors, theta, phi, normal):
+ # The theta component of a vector is the vector dotted with thetahat
+
+ theta = get_sph_theta(coords, normal)
+ phi = get_sph_phi(coords, normal)
+
+ (xprime, yprime, zprime) = get_ortho_basis(normal)
+
+ tile_shape = list(vectors.shape)[:-1] + [1]
+ Jx = np.tile(xprime,tile_shape)
+ Jy = np.tile(yprime,tile_shape)
+ Jz = np.tile(zprime,tile_shape)
+
+ thetahat = Jx*np.cos(theta)*np.cos(phi) + \
+ Jy*np.cos(theta)*np.sin(theta) - \
+ Jz*np.sin(theta)
+
+ return np.sum(vectors*thetahat, axis=-1)
https://bitbucket.org/yt_analysis/yt/changeset/10a6246dd7cd/
changeset: 10a6246dd7cd
branch: yt
user: ngoldbaum
date: 2012-10-09 03:54:00
summary: Merging.
affected #: 1 file
diff -r 3b427b92d9dab845eff0d6bd6a3105c90735c014 -r 10a6246dd7cd27f2fadcff8ec92a3c93a5ffc6d0 yt/visualization/plot_window.py
--- a/yt/visualization/plot_window.py
+++ b/yt/visualization/plot_window.py
@@ -845,17 +845,21 @@
>>> slc.save(mpl_kwargs={'bbox_inches':'tight'})
"""
+ names = []
+ if mpl_kwargs is None: mpl_kwargs = {}
if name == None:
name = str(self.pf)
elif name.endswith('.png'):
- return v.save(name)
- if mpl_kwargs is None: mpl_kwargs = {}
+ for k, v in self.plots.iteritems():
+ names.append(v.save(name,mpl_kwargs))
+ return names
axis = axis_names[self.data_source.axis]
weight = None
type = self._plot_type
if type in ['Projection','OffAxisProjection']:
weight = self.data_source.weight_field
- names = []
+ if 'Cutting' in self.data_source.__class__.__name__:
+ type = 'OffAxisSlice'
for k, v in self.plots.iteritems():
if axis:
n = "%s_%s_%s_%s" % (name, type, axis, k)
https://bitbucket.org/yt_analysis/yt/changeset/d32d3d12f91e/
changeset: d32d3d12f91e
branch: yt
user: ngoldbaum
date: 2012-10-11 21:37:14
summary: Merging.
affected #: 4 files
diff -r 10a6246dd7cd27f2fadcff8ec92a3c93a5ffc6d0 -r d32d3d12f91ef1c0c5450317a269d0480d2ccdeb doc/install_script.sh
--- a/doc/install_script.sh
+++ b/doc/install_script.sh
@@ -220,11 +220,24 @@
echo " * libncurses5-dev"
echo " * zip"
echo " * uuid-dev"
+ echo " * libfreetype6-dev"
+ echo " * tk-dev"
echo
echo "You can accomplish this by executing:"
echo
- echo "$ sudo apt-get install libssl-dev build-essential libncurses5 libncurses5-dev zip uuid-dev"
+ echo "$ sudo apt-get install libssl-dev build-essential libncurses5 libncurses5-dev zip uuid-dev libfreetype6-dev tk-dev"
echo
+ echo
+ echo " Additionally, if you want to put yt's lib dir in your LD_LIBRARY_PATH"
+ echo " so you can use yt without the activate script, you might "
+ echo " want to consider turning off LIBZ and FREETYPE in this"
+ echo " install script by editing this file and setting"
+ echo
+ echo " INST_ZLIB=0"
+ echo " INST_FTYPE=0"
+ echo
+ echo " to avoid conflicts with other command-line programs "
+ echo " (like eog and evince, for example)."
fi
if [ ! -z "${CFLAGS}" ]
then
diff -r 10a6246dd7cd27f2fadcff8ec92a3c93a5ffc6d0 -r d32d3d12f91ef1c0c5450317a269d0480d2ccdeb yt/data_objects/data_containers.py
--- a/yt/data_objects/data_containers.py
+++ b/yt/data_objects/data_containers.py
@@ -855,6 +855,22 @@
for field in temp_data.keys():
self[field] = temp_data[field]
+ def _get_pw(self, fields, center, width, origin, axes_unit, plot_type):
+ axis = self.axis
+ if fields == None:
+ if self.fields == None:
+ raise SyntaxError("The fields keyword argument must be set")
+ else:
+ self.fields = ensure_list(fields)
+ from yt.visualization.plot_window import \
+ GetBoundsAndCenter, PWViewerMPL
+ from yt.visualization.fixed_resolution import FixedResolutionBuffer
+ (bounds, center) = GetBoundsAndCenter(axis, center, width, self.pf)
+ pw = PWViewerMPL(self, bounds, origin=origin, frb_generator=FixedResolutionBuffer,
+ plot_type=plot_type)
+ pw.set_axes_unit(axes_unit)
+ return pw
+
def to_frb(self, width, resolution, center=None, height=None):
r"""This function returns a FixedResolutionBuffer generated from this
object.
@@ -916,26 +932,6 @@
frb = FixedResolutionBuffer(self, bounds, resolution)
return frb
- def to_pw(self):
- r"""Create a :class:`~yt.visualization.plot_window.PlotWindow` from this
- object.
-
- This is a bare-bones mechanism of creating a plot window from this
- object, which can then be moved around, zoomed, and on and on. All
- behavior of the plot window is relegated to that routine.
- """
- axis = self.axis
- center = self.get_field_parameter("center")
- if center is None:
- center = (self.pf.domain_right_edge
- + self.pf.domain_left_edge)/2.0
- width = (1.0, 'unitary')
- from yt.visualization.plot_window import \
- PWViewerMPL, GetBoundsAndCenter
- (bounds, center) = GetBoundsAndCenter(axis, center, width, self.pf)
- pw = PWViewerMPL(self, bounds)
- return pw
-
def interpolate_discretize(self, LE, RE, field, side, log_spacing=True):
"""
This returns a uniform grid of points between *LE* and *RE*,
@@ -1193,6 +1189,18 @@
def hub_upload(self):
self._mrep.upload()
+ def to_pw(self, fields=None, center='c', width=None, axes_unit=None,
+ origin='center-window'):
+ r"""Create a :class:`~yt.visualization.plot_window.PWViewerMPL` from this
+ object.
+
+ This is a bare-bones mechanism of creating a plot window from this
+ object, which can then be moved around, zoomed, and on and on. All
+ behavior of the plot window is relegated to that routine.
+ """
+ pw = self._get_pw(fields, center, width, origin, axes_unit, 'Slice')
+ return pw
+
class AMRCuttingPlaneBase(AMR2DData):
_plane = None
_top_node = "/CuttingPlanes"
@@ -1355,6 +1363,30 @@
return "%s/c%s_L%s" % \
(self._top_node, cen_name, L_name)
+ def to_pw(self, fields=None, center='c', width=None, axes_unit=None):
+ r"""Create a :class:`~yt.visualization.plot_window.PWViewerMPL` from this
+ object.
+
+ This is a bare-bones mechanism of creating a plot window from this
+ object, which can then be moved around, zoomed, and on and on. All
+ behavior of the plot window is relegated to that routine.
+ """
+ normal = self.normal
+ center = self.center
+ if fields == None:
+ if self.fields == None:
+ raise SyntaxError("The fields keyword argument must be set")
+ else:
+ self.fields = ensure_list(fields)
+ from yt.visualization.plot_window import \
+ GetOffAxisBoundsAndCenter, PWViewerMPL
+ from yt.visualization.fixed_resolution import ObliqueFixedResolutionBuffer
+ (bounds, center_rot) = GetOffAxisBoundsAndCenter(normal, center, width, self.pf)
+ pw = PWViewerMPL(self, bounds, origin='center-window', periodic=False, oblique=True,
+ frb_generator=ObliqueFixedResolutionBuffer, plot_type='OffAxisSlice')
+ pw.set_axes_unit(axes_unit)
+ return pw
+
def to_frb(self, width, resolution, height=None):
r"""This function returns an ObliqueFixedResolutionBuffer generated
from this object.
@@ -1762,6 +1794,18 @@
convs[:] = 1.0
return dls, convs
+ def to_pw(self, fields=None, center='c', width=None, axes_unit=None,
+ origin='center-window'):
+ r"""Create a :class:`~yt.visualization.plot_window.PWViewerMPL` from this
+ object.
+
+ This is a bare-bones mechanism of creating a plot window from this
+ object, which can then be moved around, zoomed, and on and on. All
+ behavior of the plot window is relegated to that routine.
+ """
+ pw = self._get_pw(fields, center, width, origin, axes_unit, 'Projection')
+ return pw
+
def get_data(self, fields = None):
if fields is None: fields = ensure_list(self.fields)[:]
else: fields = ensure_list(fields)
@@ -2254,6 +2298,18 @@
def add_fields(self, fields, weight = "CellMassMsun"):
pass
+ def to_pw(self, fields=None, center='c', width=None, axes_unit=None,
+ origin='center-window'):
+ r"""Create a :class:`~yt.visualization.plot_window.PWViewerMPL` from this
+ object.
+
+ This is a bare-bones mechanism of creating a plot window from this
+ object, which can then be moved around, zoomed, and on and on. All
+ behavior of the plot window is relegated to that routine.
+ """
+ pw = self._get_pw(fields, center, width, origin, axes_unit, 'Projection')
+ return pw
+
def _project_grid(self, grid, fields, zero_out):
# We split this next bit into two sections to try to limit the IO load
# on the system. This way, we perserve grid state (@restore_grid_state
diff -r 10a6246dd7cd27f2fadcff8ec92a3c93a5ffc6d0 -r d32d3d12f91ef1c0c5450317a269d0480d2ccdeb yt/visualization/image_writer.py
--- a/yt/visualization/image_writer.py
+++ b/yt/visualization/image_writer.py
@@ -379,7 +379,7 @@
take_log=True)
"""
import matplotlib
- from ._mpl_imports import *
+ from ._mpl_imports import FigureCanvasAgg, FigureCanvasPdf, FigureCanvasPS
# If this is rendered as log, then apply now.
if take_log:
@@ -420,21 +420,22 @@
else:
dpi = None
- if filename[-4:] == '.png':
- suffix = ''
- else:
+ suffix = os.path.splitext(filename)[1]
+
+ if suffix == '':
suffix = '.png'
filename = "%s%s" % (filename, suffix)
- mylog.info("Saving plot %s", fn)
+ mylog.info("Saving plot %s", filename)
if suffix == ".png":
canvas = FigureCanvasAgg(fig)
elif suffix == ".pdf":
canvas = FigureCanvasPdf(fig)
elif suffix in (".eps", ".ps"):
- canvas = FigureCanvasPS
+ canvas = FigureCanvasPS(fig)
else:
mylog.warning("Unknown suffix %s, defaulting to Agg", suffix)
canvas = FigureCanvasAgg(fig)
+
canvas.print_figure(filename)
return filename
diff -r 10a6246dd7cd27f2fadcff8ec92a3c93a5ffc6d0 -r d32d3d12f91ef1c0c5450317a269d0480d2ccdeb yt/visualization/plot_window.py
--- a/yt/visualization/plot_window.py
+++ b/yt/visualization/plot_window.py
@@ -695,6 +695,15 @@
"""
_current_field = None
+ _frb_generator = None
+ _plot_type = None
+
+ def __init__(self, *args, **kwargs):
+ if self._frb_generator == None:
+ self._frb_generator = kwargs.pop("frb_generator")
+ if self._plot_type == None:
+ self._plot_type = kwargs.pop("plot_type")
+ PWViewer.__init__(self, *args, **kwargs)
def _setup_plots(self):
if self._current_field is not None:
@@ -849,7 +858,8 @@
if mpl_kwargs is None: mpl_kwargs = {}
if name == None:
name = str(self.pf)
- elif name.endswith('.png'):
+ suffix = os.path.splitext(name)[1]
+ if suffix != '':
for k, v in self.plots.iteritems():
names.append(v.save(name,mpl_kwargs))
return names
@@ -1224,6 +1234,7 @@
_ext_widget_id = None
_current_field = None
_widget_name = "plot_window"
+ _frb_generator = FixedResolutionBuffer
def _setup_plots(self):
from yt.gui.reason.bottle_mods import PayloadHandler
@@ -1401,24 +1412,25 @@
self.cax = self.figure.add_axes(caxrect)
def save(self, name, mpl_kwargs, canvas = None):
- if name[-4:] == '.png':
- suffix = ''
+ suffix = os.path.splitext(name)[1]
+
+ if suffix == '':
+ suffix = '.png'
+ name = "%s%s" % (name, suffix)
+ mylog.info("Saving plot %s", name)
+ if suffix == ".png":
+ canvas = FigureCanvasAgg(self.figure)
+ elif suffix == ".pdf":
+ canvas = FigureCanvasPdf(self.figure)
+ elif suffix in (".eps", ".ps"):
+ canvas = FigureCanvasPS(self.figure)
else:
- suffix = '.png'
- fn = "%s%s" % (name, suffix)
- mylog.info("Saving plot %s", fn)
- if canvas is None:
- if suffix == ".png":
- canvas = FigureCanvasAgg(self.figure)
- elif suffix == ".pdf":
- canvas = FigureCanvasPdf(self.figure)
- elif suffix in (".eps", ".ps"):
- canvas = FigureCanvasPS
- else:
- mylog.warning("Unknown suffix %s, defaulting to Agg", suffix)
- canvas = FigureCanvasAgg(self.figure)
- canvas.print_figure(fn,**mpl_kwargs)
- return fn
+ mylog.warning("Unknown suffix %s, defaulting to Agg", suffix)
+ canvas = FigureCanvasAgg(self.figure)
+
+
+ canvas.print_figure(name,**mpl_kwargs)
+ return name
def _get_best_layout(self, size):
aspect = 1.0*size[0]/size[1]
https://bitbucket.org/yt_analysis/yt/changeset/6b249ecbf0e6/
changeset: 6b249ecbf0e6
branch: yt
user: ngoldbaum
date: 2012-10-11 21:53:06
summary: Need to check to see if bulk_velocity is defined.
affected #: 1 file
diff -r d32d3d12f91ef1c0c5450317a269d0480d2ccdeb -r 6b249ecbf0e6c643d2054fef8cf5c6657c740b55 yt/utilities/lib/geometry_utils.pyx
--- a/yt/utilities/lib/geometry_utils.pyx
+++ b/yt/utilities/lib/geometry_utils.pyx
@@ -355,6 +355,8 @@
cdef np.float64_t bv[3]
cdef int i, j, k
bulk_velocity = data.get_field_parameter("bulk_velocity")
+ if bulk_velocity == None:
+ bulk_velocity = np.zeros(3)
bv[0] = bulk_velocity[0]; bv[1] = bulk_velocity[1]; bv[2] = bulk_velocity[2]
if len(data['x-velocity'].shape) == 1:
# One dimensional data
https://bitbucket.org/yt_analysis/yt/changeset/ab186e76a23a/
changeset: ab186e76a23a
branch: yt
user: ngoldbaum
date: 2012-10-15 20:55:46
summary: Merging.
affected #: 18 files
diff -r 6b249ecbf0e6c643d2054fef8cf5c6657c740b55 -r ab186e76a23ad019537af4b9b3e05bac62c3284e MANIFEST.in
--- a/MANIFEST.in
+++ b/MANIFEST.in
@@ -1,3 +1,3 @@
include distribute_setup.py
recursive-include yt/gui/reason/html *.html *.png *.ico *.js
-recursive-include yt *.pyx *.pxd *.hh *.h README*
+recursive-include yt *.pyx *.pxd *.hh *.h README* CREDITS FUNDING LICENSE
diff -r 6b249ecbf0e6c643d2054fef8cf5c6657c740b55 -r ab186e76a23ad019537af4b9b3e05bac62c3284e nose.cfg
--- /dev/null
+++ b/nose.cfg
@@ -0,0 +1,4 @@
+[nosetests]
+detailed-errors=1
+where=yt
+exclude=answer_testing
diff -r 6b249ecbf0e6c643d2054fef8cf5c6657c740b55 -r ab186e76a23ad019537af4b9b3e05bac62c3284e setup.cfg
--- a/setup.cfg
+++ b/setup.cfg
@@ -1,3 +1,9 @@
[egg_info]
#tag_build = .dev
#tag_svn_revision = 1
+
+[nosetests]
+detailed-errors=1
+where=yt
+exclude=answer_testing
+with-xunit=1
diff -r 6b249ecbf0e6c643d2054fef8cf5c6657c740b55 -r ab186e76a23ad019537af4b9b3e05bac62c3284e yt/analysis_modules/halo_finding/parallel_hop/parallel_hop_interface.py
--- a/yt/analysis_modules/halo_finding/parallel_hop/parallel_hop_interface.py
+++ b/yt/analysis_modules/halo_finding/parallel_hop/parallel_hop_interface.py
@@ -31,7 +31,7 @@
from yt.funcs import *
from yt.utilities.performance_counters import yt_counters, time_function
try:
- from yt.utilities.kdtree import \
+ from yt.utilities.kdtree.api import \
chainHOP_tags_dens, \
create_tree, fKD, find_nn_nearest_neighbors, \
free_tree, find_chunk_nearest_neighbors
diff -r 6b249ecbf0e6c643d2054fef8cf5c6657c740b55 -r ab186e76a23ad019537af4b9b3e05bac62c3284e yt/analysis_modules/two_point_functions/two_point_functions.py
--- a/yt/analysis_modules/two_point_functions/two_point_functions.py
+++ b/yt/analysis_modules/two_point_functions/two_point_functions.py
@@ -30,7 +30,7 @@
from yt.utilities.parallel_tools.parallel_analysis_interface import ParallelAnalysisInterface, parallel_blocking_call, parallel_root_only
try:
- from yt.utilities.kdtree import *
+ from yt.utilities.kdtree.api import *
except ImportError:
mylog.debug("The Fortran kD-Tree did not import correctly.")
diff -r 6b249ecbf0e6c643d2054fef8cf5c6657c740b55 -r ab186e76a23ad019537af4b9b3e05bac62c3284e yt/config.py
--- a/yt/config.py
+++ b/yt/config.py
@@ -38,6 +38,7 @@
inline = 'False',
numthreads = '-1',
__withinreason = 'False',
+ __withintesting = 'False',
__parallel = 'False',
__global_parallel_rank = '0',
__global_parallel_size = '1',
diff -r 6b249ecbf0e6c643d2054fef8cf5c6657c740b55 -r ab186e76a23ad019537af4b9b3e05bac62c3284e yt/data_objects/data_containers.py
--- a/yt/data_objects/data_containers.py
+++ b/yt/data_objects/data_containers.py
@@ -3658,7 +3658,7 @@
class AMRCoveringGridBase(AMR3DData):
_spatial = True
_type_name = "covering_grid"
- _con_args = ('level', 'left_edge', 'right_edge', 'ActiveDimensions')
+ _con_args = ('level', 'left_edge', 'ActiveDimensions')
def __init__(self, level, left_edge, dims, fields = None,
pf = None, num_ghost_zones = 0, use_pbar = True, **kwargs):
"""A 3D region with all data extracted to a single, specified
@@ -3685,8 +3685,9 @@
fields=fields, pf=pf, **kwargs)
self.left_edge = np.array(left_edge)
self.level = level
- self.dds = self.pf.h.select_grids(self.level)[0].dds.copy()
- self.ActiveDimensions = np.array(dims,dtype='int32')
+ rdx = self.pf.domain_dimensions*self.pf.refine_by**level
+ self.dds = self.pf.domain_width/rdx.astype("float64")
+ self.ActiveDimensions = np.array(dims, dtype='int32')
self.right_edge = self.left_edge + self.ActiveDimensions*self.dds
self._num_ghost_zones = num_ghost_zones
self._use_pbar = use_pbar
diff -r 6b249ecbf0e6c643d2054fef8cf5c6657c740b55 -r ab186e76a23ad019537af4b9b3e05bac62c3284e yt/data_objects/tests/test_covering_grid.py
--- /dev/null
+++ b/yt/data_objects/tests/test_covering_grid.py
@@ -0,0 +1,27 @@
+from yt.testing import *
+from yt.data_objects.profiles import \
+ BinnedProfile1D, BinnedProfile2D, BinnedProfile3D
+
+def setup():
+ from yt.config import ytcfg
+ ytcfg["yt","__withintesting"] = "True"
+
+def test_covering_grid():
+ # We decompose in different ways
+ for level in [0, 1, 2]:
+ for nprocs in [1, 2, 4, 8]:
+ pf = fake_random_pf(16, nprocs = nprocs)
+ dn = pf.refine_by**level
+ cg = pf.h.covering_grid(level, [0.0, 0.0, 0.0],
+ dn * pf.domain_dimensions)
+ yield assert_equal, cg["Ones"].max(), 1.0
+ yield assert_equal, cg["Ones"].min(), 1.0
+ yield assert_equal, cg["CellVolume"].sum(), pf.domain_width.prod()
+ for g in pf.h.grids:
+ di = g.get_global_startindex()
+ dd = g.ActiveDimensions
+ for i in range(dn):
+ f = cg["Density"][dn*di[0]+i:dn*(di[0]+dd[0])+i:dn,
+ dn*di[1]+i:dn*(di[1]+dd[1])+i:dn,
+ dn*di[2]+i:dn*(di[2]+dd[2])+i:dn]
+ yield assert_equal, f, g["Density"]
diff -r 6b249ecbf0e6c643d2054fef8cf5c6657c740b55 -r ab186e76a23ad019537af4b9b3e05bac62c3284e yt/data_objects/tests/test_profiles.py
--- /dev/null
+++ b/yt/data_objects/tests/test_profiles.py
@@ -0,0 +1,74 @@
+from yt.testing import *
+from yt.data_objects.profiles import \
+ BinnedProfile1D, BinnedProfile2D, BinnedProfile3D
+
+_fields = ("Density", "Temperature", "Dinosaurs", "Tribbles")
+
+def test_profiles():
+ pf = fake_random_pf(64, nprocs = 8, fields = _fields)
+ nv = pf.domain_dimensions.prod()
+ dd = pf.h.all_data()
+ (rmi, rma), (tmi, tma), (dmi, dma) = dd.quantities["Extrema"](
+ ["Density", "Temperature", "Dinosaurs"])
+ rt, tt, dt = dd.quantities["TotalQuantity"](
+ ["Density", "Temperature", "Dinosaurs"])
+ # First we look at the
+ for nb in [8, 16, 32, 64]:
+ for lr in [True, False]:
+ # We log all the fields or don't log 'em all. No need to do them
+ # individually.
+ for lf in [True, False]:
+ # We have the min and the max, but to avoid cutting them off
+ # since we aren't doing end-collect, we cut a bit off the edges
+ for ec, e1, e2 in [(False, 0.9, 1.1), (True, 1.0, 1.0)]:
+ p1d = BinnedProfile1D(dd,
+ nb, "Density", rmi*e1, rma*e2, lf,
+ lr, end_collect=ec)
+ p1d.add_fields(["Ones", "Temperature"], weight=None)
+ yield assert_equal, p1d["Ones"].sum(), nv
+ yield assert_rel_equal, tt, p1d["Temperature"].sum(), 7
+
+ p2d = BinnedProfile2D(dd,
+ nb, "Density", rmi*e1, rma*e2, lf,
+ nb, "Temperature", tmi*e1, tma*e2, lf,
+ lr, end_collect=ec)
+ p2d.add_fields(["Ones", "Temperature"], weight=None)
+ yield assert_equal, p2d["Ones"].sum(), nv
+ yield assert_rel_equal, tt, p2d["Temperature"].sum(), 7
+
+ p3d = BinnedProfile3D(dd,
+ nb, "Density", rmi*e1, rma*e2, lf,
+ nb, "Temperature", tmi*e1, tma*e2, lf,
+ nb, "Dinosaurs", dmi*e1, dma*e2, lf,
+ lr, end_collect=ec)
+ p3d.add_fields(["Ones", "Temperature"], weight=None)
+ yield assert_equal, p3d["Ones"].sum(), nv
+ yield assert_rel_equal, tt, p3d["Temperature"].sum(), 7
+
+ p1d = BinnedProfile1D(dd, nb, "x", 0.0, 1.0, log_space=False)
+ p1d.add_fields("Ones", weight=None)
+ av = nv / nb
+ yield assert_equal, p1d["Ones"][:-1], np.ones(nb)*av
+ # We re-bin ones with a weight now
+ p1d.add_fields(["Ones"], weight="Temperature")
+ yield assert_equal, p1d["Ones"][:-1], np.ones(nb)
+
+ p2d = BinnedProfile2D(dd, nb, "x", 0.0, 1.0, False,
+ nb, "y", 0.0, 1.0, False)
+ p2d.add_fields("Ones", weight=None)
+ av = nv / nb**2
+ yield assert_equal, p2d["Ones"][:-1,:-1], np.ones((nb, nb))*av
+ # We re-bin ones with a weight now
+ p2d.add_fields(["Ones"], weight="Temperature")
+ yield assert_equal, p2d["Ones"][:-1,:-1], np.ones((nb, nb))
+
+ p3d = BinnedProfile3D(dd, nb, "x", 0.0, 1.0, False,
+ nb, "y", 0.0, 1.0, False,
+ nb, "z", 0.0, 1.0, False)
+ p3d.add_fields("Ones", weight=None)
+ av = nv / nb**3
+ yield assert_equal, p3d["Ones"][:-1,:-1,:-1], np.ones((nb, nb, nb))*av
+ # We re-bin ones with a weight now
+ p3d.add_fields(["Ones"], weight="Temperature")
+ yield assert_equal, p3d["Ones"][:-1,:-1,:-1], np.ones((nb,nb,nb))
+
diff -r 6b249ecbf0e6c643d2054fef8cf5c6657c740b55 -r ab186e76a23ad019537af4b9b3e05bac62c3284e yt/data_objects/tests/test_projection.py
--- /dev/null
+++ b/yt/data_objects/tests/test_projection.py
@@ -0,0 +1,39 @@
+from yt.testing import *
+from yt.data_objects.profiles import \
+ BinnedProfile1D, BinnedProfile2D, BinnedProfile3D
+
+def setup():
+ from yt.config import ytcfg
+ ytcfg["yt","__withintesting"] = "True"
+
+def test_projection():
+ for nprocs in [8, 1]:
+ # We want to test both 1 proc and 8 procs, to make sure that
+ # parallelism isn't broken
+ pf = fake_random_pf(64, nprocs = nprocs)
+ dims = pf.domain_dimensions
+ xn, yn, zn = pf.domain_dimensions
+ xi, yi, zi = pf.domain_left_edge + 1.0/(pf.domain_dimensions * 2)
+ xf, yf, zf = pf.domain_right_edge - 1.0/(pf.domain_dimensions * 2)
+ dd = pf.h.all_data()
+ rho_tot = dd.quantities["TotalQuantity"]("Density")[0]
+ coords = np.mgrid[xi:xf:xn*1j, yi:yf:yn*1j, zi:zf:zn*1j]
+ uc = [np.unique(c) for c in coords]
+ # Some simple projection tests with single grids
+ for ax, an in enumerate("xyz"):
+ xax = x_dict[ax]
+ yax = y_dict[ax]
+ for wf in ["Density", None]:
+ proj = pf.h.proj(ax, ["Ones", "Density"], weight_field = wf)
+ yield assert_equal, proj["Ones"].sum(), proj["Ones"].size
+ yield assert_equal, proj["Ones"].min(), 1.0
+ yield assert_equal, proj["Ones"].max(), 1.0
+ yield assert_equal, np.unique(proj["px"]), uc[xax]
+ yield assert_equal, np.unique(proj["py"]), uc[yax]
+ yield assert_equal, np.unique(proj["pdx"]), 1.0/(dims[xax]*2.0)
+ yield assert_equal, np.unique(proj["pdy"]), 1.0/(dims[yax]*2.0)
+ # wf == None
+ yield assert_equal, wf, None
+ v1 = proj["Density"].sum()
+ v2 = (dd["Density"] * dd["d%s" % an]).sum()
+ yield assert_rel_equal, v1, v2, 10
diff -r 6b249ecbf0e6c643d2054fef8cf5c6657c740b55 -r ab186e76a23ad019537af4b9b3e05bac62c3284e yt/frontends/stream/data_structures.py
--- a/yt/frontends/stream/data_structures.py
+++ b/yt/frontends/stream/data_structures.py
@@ -353,7 +353,8 @@
psize = get_psize(np.array(data[key].shape), nprocs)
grid_left_edges, grid_right_edges, temp[key] = \
decompose_array(data[key], psize, bbox)
- grid_dimensions = np.array([grid.shape for grid in temp[key]])
+ grid_dimensions = np.array([grid.shape for grid in temp[key]],
+ dtype="int32")
for gid in range(nprocs):
new_data[gid] = {}
for key in temp.keys():
@@ -364,7 +365,7 @@
sfh.update({0:data})
grid_left_edges = domain_left_edge
grid_right_edges = domain_right_edge
- grid_dimensions = domain_dimensions.reshape(nprocs,3)
+ grid_dimensions = domain_dimensions.reshape(nprocs,3).astype("int32")
handler = StreamHandler(
grid_left_edges,
diff -r 6b249ecbf0e6c643d2054fef8cf5c6657c740b55 -r ab186e76a23ad019537af4b9b3e05bac62c3284e yt/funcs.py
--- a/yt/funcs.py
+++ b/yt/funcs.py
@@ -310,7 +310,8 @@
maxval = max(maxval, 1)
from yt.config import ytcfg
if ytcfg.getboolean("yt", "suppressStreamLogging") or \
- ytcfg.getboolean("yt", "ipython_notebook"):
+ ytcfg.getboolean("yt", "ipython_notebook") or \
+ ytcfg.getboolean("yt", "__withintesting"):
return DummyProgressBar()
elif ytcfg.getboolean("yt", "__withinreason"):
from yt.gui.reason.extdirect_repl import ExtProgressBar
diff -r 6b249ecbf0e6c643d2054fef8cf5c6657c740b55 -r ab186e76a23ad019537af4b9b3e05bac62c3284e yt/testing.py
--- a/yt/testing.py
+++ b/yt/testing.py
@@ -24,7 +24,12 @@
import numpy as np
from yt.funcs import *
-from numpy.testing import assert_array_equal
+from numpy.testing import assert_array_equal, assert_almost_equal, \
+ assert_approx_equal, assert_array_almost_equal, assert_equal, \
+ assert_string_equal
+
+def assert_rel_equal(a1, a2, decimels):
+ return assert_almost_equal(a1/a2, 1.0, decimels)
def amrspace(extent, levels=7, cells=8):
"""Creates two numpy arrays representing the left and right bounds of
@@ -127,7 +132,8 @@
return left, right, level
-def fake_random_pf(ndims, peak_value = 1.0, fields = ("Density",), negative = False):
+def fake_random_pf(ndims, peak_value = 1.0, fields = ("Density",),
+ negative = False, nprocs = 1):
from yt.frontends.stream.api import load_uniform_grid
if not iterable(ndims):
ndims = [ndims, ndims, ndims]
@@ -139,5 +145,5 @@
offset = 0.0
data = dict((field, (np.random.random(ndims) - offset) * peak_value)
for field in fields)
- ug = load_uniform_grid(data, ndims, 1.0)
+ ug = load_uniform_grid(data, ndims, 1.0, nprocs = nprocs)
return ug
diff -r 6b249ecbf0e6c643d2054fef8cf5c6657c740b55 -r ab186e76a23ad019537af4b9b3e05bac62c3284e yt/utilities/kdtree/__init__.py
--- a/yt/utilities/kdtree/__init__.py
+++ b/yt/utilities/kdtree/__init__.py
@@ -1,1 +0,0 @@
-from fKDpy import *
diff -r 6b249ecbf0e6c643d2054fef8cf5c6657c740b55 -r ab186e76a23ad019537af4b9b3e05bac62c3284e yt/utilities/kdtree/api.py
--- /dev/null
+++ b/yt/utilities/kdtree/api.py
@@ -0,0 +1,9 @@
+from fKDpy import \
+ chainHOP_tags_dens, \
+ create_tree, \
+ fKD, \
+ find_nn_nearest_neighbors, \
+ free_tree, \
+ find_chunk_nearest_neighbors
+
+
diff -r 6b249ecbf0e6c643d2054fef8cf5c6657c740b55 -r ab186e76a23ad019537af4b9b3e05bac62c3284e yt/utilities/kdtree/test.py
--- a/yt/utilities/kdtree/test.py
+++ /dev/null
@@ -1,58 +0,0 @@
-from Forthon import *
-from fKDpy import *
-import numpy,random
-
-n = 32768
-
-
-fKD.tags = fzeros((64),'i')
-fKD.dist = fzeros((64),'d')
-fKD.pos = fzeros((3,n),'d')
-fKD.nn = 64
-fKD.nparts = n
-fKD.sort = True
-fKD.rearrange = True
-fKD.qv = numpy.array([16./32, 16./32, 16./32])
-
-fp = open('parts.txt','r')
-xpos = []
-ypos = []
-zpos = []
-line = fp.readline()
-while line:
- line = line.split()
- xpos.append(float(line[0]))
- ypos.append(float(line[1]))
- zpos.append(float(line[2]))
- line= fp.readline()
-
-fp.close()
-
-
-for k in range(32):
- for j in range(32):
- for i in range(32):
- fKD.pos[0][i + j*32 + k*1024] = float(i)/32 + 1./64 + 0.0001*random.random()
- fKD.pos[1][i + j*32 + k*1024] = float(j)/32 + 1./64 + 0.0001*random.random()
- fKD.pos[2][i + j*32 + k*1024] = float(k)/32 + 1./64 + 0.0001*random.random()
-
-
-
-#print fKD.pos[0][0],fKD.pos[1][0],fKD.pos[2][0]
-
-create_tree()
-
-
-find_nn_nearest_neighbors()
-
-#print 'next'
-
-#fKD.qv = numpy.array([0., 0., 0.])
-
-#find_nn_nearest_neighbors()
-
-
-#print (fKD.tags - 1)
-#print fKD.dist
-
-free_tree()
diff -r 6b249ecbf0e6c643d2054fef8cf5c6657c740b55 -r ab186e76a23ad019537af4b9b3e05bac62c3284e yt/utilities/tests/test_kdtrees.py
--- /dev/null
+++ b/yt/utilities/tests/test_kdtrees.py
@@ -0,0 +1,93 @@
+"""
+Unit test the kD trees in yt.
+
+Author: Stephen Skory <s at skory.us>
+Affiliation: U of Colorado
+Homepage: http://yt-project.org/
+License:
+ Copyright (C) 2008-2011 Stephen Skory. All Rights Reserved.
+
+ This file is part of yt.
+
+ yt is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 3 of the License, or
+ (at your option) any later version.
+
+ This program is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with this program. If not, see <http://www.gnu.org/licenses/>.
+"""
+
+from yt.testing import *
+
+try:
+ from yt.utilities.kdtree.api import *
+except ImportError:
+ mylog.debug("The Fortran kD-Tree did not import correctly.")
+
+from yt.utilities.spatial import cKDTree
+
+def setup():
+ pass
+
+def test_fortran_tree():
+ r"""This test makes sure that the fortran kdtree is finding the correct
+ nearest neighbors.
+ """
+ # Four points.
+ try:
+ fKD.pos = np.empty((3, 4), dtype='float64', order='F')
+ except NameError:
+ return
+ # Make four points by hand that, in particular, will allow us to test
+ # the periodicity of the kdtree.
+ points = np.array([0.01, 0.5, 0.98, 0.99])
+ fKD.pos[0, :] = points
+ fKD.pos[1, :] = points
+ fKD.pos[2, :] = points
+ fKD.qv = np.empty(3, dtype='float64')
+ fKD.dist = np.empty(4, dtype='float64')
+ fKD.tags = np.empty(4, dtype='int64')
+ fKD.nn = 4
+ fKD.sort = True
+ create_tree(0)
+ # Now we check to make sure that we find the correct nearest neighbors,
+ # which get stored in dist and tags.
+ fKD.qv[:] = 0.999
+ find_nn_nearest_neighbors()
+ # Fix fortran counting.
+ fKD.tags -= 1
+ # Clean up before the tests.
+ free_tree(0)
+ # What the answers should be.
+ dist = np.array([2.43e-04, 3.63e-04, 1.083e-03, 7.47003e-01])
+ tags = np.array([3, 0, 2, 1], dtype='int64')
+ assert_array_almost_equal(fKD.dist, dist)
+ assert_array_equal(fKD.tags, tags)
+
+def test_cython_tree():
+ r"""This test makes sure that the cython kdtree is finding the correct
+ nearest neighbors.
+ """
+ # Four points.
+ pos = np.empty((4, 3), dtype='float64')
+ # Make four points by hand that, in particular, will allow us to test
+ # the periodicity of the kdtree.
+ points = np.array([0.01, 0.5, 0.98, 0.99])
+ pos[:, 0] = points
+ pos[:, 1] = points
+ pos[:, 2] = points
+ kdtree = cKDTree(pos, leafsize = 2)
+ qv = np.array([0.999]*3)
+ res = kdtree.query(qv, 4, period=[1.,1.,1])
+ # What the answers should be.
+ dist = np.array([2.43e-04, 3.63e-04, 1.083e-03, 7.47003e-01])
+ tags = np.array([3, 0, 2, 1], dtype='int64')
+ assert_array_almost_equal(res[0], dist)
+ assert_array_equal(res[1], tags)
+
diff -r 6b249ecbf0e6c643d2054fef8cf5c6657c740b55 -r ab186e76a23ad019537af4b9b3e05bac62c3284e yt/visualization/plot_modifications.py
--- a/yt/visualization/plot_modifications.py
+++ b/yt/visualization/plot_modifications.py
@@ -364,39 +364,30 @@
class StreamlineCallback(PlotCallback):
_type_name = "streamlines"
- def __init__(self, field_x, field_y, factor=6.0, nx=16, ny=16,
- xstart=(0,1), ystart=(0,1), nsample=256,
- start_at_xedge=False, start_at_yedge=False,
- plot_args=None):
+ def __init__(self, field_x, field_y, factor = 16,
+ density = 1, arrowsize = 1, arrowstyle = None,
+ color = None, normalize = False):
"""
- annotate_streamlines(field_x, field_y, factor=6.0, nx=16, ny=16,
- xstart=(0,1), ystart=(0,1), nsample=256,
- start_at_xedge=False, start_at_yedge=False,
- plot_args=None):
+ annotate_streamlines(field_x, field_y, factor = 16, density = 1,
+ arrowsize = 1, arrowstyle = None,
+ color = None, normalize = False):
Add streamlines to any plot, using the *field_x* and *field_y*
- from the associated data, using *nx* and *ny* starting points
- that are bounded by *xstart* and *ystart*. To begin
- streamlines from the left edge of the plot, set
- *start_at_xedge* to True; for the bottom edge, use
- *start_at_yedge*. A line with the qmean vector magnitude will
- cover 1.0/*factor* of the image.
+ from the associated data, skipping every *factor* datapoints like
+ 'quiver'. *density* is the index of the amount of the streamlines.
"""
PlotCallback.__init__(self)
self.field_x = field_x
self.field_y = field_y
- self.xstart = xstart
- self.ystart = ystart
- self.nsample = nsample
+ self.bv_x = self.bv_y = 0
self.factor = factor
- if start_at_xedge:
- self.data_size = (1,ny)
- elif start_at_yedge:
- self.data_size = (nx,1)
- else:
- self.data_size = (nx,ny)
- if plot_args is None: plot_args = {'color':'k', 'linestyle':'-'}
- self.plot_args = plot_args
+ self.dens = density
+ self.arrowsize = arrowsize
+ if arrowstyle is None : arrowstyle='-|>'
+ self.arrowstyle = arrowstyle
+ if color is None : color = "#000000"
+ self.color = color
+ self.normalize = normalize
def __call__(self, plot):
x0, x1 = plot.xlim
@@ -404,43 +395,31 @@
xx0, xx1 = plot._axes.get_xlim()
yy0, yy1 = plot._axes.get_ylim()
plot._axes.hold(True)
- nx = plot.image._A.shape[0]
- ny = plot.image._A.shape[1]
+ nx = plot.image._A.shape[0] / self.factor
+ ny = plot.image._A.shape[1] / self.factor
pixX = _MPL.Pixelize(plot.data['px'],
plot.data['py'],
plot.data['pdx'],
plot.data['pdy'],
- plot.data[self.field_x],
+ plot.data[self.field_x] - self.bv_x,
int(nx), int(ny),
- (x0, x1, y0, y1),)
+ (x0, x1, y0, y1),).transpose()
pixY = _MPL.Pixelize(plot.data['px'],
plot.data['py'],
plot.data['pdx'],
plot.data['pdy'],
- plot.data[self.field_y],
+ plot.data[self.field_y] - self.bv_y,
int(nx), int(ny),
- (x0, x1, y0, y1),)
- r0 = np.mgrid[self.xstart[0]*nx:self.xstart[1]*nx:self.data_size[0]*1j,
- self.ystart[0]*ny:self.ystart[1]*ny:self.data_size[1]*1j]
- lines = np.zeros((self.nsample, 2, self.data_size[0], self.data_size[1]))
- lines[0,:,:,:] = r0
- mag = np.sqrt(pixX**2 + pixY**2)
- scale = np.sqrt(nx*ny) / (self.factor * mag.mean())
- dt = 1.0 / (self.nsample-1)
- for i in range(1,self.nsample):
- xt = lines[i-1,0,:,:]
- yt = lines[i-1,1,:,:]
- ix = np.maximum(np.minimum((xt).astype('int'), nx-1), 0)
- iy = np.maximum(np.minimum((yt).astype('int'), ny-1), 0)
- lines[i,0,:,:] = xt + dt * pixX[ix,iy] * scale
- lines[i,1,:,:] = yt + dt * pixY[ix,iy] * scale
- # scale into data units
- lines[:,0,:,:] = lines[:,0,:,:] * (xx1 - xx0) / nx + xx0
- lines[:,1,:,:] = lines[:,1,:,:] * (yy1 - yy0) / ny + yy0
- for i in range(self.data_size[0]):
- for j in range(self.data_size[1]):
- plot._axes.plot(lines[:,0,i,j], lines[:,1,i,j],
- **self.plot_args)
+ (x0, x1, y0, y1),).transpose()
+ X,Y = (na.linspace(xx0,xx1,nx,endpoint=True),
+ na.linspace(yy0,yy1,ny,endpoint=True))
+ if self.normalize:
+ nn = na.sqrt(pixX**2 + pixY**2)
+ pixX /= nn
+ pixY /= nn
+ plot._axes.streamplot(X,Y, pixX, pixY, density=self.dens,
+ arrowsize=self.arrowsize, arrowstyle=self.arrowstyle,
+ color=self.color, norm=self.normalize)
plot._axes.set_xlim(xx0,xx1)
plot._axes.set_ylim(yy0,yy1)
plot._axes.hold(False)
https://bitbucket.org/yt_analysis/yt/changeset/91ffdf886339/
changeset: 91ffdf886339
branch: yt
user: ngoldbaum
date: 2012-10-16 01:06:32
summary: Adding tests for the coordinate conversion code.
affected #: 2 files
diff -r ab186e76a23ad019537af4b9b3e05bac62c3284e -r 91ffdf8863397a75a1b01134f689ce36e139af03 yt/testing.py
--- a/yt/testing.py
+++ b/yt/testing.py
@@ -26,7 +26,7 @@
from yt.funcs import *
from numpy.testing import assert_array_equal, assert_almost_equal, \
assert_approx_equal, assert_array_almost_equal, assert_equal, \
- assert_string_equal
+ assert_array_less, assert_string_equal
def assert_rel_equal(a1, a2, decimels):
return assert_almost_equal(a1/a2, 1.0, decimels)
@@ -139,11 +139,16 @@
ndims = [ndims, ndims, ndims]
else:
assert(len(ndims) == 3)
- if negative:
- offset = 0.5
- else:
- offset = 0.0
+ if not iterable(negative):
+ negative = [negative for f in fields]
+ assert(len(fields) == len(negative))
+ offsets = []
+ for n in negative:
+ if n:
+ offsets.append(0.5)
+ else:
+ offsets.append(0.0)
data = dict((field, (np.random.random(ndims) - offset) * peak_value)
- for field in fields)
+ for field,offset in zip(fields,offsets))
ug = load_uniform_grid(data, ndims, 1.0, nprocs = nprocs)
return ug
diff -r ab186e76a23ad019537af4b9b3e05bac62c3284e -r 91ffdf8863397a75a1b01134f689ce36e139af03 yt/utilities/math_utils.py
--- a/yt/utilities/math_utils.py
+++ b/yt/utilities/math_utils.py
@@ -717,7 +717,7 @@
tile_shape = list(coords.shape)[:-1] + [1]
Jx = np.tile(xprime,tile_shape)
Jy = np.tile(yprime,tile_shape)
-
+
Px = np.sum(Jx*coords,axis=-1)
Py = np.sum(Jy*coords,axis=-1)
@@ -802,8 +802,6 @@
def get_sph_phi_component(vectors, phi, normal):
# The phi component of a vector is the vector dotted with phihat
- phi = get_sph_phi(coords, normal)
-
(xprime, yprime, zprime) = get_ortho_basis(normal)
tile_shape = list(vectors.shape)[:-1] + [1]
@@ -817,18 +815,15 @@
def get_sph_theta_component(vectors, theta, phi, normal):
# The theta component of a vector is the vector dotted with thetahat
- theta = get_sph_theta(coords, normal)
- phi = get_sph_phi(coords, normal)
-
(xprime, yprime, zprime) = get_ortho_basis(normal)
tile_shape = list(vectors.shape)[:-1] + [1]
Jx = np.tile(xprime,tile_shape)
Jy = np.tile(yprime,tile_shape)
Jz = np.tile(zprime,tile_shape)
-
+
thetahat = Jx*np.cos(theta)*np.cos(phi) + \
- Jy*np.cos(theta)*np.sin(theta) - \
+ Jy*np.cos(theta)*np.sin(phi) - \
Jz*np.sin(theta)
return np.sum(vectors*thetahat, axis=-1)
https://bitbucket.org/yt_analysis/yt/changeset/018cc61ee15a/
changeset: 018cc61ee15a
branch: yt
user: ngoldbaum
date: 2012-10-16 01:07:48
summary: Forgot to hg add.
affected #: 2 files
diff -r 91ffdf8863397a75a1b01134f689ce36e139af03 -r 018cc61ee15ac835460585283838f1a06461ca60 yt/utilities/lib/tests/test_geometry_utils.py
--- /dev/null
+++ b/yt/utilities/lib/tests/test_geometry_utils.py
@@ -0,0 +1,30 @@
+from yt.testing import *
+from yt.utilities.lib import obtain_rvec, obtain_rv_vec
+
+_fields = ("Density", "x-velocity", "y-velocity", "z-velocity")
+
+def test_obtain_rvec():
+ pf = fake_random_pf(64, nprocs=8, fields=_fields,
+ negative = [False, True, True, True])
+
+ dd = pf.h.sphere((0.5,0.5,0.5), 0.2)
+
+ coords = obtain_rvec(dd)
+
+ r = np.sqrt(np.sum(coords*coords,axis=0))
+
+ assert_array_less(r.max(), 0.2)
+
+ assert_array_less(0.0, r.min())
+
+def test_obtain_rv_vec():
+ pf = fake_random_pf(64, nprocs=8, fields=_fields,
+ negative = [False, True, True, True])
+
+ dd = pf.h.all_data()
+
+ vels = obtain_rv_vec(dd)
+
+ assert_array_equal(vels[0,:], dd['x-velocity'])
+ assert_array_equal(vels[1,:], dd['y-velocity'])
+ assert_array_equal(vels[2,:], dd['z-velocity'])
diff -r 91ffdf8863397a75a1b01134f689ce36e139af03 -r 018cc61ee15ac835460585283838f1a06461ca60 yt/utilities/tests/test_coordinate_conversions.py
--- /dev/null
+++ b/yt/utilities/tests/test_coordinate_conversions.py
@@ -0,0 +1,128 @@
+from yt.testing import *
+from yt.utilities.math_utils import \
+ get_sph_r_component, \
+ get_sph_theta_component, \
+ get_sph_phi_component, \
+ get_cyl_r_component, \
+ get_cyl_z_component, \
+ get_cyl_theta_component, \
+ get_cyl_r, get_cyl_theta, \
+ get_cyl_z, get_sph_r, \
+ get_sph_theta, get_sph_phi
+
+# Randomly generated coordinates in the domain [[-1,1],[-1,1],-1,1]]
+coords = np.array([[-0.41503037, -0.22102472, -0.55774212],
+ [ 0.73828247, -0.17913899, 0.64076921],
+ [ 0.08922066, -0.94254844, -0.61774511],
+ [ 0.10173242, -0.95789145, 0.16294352],
+ [ 0.73186508, -0.3109153 , 0.75728738],
+ [ 0.8757989 , -0.41475119, -0.57039201],
+ [ 0.58040762, 0.81969082, 0.46759728],
+ [-0.89983356, -0.9853683 , -0.38355343]])
+
+def test_spherical_coordinate_conversion():
+ normal = [0, 0, 1]
+ real_r = [ 0.72950559, 0.99384957, 1.13047198, 0.97696269,
+ 1.09807968, 1.12445067, 1.10788685, 1.38843954]
+ real_theta = [ 2.44113629, 0.87012028, 2.14891444, 1.4032274 ,
+ 0.80979483, 2.10280198, 1.13507735, 1.85068416]
+ real_phi = [-2.65224483, -0.23804243, -1.47641858, -1.46498842,
+ -0.40172325, -0.4422801 , 0.95466734, -2.31085392]
+
+ calc_r = get_sph_r(coords)
+ calc_theta = get_sph_theta(coords, normal)
+ calc_phi = get_sph_phi(coords, normal)
+
+ assert_array_almost_equal(calc_r, real_r)
+ assert_array_almost_equal(calc_theta, real_theta)
+ assert_array_almost_equal(calc_phi, real_phi)
+
+ normal = [1, 0, 0]
+ real_theta = [ 2.17598842, 0.73347681, 1.49179079, 1.46647589,
+ 0.8412984 , 0.67793705, 1.0193883 , 2.27586987]
+ real_phi = [-0.37729951, -2.86898397, -0.99063518, -1.73928995,
+ -2.75201227,-0.62870527, 2.08920872, -1.19959244]
+
+ calc_theta = get_sph_theta(coords, normal)
+ calc_phi = get_sph_phi(coords, normal)
+
+ assert_array_almost_equal(calc_theta, real_theta)
+ assert_array_almost_equal(calc_phi, real_phi)
+
+def test_cylindrical_coordiante_conversion():
+ normal = [0, 0, 1]
+ real_r = [ 0.47021498, 0.75970506, 0.94676179, 0.96327853,
+ 0.79516968, 0.96904193, 1.00437346, 1.3344104 ]
+ real_theta = [-2.65224483, -0.23804243, -1.47641858, -1.46498842,
+ -0.40172325, -0.4422801 , 0.95466734, -2.31085392]
+ real_z = [-0.55774212, 0.64076921, -0.61774511, 0.16294352,
+ 0.75728738, -0.57039201, 0.46759728, -0.38355343]
+
+ calc_r = get_cyl_r(coords, normal)
+ calc_theta = get_cyl_theta(coords, normal)
+ calc_z = get_cyl_z(coords, normal)
+
+ assert_array_almost_equal(calc_r, real_r)
+ assert_array_almost_equal(calc_theta, real_theta)
+ assert_array_almost_equal(calc_z, real_z)
+
+ normal = [1, 0, 0]
+ real_r = [ 0.59994016, 0.66533898, 1.12694569, 0.97165149,
+ 0.81862843, 0.70524152, 0.94368441, 1.05738542]
+ real_theta = [-0.37729951, -2.86898397, -0.99063518, -1.73928995,
+ -2.75201227, -0.62870527, 2.08920872, -1.19959244]
+ real_z = [-0.41503037, 0.73828247, 0.08922066, 0.10173242,
+ 0.73186508, 0.8757989 , 0.58040762, -0.89983356]
+
+ calc_r = get_cyl_r(coords, normal)
+ calc_theta = get_cyl_theta(coords, normal)
+ calc_z = get_cyl_z(coords, normal)
+
+ assert_array_almost_equal(calc_r, real_r)
+ assert_array_almost_equal(calc_theta, real_theta)
+ assert_array_almost_equal(calc_z, real_z)
+
+def test_spherical_coordinate_projections():
+ normal = [0, 0, 1]
+ theta = get_sph_theta(coords, normal)
+ theta_pass = np.tile(theta, (3, 1)).T
+ phi = get_sph_phi(coords, normal)
+ phi_pass = np.tile(phi, (3, 1)).T
+ zero = np.tile(0,coords.shape[0])
+
+ # Purely radial field
+ vecs = np.array([np.sin(theta)*np.cos(phi), np.sin(theta)*np.sin(phi), np.cos(theta)]).T
+ assert_array_almost_equal(zero, get_sph_theta_component(vecs, theta_pass, phi_pass, normal))
+ assert_array_almost_equal(zero, get_sph_phi_component(vecs, phi_pass, normal))
+
+ # Purely toroidal field
+ vecs = np.array([-np.sin(phi), np.cos(phi), zero]).T
+ assert_array_almost_equal(zero, get_sph_theta_component(vecs, theta_pass, phi_pass, normal))
+ assert_array_almost_equal(zero, get_sph_r_component(vecs, theta_pass, phi_pass, normal))
+
+ # Purely poloidal field
+ vecs = np.array([np.cos(theta)*np.cos(phi), np.cos(theta)*np.sin(phi), -np.sin(theta)]).T
+ assert_array_almost_equal(zero, get_sph_phi_component(vecs, phi_pass, normal))
+ assert_array_almost_equal(zero, get_sph_r_component(vecs, theta_pass, phi_pass, normal))
+
+def test_cylindrical_coordinate_projections():
+ normal = [0, 0, 1]
+ theta = get_cyl_theta(coords, normal)
+ theta_pass = np.tile(theta, (3, 1)).T
+ z = get_cyl_z(coords, normal)
+ zero = np.tile(0, coords.shape[0])
+
+ # Purely radial field
+ vecs = np.array([np.cos(theta), np.sin(theta), zero]).T
+ assert_array_almost_equal(zero, get_cyl_theta_component(vecs, theta_pass, normal))
+ assert_array_almost_equal(zero, get_cyl_z_component(vecs, normal))
+
+ # Purely toroidal field
+ vecs = np.array([-np.sin(theta), np.cos(theta), zero]).T
+ assert_array_almost_equal(zero, get_cyl_z_component(vecs, normal))
+ assert_array_almost_equal(zero, get_cyl_r_component(vecs, theta_pass, normal))
+
+ # Purely z field
+ vecs = np.array([zero, zero, z]).T
+ assert_array_almost_equal(zero, get_cyl_theta_component(vecs, theta_pass, normal))
+ assert_array_almost_equal(zero, get_cyl_r_component(vecs, theta_pass, normal))
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