Implement Stepanov-Sinha formalism

udkm1Dsim/__init__.py

@@ -9,10 +9,10 @@
 from .simulations.heat import Heat
 from .simulations.phonons import Phonon, PhononNum, PhononAna
 from .simulations.magnetization import Magnetization, LLB
-from .simulations.xrays import Xray, XrayKin, XrayDyn, XrayDynMag
+from .simulations.xrays import Xray, XrayKin, XrayDyn, XrayDynMag, XrayStepanovSinha
 
 __all__ = ['Atom', 'AtomMixed', 'Layer', 'AmorphousLayer', 'UnitCell', 'Structure',
            'Simulation', 'Heat', 'Phonon', 'PhononNum', 'PhononAna', 'Magnetization', 'LLB',
-           'Xray', 'XrayKin', 'XrayDyn', 'XrayDynMag', 'u', 'Q_']
+           'Xray', 'XrayKin', 'XrayDyn', 'XrayDynMag', 'XrayStepanovSinha', 'u', 'Q_']
 
 __version__ = '1.5.6'

udkm1Dsim/simulations/xrays.py

@@ -22,7 +22,7 @@
 # OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE
 # OR OTHER DEALINGS IN THE SOFTWARE.
 
-__all__ = ['Xray', 'XrayKin', 'XrayDyn', 'XrayDynMag']
+__all__ = ['Xray', 'XrayKin', 'XrayDyn', 'XrayDynMag', 'XrayStepanovSinha']
 
 __docformat__ = 'restructuredtext'
 
@@ -1502,7 +1502,7 @@ def calc_reflectivity_from_matrix(M):
 class XrayDynMag(Xray):
     r"""XrayDynMag
 
-    Dynamical magnetic X-ray scattering simulations.
+    Dynamical magnetic X-ray scattering simulations by Elzo et al. [10]_.
 
     Adapted from Elzo et.al. [10]_ and initially realized in `Project Dyna
     <http://dyna.neel.cnrs.fr>`_.
@@ -2733,3 +2733,250 @@ def calc_roughness_matrix(roughness, k_z, last_k_z):
         W[:, :, 3, 3] = rugosn
 
         return W
+
+
+class XrayStepanovSinha(Xray):
+    r"""XrayStepanovSinha
+
+    Dynamical magnetic X-ray scattering simulations by
+    Stepanov & Sinha [12]_.
+
+    In collaboration with Samuel Flewett
+    (`@sflewett <https://github.com/sflewett>`_)
+
+    Args:
+        S (Structure): sample to do simulations with.
+        force_recalc (boolean): force recalculation of results.
+
+    Keyword Args:
+        save_data (boolean): true to save simulation results.
+        cache_dir (str): path to cached data.
+        disp_messages (boolean): true to display messages from within the
+            simulations.
+        progress_bar (boolean): enable tqdm progress bar.
+
+    Attributes:
+        S (Structure): sample structure to calculate simulations on.
+        force_recalc (boolean): force recalculation of results.
+        save_data (boolean): true to save simulation results.
+        cache_dir (str): path to cached data.
+        disp_messages (boolean): true to display messages from within the
+            simulations.
+        progress_bar (boolean): enable tqdm progress bar.
+        energy (ndarray[float]): photon energies :math:`E` of scattering light
+        wl (ndarray[float]): wavelengths :math:`\lambda` of scattering light
+        k (ndarray[float]): wavenumber :math:`k` of scattering light
+        theta (ndarray[float]): incidence angles :math:`\theta` of scattering
+            light
+        qz (ndarray[float]): scattering vector :math:`q_z` of scattering light
+        polarizations (dict): polarization states and according names.
+        pol_in_state (int): incoming polarization state as defined in
+            polarizations dict.
+        pol_out_state (int): outgoing polarization state as defined in
+            polarizations dict.
+        pol_in (float): incoming polarization factor (can be a complex ndarray).
+        pol_out (float): outgoing polarization factor (can be a complex ndarray).
+        last_atom_ref_trans_matrices (list): remember last result of
+           atom ref_trans_matrices to speed up calculation.
+
+    References:
+
+        .. [12] S. A. Stepanov and S. K. Sinha,
+           *X-Ray Resonant Reflection from Magnetic Multilayers: Recursion
+           Matrix Algorithm*,
+           `Phys. Rev. B 61, 15302 (2000).
+           <https://doi.org/10.1103/PhysRevB.61.15302>`_
+
+    """
+
+    def __init__(self, S, force_recalc, **kwargs):
+        super().__init__(S, force_recalc, **kwargs)
+
+    def __str__(self):
+        """String representation of this class"""
+        class_str = 'Stepanov-Sinha Magnetic X-Ray Diffraction simulation properties:\n\n'
+        class_str += super().__str__()
+        return class_str
+
+    def get_hash(self, **kwargs):
+        """get_hash
+
+        Calculates an unique hash given by the energy :math:`E`, :math:`q_z`
+        range, polarization states as well as the sample structure hash for
+        relevant x-ray and magnetic parameters. Optionally, part of the
+        ``strain_map`` and ``magnetization_map`` are used.
+
+        Args:
+            **kwargs (ndarray[float]): spatio-temporal strain and magnetization
+                profile.
+
+        Returns:
+            hash (str): unique hash.
+
+        """
+        param = [self.pol_in_state, self.pol_out_state, self._qz, self._energy]
+
+        if 'strain_map' in kwargs:
+            strain_map = kwargs.get('strain_map')
+            if np.size(strain_map) > 1e6:
+                strain_map = strain_map.flatten()[0:1000000]
+            param.append(strain_map)
+        if 'magnetization_map' in kwargs:
+            magnetization_map = kwargs.get('magnetization_map')
+            if np.size(magnetization_map) > 1e6:
+                magnetization_map = magnetization_map.flatten()[0:1000000]
+            param.append(magnetization_map)
+
+        return self.S.get_hash(types=['xray', 'magnetic']) + '_' + make_hash_md5(param)
+
+    def set_incoming_polarization(self, pol_in_state):
+        """set_incoming_polarization
+
+        Sets the incoming polarization factor for circular +, circular -, sigma,
+        pi, and unpolarized polarization.
+
+        Args:
+            pol_in_state (int): incoming polarization state id.
+
+        """
+        pass
+
+    def set_outgoing_polarization(self, pol_out_state):
+        """set_outgoing_polarization
+
+        Sets the outgoing polarization factor for circular +, circular -, sigma,
+        pi, and unpolarized polarization.
+
+        Args:
+            pol_out_state (int): outgoing polarization state id.
+
+        """
+        pass
+
+    def homogeneous_reflectivity(self, *args):
+        r"""homogeneous_reflectivity
+
+        Calculates the reflectivity :math:`R` of the whole sample structure
+        allowing only for homogeneous strain and magnetization.
+
+        Args:
+            args (ndarray[float], optional): strains and magnetization for each
+                sub-structure.
+
+        Returns:
+            (tuple):
+            - *R (ndarray[float])* - homogeneous reflectivity.
+            - *R_phi (ndarray[float])* - homogeneous reflectivity for opposite
+              magnetization.
+
+        """
+        pass
+
+    def calculate_chi(self, atom, density, *args):
+        r"""calculate_chi
+
+        Calculates the dielectric susceptibility tensor :math:`\chi_{ij}` of a
+        layer following Eq. (17) of Ref. [12]_:
+
+        .. math::
+
+            \chi_{ij} = (\chi_0 + A) \delta_{ij} - i B \varepsilon_{ijk} M_k
+                + C M_i M_k
+
+        with the :math:`\chi_0` being the mean dielectric susceptibility.
+
+        Add more details from the paper here!
+
+        Args:
+            atom (Atom, AtomMixed): atom or mixed atom.
+            density (float): density around the atom [kg/m³].
+
+        Returns:
+            (tuple):
+            - *chi (ndarray[complex])* - dielectric susceptibility tensor.
+            - *chi_zero (ndarray[complex])* - mean dielectric susceptibility.
+
+        """
+        try:
+            magnetization = args[0]
+            mag_amplitude = magnetization[0]
+            mag_phi = magnetization[1]
+            mag_gamma = magnetization[2]
+        except IndexError:
+            # here we catch magnetizations with only one instead of three
+            # elements
+            try:
+                mag_amplitude = atom.mag_amplitude
+            except AttributeError:
+                mag_amplitude = 0
+            try:
+                mag_phi = atom._mag_phi
+            except AttributeError:
+                mag_phi = 0
+            try:
+                mag_gamma = atom._mag_gamma
+            except AttributeError:
+                mag_gamma = 0
+
+        M = len(self._energy)  # number of energies
+        N = np.shape(self._qz)[1]  # number of q_z
+
+        chi = np.zeros([M, N, 3, 3], dtype=np.cfloat)
+        chi_zero = np.zeros([M, N], dtype=np.cfloat)
+
+        energy = self._energy
+        wl = self._wl
+
+        try:
+            cf = atom.get_atomic_form_factor(energy)
+        except AttributeError:
+            cf = np.zeros_like(energy, dtype=np.cfloat)
+        try:
+            mf = atom.get_magnetic_form_factor(energy)
+        except AttributeError:
+            mf = np.zeros_like(energy, dtype=np.cfloat)
+
+        try:
+            molar_density = density/(atom.mass_number_a/1000)*6.0222e23
+        except AttributeError:
+            molar_density = 0
+        r0 = 2.82e-15  # classical electron radius
+        multiplier = wl**2*r0/np.pi
+
+        chi_zero = molar_density*multiplier*cf
+
+        B = molar_density*multiplier*mf
+        C = 0
+        # This needs to be set in the case that a quadratic term is present
+        # And should be set in the initialization routine as a global parameter
+
+        # magnetization is possibly in cartesian coordinates
+        # need to be converted
+
+        # m1=np.array(M[0,...])
+        # m2=np.array(M[1,...])
+        # m3=np.array(M[2,...])
+        # empty=np.zeros(m1.shape)
+        # ones=empty+1
+        # delta=np.array([[ones,empty,empty],
+        #                 [empty,ones,empty],
+        #                 [empty,empty,ones]])
+        # m_epsilon1=np.array([[empty,empty,empty],
+        #         [empty,empty,m1],
+        #         [empty,-m1,empty]])
+        # m_epsilon2=np.array([[empty,empty,-m2],
+        #         [empty,empty,empty],
+        #         [m2,empty,empty]])
+        # m_epsilon3=np.array([[empty,m3,empty],
+        #         [-m3,empty,empty],
+        #         [empty,empty,empty]])
+        # temp=m_epsilon1+m_epsilon2+m_epsilon3
+        # MM = np.array([[M[j,...]*M[i,...] for i in range(3)] for j in range(3)])
+        # chi = (chi_zero)*delta\
+        #     -complex(0,1)*B*temp+C*MM
+        # self.chi=chi
+        # self.chi_zero=chi_zero*ones
+
+        chi = mag_amplitude*mag_gamma*mag_phi*B*C
+
+        return chi, chi_zero