Functional Layers

Reflectometry layers can be arbitrary functions. This is a rather arbitrary example, with a sinusoidal nuclear profile and an exponential magnetic profile. A simulated dataset is generated from the model.

import numpy as np
from numpy import sin, pi, log, exp, hstack
from bumps.util import push_seed

from refl1d.names import *

FunctionalProfile and FunctionalMagnetism are already available from refl1d.names, but a couple of aliases make them a little easier to access.

from refl1d.flayer import FunctionalProfile as FP
from refl1d.flayer import FunctionalMagnetism as FM

Define the nuclear profile function.

The first parameter to the function is z, which is the points within the layer at which to evaluate the function. The z steps are controlled by the dz parameter to the Experiment definition, which defaults to \(\min(5, \tfrac{1}{10} \tfrac{2 \pi}{Q_\text{max}})\) in angstroms. The remaining parameters become fittable parameters in the model.

The returned value is a complex number whose real part is rho and whose imaginary part is irho. This example is for neutron reflectometry which for the most part does not have a strong absorption cross section.

def nuc(z, period, phase):
    """Nuclear profile"""
    return sin(2*pi*(z/period + phase))

Define the magnetic profile. Like the nuclear profile, the first parameter is z and the remaining parameters become fittable parameters. The returned value is rhoM or the pair rhoM, thetaM, with thetaM defaulting to 0 if it is not returned. Either rhoM or thetaM can be constant.

def mag(z, z1, z2, M1, M2, M3):
    r"""Magnetic profile

    Return the following function:

    .. math::

        f(z) = \left{ \begin{array}{ll}
            C & \mbox{if } z < z_1 \\
            re^{kz} & \mbox{if } z_1 \leq z \leq z_2 \\
            az+b & \mbox{if } z > z_2
            \end{array} \right.

    where :math:`C = M_1`, :math:`r,k` are set such that :math:`re^{kz_1} = M_1` and
    :math:`re^{kz_2} = M_2`, and :math:`a,b` are set such that :math:`az_2 + b = M_2`
    and :math:`az_{\rm end} + b = M_3`.
    """
    # Make sure z1 > z2, swapping if they are different.  Note that in the
    # posterior probability this will set P(z1, z2)=P(z2, z1) always.
    if z1 > z2:
        z1, z2 = z2, z1
    C = M1
    k = (log(M2) - log(M1)) / (z2 - z1)
    r = M1/exp(k*z1)
    a = (M3 - M2) / (z[-1] - z2)
    b = M2 - a*z2

    part1 = z[z < z1]*0+C
    part2 = r*exp(k*z[(z >= z1)&(z <= z2)])
    part3 = a*z[z > z2] + b
    return hstack((part1, part2, part3))

Use these functions to define the functional layer.

flayer = FP(100, 0, name="sin", profile=nuc, period=10, phase=0.2,
            magnetism=FM(profile=mag, M1=1, M2=4, M3=5, z1=10, z2=40))

The functional layer is a normal layer which can be stacked into the model. flayer.start and flayer.end are materials objects whose rho/irho values correspond to the complex rho + j irho value returned by the function at the start and end of the layer. Similarly, magnetism.start and magnetism.end return a magnetic layer defined by the start and end of the magnetic profile.

sample = (silicon(0, 5)
          | flayer
          | flayer.end(35, 15, magnetism=flayer.magnetism.end)
          | air)

Need to be able to compute the thickness of the functional magnetic layer, which unfortunately requires the layer stack and an index. The index can be layer number, layer name, or if there are multiple layers with the same name, (layer name, k), where the magnetism is attached to the kth layer.

flayer.magnetism.set_anchor(sample, 'sin')

Set the fittable parameters. Note that the parameters to the function after the first parameter z become fittable parameters.

sample['sin'].period.range(0, 100)
sample['sin'].phase.range(0, 1)
sample['sin'].thickness.range(0, 1000)
sample['sin'].magnetism.M1.range(0, 10)
sample['sin'].magnetism.M2.range(0, 10)
sample['sin'].magnetism.M3.range(0, 10)
sample['sin'].magnetism.z1.range(0, 100)
sample['sin'].magnetism.z2.range(0, 100)

Define the model. Since this is a simulation, we need to define the incident beam in terms of angles, wavelengths and dispersion. This gets attached to the model forming an experiment. Finally, we simulate data for the experiment with 5% dR/R. We set the seed for the simulation so that the result is reproducible. We could instead set the seed to None so that it pulls a random seed from entropy.

T = np.linspace(0, 5, 100)
xs = [NeutronProbe(T=T, dT=0.01, L=4.75, dL=0.0475, name=name)
      for name in ("--", "-+", "+-", "++")]
probe = PolarizedNeutronProbe(xs)
M = Experiment(probe=probe, sample=sample, dz=0.1)
with push_seed(1):
    M.simulate_data(5)

problem = FitProblem(M)