D03 crefins implemented.

cosmic_dustbox/crefin.py

 from scipy import interpolate as _interp
 import astropy.units as _u
 import numpy as _np
+import os as _os
 
 
 class Crefin(object):
@@ -11,25 +12,128 @@ class Crefin(object):
 
     def __call__(self, a, wave):
         return self.f(
-            a.to(_u.micron, equivalencies=_u.spectral()),
-            wave.to(_u.micron, equivalencies=_u.spectral()))
+            a.to(_u.micron, equivalencies=_u.spectral()).value,
+            wave.to(_u.micron, equivalencies=_u.spectral()).value)
 
     @classmethod
-    def fromData(cls, data, a_bounds_error=False, lam_bounds_error=True):
-        if len(data['a']) == 1:
-            interpolation = _interp.interp1d(
-                data['lam'], data['n'], bounds_error=lam_bounds_error)
+    def fromData(cls, a, lam, n, bounds_error=True, a_bounds_error=False):
+        """
+        Interpolate complex refractive index from given data.
 
-            def f(a, lam):
-                n = interpolation(lam)
-                return _np.broadcast_to(n, (len(a), len(lam)))
-        elif len(data['a']) == 2:
-            interpolation = _interp.interp2d(
-                data['a'], data['lam'], data['n'], bounds_error=True)
+        Parameters
+        ----------
+        a, lam : array
+            1-D arrays of grain sizes and wavelengths in micron. Must have
+            length >= 2. `a` can also be None if there is no size dependence or
+            it is not known
+
+        n : array
+            Array of the complex refractive index. Note: if `n` is 2-D, its
+            first index is wavelength and its second index is grain size. For
+            more details see the documentation of scipy.interpolate.interp2d.
+
+        bounds_error : bool, optional
+            Whether to raise out of bounds errors when interpolating.
+
+        a_bounds_error : bool, optional
+            If `bounds_error` is True, this can be used to specify if bounds
+            errors on grain size should be ignored by setting it False. This is
+            useful because often `n` is not known at many grain sizes. If
+            `bounds_error` is False or if `a` is None, this is ignored.
+
+        Note
+        ----
+        scipy.interpolate.interp2d implicitly sorts its `x` and `y` arguments
+        and returns an array that's sorted accordingly. Here we shuffle the
+        rows and columns of the returned array around so that their
+        row-(column-)order corresponds to the order of the input arrays. This
+        is highly inefficient but convenient. Like this the interpolation
+        behaves exactly as a 2D function evaluated on a meshgrid of `x` and `y`
+        would.
+        """
+        if a is None:
+            interpolation = _interp.interp1d(lam, n, bounds_error=bounds_error)
 
             def f(a, lam):
-                r = _np.clip(a, data['a'][0], data['a'][1])
-                return interpolation(a, lam)
+                return _np.broadcast_to(interpolation(lam), (len(a), len(lam)))
         else:
-            raise ValueError("Currently no method implemented!")
+            rinterpolation = _interp.interp2d(a, lam, n.real,
+                bounds_error=bounds_error)
+            iinterpolation = _interp.interp2d(a, lam, n.imag,
+                bounds_error=bounds_error)
+
+            if bounds_error and not a_bounds_error:
+                def f(aa, lam):
+                    r = _np.clip(aa, _np.amin(a), _np.amax(a))
+                    xind = _np.argsort(r)
+                    yind = _np.argsort(lam)
+                    return (
+                        rinterpolation(r, lam)
+                        + 1j*iinterpolation(r, lam))[yind][:, xind].T
+            else:
+                def f(a, lam):
+                    xind = _np.argsort(r)
+                    yind = _np.argsort(lam)
+                    return (
+                        rinterpolation(a, lam)
+                        + 1j*iinterpolation(a, lam))[yind][:, xind].T
         return cls(f)
+
+
+class SGPAHCrefin(Crefin):
+
+    @classmethod
+    def parseCrefinFile(cls, path):
+        with open(path) as f:
+            data = _np.loadtxt(f, skiprows=5)
+            f.seek(0)
+            next(f)
+            size = float(f.readline().split('=')[0])
+        return size, data
+
+    @classmethod
+    def fromFiles(cls, paths):
+        a = []
+        n = []
+        for j, path in enumerate(paths):
+            size, data = cls.parseCrefinFile(path)
+            if size in a:
+                raise ValueError("Grain size "
+                                 + str(size) + "provided twice.")
+            else:
+                a.append(size)
+            if j == 0:
+                lam = data[:, 0]
+            else:
+                if not _np.array_equal(lam, data[:, 0]):
+                    raise ValueError(
+                        "Wavelengths in file " + path
+                        + "differ from those in the other files visited "
+                        "so far.")
+            nt = data[:, 3] + 1 + 1j*data[:, 4]
+            n.append(nt)
+        if len(a) > 1:
+            a = _np.array(a)
+            n = _np.array(n).T
+        elif len(a) == 1:
+            a = None
+            n = _np.array(n)
+        return cls.fromData(a, lam, n, bounds_error=True, a_bounds_error=False)
+
+
+cpaD03 = SGPAHCrefin.fromFiles([
+    _os.path.join(_os.path.dirname(__file__),
+                  'data/crefin/callindex.out_CpaD03_0.01'),
+    _os.path.join(_os.path.dirname(__file__),
+                  'data/crefin/callindex.out_CpaD03_0.10')
+])
+cpeD03 = SGPAHCrefin.fromFiles([
+    _os.path.join(_os.path.dirname(__file__),
+                  'data/crefin/callindex.out_CpeD03_0.01'),
+    _os.path.join(_os.path.dirname(__file__),
+                  'data/crefin/callindex.out_CpeD03_0.10')
+])
+silD03 = SGPAHCrefin.fromFiles([
+    _os.path.join(_os.path.dirname(__file__),
+                  'data/crefin/callindex.out_silD03')
+])

cosmic_dustbox/tests/test_crefin.py

@@ -2,28 +2,198 @@ from unittest import TestCase
 from .. import crefin
 import numpy as np
 import astropy.units as u
+import os
 
 
 class TestFromData(TestCase):
 
     @classmethod
     def setUpClass(cls):
-        cls.data = {
-            'a': [1e-9],
-            'lam': np.linspace(1, 10, num=10),
-            'n': np.random.random(10) + 1j * np.random.random(10),
-        }
+        cls.a1 = None
+        cls.lam1 = np.linspace(1, 10, num=10)
+        cls.n1= np.random.random(10) + 1j * np.random.random(10)
+
+        cls.a2 = np.array([1e-9, 1e-8])
+        cls.lam2 = np.linspace(1, 10, num=10)
+        cls.n2 = (
+            np.random.random(20) + 1j * np.random.random(20)).reshape(10, 2)
         return
 
     @classmethod
     def tearDownClass(cls):
         return
 
-    def test_smoke_test(self):
-        a = crefin.Crefin.fromData(self.__class__.data)
+    def test_no_a(self):
+        c = self.__class__
+        a = crefin.Crefin.fromData(c.a1, c.lam1, c.n1)
+        return
+
+    def test_normal(self):
+        c = self.__class__
+        a = crefin.Crefin.fromData(c.a2, c.lam2, c.n2)
+        return
+
+    def test_normal_no_bounds(self):
+        c = self.__class__
+        a = crefin.Crefin.fromData(c.a2, c.lam2, c.n2, bounds_error=False)
+        return
+
+
+class TestSGPAHCrefin(TestCase):
+
+    @classmethod
+    def setUpClass(cls):
+        cls.basepath = os.path.join(
+            os.path.dirname(os.path.abspath(__file__)), '../data/crefin')
+        return
+
+    @classmethod
+    def tearDownClass(cls):
+        return
+
+    def test_parseCrefinFile_asil(self):
+        path = os.path.join(self.__class__.basepath, 'callindex.out_silD03')
+        size, data = crefin.SGPAHCrefin.parseCrefinFile(path)
+        self.assertAlmostEqual(0.1, size)
+        self.assertEqual(data.shape[0], 837)
+        self.assertEqual(data.shape[1], 5)
+        return
+
+    def test_parseCrefinFile_cpar(self):
+        path = os.path.join(self.__class__.basepath,
+                            'callindex.out_CpaD03_0.01')
+        size, data = crefin.SGPAHCrefin.parseCrefinFile(path)
+        self.assertAlmostEqual(0.01, size)
+        self.assertEqual(data.shape[0], 386)
+        self.assertEqual(data.shape[1], 5)
+        return
+
+    def test_parseCrefinFile_cperp(self):
+        path = os.path.join(self.__class__.basepath,
+                            'callindex.out_CpeD03_0.01')
+        size, data = crefin.SGPAHCrefin.parseCrefinFile(path)
+        self.assertAlmostEqual(0.01, size)
+        self.assertEqual(data.shape[0], 383)
+        self.assertEqual(data.shape[1], 5)
+        return
+
+    def test_fromFiles_asil(self):
+        path = os.path.join(self.__class__.basepath, 'callindex.out_silD03')
+        c = crefin.SGPAHCrefin.fromFiles([path])
+        a = np.array([0.1])*u.micron
+        wave = np.array([1.239840e5, 5.481167e2, 6.199200e-5])*u.micron
+        nre_expect = [2.4350+1, 2.4040+1, -1.9450e-6+1]
+        nim_expect = [1.1190e-3, 8.2550e-2, 1.7830e-8]
+        n = c(a, wave)
+
+        self.assertEqual(n.shape[0], len(a))
+        self.assertEqual(n.shape[1], len(wave))
+        for j, r in enumerate(n[0, :]):
+            self.assertAlmostEqual(r.real, nre_expect[j])
+            self.assertAlmostEqual(r.imag, nim_expect[j])
+        return
+
+    def test_fromFiles_asil_multiple_a(self):
+        path = os.path.join(self.__class__.basepath, 'callindex.out_silD03')
+        c = crefin.SGPAHCrefin.fromFiles([path])
+        a = np.array([1.0, 0.1, 0.01])*u.micron
+        wave = np.array([1.239840e5, 5.481167e2, 6.199200e-5])*u.micron
+        nre_expect = [2.4350+1, 2.4040+1, -1.9450e-6+1]
+        nim_expect = [1.1190e-3, 8.2550e-2, 1.7830e-8]
+        n = c(a, wave)
+
+        self.assertEqual(n.shape[0], len(a))
+        self.assertEqual(n.shape[1], len(wave))
+        for i in range(len(a)):
+            for j, r in enumerate(n[i, :]):
+                self.assertAlmostEqual(r.real, nre_expect[j])
+                self.assertAlmostEqual(r.imag, nim_expect[j])
+        return
+
+    def test_fromFiles_cperp(self):
+        paths = [
+            os.path.join(self.__class__.basepath, 'callindex.out_CpeD03_0.01'),
+            os.path.join(self.__class__.basepath, 'callindex.out_CpeD03_0.10')
+        ]
+        c = crefin.SGPAHCrefin.fromFiles(paths)
+        a = np.array([1.0, 0.1, 0.01, 0.001])*u.micron
+        wave = np.array([1.239840E+05, 6.322489E+00, 6.199200E-05])*u.micron
+        nre_expect = np.array([
+            [1.0596E+03, 4.5430E+00, 1.5405E-10],
+            [1.0596E+03, 4.5430E+00, 1.5405E-10],
+            [3.4137E+02, 4.5568E+00, 1.5405E-10],
+            [3.4137E+02, 4.5568E+00, 1.5405E-10],
+        ])+1
+        nim_expect = [
+            [1.0636E+03, 4.2997E+00, 1.4471E-17],
+            [1.0636E+03, 4.2997E+00, 1.4471E-17],
+            [3.4149E+02, 4.3122E+00, 1.3963E-16],
+            [3.4149E+02, 4.3122E+00, 1.3963E-16],
+        ]
+        n = c(a, wave)
+
+        self.assertEqual(n.shape[0], len(a))
+        self.assertEqual(n.shape[1], len(wave))
+        for i in [0, 1]:
+            for j, r in enumerate(n[i, :]):
+                self.assertAlmostEqual(r.real, nre_expect[i][j],
+                                       msg='Failed at i='+str(i)+' j='+str(j))
+                self.assertAlmostEqual(r.imag, nim_expect[i][j])
+        return
+
+    def test_fromFiles_cpar(self):
+        paths = [
+            os.path.join(self.__class__.basepath, 'callindex.out_CpaD03_0.01'),
+            os.path.join(self.__class__.basepath, 'callindex.out_CpaD03_0.10')
+        ]
+        c = crefin.SGPAHCrefin.fromFiles(paths)
+        a = np.array([1.0, 0.1, 0.01, 0.001])*u.micron
+        wave = np.array([1.239840E+05, 6.199200E+00, 6.199200E-05])*u.micron
+        nre_expect = np.array([
+            [1.5051E+02, 1.1266E+00, 2.3165E-10],
+            [1.5051E+02, 1.1266E+00, 2.3165E-10],
+            [1.4329E+02, 1.0261E+00, 2.3165E-10],
+            [1.4329E+02, 1.0261E+00, 2.3165E-10],
+        ])+1
+        nim_expect = [
+            [1.5156E+02, 2.3038E-02, 1.4161E-17],
+            [1.5156E+02, 2.3038E-02, 1.4161E-17],
+            [1.4434E+02, 3.0850E+00, 1.5614E-17],
+            [1.4434E+02, 3.0850E+00, 1.5614E-17],
+        ]
+        n = c(a, wave)
+
+        self.assertEqual(n.shape[0], len(a))
+        self.assertEqual(n.shape[1], len(wave))
+        for i in [0, 1]:
+            for j, r in enumerate(n[i, :]):
+                self.assertAlmostEqual(r.real, nre_expect[i][j],
+                                       msg='Failed at i='+str(i)+' j='+str(j))
+                self.assertAlmostEqual(r.imag, nim_expect[i][j])
+        return
+
+    def test_fromFiles_cpar(self):
+        paths = [
+            os.path.join(self.__class__.basepath, 'callindex.out_CpaD03_0.01'),
+            os.path.join(self.__class__.basepath, 'callindex.out_CpaD03_0.10')
+        ]
+        a = crefin.SGPAHCrefin.fromFiles(paths)
+        return
+
+    def test_fromFiles_diff_wavelengths(self):
+        paths = [
+            os.path.join(self.__class__.basepath, 'callindex.out_CpeD03_0.01'),
+            os.path.join(self.__class__.basepath, 'callindex.out_CpaD03_0.10')
+        ]
+        with self.assertRaises(ValueError) as context:
+            a = crefin.SGPAHCrefin.fromFiles(paths)
         return
 
-    def test_interp(self):
-        a = crefin.Crefin.fromData(self.__class__.data)
-        a(np.array([1.0])*u.micron, np.array([1.0])*u.micron)
+    def test_fromFiles_same_size(self):
+        paths = [
+            os.path.join(self.__class__.basepath, 'callindex.out_CpeD03_0.01'),
+            os.path.join(self.__class__.basepath, 'callindex.out_CpeD03_0.01')
+        ]
+        with self.assertRaises(ValueError) as context:
+            a = crefin.SGPAHCrefin.fromFiles(paths)
         return