llb_params.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "f5343db9",
   "metadata": {
    "tags": []
   },
   "source": [
    "# LLB parameter implementation"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "218eae7f-8e68-4a6f-9928-fd9f03fefbd2",
   "metadata": {},
   "source": [
    "This code computes magnetization dynamics of a one-dimensional sample consisting of an arbitrary amaount of different mataerials with different parameters within. The equations for LLB and all said materials can be found below in documentation. Importantly, temperature dynamics for optical laser excitement are not included here, rather a temperature map on the same sample dimensions and an arbitrary amount of timesteps has to be imported. \n",
    "\n",
    "The code can be roughly seperated in four parts: \n",
    "\n",
    "The static part is merely the definition of the materials and samplestrucutre (things that only have to be computed once in the whole simulation): \n",
    "\n",
    "I. Defining the parameters for each consituent mateiral of the sample (section 1.1)\n",
    "\n",
    "II. Creating a 1d sample of these constituents and mapping all parameters and functions that rely on them on the scale of the created sample structure (sections 1.2 and 1.3)\n",
    "\n",
    "The dynamical part includes all temperature and magnetization dependence (things that have to be computed at every timestep):\n",
    "\n",
    "III. Defining the temperature and magnetization dependence of the parameters (sections 2.2 and 2.3)\n",
    "\n",
    "IV. Joining all predefined functions to run the simulation and create an output (section 2.5)\n",
    "\n",
    "More detailed information on the steps and explanations for the subsections in between that have not been mentioned so far can be found below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "89486d44",
   "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"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ba92f4df-e74e-4bd3-869b-cf731827d56b",
   "metadata": {},
   "source": [
    "For the attempt to solve the differential euqation with the integrated solver odeint we import two more modules:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "4a803c7a-0ff3-4ee3-ac32-0bec13a51a4c",
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.integrate import solve_ivp\n",
    "from scipy.interpolate import UnivariateSpline"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2cd3d85e",
   "metadata": {
    "tags": []
   },
   "source": [
    "## (Magnetic) material class"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "96d775b3-dd8e-4870-99b6-ef980eb4f5c4",
   "metadata": {},
   "source": [
    "I define a class that should compute and hold all the relevant parameters needed for the LLB computation. This class holds only material-specific parameters, while interactions between different materials should be defined and imported from the Structure class."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "8b03fc2f",
   "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=create_mean_mag_map(self)         # 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"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "e2e61949",
   "metadata": {},
   "outputs": [],
   "source": [
    "    def __str__(self):\n",
    "        return self.name"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "10bc9ddc",
   "metadata": {
    "tags": []
   },
   "source": [
    "### Creation of the mean magnetization map"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "342d2991-e019-4cf8-8ade-8b096811fc81",
   "metadata": {},
   "source": [
    "Here a temperature map of the mean field magnetization $m_{eq}(T)$ is created.\n",
    "This means solving the self consistent equation\n",
    "\n",
    "\\begin{eqnarray}\n",
    "    m_{eq}(T)&=&B_S(m_{eq}, T) \\label{meq_def} \\\\\n",
    "    B_S(m, T)&=& \\frac{2 S+1}{2S} \\coth{(\\frac{2S+1}{2S}\\frac{\\text{self.}Jm}{k_BT})}-\\frac{1}{2S} \\coth{(\\frac{1}{2S} \\frac{\\text{self.}Jm}{k_B T})} \\label{Brillouin},\n",
    "\\end{eqnarray}\n",
    "\n",
    "where\n",
    "\n",
    "\\begin{align}\n",
    "\\text{self.}J=3\\frac{S}{S+1}k_B \\ \\text{self.}T_C\n",
    "\\end{align}\n",
    "\n",
    "is the mean field exchange coupling constant for effective spin self.$S$ and Curie temperature self.$T_C$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "1778a430",
   "metadata": {},
   "outputs": [],
   "source": [
    "    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)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5cb39c84-ab8a-4fc1-85b3-c69b072fc76a",
   "metadata": {
    "tags": []
   },
   "source": [
    "### Creation of the map of derivative of the Brillouin function"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1f59aceb",
   "metadata": {},
   "source": [
    "To later define the longitudinal suszeptibility, we need to define the derivative of the Brillouin function.\n",
    "\n",
    "\\begin{align}\n",
    "\\frac{dB}{dx} = \\frac{1}{4S^2\\sinh^2(x/2S)}-\\frac{(2S+1)^2}{4S^2\\sinh^2(\\frac{(2S+1)x}{2S})}, \n",
    "\\end{align}\n",
    "\n",
    "where $x=\\frac{Jm}{k_BT}$.\n",
    "\n",
    "The derivative will be evalueated at equilibrium $m=m_{eq}(T)$. In the following bit we will define the material-dependent prefactors of both terms."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "add9cc0f-f1ac-4ba7-a7cf-5f6233cb919c",
   "metadata": {},
   "outputs": [],
   "source": [
    "    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"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d918d494-ce90-4182-b83c-e322173fe366",
   "metadata": {
    "tags": []
   },
   "source": [
    "### Calling the mean field magnetization"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5325cc81",
   "metadata": {},
   "source": [
    "In the __init()__ function, the maps of mean field magnetization and derivative of Brillouin function are called and saved with create_ functions. To call use these maps for any array of temperatures we define the get_ functions. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "134b257a",
   "metadata": {},
   "outputs": [],
   "source": [
    "    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"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "14308b33",
   "metadata": {},
   "source": [
    "### Material dependent parameters"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "be18c309-1070-473a-9a01-b5106f557f1b",
   "metadata": {},
   "source": [
    "#### Transverse and longitudinal damping parameters"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d7fdf949-e4dd-406c-b018-369b7af8a1b7",
   "metadata": {},
   "source": [
    "Here we define the material dependent parameters of longitudinal and transversal damping parameters. Their dependence on temperature and mean magnetization will be added later when the sample structure is implemented:\n",
    "\n",
    "\\begin{eqnarray}\n",
    "\\alpha_{\\parallel} &=&_{T<T_c} \\frac{2\\lambda}{S+1}\\frac{1}{\\sinh(2q_s)}\\\\\n",
    "\\alpha_{\\perp} &=&_{T<T_c} \\frac{\\lambda}{m_{eq}(T)}(\\frac{\\tanh(q_s)}{q_s}-\\frac{T}{3T_C}),\n",
    "\\end{eqnarray}\n",
    "\n",
    "\\begin{align}\n",
    "\\alpha_{\\parallel, \\perp} =_{T>T_c} \\frac{2 \\lambda}{3} \\frac{T}{T_C}\n",
    "\\end{align}\n",
    "\n",
    "where\n",
    "\n",
    "\\begin{align}\n",
    "q_s=\\frac{3 T_C m_{eq}(T)}{(2S+1)T}\n",
    "\\end{align}\n",
    "\n",
    "$m_{eq}(T)$ is an argument of the function as get_mean_mag() will be called at every timestep before calling the functions for temperature dependent parameters.\n",
    "As $q_s$ is temperature dependent, we will only compute this part later if a sample structure and the according mean field magnetization profile are created."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "889be63e",
   "metadata": {},
   "outputs": [],
   "source": [
    "    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",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1928def4-bb74-4bb9-a467-64db7f97ef0b",
   "metadata": {},
   "source": [
    "#### Longitudinal Susceptibility"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c16e6ad6-5b6a-4cb1-aae8-e251f1ab2e31",
   "metadata": {},
   "source": [
    "The longitudinal susceptibility in MFA is expressed as\n",
    "\n",
    "\\begin{align}\n",
    "\\chi_{\\parallel}=_{T<T_C}\\frac{\\beta \\mu_{\\rm{at}} B_S'(m_{eq},T)}{1-\\beta J B_S'(m_{eq}, T)}\n",
    "\\end{align}\n",
    "\n",
    "\\begin{align}\n",
    "\\chi_{\\parallel}=_{T>T_C} \\frac{\\mu_{\\rm{at}}T_C}{J(T-T_C)}\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "b7e0416f-4bdb-4e82-aec8-4633d01a0fe2",
   "metadata": {},
   "outputs": [],
   "source": [
    "    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"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f091dc70-9c9e-46fd-870b-debf6de1cdbe",
   "metadata": {},
   "source": [
    "#### Uniaxial Anisotropy strength"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b6a0fc17-5fa1-43b0-8eeb-605110d77b60",
   "metadata": {},
   "source": [
    "The uniaxial anisotropy field along of grain i along the easy (m,n,o)-axis is expressed in the form.\n",
    "\n",
    "\\begin{align}\n",
    "\\mathbf{H}_{i,A_m}=-\\frac{2 K_i(T)}{M_s(T)} \\left( m_{i,n} \\mathbf{e}_n+ m_{i,o} \\mathbf{e}_o \\right),\n",
    "\\end{align}\n",
    "\n",
    "where $K_i(T)=K_0 \\ m_{eq}(T)^{\\kappa}$ scales in a power law with the equilibrium magnetization, $\\kappa$ being a material parameter. Here, only the temperature-indipendent but material_dependent part of the prefactor is computed. This is done, so that later when the sample-structure und temperature profile are defined, the corresponding field for the whole sample can be computed with minimal computational effort (I hope and think so at least). What we thus do here is to compute\n",
    "\n",
    "\\begin{align}\n",
    "    E_{0,anis}=-2K_0\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "8428f5a6-85b8-4de6-bc53-db2fc0857b40",
   "metadata": {},
   "outputs": [],
   "source": [
    "    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": "2426e362-fff0-45b2-972a-3c22edde7600",
   "metadata": {},
   "source": [
    "#### Exchange stiffness"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d05f2508-8762-4fca-a5a1-83bfc380116f",
   "metadata": {},
   "source": [
    "The exchange field is expressed by\n",
    "\n",
    "\\begin{align}\n",
    "H_{i,ex}=\\frac{2A(T)}{M_s(0) m_{eq}^2 \\Delta^2} \\sum_{j \\ neighb \\ i} (m_j-m_i),\n",
    "\\end{align}\n",
    "\n",
    "where $A(T)=A_0 m_{eq}(T)^2$ is the temperature dependent exchange stiffness and $\\Delta$ is the dimension of the cubic magnetic grains. At this point there is nothing to compute here, this block serves only for documentation. \\textcolor{red}{Whats the best way to define grain size in the code?}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "654a0d3a",
   "metadata": {
    "tags": []
   },
   "source": [
    "## Creating some sample structure"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fead8e22-382e-4626-9039-4f0b78a432d9",
   "metadata": {},
   "source": [
    "Apart from defining all the necessary fields and material parameters, the most important part is to properly define the sample structure and simplify computation of the magnetic fields for any given sample. For this, we define a rather complicated sample-composition in the following function get_sample(), consisting of three different 'materials' (defined on the grainsize of the micromagnetic simulation, so some nm in thickness) with different magnetic properties and respective interface properties, that need to be defined in the following. Apart from only the sample-structure, we retrieve also the magnetization coordinates of every magnetic grain, and the positions (indices of the 1d array that defines the sample of magnetic grains) of each of the material. This will help to later call the functions within the material class on an array that contains all magnetic grains of each material.\n",
    "\n",
    "Example: Sample consists of three materials in the following order:\n",
    "\n",
    "\\begin{align}\n",
    "[mat1, mat2, mat3, mat3, mat1]\n",
    "\\end{align}\n",
    "\n",
    "material_grain_indices is a list of M lists (M=number of different materials in the sample), with each of the M lists containing the indices of all positions of this material, so:\n",
    "\n",
    "\\begin{align}\n",
    "\\text{material_grain_indices}=[[0,4],[1],[2,3]]\n",
    "\\end{align}\n",
    "\n",
    "sample_sorter is an index array that brings the flattened version of material_grain_indices back into the order of the sample-structure. This will be used to read out at every timestep the mean magnetization of all grains of each material with a given temperature profile (which minimizers the computation of mmag, only once per material) and then flatten it back into the sample structure to perform time-efficient numpy operations with it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "298b2eb8",
   "metadata": {
    "tags": []
   },
   "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.005, 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_1, mat_2, mat_3]\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(Nickel_1+Cobalt+Iron+Nickel_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(60)\n",
    "    m_phi = np.zeros(60)\n",
    "    m_gamma = np.zeros(60)\n",
    "    return materials, sample, m_amp, m_phi, m_gamma, material_grain_indices, sample_sorter, mat_locator"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f438af9e",
   "metadata": {
    "tags": []
   },
   "source": [
    "### Creating magnetization vectors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "7ccd4597",
   "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": "f4ed8d47-aa5d-461f-b142-9a0419edfb83",
   "metadata": {
    "tags": []
   },
   "source": [
    "### Plotting functions for the $m_{eq}(T)$ maps"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "66b79a15-6458-4aa8-afad-b46c3bf9c2c5",
   "metadata": {},
   "source": [
    "This should really work for any sample preperation. Just to visualize the output of mean field magnetization for different parameters, being spin and critical temperature"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "76af2264",
   "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": "a4b3935d-4f1b-437c-a354-114017329d0c",
   "metadata": {
    "tags": []
   },
   "source": [
    "## Embedding material dependent parameters in the smple structure"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6f8e5a6a-ebff-42b1-8b4f-7cccf4ca0b79",
   "metadata": {},
   "source": [
    "To hopefully reduce the computation time we now define all material parameters and functions on the the sample structure. This minimizes the computations to be made in the dynamical simulation, saving everything possible in numpy-array-format to maximally make use of numpy's quickness."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2603bce5-7be3-4294-9d20-5b4f5eee2835",
   "metadata": {},
   "source": [
    "### Exchange coupling"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3f9e0cb8-64bb-477e-80a0-a46dcd64d8fd",
   "metadata": {},
   "source": [
    "The follwoing peace of code defines a symmetric matrix of dimension len(materials)xlen(materials) (here three, indices a,b,c) between potential neighbouring grains with indices i,j\n",
    "\n",
    "\\begin{align}\n",
    "    \\mathcal{J}_{ij}= \\left( \\begin{matrix}J_{aa} \\ J_{ab} \\ J_{ac} \\\\ J_{ba} \\ J_{bb} \\ J_{cb} \\\\ J_{ca} \\ J_{cb} \\ J_{cc}\\end{matrix} \\right),\n",
    "\\end{align}\n",
    "\n",
    "where $J_{lk}=J_{kl}$ for $k,l \\in [a,b,c]$ \n",
    "and $J_{ll}=3 \\frac{l.S}{l.S+1} k_B \\ l.T_C$ is just the mean field exchange coupling constant for material $l$.\n",
    "\n",
    "Now one can fill the array exch_coup_arr with the respective coupling constants of neighbouring grains as defined in sample. \n",
    "For a sample of five grains\n",
    "\n",
    "$[mat1, mat2, mat3, mat3, mat1]$\n",
    "\n",
    "this would produce\n",
    "\n",
    "$[[0, J_{12}],[J_{21}, J_{23}],[J_{32}, J_{33}],[J_{33}, J_{31}],[J_{31}, 0]]$\n",
    "\n",
    "This function just needs to be called once to create the proper interaction array for the simulated sample."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "051c88b7-08e1-4ab5-889f-eb2941dfd153",
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_exch_coup_sample(materials, sample, mat_loc):\n",
    "    # This function takes as input parameters:\n",
    "    #   (i) the 1d-list of magnetic unique materials in the sample (size M)\n",
    "    #   (ii) the 1d numpy array of the sample, consisting of a material (from class material) at each grain (size N)\n",
    "    # As an output we get a 2d numpy array of dimension Nx2 for coupling each site with its 2 neighbours in the linear chain of grains.\n",
    "\n",
    "    # Define a matrix J_mat of dimension len(materials)xlen(materials) with the exchange coupling constants of mat_i and mat_j at J_mat[i][j]=J_mat[j][i]\n",
    "    J_mat=np.zeros((len(materials), len(materials)))\n",
    "    # fill the diagonal with the mean field exchange constant of each material:\n",
    "    for i, mat in enumerate(materials):\n",
    "        J_mat[i][i]=mat.J\n",
    "    # define the off-diagonals, namely some values for exchange coupling constants of different materials:\n",
    "    J_mat[0][1]=1e-20\n",
    "    J_mat[1][2]=5e-20\n",
    "    J_mat[0][2]=1e-19\n",
    "    # symmetrize the matrix so that also elements [i][j] with i>j can be read out:\n",
    "    for i in range(1,len(materials)):\n",
    "        for j in range(i):\n",
    "            J_mat[i][j]=J_mat[j][i]\n",
    "\n",
    "    # Now we can assign the coupling of each grain to its nearest neighbours by filling the output array with the respective matrix entry:\n",
    "    # Let's define the output array:\n",
    "    ex_coup_arr=np.zeros((len(sample),2))\n",
    "    \n",
    "    # This list can assign the proper matrix elements to the output matrix\n",
    "    for i, grain in enumerate(sample):\n",
    "        if i>0:\n",
    "            ex_coup_arr[i][0]=J_mat[mat_loc[i]][mat_loc[i-1]]\n",
    "        if i<len(sample)-1:\n",
    "            ex_coup_arr[i][1]=J_mat[mat_loc[i]][mat_loc[i+1]]\n",
    "    return ex_coup_arr"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "944f039e-2755-45a6-a7cd-55ddf0255fde",
   "metadata": {},
   "source": [
    "### Exchange stiffness"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0f8acf90-61b9-4d21-b66a-a2be229171fa",
   "metadata": {},
   "source": [
    "Exactly the same procedure for the exchange stiffness. Here the lateral size of magnetic grains is still missing!!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "e19285c3-91c2-42df-83cd-6c1d97e980a2",
   "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])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "862437fb-61f1-4609-96d5-34298e3309ef",
   "metadata": {},
   "source": [
    "### Spin, $T_C$, J, $\\lambda$, $\\mu_{\\rm{at}}, M_S, \\Delta$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "b14b5b33-6e53-4f07-9d09-f01db3d91cc9",
   "metadata": {},
   "outputs": [],
   "source": [
    "def S_sample(sample):\n",
    "    return np.array([mat.S for mat in sample])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "213ad972-6435-4851-ae2a-7e02053f7572",
   "metadata": {},
   "outputs": [],
   "source": [
    "def Tc_sample(sample):\n",
    "    return np.array([mat.Tc for mat in sample])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "2365477a-75e1-4868-9a74-72e2bddc43a1",
   "metadata": {},
   "outputs": [],
   "source": [
    "def J_sample(sample):\n",
    "    return np.array([mat.J for mat in sample])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "9f8fa53d-a7d8-41e6-a2ff-a482ef00da41",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def lamda_sample(sample):\n",
    "    return np.array([mat.lamda for mat in sample])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "43c44310-521b-49bf-86a7-1a6255b5f311",
   "metadata": {},
   "outputs": [],
   "source": [
    "def muat_sample(sample):\n",
    "    return np.array([mat.muat for mat in sample])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "01989706-f113-47c2-957c-2d699a98228a",
   "metadata": {},
   "outputs": [],
   "source": [
    "def Ms_sample(sample):\n",
    "    return np.array([mat.Ms for mat in sample])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "e8ba4261-20c0-4222-9e9c-13798050445b",
   "metadata": {},
   "outputs": [],
   "source": [
    "def Delta2_sample(sample):\n",
    "    return np.array([mat.Delta**2 for mat in sample])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e492a9d0-c58d-4ea0-ab26-98cce3b71d94",
   "metadata": {},
   "source": [
    "### Anisotropy strength, exponent and axis"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0344d57e-e880-44a4-a62b-44c6d7c86fce",
   "metadata": {},
   "source": [
    "We allow for different easy axis of different materials within the sample. This information is saved in the arrays ani_axis_mask_sam (to filter out the easy axis, necessary to compute the scalar anisotropy strength) and ani_axis_inv_sam (to filter out the hard axis, used to define the direction of the anisotropy field later)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "f31e1958-c675-4da9-8e9c-16246d3d0bc7",
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_ani_sample(sample, Ms_sam):\n",
    "    ani_sam=np.divide(np.array([anisotropy(mat) 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"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5d4b2f8e-5e4e-49c1-851e-88626ada7545",
   "metadata": {},
   "source": [
    "### Damping paramteres"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "22755a67-4df2-41e8-8f70-f292fa4f4b3b",
   "metadata": {},
   "outputs": [],
   "source": [
    "def alpha_par_sample(sample):\n",
    "    return np.array([alpha_par(mat) for mat in sample])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "628ae266-6136-49c8-9338-e8852e9a4643",
   "metadata": {},
   "outputs": [],
   "source": [
    "def qs_sample(sample):\n",
    "    return np.array([qs(mat) for mat in sample])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c9f59f82-0907-4e3d-90d8-38272eb181cd",
   "metadata": {},
   "source": [
    "### Longitudinal Susceptibility"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "0684cda3-cd6c-4e73-b404-08e4f3c41017",
   "metadata": {},
   "outputs": [],
   "source": [
    "def dbrillouin_t1_sample(sample):\n",
    "    return np.array([dbrillouin_t1(mat) for mat in sample])\n",
    "\n",
    "def dbrillouin_t2_sample(sample):\n",
    "    return np.array([dbrillouin_t2(mat) for mat in sample])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "dcb22903-6954-4f0a-ae69-2609a0fde6c9",
   "metadata": {},
   "outputs": [],
   "source": [
    "def chi_par_num_sample(sample):\n",
    "    return np.array([chi_par_num(mat) for mat in sample])\n",
    "\n",
    "def chi_par_denomm1_sample(sample):\n",
    "    return np.array([chi_par_denomm1(mat) for mat in sample])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "00fdaad6-56ae-4e4e-a078-b95269ae40ed",
   "metadata": {
    "tags": []
   },
   "source": [
    "# LLB 1d dynamical simulation"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "23ee768d-b14f-4569-a8c1-29f1f8e27ffa",
   "metadata": {},
   "source": [
    "## Overview of dynamical simulation"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "27c6c755-0ab2-4e10-a981-279cc83e950f",
   "metadata": {},
   "source": [
    "Now comes the actual dynamical qLLB simulation. We will import a temperature profile dependent on time and 1d-space, and evaluate the change of magnetization vectors of every grain within the sample at every timestep. The dynamical LLB equations are expressed as follows:\n",
    "\n",
    "\\begin{align}\n",
    "    \\frac{1}{\\gamma}\\frac{d\\mathbf{m}}{dt}=-\\mathbf{m} \\times \\mathbf{H}_{eff}- \\frac{\\alpha_{\\perp}}{m^2}\\mathbf{m} \\times (\\mathbf{m} \\times \\mathbf{H}_{eff}) + \\frac{\\alpha_{\\parallel}}{m^2} (\\mathbf{m}  \\cdot \\mathbf{H}_{eff}) \\cdot \\mathbf{m},\n",
    "\\end{align}\n",
    "\n",
    "with the above defined damping parameters. The three terms describe i) precession at Lamor frequency, ii) transversal damping (conserving the macrospin length) and iii) longitudinal damping (changing macrospin length due to incoherent atomistic spin excitations within the grainsize the macrospin is defined on). The effective magnetic field is the sum of all relevant magnetic interactions:\n",
    "\n",
    "\\begin{align}\n",
    "    \\mathbf{H}_{eff}= \\mathbf{H}_{ext}+\\mathbf{H}_A+\\mathbf{H}_{ex}+\\mathbf{H}_{th},\n",
    "\\end{align}\n",
    "\n",
    "where $\\mathbf{H}_{ext}$ is an external magnetic field, and $\\mathbf{H}_{th}$ is a thermal field, defined as:\n",
    "\n",
    "\\begin{align}\n",
    "    \\mathbf{H}_{th}=_{T<T_C} \\frac{1}{2\\chi_{\\parallel}}(1-\\frac{m^2}{m_e^2})\\mathbf{m}\n",
    "\\end{align}\n",
    "\n",
    "\\begin{align}\n",
    "    \\mathbf{H}_{th}=_{T>T_C} -\\frac{1}{\\chi_{\\parallel}}(1+\\frac{3}{5} \\frac{T_C}{T-T_C}m^2)\\mathbf{m}\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bee49dcc-6a60-43aa-816b-88aa96db792c",
   "metadata": {},
   "source": [
    "## Temperature dependence of parameters"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f15c7868-ab4b-42a9-b872-64ba7124fb45",
   "metadata": {},
   "source": [
    "are all temperature dependent functions. Expecially the mean field magnetization profile for the sample at a given temperature profile, and the corresponding effective magnetic field."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d70ad2e9-54a5-4ed1-816e-a3cfaa19498a",
   "metadata": {},
   "source": [
    "### Profile of the mean field magnetization"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2dce6a7e-1320-4f64-8f50-245af5944b6e",
   "metadata": {},
   "source": [
    "Now we make use of the material-specific seperation of the sample in terms of mat_gr_ind and construct a mean mag map for a given (1d) array of temperatures and an arbitrary sample composition."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "f52ff68f-10ad-4ace-93c6-a177940970d3",
   "metadata": {},
   "outputs": [],
   "source": [
    "def split_sample_T(T, tc_mask, mat_gr_ind, materials):\n",
    "    T_sep_red=[np.array([T[i] for i in mat_ind])/materials[j].Tc for j, mat_ind in enumerate(mat_gr_ind)]\n",
    "    tc_mask_sep=[np.array([tc_mask[i] for i in mat_ind]) for mat_ind in mat_gr_ind]\n",
    "    return T_sep_red, tc_mask_sep"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "32368cac-a3fe-4dff-b313-251166991c98",
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_mean_mag_sample_Ts(tes, under_tc, sample):\n",
    "    mmag_sam_T=np.array([[get_mean_mag(mat, t[i], under_tc[i]) for mat in sample] for i in range(len(tes))])\n",
    "    return mmag_sam_T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "id": "fbcae98f-9835-4024-9ffd-9902bf2fa5fa",
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_mean_mag_sample_T(mat_gr_ind_flat, materials, T_sep_red, tc_mask_sep):\n",
    "    Tc_vals = np.array([mat.Tc for mat in materials])\n",
    "    tc_mask_sep_norm = [np.array(tc_mask) for tc_mask in tc_mask_sep]\n",
    "    mean_mags = np.array([get_mean_mag(mat, T, tc_mask) for mat, T, tc_mask in zip(materials, T_sep_red, tc_mask_sep_norm)])\n",
    "    mmag_sam_T_flat = np.concatenate(mean_mags)[mat_gr_ind_flat]\n",
    "    return mmag_sam_T_flat"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "93217089-6673-4bb6-abe8-5ee3a3176ca7",
   "metadata": {},
   "source": [
    "Now we can make use of numpy library to efficiently compute temperature dependence of effective field, damping parameters and susceptibility."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1b6cb406-0acf-4b4e-a5a5-fd0dbd9c6d78",
   "metadata": {},
   "source": [
    "### Anisotropy strength"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3e24c535-bd6a-416c-b694-c3e524d33cf0",
   "metadata": {},
   "source": [
    "### Exchange stiffness"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "id": "c30e873a-0847-4bf4-9955-2fb989900bac",
   "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))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c92c1ff5-5131-494b-84bb-0964dd1bf47a",
   "metadata": {},
   "source": [
    "scales with $m_{eq}^2$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "id": "4507cc94-e973-45ba-8646-a22a496c4761",
   "metadata": {},
   "outputs": [],
   "source": [
    "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)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9da5d716-8b05-4fe5-8f75-5db270631df4",
   "metadata": {},
   "source": [
    "### Damping parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "id": "888f3be7-b667-46ea-9add-0fcc69404bd7",
   "metadata": {},
   "outputs": [],
   "source": [
    "def qs_sample_T(qs_sam, mmag_sam_T, T):\n",
    "    return qs_sam*mmag_sam_T/T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "id": "1bd31049-2719-417c-a219-e0cd0c11f288",
   "metadata": {},
   "outputs": [],
   "source": [
    "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"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "id": "8c9a3bfd-21d4-4916-b32a-10649b2d04b2",
   "metadata": {},
   "outputs": [],
   "source": [
    "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"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c20fd020-509f-4127-8b41-be76149e7f4d",
   "metadata": {},
   "source": [
    "### Longitudinal susceptibility"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "42240ae3-792e-4aa8-8a7f-f394855ae7d5",
   "metadata": {},
   "source": [
    "The longitudinal susceptibility was split apart into numerator and denomenator above, the Derivative of Brillouin function was taken apart into two terms. Now, as we have a temperature-/ and mean_mag-profile, we can define the temperature dependence of the respective terms and merge them. We start by computing the ratio eta that is argument of the Brillouin function, then compute the full derivative of Brillouin function and ultimately compute the full longitudinal susceptibility."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "id": "3e53fbde-1576-41dc-8f3e-b828754b94e6",
   "metadata": {},
   "outputs": [],
   "source": [
    "def eta_sample_T(mmag_sam_T, J_sam, T):\n",
    "    return J_sam*mmag_sam_T/sp.k/T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "id": "ccdd19b6-1837-4773-899f-3de0c2977bff",
   "metadata": {},
   "outputs": [],
   "source": [
    "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"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "id": "28fde98f-fd77-4399-bd39-48ee0e94303a",
   "metadata": {},
   "outputs": [],
   "source": [
    "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": "e7f1537c-a6f5-404d-85fd-d9b75b767a56",
   "metadata": {},
   "source": [
    "## Magnetization dependent dynamical functions"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5587e0d7-d916-422b-aa69-38d809367c4d",
   "metadata": {},
   "source": [
    "### Anisotropy field"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f98fca24-c47f-423e-adf1-48b0c54f8187",
   "metadata": {},
   "source": [
    "For the anisotropy field we have defined the unique easy axis in the materials class. Note that for different grains the easy axis can be oriented differently. As an output we get the anisotropy field as an array of dimension (len(sample) x 3), just as the magnetization is also saved."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "id": "74410e98-c60b-45c7-8455-aef7dc2cb1a1",
   "metadata": {},
   "outputs": [],
   "source": [
    "def anis_field(anis_sam_T, m, ani_perp_sam):\n",
    "    return anis_sam_T[:, np.newaxis]*(m*ani_perp_sam)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1f3c94c7-315f-4ca0-9bb1-062c8fce53fc",
   "metadata": {},
   "source": [
    "### Exchange field"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "15b2ba80-10d0-4ebd-b49c-5d7a1f85aba9",
   "metadata": {},
   "source": [
    "As an input we alreadz get the difference vectors of neighbouring magnetic grains, so that we just have to multiply is with the corresponding coupling parameter from the (2 x len(sample)) matrix that holds the now temperature dependent coupling constants. diff_up corresponds here to $m_{i-1}-m_{i}$, diff_down to $m_{i+1}-m_{i}$ at each position $i$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "id": "219e7b8f-09b5-4001-9470-2d2cc5026490",
   "metadata": {},
   "outputs": [],
   "source": [
    "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"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0d114768-0ac2-4cdd-90fc-007daa906f7c",
   "metadata": {},
   "source": [
    "### Thermal field"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3f2844f9-871d-461e-a91e-c8c4ce63c984",
   "metadata": {},
   "source": [
    "In defining the thermal fields we must differentiate where the temperature exceeds T_c at any grain and where not."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "id": "73039fac-4327-410f-8488-924ddda4344b",
   "metadata": {},
   "outputs": [],
   "source": [
    "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": "69ccfb25-4007-4f50-878e-d4c71d18fbf4",
   "metadata": {},
   "source": [
    "Next steps: 1. add up fields for effective field, 2. import temperatures and create map, 3. run test, 4. debug, 5. physical check "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3d000b2b-e898-4f7d-b233-b66bada68502",
   "metadata": {},
   "source": [
    "## Import temperature map"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8a308795-14fc-4288-b33d-487eebab3b65",
   "metadata": {},
   "source": [
    "Now I import a precomputed temperature map of sample of Nickel, Cobalt, Iron Nickel. The sample is 60 grains long in total and N timesteps were computed. In the following block I create an electron tempererature map of dimension Nx60. Also, we define the timestep of simulations, which is 1 as for this simulation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "id": "280d391c-f30c-4680-8c56-f203dc76401d",
   "metadata": {},
   "outputs": [],
   "source": [
    "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=np.array([[tes[i,0]] for i in range(len(delay))])*0.7"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "id": "efd4fe1d-55dd-4bd3-b0b0-9620295fa747",
   "metadata": {},
   "outputs": [],
   "source": [
    "#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])\n",
    "#new_delay=np.arange(0,delay[-1], 1e-4)#[:10001]\n",
    "\n",
    "#tes=float_dat[:,5]\n",
    "#new_tes=np.array(ip.interp1d(delay, tes)(new_delay))#[:10001]\n",
    "#new_tes = np.reshape(new_tes,(-1,1))\n",
    "#print(len(new_tes))\n",
    "\n",
    "#mxs= float_dat[:,1]\n",
    "#mys=float_dat[:,2]\n",
    "#mzs=float_dat[:,3]\n",
    "#m2s=float_dat[:,4]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "deecdac2-d76a-4ddf-9081-1d416da48f16",
   "metadata": {},
   "source": [
    "## Run dynamical simulation"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bdf59ec0-c4d1-4596-ac38-183ab01aec30",
   "metadata": {},
   "source": [
    "Now we'll define a function that loops over all N timesteps in the temperature map and creates the corresponding magnetization map as an output. I see now that it is way prettier to define a class for both the sample and the simulation as well, but I won't go back to rewrite this part now, this can be adapted in your code accordingly. Because we'll use Heun method for the computation of each mag increment, I will define a function to compute these increments. Because neither sample structure nor simulation are implemented as a class here, this function takes a lot of arguments: Let's start with the function to compute the increments:"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1da724cb-9809-41ca-8981-d118f5a8c430",
   "metadata": {},
   "source": [
    "### Function to compute magnetization increments"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "id": "812947f7-c523-41c0-95b7-5f100e771d90",
   "metadata": {},
   "outputs": [],
   "source": [
    "def mag_incr(materials, sample, m_amp, m_phi, m_gamma, mat_gr_ind, mat_gr_ind_flat, mat_loc, Ms_sam, ex_stiff_sam, S_sam, Tc_sam, J_sam, lamda_sam, muat_sam, K0_sam, kappa_ani_sam, ani_perp_sam, alpha_par_sam, qs_sam, dbrillouin_t1_sam, dbrillouin_t2_sam, chi_par_num_sam, chi_par_denomm1_sam, m, Te, H_ext, gamma, dt):\n",
    "    #some operations of the magnetization we need:\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",
    "    # at every timestep, we have to calculate the temperature dependent parameters, so let's call all the functions defined before\n",
    "    t_reduced = np.divide(Te, Tc_sam)\n",
    "    under_tc=t_reduced<1.\n",
    "    over_tc=~under_tc\n",
    "    Temp_sep, tc_mask_sep=split_sample_T(Te, under_tc, mat_gr_ind, materials)\n",
    "    mmag_sam_T=get_mean_mag_sample_T(mat_gr_ind_flat, materials, Temp_sep, tc_mask_sep)\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",
    "    # from all these we can define 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",
    "    # and 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=gamma*dt*(-m_rot -trans_damp + long_damp)\n",
    "    \n",
    "    return dm"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7744c820-1ded-47bc-bd4b-e81b1a5b539f",
   "metadata": {},
   "source": [
    "And then the main function, to be called to run the simulation. It calls all the static funcitons to set up parameters and constants on the sample structure and then in Heun method fashion calls the function mag_inr() accordingly. It returns the finished magnetization map for all layers and all times."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1a45f8ba-8859-42c5-859a-01badf49bbf0",
   "metadata": {},
   "source": [
    "### Main function computing the magnetization map"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "id": "6abd0a52-8470-48c4-86ac-45f8b50da404",
   "metadata": {},
   "outputs": [],
   "source": [
    "def run_LLB(tes, dt):\n",
    "    starttime=time.time()\n",
    "    # let's define some constants for simulation:\n",
    "    gamma=1.76e11 # gyromagnetic ratio in (Ts)^{-1}\n",
    "    H_ext=np.array([[0,0,0] for _ in range(60)]) #external field in T\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",
    "    Ms_sam=Ms_sample(sample)\n",
    "    #exch_coup_const_sam=get_exch_coup_sample(materials, sample, mat_loc)\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",
    "    # initialize the starting magnetization\n",
    "    m=get_mag(np.array([m_amp, m_phi, m_gamma]))\n",
    "    #m=np.array([[0.8,0.,0.6]])\n",
    "    mag_map=[m]\n",
    "    \n",
    "    initime=time.time()\n",
    "    print('Magnetization parameters initialized. Time spent:' , str(initime-starttime) , 's')\n",
    "    \n",
    "    for i, Te in enumerate(tes):\n",
    "        dm1=mag_incr(materials, sample, m_amp, m_phi, m_gamma, mat_gr_ind, mat_gr_ind_flat, mat_loc, Ms_sam, ex_stiff_sam, S_sam, Tc_sam, J_sam, lamda_sam, muat_sam, K0_sam, kappa_ani_sam, ani_perp_sam, alpha_par_sam, qs_sam, dbrillouin_t1_sam, dbrillouin_t2_sam, chi_par_num_sam, chi_par_denomm1_sam, m, Te, H_ext, gamma, dt)\n",
    "        m_heun=m+dm1\n",
    "        dm2=mag_incr(materials, sample, m_amp, m_phi, m_gamma, mat_gr_ind, mat_gr_ind_flat, mat_loc, Ms_sam, ex_stiff_sam, S_sam, Tc_sam, J_sam, lamda_sam, muat_sam, K0_sam, kappa_ani_sam, ani_perp_sam, alpha_par_sam, qs_sam, dbrillouin_t1_sam, dbrillouin_t2_sam, chi_par_num_sam, chi_par_denomm1_sam, m_heun, Te, H_ext, gamma, dt)\n",
    "        newmag=m+np.divide((dm1+dm2),2)\n",
    "        m=newmag\n",
    "        mag_map.append(m)\n",
    "    endtime=time.time()\n",
    "    print('Magnetization map created. Total time spent:' , str(endtime-starttime) , 's')\n",
    "    return np.array(mag_map)"
   ]
  }