[yt-users] camera and boolean object
M Ryan Joung
moo at astro.columbia.edu
Tue Oct 23 07:05:12 PDT 2012
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 at skory.us
> http://stephenskory.com/
> 510.621.3687 (google voice)
> _______________________________________________
> yt-users mailing list
> yt-users at lists.spacepope.org
> http://lists.spacepope.org/listinfo.cgi/yt-users-spacepope.org
>
More information about the yt-users
mailing list