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

# # Hysteresis simulations
# 
# In this tutorial, we summarise some of the convenience methods offered by Ubermag that can be used for simulating hysteresis. Let us first define a simple system object. For details on how to define simulations, please refer to other tutorials.

# In[1]:


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

region = df.Region(p1=(-50e-9, -50e-9, -50e-9), p2=(50e-9, 50e-9, 50e-9))
mesh = df.Mesh(region=region, cell=(5e-9, 5e-9, 5e-9))

system = mm.System(name="hysteresis")
system.energy = (
    mm.Exchange(A=1e-12)
    + mm.UniaxialAnisotropy(K=4e5, u=(0, 0, 1))
    + mm.DMI(D=1e-3, crystalclass="T")
)  # + mm.Demag()


def Ms_fun(point):
    x, y, z = point
    if x**2 + y**2 + z**2 <= 50e-9**2:
        return 1e6
    else:
        return 0


system.m = df.Field(mesh, nvdim=3, value=(0, 0, -1), norm=Ms_fun, valid="norm")


# Now that we have `system` object we can simulate hysteresis using `HysteresisDriver`. Like all other drivers, `HysteresisDriver` has `drive` method. As input arguments it takes `system` object (as usual) and:
# 
# - `Hmin` - the starting value of magnetic field
# - `Hmin` - the end value of magnetic field
# - `n` - the number of points between `Hmin` and `Hmax`
# 
# For instance, let us say we want to simulate hysteresis between $-1\,\text{T}$ and $1\,\text{T}$ applied along the $z$-direction in steps of $0.2\,\text{T}$. Accordingly, `Hmin` and `Hmax` are: 

# In[2]:


Hmin = (0, 0, -1 / mm.consts.mu0)
Hmax = (0, 0, 1 / mm.consts.mu0)


# The number of steps `n` should be 21, so that the values of magnetic field are: $-1\,\text{T}$, $-0.9\,\text{T}$, $-0.8\,\text{T}$, ..., $0.9\,\text{T}$, $1\,\text{T}$, $0.9\,\text{T}$, $0.8\,\text{T}$, ..., $-0.9\,\text{T}$, $-1\,\text{T}$.

# In[3]:


n = 21


# We can now create `HysteresisDriver` and drive the system.

# In[4]:


hd = mc.HysteresisDriver()
hd.drive(system, Hmin=Hmin, Hmax=Hmax, n=n)


# After the simulation is complete, we can have a look at the last magnetisation step:

# In[5]:


system.m.sel("y").mpl()


# Similarly, table can be viewed:

# In[6]:


system.table.data.head()


# From the table, we can see at what external magnetic field the system was relaxed:

# In[7]:


system.table.data["B_hysteresis"]


# The units of $B$ are:

# In[8]:


system.table.units["B_hysteresis"]


# ## Plotting hysteresis loop

# Using `Table` object, we can plot different values as a function of external magnetic field. For that, as usual, we use `mpl` method. We have to specify to that method what do we want to have on the $y$-axis.

# In[9]:


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


# This does not look like a hysteresis loop as we expected. The reason is that on the horizontal axis we have the magnitude `B` by default, which is always positive. We can change that by passing `Bz` for `x`:

# In[10]:


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


# `Table.mpl` is based on `matplotlib.pyplot.plot`, so any keyword argument accepted by it can be passed:

# In[11]:


system.table.mpl(
    x="Bz_hysteresis", y=["mz"], marker="o", linewidth=2, linestyle="dashed"
)


# ## `micromagneticdata` analysis
# 
# We can also analyse magnetisation fields at different values of external magnetic field as well as building interactive plots using `micromagneticdata` package.

# In[12]:


import micromagneticdata as md


# We can create a data object by passing system's name.

# In[13]:


data = md.Data(name=system.name)


# Now, we can have a look at all drives we did so far:

# In[14]:


data.info


# There is only one drive with index `0`. We can get if by indexing the data object:

# In[15]:


drive = data[0]


# The number of steps saved is:

# In[16]:


drive.n


# We can get an individual step by passing a value between 0 and 40 as an index:

# In[17]:


drive[0].sel("x").mpl()


# We can also interactively plot all individual steps by using `drive.hv`. For more details on creating interactive plots, please refer to other tutorials.

# In[18]:


drive.hv(kdims=["x", "y"])


# ## Multi-step Hysteresis Drivers
# Ubermag also provides a flexible version of a hysteresis driver which enables a hysteresis loop to be broken down into multiple curves.
# A list of hysteresis loop segments can be provided as an argument `Hsteps`. Each element of the list has the form `(Hstart, Hend, n)`.
# 
# To see the power of this technique, we will first initialise the field and then perform a full hysteresis loop including the virgin curve.

# In[19]:


system.m = df.Field(
    mesh,
    nvdim=3,
    value=lambda x: (0, 0, -1) if x[0] > 0 else (0, 0, 1),
    norm=Ms_fun,
    valid="norm",
)


# In[20]:


hd = mc.HysteresisDriver()
hd.drive(system, Hsteps=[[(0, 0, 0), Hmax, 10], [Hmax, Hmin, 10], [Hmin, Hmax, 20]])


# Finally we can plot the full curve.

# In[21]:


system.table.mpl(
    x="Bz_hysteresis", y=["mz"], marker="o", linewidth=2, linestyle="dashed"
)