Hi Stephen, Sure, here it is. (I should give credit to Matt for a large chunk of the script.) Thanks! -Ryan ------------------------------------------------------------------------------------------------ from yt.mods import * from yt.utilities.math_utils import ortho_find import yt.visualization.volume_rendering.camera as camera from yt.utilities.lib import healpix_aitoff_proj, pixelize_healpix from itertools import izip, chain, repeat import h5py pf = load("~/enzo/data/cosmo/r0058/redshift0058") com = na.array([ 0.40329299, 0.47176085, 0.46132502]) L = na.array([ 0.25527702, -0.79298414, -0.55318153]) # We now have our angular momentum vector. v1, v2, v3 = ortho_find(L) # v1 = L, v2 & v3 = are two vectors orthogonal to L c = v2 * 8.0/pf['kpc'] + com print "L =", L print "v1, v2, v3 =", v1, v2, v3 sp = pf.h.sphere(c, (1.0, 'kpc')) stars = sp["particle_type"] == 2 xv, yv, zv = [sp["particle_velocity_%s" % ax] for ax in 'xyz'] bv = na.array([na.average(xv), na.average(yv), na.average(zv)]) rotation = na.array([v2,v3,v1]) sphvir = pf.h.sphere(c, (100.0, 'kpc')) rdisk = 18./pf['kpc'] hdisk = 4./pf['kpc'] disk = pf.h.disk(c, v1, rdisk, hdisk) halo = pf.h.boolean([sphvir, "NOT", disk]) # exclude the disk from the sphere for g in pf.h.grids: g.set_field_parameter("center", c) g.set_field_parameter("bulk_velocity", bv) g.set_field_parameter("normal",np.array([0,0,1])) NP = 256 resolution = 512 rr = 10 ntheta = 720 nphi = 360 field_names = [] for vi in na.mgrid[-250:250:10]: fname = "VelocityTransmittedHI%+03i%+03i" % (vi, vi+10) print fname def _get(vel): @derived_field(name=fname) def VelocityTransmittedHI(field, data): return data["HI_Density"] * ( (((vel * 1e5) < data["RadialVelocity"]) &(data["RadialVelocity"] < (vel + 10) * 1e5))).astype("float64") field_names.append(fname) _get(vi) for field in field_names: f = h5py.File("hi_maps_pixel.h5", "a") if field in f["/"]: continue allsky = camera.allsky_projection(halo, c, #pf, c, 100.0/pf['kpc'], NP, field, inner_radius = rr, rotation = rotation) print "Computed", field, allsky.min(), allsky.max() img, count = pixelize_healpix(NP, allsky, ntheta, nphi, na.eye(3)) print "Saving", field, img.min(), img.max() f.create_dataset("%sPixel" % field, data=img) f.flush() f.close() ------------------------------------------------------------------------------------------------ On Oct 23, 2012, at 9:28 AM, Stephen Skory wrote:
Hi Ryan,
I have a question about volume rendering. I'd like to run camera.allsky_projection on a data object with a particular shape, in this case a sphere with a cylindrical hole in the center. I was able to define such an object using the boolean operator "NOT". But feeding it to camera.allsky_projection leads to an error, with the traceback below.
I think I might have an idea on how to address this, but since I'd like to A) test the fix myself before pushing it and B) have a better understanding of what you're doing in the first place, would you mind sharing the script you're using? That would be most helpful, thanks!
-- Stephen Skory s@skory.us http://stephenskory.com/ 510.621.3687 (google voice) _______________________________________________ yt-users mailing list yt-users@lists.spacepope.org http://lists.spacepope.org/listinfo.cgi/yt-users-spacepope.org