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

# # FMR standard problem
# 
# ## Problem specification
# 
# We choose a cuboidal thin film permalloy sample measuring $120 \times 120 \times 10 \,\text{nm}^{3}$. The choice of a cuboid is important as it ensures that the finite difference method employed by OOMMF does not introduce errors due to irregular boundaries that cannot be discretized well. We choose the thin film geometry to be thin enough so that the variation of magnetization dynamics along the out-of-film direction can be neglected. Material parameters based on permalloy are:
# 
# - exchange energy constant $A = 1.3 \times 10^{-11} \,\text{J/m}$,
# - magnetisation saturation $M_\text{s} = 8 \times 10^{5} \,\text{A/m}$,
# - Gilbert damping $\alpha = 0.008$.
# 
# An external magnetic bias field with magnitude $80 \,\text{kA/m}$ is applied along the direction $e = (1, 0.715, 0)$. We choose the external magnetic field direction slightly off the sample diagonal in order to break the system’s symmetry and thus avoid degenerate eigenmodes. First, we initialize the system with a uniform out-of-plane magnetization $m_{0} = (0, 0, 1)$. The system is allowed to relax for $5 \,\text{ns}$, which was found to be sufficient time to obtain a well-converged equilibrium magnetization configuration. We refer to this stage of simulation as the relaxation stage, and its final relaxed magnetization configuration is saved to serve as the initial configuration for the next dynamic stage. Because we want to use a well defined method that is supported by all simulation tools, we minimize the system’s energy by integrating the LLG equation with a large, quasistatic Gilbert damping $\alpha = 1$ for $5 \,\text{ns}$. In the next step (dynamic stage), a simulation is started using the equilibrium magnetisation configuration from the relaxation stage as the initial configuration. Now, the direction of an external magnetic field is altered to $e = (1, 0.7, 0)$. This simulation stage runs for $T = 20 \,\text{ns}$ while the (average and spatially resolved) magnetization $M(t)$ is recorded every $\Delta t = 5 \,\text{ps}$. The Gilbert damping in this dynamic simulation stage is $\alpha = 0.008$.
# 
# Details of this standard problem specification can be found in Ref. 1.
# 
# ## Relaxation stage
# 
# Firstly, all required modules are imported.

# In[1]:


import oommfc as oc
import discretisedfield as df


# Now, we specify all simulation parameters.

# In[2]:


import numpy as np

lx = ly = 120e-9  # x and y dimensions of the sample(m)
lz = 10e-9  # sample thickness (m)
dx = dy = dz = 5e-9  # discretisation in x, y, and z directions (m)

Ms = 8e5  # saturation magnetisation (A/m)
A = 1.3e-11  # exchange energy constant (J/m)
H = 8e4 * np.array([0.81345856316858023, 0.58162287266553481, 0.0])
alpha = 0.008  # Gilbert damping
gamma = 2.211e5


# Now, the system object can be created and mesh, magnetisation, hamiltonian, and dynamics are specified.

# In[3]:


mesh = oc.Mesh(p1=(0, 0, 0), p2=(lx, ly, lz), cell=(dx, dy, dz))

system = oc.System(name="stdprobfmr")

system.hamiltonian = oc.Exchange(A) + oc.Demag() + oc.Zeeman(H)
system.dynamics = oc.Precession(gamma) + oc.Damping(alpha)
system.m = df.Field(mesh, value=(0, 0, 1), norm=Ms)


# Finally, the system is relaxed.

# In[4]:


md = oc.MinDriver()
md.drive(system)


# We can now load the relaxed state to the Field object and plot the $z$ slice of magnetisation.

# In[5]:


get_ipython().run_line_magic('matplotlib', 'inline')
system.m.plot_plane("z")


# ## Dynamic stage

# In the dynamic stage, we use the relaxed state from the relaxation stage.

# In[6]:


# Change external magnetic field.
H = 8e4 * np.array([0.81923192051904048, 0.57346234436332832, 0.0])
system.hamiltonian.zeeman.H = H


# Finally, we run the multiple stage simulation using `TimeDriver`.

# In[7]:


T = 20e-9
n = 4000

td = oc.TimeDriver()
td.drive(system, t=T, n=n)


# ## Postprocessing

# From the obtained vector field samples, we can compute the average of magnetisation $y$ component and plot its time evolution. 

# In[8]:


import matplotlib.pyplot as plt

t = system.dt['t'].as_matrix()
my = system.dt['mx'].as_matrix()

# Plot <my> time evolution.
plt.figure(figsize=(8, 6))
plt.plot(t, my)
plt.xlabel('t (ns)')
plt.ylabel('my average')
plt.grid()


# From the $<m_{y}>$ time evolution, we can compute and plot its Fourier transform.

# In[10]:


import scipy.fftpack

psd = np.log10(np.abs(scipy.fftpack.fft(my))**2)
f_axis = scipy.fftpack.fftfreq(4000, d=20e-9/4000)

plt.plot(f_axis/1e9, psd)
plt.xlim([6, 12])
plt.xlabel('f (GHz)')
plt.ylabel('Psa (a.u.)')
plt.grid()


# ## References
# [1] A. Baker, M. Beg, G. Ashton, M. Albert, D. Chernyshenko, W. Wang, S. Zhang, M.-A. Bisotti, M. Franchin, C.L. Hu, R. Stamps, T. Hesjedal, and H. Fangohr, *J. Magn. Magn. Mater.* **421**, 428 (2017).