# This program is public domain
# Author: Paul Kienzle
r"""
Scattering length density profile.
In order to render a reflectometry model, the theory function calculator
renders each layer in the model for each energy in the probe. For slab
layers this is easy: just accumulate the slabs, with the 1-\ $\sigma$ Gaussian
interface width between the slabs. For freeform or functional layers,
this is more complicated. The rendering needs to chop each layer into
microslabs and evaluate the profile at each of these slabs.
Example
-------
This example sets up a model which uses tanh to transition from
silicon to gold in 20 |Ang| with 2 |Ang| steps.
First define the profile, and put in the substrate:
>>> S = Microslabs(nprobe=1, dz=2)
>>> S.clear()
>>> S.append(w=0, rho=2.07)
Next add the interface. This uses :meth:`microslabs` to select
the points at which the interface is evaluated, much like you
would do when defining your own special layer type. Note that the
points Pz are in the center of the micro slabs. The width of the
final slab may be different. You do not need to use fixed width
microslabs if you can more efficiently represent the profile with
a smaller number of variable width slabs, but :meth:`contract_profile`
serves the same purpose with less work on your part.
>>> from numpy import tanh
>>> Pw, Pz = S.microslabs(20)
>>> print("widths = %s ..."%(" ".join("%g"%v for v in Pw[:5])))
widths = 2 2 2 2 2 ...
>>> print("centers = %s ..."%(" ".join("%g"%v for v in Pz[:5])))
centers = 1 3 5 7 9 ...
>>> rho = (1-tanh((Pz-10)/5))/2*(2.07-4.5)+4.5
>>> S.extend(w=Pw, rho=[rho])
Finally, add the incident medium and see the results. Note that *rho*
is a matrix, with one column for each incident energy. We are only
using one energy so we only show the first column.
>>> S.append(w=0, rho=4.5)
>>> print("width = %s ..."%(" ".join("%g"%v for v in S.w[:5])))
width = 0 2 2 2 2 ...
>>> print("rho = %s ..."%(" ".join("%.2f"%v for v in S.rho[0, :5])))
rho = 2.07 2.13 2.21 2.36 2.63 ...
Since *irho* and *sigma* were not specified, they will be zero.
>>> print("sigma = %s ..."%(" ".join("%g"%v for v in S.sigma[:5])))
sigma = 0 0 0 0 0 ...
>>> print("irho = %s ..."%(" ".join("%g"%v for v in S.irho[0, :5])))
irho = 0 0 0 0 0 ...
"""
from __future__ import division, print_function
import numpy as np
from numpy import inf, nan, isnan
from scipy.special import erf
from .reflectivity import BASE_GUIDE_ANGLE as DEFAULT_THETA_M
[docs]
class Microslabs(object):
"""
Manage the micro slab representation of a model.
In order to compute reflectivity, we need a series of slabs with thickness,
roughness and scattering potential for each slab. Because scattering
potentials are probe dependent we store an array of potentials for each
probe value.
Some slab models use non-uniform layers, and so need the additional
parameter of dz for the step size within the layer.
The space for the slabs is saved even after reset, in preparation for a
new set of slabs from different fitting parameters.
"""
def __init__(self, nprobe, dz=1):
self._num_slabs = 0
# _slabs contains the 1D objects w, sigma of len n
# _slabs_rho contains 2D objects rho, irho, with one for each wavelength
# rhoM, thetaM will contain magnetic moment and angle
self._slabs = np.empty(shape=(0, 2))
self._slabs_rho = np.empty(shape=(0, nprobe, 2))
self.rhoM = None # type: np.ndarray
self.thetaM = None # type: np.ndarray
self._slabs_mag = np.empty(shape=(0, nprobe, 2))
self.dz = dz
self._magnetic_sections = []
self._z_left = self._z_right = 0.
self._z_offset = 0.
[docs]
def microslabs(self, thickness=0):
"""
Return a set of microslabs for a layer of the given *thickness*.
The step size slabs.dz was defined when the Microslabs
object was created.
This is a convenience function. Layer definitions can choose
their own slices so long as the step size is approximately
slabs.dz in the varying region.
:Parameters:
*thickness* : float | A
Layer thickness
:Returns:
*widths*: vector | A
Microslab widths
*centers*: vector | A
Microslab centers
"""
# TODO: force dz onto a common boundary to avoid remeshing
# in the smooth profile function
edges = np.arange(0, thickness + self.dz, self.dz, dtype='d')
edges[-1] = thickness
centers = (edges[1:] + edges[:-1]) / 2
widths = edges[1:] - edges[:-1]
return widths, centers
[docs]
def clear(self):
"""
Reset the slab model so that none are present.
"""
self._num_slabs = 0
self._magnetic_sections = []
def __len__(self):
return self._num_slabs
[docs]
def repeat(self, start=0, count=1, interface=0):
"""
Extend the model so that there are *count* versions of the slabs
from *start* to the final slab.
This is equivalent to L.extend(L[start:]*(count-1)) for list L.
"""
# For now use the dumb implementation; a better implementation
# would remember the repeats and pre-calculate the matrix product
# for the repeating region, saving much work later. This has
# to work in conjunction with interfaces and with magnetic profiles.
repeats = count - 1
end = len(self)
length = end - start
fromidx = slice(start, end)
toidx = slice(end, end + repeats * length)
self._reserve(repeats * length)
self._slabs[toidx] = np.tile(self._slabs[fromidx], [repeats, 1])
self._slabs_rho[toidx] = np.tile(self._slabs_rho[fromidx], [repeats, 1, 1])
self._num_slabs += repeats * length
# Replace interface on the top
self._slabs[self._num_slabs - 1, 1] = interface
if self._magnetic_sections:
raise NotImplementedError("Repeated magnetic layers not implemented")
def _reserve(self, nadd):
"""
Reserve space for at least *nadd* slabs.
"""
ns, nl, _ = self._slabs_rho.shape
if ns < self._num_slabs + nadd:
new_ns = self._num_slabs + nadd + 50
self._slabs = self._slabs.copy()
self._slabs.resize((new_ns, 4), refcheck=False)
self._slabs_rho = self._slabs_rho.copy()
self._slabs_rho.resize((new_ns, nl, 2), refcheck=False)
[docs]
def extend(self, w=0, sigma=0, rho=0, irho=0):
"""
Extend the micro slab model with the given layers.
"""
nadd = len(w)
self._reserve(nadd)
idx = slice(self._num_slabs, self._num_slabs + nadd)
self._num_slabs += nadd
self._slabs[idx, 0] = w
self._slabs[idx, 1] = sigma
self._slabs_rho[idx, :, 0] = np.asarray(rho).T
self._slabs_rho[idx, :, 1] = np.asarray(irho).T
[docs]
def append(self, w=0, sigma=0, rho=0, irho=0):
"""
Extend the micro slab model with a single layer.
"""
#self.extend(w=[w], sigma=[sigma], rho=[rho], irho=[irho])
#return
self._reserve(1)
self._slabs[self._num_slabs, 0] = w
self._slabs[self._num_slabs, 1] = sigma
self._slabs_rho[self._num_slabs, :, 0] = rho
self._slabs_rho[self._num_slabs, :, 1] = irho
self._num_slabs += 1
[docs]
def add_magnetism(self, anchor, w, rhoM=0, thetaM=DEFAULT_THETA_M, sigma=0):
"""
Add magnetic layers.
Note that *sigma* is a pair *(interface_below, interface_above)*
representing the magnetic roughness, which may be different
from the nuclear roughness at the layer boundaries.
"""
w = np.asarray(w, 'd')
if np.isscalar(sigma):
sigma = (sigma, sigma)
self._magnetic_sections.append((np.vstack((w, rhoM, thetaM)),
anchor, sigma))
[docs]
def thickness(self):
"""
Total thickness of the profile.
Note that thickness includes the thickness of the substrate and
surface layers. Normally these will be zero, but the contract
profile operation may result in large values for either.
"""
return np.sum(self._slabs[1:self._num_slabs, 0])
@property
def w(self):
"Thickness (A)"
return self._slabs[:self._num_slabs, 0]
@property
def sigma(self):
"rms roughness (A)"
return self._slabs[:self._num_slabs - 1, 1]
@property
def surface_sigma(self):
"roughness for the current top layer, or nan if substrate"
return self._slabs[self._num_slabs - 1, 1] if self._num_slabs > 0 else nan
@property
def rho(self):
"Scattering length density (10^-6 number density)"
return self._slabs_rho[:self._num_slabs, :, 0].T
@property
def irho(self):
"Absorption (10^-6 number density)"
return self._slabs_rho[:self._num_slabs, :, 1].T
@property
def ismagnetic(self):
"True if there are magnetic materials in any slab"
return self._magnetic_sections != []
[docs]
def limited_sigma(self, limit=0):
"""
Limit the roughness by some fraction of layer thickness.
This function should be called before :meth:`finalize`, but after
all slabs have been added to the profile.
*limit* is the number of times sigma has to fit in the layers
on either side of the interface. The returned sigma is
truncated to min(wlow, whigh)/*limit* where wlow is the thickness
of the layer below the interface, and whigh is the thickness above
the interface. A *limit* value of 0 returns the original sigma.
Although a gaussian inteface extends to infinity, in practice
setting a *limit* of 3 allows the layer to reach its bulk value,
with no cross talk between the interfaces. For very large
roughnesses, the blending algorithm allows the sld beyond
the interface to bleed through the entire layer and into the
next. In this case the roughness should be the same on both
sides of the layer to avoid artifacts at the interface.
Magnetic roughness is ignored for now.
"""
self.sigma[:] = compute_limited_sigma(self.w, self.sigma, limit)
[docs]
def finalize(self, step_interfaces, dA):
"""
Rendering complete.
Call this method after the microslab model has been constructed,
so any post-rendering processes can be completed.
In addition to clearing any width from the substrate and
the surface surround, this will align magnetic and nuclear slabs,
convert interfaces to step interfaces if desired, and merge slabs
with similar scattering potentials to reduce computation time.
*step_interfaces* is True if interfaces should be rendered using
slabs.
*dA* is the tolerance to use when deciding if similar layers can
be merged.
"""
if self.ismagnetic:
self._align_magnetic_and_nuclear()
self._set_z_range()
# render step interfaces
if step_interfaces:
self._render_interfaces()
if self.ismagnetic:
self._contract_magnetic(dA)
else:
self._contract_profile(dA)
def _set_z_range(self):
"""
Make sure z-range includes 3-sigma around every interface.
"""
self.w[0] = self.w[-1] = 0
offset = np.cumsum(self.w[:-1])
self._z_left = min(-10, np.min(offset - 3*self.sigma))
self._z_right = max(offset[-1]+10, np.max(offset + 3*self.sigma))
def _align_magnetic_and_nuclear(self):
"""
Add magnetic information to the nuclear slabs, introducing new
slabs as necessary where magnetic and nuclear do not match.
"""
from .refllib import align_magnetic
# Nuclear profile (one wavelength only)
#if self.rho.shape[0] != 1:
# raise ValueError("wavelength-dependent magnetism not supported")
w, sigma, rho, irho = self.w, self.sigma, self.rho[0], self.irho[0]
# Fill in gaps for magnetic profile
wM, sigmaM, rhoM, thetaM = self._join_magnetic_sections(gap_size=1e-6)
# Align nuclear and magnetic
w, sigma, rho, irho, wM, sigmaM, rhoM, thetaM = [
np.ascontiguousarray(v, 'd')
for v in (w, sigma, rho, irho, wM, sigmaM, rhoM, thetaM)
]
output = np.empty((len(w)+len(wM), 6), 'd')
n = align_magnetic(w, sigma, rho, irho, wM, sigmaM, rhoM, thetaM, output)
# Store the resulting profile
self._reserve(n - self._num_slabs) # make sure there is space
self._num_slabs = n
self.w[:] = output[:n, 0]
self.sigma[:] = output[:n-1, 1]
self.rho[0][:] = output[:n, 2]
self.irho[0][:] = output[:n, 3]
self.rhoM = output[:n, 4]
self.thetaM = output[:n, 5]
def _render_interfaces(self):
"""
Convert all interfaces into step interfaces by sampling the
analytic version of the smoothed profile at intervals of dz.
The interface effects are limited to the surrounding layers.
Use of contract_profile afterward is strongly recommended, for
better performance on models with large sections of constant
scattering potential.
"""
z = np.arange(self._z_left, self._z_right + 0.5*self.dz, self.dz)
n_slabs = len(z)
n_profiles = self.rho.shape[0]
offsets = np.cumsum(self.w[:-1]) # assumes w[0] == 0 in _set_z_range
# generate profiles
rho = np.empty((n_profiles, n_slabs), 'd')
irho = np.empty((n_profiles, n_slabs), 'd')
for k in range(n_profiles):
# Gd support: cycle through wavelength dependent rho/irho
rho[k] = build_profile(z, offsets, self.sigma, self.rho[k])
irho[k] = build_profile(z, offsets, self.sigma, self.irho[k])
if self.ismagnetic:
rhoM = build_profile(z, offsets, self.sigma, self.rhoM)
thetaM = build_profile(z, offsets, self.sigma, self.thetaM)
w = self.dz * np.ones(n_slabs)
w[0] = w[-1] = 0.
# update slabs
self._reserve(n_slabs - self._num_slabs)
self._num_slabs = n_slabs
self.w[:] = w
self.sigma[:] = 0
self.rho[:,:] = rho
self.irho[:,:] = irho
if self.ismagnetic:
self.rhoM = rhoM
self.thetaM = thetaM
self._z_offset = self._z_left
def _contract_profile(self, dA):
from .refllib import contract_by_area
if dA is None:
return
# TODO: need a separate implementation for multiple wavelengths
if self.rho.shape[0] > 1:
# Note: should at least check for duplicates otherwise thick
# layers will get extremely slow
return
w, sigma, rho, irho = [
np.ascontiguousarray(v, 'd')
for v in (self.w, self.sigma, self.rho[0], self.irho[0])
]
#print "final sld before contract", rho[-1]
n = contract_by_area(w, sigma, rho, irho, dA)
self._num_slabs = n
self.w[:] = w[:n]
self.rho[0, :] = rho[:n]
self.irho[0, :] = irho[:n]
self.sigma[:] = sigma[:n-1]
#print "final sld after contract", rho[n-1], self.rho[0][n-1], n
def _contract_magnetic(self, dA):
from .refllib import contract_mag
if dA is None:
return
# TODO: need a separate implementation for multiple wavelengths
if self.rho.shape[0] > 1:
# Note: should at least check for duplicates otherwise thick
# layers will get extremely slow
return
w, sigma, rho, irho, rhoM, thetaM = \
[np.ascontiguousarray(v, 'd')
for v in (self.w, self.sigma, self.rho[0], self.irho[0], self.rhoM, self.thetaM)]
#print "final sld before contract", rho[-1]
n = contract_mag(w, sigma, rho, irho, rhoM, thetaM, dA)
self._num_slabs = n
self.w[:] = w[:n]
self.rho[0][:] = rho[:n]
self.irho[0][:] = irho[:n]
self.rhoM = rhoM[:n]
self.thetaM = thetaM[:n]
self.sigma[:] = sigma[:n-1]
#self.sigma[:] = 0
#print "final sld after contract", rho[n-1], self.rho[0][n-1], n
def _DEAD_apply_smoothness(self, dA, smoothness=0.3):
"""
Set a guassian interface for layers which have been coalesced using
the contract_profile function.
Note that we guess which interfaces this applies to after the fact
using criteria similar to those used to coalesce the microslabs
into layers, and so it may apply to layers which are close in
scattering length density and have zero sigma, but which were
distinct in the original model. The displayed profile will show
the profile used to calculate the reflectivity, so even though
this behaviour is different from what the user intended, the
result will not be misleading.
In a detailed example of a tethered polymer model, smoothness
was found to be worse than no smoothness, so this function has
been removed from the execution stream.
"""
if dA is None or smoothness == 0:
return
# TODO: refine this so that it can look forward as well as back
# TODO: also need to avoid changing explicit sigma=0...
w, sigma, rho, irho = self.w, self.sigma, self.rho[0], self.irho[0]
step = (abs(np.diff(rho)) + abs(np.diff(irho)))
step[:-1] *= w[1:-1] # compute dA of step; substrate uses w=1
idx = np.nonzero(sigma == 0)[0]
fix = step[idx] < 3 * dA
sigma[idx[fix]] = w[idx[fix]] * smoothness
[docs]
def step_profile(self):
"""
Return a step profile representation of the microslab structure.
Nevot-Croce interfaces are not represented.
"""
rho = np.vstack([self.rho[0, :], self.rho[0, :]]).T.flatten()
irho = np.vstack([self.irho[0, :], self.irho[0, :]]).T.flatten()
if len(self.w) > 2:
offsets = np.cumsum(self.w[0:-1])
z = np.vstack([np.hstack([-10, offsets]),
np.hstack([offsets, offsets[-1]+10])]).T.flatten()
else:
z = np.array([-10, 0, 0, 10])
return z+self._z_offset, rho, irho
[docs]
def magnetic_step_profile(self):
"""
Return a step profile representation of the microslab structure.
Nevot-Croce interfaces are not represented.
"""
z, rho, irho = self.step_profile()
rhoM = np.vstack([self.rhoM, self.rhoM]).T.flatten()
thetaM = np.vstack([self.thetaM, self.thetaM]).T.flatten()
return z+self._z_offset, rho, irho, rhoM, thetaM
[docs]
def smooth_profile(self, dz=0.1):
"""
Return a smooth profile representation of the microslab structure
Nevot-Croce roughness is approximately represented, though the
calculation is incorrect for layers with large roughness compared
to the thickness.
The returned profile has uniform step size *dz*.
"""
z = np.arange(self._z_left, self._z_right + 0.5*dz, dz)
offsets = np.cumsum(self.w) + self._z_offset
irho = build_profile(z, offsets, self.sigma, self.irho[0])
rho = build_profile(z, offsets, self.sigma, self.rho[0])
return z, rho, irho
[docs]
def magnetic_smooth_profile(self, dz=0.1):
"""
Return a profile representation of the magnetic microslab structure.
"""
z = np.arange(self._z_left, self._z_right + 0.5*dz, dz)
offsets = np.cumsum(self.w) + self._z_offset
irho = build_profile(z, offsets, self.sigma, self.irho[0])
rho = build_profile(z, offsets, self.sigma, self.rho[0])
rhoM = build_profile(z, offsets, self.sigma, self.rhoM)
thetaM = build_profile(z, offsets, self.sigma, self.thetaM)
return z, rho, irho, rhoM, thetaM
def _join_magnetic_sections(self, gap_size):
"""
Convert anchored magnetic sections into coarse magnetic slabs.
"""
# Find the magnetic blocks
blocks, offsets, sigmas = zip(*self._magnetic_sections)
#print "blocks", blocks
#print "offsets", offsets
#print "sigmas", sigmas
# * Splice the blocks together with rhoM=0 in the gaps and
# thetaM=(thetaM_below+thetaM_above)/2 in the gaps.
# * Result is:
# slices = [(thickness, rhoM, thetaM), (thickness, rhoM, thetaM), ...]
# * Initialize slices with the magnetism of the substrate, which will
# be [thickness=0, rhoM=0, thetaM=first thetaM] unless substrate
# magnetism has been specified
substrate_magnetism = isnan(sigmas[0][0])
if substrate_magnetism:
slices = [[[], [], []]]
else:
slices = [[[0], [0], [blocks[0][2, 0]]]]
interfaces = []
pos = 0
for i, B in enumerate(blocks):
anchor = offsets[i]
w = anchor - pos
if w >= gap_size: # Big gap, so need spacer
# Target average theta between blocks.
if i == 0:
thetaM = B[2, 0]
interfaces.append(0)
else:
thetaM = (B[2, 0] + blocks[i - 1][2, -1]) / 2.
interfaces.append(sigmas[i - 1][1])
slices.append([[w], [0], [thetaM]])
interfaces.append(sigmas[i][0])
elif w >= -1e-6:
# Small gap, so add it to the start of the next block
B[0, 0] += w
anchor -= w
if i == 0:
if not substrate_magnetism:
interfaces.append(sigmas[0][0])
else:
# Use interface_above between blocks which are connected,
# ignoring interface_below.
interfaces.append(sigmas[i - 1][1])
else:
# negative gap should never happen
raise ValueError("Overlapping magnetic layers at %d" % i)
slices.append(B)
nslabs = len(B[0, :])
interfaces.extend([0] * (nslabs - 1))
width = np.sum(B[0, :])
pos = anchor + width
# Add the final slice
w = self.thickness() - pos
theta = blocks[-1][2, -1]
slices.append([[w], [0], [theta]])
interfaces.append(sigmas[-1][1])
wM, rhoM, thetaM = [np.hstack(v) for v in zip(*slices)]
sigmaM = np.array(interfaces)
#print "result", wM, rhoM, thetaM, sigmaM
return wM, sigmaM, rhoM, thetaM
[docs]
def compute_limited_sigma(thickness, roughness, limit):
# Limit roughness to the depths of the surrounding layers. Roughness
# of the first and last layers interfaces is limited only by the
# depth of the first and last layers. We must check explicitly for
# a pure substrate system since that has no limits on roughness.
if limit > 0 and len(thickness) > 2:
s = np.min((thickness[:-1], thickness[1:]), axis=0) / limit
s[+0] = thickness[+1] / limit
s[-1] = thickness[-2] / limit
roughness = np.where(roughness < s, roughness, s)
return roughness
[docs]
def build_profile(z, offset, roughness, value):
"""
Convert a step profile to a smooth profile.
*z* calculation points
*offset* offset for each interface
*roughness* roughness of each interface
*value* target value for each slab
*max_rough* limit the roughness to a fraction of the layer thickness
"""
contrast = np.diff(value)
result = np.zeros_like(z) + value[0]
for offset_k, sigma_k, contrast_k in zip(offset, roughness, contrast):
delta = contrast_k * blend(z, sigma_k, offset_k)
result += delta
return result
SQRT1_2 = 1. / np.sqrt(2.0)
[docs]
def blend(z, sigma, offset):
"""
blend function
Given a Gaussian roughness value, compute the portion of the neighboring
profile you expect to find in the current profile at depth z.
"""
if sigma <= 0.0:
return 1.0 * (z >= offset)
else:
return 0.5 * erf(SQRT1_2 * (z - offset) / sigma) + 0.5