#!/usr/bin/env python
# coding: utf-8

# # Dynamics equation
# 
# The dynamics of magnetisation field $\mathbf{m}$, without external excitations (e.g. spin-polarised current) is governed by the Landau-Lifshitz-Gilbert (LLG) equation
# 
# $$\frac{d\mathbf{m}}{dt} = \underbrace{-\gamma_{0}(\mathbf{m} \times \mathbf{H}_\text{eff})}_\text{precession} + \underbrace{\alpha\left(\mathbf{m} \times \frac{d\mathbf{m}}{dt}\right)}_\text{damping},$$
# 
# where $\gamma_{0} = \mu_{0}\gamma$ is the gyromagnetic ratio, $\alpha$ is the Gilbert damping, and $\mathbf{H}_\text{eff} = -\frac{1}{\mu_{0}M_\text{s}}\frac{\delta w(\mathbf{m})}{\delta \mathbf{m}}$ is the effective field. It consists of two terms: precession and damping. In this tutorial, we will explore some basic properties of this equation to understand how to define it in simulations.
# 
# ## Macrospin
# 
# We will study the simplest "zero-dimensional" case - macrospin. In the first step, after we import necessary modules and create the mesh which consists of a single discretisation cell.

# In[1]:


import oommfc as mc
import discretisedfield as df
import micromagneticmodel as mm

# Define a macrospin mesh (i.e. one discretisation cell).
p1 = (0, 0, 0)  # first point of the mesh domain (m)
p2 = (1e-9, 1e-9, 1e-9)  # second point of the mesh domain (m)
n = (1, 1, 1)  # discretisation cell size (m)

region = df.Region(p1=p1, p2=p2)
mesh = df.Mesh(region=region, n=n)


# Now, we can create a micromagnetic system object.

# In[2]:


system = mm.System(name="macrospin")


# Let us assume we have a simple Hamiltonian which consists of only Zeeman energy term
# 
# $$w = -\mu_{0}M_\text{s}\mathbf{m}\cdot\mathbf{H},$$
# 
# where $M_\text{s}$ is the saturation magnetisation, $\mu_{0}$ is the magnetic constant, and $\mathbf{H}$ is the external magnetic field. We apply the external magnetic field with magnitude $H = 2 \times 10^{6} \,\text{A}\,\text{m}^{-1}$ in the positive $z$ direction.

# In[3]:


H = (0, 0, 2e6)  # external magnetic field (A/m)
system.energy = mm.Zeeman(H=H)


# ## Dynamics simulation
# 
# In the next step we can define the system's dynamics. Let us assume we have $\gamma_{0} = 2.211 \times 10^{5} \,\text{m}\,\text{A}^{-1}\,\text{s}^{-1}$ and $\alpha=0.1$.

# In[4]:


gamma0 = 2.211e5  # gyromagnetic ratio (m/As)
alpha = 0.1  # Gilbert damping

system.dynamics = mm.Precession(gamma0=gamma0) + mm.Damping(alpha=alpha)


# To check what is our dynamics equation:

# In[5]:


system.dynamics


# Before we start running time evolution simulations, we need to initialise the magnetisation. In this case, our magnetisation is pointing in the positive $x$ direction with $M_\text{s} = 8 \times 10^{6} \,\text{A}\,\text{m}^{-1}$. The magnetisation is defined using `Field` class from the `discretisedfield` package we imported earlier.

# In[6]:


initial_m = (1, 0, 0)  # vector in x direction
Ms = 8e6  # magnetisation saturation (A/m)

system.m = df.Field(mesh, nvdim=3, value=initial_m, norm=Ms)


# Now, we can run the time evolution using `TimeDriver` for $t=0.1 \,\text{ns}$ and save the magnetisation configuration in $n=200$ steps.

# In[7]:


td = mc.TimeDriver()
td.drive(system, t=0.1e-9, n=200)


# ## Simulation results
# 
# How different system parameters vary with time, we can inspect by showing the system's datatable.

# In[8]:


system.table.data


# However, in our case it is much more informative if we plot the time evolution of magnetisation $z$ component $m_{z}(t)$.

# In[9]:


system.table.mpl(y=["mz"])


# Similarly, we can plot all three magnetisation components

# In[10]:


system.table.mpl(y=["mx", "my", "mz"])


# We can see that after some time the macrospin aligns parallel to the external magnetic field in the $z$ direction. We can also have a look at the dynamics using an interactive plot.

# ### Exercise 1
# 
# Modify Gilbert damping and set it to $\alpha=0.005$ and see what happens with the dynamics.
# 
# **Solution**

# In[11]:


system.dynamics.damping.alpha = 0.005
system.m = df.Field(mesh, nvdim=3, value=initial_m, norm=Ms)

td.drive(system, t=0.1e-9, n=200)

system.table.mpl(y=["mx", "my", "mz"])


# ### Exercise 2
# 
# Repeat the simulation with $\alpha=0.1$ and $\mathbf{H} = (0, 0, -2\times 10^{6})\,\text{A}\,\text{m}^{-1}$.
# 
# **Solution**

# In[12]:


system.energy.zeeman.H = (0, 0, -2e6)
system.dynamics.damping.alpha = 0.1
system.m = df.Field(mesh, nvdim=3, value=initial_m, norm=Ms)

td.drive(system, t=0.1e-9, n=200)

system.table.mpl(y=["mx", "my", "mz"])


# ### Exercise 3
# 
# Keep using $\alpha=0.1$. Change the field from `H = (0, 0, -2e6)` to `H = (0, -1.41e6, -1.41e6)`, and plot
# $m_x(t)$, $m_y(t)$ and $m_z(t)$ as above. Can you explain the (initially non-intuitive) output?

# In[13]:


system.energy.zeeman.H = (0, -1.41e6, -1.41e6)

td.drive(system, t=0.1e-9, n=200)

system.table.mpl(y=["mx", "my", "mz"])