[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