PyGTM integration

docs/source/examples.rst

@@ -8,5 +8,6 @@ Examples
     Microscopic 3-Temperature-Model <examples/m3tm.ipynb>
     Landau-Lifshitz-Bloch <examples/LLB.ipynb>
     Phonons <examples/phonons.ipynb>
+    General Transfer Matrix Formalism <examples/gtm.ipynb>
     Dynamical Xray Scattering <examples/dynamical_xray.ipynb>
     Dynamical Magnetic Xray Scattering <examples/dynamical_magnetic_xray.ipynb>

docs/source/examples/Au_Johnson_nk.csv

@@ -0,0 +1,50 @@
+Wavelength. µm	n	k
+0.1879	1.28	1.188
+0.1916	1.32	1.203
+0.1953	1.34	1.226
+0.1993	1.33	1.251
+0.2033	1.33	1.277
+0.2073	1.3	1.304
+0.2119	1.3	1.35
+0.2164	1.3	1.387
+0.2214	1.3	1.427
+0.2262	1.31	1.46
+0.2313	1.3	1.497
+0.2371	1.32	1.536
+0.2426	1.32	1.577
+0.249	1.33	1.631
+0.2551	1.33	1.688
+0.2616	1.35	1.749
+0.2689	1.38	1.803
+0.2761	1.43	1.847
+0.2844	1.47	1.869
+0.2924	1.49	1.878
+0.3009	1.53	1.889
+0.3107	1.53	1.893
+0.3204	1.54	1.898
+0.3315	1.48	1.883
+0.3425	1.48	1.871
+0.3542	1.5	1.866
+0.3679	1.48	1.895
+0.3815	1.46	1.933
+0.3974	1.47	1.952
+0.4133	1.46	1.958
+0.4305	1.45	1.948
+0.4509	1.38	1.914
+0.4714	1.31	1.849
+0.4959	1.04	1.833
+0.5209	0.62	2.081
+0.5486	0.43	2.455
+0.5821	0.29	2.863
+0.6168	0.21	3.272
+0.6595	0.14	3.697
+0.7045	0.13	4.103
+0.756	0.14	4.542
+0.8211	0.16	5.083
+0.892	0.17	5.663
+0.984	0.22	6.35
+1.088	0.27	7.15
+1.216	0.35	8.145
+1.393	0.43	9.519
+1.61	0.56	11.21
+1.937	0.92	13.78

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.scattering import Scattering, GTM, XrayKin, XrayDyn, XrayDynMag
 
 __all__ = ['Atom', 'AtomMixed', 'Layer', 'AmorphousLayer', 'UnitCell', 'Structure',
            'Simulation', 'Heat', 'Phonon', 'PhononNum', 'PhononAna', 'Magnetization', 'LLB',
-           'Xray', 'XrayKin', 'XrayDyn', 'XrayDynMag', 'u', 'Q_']
+           'Scattering', 'GTM', 'XrayKin', 'XrayDyn', 'XrayDynMag', 'u', 'Q_']
 
 __version__ = '2.0.3'

udkm1Dsim/structures/layers.py

@@ -54,8 +54,10 @@
         sound_vel (float): sound velocity.
         phonon_damping (float): phonon damping.
         opt_pen_depth (float): optical penetration depth.
-        opt_ref_index (float): refractive index.
-        opt_ref_index_per_strain (float): change of refractive index per strain.
+        opt_ref_index (complex): refractive index.
+        opt_ref_index_per_strain (complex): change of refractive index per strain.
+        epsilon (complex): permittivity.
+        euler_angles (float): euler angles theta, phi, psi of the permittivity.
         heat_capacity (float): heat capacity.
         therm_cond (float): thermal conductivity.
         lin_therm_exp (float): linear thermal expansion.
@@ -82,6 +84,10 @@
         opt_ref_index_per_strain (ndarray[float]): optical refractive
            index change per strain - real and imagenary part
            :math:`\frac{d n}{d \eta} + i\frac{d \kappa}{d \eta}`.
+        epsilon (list[@lambda]): list of photon frequency [Hz] dependent
+           permittivity elements xx, yy, zz.
+        euler_angles (list[float]): euler angles theta, phi, psi of the permittivity.
+        euler_matrix (ndarray[complex]): euler matrix for permittivity calculations.
         therm_cond (list[@lambda]): list of T-dependent thermal conductivity
             [W/(m K)].
         lin_therm_exp (list[@lambda]): list of T-dependent linear thermal
@@ -121,6 +127,9 @@
         self.opt_pen_depth = kwargs.get('opt_pen_depth', 0.0*u.nm)
         self.opt_ref_index = kwargs.get('opt_ref_index', 0.0+0.0j)
         self.opt_ref_index_per_strain = kwargs.get('opt_ref_index_per_strain', 0.0+0.0j)
+        self.epsilon = kwargs.get('epsilon', 0.0)
+        self.euler_matrix = np.identity(3, dtype=np.complex128)
+        self.euler_angles = kwargs.get('euler_angles', [0.0, 0.0, 0.0]*u.deg)
         self.heat_capacity = kwargs.get('heat_capacity', 0.0)
         self.therm_cond = kwargs.get('therm_cond', 0.0)
         self.lin_therm_exp = kwargs.get('lin_therm_exp', 0.0)
@@ -161,6 +170,8 @@
                       self.opt_ref_index)],
                   ['opt. ref. index/strain', '{0.real:.4f} + {0.imag:.4f}i'.format(
                       self.opt_ref_index_per_strain)],
+                  ['epsilon', self.epsilon],
+                  ['euler angles', self.euler_angles],
                   ['thermal conduct.', ' W/(m K)\n'.join(self.therm_cond_str) + ' W/(m K)'],
                   ['linear thermal expansion', '\n'.join(self.lin_therm_exp_str)],
                   ['heat capacity', ' J/(kg K)\n'.join(self.heat_capacity_str) + ' J/(kg K)'],
@@ -271,7 +282,8 @@
                                'xray': ['num_atoms', '_area', '_mass', '_deb_wal_fac',
                                         '_thickness'],
                                'optical': ['_c_axis', '_opt_pen_depth', 'opt_ref_index',
-                                           'opt_ref_index_per_strain'],
+                                           'opt_ref_index_per_strain', 'epsilon', 'euler_angles',
+                                           '_thickness'],
                                'magnetic': ['_thickness', 'magnetization', 'eff_spin',
                                             '_curie_temp', '_aniso_exponents', '_anisotropy',
                                             '_exch_stiffness', '_mag_saturation', 'lamda'],
@@ -441,6 +453,68 @@
     def opt_pen_depth(self, opt_pen_depth):
         self._opt_pen_depth = opt_pen_depth.to_base_units().magnitude
 
+    @property
+    def epsilon(self):
+        return self._epsilon
+
+    @epsilon.setter
+    def epsilon(self, epsilon):
+        if not type(epsilon) is list:
+            epsilon = [epsilon]
+        self._epsilon = [0, 0, 0]
+        if epsilon[0] != 0:
+            if isfunction(epsilon[0]):
+                self._epsilon[0] = epsilon[0]
+            else:
+                self._epsilon[0] = lambda f: epsilon[0]*np.ones_like(f)
+        else:
+            self._epsilon[0] = lambda f: (1.0 + 0.0j)*np.ones_like(f)
+        try:
+            if isfunction(epsilon[1]):
+                self._epsilon[1] = epsilon[1]
+            else:
+                self._epsilon[1] = lambda f: epsilon[1]*np.ones_like(f)
+        except IndexError:
+            self._epsilon[1] = self._epsilon[0]
+        try:
+            if isfunction(epsilon[2]):
+                self._epsilon[2] = epsilon[2]
+            else:
+                self._epsilon[2] = lambda f: epsilon[2]*np.ones_like(f)
+        except IndexError:
+            self._epsilon[2] = self._epsilon[0]
+
+    @property
+    def euler_angles(self):
+        res = []
+        for eu in self._euler_angles:
+            res.append(Q_(eu, u.rad).to('deg'))
+        return res
+
+    @euler_angles.setter
+    def euler_angles(self, euler_angles):
+        self._euler_angles = []
+        for eu in euler_angles:
+            self._euler_angles.append(eu.to_base_units().magnitude)
+
+        theta = self._euler_angles[0]
+        phi = self._euler_angles[1]
+        psi = self._euler_angles[2]
+
+        self.euler_matrix[0, 0] = np.cos(psi) * np.cos(phi) \
+            - np.cos(theta) * np.sin(phi) * np.sin(psi)
+        self.euler_matrix[0, 1] = -np.sin(psi) * np.cos(phi) \
+            - np.cos(theta) * np.sin(phi) * np.cos(psi)
+        self.euler_matrix[0, 2] = np.sin(theta) * np.sin(phi)
+        self.euler_matrix[1, 0] = np.cos(psi) * np.sin(phi) \
+            + np.cos(theta) * np.cos(phi) * np.sin(psi)
+        self.euler_matrix[1, 1] = -np.sin(psi) * np.sin(phi) \
+            + np.cos(theta) * np.cos(phi) * np.cos(psi)
+        self.euler_matrix[1, 2] = -np.sin(theta) * np.cos(phi)
+        self.euler_matrix[2, 0] = np.sin(theta) * np.sin(psi)
+        self.euler_matrix[2, 1] = np.sin(theta) * np.cos(psi)
+        self.euler_matrix[2, 2] = np.cos(theta)
+
     @property
     def roughness(self):
         return Q_(self._roughness, u.meter).to('nm')
@@ -612,6 +686,14 @@
     def mag_saturation(self, mag_saturation):
         self._mag_saturation = float(mag_saturation.to_base_units().magnitude)
 
+    def get_epsilon_matrix(self, f):
+        epsilon_matrix = np.zeros((len(f), 3, 3), dtype=np.complex128)
+        epsilon_matrix[:, 0, 0] = self._epsilon[0](f)
+        epsilon_matrix[:, 1, 1] = self._epsilon[1](f)
+        epsilon_matrix[:, 2, 2] = self._epsilon[2](f)
+        return np.matmul(np.linalg.inv(self.euler_matrix),
+                         np.matmul(epsilon_matrix, self.euler_matrix))
+
 
 class AmorphousLayer(Layer):
     r"""AmorphousLayer