[yt-svn] commit/yt: ngoldbaum: Merged in MatthewTurk/yt/yt-3.0 (pull request #800)
commits-noreply at bitbucket.org
commits-noreply at bitbucket.org
Wed Apr 9 08:16:38 PDT 2014
1 new commit in yt:
https://bitbucket.org/yt_analysis/yt/commits/b7bc0683be9a/
Changeset: b7bc0683be9a
Branch: yt-3.0
User: ngoldbaum
Date: 2014-04-09 17:16:26
Summary: Merged in MatthewTurk/yt/yt-3.0 (pull request #800)
Spherical slicing and coordinates, making bounds for Cylindrical in Cartesian
Affected #: 8 files
diff -r bf628d9eec62eb47e086acc690d045f89b31ce4e -r b7bc0683be9ad0e96cf6a0fa2b570c92b963d20c yt/data_objects/static_output.py
--- a/yt/data_objects/static_output.py
+++ b/yt/data_objects/static_output.py
@@ -50,6 +50,8 @@
PolarCoordinateHandler
from yt.geometry.cylindrical_coordinates import \
CylindricalCoordinateHandler
+from yt.geometry.spherical_coordinates import \
+ SphericalCoordinateHandler
# We want to support the movie format in the future.
# When such a thing comes to pass, I'll move all the stuff that is contant up
@@ -365,6 +367,8 @@
self.coordinates = CylindricalCoordinateHandler(self)
elif self.geometry == "polar":
self.coordinates = PolarCoordinateHandler(self)
+ elif self.geometry == "spherical":
+ self.coordinates = SphericalCoordinateHandler(self)
else:
raise YTGeometryNotSupported(self.geometry)
diff -r bf628d9eec62eb47e086acc690d045f89b31ce4e -r b7bc0683be9ad0e96cf6a0fa2b570c92b963d20c yt/frontends/gdf/data_structures.py
--- a/yt/frontends/gdf/data_structures.py
+++ b/yt/frontends/gdf/data_structures.py
@@ -163,7 +163,8 @@
_field_info_class = GDFFieldInfo
def __init__(self, filename, dataset_type='grid_data_format',
- storage_filename=None):
+ storage_filename=None, geometry = 'cartesian'):
+ self.geometry = geometry
self.fluid_types += ("gdf",)
Dataset.__init__(self, filename, dataset_type)
self.storage_filename = storage_filename
diff -r bf628d9eec62eb47e086acc690d045f89b31ce4e -r b7bc0683be9ad0e96cf6a0fa2b570c92b963d20c yt/geometry/cartesian_coordinates.py
--- a/yt/geometry/cartesian_coordinates.py
+++ b/yt/geometry/cartesian_coordinates.py
@@ -62,7 +62,8 @@
period = self.period[:2].copy() # dummy here
period[0] = self.period[self.x_axis[dim]]
period[1] = self.period[self.y_axis[dim]]
- period = period.in_units("code_length").d
+ if hasattr(period, 'in_units'):
+ period = period.in_units("code_length").d
buff = _MPL.Pixelize(data_source['px'], data_source['py'],
data_source['pdx'], data_source['pdy'],
data_source[field], size[0], size[1],
diff -r bf628d9eec62eb47e086acc690d045f89b31ce4e -r b7bc0683be9ad0e96cf6a0fa2b570c92b963d20c yt/geometry/cylindrical_coordinates.py
--- a/yt/geometry/cylindrical_coordinates.py
+++ b/yt/geometry/cylindrical_coordinates.py
@@ -21,6 +21,8 @@
_unknown_coord, \
_get_coord_fields
import yt.visualization._MPL as _MPL
+from yt.utilities.lib.misc_utilities import \
+ pixelize_cylinder
#
# Cylindrical fields
#
@@ -71,11 +73,12 @@
units = "code_length**3")
- def pixelize(self, dimension, data_source, field, bounds, size, antialias = True):
+ def pixelize(self, dimension, data_source, field, bounds, size,
+ antialias = True, periodic = True):
ax_name = self.axis_name[dimension]
if ax_name in ('r', 'theta'):
return self._ortho_pixelize(data_source, field, bounds, size,
- antialias)
+ antialias, dimension, periodic)
elif ax_name == "z":
return self._cyl_pixelize(data_source, field, bounds, size,
antialias)
@@ -83,20 +86,26 @@
# Pixelizing along a cylindrical surface is a bit tricky
raise NotImplementedError
- def _ortho_pixelize(self, data_source, field, bounds, size, antialias):
+ def _ortho_pixelize(self, data_source, field, bounds, size, antialias,
+ dim, periodic):
+ period = self.period[:2].copy() # dummy here
+ period[0] = self.period[self.x_axis[dim]]
+ period[1] = self.period[self.y_axis[dim]]
+ if hasattr(period, 'in_units'):
+ period = period.in_units("code_length").d
buff = _MPL.Pixelize(data_source['px'], data_source['py'],
data_source['pdx'], data_source['pdy'],
data_source[field], size[0], size[1],
bounds, int(antialias),
- True, self.period).transpose()
+ period, int(periodic)).transpose()
return buff
def _cyl_pixelize(self, data_source, field, bounds, size, antialias):
buff = pixelize_cylinder(data_source['r'],
data_source['dr'],
data_source['theta'],
- data_source['dtheta'],
- size[0], data_source[field], bounds[0])
+ data_source['dtheta']/2.0, # half-widths
+ size, data_source[field], bounds)
return buff
axis_name = { 0 : 'r', 1 : 'z', 2 : 'theta',
diff -r bf628d9eec62eb47e086acc690d045f89b31ce4e -r b7bc0683be9ad0e96cf6a0fa2b570c92b963d20c yt/geometry/polar_coordinates.py
--- a/yt/geometry/polar_coordinates.py
+++ b/yt/geometry/polar_coordinates.py
@@ -91,8 +91,8 @@
buff = pixelize_cylinder(data_source['r'],
data_source['dr'],
data_source['theta'],
- data_source['dtheta'],
- size[0], data_source[field], bounds[1])
+ data_source['dtheta'] / 2.0, # half-widths
+ size, data_source[field], bounds)
return buff
axis_name = { 0 : 'r', 1 : 'theta', 2 : 'z',
diff -r bf628d9eec62eb47e086acc690d045f89b31ce4e -r b7bc0683be9ad0e96cf6a0fa2b570c92b963d20c yt/geometry/spherical_coordinates.py
--- /dev/null
+++ b/yt/geometry/spherical_coordinates.py
@@ -0,0 +1,162 @@
+"""
+Spherical fields
+
+
+
+
+"""
+
+#-----------------------------------------------------------------------------
+# Copyright (c) 2013, yt Development Team.
+#
+# Distributed under the terms of the Modified BSD License.
+#
+# The full license is in the file COPYING.txt, distributed with this software.
+#-----------------------------------------------------------------------------
+
+import numpy as np
+from .coordinate_handler import \
+ CoordinateHandler, \
+ _unknown_coord, \
+ _get_coord_fields
+import yt.visualization._MPL as _MPL
+from yt.utilities.lib.misc_utilities import \
+ pixelize_cylinder, pixelize_aitoff
+
+class SphericalCoordinateHandler(CoordinateHandler):
+
+ def __init__(self, pf, ordering = 'rtp'):
+ if ordering != 'rtp': raise NotImplementedError
+ super(SphericalCoordinateHandler, self).__init__(pf)
+
+ def setup_fields(self, registry):
+ # return the fields for r, z, theta
+ registry.add_field(("index", "dx"), function=_unknown_coord)
+ registry.add_field(("index", "dy"), function=_unknown_coord)
+ registry.add_field(("index", "dz"), function=_unknown_coord)
+ registry.add_field(("index", "x"), function=_unknown_coord)
+ registry.add_field(("index", "y"), function=_unknown_coord)
+ registry.add_field(("index", "z"), function=_unknown_coord)
+ f1, f2 = _get_coord_fields(0)
+ registry.add_field(("index", "dr"), function = f1,
+ display_field = False,
+ units = "code_length")
+ registry.add_field(("index", "r"), function = f2,
+ display_field = False,
+ units = "code_length")
+
+ f1, f2 = _get_coord_fields(1, "")
+ registry.add_field(("index", "dtheta"), function = f1,
+ display_field = False,
+ units = "")
+ registry.add_field(("index", "theta"), function = f2,
+ display_field = False,
+ units = "")
+
+ f1, f2 = _get_coord_fields(2, "")
+ registry.add_field(("index", "dphi"), function = f1,
+ display_field = False,
+ units = "")
+ registry.add_field(("index", "phi"), function = f2,
+ display_field = False,
+ units = "")
+
+ def _SphericalVolume(field, data):
+ # r**2 sin theta dr dtheta dphi
+ vol = data["index", "r"]**2.0
+ vol *= data["index", "dr"]
+ vol *= np.sin(data["index", "theta"])
+ vol *= data["index", "dtheta"]
+ vol *= data["index", "dphi"]
+ return vol
+ registry.add_field(("index", "cell_volume"),
+ function=_SphericalVolume,
+ units = "code_length**3")
+
+ def pixelize(self, dimension, data_source, field, bounds, size,
+ antialias = True, periodic = True):
+ if dimension == 0:
+ return self._ortho_pixelize(data_source, field, bounds, size,
+ antialias, dimension, periodic)
+ elif dimension in (1, 2):
+ return self._cyl_pixelize(data_source, field, bounds, size,
+ antialias, dimension)
+ else:
+ raise NotImplementedError
+
+ def _ortho_pixelize(self, data_source, field, bounds, size, antialias,
+ dim, periodic):
+ # We should be using fcoords
+ period = self.period[:2].copy() # dummy here
+ period[0] = self.period[self.x_axis[dim]]
+ period[1] = self.period[self.y_axis[dim]]
+ period = period.in_units("code_length").d
+ buff = _MPL.Pixelize(data_source['px'], data_source['py'],
+ data_source['pdx'], data_source['pdy'],
+ data_source[field], size[0], size[1],
+ bounds, int(antialias),
+ period, int(periodic)).transpose()
+ return buff
+
+ def _cyl_pixelize(self, data_source, field, bounds, size, antialias,
+ dimension):
+ if dimension == 1:
+ buff = pixelize_cylinder(data_source['r'],
+ data_source['dr'],
+ data_source['phi'],
+ data_source['dphi'] / 2.0, # half-widths
+ size, data_source[field], bounds)
+ elif dimension == 2:
+ buff = pixelize_cylinder(data_source['r'],
+ data_source['dr'],
+ data_source['theta'],
+ data_source['dtheta'] / 2.0, # half-widths
+ size, data_source[field], bounds)
+ buff = pixelize_cylinder(data_source['r'],
+ data_source['dr'],
+ 2.0*np.pi - data_source['theta'],
+ data_source['dtheta'] / 2.0, # half-widths
+ size, data_source[field], bounds,
+ input_img = buff)
+ else:
+ raise RuntimeError
+ return buff
+
+
+ def convert_from_cartesian(self, coord):
+ raise NotImplementedError
+
+ def convert_to_cartesian(self, coord):
+ raise NotImplementedError
+
+ def convert_to_cylindrical(self, coord):
+ raise NotImplementedError
+
+ def convert_from_cylindrical(self, coord):
+ raise NotImplementedError
+
+ def convert_to_spherical(self, coord):
+ raise NotImplementedError
+
+ def convert_from_spherical(self, coord):
+ raise NotImplementedError
+
+ # Despite being mutables, we uses these here to be clear about how these
+ # are generated and to ensure that they are not re-generated unnecessarily
+ axis_name = { 0 : 'r', 1 : 'theta', 2 : 'phi',
+ 'r' : 'r', 'theta' : 'theta', 'phi' : 'phi',
+ 'R' : 'r', 'Theta' : 'theta', 'Phi' : 'phi'}
+
+ axis_id = { 'r' : 0, 'theta' : 1, 'phi' : 2,
+ 0 : 0, 1 : 1, 2 : 2}
+
+ x_axis = { 'r' : 1, 'theta' : 0, 'phi' : 0,
+ 0 : 1, 1 : 0, 2 : 0}
+
+ y_axis = { 'r' : 2, 'theta' : 2, 'phi' : 1,
+ 0 : 2, 1 : 2, 2 : 1}
+
+ @property
+ def period(self):
+ return self.pf.domain_width
+
diff -r bf628d9eec62eb47e086acc690d045f89b31ce4e -r b7bc0683be9ad0e96cf6a0fa2b570c92b963d20c yt/utilities/lib/misc_utilities.pyx
--- a/yt/utilities/lib/misc_utilities.pyx
+++ b/yt/utilities/lib/misc_utilities.pyx
@@ -477,23 +477,27 @@
np.ndarray[np.float64_t, ndim=1] dradius,
np.ndarray[np.float64_t, ndim=1] theta,
np.ndarray[np.float64_t, ndim=1] dtheta,
- int buff_size,
+ buff_size,
np.ndarray[np.float64_t, ndim=1] field,
- np.float64_t rmax=-1.0) :
+ extents, input_img = None):
cdef np.ndarray[np.float64_t, ndim=2] img
cdef np.float64_t x, y, dx, dy, r0, theta0
+ cdef np.float64_t rmax, x0, y0, x1, y1
cdef np.float64_t r_i, theta_i, dr_i, dtheta_i, dthetamin
cdef int i, pi, pj
- if rmax < 0.0 :
- imax = radius.argmax()
- rmax = radius[imax] + dradius[imax]
+ imax = radius.argmax()
+ rmax = radius[imax] + dradius[imax]
- img = np.zeros((buff_size, buff_size))
- extents = [-rmax, rmax] * 2
- dx = (extents[1] - extents[0]) / img.shape[0]
- dy = (extents[3] - extents[2]) / img.shape[1]
+ if input_img is None:
+ img = np.zeros((buff_size[0], buff_size[1]))
+ img[:] = np.nan
+ else:
+ img = input_img
+ x0, x1, y0, y1 = extents
+ dx = (x1 - x0) / img.shape[0]
+ dy = (y1 - y0) / img.shape[1]
dthetamin = dx / rmax
@@ -513,14 +517,73 @@
continue
x = r_i * math.cos(theta_i)
y = r_i * math.sin(theta_i)
- pi = <int>((x + rmax)/dx)
- pj = <int>((y + rmax)/dy)
- img[pi, pj] = field[i]
+ pi = <int>((x - x0)/dx)
+ pj = <int>((y - y0)/dy)
+ if pi >= 0 and pi < img.shape[0] and \
+ pj >= 0 and pj < img.shape[1]:
+ if img[pi, pj] != img[pi, pj]:
+ img[pi, pj] = 0.0
+ img[pi, pj] = field[i]
r_i += 0.5*dx
theta_i += dthetamin
return img
+ at cython.cdivision(True)
+ at cython.boundscheck(False)
+ at cython.wraparound(False)
+def pixelize_aitoff(np.ndarray[np.float64_t, ndim=1] theta,
+ np.ndarray[np.float64_t, ndim=1] dtheta,
+ np.ndarray[np.float64_t, ndim=1] phi,
+ np.ndarray[np.float64_t, ndim=1] dphi,
+ buff_size,
+ np.ndarray[np.float64_t, ndim=1] field,
+ extents, input_img = None):
+
+ cdef np.ndarray[np.float64_t, ndim=2] img
+ cdef int i, j, nf, fi
+ cdef np.float64_t x, y, z, zb
+ cdef np.float64_t dx, dy, inside
+ cdef np.float64_t theta1, dtheta1, phi1, dphi1
+ cdef np.float64_t theta0, phi0
+ cdef np.float64_t PI = np.pi
+ cdef np.float64_t s2 = math.sqrt(2.0)
+ nf = field.shape[0]
+
+ if input_img is None:
+ img = np.zeros((buff_size[0], buff_size[1]))
+ img[:] = np.nan
+ else:
+ img = input_img
+ dx = 2.0 / (img.shape[0] - 1)
+ dy = 2.0 / (img.shape[1] - 1)
+ for i in range(img.shape[0]):
+ x = (-1.0 + i*dx)*s2*2.0
+ for j in range(img.shape[1]):
+ y = (-1.0 + j * dy)*s2
+ zb = (x*x/8.0 + y*y/2.0 - 1.0)
+ if zb > 0: continue
+ z = (1.0 - (x/4.0)**2.0 - (y/2.0)**2.0)
+ z = z**0.5
+ # Longitude
+ phi0 = (2.0*math.atan(z*x/(2.0 * (2.0*z*z-1.0))) + PI)
+ # Latitude
+ # We shift it into co-latitude
+ theta0 = (math.asin(z*y) + PI/2.0)
+ # Now we just need to figure out which pixel contributes.
+ # We do not have a fast search.
+ for fi in range(nf):
+ theta1 = theta[fi]
+ dtheta1 = dtheta[fi]
+ if not (theta1 - dtheta1 <= theta0 <= theta1 + dtheta1):
+ continue
+ phi1 = phi[fi]
+ dphi1 = dphi[fi]
+ if not (phi1 - dphi1 <= phi0 <= phi1 + dphi1):
+ continue
+ img[i, j] = field[fi]
+ return img
+
#@cython.cdivision(True)
#@cython.boundscheck(False)
#@cython.wraparound(False)
diff -r bf628d9eec62eb47e086acc690d045f89b31ce4e -r b7bc0683be9ad0e96cf6a0fa2b570c92b963d20c yt/visualization/plot_window.py
--- a/yt/visualization/plot_window.py
+++ b/yt/visualization/plot_window.py
@@ -165,8 +165,20 @@
width = get_sanitized_width(axis, width, None, pf)
center = get_sanitized_center(center, pf)
elif pf.geometry in ("polar", "cylindrical"):
+ # Set our default width to be the full domain
width = [pf.domain_right_edge[0]*2.0, pf.domain_right_edge[0]*2.0]
center = pf.arr([0.0, 0.0, 0.0], "code_length")
+ elif pf.geometry == "spherical":
+ if axis == 0:
+ width = pf.domain_width[1], pf.domain_width[2]
+ center = 0.5*(pf.domain_left_edge +
+ pf.domain_right_edge).in_units("code_length")
+ else:
+ # Our default width here is the full domain
+ width = [pf.domain_right_edge[0]*2.0, pf.domain_right_edge[0]*2.0]
+ center = pf.arr([0.0, 0.0, 0.0], "code_length")
+ else:
+ raise NotImplementedError
bounds = (center[x_dict[axis]]-width[0] / 2,
center[x_dict[axis]]+width[0] / 2,
center[y_dict[axis]]-width[1] / 2,
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