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

# # Tutorial 05: Cubic anisotropy energy term
# 
# Cubic anisotropy energy density is computed as
# 
# $$w_{ca} = -K_{1} [(\mathbf{m} \cdot
#            \mathbf{u}_{1})^{2}(\mathbf{m} \cdot
#            \mathbf{u}_{2})^{2} + (\mathbf{m} \cdot
#            \mathbf{u}_{2})^{2}(\mathbf{m} \cdot
#            \mathbf{u}_{3})^{2} + (\mathbf{m} \cdot
#            \mathbf{u}_{1})^{2}(\mathbf{m} \cdot
#            \mathbf{u}_{3})^{2}]$$
# 
# where $\mathbf{m}$ is the normalised ($|\mathbf{m}|=1$) magnetisation, $K_{1}$ is the cubic anisotropy constant, and $\mathbf{u}_{1}$ and $\mathbf{u}_{2}$ are the anisotropy axes. Cubic anisotropy energy term tends to align all magnetic moments parallel or antiparallel to one of the three anisotropy axes.
# 
# In `oommfc`, $\mathbf{m}$ is a part of the magnetisation field `system.m`. Therefore, only cubic anisotropy constant $K_{1}$ and axes $\mathbf{u}_{1}$ and $\mathbf{u}_{2}$ should be provided as input parameters to uniquely define the energy term. All parameters can be constant in space or spatially varying.
# 
# ## Spatially constant $K_{1}$ and $\mathbf{u}$
# 
# Let us start by assembling a simple simple simulation where neither $K$ nor $\mathbf{u}$ vary in space. The sample is a "one-dimensional" chain of magnetic moments.

# In[1]:


import discretisedfield as df
import micromagneticmodel as mm

import oommfc as oc

p1 = (-10e-9, 0, 0)
p2 = (10e-9, 1e-9, 1e-9)
cell = (1e-9, 1e-9, 1e-9)

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


# The mesh is

# In[2]:


# NBVAL_IGNORE_OUTPUT
mesh.k3d()


# The system has a Hamiltonian, which consists of only cubic anisotropy energy term.

# In[3]:


K = 1e5  # cubic anisotropy constant (J/m**3)
u1 = (0, 0, 1)  # cubic anisotropy axis
u2 = (0, 1, 0)  # cubic anisotropy axis
system = mm.System(name='cubicanisotropy_constant_K_u')
system.energy = mm.CubicAnisotropy(K=K, u1=u1, u2=u2)


# We are going to minimise the system's energy using `oommfc.MinDriver` later. Therefore, we do not have to define the system's dynamics equation. Finally, we need to define the system's magnetisation (`system.m`). We are going to make it random with $M_\text{s}=8\times10^{5} \,\text{Am}^{-1}$

# In[4]:


import random

import discretisedfield as df

Ms = 8e5  # saturation magnetisation (A/m)

def m_fun(pos):
    return [2*random.random()-1 for i in range(3)]

system.m = df.Field(mesh, dim=3, value=m_fun, norm=Ms)


# The magnetisation, we set is

# In[5]:


# NBVAL_IGNORE_OUTPUT
system.m.k3d.vector(color_field=system.m.z)


# Now, we can minimise the system's energy by using `oommfc.MinDriver`.

# In[6]:


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


# We expect that now all magnetic moments are aligned parallel or antiparallel to the anisotropy axes.

# In[7]:


# NBVAL_IGNORE_OUTPUT
system.m.k3d.vector(color_field=system.m.z)


# ## Spatially varying $\mathbf{K}_{1}$
# 
# There are two different ways how a parameter can be made spatially varying, by using:
# 1. Dictionary
# 2. `discretisedfield.Field`
# 
# ### Dictionary
# 
# In order to define a parameter using a dictionary, regions must be defined in the mesh. Regions are defined as a dictionary, whose keys are the strings and values are `discretisedfield.Region` objects, which take two corner points of the region as input parameters. 

# In[8]:


p1 = (-10e-9, 0, 0)
p2 = (10e-9, 1e-9, 1e-9)
cell = (1e-9, 1e-9, 1e-9)
subregions = {'region1': df.Region(p1=(-10e-9, 0, 0), p2=(0, 1e-9, 1e-9)),
              'region2': df.Region(p1=(0, 0, 0), p2=(10e-9, 1e-9, 1e-9))}
region = df.Region(p1=p1, p2=p2)
mesh = df.Mesh(region=region, cell=cell, subregions=subregions)


# In[9]:


# NBVAL_IGNORE_OUTPUT
mesh.k3d_subregions()


# Let us say that there is no cubic anisotropy ($K=0$) in region 1, whereas in region 2 it is $K=10^{5} \,\text{Jm}^{-3}$. `K` is now defined as a dictionary:

# In[10]:


K = {'region1': 0, 'region2': 1e5}


# The system object is

# In[11]:


system = mm.System(name='cubicanisotropy_dict_K')
system.energy = mm.CubicAnisotropy(K=K, u1=u1, u2=u2)
system.m = df.Field(mesh, dim=3, value=m_fun, norm=Ms)


# Its magnetisation is

# In[12]:


# NBVAL_IGNORE_OUTPUT
system.m.k3d.vector(color_field=system.m.z)


# After we minimise the energy

# In[13]:


# NBVAL_IGNORE_OUTPUT
md.drive(system)


# The magnetisation is as we expected.

# In[14]:


# NBVAL_IGNORE_OUTPUT
system.m.k3d.vector(color_field=system.m.z)


# ### `discretisedfield.Field`
# 
# Let us define the spatailly varying uniaxial anisotropy, so that
# 
# $K(x, y, z) = \left\{
# \begin{array}{ll}
# 0 & x \le 0 \\
# 1e5 & x > 0 \\
# \end{array}
# \right. $
# 
# The value of `K` for the spatially varying anisotropy is set using a Python function.

# In[15]:


def K_fun(pos):
    x, y, z = pos
    if x <= 0:
        return 0
    else:
        return 1e5


# The uniaxial anisotropy parameters are

# In[16]:


K = df.Field(mesh, dim=1, value=K_fun)


# The system is

# In[17]:


system = mm.System(name='cubicanisotropy_field_u')
system.energy = mm.CubicAnisotropy(K=K, u1=u1, u2=u2)
system.m = df.Field(mesh, dim=3, value=m_fun, norm=Ms)


# and its magnetisation is

# In[18]:


# NBVAL_IGNORE_OUTPUT
system.m.k3d.vector(color_field=system.m.z)


# After the energy minimisation, the magnetisation is:

# In[19]:


# NBVAL_IGNORE_OUTPUT
md.drive(system)
system.m.k3d.vector(color_field=system.m.z)