[Yt-svn] yt-commit r357 - in trunk/yt: lagos raven
mturk at wrangler.dreamhost.com
mturk at wrangler.dreamhost.com
Fri Jan 11 09:18:42 PST 2008
Author: mturk
Date: Fri Jan 11 09:18:23 2008
New Revision: 357
URL: http://yt.spacepope.org/changeset/357
Log:
Fixes for cutting plane -- previously I had not normalized the plane-coord x
and y vectors, so it was looking all distorted. Also sped things up a bit.
Fixed the problem ( Ticket #59 ) of log scales on phase plots. I thought I was
being a Clever Dan before, but I wasn't, I was just being dumb.
Finally, added some properties to the Enzo3DData types that make them behave
more like hierarchies. This is in preparation for decentralizing hierarchy
support. (First, optimization.)
There are probably some changes scattered throughout related to the upcoming clumpfinder.
Modified:
trunk/yt/lagos/BaseDataTypes.py
trunk/yt/lagos/BaseGridType.py
trunk/yt/lagos/DataReadingFuncs.py
trunk/yt/lagos/DerivedFields.py
trunk/yt/lagos/HierarchyType.py
trunk/yt/lagos/WeaveStrings.py
trunk/yt/raven/PlotCollection.py
trunk/yt/raven/PlotTypes.py
Modified: trunk/yt/lagos/BaseDataTypes.py
==============================================================================
--- trunk/yt/lagos/BaseDataTypes.py (original)
+++ trunk/yt/lagos/BaseDataTypes.py Fri Jan 11 09:18:23 2008
@@ -74,6 +74,9 @@
def set_field_parameter(self, name, val):
self.field_parameters[name] = val
+ def has_field_parameter(self, name):
+ return self.field_parameters.has_key(name)
+
def convert(self, datatype):
return self.hierarchy[datatype]
@@ -372,6 +375,7 @@
def __init__(self, normal, center, fields = None):
Enzo2DData.__init__(self, 4, fields)
self.center = center
+ self.set_field_parameter('center',center)
self._cut_masks = {}
# Let's set up our plane equation
# ax + by + cz + d = 0
@@ -382,9 +386,9 @@
_t = na.cross(self._norm_vec, vecs).sum(axis=1)
ax = nd.maximum_position(_t)
self._x_vec = na.cross(vecs[ax,:], self._norm_vec).ravel()
- self._x_vec /= na.dot(self._x_vec, self._x_vec)
+ self._x_vec /= na.sqrt(na.dot(self._x_vec, self._x_vec))
self._y_vec = na.cross(self._norm_vec, self._x_vec).ravel()
- self._y_vec /= na.dot(self._y_vec, self._y_vec)
+ self._y_vec /= na.sqrt(na.dot(self._y_vec, self._y_vec))
self._rot_mat = na.array([self._x_vec,self._y_vec,self._norm_vec])
self._inv_mat = na.linalg.pinv(self._rot_mat)
self._refresh_data()
@@ -412,12 +416,15 @@
def _get_cut_mask(self, grid):
if self._cut_masks.has_key(grid.id):
return self._cut_masks[grid.id]
- # This is slow.
- D = na.sum(self._norm_vec.reshape((3,1,1,1)) * \
- na.array([grid['x'],grid['y'],grid['z']]), axis=0) + self._d
- diag_dist = na.sqrt(grid['dx'][0]**2.0
- + grid['dy'][0]**2.0
- + grid['dz'][0]**2.0)
+ # This is slow. Suggestions for improvement would be great...
+ ss = grid.ActiveDimensions
+ D = na.ones(ss) * self._d
+ D += (grid['x'][:,0,0] * self._norm_vec[0]).reshape(ss[0],1,1)
+ D += (grid['y'][0,:,0] * self._norm_vec[1]).reshape(1,ss[1],1)
+ D += (grid['z'][0,0,:] * self._norm_vec[2]).reshape(1,1,ss[2])
+ diag_dist = na.sqrt(grid.dx**2.0
+ + grid.dy**2.0
+ + grid.dz**2.0)
cm = na.where(na.abs(D) <= 0.5*diag_dist)
self._cut_masks[grid.id] = cm
return cm
@@ -886,7 +893,9 @@
new_field = na.ones(grid.ActiveDimensions, dtype=dtype) * default_val
np = pointI[0].ravel().size
new_field[pointI] = self[field][i:i+np]
- grid[field] = new_field.copy()
+ if grid.data.has_key(field): del grid.data[field]
+ grid[field] = new_field
+ i += np
def _generate_field(self, field):
if fieldInfo.has_key(field):
@@ -979,6 +988,22 @@
gridRightEdge = property(__get_gridRightEdge, __set_gridRightEdge,
__del_gridRightEdge)
+ def __get_gridLevels(self):
+ if self.__gridLevels == None:
+ self.__gridLevels = na.array([g.Level for g in self._grids])
+ return self.__gridLevels
+
+ def __del_gridLevels(self):
+ del self.__gridLevels
+ self.__gridLevels = None
+
+ def __set_gridLevels(self, val):
+ self.__gridLevels = val
+
+ __gridLevels = None
+ gridLevels = property(__get_gridLevels, __set_gridLevels,
+ __del_gridLevels)
+
class ExtractedRegionBase(Enzo3DData):
def __init__(self, base_region, indices):
cen = base_region.get_field_parameter("center")
@@ -1082,7 +1107,11 @@
self._refresh_data()
def _get_list_of_grids(self, field = None):
- self._grids,ind = self.hierarchy.find_sphere_grids(self.center, self.radius)
+ grids,ind = self.hierarchy.find_sphere_grids(self.center, self.radius)
+ # Now we sort by level
+ grids = grids.tolist()
+ grids.sort(key=lambda x: (x.Level, x.LeftEdge[0], x.LeftEdge[1], x.LeftEdge[2]))
+ self._grids = na.array(grids)
@restore_grid_state
def _get_cut_mask(self, grid, field=None):
@@ -1118,27 +1147,27 @@
def _get_list_of_grids(self):
grids, ind = self.pf.hierarchy.get_box_grids(self.left_edge, self.right_edge)
level_ind = na.where(self.pf.hierarchy.gridLevels.ravel()[ind] <= self.level)
- sort_ind = na.flipud(na.argsort(self.pf.h.gridLevels.ravel()[ind][level_ind]))
+ sort_ind = na.argsort(self.pf.h.gridLevels.ravel()[ind][level_ind])
self._grids = self.pf.hierarchy.grids[ind][level_ind][(sort_ind,)]
def __setup_weave_dict(self, grid):
return {
- 'nxc': int(grid.ActiveDimensions[0]),
- 'nyc': int(grid.ActiveDimensions[1]),
- 'nzc': int(grid.ActiveDimensions[2]),
- 'leftEdgeCoarse': na.array(grid.LeftEdge),
+ 'nx_g': int(grid.ActiveDimensions[0]),
+ 'ny_g': int(grid.ActiveDimensions[1]),
+ 'nz_g': int(grid.ActiveDimensions[2]),
+ 'leftEdgeGrid': na.array(grid.LeftEdge),
'rf': int(grid.dx / self.dx),
- 'dx_c': float(grid.dx),
- 'dy_c': float(grid.dy),
- 'dz_c': float(grid.dz),
- 'dx_f': float(self.dx),
- 'dy_f': float(self.dy),
- 'dz_f': float(self.dz),
- 'cubeLeftEdge': na.array(self.left_edge),
+ 'dx_g': float(grid.dx),
+ 'dy_g': float(grid.dy),
+ 'dz_g': float(grid.dz),
+ 'dx_m': float(self.dx),
+ 'dy_m': float(self.dy),
+ 'dz_m': float(self.dz),
+ 'leftEdgeCube': na.array(self.left_edge),
'cubeRightEdge': na.array(self.right_edge),
- 'nxf': int(self.ActiveDimensions[0]),
- 'nyf': int(self.ActiveDimensions[1]),
- 'nzf': int(self.ActiveDimensions[2]),
+ 'nx_m': int(self.ActiveDimensions[0]),
+ 'ny_m': int(self.ActiveDimensions[1]),
+ 'nz_m': int(self.ActiveDimensions[2]),
'childMask' : grid.child_mask
}
@@ -1164,7 +1193,8 @@
for grid in self._grids:
self._get_data_from_grid(grid, field)
if not na.any(self[field] == -999): break
- if na.any(self[field] == -999) and len(self._grids) > 1:
+ if na.any(self[field] == -999) and self.dx < self.hierarchy.grids[0].dx:
+ print na.where(self[field]==-999)[0].size
raise KeyError
def flush_data(self, field=None):
@@ -1178,7 +1208,6 @@
mylog.info("Flushing field %s to %s possible grids",
field, len(self._grids))
grid_list = self._grids.tolist()
- grid_list.reverse()
for grid in grid_list:
self._flush_data_to_grid(grid, field)
@@ -1200,9 +1229,10 @@
locals_dict = self.__setup_weave_dict(grid)
locals_dict['fieldData'] = grid[field]
locals_dict['cubeData'] = self[field]
- locals_dict['lastLevel'] = int(grid.Level == self.level)
+ locals_dict['lastLevel'] = 1
weave.inline(DataCubeReplaceData,
locals_dict.keys(), local_dict=locals_dict,
compiler='gcc',
type_converters=converters.blitz,
auto_downcast=0, verbose=2)
+ grid[field] = locals_dict['fieldData'][:]
Modified: trunk/yt/lagos/BaseGridType.py
==============================================================================
--- trunk/yt/lagos/BaseGridType.py (original)
+++ trunk/yt/lagos/BaseGridType.py Fri Jan 11 09:18:23 2008
@@ -446,13 +446,16 @@
child_mask = property(fget=_get_child_mask, fdel=_del_child_mask)
child_indices = property(fget=_get_child_indices, fdel = _del_child_indices)
- def retrieve_ghost_zones(self, n_zones, fields):
+ def retrieve_ghost_zones(self, n_zones, fields, all_levels=False):
# We will attempt this by creating a datacube that is exactly bigger
# than the grid by nZones*dx in each direction
new_left_edge = self.LeftEdge - n_zones * self.dx
new_right_edge = self.RightEdge + n_zones * self.dx
# Something different needs to be done for the root grid, though
- cube = self.hierarchy.covering_grid(self.Level,
+ level = self.Level
+ if all_levels:
+ level = self.hierarchy.max_level + 1
+ cube = self.hierarchy.covering_grid(level,
new_left_edge, new_right_edge,
self.ActiveDimensions + 2*n_zones, fields,
num_ghost_zones=n_zones)
Modified: trunk/yt/lagos/DataReadingFuncs.py
==============================================================================
--- trunk/yt/lagos/DataReadingFuncs.py (original)
+++ trunk/yt/lagos/DataReadingFuncs.py Fri Jan 11 09:18:23 2008
@@ -62,7 +62,7 @@
"""
sets = SD.SD(self.filename).datasets()
for set in sets:
- self.readDataFast(set)
+ self[set] = self.readDataFast(set)
def readDataHDF5(self, field):
"""
Modified: trunk/yt/lagos/DerivedFields.py
==============================================================================
--- trunk/yt/lagos/DerivedFields.py (original)
+++ trunk/yt/lagos/DerivedFields.py Fri Jan 11 09:18:23 2008
@@ -142,6 +142,7 @@
doesnt_have = []
for p in self.prop:
if not hasattr(data,p):
+ print type(data), p
doesnt_have.append(p)
if len(doesnt_have) > 0:
raise NeedsProperty(doesnt_have)
@@ -218,13 +219,13 @@
def _GridLevel(field, data):
return na.ones(data["Density"].shape)*(data.Level)
-add_field("GridLevel", validators=[ValidateProperty('Level'),
+add_field("GridLevel", validators=[#ValidateProperty('Level'),
ValidateSpatial(0)])
def _GridIndices(field, data):
return na.ones(data["Density"].shape)*(data.id-1)
-add_field("GridIndices", validators=[ValidateProperty('id'),
- ValidateSpatial(0)])
+add_field("GridIndices", validators=[#ValidateProperty('id'),
+ ValidateSpatial(0)], take_log=False)
def _OnesOverDx(field, data):
return na.ones(data["Density"].shape,
@@ -468,6 +469,39 @@
def _Contours(field, data):
return na.ones(data["Density"].shape)*-1
add_field("Contours", validators=[ValidateSpatial(0)], take_log=False)
+add_field("tempContours", function=_Contours, validators=[ValidateSpatial(0)], take_log=False)
+
+def _SpecificAngularMomentum(field, data):
+ """
+ Calculate the angular velocity. Returns a vector for each cell.
+ """
+ if data.has_field_parameter("bulk_velocity"):
+ bv = data.get_field_parameter("bulk_velocity")
+ xv = data["x-velocity"] - bv[0]
+ yv = data["y-velocity"] - bv[1]
+ zv = data["z-velocity"] - bv[2]
+ else:
+ xv = self["x-velocity"]
+ yv = self["y-velocity"]
+ zv = self["z-velocity"]
+
+ center = data.get_field_parameter('center')
+ coords = na.array([data['x'],data['y'],data['z']])
+ r_vec = coords - na.reshape(center,(3,1))
+ v_vec = na.array([xv,yv,zv])
+ tr = na.zeros(v_vec.shape, 'float64')
+ tr[0,:] = (r_vec[1,:] * v_vec[2,:]) - \
+ (r_vec[2,:] * v_vec[1,:])
+ tr[1,:] = (r_vec[0,:] * v_vec[2,:]) - \
+ (r_vec[2,:] * v_vec[0,:])
+ tr[2,:] = (r_vec[0,:] * v_vec[1,:]) - \
+ (r_vec[1,:] * v_vec[0,:])
+ return tr
+def _convertSpecificAngularMomentum(data):
+ return data.convert("cm")
+add_field("SpecificAngularMomentum",
+ convert_function=_convertSpecificAngularMomentum,
+ units=r"\rm{cm}^2/\rm{s}", validators=[ValidateParameter('center')])
def _Radius(field, data):
center = data.get_field_parameter("center")
@@ -523,16 +557,22 @@
+ (data['y']-center[1])*data["y-velocity"]
+ (data['z']-center[2])*data["z-velocity"])/data["RadiusCode"]
return new_field
+def _RadialVelocityABS(field, data):
+ return na.abs(_RadialVelocity(field, data))
def _ConvertRadialVelocity(data):
- return (data.convert("x-velocity") / 1e5)
-def _ConvertRadialVelocityCGS(data):
return (data.convert("x-velocity"))
+def _ConvertRadialVelocityKMS(data):
+ return (data.convert("x-velocity") / 1e5)
add_field("RadialVelocity", function=_RadialVelocity,
- convert_function=_ConvertRadialVelocity, units=r"\rm{km}/\rm{s}",
+ convert_function=_ConvertRadialVelocity, units=r"\rm{cm}/\rm{s}",
+ validators=[ValidateParameter("center"),
+ ValidateParameter("bulk_velocity")])
+add_field("RadialVelocityABS", function=_RadialVelocityABS,
+ convert_function=_ConvertRadialVelocity, units=r"\rm{cm}/\rm{s}",
validators=[ValidateParameter("center"),
ValidateParameter("bulk_velocity")])
-add_field("RadialVelocityCGS", function=_RadialVelocity,
- convert_function=_ConvertRadialVelocityCGS, units=r"\rm{cm}/\rm{s}",
+add_field("RadialVelocityKMS", function=_RadialVelocity,
+ convert_function=_ConvertRadialVelocityKMS, units=r"\rm{km}/\rm{s}",
validators=[ValidateParameter("center"),
ValidateParameter("bulk_velocity")])
@@ -563,11 +603,11 @@
def _convertEnergy(data):
return data.convert("x-velocity")**2.0
-fieldInfo["Gas_Energy"].units = r"\rm{ergs}/\rm{g}"
+fieldInfo["Gas_Energy"]._units = r"\rm{ergs}/\rm{g}"
fieldInfo["Gas_Energy"]._convert_function = _convertEnergy
-fieldInfo["Total_Energy"].units = r"\rm{ergs}/\rm{g}"
+fieldInfo["Total_Energy"]._units = r"\rm{ergs}/\rm{g}"
fieldInfo["Total_Energy"]._convert_function = _convertEnergy
-fieldInfo["Temperature"].units = r"\rm{K}"
+fieldInfo["Temperature"]._units = r"\rm{K}"
def _convertVelocity(data):
return data.convert("x-velocity")
Modified: trunk/yt/lagos/HierarchyType.py
==============================================================================
--- trunk/yt/lagos/HierarchyType.py (original)
+++ trunk/yt/lagos/HierarchyType.py Fri Jan 11 09:18:23 2008
@@ -377,6 +377,7 @@
# For multiple grids on the root level
if v == -1: self.gridReverseTree[i] = None
self.maxLevel = self.gridLevels.max()
+ self.max_level = self.maxLevel
# Now we do things that we need all the grids to do
#self.fieldList = self.grids[0].getFields()
# The rest of this can probably be done with list comprehensions, but
Modified: trunk/yt/lagos/WeaveStrings.py
==============================================================================
--- trunk/yt/lagos/WeaveStrings.py (original)
+++ trunk/yt/lagos/WeaveStrings.py Fri Jan 11 09:18:23 2008
@@ -123,48 +123,46 @@
_baseDataCubeRefineCoarseData = \
r"""
-// For every cell in fieldData, there are (dx_c / dx_f)**3.0 cells that it maps
-// to in our cube, where dx_f is the fine resolution and dx_c is the coarse
-// resolution.
#define MIN(A,B) ((A) < (B) ? (A) : (B))
#define MAX(A,B) ((A) < (B) ? (B) : (A))
-int bx_f, by_f, bz_f, xc, yc, zc, xf, yf, zf;
-int ixf, iyf, izf;
-int max_bx_f, max_by_f, max_bz_f;
+int bx_m, by_m, bz_m, xg, yg, zg, xm, ym, zm;
+int ixm, iym, izm;
+int max_bx_m, max_by_m, max_bz_m;
-for (xc = 0; xc < nxc; xc++)
- for (yc = 0; yc < nyc; yc++)
- for (zc = 0; zc < nzc; zc++) {
- if (!lastLevel)
- if (childMask(xc, yc, zc) == 0) continue;
- // Now we have a cell. What are the left edges?
- if ((leftEdgeCoarse(0)+dx_c*xc > cubeRightEdge(0))
- || (leftEdgeCoarse(1)+dy_c*yc > cubeRightEdge(1))
- || (leftEdgeCoarse(2)+dz_c*zc > cubeRightEdge(2))) continue;
- if ((leftEdgeCoarse(0)+dx_c*(xc+1) < cubeLeftEdge(0))
- || (leftEdgeCoarse(1)+dy_c*(yc+1) < cubeLeftEdge(1))
- || (leftEdgeCoarse(2)+dz_c*(zc+1) < cubeLeftEdge(2))) continue;
- bx_f = MAX(floorl((leftEdgeCoarse(0)+dx_c*xc - cubeLeftEdge(0))/dx_f),0);
- by_f = MAX(floorl((leftEdgeCoarse(1)+dy_c*yc - cubeLeftEdge(1))/dy_f),0);
- bz_f = MAX(floorl((leftEdgeCoarse(2)+dz_c*zc - cubeLeftEdge(2))/dz_f),0);
- max_bx_f = MIN(floorl((leftEdgeCoarse(0)+dx_c*(xc+1) - cubeLeftEdge(0))/dx_f),nxf);
- max_by_f = MIN(floorl((leftEdgeCoarse(1)+dy_c*(yc+1) - cubeLeftEdge(1))/dy_f),nyf);
- max_bz_f = MIN(floorl((leftEdgeCoarse(2)+dz_c*(zc+1) - cubeLeftEdge(2))/dz_f),nzf);
- for (xf = bx_f; xf < max_bx_f ; xf++) {
- if (xf < 0 || xf > nxf) continue;
- for (yf = by_f; yf < max_by_f ; yf++) {
- if (yf < 0 || yf > nyf) continue;
- for (zf = bz_f; zf < max_bz_f ; zf++) {
- if (zf < 0 || zf > nzf) continue;
- %s
- }
+for (xg = 0; xg < nx_g; xg++) {
+ if (leftEdgeGrid(0)+dx_g*xg > cubeRightEdge(0)) continue;
+ if (leftEdgeGrid(0)+dx_g*(xg+1) < leftEdgeCube(0)) continue;
+ bx_m = MAX(ceill((leftEdgeGrid(0)+dx_g*xg - leftEdgeCube(0))/dx_m),0);
+ max_bx_m = MIN(floorl((leftEdgeGrid(0)+dx_g*(xg+1) - leftEdgeCube(0))/dx_m),nx_m);
+ for (yg = 0; yg < ny_g; yg++) {
+ if (leftEdgeGrid(1)+dy_g*yg > cubeRightEdge(1)) continue;
+ if (leftEdgeGrid(1)+dy_g*(yg+1) < leftEdgeCube(1)) continue;
+ by_m = MAX(ceill((leftEdgeGrid(1)+dy_g*yg - leftEdgeCube(1))/dy_m),0);
+ max_by_m = MIN(floorl((leftEdgeGrid(1)+dy_g*(yg+1) - leftEdgeCube(1))/dy_m),ny_m);
+ for (zg = 0; zg < nz_g; zg++) {
+ if (!lastLevel)
+ if (childMask(xg, yg, zg) == 0) continue;
+ if (leftEdgeGrid(2)+dz_g*zg > cubeRightEdge(2)) continue;
+ if (leftEdgeGrid(2)+dz_g*(zg+1) < leftEdgeCube(2)) continue;
+ bz_m = MAX(ceill((leftEdgeGrid(2)+dz_g*zg - leftEdgeCube(2))/dz_m),0);
+ max_bz_m = MIN(floorl((leftEdgeGrid(2)+dz_g*(zg+1) - leftEdgeCube(2))/dz_m),nz_m);
+ for (xm = bx_m; xm < max_bx_m ; xm++) {
+ //if (xm < 0 || xm > nx_m) continue;
+ for (ym = by_m; ym < max_by_m ; ym++) {
+ //if (ym < 0 || ym > ny_m) continue;
+ for (zm = bz_m; zm < max_bz_m ; zm++) {
+ //if (zm < 0 || zm > nz_m) continue;
+ %s
}
}
+ }
}
+ }
+}
"""
DataCubeRefineCoarseData = _baseDataCubeRefineCoarseData % \
- "cubeData(xf,yf,zf) = fieldData(xc,yc,zc);"
+ ("cubeData(xm,ym,zm) = fieldData(xg,yg,zg);")
DataCubeReplaceData = _baseDataCubeRefineCoarseData % \
- "fieldData(xc,yc,zc) = cubeData(xf,yf,zf);"
+ ("fieldData(xg,yg,zg) = cubeData(xm,ym,zm);")
Modified: trunk/yt/raven/PlotCollection.py
==============================================================================
--- trunk/yt/raven/PlotCollection.py (original)
+++ trunk/yt/raven/PlotCollection.py Fri Jan 11 09:18:23 2008
@@ -36,7 +36,7 @@
self.plots = []
self._run_id = deliverator_id
self.pf = pf
- if not center:
+ if center == None:
v,self.c = pf.h.findMax("Density") # @todo: ensure no caching
else:
self.c = center
@@ -117,10 +117,13 @@
def add_cutting_plane(self, field, normal,
center=None, use_colorbar=True,
- figure = None, axes = None, fig_size=None):
+ figure = None, axes = None, fig_size=None, obj=None):
if center == None:
center = self.c
- cp = self.pf.hierarchy.cutting(normal, center, field)
+ if not obj:
+ cp = self.pf.hierarchy.cutting(normal, center, field)
+ else:
+ cp = obj
p = self._add_plot(PlotTypes.CuttingPlanePlot(cp, field,
use_colorbar=use_colorbar, axes=axes, figure=figure,
size=fig_size))
Modified: trunk/yt/raven/PlotTypes.py
==============================================================================
--- trunk/yt/raven/PlotTypes.py (original)
+++ trunk/yt/raven/PlotTypes.py Fri Jan 11 09:18:23 2008
@@ -412,17 +412,20 @@
pxs, pys, pzs = self.data['px'], self.data['py'], self.data['pz']
xs, ys, zs = self.data['x'], self.data['y'], self.data['z']
dxs, dys, dzs = self.data['pdx'], self.data['pdy'], self.data['pdz']
- ds = self.data[self.axis_names['Z']]
- nx = ds.size
+ field = self.axis_names['Z']
+ ds = self.data[field]
+ indices = na.argsort(dxs)[::-1]
+ nx = indices.size
inv_mat = self.data._inv_mat
center = na.array(self.data.center)
buff = na.zeros((width,height), dtype='float64')
+ count = na.zeros((width,height), dtype='float64')
weave.inline(_pixelize_cp,
['pxs','pys','pzs','xs','ys','zs','dxs','dys','dzs',
- 'buff','ds','nx','inv_mat','width','height',
- 'px_min','px_max','py_min','py_max', 'center'],
+ 'buff','ds','nx','inv_mat','width','height','count',
+ 'px_min','px_max','py_min','py_max', 'center', 'indices'],
compiler='gcc', type_converters=converters.blitz,
- auto_downcast = 0)
+ auto_downcast = 0, verbose=2)
return buff
def _refresh_display_width(self, width=None):
@@ -457,11 +460,9 @@
self.axis_names["Y"] = fields[1]
self.axis_names["Z"] = fields[2]
- log_field, self.x_v, self.x_bins = self.setup_bins(self.fields[0],
- self._axes.set_xscale)
- log_field, self.y_v, self.y_bins = self.setup_bins(self.fields[1],
- self._axes.set_yscale)
- self.log_z, self.z_v, self.z_bins = self.setup_bins(self.fields[2])
+ self._log_x, self.x_v, self.x_bins = self.setup_bins(self.fields[0])
+ self._log_y, self.y_v, self.y_bins = self.setup_bins(self.fields[1])
+ self._log_z, self.z_v, self.z_bins = self.setup_bins(self.fields[2])
self.colorbar = None
@@ -494,19 +495,17 @@
def switch_x(self, field):
self.fields[0] = field
self.axis_names["X"] = field
- log_field, self.x_v, self.x_bins = self.setup_bins(self.fields[0],
- self._axes.set_xscale)
+ self._log_x, self.x_v, self.x_bins = self.setup_bins(self.fields[0])
def switch_y(self, field):
self.fields[1] = field
self.axis_names["Y"] = field
- log_field, self.y_v, self.y_bins = self.setup_bins(self.fields[1],
- self._axes.set_yscale)
+ self._log_y, self.y_v, self.y_bins = self.setup_bins(self.fields[1])
def switch_z(self, field):
self.fields[2] = field
self.axis_names["Z"] = field
- self.log_z, self.z_v, self.z_bins = self.setup_bins(self.fields[2])
+ self._log_z, self.z_v, self.z_bins = self.setup_bins(self.fields[2])
def switch_weight(self, weight):
if weight == "": weight=None
@@ -562,7 +561,7 @@
vmin = na.nanmin(vals[vit])
vmax = na.nanmax(vals[vit])
vals[vi] = 0.0
- if self.log_z:
+ if self._log_z:
# We want smallest non-zero vmin
vmin=vals[vals>0.0].min()
self.norm=matplotlib.colors.LogNorm(vmin=vmin, vmax=vmax,
@@ -579,16 +578,12 @@
self.cmap.set_under("w")
self.cmap.set_over("w")
self._axes.clear()
- xs = self._axes.get_xscale()
- ys = self._axes.get_yscale()
- self._axes.set_xscale("linear")
- self._axes.set_yscale("linear")
self.image = self._axes.pcolormesh(self.x_bins, self.y_bins, \
vals,shading='flat', \
norm=self.norm, cmap=self.cmap)
+ self._axes.set_xscale("log" if self._log_x else "linear")
+ self._axes.set_yscale("log" if self._log_y else "linear")
self.vals = vals
- self._axes.set_xscale(xs)
- self._axes.set_yscale(ys)
#self.ticker = matplotlib.ticker.LogLocator(subs=[0.25, 0.5, 0.75, 1])
if self.colorbar == None:
@@ -707,9 +702,11 @@
_pixelize_cp = r"""
-long double md, lxpx, rxpx, lypx, rypx;
-long double lrx, lry, lrz, rrx, rry, rrz;
-int lc, lr, rc, rr;
+long double md, cxpx, cypx;
+long double cx, cy, cz;
+long double lrx, lry, lrz;
+long double rrx, rry, rrz;
+int lc, lr, rc, rr, p;
long double px_dx, px_dy, px_dz, overlap1, overlap2, overlap3;
px_dx = (px_max-px_min)/height;
@@ -720,35 +717,31 @@
#define max(X,Y) ((X) > (Y) ? (X) : (Y))
using namespace std;
-for(int i=0; i<width; i++) for(int j=0; j<height; j++) buff(i,j)=0.0;
+//for(int i=0; i<width; i++) for(int j=0; j<height; j++) count(i,j)=buff(i,j)=0.0;
-for(int p=0; p<nx; p++)
+for(int pp=0; pp<nx; pp++)
{
+ p = indices(pp);
// Any point we want to plot is at most this far from the center
- md = sqrt(dxs(p)*dxs(p) + dys(p)*dys(p) + dzs(p)*dzs(p));
+ md = 2.0*sqrtl(dxs(p)*dxs(p) + dys(p)*dys(p) + dzs(p)*dzs(p));
if(((pxs(p)+md<px_min) ||
(pxs(p)-md>px_max)) ||
((pys(p)+md<py_min) ||
(pys(p)-md>py_max))) continue;
- lc = max(((pxs(p)-md-px_min)/px_dx),0);
- lr = max(((pys(p)-md-py_min)/px_dy),0);
- rc = min(((pxs(p)+md-px_min)/px_dx),height);
- rr = min(((pys(p)+md-py_min)/px_dy),width);
+ lc = max(floorl((pxs(p)-md-px_min)/px_dx),0);
+ lr = max(floorl((pys(p)-md-py_min)/px_dy),0);
+ rc = min(ceill((pxs(p)+md-px_min)/px_dx),height);
+ rr = min(ceill((pys(p)+md-py_min)/px_dy),width);
for (int i=lr;i<rr;i++) {
- lypx = px_dy * i + py_min;
- rypx = px_dy * (i+1) + py_min;
+ cypx = px_dy * (i+0.5) + py_min;
for (int j=lc;j<rc;j++) {
- lxpx = px_dx * j + px_min;
- rxpx = px_dx * (j+1) + px_min;
- lrx = inv_mat(0,0)*lxpx + inv_mat(0,1)*lypx + center(0);
- lry = inv_mat(1,0)*lxpx + inv_mat(1,1)*lypx + center(1);
- lrz = inv_mat(2,0)*lxpx + inv_mat(2,1)*lypx + center(2);
- rrx = inv_mat(0,0)*rxpx + inv_mat(0,1)*rypx + center(0);
- rry = inv_mat(1,0)*rxpx + inv_mat(1,1)*rypx + center(1);
- rrz = inv_mat(2,0)*rxpx + inv_mat(2,1)*rypx + center(2);
- if(((abs(xs(p)-lrx)>1.01*dxs(p)) && (abs(xs(p)-rrx)>1.01*dxs(p))) ||
- ((abs(ys(p)-lry)>1.01*dys(p)) && (abs(ys(p)-rry)>1.01*dys(p))) ||
- ((abs(zs(p)-lrz)>1.01*dzs(p)) && (abs(zs(p)-rrz)>1.01*dzs(p)))) continue;
+ cxpx = px_dx * (j+0.5) + px_min;
+ cx = inv_mat(0,0)*cxpx + inv_mat(0,1)*cypx + center(0);
+ cy = inv_mat(1,0)*cxpx + inv_mat(1,1)*cypx + center(1);
+ cz = inv_mat(2,0)*cxpx + inv_mat(2,1)*cypx + center(2);
+ if( ((xs(p)-cx)>1.01*dxs(p)) || ((xs(p)-cx)<(-1.01*dxs(p)))
+ || ((ys(p)-cy)>1.01*dys(p)) || ((ys(p)-cy)<(-1.01*dys(p)))
+ || ((zs(p)-cz)>1.01*dzs(p)) || ((zs(p)-cz)<(-1.01*dzs(p))) ) continue;
buff(i,j) = ds(p);
}
}
More information about the yt-svn
mailing list