...
 
Commits (5)
......@@ -6,7 +6,7 @@ import numpy as _np
class Crefin(object):
"""
Base class used to represent the complex refractive index of a certain
grain material possibly along a certain axis.
grain material, possibly along a certain axis.
This class implements a `__call__` method which computes the crefin on a
grid of grain size and photon energy (or equivalent).
......@@ -97,43 +97,74 @@ class Crefin(object):
interpolates the crefin on the given data.
"""
if a is None:
interpolation = _interp.interp1d(lam, n, bounds_error=bounds_error)
interpolator = _interp.interp1d(lam, n, bounds_error=bounds_error)
def f(a, lam):
return _np.broadcast_to(interpolation(lam), (len(a), len(lam)))
return _np.broadcast_to(interpolator(lam), (len(a), len(lam)))
else:
rinterpolation = _interp.interp2d(lam, a, n.real,
bounds_error=bounds_error)
iinterpolation = _interp.interp2d(lam, a, n.imag,
bounds_error=bounds_error)
# 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 done for convenience, however it is inefficient
# and might be improved in the future.
interpolator = cls._interp2d_imag_sorted(a, lam, n, bounds_error)
if bounds_error and not a_bounds_error:
def f(aa, lam):
r = _np.clip(aa, _np.amin(a), _np.amax(a))
n = rinterpolation(lam, r) + 1j*iinterpolation(lam, r)
if _np.prod(r.shape)*_np.prod(lam.shape) > 1:
rowind_unsorted = _np.argsort(r)
colind_unsorted = _np.argsort(lam)
return n[rowind_unsorted][:, colind_unsorted]
else:
return _np.atleast_2d(n)
# clip a values to prevent bounds errors
aaa = _np.clip(aa, _np.amin(a), _np.amax(a))
return interpolator(aaa, lam)
else:
def f(a, lam):
n = rinterpolation(lam, a) + 1j*iinterpolation(lam, a)
if _np.prod(a.shape)*_np.prod(lam.shape) > 1:
rowind_unsorted = _np.argsort(lam)
colind_unsorted = _np.argsort(a)
return n[rowind_unsorted][:, colind_unsorted]
else:
return _np.atleast_2d(n)
def f(aa, lam):
return interpolator(aa, lam)
return cls(f)
@classmethod
def _interp2d_imag_sorted(cls, x, y, z, bounds_error=False):
"""
2-D interpolation of complex function with input-sorted output.
Returns interpolator for R**2 -> C function. When the interpolator is
called, it sorts the output array according to how the input arrays for
`x` and `y` are sorted.
Parameters
----------
x, y : array
1-D arrays of `x` and `y` values at which function is given.
z : array
2-D array of complex values of function evaluated at `x` and
`y`. `x` corresponds to first index and `y` to second.
bounds_error : bool, optional
Whether to extrapolate out of bounds values (False) or to raise an
error (True).
Returns
-------
f : function
A function that takes two array arguments of `x` and `y` values and
returns interpolated function values in 2-D array.
Notes
-----
scipy.interpolate.interp2d is used internally here. It 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 done for convenience, however it is inefficient
and might be improved in the future.
"""
rinterpolation = _interp.interp2d(y, x, z.real,
bounds_error=bounds_error)
iinterpolation = _interp.interp2d(y, x, z.imag,
bounds_error=bounds_error)
def interpolator(xx, yy):
n = _np.atleast_2d(
rinterpolation(yy, xx) + 1j*iinterpolation(yy, xx)
)
rowind_unsorted = _np.argsort(xx)
colind_unsorted = _np.argsort(yy)
return n[rowind_unsorted][:, colind_unsorted]
return interpolator
class SGPAHCrefin(Crefin):
"""
......@@ -145,7 +176,7 @@ class SGPAHCrefin(Crefin):
with open(path) as f:
data = _np.loadtxt(f, skiprows=6)
f.seek(0)
# skip tow lines
# skip two lines
[next(f) for i in range(2)]
size = float(f.readline().split('=')[0])
return size, data
......
......@@ -24,7 +24,7 @@ class SphericalGS(GrainSpecies):
return
def volume(self, sizes):
return 4*_np.pi/3 * size**3
return 4*_np.pi/3 * sizes**3
def geomCS(self, sizes):
return _np.pi * sizes**2
......
......@@ -21,6 +21,18 @@ class Testsdists(TestCase):
D03.crefins(name)
return
@data(
'cpa',
'cpe',
'sil',
)
def test_call(self, name):
crefin = D03.crefins(name)
a = np.array([0.1])*u.micron
energies = np.array([10])*u.eV
crefin(a, energies)
return
@ddt
class Testgspecies(TestCase):
......