Commits

John Wise committed 34c9827

Adding two rendering scripts for the Planck transfer function and one
that colors the gas by entropy.

Comments (0)

Files changed (2)

entropy_colored.py

+import os, sys
+from matplotlib.cm import get_cmap
+from yt.mods import *
+from yt.utilities.amr_kdtree.api import AMRKDTree
+import yt.visualization.volume_rendering.api as vr
+import yt.visualization.image_writer as iw
+
+serial = False
+
+fn = sys.argv[-1]
+if not os.path.exists(fn):
+    print "%s not found"
+    sys.exit()
+
+#pf = load('RD0009/RedshiftOutput0009')
+pf = load(fn)
+
+N = 800 # Pixels (1024^2)
+c = na.array([0.563764, 0.540718, 0.509558])
+initial_width = 5.0 / pf["kpc"]
+final_width = 1.0 / pf["kpc"]
+Nframes = 100
+Nrotate = 100
+#cc = c - 0.5*radius/3.0**0.5
+L = [1.0, 1.0, 1.0]
+
+Nc = 256
+colormap = "idl05"
+#colormap = "RdYlBu"
+#colormap = "kamae"
+wfield="Density"
+#wfield="DensitySquared"
+wf_log=True
+field='Entropy'
+#field='Temperature'
+
+def _DensitySquared(field, data):
+    return (data["Density"]/1.673e-24)**2
+add_field("DensitySquared", function=_DensitySquared,
+          take_log = True)
+
+def _RelativeDensity(field, data):
+    dd = na.maximum(na.minimum(data["Density"], dma), dmi)
+    return (dd/dma)**2
+add_field("RelativeDensity", function=_RelativeDensity,
+          take_log = False)
+
+#reg = pf.h.region(cc, cc-radius, cc+radius)
+#tmi,tma = pf.h.all_data().quantities['Extrema'](field)[0]
+#tmi,tma = na.log10(tmi), na.log10(tma)
+#print tmi, tma
+#tmi = max(tmi, tma-3)
+#tma = tma+1
+#tmi = tmi-1
+
+#dmi,dma = pf.h.all_data().quantities['Extrema'](wfield)[0]
+#print dmi, dma
+#dmi,dma = na.log10(dmi), na.log10(dma)
+#dmi = max(dmi, dma-2)
+#dma = dma+1
+#dmi = dmi-1
+
+# Good for density/temperature
+#tmi,tma = 1.5, 5.0
+#dmi,dma = -29.0, -25.0
+# Good for density^2/entropy
+if field == "Entropy":
+    tmi,tma = -15.0, -7.0
+elif field == "Temperature":
+    tmi,tma = 1.5, 5.0
+#dmi,dma = -5.0, 1.0
+dmi,dma = -26.0, -23.0
+print 'Density:', dmi, dma
+print 'Temperature:', tmi, tma
+#tmi,tma = -11.0,-7.0
+#dmi,dma = -2.0, +4.0
+
+# Construct transfer function
+ttf = [vr.TransferFunction((tmi, tma)),
+       vr.TransferFunction((tmi, tma)),
+       vr.TransferFunction((tmi, tma))]
+xx = na.linspace(0.0, 1.0, Nc)
+cmap = get_cmap(colormap)
+for i in range(Nc):
+    ttf[0].y[i],ttf[1].y[i],ttf[2].y[i], alpha = cmap(xx[i])
+mv = vr.MultiVariateTransferFunction()
+
+for i in range(3):
+    mv.add_field_table(ttf[i], 0, weight_table_id=3)
+    mv.link_channels(i,i)
+#dtf = vr.TransferFunction((rdmi, rdma))
+#dtf.add_gaussian(rdma, 0.2*(rdma-rdmi), 1.0)
+dtf = vr.TransferFunction((dmi, dma))
+
+# Linear
+dtf.y = 0.1 * (dtf.x - dtf.x_bounds[0]) / (dtf.x_bounds[1] - dtf.x_bounds[0])
+dtf.y = dtf.y**0.5
+
+dtfg = vr.TransferFunction((dmi, dma))
+dtfg.add_gaussian(-25.5, 0.05, 1.0)
+dtfg.add_gaussian(-25.0, 0.05, 1.0)
+dtfg.add_gaussian(-24.5, 0.05, 1.0)
+dtfg.add_gaussian(-24.0, 0.05, 1.0)
+dtfg.add_gaussian(-23.5, 0.05, 1.0)
+dtfg.add_gaussian(-23.0, 0.05, 1.0)
+dtf.y += dtfg.y
+dtf.y = na.minimum(dtf.y, 1.0)
+
+# 1/x
+#alpha0 = 0.05
+#dtf.y = 1.0 / ((dtf.x_bounds[1] - dtf.x) / (dtf.x_bounds[1] - dtf.x_bounds[0]) /alpha0 + 1)
+
+#dtf.add_gaussian(dma, 1.0*(dma-dmi), 0.5)
+#dtf.add_gaussian(dma, 0.2*(dma-dmi), 1.0)
+#Ng = 8
+#width = 0.01 * (dma-dmi)/Ng
+#vv = na.mgrid[dmi:dma:Ng*1j]
+#alpha = na.linspace(0.001,1,Ng)
+#for i in range(Ng):
+#    dtf.add_gaussian(vv[i], width, alpha[i])
+#dtf.plot('test.png')
+#sys.exit()
+mv.add_field_table(dtf, 1)
+mv.link_channels(4, [3,4,5])
+
+
+# Create the camera object. Use the keyword: no_ghost=True if a lot of time is
+# spent creating vertex-centered data. In this case I'm running with 8
+# processors, and am splitting the image plane into 4 pieces and using 2
+# processors on each piece.
+cam = vr.Camera(c, L, initial_width, N,
+                transfer_function = mv,
+                fields = [field, wfield],
+                log_fields = [True, wf_log],
+                volume = None,
+                no_ghost = False,
+                sub_samples = 5, 
+                pf=pf,
+                le=c-0.5*initial_width,
+                re=c+0.5*initial_width)
+
+
+if os.path.exists("%s_kd_bricks.h5" % pf):
+    cam.volume.load_kd_bricks("%s.bricks" % pf)
+else:
+    cam.snapshot("first-vr.png")
+    cam.volume.store_kd_bricks(fn="%s.bricks" % pf)
+
+cam.snapshot("entropy-%s-second-vr.png" % pf)
+sys.exit()
+
+frame = 0
+for i, snapshot in enumerate(cam.zoomin(initial_width/final_width, Nframes)):
+    write_bitmap(snapshot, "entropy-frames/test_%04i.png" % frame)
+    frame += 1
+
+for i, snapshot in enumerate(cam.rotation(2*na.pi, Nrotate)):
+    write_bitmap(snapshot, "entropy-frames/test_%04i.png" % frame)
+    frame += 1
+    
+import os, sys
+from matplotlib.cm import get_cmap
+from yt.mods import *
+from yt.utilities.amr_kdtree.api import AMRKDTree
+import yt.visualization.volume_rendering.api as vr
+import yt.visualization.image_writer as iw
+
+serial = False
+
+fn = sys.argv[-1]
+if not os.path.exists(fn):
+    print "%s not found"
+    sys.exit()
+
+#pf = load('RD0009/RedshiftOutput0009')
+pf = load(fn)
+
+N = 800 # Pixels (1024^2)
+#c = na.array([0.5,0.5,0.5])
+c = na.array([0.563764, 0.540718, 0.509558])
+initial_width = 5.0 / pf["kpc"]
+final_width = 1.0 / pf["kpc"]
+Nframes = 100
+Nrotate = 100
+#cc = c - 0.5*radius/3.0**0.5
+L = [1.0, 1.0, 1.0]
+
+tmi,tma = 30, 1e7
+ltmi, ltma = na.log10(tmi), na.log10(tma)
+dmi,dma = 3e-29, 6e-26
+
+def _RelativeDensity(field, data):
+    dd = na.maximum(na.minimum(data["Density"], dma), dmi)
+    return (dd/dma)
+add_field("RelativeDensity", function=_RelativeDensity,
+          take_log = False)
+def _ClipTemperature(field, data):
+    return na.maximum(na.minimum(data["Temperature"], tma), tmi)
+add_field("ClipTemperature", function=_ClipTemperature,
+          take_log = True)
+
+planck = vr.PlanckTransferFunction( (ltmi-0.1, ltma+0.1), (dmi, dma), 
+                                    nbins=1024, anorm=1e6 )
+
+# Create the camera object. Use the keyword: no_ghost=True if a lot of time is
+# spent creating vertex-centered data. In this case I'm running with 8
+# processors, and am splitting the image plane into 4 pieces and using 2
+# processors on each piece.
+cam = vr.Camera(c, L, initial_width, N,
+                transfer_function = planck,
+                fields = ["ClipTemperature", "RelativeDensity"],
+                log_fields = [True, False],
+                volume = None,
+                no_ghost = False,
+                sub_samples = 5, 
+                pf=pf,
+                le=c-0.5*initial_width,
+                re=c+0.5*initial_width)
+
+
+if os.path.exists("%s_kd_bricks.h5" % pf):
+    cam.volume.load_kd_bricks()
+else:
+    cam.snapshot("first-vr.png")
+    cam.volume.store_kd_bricks()
+
+#cam.snapshot("preal-%s-second-vr.png" % pf)
+sys.exit()
+frame = 0
+for i, snapshot in enumerate(cam.zoomin(initial_width/final_width, Nframes)):
+    write_bitmap(snapshot, "frames/test_%04i.png" % frame)
+    frame += 1
+
+for i, snapshot in enumerate(cam.rotation(2*na.pi, Nrotate)):
+    write_bitmap(snapshot, "frames/test_%04i.png" % frame)
+    frame += 1