[Yt-svn] yt: Adding "scattering" for absorption in the Planck function. ...

[hg at spacepope.org]

hg Repository: yt
details:   yt/rev/3da92a5ea545
changeset: 1515:3da92a5ea545
user:      Matthew Turk <matthewturk at gmail.com>
date:
Fri Apr 02 08:02:09 2010 -0700
description:
Adding "scattering" for absorption in the Planck function.  Really, it's a huge
hack and because basically everything is done unit free it doesn't make that
much sense just now.  Maybe we'll fix it some day to be unitful, and thus,
at least self-consistent.

diffstat:

 yt/_amr_utils/VolumeIntegrator.pyx                 |   4 +-
 yt/extensions/volume_rendering/TransferFunction.py |  42 ++++++++++++---------
 yt/extensions/volume_rendering/UBVRI.py            |   5 ++
 3 files changed, 32 insertions(+), 19 deletions(-)

diffs (94 lines):

diff -r fab318774483 -r 3da92a5ea545 yt/_amr_utils/VolumeIntegrator.pyx
--- a/yt/_amr_utils/VolumeIntegrator.pyx	Thu Apr 01 22:30:30 2010 -0700
+++ b/yt/_amr_utils/VolumeIntegrator.pyx	Fri Apr 02 08:02:09 2010 -0700
@@ -195,6 +195,8 @@
             if fid != -1: self.istorage[i] *= self.istorage[fid]
         for i in range(6):
             trgba[i] = self.istorage[self.field_table_ids[i]]
+            #print i, trgba[i],
+        #print
         # A few words on opacity.  We're going to be integrating equation 1.23
         # from Rybicki & Lightman.  dI_\nu / ds = -\alpha_\nu I_\nu + j_\nu
         # \alpha_nu = \kappa \rho , but we leave that up to the input
@@ -206,7 +208,7 @@
             # This is the new way: alpha corresponds to opacity of a given
             # slice.  Previously it was ill-defined, but represented some
             # measure of emissivity.
-            ta = (1.0 - dt*trgba[i+3])
+            ta = fmax((1.0 - dt*trgba[i+3]), 0.0)
             rgba[i  ] = dt*trgba[i  ] + ta * rgba[i  ]
             rgba[i+3] = dt*trgba[i+3] + ta * rgba[i+3]
             # This is the old way:
diff -r fab318774483 -r 3da92a5ea545 yt/extensions/volume_rendering/TransferFunction.py
--- a/yt/extensions/volume_rendering/TransferFunction.py	Thu Apr 01 22:30:30 2010 -0700
+++ b/yt/extensions/volume_rendering/TransferFunction.py	Fri Apr 02 08:02:09 2010 -0700
@@ -179,8 +179,22 @@
 
 class PlanckTransferFunction(MultiVariateTransferFunction):
     def __init__(self, T_bounds, rho_bounds, nbins=256,
-                 red='R', green='V', blue='B'):
+                 red='R', green='V', blue='B',
+                 anorm = 1e6):
+        """
+        This sets up a planck function for multivariate emission and
+        absorption.  We assume that the emission is black body, which is then
+        convolved with appropriate Johnson filters for *red*, *green* and
+        *blue*.  *T_bounds* and *rho_bounds* define the limits of tabulated
+        emission and absorption functions.  *anorm* is a "fudge factor" that
+        defines the somewhat arbitrary normalization to the scattering
+        approximation: because everything is done largely unit-free, and is
+        really not terribly accurate anyway, feel free to adjust this to change
+        the relative amount of reddenning.  Maybe in some future version this
+        will be unitful.
+        """
         MultiVariateTransferFunction.__init__(self)
+        mscat = -1
         from UBVRI import johnson_filters
         for i, f in enumerate([red, green, blue]):
             jf = johnson_filters[f]
@@ -188,24 +202,16 @@
             tf.add_filtered_planck(jf['wavelen'], jf['trans'])
             self.add_field_table(tf, 0, 1)
             self.link_channels(i, i) # 0 => 0, 1 => 1, 2 => 2
+            mscat = max(mscat, jf["Lchar"]**-4)
 
-        # These are the alpha: r_alpha
-        tf = TransferFunction(rho_bounds)
-        tf.add_line( (rho_bounds[0], 0.1), (rho_bounds[1], 0.1) )
-        self.add_field_table(tf, 1)
-        self.link_channels(3, 3)
-
-        # These are the alpha: g_alpha
-        tf = TransferFunction(rho_bounds)
-        tf.add_line( (rho_bounds[0], 0.1), (rho_bounds[1], 0.1) )
-        self.add_field_table(tf, 1)
-        self.link_channels(4, 4)
-
-        # These are the alpha: b_alpha
-        tf = TransferFunction(rho_bounds)
-        tf.add_line( (rho_bounds[0], 0.1), (rho_bounds[1], 0.1) )
-        self.add_field_table(tf, 1)
-        self.link_channels(5, 5)
+        for i, f in enumerate([red, green, blue]):
+            # Now we set up the scattering
+            scat = (johnson_filters[f]["Lchar"]**-4 / mscat)*1e6
+            tf = TransferFunction(rho_bounds)
+            print "Adding: %s with relative scattering %s" % (f, scat)
+            tf.y *= 0.0; tf.y += scat
+            self.add_field_table(tf, 1, weight_field_id = 1)
+            self.link_channels(i+3, i+3)
 
         self._normalize()
 
diff -r fab318774483 -r 3da92a5ea545 yt/extensions/volume_rendering/UBVRI.py
--- a/yt/extensions/volume_rendering/UBVRI.py	Thu Apr 01 22:30:30 2010 -0700
+++ b/yt/extensions/volume_rendering/UBVRI.py	Fri Apr 02 08:02:09 2010 -0700
@@ -98,3 +98,8 @@
         dtype='float64'),
       ),
     )
+
+for filter, vals in johnson_filters.items():
+    wavelen = vals["wavelen"]
+    trans = vals["trans"]
+    vals["Lchar"] = wavelen[na.argmax(trans)]