LLB_te_interpolation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "2c5776aa-73a0-410c-8b70-8abcbd65079f",
   "metadata": {},
   "source": [
    "# LLB implementation with dynamical evaluation of constants"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f430c913-357e-4852-985b-d7c042270dae",
   "metadata": {},
   "source": [
    "In this attempt, I will set up the simulation to only precalulate and interpolate the temperature map and mean field magnetization. All other temperature dependent paramters will be evaluated dynamically from the interpolated temperatures and mean field amgnetizations. To determin the under/obver T_c behavoiur at each timestep and in each grain we precondruct a boolean-array-map saving only under/over tc. To read this information out at desired timesteps we use finderb "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "47bb1c1e-c1ee-4a8e-b707-790b2a871386",
   "metadata": {},
   "source": [
    "## Import packages"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "4157cfbb-130a-42c9-af17-f14a468c1ed5",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy import constants as sp\n",
    "from scipy import optimize as op\n",
    "from scipy import interpolate as ip\n",
    "from matplotlib import pyplot as plt\n",
    "import itertools\n",
    "import time\n",
    "from scipy.integrate import solve_ivp"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "990f4903-2a45-4031-b871-064eeff59a95",
   "metadata": {},
   "source": [
    "## Define finderb"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3f6bbef9-3509-44a9-8909-d163b699de87",
   "metadata": {},
   "source": [
    "Just as used in udkm1dsim. Is used here only to find the proper timestep to select over/under Tc behaviour."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "b1f4ca14-d847-4433-8eac-73e52d7cbef6",
   "metadata": {},
   "outputs": [],
   "source": [
    "def finderb(key, array):\n",
    "    \"\"\"finderb\n",
    "\n",
    "    Binary search algorithm for sorted array. Searches for the first index\n",
    "    ``i`` of array where ``key`` >= ``array[i]``. ``key`` can be a scalar or\n",
    "    a np.ndarray of keys. ``array`` must be a sorted np.ndarray.\n",
    "\n",
    "    Author: André Bojahr.\n",
    "    Licence: BSD.\n",
    "\n",
    "    Args:\n",
    "        key (float, ndarray[float]): single or multiple sorted keys.\n",
    "        array (ndarray[float]): sorted array.\n",
    "\n",
    "    Returns:\n",
    "        i (ndarray[float]): position indices for each key in the array.\n",
    "\n",
    "    \"\"\"\n",
    "    key = np.array(key, ndmin=1)\n",
    "    n = len(key)\n",
    "    i = np.zeros([n], dtype=int)\n",
    "\n",
    "    for m in range(n):\n",
    "        i[m] = finderb_nest(key[m], array)\n",
    "    return i"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "9093b1b9-4d0a-4108-9b0f-f223a8bb8314",
   "metadata": {},
   "outputs": [],
   "source": [
    "def finderb_nest(key, array):\n",
    "    \"\"\"finderb_nest\n",
    "\n",
    "    Nested sub-function of :func:`.finderb` for one single key.\n",
    "\n",
    "    Author: André Bojahr.\n",
    "    Licence: BSD.\n",
    "\n",
    "    Args:\n",
    "        key (float): single key.\n",
    "        array (ndarray[float]): sorted array.\n",
    "\n",
    "    Returns:\n",
    "        a (float): position index of key in the array.\n",
    "\n",
    "    \"\"\"\n",
    "    a = 0  # start of intervall\n",
    "    b = len(array)  # end of intervall\n",
    "\n",
    "    # if the key is smaller than the first element of the\n",
    "    # vector we return 1\n",
    "    if key < array[0]:\n",
    "        return 0\n",
    "\n",
    "    while (b-a) > 1:  # loop until the intervall is larger than 1\n",
    "        c = int(np.floor((a+b)/2))  # center of intervall\n",
    "        if key < array[c]:\n",
    "            # the key is in the left half-intervall\n",
    "            b = c\n",
    "        else:\n",
    "            # the key is in the right half-intervall\n",
    "            a = c\n",
    "\n",
    "    return a"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a9ebc94e-561e-498b-bdb3-538d582c2ccb",
   "metadata": {},
   "source": [
    "## Define material class"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "456cd41e-b032-4187-bdf0-2ea4dd4a3ef1",
   "metadata": {},
   "outputs": [],
   "source": [
    "class material():\n",
    "    def __init__(self, name, S, Tc, lamda, muat, kappa_anis, anis_axis, K_0, A_0, Ms, Delta):\n",
    "        self.name=name                                      # name of the material used for the string representation of the class\n",
    "        self.S=S                                            # effective spin\n",
    "        self.Tc=Tc                                          # Curie temperature\n",
    "        self.J=3*self.S/(self.S+1)*sp.k*self.Tc             # mean field exchange coupling constant\n",
    "        self.mean_mag_map=self.create_mean_mag_map()         # creates the mean magnetization map over temperature as an interpolation function\n",
    "        self.lamda=lamda                                    # intrinsic coupling to bath parameter\n",
    "        self.muat=muat                                      # atomic magnetic moment in units of mu_Bohr\n",
    "        self.kappa_anis=kappa_anis                          # exponent for the temperature dependence of uniaxial anisotropy\n",
    "        self.anis_axis=anis_axis                            # uniaxials anisotropy axis (x:0, y:1, z:2) other anisotropies are not yet implemented \n",
    "        self.K_0=K_0                                        # value for the anisotropy at T=0 K in units of J/m^3\n",
    "        self.A_0=A_0                                        # value for the exchange stiffness at T=0 K in units of J/m\n",
    "        self.Ms=Ms                                          # value for the saturation magnetization at 0K in J/T/m^3\n",
    "        self.Delta=Delta                                    # length of the grain in depth direction in m\n",
    "        \n",
    "    def __str__(self):\n",
    "        return self.name\n",
    "    \n",
    "    def create_mean_mag_map(self):\n",
    "        # This function computes the mean field mean magnetization map by solving the self-consistent equation m=B(m, T)\n",
    "        # As an output we get an interpolation function of the mean field magnetization at any temperature T<=T_c (this can of course be extended to T>T_c with zeros).\n",
    "\n",
    "        # Start by defining a unity function m=m:\n",
    "        def mag(m):\n",
    "            return m\n",
    "\n",
    "        # Define the Brillouin function as a function of scalars, as fsolve takes functions of scalars:\n",
    "        def Brillouin(m, T):\n",
    "            # This function takes input parameters\n",
    "            #   (i) magnetization amplitude m_amp_grid (scalar)\n",
    "            #   (ii) (electron) temperature (scalar)\n",
    "            # As an output we get the Brillouin function evaluated at (i), (ii) (scalar)\n",
    "\n",
    "            eta = self.J * m / sp.k / T /self.Tc\n",
    "            c1 = (2 * self.S + 1) / (2 * self.S)\n",
    "            c2 = 1 / (2 * self.S)\n",
    "            bri_func = c1 / np.tanh(c1 * eta) - c2 / np.tanh(c2 * eta)\n",
    "            return bri_func\n",
    "\n",
    "        # Then we also need a temperature grid. I'll make it course grained for low temperatures (<0.8*Tc) (small slope) and fine grained for large temperatures (large slope):\n",
    "        temp_grid=np.array(list(np.arange(0, 0.8, 1e-3))+list(np.arange(0.8, 1+1e-5, 1e-5)))\n",
    "\n",
    "        # I will define the list of m_eq(T) here and append the solutions of m=B(m, T). It will have the length len(temp_grid) at the end.\n",
    "        meq_list=[1.]\n",
    "\n",
    "        # Define a function to find the intersection of m and B(m, T) for given T with scipy:\n",
    "        def find_intersection_sp(m, Bm, m0):\n",
    "            return op.fsolve(lambda x: m(x) - Bm(x), m0)\n",
    "\n",
    "        # Find meq for every temperature, starting point for the search being (1-T/Tc)^(1/2), fill the list\n",
    "        for i,T in enumerate(temp_grid[1:]):\n",
    "            # Redefine the Brillouin function to set the temperature parameter (I did not find a more elegant solution to this):\n",
    "            def Brillouin_2(m):\n",
    "                return Brillouin(m, T)\n",
    "            # Get meq:\n",
    "            meq=find_intersection_sp(mag, Brillouin_2, np.sqrt(1-T))\n",
    "            if meq[0]<0:            # This is a comletely unwarranted fix for values of meq<0 at temperatures very close to Tc, that fsolve produces. It seems to work though, as the interpolated function plotted by plot_mean_mags() seems clean.\n",
    "                meq[0]*=-1\n",
    "            # Append it to list me(T)\n",
    "            meq_list.append(meq[0])\n",
    "        meq_list[-1]=0              # This fixes slight computational errors to fix m_eq(Tc)=0 (it produces something like m_eq[-1]=1e-7)\n",
    "        return ip.interp1d(temp_grid, meq_list)\n",
    "    \n",
    "    def dbrillouin_t1(self):\n",
    "        return 1/4/self.S**2\n",
    "    \n",
    "    def dbrillouin_t2(self):\n",
    "        return (2*self.S+1)**2/4/self.S**2\n",
    "    \n",
    "    def get_mean_mag(self, T, tc_mask):\n",
    "        # After creating the map, this function can be called to give m_eq at any temperature\n",
    "        # The function takes a 1d-array of temperatures as an input (temperature map at each timestep) and returns an array with the respective mean field equilibrium magnetization\n",
    "        meq=np.zeros(T.shape)\n",
    "        meq[tc_mask]=self.mean_mag_map(T[tc_mask])\n",
    "        return meq\n",
    "    \n",
    "    def alpha_par(self):\n",
    "        # This funtion computes the longitudinal damping parameter alpha_parallel\n",
    "        return 2*self.lamda/(self.S+1)\n",
    "\n",
    "    def qs(self):\n",
    "        # This function computes the first term of the transverse damping parameter alpha_transverse\n",
    "        qs=3*self.Tc/(2*self.S+1)\n",
    "        return qs\n",
    "    \n",
    "    def chi_par_num(self):\n",
    "        return 1/sp.k*self.muat*9.274e-24\n",
    "    \n",
    "    def chi_par_denomm1(self):\n",
    "        return self.J/sp.k\n",
    "    \n",
    "    def anisotropy(self):\n",
    "        #This takes mean field magnetization (1d-array of length N (number of grains)), magnetization vectors (dimension 3xN), magnetization amplitudes (length N) and easy axis ([0,1,2] corresponding to [x,y,z])\n",
    "        return -2*self.K_0"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3fbd4798-2832-4bd5-b4fc-b91db15a7e97",
   "metadata": {},
   "source": [
    "## Create a sample structure"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "023cfecb-0fa2-42cd-9e14-8843c9fc0fd2",
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_sample():\n",
    "    # This is a dummy function that should definitely be replaced by outputs from your code. It does not take any input parameters as I define everything here.\n",
    "    # As an output we get\n",
    "    #   (i) a 1d list of M materials within the sample (materials on the scale of the grainsize of the macrospins)\n",
    "    #   (ii) a 1d numpy array of the actual sample consisting of stacked layers of the M materials\n",
    "    #   (iii-v) magnetization amplitudes and angles\n",
    "    \n",
    "    # Define define three dummy materials with different parameters:\n",
    "    mat_1 = material('Nickel', 0.5, 630., 0.005, 0.393, 3, 2,  0.45e6, -1e-11, 500e3, 1e-9)\n",
    "    mat_2 = material('Cobalt', 1e6, 1480., 0.1, 0.393, 3, 2, 0.45e6, -1e-11, 1400e3, 1e-9)\n",
    "    mat_3 = material('Iron', 2., 1024., 0.005, 2.2, 3, 2, 0.45e6, -1e-11, 200e3, 1e-9)\n",
    "    #FGT = material ('FGT', 0.5, 220., 0.01, 2.2, 3, 2, 0.45e6, 1e-11, 200e-13, 1e-9)\n",
    "    #FGT2 = material ('FGT2', 2., 220., 0.01, 2.2, 3, 2, 0.45e6, 1e-11, 200e-13, 1e-9)\n",
    "    \n",
    "    materials=[mat_2]\n",
    "    \n",
    "    Nickel_1=[mat_1 for _ in range(10)]\n",
    "    Cobalt=[mat_2 for _ in range(15)]\n",
    "    Iron=[mat_3 for _ in range(10)]\n",
    "    Nickel_2=[mat_1 for _ in range(25)]\n",
    "    \n",
    "    sample=np.array([mat_2])\n",
    "    \n",
    "    #The following constructs a list of lists, containing in list[i] a list of indices of material i in the sample_structure. This will help compute the mean field magnetization only once for every material at each timestep. \n",
    "    material_grain_indices=[]\n",
    "    for mat in materials:\n",
    "        material_grain_indices.append([i for i in range(len(sample)) if sample[i]==mat])\n",
    "    material_grain_indices_flat=[index for mat_list in material_grain_indices for index in mat_list]\n",
    "    sample_sorter=np.array([material_grain_indices_flat.index(i) for i in np.arange(len(sample))])\n",
    "    \n",
    "    #The following list locates which material is positioned at which grain of the sample. THis will later be used to define an array of material paramters for the whole sample \n",
    "    mat_locator=[materials.index(grain) for grain in sample]\n",
    "\n",
    "    #Define initial magnetization on the whole sample (for simplicity uniform) and fully magnetized along the z-axis\n",
    "    m_amp = np.ones(len(sample))\n",
    "    m_phi = np.zeros(len(sample))\n",
    "    m_gamma = np.zeros(len(sample))\n",
    "    return materials, sample, m_amp, m_phi, m_gamma, material_grain_indices, sample_sorter, mat_locator"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "048bf70a-d42d-4262-82e4-828aec4b8dfc",
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_mean_mags(materials):\n",
    "    #define a temperature grid:\n",
    "    temps=np.arange(0,2+1e-4, 1e-4)\n",
    "    tc_mask=temps<1.\n",
    "    temps[-1]=1.\n",
    "    for i,m in enumerate(materials):\n",
    "        mmag=get_mean_mag(m, temps, tc_mask)\n",
    "        label=str(m.name)\n",
    "        plt.plot(temps*m.Tc, mmag, label=label)\n",
    "\n",
    "    plt.xlabel(r'Temperature [K]', fontsize=16)\n",
    "    plt.ylabel(r'$m_{\\rm{eq}}$', fontsize=16)\n",
    "    plt.legend(fontsize=14)\n",
    "    plt.title(r'$m_{\\rm{eq}}$ for all materials in sample', fontsize=18)\n",
    "    plt.savefig('plots/meqtest.pdf')\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8be7b3c6-9e88-438e-a638-98016bdb9be7",
   "metadata": {},
   "source": [
    "## Define parameters on sample structure"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "53c7455a-2e4d-4b84-9059-23e93ca45576",
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_ex_stiff_sample(materials, sample, mat_loc, Ms_sam, Delta2_sam):\n",
    "    #This computes a grid for the exchange stiffness in analogous fashion to get_exch_coup_sam()\n",
    "    A_mat=np.zeros((len(materials), len(materials)))\n",
    "    for i, mat in enumerate(materials):\n",
    "        A_mat[i][i]=mat.A_0\n",
    "        \n",
    "    #A_mat[0][1]=-1e-11\n",
    "    #A_mat[1][2]=-5e-11\n",
    "    #A_mat[0][2]=-2.5e-11\n",
    "    \n",
    "    for i in range(1, len(materials)):\n",
    "        for j in range(i):\n",
    "            A_mat[i][j]=A_mat[j][i]\n",
    "            \n",
    "    ex_stiff_arr=np.ones((len(sample),2))*A_mat[0][0]\n",
    "    \n",
    "    for i, grain in enumerate(sample):\n",
    "        if i>0:\n",
    "            ex_stiff_arr[i][0]=A_mat[mat_loc[i]][mat_loc[i-1]]\n",
    "        if i<len(sample)-1:\n",
    "            ex_stiff_arr[i][1]=A_mat[mat_loc[i]][mat_loc[i+1]]\n",
    "    return np.divide(ex_stiff_arr, np.multiply(Ms_sam, Delta2_sam)[:,np.newaxis])\n",
    "\n",
    "def S_sample(sample):\n",
    "    return np.array([mat.S for mat in sample])\n",
    "\n",
    "def Tc_sample(sample):\n",
    "    return np.array([mat.Tc for mat in sample])\n",
    "\n",
    "def J_sample(sample):\n",
    "    return np.array([mat.J for mat in sample])\n",
    "\n",
    "def lamda_sample(sample):\n",
    "    return np.array([mat.lamda for mat in sample])\n",
    "\n",
    "def muat_sample(sample):\n",
    "    return np.array([mat.muat for mat in sample])\n",
    "\n",
    "def Ms_sample(sample):\n",
    "    return np.array([mat.Ms for mat in sample])\n",
    "\n",
    "def Delta2_sample(sample):\n",
    "    return np.array([mat.Delta**2 for mat in sample])\n",
    "\n",
    "def get_ani_sample(sample, Ms_sam):\n",
    "    ani_sam=np.divide(np.array([mat.anisotropy() for mat in sample]), Ms_sam)\n",
    "    kappa_ani_sam=np.array([mat.kappa_anis for mat in sample])\n",
    "    ani_perp_sam= np.ones((len(sample), 3))\n",
    "    for i,mat in enumerate(sample):\n",
    "        ani_perp_sam[i, mat.anis_axis]=0\n",
    "    return ani_sam, kappa_ani_sam, ani_perp_sam\n",
    "\n",
    "def alpha_par_sample(sample):\n",
    "    return np.array([mat.alpha_par() for mat in sample])\n",
    "\n",
    "def qs_sample(sample):\n",
    "    return np.array([mat.qs() for mat in sample])\n",
    "\n",
    "def dbrillouin_t1_sample(sample):\n",
    "    return np.array([mat.dbrillouin_t1() for mat in sample])\n",
    "\n",
    "def dbrillouin_t2_sample(sample):\n",
    "    return np.array([mat.dbrillouin_t2() for mat in sample])\n",
    "\n",
    "def chi_par_num_sample(sample):\n",
    "    return np.array([mat.chi_par_num() for mat in sample])\n",
    "\n",
    "def chi_par_denomm1_sample(sample):\n",
    "    return np.array([mat.chi_par_denomm1() for mat in sample])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f17b5232-5984-4712-aa24-5b2714de1401",
   "metadata": {},
   "source": [
    "## Functions to import temperature map and create corresponding mean_mag_map"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "b5f4ab3c-4a7c-4365-bc6b-e9ea8003efca",
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_tes():\n",
    "    delay=np.load('temp_test/delays.npy')\n",
    "    teNi1=np.load('temp_test/tesNickel0.npy')\n",
    "    teCo2=np.load('temp_test/tesCobalt1.npy')\n",
    "    teFe3=np.load('temp_test/tesIron2.npy')\n",
    "    teNi4=np.load('temp_test/tesNickel3.npy')\n",
    "    tes=np.append(teNi1, teCo2, axis=1)\n",
    "    tes=np.append(tes, teFe3, axis=1)\n",
    "    tes=np.append(tes, teNi4, axis=1)\n",
    "    tes=teCo2\n",
    "    return delay, tes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "d727a0f1-d75f-4866-ba82-6a0ed1eb47e0",
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_tes_paul():    \n",
    "    paul_dat=open('temp_test/paul_data.txt', 'r').readlines()\n",
    "    sep_dat=[line.split() for line in paul_dat]\n",
    "    float_dat=np.array([[float(num) for num in line] for line in sep_dat])\n",
    "\n",
    "    delay=np.array(float_dat[:,0])*1e-12\n",
    "    new_delay=np.arange(0,delay[-1], 1e-16)\n",
    "\n",
    "    tes=float_dat[:,5]\n",
    "    new_tes=np.array(ip.interp1d(delay, tes)(new_delay))\n",
    "    new_tes = np.reshape(new_tes,(-1,1))\n",
    "    \n",
    "    mxs= float_dat[:,1]\n",
    "    mys=float_dat[:,2]\n",
    "    mzs=float_dat[:,3]\n",
    "    m2s=float_dat[:,4]\n",
    "    return new_delay, new_tes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "c97f41c8-149a-4f69-9d7c-37d51230f012",
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_ms_paul():    \n",
    "    paul_dat=open('temp_test/paul_data.txt', 'r').readlines()\n",
    "    sep_dat=[line.split() for line in paul_dat]\n",
    "    float_dat=np.array([[float(num) for num in line] for line in sep_dat])\n",
    "\n",
    "    delay=np.array(float_dat[:,0])*1e-12\n",
    "    new_delay=np.arange(0,delay[-1], 1e-16)\n",
    "\n",
    "    tes=float_dat[:,5]\n",
    "    new_tes=np.array(ip.interp1d(delay, tes)(new_delay))\n",
    "    new_tes = np.reshape(new_tes,(-1,1))\n",
    "    \n",
    "    mxs= float_dat[:,1]\n",
    "    mys=float_dat[:,2]\n",
    "    mzs=float_dat[:,3]\n",
    "    m2s=float_dat[:,4]\n",
    "    return delay, mxs, mys, mzs, m2s"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "9e5a190f-22e8-4d44-b9dc-0fc671d68236",
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_tc_mask(te_red, Tc_sam):\n",
    "    under_tc=te_red<1.\n",
    "    return under_tc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "688f9f64-a7a0-4e78-9f52-8822a80c3273",
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_mean_mag_sample_Ts(te_red, under_tc, mat_gr_ind, mat_gr_ind_flat, materials):\n",
    "    mmag_sam_T=[[materials[i].get_mean_mag(te_red[:,j], under_tc[:,j]) for j in mat_gr_ind[i]] for i in range(len(mat_gr_ind))]\n",
    "    mmag_sam_T_flat=np.concatenate(mmag_sam_T)[mat_gr_ind_flat]\n",
    "    return mmag_sam_T_flat"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3bf00fcc-4c46-4107-a0e0-65573993d6bd",
   "metadata": {},
   "source": [
    "## Transform initial polar magnetization input into cartesian coordinates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "7b722da9-630d-4ae2-a460-a71293de144e",
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_mag(polar_dat):\n",
    "    # This function takes as input parameters the amplitude and angles (A, gamma, phi) and puts out a numpy array of dimension 3xlen(sample)\n",
    "    # with 3 magnetization components for len(sample) grains\n",
    "    amp=polar_dat[0,:]\n",
    "    gamma=polar_dat[1,:]\n",
    "    phi=polar_dat[2,:]\n",
    "    sin_phi=np.sin(phi)\n",
    "    \n",
    "    mx=amp*sin_phi*np.cos(gamma)\n",
    "    my=amp*sin_phi*np.sin(gamma)\n",
    "    mz=amp*np.cos(phi)\n",
    "    \n",
    "    return np.array([mx,my,mz]).T"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8f378ac8-c172-4193-b09d-8e688059f1d9",
   "metadata": {},
   "source": [
    "## Setup of temperature dynamic parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "cf288009-7a3b-46bf-8cbe-bbaea4b13a8b",
   "metadata": {},
   "outputs": [],
   "source": [
    "def ani_sample_T(mmag_sam_T, K0_sam, kappa_ani_sam):\n",
    "    return np.multiply(K0_sam,np.power(mmag_sam_T,kappa_ani_sam-2))\n",
    "\n",
    "def ex_stiff_sample_T(mmag_sam_T, ex_stiff_sam):\n",
    "    return np.multiply(np.power(mmag_sam_T[:, np.newaxis],2-2),ex_stiff_sam)\n",
    "\n",
    "def qs_sample_T(qs_sam, mmag_sam_T, T):\n",
    "    return qs_sam*mmag_sam_T/T\n",
    "\n",
    "def alpha_par_sample_T(mmag_sam_T, T, alpha_par_sam, qs_sam_T, Tc_sam, under_tc, over_tc, lambda_sam):\n",
    "    apsT=np.zeros(len(T))\n",
    "    apsT[under_tc]=alpha_par_sam[under_tc]/np.sinh(2*qs_sam_T[under_tc])\n",
    "    apsT[over_tc]=lambda_sam[over_tc]*2/3*np.divide(T[over_tc], Tc_sam[over_tc])\n",
    "    return apsT\n",
    "\n",
    "def alpha_trans_sample_T(mmag_sam_T, lamda_sam, T, qs_sam_T, Tc_sam, under_tc, over_tc, lambda_sam):\n",
    "    atsT=np.zeros(len(T))\n",
    "    atsT[under_tc]=np.multiply(lambda_sam[under_tc], (np.divide(np.tanh(qs_sam_T[under_tc]), qs_sam_T[under_tc])-np.divide(T[under_tc],3*Tc_sam[under_tc])))\n",
    "    atsT[over_tc]=lambda_sam[over_tc]*2/3*np.divide(T[over_tc], Tc_sam[over_tc])\n",
    "    return atsT\n",
    "\n",
    "def eta_sample_T(mmag_sam_T, J_sam, T):\n",
    "    return J_sam*mmag_sam_T/sp.k/T\n",
    "\n",
    "def dbrillouin_sample_T(eta_sam_T, S_sam, dbrillouin_t1_sam, dbrillouin_t2_sam):\n",
    "    two_S_sam=2*S_sam\n",
    "    x1=np.divide(eta_sam_T,two_S_sam)\n",
    "    x2=np.divide(np.multiply(eta_sam_T,(two_S_sam+1)),two_S_sam)\n",
    "    sinh_func=1/np.sinh(np.array([x1,x2]))**2\n",
    "    dbrillouin_sam_T=dbrillouin_t1_sam*sinh_func[0]-dbrillouin_t2_sam*sinh_func[1]\n",
    "    return dbrillouin_sam_T\n",
    "\n",
    "def chi_par_sample_T(chi_par_num_sam, chi_par_denomm1_sam, dbrillouin_sam_T, T, under_tc, over_tc, muat_sam, Tc_sam, J_sam):\n",
    "    cpsT=np.zeros(len(T))\n",
    "    cpsT[under_tc]=np.multiply(chi_par_num_sam[under_tc], np.divide(dbrillouin_sam_T, T[under_tc]-np.multiply(chi_par_denomm1_sam[under_tc], dbrillouin_sam_T)))\n",
    "    cpsT[over_tc]=np.divide(np.multiply(muat_sam[over_tc]*9.274e-24, Tc_sam[over_tc]), J_sam[over_tc]*(T[over_tc]-Tc_sam[over_tc]+1e-1))\n",
    "    return cpsT"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8d528a61-6fef-4e48-8137-f64383df95e2",
   "metadata": {},
   "source": [
    "## Components of effective field (magnetization dependent)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "8efba811-800b-4b1e-b4ee-1fbe9c6dfe11",
   "metadata": {},
   "outputs": [],
   "source": [
    "def anis_field(anis_sam_T, m, ani_perp_sam):\n",
    "    return anis_sam_T[:, np.newaxis]*(m*ani_perp_sam)\n",
    "\n",
    "def ex_field(ex_stiff_sam_T, m_diff_up, m_diff_down):\n",
    "    return ex_stiff_sam_T[:,0][:,np.newaxis]*m_diff_up+ex_stiff_sam_T[:,1][:,np.newaxis]*m_diff_down\n",
    "\n",
    "def th_field(m, m_squared, mmag_sam_T, T, Tc_sam, chi_par_sam_T, under_tc, over_tc):\n",
    "    factor = 1/chi_par_sam_T\n",
    "    H_th = np.zeros(len(T))\n",
    "    H_th[under_tc] = (1-m_squared[under_tc]/mmag_sam_T[under_tc]**2)*factor[under_tc]/2\n",
    "    H_th[over_tc] = -(1+3/5*Tc_sam[over_tc]/(T[over_tc]-Tc_sam[over_tc]+1e-1))*m_squared[over_tc]*factor[over_tc]\n",
    "    return H_th[:, np.newaxis]*m  "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "69fe7c0f-6ed1-4742-8b61-071ed5108213",
   "metadata": {},
   "source": [
    "## Main function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "42300e13-f396-459a-a102-260139a1a976",
   "metadata": {},
   "outputs": [],
   "source": [
    "def integrate_LLB():\n",
    "    starttime=time.time()\n",
    "    gamma=1.76e11 # gyromagnetic ratio in (Ts)^{-1}\n",
    "\n",
    "    # import temperature map\n",
    "    times, tes=get_tes_paul()\n",
    "    N=len(times)\n",
    "    \n",
    "    # load a sample and call the functions to get all parameters on the sample structure:\n",
    "    materials, sample, m_amp, m_phi, m_gamma, mat_gr_ind, mat_gr_ind_flat, mat_loc=get_sample()\n",
    "    \n",
    "    H_ext=np.array([[0,0,0] for _ in sample]) #external field in T\n",
    "    \n",
    "    \n",
    "    # load initial magnetization\n",
    "    m0=get_mag(np.array([m_amp, m_phi, m_gamma])).flatten()\n",
    "    m0=np.array([0.8, 0., 0.6])\n",
    "    \n",
    "    # Create static parameters on sample structure\n",
    "    Ms_sam=Ms_sample(sample)\n",
    "    Delta2_sam=Delta2_sample(sample)\n",
    "    S_sam=S_sample(sample)\n",
    "    Tc_sam=Tc_sample(sample)\n",
    "    J_sam=J_sample(sample)\n",
    "    lamda_sam=lamda_sample(sample)\n",
    "    muat_sam=muat_sample(sample)\n",
    "    ex_stiff_sam=get_ex_stiff_sample(materials, sample, mat_loc, Ms_sam, Delta2_sam)\n",
    "    K0_sam, kappa_ani_sam, ani_perp_sam=get_ani_sample(sample, Ms_sam)\n",
    "    alpha_par_sam=alpha_par_sample(sample)\n",
    "    qs_sam=qs_sample(sample)\n",
    "    dbrillouin_t1_sam=dbrillouin_t1_sample(sample)\n",
    "    dbrillouin_t2_sam=dbrillouin_t2_sample(sample)\n",
    "    chi_par_num_sam=chi_par_num_sample(sample)\n",
    "    chi_par_denomm1_sam=chi_par_denomm1_sample(sample)\n",
    "    \n",
    "    # create boolean array seperating temperature values under and over the respective Curie temperatures of materials in the sample:\n",
    "    Te_red=np.divide(tes, Tc_sam)\n",
    "    utc=get_tc_mask(Te_red, Tc_sam)\n",
    "    otc=~utc\n",
    "\n",
    "    # create the mean magnetization map for the imported temperature map:\n",
    "    mms=get_mean_mag_sample_Ts(Te_red, utc, mat_gr_ind, mat_gr_ind_flat, materials)\n",
    "    \n",
    "    # timecheck incoming\n",
    "    inichecktime=time.time()\n",
    "    initime=inichecktime-starttime\n",
    "    print('Time used to import temperature map and initialize mean mag map and set up static parameters on sample strucure:' , str(initime) , 's')\n",
    "    \n",
    "    # define and call interpolation functions for temperature map and mean_magnetization\n",
    "    def mean_mag_interpol(t):\n",
    "        return ip.interp1d(t,mms)\n",
    "    \n",
    "    def te_map(t):\n",
    "        return ip.interp1d(t, tes.T)\n",
    "    \n",
    "    mmag_ipl=mean_mag_interpol(times)\n",
    "    tes_ipl=te_map(times)\n",
    "    \n",
    "    #another timecheck\n",
    "    inichecktime2=time.time()\n",
    "    ini2time=inichecktime2-inichecktime\n",
    "    print('Creation of interpolation functions of (implicitly) time dependent parameters took' , str(ini2time), 's')\n",
    "    \n",
    "    def LLB(t, m_flat, mmag_ipl, tes_ipl):\n",
    "        \n",
    "        mmag_sam_T=mmag_ipl(t)\n",
    "        Te=tes_ipl(t)\n",
    "        \n",
    "        # find the closest timestep and select the corresponding boolean array for Tc_selection\n",
    "        t_index=finderb(t, times)[0]\n",
    "        under_tc=utc[t_index, :]\n",
    "        over_tc=~under_tc\n",
    "        \n",
    "        # unflatten the magnetization input and compute squares and\n",
    "        m=m_flat.reshape(len(sample), 3)\n",
    "        m_squared=np.sum(np.power(m,2), axis=-1)\n",
    "        m_diff_down=np.concatenate((np.diff(m, axis=0), np.zeros((1, 3))), axis=0)\n",
    "        m_diff_up=-np.roll(m_diff_down, 1)\n",
    "        \n",
    "        # evaluate the temperature dependence of parameters\n",
    "        anis_sam_T=ani_sample_T(mmag_sam_T, K0_sam, kappa_ani_sam)\n",
    "        ex_stiff_sam_T=ex_stiff_sample_T(mmag_sam_T, ex_stiff_sam)\n",
    "        qs_sam_T=qs_sample_T(qs_sam, mmag_sam_T, Te)\n",
    "        eta_sam_T=eta_sample_T(mmag_sam_T, J_sam, Te)\n",
    "        dbrillouin_sam_T=dbrillouin_sample_T(eta_sam_T[under_tc], S_sam[under_tc], dbrillouin_t1_sam[under_tc], dbrillouin_t2_sam[under_tc])\n",
    "        chi_par_sam_T=chi_par_sample_T(chi_par_num_sam, chi_par_denomm1_sam, dbrillouin_sam_T, Te, under_tc, over_tc, muat_sam, Tc_sam, J_sam)\n",
    "\n",
    "\n",
    "        # and the effective field\n",
    "        H_ex=ex_field(ex_stiff_sam_T, m_diff_up, m_diff_down)\n",
    "        H_ani=anis_field(anis_sam_T, m, ani_perp_sam)\n",
    "        H_th=th_field(m, m_squared, mmag_sam_T, Te, Tc_sam, chi_par_sam_T, under_tc, over_tc)\n",
    "\n",
    "        H_eff=H_ani+H_ex+H_th+H_ext\n",
    "\n",
    "        # only missing the damping parameters:\n",
    "        alpha_par_sam_T=alpha_par_sample_T(mmag_sam_T, Te, alpha_par_sam, qs_sam_T, Tc_sam, under_tc, over_tc, lamda_sam)\n",
    "        alpha_trans_sam_T=alpha_trans_sample_T(mmag_sam_T, lamda_sam, Te, qs_sam_T, Tc_sam, under_tc, over_tc, lamda_sam)\n",
    "\n",
    "        # we will precompute the prefactors that are not of dimension len(sample)x3:\n",
    "\n",
    "        pref_trans=np.divide(alpha_trans_sam_T, m_squared)\n",
    "        pref_long=np.multiply(np.divide(alpha_par_sam_T, m_squared),np.einsum('ij,ij->i', m, H_eff))\n",
    "\n",
    "        # and all the cross products:\n",
    "        m_rot=np.cross(m,H_eff) #precessional term\n",
    "        m_trans=np.cross(m,m_rot) #transverse damping term\n",
    "\n",
    "        trans_damp= np.multiply(pref_trans[:,np.newaxis],m_trans)\n",
    "        long_damp= np.multiply(pref_long[:,np.newaxis], m)\n",
    "\n",
    "        # Now compute the magnetization increment\n",
    "        dm_dt=gamma*(-m_rot -trans_damp + long_damp)\n",
    "        return dm_dt.flatten()\n",
    "    \n",
    "    #Now call the solver with arguments of function, initial condition, time-array and parameterlist (the gathering of interpolation functions)\n",
    "    solved_LLB=solve_ivp(fun=lambda t,m: LLB(t, m, mmag_ipl, tes_ipl), t_span=(0,times[-1]), y0=m0, method='RK45')\n",
    "    \n",
    "    \n",
    "    #last timecheck\n",
    "    endtime=time.time()\n",
    "    dyntime=endtime-inichecktime\n",
    "    alltime=endtime-starttime\n",
    "    print('Dynamical calculation took' , str(dyntime) , 's')\n",
    "    print('The whole simulation took' , str(alltime), 's')\n",
    "    \n",
    "    #this is the mag map output:\n",
    "    return solved_LLB"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "be358376-f699-4d22-a670-2681a9691076",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\tgrie\\AppData\\Local\\Packages\\PythonSoftwareFoundation.Python.3.8_qbz5n2kfra8p0\\LocalCache\\local-packages\\Python38\\site-packages\\scipy\\optimize\\minpack.py:175: RuntimeWarning: The iteration is not making good progress, as measured by the \n",
      "  improvement from the last ten iterations.\n",
      "  warnings.warn(msg, RuntimeWarning)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Time used to import temperature map and initialize mean mag map and set up static parameters on sample strucure: 11.268680095672607 s\n",
      "Creation of interpolation functions of (implicitly) time dependent parameters took 0.8874285221099854 s\n",
      "Dynamical calculation took 410.8719754219055 s\n",
      "The whole simulation took 422.1406555175781 s\n"
     ]
    }
   ],
   "source": [
    "mag_map=integrate_LLB()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "e0b8f55c-458a-4931-b1ab-b372464651e6",
   "metadata": {},
   "outputs": [],
   "source": [
    "sim_delay=mag_map.t\n",
    "mm_rs=mag_map.y.T.flatten().reshape(len(mag_map.t), 1,3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "90646f35-9593-4ad6-8ed9-72972c7fa2cd",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAAEMCAYAAAAh7MZPAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/av/WaAAAACXBIWXMAAAsTAAALEwEAmpwYAAAclElEQVR4nO3de9wdVX3v8c83CRDut0SIJJCgqRqpcnmEUlrlaNCACPbYloucQkHzUkpFUSs0fVFE+zpcDh5boUKKlCgot6qkEEVNoVoUTLiTQEgMCAEUwrFcRAiB3/lj1iaTzb7MfvZ9nu/79ZrXnlkze+b3rOz8nvWsmb2WIgIzMyuncf0OwMzMusdJ3sysxJzkzcxKzEnezKzEnOTNzErMSd7MrMQm9DuATpu0tWL6ZLJfX+MA1VmvbIu0Mb5qZ/WB43MnUJP16oU629Qpo8Yx1WW1jq9X1op23z/E7rgNpgFbApsB48YDm6RlfG6p9cHKfwag/meAqjJq7IPm/8719tU7ttXyRsrzGbntttvXRsTkds4xZ877Yu3atUWvd0NEzGnneq0qXZKfviMs/SKwBTAxLfn16rIJ44CtcsvE9LpZbnti1fqEtD+/XXmtXp9Alhgm1Fiosd7otdF6o7JG5aM9roS2FMwD9gN2B7baFtgZ2AnYLi2Vz0n+37/WZ6B6qXwGoLP/5s22R1veSK/e033SJr9s9xxr165l6dKbC15v80ntXq9Vg1nzZmZDI4D1/Q6iLid5M7O2OMmbmZWYk7yZWYk5yZuZlVgAL/Q7iLqc5M3M2hLAy/0Ooi4neTOztri7xsysxJzkzcxKzknezKykXsE3Xs3MSsvdNWZmJeYkb2ZWck7yZmYl5Za8mVmJOcmbmZXYYD9d05fp/yTtIOmHklam1+0bHLuNpDWSzu9ljGZmxa0vuPRev+Z4PRVYHBEzgcVpu54vAD/uSVRmZi2rdNc4yecdDixI6wuAD9Y6SNI+ZHOv/aA3YZmZtcpJvpadIuLxtP4rskS+EUnjgPOAzzQ7maS5kpZKWvrkc50N1MysscoolEWW3uvajVdJPyKbBbnavPxGRISkqHHcicCiiFgjNZ4dPiLmA/MBRqbXPJeZWZeM0fHkI2J2vX2Sfi1pSkQ8LmkK8ESNw/YH/ljSicBWwKaSnouIRv33ZmY95kcoa1kIHAuclV6vrT4gIj5cWZd0HDDiBG9mg2ewk3y/+uTPAg6StBKYnbaRNCLp4j7FZGY2CoN947UvLfmIeAp4T43ypcBHapRfClza9cDMzFrmlryZWcl1riUvaY6kFZJWSXpNF7WkXSXdKOkOSXdLOqTR+TysgZlZWzo3rIGk8cAFwEHAGmCJpIURsTx32N8BV0XEVyXNAhYB0+ud0y15M7O2dLRPfl9gVUSsjoh1wBVkXx6tvuA2aX1b4LFGJ3RL3sysLS31yU+StDS3PT99z6diF+CR3PYaYL+qc5wB/EDSXwNbkj28UpeTvJlZ2won+bURMdLmxY4CLo2I8yTtD3xD0h4R8Uqtg53kzcza0tGnax4FpuW2p6ayvBOAOQAR8TNJE4FJ1P5Sqfvkzcza09E++SXATEkzJG0KHEn25dG8h0mPoEt6CzAReLLeCd2SNzNrSwAvduZMEeslnQTcAIwHLomIZZLOBJZGxELg08C/SPpUuvhxEVF3zC4neTOztnT2y1ARsYjssch82em59eXAAUXP5yRvZtaWwf7Gq5O8mVlbnOTNzErMSd7MrMTG6KQhZmZjg1vyZmYl5iRvZlZiTvJmZiXmJG9mVnJO8mZmJdW5SUO6wUnezKwt7q4xMysxJ3kzs5J7ud8B1OUkb2bWFrfkzcxKzEnezKzEPHaNmVmJuSVvZlZiTvJmZiXmJG9mVm7hRyjNzMopgJf6HUR9TvJmZu0IBvm7UE7yZmZte6XfAdTnJG9m1g635M3MSm6AW/Lj+nFRSTtI+qGklel1+xrH7CnpZ5KWSbpb0hH9iNXMrKFKS77I0gd9SfLAqcDiiJgJLE7b1Z4H/iIi3grMAb4sabvehWhmVkDl6ZoiSx/0K8kfDixI6wuAD1YfEBEPRMTKtP4Y8AQwuVcBmpkVNsAt+X71ye8UEY+n9V8BOzU6WNK+wKbAL7odmJlZS8bqjVdJPwJ2rrFrXn4jIkJSNDjPFOAbwLERUfP2hqS5wFyAXXccdchmZqMzwDdeu5bkI2J2vX2Sfi1pSkQ8npL4E3WO2wa4HpgXEbc0uNZ8YD7AyPT6vzDMzDpuwFvy/eqTXwgcm9aPBa6tPkDSpsB3gK9HxDU9jM3MrDWvFFz6oF9J/izgIEkrgdlpG0kjki5Ox/w58E7gOEl3pmXPvkRrZlbPK8C6gksf9OXGa0Q8BbynRvlS4CNp/TLgsh6HZmbWurHYJ29mNiYMeJ+8k7yZWbuc5M3MSipwd42ZWWl50hAzs5Ib4O6afj1CaWZWDh0ehVLSHEkrJK2SVGvwRiT9uaTlaZTebzY6X9OWvKTJwEeB6fnjI+L4YiGbmZVch/rkJY0HLgAOAtYASyQtjIjluWNmAqcBB0TEbyS9rtE5i3TXXAv8BPgRA/1HiZlZH3T2Ecp9gVURsRpA0hVko/Yuzx3zUeCCiPgNQETUHBamokiS3yIiPje6eM3MxoDiSX6SpKW57flp7K2KXYBHcttrgP2qzvF7AJJuBsYDZ0TE9+tdsEiSv07SIRGxqMCxZmZjS2tP16yNiJE2rzgBmAkcCEwFfizp9yPiv2sdXOTG68lkif4FSc+m5Zk2gzQzK4fKc/KdGaDsUWBabntqKstbAyyMiJci4kHgAbKkX1PTJB8RW0fEuIiYmNa3johtCoVrZjYWdO7pmiXATEkz0ki8R5KN2pv3XbJWPJImkXXfrK53wkLPyUs6jGxESICbIuK6QuGamZVdB2+8RsR6SScBN5D1t18SEcsknQksjYiFad97JS1PV/5sGvSxpiKPUJ4FvAO4PBWdLOmAiDitzZ/HzKwcOjisQbr/uaiq7PTcegCnpKWpIi35Q4A9K1PvSVoA3EH2nKaZ2dhWGU9+QBX9xut2ufVtuxCHmdnwGuCZoYq05P83cIekGwGR9c3X/KqtmdmYM+zjyUfEtyTdRNYvD/C5iPhVV6MyMxsmwzjUsKQ3R8T9kvZORWvS6+slvT4ibu9+eGZmA26IW/KnAHOB82rsC+DdXYnIzGzYDGOSj4i5afXgiHghv0/SxK5GZWY2LAZ80pAiT9f8tGCZmdnY09lhDTquUZ/8zmQjom0uaS+yJ2sAtgG26EFsZmbDYRi7a4D3AceRDZDzpVz5s8DfdjEmM7PhMaw3XiNiAbBA0oci4t96GJOZ2XAZ0kcoj4mIy4Dpkl4zRkJEfKnG28zMxpZhbckDW6bXrXoRiJnZUBrwp2sadddclF4/37twzMyGzIC35Js+QinpHEnbSNpE0mJJT0o6phfBmZkNhQF+hLLIc/LvjYhngEOBh4A3Ap/tZlBmZkOj0pLvzMxQHVdkFMrKMe8Hro6IpyU1Ot7MbGwZ4O6aIkn+Okn3A78DPi5pMvBCk/eYmY0NA37jtchE3qcCfwiMRMRLwG+Bw7sdmJnZUBj27hpJmwDHAO9M3TT/CVzY5bjMzIbHMH4ZKuerwCbAP6ft/5XKPtKtoMzMhsaAP0JZJMm/IyLentv+D0l3deLikuYA/wiMBy6OiLOq9m8GfB3YB3gKOCIiHurEtc3MOmaAW/JFHqF8WdIbKhuSdqcDv7ckjQcuAA4GZgFHSZpVddgJwG8i4o3A/wXObve6ZmYdNex98mTPxN8oaTXZcMO7AX/ZgWvvC6yKiNUAkq4gu6G7PHfM4cAZaf0a4HxJiojowPXNzNo34E/XFJnIe7GkmcCbUtGKiHixA9feBXgkt70G2K/eMRGxXtLTwI7A2vxBkuaSTVXIrjt2IDIzs1YMc598murvROCPyH5n/UTShdVTAvZTRMwH5gOMTJdb+WbWOwN+47VIn/zXgbcCXwHOT+vf6MC1HwWm5banprKax0iaAGxLdgPWzGxwDPDYNUX65PeIiPwN0RslLa97dHFLgJmSZpAl8yOBo6uOWQgcC/wM+FPgP9wfb2YDpQQt+dsl/UFlQ9J+wNJ2LxwR64GTgBuA+4CrImKZpDMlHZYO+xqwo6RVwCnAqe1e18yso4Z1Iu+cfYCfSno4be8KrJB0DxAR8bbRXjwiFgGLqspOz62/APzZaM9vZtZ1AazrdxD1FUnyc7oehZnZMBvgL0MVeYTyl70IxMxsKA14n3yRlryZmTXiJG9mVlKVG68DyknezKwdAz6sQZGJvP+npJWSnpb0jKRnJT3Ti+DMzIbCkA9Qdg7wgYi4r9vBmJkNnRLceP21E7yZWQND3ie/VNKVwHeBV0efjIhvdysoM7OhMeAt+SLDGmwDPA+8F/hAWg7tZlBmZkOjw5OGSJojaYWkVZLqDuUi6UOSQtJIo/MV+TJUJyYIMTMrpw4+XZObMe8gsjk2lkhaGBHLq47bGjgZuLXZOesmeUl/ExHnSPoK2Y+xkYj4RIvxm5mVU+f65IvMmAfwBbLpUD/b7ISNWvKVm61tjzhpZlZarfXJT5KUz6nz06RHFU1nzJO0NzAtIq6XNPokHxH/nl4XFInczGzMKp7k10ZEwz70RiSNA74EHFf0PUWm/5sMfA6YBUyslEfEu1sP0cysZDo7rEGzGfO2BvYAbpIEsDOwUNJhEVGz16XI0zWXk3XdzAA+DzxENquTmZlBJ5+ueXXGPEmbks2Yt7CyMyKejohJETE9IqYDtwB1EzwUS/I7RsTXgJci4j8j4njArXgzM9jwdE2Rpdmpis2Y15IiX4aqhPa4pPcDjwE7jOZiZmal0+EvQzWbMa+q/MBm5yuS5L8oaVvg08BXyL4c9akC7zMzGxuGdViD9GD+zIi4Dnga+B89icrMbFgM87AGEfEycFSPYjEzG05DPtTwzZLOB64EflspjIjbuxaVmdmwGPBJQ4ok+T3T65m5ssBP2JiZDXx3TZEBytwPb2bWyLDeeAWQdEqN4qeB2yLizo5HZGY2ZAa4IV+ou2YkLf+etg8F7gY+JunqiDinW8GZmQ26Ae+tKZTkpwJ7R8RzAJL+HrgeeCdwG9kcsGZmY9YA99YUSvKvIzftH9l95J0i4neSXqzzHjOzMWHAH64plOQvB26VdG3a/gDwTUlb8tqB7M3MxpSh766JiC9I+h5wQCr6WG7Esw93LTIzsyEw9EkeICV1zxBlZlbDsPfJm5lZHaVoyZuZWX1O8mZmJTXoT9cUmRmqayTNkbRC0ipJp9bYf4qk5ZLulrRY0m79iNPMrJ7KFK9Fln7oW5JPY9VfABxMNkn4UZJmVR12BzASEW8DrsFfvDKzATTAIw33tSW/L7AqIlZHxDrgCuDw/AERcWNEPJ82byH79q2Z2cCo3Hh1kn+tXYBHcttrUlk9JwDf62pEZmajMMjdNUNx41XSMWSDpL2rzv65wFyAXXfsYWBmNua9AqzrdxAN9LMl/ygwLbc9NZVtRNJsYB5wWETUHCsnIuZHxEhEjEzeqiuxmpnV5ZZ8bUuAmZJmkCX3I4Gj8wdI2gu4CJgTEU/0PkQzs8b8Zag6ImK9pJOAG4DxwCURsUzSmcDSiFgInAtsBVwtCeDhiDisXzGbmdXiYQ3qiIhFwKKqstNz67N7HpSZWQvckjczKzkneTOzkhr0YQ2c5M3M2uDuGjOzkvONVzOzknJL3sys5NySNzMrqWCwhzVwkjcza0NlPPlB5SRvZtYm98mbmZWUb7yamZWcu2vMzErKLXkzsxLzsAZmZiU26C35fs4MZWZWCp2cGUrSHEkrJK2SdGqN/adIWi7pbkmLJe3W6HxO8mZmbai05IsszUgaD1wAHAzMAo6SNKvqsDuAkYh4G3ANcE6jczrJm5m1qVNJHtgXWBURqyNiHXAFcHj+gIi4MSKeT5u3kM2PXZeTvJlZGyrfeC3YXTNJ0tLcMrfqdLsAj+S216Syek4AvtcoPt94NTNrQ4tP16yNiJFOXFfSMcAI8K5GxznJm5m1qYNP1zwKTMttT01lG5E0G5gHvCsiXmx0Qid5M7M2dPgRyiXATEkzyJL7kcDR+QMk7QVcBMyJiCeandBJ3sysTZ0a1iAi1ks6CbgBGA9cEhHLJJ0JLI2IhcC5wFbA1ZIAHo6Iw+qd00nezKwNnf4yVEQsAhZVlZ2eW5/dyvmc5M3M2uBhDczMSmzQhzVwkjcza5OHGjYzKym35M3MSs5J3syspDyRt5lZiQWwrt9BNOAkb2bWJrfkzcxKyjdezcxKzi15M7OSckvezKzkBjnJ93VmqGYT1uaO+5CkkNSRwfbNzDqlMnZNkaUf+pbkC05Yi6StgZOBW3sboZlZc52cyLsb+tmSbzphbfIF4GzghV4GZ2ZWRItzvPZcP5N80wlrJe0NTIuI6xudSNLcysS4Tz7X+UDNzBpxS34UJI0DvgR8utmxETE/IkYiYmTyVt2PzcysYtBb8v18uqbZhLVbA3sAN6UprnYGFko6LCKW9ixKM7MGPKxBfQ0nrI2Ip4FJlW1JNwGfcYI3s0EzyF+G6lt3TUSsByoT1t4HXFWZsFZS3UlpzcwGyaA/XdPXL0M1m7C2qvzAXsRkZtaqQf4ylL/xambWBo8nb2ZWcm7Jm5mVVGVYg0HlJG9m1gaPQmlmVnLukzczKym35M3MSsxJ3sys5NxdY2ZWUn66xsysxNxdY2ZWck7yZmYl5WENzMxKzi15M7OS8o1XM7MS841XM7OSG+Q+eUVEv2PoKEnPAiv6HUcNk4C1/Q6iBsfVmkGMaxBjguGIa7eImNzOySR9n9xUpU2sjYg57VyvVWVM8ksjYqTfcVRzXK1xXMUNYkzguAZF3+Z4NTOz7nOSNzMrsTIm+fn9DqAOx9Uax1XcIMYEjmsglK5P3szMNihjS97MzBIneTOzEitVkpc0R9IKSaskndqF80+TdKOk5ZKWSTo5le8g6YeSVqbX7VO5JP1TiuduSXvnznVsOn6lpGNz5ftIuie9558kqYX4xku6Q9J1aXuGpFvTua6UtGkq3yxtr0r7p+fOcVoqXyHpfbnyUdWtpO0kXSPpfkn3Sdp/EOpL0qfSv+G9kr4laWI/6kvSJZKekHRvrqzr9VPvGg1iOjf9G94t6TuSthttHYymnuvFldv3aUkhaVIv62ooREQpFmA88Atgd2BT4C5gVoevMQXYO61vDTwAzALOAU5N5acCZ6f1Q4DvAQL+ALg1le8ArE6v26f17dO+n6djld57cAvxnQJ8E7gubV8FHJnWLwQ+ntZPBC5M60cCV6b1WaneNgNmpPoc307dAguAj6T1TYHt+l1fwC7Ag8DmuXo6rh/1BbwT2Bu4N1fW9fqpd40GMb0XmJDWz87F1HIdtFrPjeJK5dOAG4BfApN6WVfDsPQ9gI79ILA/cENu+zTgtC5f81rgILJv2E5JZVOAFWn9IuCo3PEr0v6jgIty5RelsinA/bnyjY5rEstUYDHwbuC69EFdm/uP+Wr9pP8Q+6f1Cek4VddZ5bjR1i2wLVkyVVV5X+uLLMk/kv6jT0j19b5+1RcwnY0Tatfrp9416sVUte9PgMtr/WzN6mA0n8tmcQHXAG8HHmJDku9ZXQ36Uqbumsp/3Io1qawr0p+SewG3AjtFxONp16+AnZrE1Kh8TY3yIr4M/A0bhtHYEfjviFhf41yvXj/tfzod32q8zcwAngT+VVk30sWStqTP9RURjwL/B3gYeJzs57+N/tdXRS/qp941ijierKU7mphG87msS9LhwKMRcVfVrkGpq74rU5LvGUlbAf8GfDIinsnvi+zXfU+fS5V0KPBERNzWy+sWMIHsz+uvRsRewG/J/tx9VZ/qa3vgcLJfQq8HtgR6Op5IUb2on1auIWkesB64vJsxFYxlC+BvgdN7dc1+fF7bVaYk/yhZ31zF1FTWUZI2IUvwl0fEt1PxryVNSfunAE80ialR+dRR/AwHAIdJegi4gqzL5h+B7SRVRhrNn+vV66f92wJPjSLeZtYAayLi1rR9DVnS73d9zQYejIgnI+Il4Ntkddjv+qroRf3Uu0Zdko4DDgU+nJLdaGJ6itbruZ43kP2ivit99qcCt0vaeRRxdbSuBkq/+4s6tZC1GleT/aNXbvS8tcPXEPB14MtV5eey8Y2Zc9L6+9n45s/PU/kOZH3V26flQWCHtK/65s8hLcZ4IBtuvF7Nxje4Tkzrf8XGN7iuSutvZeObaKvJbqCNum6BnwBvSutnpLrqa30B+wHLgC3S+xYAf92v+uK1ffJdr59612gQ0xxgOTC56riW66DVem4UV9W+h9jQJ9+zuhr0pe8BdPSHye6oP0B2V39eF87/R2R/qt0N3JmWQ8j6DRcDK4Ef5T40Ai5I8dwDjOTOdTywKi1/mSsfAe5N7zmfqhtPBWI8kA1Jfvf0wV2V/mNtlsonpu1Vaf/uuffPS9deQe5JldHWLbAnsDTV2XfTf6y+1xfweeD+9N5vkCWpntcX8C2y+wIvkf3lc0Iv6qfeNRrEtIqsL/vOtFw42joYTT3Xi6tq/0NsSPI9qathWDysgZlZiZWpT97MzKo4yZuZlZiTvJlZiTnJm5mVmJO8mVmJOcmbmZWYk7xZF0naUtLSNPSEWc85yVtTysaEPzG3/dMeXPMMSZ8pcNxzHbjWdEm/k3Rnrqzm2OX1xkhv4HNkQ+uONrbNJd0paV1lrHSzVjjJWxHbkY31DUBE/GH/QumaX0TEnrntS6katEzSeLJvUR5MNo76UZJmpX2/L+m6quUIsqEARj3WSUT8LsX12GjPYWObk7wVcRbwhtSiPLfSek4t4PslXSrpAUmXS5ot6eY0i86+lRNIOkbSz9M5LkoJcyOS5qXz/Bfwpqp9Rd7/XUm3KZvxaW4qO1PSJ3PH/IPSjF6NRMSPgf9XVbwvsCoiVkfEOrLB4A5Px98TEYfmF+BtZGOhHA18VNK4XJ1drmymrGskbZG6da6XdJey2aqOaBajWRFO8lbEqaSWbkR8tmrfG4HzgDen5WiyMX4+QzYMLJLeAhwBHJBapS8DH86fRNI+ZINS7Uk25sk7cvuavj85PiL2IRuD5BOSdgQuAf4inWdcusZlo6gDaHGc+IiYFxGfJJut618iojLW/5uAf46ItwDPkP2VNAd4LCLeHhF7AN8fZYxmG5nQ/BCzhh6MiHsAJC0DFkdESLqHbMRAgPcA+wBL0rSZm/PaLow/Br4TEc+ncy3M7SvyfsgS+5+k9WnAzIi4RdJTkvYim+zhjohoNHxtx0XEpVVFj0TEzWn9MuATwELgPElnkw0w95Mehmgl5iRv7Xoxt/5KbvsVNny+BCyIiNNGeY2m75d0INk48ftHxPOSbiIb1RDgYrI5XHcma9mPVqfmLKgeFTAi4gFlk00fAnxR0uKIOHOUcZq9yt01VsSzZBOXj9Zi4E8lvQ5A0g6Sdqs65sfAB9PTJFsDH2jx/dsCv0kJ/s1kfeEV3yHrDnkH2Tyio7UEmClphqRNybp+FjZ5Ty27Sto/rR8N/Jek1wPPR8RlZOOX791GnGavckvemoqIp9LN1HvZMLdnK+9fLunvgB+kfvGXyCaI+GXumNslXUk2ucQTZAm18PvJ+rA/Juk+snHNb8m9f52kG8nmFn25SMySvkU2Nv8kSWuAv4+Ir0k6iewXxXjgkohY1mJ1kOL7K0mXkD1981Wy7qpzJb2Sfr6Pj+K8Zq/h8eSt9NIvhtuBP4uIlTX2TyfrB9+jB7GM6lrKprcbiYi13YjLysvdNVZq6Tn2VWQ3hF+T4JOXgW3zX4YaFJUvQwGbkN3nMGuJW/JmZiXmlryZWYk5yZuZlZiTvJlZiTnJm5mVmJO8mVmJOcmbmZWYk7yZWYk5yZuZlZiTvJlZif1/wRdLGAGsjngAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "color_data = mm_rs[:, :, 2].T\n",
    "plt.imshow(color_data, cmap='hot', aspect='auto')\n",
    "plt.ylabel(r'grain position')\n",
    "plt.xlabel(r'time delay [$10^{-4}$ps]')\n",
    "plt.colorbar()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "d4ef83c6-9075-4067-941e-5147345e0bad",
   "metadata": {},
   "outputs": [],
   "source": [
    "pt, pmx, pmy, pmz, pm2=get_ms_paul()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "4da0b689-8f15-4775-96de-caaa3cdb55ef",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEGCAYAAABo25JHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/av/WaAAAACXBIWXMAAAsTAAALEwEAmpwYAAAkVUlEQVR4nO3df5xcdX3v8dc7m1+LJKSQQGUhBjAEUxKJRH5IraBwE6GFlFoFpUVK4arFVtFcwwX5JdzgzUNq9eKPqBR/UAEtjVGQaPkhFIESDElINBD5IVlUIpIESUg2y+f+cc5ulmFn58zunJk5O+/n4zEP9pw5M/M55Jz5zPe3IgIzM2tdIxodgJmZNZYTgZlZi3MiMDNrcU4EZmYtzonAzKzFjWx0ANWaOHFiTJkypdFhmJkVykMPPfS7iJjU33OFSwRTpkxh+fLljQ7DzKxQJD1V7jlXDZmZtTgnAjOzFudEYGbW4pwIzMxanBOBmVmLK1yvIbMsLlqymm/d/6uGfPZIwfqFJzXks80Gw4nABmXKglsaHULT2hnw+gtucTIokCUrOvnIjQ83OozMxo9pY9Vlc2v2fk4E1uv1F9zCTs9KXhP+/1g/rfijZMv2bmZeclvNkoETQQtqxRvHmp+vy+ps2d5ds/dyIhjmfHNZMzjh6rt47NkXGx2GleFEMMz4i98ayddfMTkRDAO++awRfN011vgxbTV7LyeCgvJNWL0nr6p9L55W+Xd431fu495f/r7RYVjKvYZaXDN/8eTxRWuN08zXWj200vXsRFAQjbopW+lmsOHz5e/rtjpOBE2uHjfmMQftyfXnHJ3751hzauYvf3+h14cTQZPK6+YU8IRvrpoZIXi5n8FjI1T/WKrRDCNp/SXfPJwImlCtk4BvuPz0lwQG2t9ojRg97uuv+TkRNJFaJgDffPUxoX0Um7Z19bu/mRx55Y/57Qs7cv2Mz77nMObN6sj1MywfTgRNohZJwF/+9acyVUDl9tdbnlVAvt6GDyeCBqtF/2zfkI2zaeurSwMD7a+nWlcxenrt4cuJoIGGcqP6pmwOe5SpGtqjgVVDtVyL4YyjJnPFvBk1eS9rXk4EDTKUJOASQPNotqqhWpUCfI21FieCBhjszeqbs/k0U9XQUJOAr6/W5URQZ4O5Wd0bo3ntO6Gdzk3b+t1fTy5h2lA4EdTRYG5W36TN7bhDJvVbH3/cIZPq8vlDmeff15b1cCKok2qTgBuDi+HOX2ysan8tDXZwmBOAlXIiqINqk4Bv1OJ4pp9qoYH214qTgNWSE0HOnASGt0a0Ecy85Laqk4CvKxvIiEYHMJw5CQx/5doC8mojuGjJ6qoXLfd1ZZU4EeTk9Rc4CbSCercRVDNQTPi6smycCHKwZEVnVUV336zFVc82gmpKmGccNdnTjVtmuSYCSXMlrZO0XtKCfp6fLOlOSSskrZJ0Yp7x1Es1k3w5CRRbubaAWrcRVJMEnrzqJE8LYVXJLRFIagOuAd4JTAdOlzS95LCLgJsiYhZwGvCFvOKpl2pvWCu2+XOmMapkFZpRI8T8OdNq9hmHXHhr5mN9Tdlg5FkiOAJYHxGPR8QO4AbglJJjAhif/r0H8EyO8eTOSaBFlc4rVON5hl7qzlbP6GvKBivPRNABPN1ne0O6r69LgTMkbQBuBT7c3xtJOlfScknLN27Mf6DOYLzvK/dlPtY37PCxaNk6ukq+qLu6g0XL1tXk/bP+uDjmoD1r8nnWmhrdWHw6cF1E7AecCHxT0qtiiojFETE7ImZPmlSfofvVyrqmgJPA8JJnY3E1VULXn3P0kD/PWleeiaAT2L/P9n7pvr7OBm4CiIj7gLHAxBxjykXWX237jBudcyRWb3k2FrtKyOolz0TwIDBV0gGSRpM0Bi8tOeZXwDsAJL2BJBE0Z91PGdWMF3jgwhNyjMQaIa/G4qw/LpwErBZySwQRsRM4D1gG/Jykd9AaSZdLOjk97GPAOZJWAt8G3h8Rg5hBpXGyjhfwDTuM5dxYXI7bBaxWVLDvXWbPnh3Lly9vdBiAf7UZHHPVHf3ONdQxoZ17F7x9UO/p68ryIOmhiJjd33ONbiwurKwNeW4XGN4aNfuok4DVkhPBIGVtyHO7wPBW68biLKWB8WPaBvXeZuVUTASSvpllXytx0d16zJ8zjfZRr/xibh/VVtORxaVWXTY3t/e21pRlPYI/6buRTh1xeD7hNL+LlqzOdJwb8lpDz1rSl3xvDZtf6uK148fyiXceMqg1prP8wPjsew6r+n3NKilbIpB0gaQXgJmStqSPF4Bnge/VLcImk3UaYA/waTF16ik0mARjVknZRBARCyNiHLAoIsanj3ERsVdEXFDHGJtG1jEDrhJqHUtWdHLBzavZvK0LgF9veYkLbl7NkhWlYycHluXa8nVleanYRhARF0jqkPQWSX/W86hHcM0my5gB9xJqLYuWrWNb1ytXDNvW1V31XEODWX/YrFYqthFIuopkVPBaoOeKD+DuHONqOlkbiN1LqLXUovtoltLDGUdNzvx+ZtXK0lj8l8C0iNiedzBF56J766nF4vVZFjLyQjOWpyzjCB4HRuUdSDOrdhF6ax316D7qcQOWtywlgq3Aw5JuB3pLBRHxj7lF1USydhd1aaA19fTiuXTpGjZt6+KPx49lQRXdR7OsY+FxA5a3LIlgKa+eNbRlZOku6gZiG6ys61iY5aliIoiIr0tqByZHRG2WXSqII6/8cabj3EDcunq6j/b0HPpN2n0UatPn3yVNq4csU0z8BfAwcFu6fZikligh/PaFHRWP8UjP1jaU7qMnXH1XTlGZVSdLY/GlJAvRbwKIiIeBA3OLqEnMvOS2TMd5pGdrG0r30ceefXHA50fWabSyWZZE0BURm0v2vZxHMM1ky/buise42G55LlW5fqGvL6uPLIlgjaT3Am2Spkr6PPDTnONqqGoWDbfW1ojZR81qLUsi+DDJDKTbSZaT3AJ8JMeYGi7LWgMuDRgkVYMLT53BhPZkqM0+48aw8NQZFasMD6gwNmVsm+uFrH6y9BraClyYPoa9ahajNxusSj81fnHliXWJwwyyzTV0MPBxYErf4yNicAuyNrksk3+5NGA9SruP/vaF7TXtPmpWD1kGlH0H+BLwVXZNOjcsuTRg1Rqo+2i5RJB1tLpZvWRJBDsj4ou5R9IEXBqwag2m+2il0eq+xqzesjQWf1/ShyS9VtKePY/cI6uzLD2F3K/bSuXZfdSsXrIkgjOB+SRdRh9KH8vzDKoRsvQUcr9uK+XuozYcZOk1dEA9AmmkLHMKuTBg/elpB7js+2t4fmsX+4wbwwUnvmHQDcWectoaIUuvoVHAB4Ge5SnvAr4cEV05xlVXWeYUesL1tlZGz5f+R258mG+fexQHTtq97LGVqiA95bQ1QpbG4i+SLEzzhXT7b9J9f59XUPWUZT54lwasEmW8SLJUQZrVW5Y2gjdHxJkRcUf6OAt4c96B1UuW+eBdGrCBLFnRyaVL1wBw2uL7M61BbNZMsiSCbkkH9WxIOpBhMp7AN6wNVc+Asue3JjWlz6YDygZzbXlaCWuULIlgPnCnpLsk/QS4A/hYvmHVR5ZFw92n2wZSzXoElQaSeVoJa5QsvYZulzQV6OkPty4itg/0GrNWUc2AsizLnpo1QpZeQ2OBDwF/SjJX1j2SvhQRL+UdXJ6mVJj9EVwasMr2ndBOZz9f+h5QZkWSpWroGyTTUH8e+H/p39/MMyizovCAMhsOsnQfPTQipvfZvlPS2rwCqgeXBqxWesYQXP79tfx+6w72HjeG/z2IAWX7jBudR3hmmWRJBD+TdFRE3A8g6UgyTjEhaS7wL0Ab8NWIuKqfY95Nsi5yACsj4r0ZY89syYpOFi1bxzObtrnIbjU3b1YHo0eO4EPX/4xvnn0k0/54XNXv8cCFJ+QQmVk2WRLB4cBPJfW0dE0G1klaDUREzOzvRZLagGuAE4ANwIOSlkbE2j7HTAUuAI6JiOcl7T2Ec+lX6Xzx/dXnljrjqMm1DsNaRJRZcibLNCZmjZIlEQx2zPsRwPqIeBxA0g3AKUDfaqVzgGsi4nmAiHh2kJ9VVn/d+yq5Yt6MWodhw1zPCIAoM3A4yzQmZo1SsbE4Ip4iWad4D2CvnkdEPJU+V04H8HSf7Q3pvr4OBg6WdK+k+9OqpFeRdK6k5ZKWb9y4sVLIrzDQvPD9mbr3a6o63mzJik4u/I9HADjz2v/2QEUrnCzdRz8FvB/4JbuWWg2gFktVjgSmAscC+wF3S5oREZv6HhQRi4HFALNnz65qspZy3fvK+fH5x1bz9tbiSqsen/VSlVZAWbqPvhs4KCKOjYjj0keWJNAJ7N9ne790X18bgKUR0RURTwCPkiSGmumve185HuJv1apmZHE5n33PYTWOyqw6WRLBI8CEQbz3g8BUSQdIGg2cBiwtOWYJSWkASRNJqooeH8RnlTVvVgcLT81W5+8h/latwSxVWcolB2u0LIlgIbBC0jJJS3selV4UETuB84BlwM+BmyJijaTLJZ2cHrYMeC4dl3AnMD8inhvcqZTnG83ykmWpyhOuvqtO0ZgNTpZeQ18HPg2sBl6u5s0j4lbg1pJ9F/f5O4Dz00duPIDM8jJ/zrRXtBHAq0cWP/bsi40IzSyzLIlga0R8LvdIzAqod2TxD9by+xd3MGn3MVx40uCXqjRrhCyJ4B5JC0nq93tnHY2In+UWVQ3NvOS2ise4NGBDMW9WB7uPGcnff2M5X3v/bGbuN6HRIZlVJUsimJX+96g++2rVfTR3W7YPizV0rMn1LFVZbkBZOe4xZM0gy3oEx9UjkEY55qA9Gx2CDQNZ1ywu5SokawYVew1J2kPS1T0jeyV9RtIe9QiuHq4/5+hGh2AFt2RFJ/O/swqAv//6co8stsLJ0n30WuAFkoFl7yaZbuJf8wyqlsaPKT+YzAPIbKh6RhY/92Iyl9DGP7xyzeL3feW+RoZnlkmWRHBQRFwSEY+nj8uAA/MOrFZWXTa332Qwtk0eQGZDVmlk8b2//H0jwjKrSpbG4m2S/jQi/gtA0jFAdTO5NdiqywY7garZwGoxstis0bIkgg8CX+/TLvA8ySR0Zi1vKGsWZymOm9VDlmmoH46INwIzgZkRMSsiVuYfmlnzG8qaxVe766g1iSy9hv6PpAkRsSUitkj6I0lX1CM4s2bXM6nhXq9J1hyeuPtoFp46I1O3UHcdtWaRpXT6zr7rA6SribmV1Sw1b1YH/5z+uv/SGYf7C94KJ0siaJM0pmdDUjswZoDjzVpO78jixoZhNihZGouvB26X1DN24CySGUnNLKV01eK+U0x4YJkVRZYpJj4taSVwfLrrUxGxLN+wzIqlvykmPv4d96mwYshSIiAibgMqT+Np1uKiT5Fg58uuKLJiyJQIzGxg967/HQDvWXw/HRPaM3UfNWsWHtNiNkRLVnTy1Xue6N3u3LSNC25ePeBrzjhqct5hmWXmRGA2RIuWrWNH9ytXcS2df6jUFfNm5BmSWVUqVg2lcwtdCrwuPV4kyw0XZuI5szx5XiEruixtBF8DPgo8BHi5L7MS5eYbMiuKLFVDmyPihxHxbEQ81/PIPTKzgpg/Zxqj2155K5XOP2TWzLIkgjslLZJ0tKQ39Txyj8ysIObN6uB/vm1XTWnHhHYWnuo2ACuOLFVDR6b/nd1nX2EWrzerh7dOncTn71jPt84+kj+dOtGjiq1QWn7xerNa2DXXUDKI7BP/vqqB0ZhVp+UXrzerhZ4ZJnoGFm/f+XLZY82azbBfvN6sHvqba8isKLK0ERwUEX/VZ/sySQ/nFI9ZoWWZXcijiq3ZtMTi9WZ5u/vRZK6hM6/9bzoqrFfsUcXWbKpdvF7A7/Hi9Wa9lqzo5Et3/bJ324PLrGiy9Bp6GHijpPHp9pa8gzIrkkXL1rG9243DVlxlE4GkMyLiW5LOL9kPQERcnXNsZoXguYas6AYqEbwm/e+4fp7zihtmKc81ZEVXNhFExJfTP/8zIu7t+1zaYGxmJHMNfeK7q1w9ZIWVZRzB5zPuexVJcyWtk7Re0oIBjvsrSSFpdrljzJrVvFkd/MPbX9+7XanXkFmzGaiN4GjgLcCkknaC8UDFqRUltQHXACcAG4AHJS2NiLUlx40D/gl4oPrwzZrDsdMmcfWPH+Wrfzub46fvw5QFtzQ6JLPMBioRjAZ2J0kW4/o8tgDvyvDeRwDrI+LxiNgB3ACc0s9xnwI+DbxURdxmTUXpJBMBnnDOCmegNoKfAD+RdF1EPDWI9+4Anu6zvYFdM5kCkE5nvX9E3CJpfrk3knQucC7A5MkelWnNp3fSuQgu/I+B1ys2azZZBpRtlbQI+BNgbM/OiBjSNNSSRgBXk2FwWkQsBhYDzJ492z2WrKm9uMML+VmxZGksvh74BXAAcBnwJPBghtd1Avv32d4v3ddjHHAocJekJ4GjgKVuMLYiq/QrpX1UllvOrL6yXJV7RcTXgK6I+ElE/B3ZFqV5EJgq6QBJo4HTgKU9T0bE5oiYGBFTImIKcD9wckQsr/40zBprV9XQwMctPHVm/sGYVSlLIuhK//trSSdJmgXsWelFEbETOA9YBvwcuCki1ki6XNLJg47YrAlp14oEAx43b1ZH/sGYVSlLG8EV6YRzHyMZPzAe+GiWN4+IW4FbS/ZdXObYY7O8p1kzumvdswB84Fs/a3AkZtXLkghWRsRmYDNwHICkP841KrMCWbKik3+5/bFGh2E2aFmqhp6Q9G1Ju/XZd2vZo81azKJl6zItTekRx9assiSC1cA9wH9JOijd54X5zFJZZh9tH9XG/DnT6hCNWfWyJIKIiC8AHwa+L+kv8OyjZr32rfBLv2NCOwtPneGGYmtaWRKBANIZSN8B/C/gkDyDMiuS+XOmMWbkK2+l9lFtvHaPscz5k324d8HbnQSsqWVJBCf2/BERvyZpMJ6bW0RmBTNvVgcfPf7g3u2eEgDAuLGjGhWWWWYVVygDTu9ZlazE3blFZVYwx0/fm6tu+wWfO30WJ79xXwAuWvII48Zm6Zhn1lheocysJnqXcAWg++XgD9t3ukRgheAVysxqoLTQ/IftOwEY7xKBFUCuK5SZtYreCSbSsvILLyUzs7hqyIogtxXKzFpJaTvaCy8lJQJXDVkRDPRzpXSFsh5ZVygzazmRNp/tSgQuEVjzy7xCmaTdImJrHWMzK4zyVUMuEVjzy9JGsK+ktSSL0yDpjZK+kG9YZsXynz//LQDn37SSY666g7vWbQRg9zEuEVjzy5IIPgvMAZ4DiIiVwJ/lGJNZoSxZ0cmiZet6tzs3beOGB38FwG6j3ZxmzS/TunkR8XTJLi/Kapbqb/bRru6kjsiJwIogS7n1aUlvAULSKOCfSFYcMzMGnn203YnACiBLieADwD8AHSSLzx+WbpsZA88+OrrNi9Vb86t4lUbE7yLifRGxT0TsHRFnRMRz9QjOrAj6m320bYQYM3LEq8YXmDWjilVDkiYB5wBT+h4fEX+XX1hmxTFvVgfPv7iDy36wFkhmH528526s3/iHBkdmlk2Wcuv3gD2A/wRu6fMws9ScQ5NlvK86dQb3Lng7e48f44ZiK4wsjcW7RcQnco/EbBjZuqOb9lFOBFYMWUoEP5B0YuXDzFpXT1NAz/zsL3V1u8eQFUaWRPBPJMlgm6Qtkl6QtCXvwMyKRL3rESTbW3d0u2rICqNi1VBE9LcwjZn1satEkGSCrTu6+aPdPM+QFUOWXkNv6mf3ZuCpiNhZ+5DMiqd00rltO3bSPtrzDFkxZLlSvwC8CVidbs8AHgH2kPTBiPhRXsGZFUbJcIHnt3Zx+89/ywELbmHfCe3MnzONebM6GhObWQVZ2gieAWZFxOERcTjJyOLHgROA/5tjbGaFEyST0G3e1sXWHd0EySR0F9y8miUrOhsdnlm/spQIDo6INT0bEbFW0iER8bhHTZollj3yGwA+ueQR2vq5L7Z1dbNo2TqXCqwpZUkEayR9Ebgh3X4PsFbSGKArt8jMCmLJik6uvGXXPIzdPQ0FJQaanM6skbJUDb0fWA98JH08nu7rAo7LJyyz4li0bB0vlUxD3Z+BJqcza6Qs3Ue3AZ9JH6U8mYq1vCy/9NtHtTF/zrQ6RGNWvSzdR6cCC4HpwNie/RFxYI5xmRXGvhPa6ewnGajP8+41ZM0sS9XQvwJfBHaSVAV9A/hWljeXNFfSOknrJS3o5/nzJa2VtErS7ZJeV03wZs1g/pxpjB356lvptCP254mrTuLeBW93ErCmliURtEfE7YAi4qmIuBQ4qdKLJLUB1wDvJClNnC5peslhK4DZETET+C7ujmoFNG9WB5/8812X9t7jxgBw1IF7NSoks6pkSQTbJY0AHpN0nqS/BHbP8LojgPUR8XhE7CDpdXRK3wMi4s6I2Jpu3g/sV0XsZk3jz2fuC8An/3w61511BMCrFqsxa1ZZJ53bDfhH4HDgb4AzM7yuA+i76P2GdF85ZwM/7O8JSedKWi5p+caNGzN8tFmd9cw1FMGO7qQH0WgnAiuILL2GHkz//ANwVh5BSDoDmA28rUwMi4HFALNnz+6/k7ZZA/UdQ7Yj7Uo6us2zj1oxZOk1NBu4EHgdr1yqcmaFl3YC+/fZ3i/dV/r+x6fv/7aI2J4hZrOm1psIXCKwgsgysvh6YD7JpHOVR83s8iAwVdIBJAngNOC9fQ+QNAv4MjA3Ip6t4r3Nmkrf2Ud3dHcDTgRWHFkSwcaIWFrtG0fETknnAcuANuDaiFgj6XJgefqei0ganr+Tzlv0q4g4udrPMmu0nnm3guhTNeREYMWQJRFcIumrwO1Ab9VNRNxc6YURcStwa8m+i/v8fXz2UM2aV98SwfbeqiFPymjFkCURnAUcAoxiV9VQABUTgVmrcGOxFVmWRPDmiPAkKWYZBNDVnXRscxuBFUWWRPBTSdMjYm3u0ZgV1A9W/hqAq374C/ZoT9YqdiKwoshypR4FPJzOGbRK0mpJq/IOzKwolqzo5OLvPdK7vXlbskzHj9b8plEhmVUlS4lgbu5RmBVYufUIPnf7Y5x2xOQGRGRWnSwji5+qRyBmRVVuPYJnNr9U50jMBseVmGZDVG7lsQ6vSGYF4URgNkTl1iPwimRWFE4EZkM0b1YHV/7ljN7t3ceMpG2EvBiNFYYTgVkNnHJYsh7B+ScczHuPnMyoNo8qtuJwIjCrgd65hgK6Xw7a5ERgxeFEYFYDfb/2u18ORoxwIrDicCIwq6EgeDmCNicCKxAnArMaUO9Sla4asuJxIjCrgV3rEcDL4aohKxYnArNainCJwAony1xDZlbBkhXJctyfu2M9u41qY8wo/8ay4vDVajZES1Z0csHNq3u3t3Z1s2lrV29yMGt2TgRmQ7Ro2Tq2dXW/Yl+k+82KwInAbIjKzj5aZr9Zs3EiMBuicrOPlttv1mycCMyGaP6cabSPeuVC9cKzj1pxuNeQ2RD1zDL6kRsfBmDsyBHstfsYzz5qheESgVkNzJvVwWtGt/F3xxzA0QftxZ6vGd3okMwyc4nArEZ2HzuSF7fvpDuSBeyPueoOntm0jX0ntDN/zjSXEKxpORGY1ci4saN4YXsXv9m8jQ3Pb+XlSPZ3btrWO87AycCakauGzGpgyYpOnnruRW5d/Rse/e0fepNAj21d3R5XYE3LicBsiHpGFnd1x4DHeVyBNSsnArMh6m9kcX88rsCalROB2RBl+aXfPqrN4wqsaTkRmA1RuV/6bRICOia0s/DUGW4otqblXkNmQzR/zjTmf2clXSUtxJ959xv95W+F4BKBWS30sw7NZd9fwwELbuGYq+7wlNTW1FwiMBuiRcvW9dtj6PmtXYDHEVjzy7VEIGmupHWS1kta0M/zYyTdmD7/gKQpecZjlofODI3FHkdgzSy3RCCpDbgGeCcwHThd0vSSw84Gno+I1wP/DHw6r3jMGi1LwjBrhDxLBEcA6yPi8YjYAdwAnFJyzCnA19O/vwu8Q/Kq32Zm9ZRnIugAnu6zvSHd1+8xEbET2AzsVfpGks6VtFzS8o0bN+YUrplZaypEr6GIWBwRsyNi9qRJkxodjpnZsJJnIugE9u+zvV+6r99jJI0E9gCeyzEmMzMrkWcieBCYKukASaOB04ClJccsBc5M/34XcEdEDDxzl1mTefKqk2p6nFm95TaOICJ2SjoPWAa0AddGxBpJlwPLI2Ip8DXgm5LWA78nSRZmheMveSuyXAeURcStwK0l+y7u8/dLwF/nGYOZmQ2sEI3FZmaWHycCM7MW50RgZtbinAjMzFqcitZbU9JG4KlBvnwi8LsahlMEPufW4HNuDUM559dFRL8jcguXCIZC0vKImN3oOOrJ59wafM6tIa9zdtWQmVmLcyIwM2txrZYIFjc6gAbwObcGn3NryOWcW6qNwMzMXq3VSgRmZlbCicDMrMUNy0Qgaa6kdZLWS1rQz/NjJN2YPv+ApCkNCLOmMpzz+ZLWSlol6XZJr2tEnLVU6Zz7HPdXkkJS4bsaZjlnSe9O/63XSPq3esdYaxmu7cmS7pS0Ir2+T2xEnLUi6VpJz0p6pMzzkvS59P/HKklvGvKHRsSwepBMef1L4EBgNLASmF5yzIeAL6V/nwbc2Oi463DOxwG7pX9/sBXOOT1uHHA3cD8wu9Fx1+HfeSqwAvijdHvvRsddh3NeDHww/Xs68GSj4x7iOf8Z8CbgkTLPnwj8EBBwFPDAUD9zOJYIjgDWR8TjEbEDuAE4peSYU4Cvp39/F3iHJNUxxlqreM4RcWdEbE037ydZMa7Isvw7A3wK+DTwUj2Dy0mWcz4HuCYingeIiGfrHGOtZTnnAManf+8BPFPH+GouIu4mWZ+lnFOAb0TifmCCpNcO5TOHYyLoAJ7us70h3dfvMRGxE9gM7FWX6PKR5Zz7OpvkF0WRVTzntMi8f0TcUs/AcpTl3/lg4GBJ90q6X9LcukWXjyznfClwhqQNJOuffLg+oTVMtfd7RbkuTGPNR9IZwGzgbY2OJU+SRgBXA+9vcCj1NpKkeuhYklLf3ZJmRMSmRgaVs9OB6yLiM5KOJln18NCIeLnRgRXFcCwRdAL799neL93X7zGSRpIUJ5+rS3T5yHLOSDoeuBA4OSK21ym2vFQ653HAocBdkp4kqUtdWvAG4yz/zhuApRHRFRFPAI+SJIaiynLOZwM3AUTEfcBYksnZhqtM93s1hmMieBCYKukASaNJGoOXlhyzFDgz/ftdwB2RtsIUVMVzljQL+DJJEih6vTFUOOeI2BwREyNiSkRMIWkXOTkiljcm3JrIcm0vISkNIGkiSVXR43WMsdaynPOvgHcASHoDSSLYWNco62sp8Ldp76GjgM0R8euhvOGwqxqKiJ2SzgOWkfQ4uDYi1ki6HFgeEUuBr5EUH9eTNMqc1riIhy7jOS8Cdge+k7aL/yoiTm5Y0EOU8ZyHlYznvAz4H5LWAt3A/IgobGk34zl/DPiKpI+SNBy/v8g/7CR9mySZT0zbPS4BRgFExJdI2kFOBNYDW4GzhvyZBf7/ZWZmNTAcq4bMzKwKTgRmZi3OicDMrMU5EZiZtTgnAjOzBqo0ydwg3u/Tkh5JH+/J8honAjNA0qWSPj7UYzJ+1nWSnpD0gQGOeWs6g2hNvhysqV0H1GQqEEknkUxYdxhwJPBxSeMHfBFOBGaNMj/tE96viLiHpK+4DXP9TTIn6SBJt0l6SNI9kg7J+HbTgbsjYmdEvAisIkOScSKwliXpQkmPSvovYFqf/RVvQknnSHpQ0kpJ/y5pN0nj0l/6o9JjxvfdHiCOv06L8Ssl3V3zE7UiWgx8OCIOBz4OfCHj61YCc9PrcSLJ9PP7V3jN8BtZbJaFpMNJRpQfRnIf/Ax4KH16MfCBiHhM0pEkN+HbS97i5oj4SvpeVwBnR8TnJd0FnEQy1cNp6XFdFcK5GJgTEZ2SJgzx1KzgJO0OvIVdswAAjEmfOxW4vJ+XdUbEnIj4kaQ3Az8lmWbjPpIR5gNyIrBW9VbgP3rWaJC0NP1v2ZuwxKFpAphAMnXHsnT/V4H/RZIIziJZH6CSe4HrJN0E3DyIc7HhZQSwKSIOK30iIm6mwjUSEVcCVwIoWaHu0SwfaGa79N6EfR5v6Oe464DzImIGcBnJRGdExL3AFEnHAm0RUbGxNyI+AFxEUoR/SFKR18awIYqILcATkv4aepemfGOW10pq67l+JM0EZgI/qvQ6JwJrVXcD8yS1SxoH/AVUdROOA36d1v+/r+S5bwD/BvxrlkAkHRQRD0TExSTF+Yp1ujZ8pJPM3QdMk7RB0tkk19TZklYCa+h/9b3+jALuSScdXAyckS6+NSBXDVlLioifSbqRpHHtWZLpjnu8D/iipItIbqwb0uP6+iTwAMkX9wMkiaHH9cAVwLczhrNI0lSSNWhv7+ezbBiLiNPLPFV1l9KIeImk51BVPPuoWY1JehdwSkT8TZnnrwN+EBHfrfA+U9LjDq15kGZ9uERgVkOSPg+8k4HHAGwGPiVpYrmxBJLeStJb6Xe1j9LslVwiMDNrcW4sNjNrcU4EZmYtzonAzKzFORGYmbU4JwIzsxb3/wEeAHt1fk302gAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#plt.plot(delay, mag_map_man[:,:,2][:-1])\n",
    "#plt.scatter(delay, mys)\n",
    "\n",
    "plt.plot(sim_delay, mm_rs[:,:,2])\n",
    "plt.scatter(pt, pmz)\n",
    "plt.xlabel(r'delay [s]')\n",
    "plt.ylabel(r'magnetization component')\n",
    "\n",
    "#plt.plot(new_delay, mag_map_man[:,0,0][:-1])\n",
    "#plt.scatter(delay, mxs)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1f0eeba9-3410-4b6c-a266-ddd31894c0ec",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}