Merge FreeMagnetismInterface feature enhancement

refl1d/magnetism.py

@@ -38,7 +38,7 @@
 """
 
 from dataclasses import dataclass
-from typing import Dict, List, Literal, Union
+from typing import Optional, Any, Union, Dict, List, Literal
 
 import numpy as np
 from bumps.mono import monospline
@@ -437,3 +437,150 @@ def __str__(self):
 
     def __repr__(self):
         return "FreeMagnetism"
+
+
+@dataclass(init=False)
+class FreeMagnetismInterface(BaseMagnetism):
+    r"""
+    Spline change in magnetism throughout layer.
+
+    Defines monotonic splines for rhoM and thetaM with shared knot positions.
+
+    *dz* is the relative $z$ step between the knots, with position $z_k$
+    defined by $z_k = w \sum_{i=0}^k \delta z_i / \sum_{i=0}^n \delta z_i$,
+    where $n$ is the number of intervals. The resulting $z$ must be monotonic,
+    with $\delta z_i \ge 0$ for all intervals.
+
+    *drhoM* gives the relative $\rho_M$ step between knots. Unlike
+    $\rho_{Mk} = \sum_{i=0}^k \delta \rho_{Mi} / \sum_{i=0}^n \delta \delta \rho_{Mi}$.
+
+    *dthetaM* gives the magnetic angle for each knot.
+
+    *name* is the base name for the various layer parameters.
+
+    *dead_above* and *dead_below* define magnetically dead layers at the
+    nuclear boundaries.  These can be negative if magnetism extends beyond
+    the nuclear boundary.
+
+    *interface_above* and *interface_below* define the magnetic interface
+    at the boundaries, if it is different from the nuclear interface.
+
+    *mbelow* and *mabove* are the rhoM parameter values of
+    the above and below layers respectively. Do not specify if
+    layers either side are not magnetic.
+
+    *tbelow* and *tabove* are the thetaM parameter values of
+    the above and below layers respectively. Do not specify if
+    layers either side are not magnetic.
+    """
+
+    name: Optional[str]
+    mbelow: Optional[Union[Parameter, float]]
+    mabove: Optional[Union[Parameter, float]]
+    tbelow: Optional[Any]
+    tabove: Optional[Any]
+    dz: List[Union[float, Parameter]]
+    drhoM: List[Union[float, Parameter]]
+    dthetaM: List[Union[float, Parameter]]
+
+    magnetic = True
+
+    def __init__(
+        self, dz=(), drhoM=(), dthetaM=(), mbelow=0, mabove=0, tbelow=None, tabove=None, name="MagInterface", **kw
+    ):
+        BaseMagnetism.__init__(self, name=name, **kw)
+
+        def parvec(vector, name, limits):
+            return [Parameter.default(p, name=name + "[%d]" % i, limits=limits) for i, p in enumerate(vector)]
+
+        self.drhoM, self.dthetaM, self.dz = [
+            parvec(v, name + " " + part, limits)
+            for v, part, limits in zip((drhoM, dthetaM, dz), ("drhoM", "dthetaM", "dz"), ((-1, 1), (-1, 1), (0, 1)))
+        ]
+        self.mbelow = Parameter.default(mbelow, name=name + " mbelow", limits=(-np.inf, np.inf))
+        self.mabove = Parameter.default(mabove, name=name + " mabove", limits=(-np.inf, np.inf))
+
+        # if only tabove or tbelow is defined then they are made equal
+        # this is to deal with the situation of a non-magnetic
+        # layer next to a magnetic one
+        if tabove is None and tbelow is None:
+            tbelow = tabove = DEFAULT_THETA_M
+        elif tbelow is None:
+            tbelow = tabove
+        elif tabove is None:
+            tabove = tbelow
+
+        self.tbelow = Parameter.default(tbelow, name=name + " tbelow", limits=(0, 360))
+        self.tabove = Parameter.default(tabove, name=name + " tabove", limits=(0, 360))
+        if len(self.dz) != len(self.drhoM):
+            raise ValueError("Need one dz for each drhoM")
+        if len(self.dthetaM) > 0 and len(self.drhoM) != len(self.dthetaM):
+            raise ValueError("Need one dthetaM for each drhoM")
+
+    def parameters(self):
+        parameters = BaseMagnetism.parameters(self)
+        parameters.update(
+            drhoM=self.drhoM,
+            dthetaM=self.dthetaM,
+            dz=self.dz,
+            mbelow=self.mbelow,
+            mabove=self.mabove,
+            tbelow=self.tbelow,
+            tabove=self.tabove,
+        )
+        return parameters
+
+    def to_dict(self):
+        result = BaseMagnetism.to_dict(self)
+        result["dz"] = to_dict(self.dz)
+        result["drhoM"] = to_dict(self.drhoM)
+        result["dthetaM"] = to_dict(self.dthetaM)
+        result["mbelow"] = to_dict(self.mbelow)
+        result["mabove"] = to_dict(self.mabove)
+        result["tbelow"] = to_dict(self.tbelow)
+        result["tabove"] = to_dict(self.tabove)
+        return result
+
+    def profile(self, Pz, thickness):
+        z = np.hstack((0, np.cumsum(np.asarray([v.value for v in self.dz], "d"))))
+        if z[-1] == 0:
+            z[-1] = 1
+        z *= thickness / z[-1]
+
+        rhoM_fraction = np.hstack((0, np.cumsum(np.asarray([v.value for v in self.drhoM], "d"))))
+        # AJC added since without the line below FreeMagnetismInterface
+        # does not initialise properly - fixes strange behaviour at drho=0 on end point
+        if rhoM_fraction[-1] == 0:
+            rhoM_fraction[-1] = 1
+
+        rhoM_fraction *= 1 / rhoM_fraction[-1]
+        PrhoM = np.clip(monospline(z, rhoM_fraction, Pz), 0, 1)
+
+        if self.dthetaM:
+            thetaM_fraction = np.hstack((0, np.cumsum(np.asarray([v.value for v in self.dthetaM], "d"))))
+            if thetaM_fraction[-1] == 0:
+                thetaM_fraction[-1] = 1
+
+            thetaM_fraction *= 1 / thetaM_fraction[-1]
+            PthetaM = np.clip(monospline(z, thetaM_fraction, Pz), 0, 1)
+        else:
+            # AJC changed from len(z) to PrhoM - since PrhoM is the length of the vector
+            # we want PthetaM to match - otherwise slabs.add_magnetism throws an error
+            PthetaM = np.linspace(0.0, 1.0, len(PrhoM))
+
+        return PrhoM, PthetaM
+
+    def render(self, probe, slabs, thickness, anchor, sigma):
+        Pw, Pz = slabs.microslabs(thickness)
+        rhoM_profile, thetaM_profile = self.profile(Pz, thickness)
+        mbelow, mabove = self.mbelow.value, self.mabove.value
+        tbelow, tabove = self.tbelow.value, self.tabove.value
+        rhoM = (1 - rhoM_profile) * mbelow + rhoM_profile * mabove
+        thetaM = (1 - thetaM_profile) * tbelow + thetaM_profile * tabove
+        slabs.add_magnetism(anchor=anchor, w=Pw, rhoM=rhoM, thetaM=thetaM, sigma=sigma)
+
+    def __str__(self):
+        return "freemagint(%d)" % (len(self.drhoM))
+
+    def __repr__(self):
+        return "FreeMagnetismInterface"

refl1d/mono.py

@@ -13,6 +13,8 @@
 from numpy import asarray, clip, cumsum, diff, hstack, inf, sort
 
 from . import util
+from .magnetism import BaseMagnetism
+from .material import Scatterer
 from .model import Layer
 
 
@@ -115,16 +117,19 @@ class FreeInterface(Layer):
     with slabs.
     """
 
-    name: Optional[str]
-    below: Optional[Any]
-    above: Optional[Any]
+    name: str
+    below: Scatterer
+    above: Scatterer
     thickness: Par
     interface: Par
-    dz: List[Union[float, Par]]
-    dp: List[Union[float, Par]]
+    dz: List[Par]
+    dp: List[Par]
+    magnetism: Optional[BaseMagnetism] = None
     # inflections: List[Any]
 
-    def __init__(self, thickness=0, interface=0, below=None, above=None, dz=None, dp=None, name="Interface"):
+    def __init__(
+        self, thickness=0, interface=0, below=None, above=None, dz=None, dp=None, name="Interface", magnetism=None
+    ):
         self.name = name
         self.below, self.above = below, above
         self.thickness = Par.default(thickness, limits=(0, inf), name=name + " thickness")
@@ -146,6 +151,7 @@ def __init__(self, thickness=0, interface=0, below=None, above=None, dz=None, dp
         self.dp = [Par.default(p, name=name + " dp[%d]" % i, limits=(0, inf)) for i, p in enumerate(dp)]
         self.inflections = Par(name=name + " inflections")
         self.inflections.equals(ParFunction(inflections, dx=self.dz, dy=self.dp))
+        self.magnetism = magnetism
 
     def parameters(self):
         return {
@@ -156,6 +162,7 @@ def parameters(self):
             "below": self.below.parameters(),
             "above": self.above.parameters(),
             "inflections": self.inflections,
+            "magnetism": self.magnetism.parameters() if self.magnetism is not None else None,
         }
 
     def to_dict(self):
@@ -175,14 +182,22 @@ def profile(self, Pz):
         if p[-1] == 0:
             p[-1] = 1
         p *= 1 / p[-1]
+        # AJC included condition as if z[-1] == 0 then z *= thickness/z[-1] == [nan]*len(z)
+        # This then ends with bumps.mono.Monospline adding an extra element as a result of
+        # line 42 in bumps.mono.Monospline
+        if z[-1] == 0:
+            z[-1] = 1
         z *= thickness / z[-1]
         profile = clip(monospline(z, p, Pz), 0, 1)
         return profile
 
     def render(self, probe, slabs):
         thickness = self.thickness.value
+
+        # TODO: why is provided if it is ignored?
         # interface ignored for FreeInterface
         # interface = self.interface.value
+
         below_rho, below_irho = self.below.sld(probe)
         above_rho, above_irho = self.above.sld(probe)
         # Pz is the center, Pw is the width

refl1d/names.py

@@ -26,7 +26,7 @@
 from .flayer import FunctionalMagnetism, FunctionalProfile
 from .instrument import Monochromatic, Pulsed
 from .magnetic import FreeMagnetic, MagneticSlab, MagneticStack, MagneticTwist
-from .magnetism import FreeMagnetism, Magnetism, MagnetismStack, MagnetismTwist
+from .magnetism import FreeMagnetism, FreeMagnetismInterface, Magnetism, MagnetismStack, MagnetismTwist
 from .material import SLD, Compound, Material, Mixture
 
 # Pull in common materials for reflectometry experiments.