<jupyter_start><jupyter_text>**Calculating Stoner-Wohlfarth Curves** The Stoner-Wohlfarth model is given in the paper "A Mechanism of Magnetic Hysteresis in Heterogeneous Alloys":> [<NAME> and <NAME>, Phil. Trans. Royal Soc. London 240 (1948) 599](https://doi.org/10.1098/rsta.1948.0007)- The SW-model considers magnetisation processes governed by **rotations**  of the magnetisation vector.- Magnetisation processes based on domains and domain wall movement are **NOT** considered.The model is therefore strictly valid for:- single-domain particles ($d < d_{c}$) or - ensembles of non-interacting, single-domain particlesBut we can also apply the SW-model in cases where *domain wall motion is not relevant*:- hard axis magnetisation curve of a single crystal with uniaxial anisotropy- magnetisation curve of a polycrystal at applied field strengths above those where domain wall movement occurs (neglecting particle-particle interactions). Consider a single-domain particle with uniaxial magnetocrystalline anisotropy.- In zero applied field, $H$, the magnetisation, $M_{0}$, will be parallel to the easy axis.- If we apply a field, $H$, perpendicular to the easy axis, $M_{0}$ will be forced to rotate away from the easy axis, towards $H$.  The angle between the field and the magnetisation, $\theta_{M}$, results from miminising the sum of the **anisotropy energy** and the **Zeeman (field) energy**: $E_{v} = K_{1}sin^{2}(\theta_{M}-\theta_{H})  + K_{2}sin^{4}(\theta_{M}-\theta_{H}) - M_{0}Hcos{\theta_{M}}$ where:$E_{v}$ = the total energy per unit volume$M_{0}$ = the spontaneous magnetisation of the material$K_{1}$, $K_{2}$ = second and fourth order magnetocrystalline anisotropy constants$\theta_{H}$ = angle between the applied field and the easy axis$\theta_{M}$ = angle between the applied field and the magnetisation$H$ = applied field strength Before we can write down the equation in code, we need to import some things we will need:<jupyter_code>import numpy
import matplotlib.pyplot as plt
from matplotlib.ticker import StrMethodFormatter
from scipy.optimize import minimize
import ipywidgets as widgets
from ipywidgets import interactive, fixed

%matplotlib inline
%config InlineBackend.figure_format = 'svg'

plt.rc('font', size=10)
plt.rc('figure', figsize=(6,4))
def get_energy_terms(thetaM, thetaH, H, M0, K1, K2):
    """
    Calculate the anisotropy energy, the Zeeman energy and the total energy, given
    the arguments.
    
    Args:
        thetaM = angle in degrees between the magnetisation and the applied field
        thetaH = angle in degrees between the easy axis and the applied field
        H = applied field in Oe
        M0 = spontaneous magnetisation in erg/Oe.cm^3 (== kA/m)
        K1 = second order anisotropy constant in erg/cm^3
        K2 = fourth order anisotropy constant in erg/cm^3
    
    Returns: 
        AE = float, the anisotropy energy per unit volume for the conditions given, unit: erg/cm^3.
        ZE = float, the Zeeman energy per unit volume for the conditions given, unit: erg/cm^3.
        TE = float, the total energy per unit volume for the conditions given, unit: erg/cm^3.
    """
    t1 = K1*numpy.sin(numpy.deg2rad(thetaM-thetaH))**2
    t2 = K2*numpy.sin(numpy.deg2rad(thetaM-thetaH))**4

    AE = t1+t2
    ZE = M0*H*numpy.cos(numpy.deg2rad(thetaM))
    TE = AE-ZE
    
    return AE, ZE, TE<jupyter_output><empty_output><jupyter_text>Example with $Y_{2}Fe_{14}B$This material has a tetragonal crystal structure and uniaxial magnetocrystalline anisotropy, where the  axis is the easy axis.Here are values for for the constants for $Y_{2}Fe_{14}B$ extracted from [Givord *et al.* Solid State Comm. 51 (1984) 857-860](https://doi.org/10.1016/0038-1098(84)91087-1).<jupyter_code>M0 = 1053 # erg/Oecm^3 == kA/m (~ 1.4 T )
K1 = 1.25e7 # erg/cm^3 
K2 = 0.05e7 # erg/cm^3<jupyter_output><empty_output><jupyter_text>Let's examine how the energy per unit volume, $E_{v}$, changes as a function of the angle between the field and the magnetisation, $\theta_{M}$ for different applied field strengths, $H$, using an interactive plot.<jupyter_code>def plot_func(thetaH, H, M0, K1, K2):
    thetaM = numpy.arange(91)
    AE, ZE, TE = get_energy_terms(thetaM, thetaH, H, M0, K1, K2)
    
    fig, ax = plt.subplots()
    
    ax.plot(thetaM, AE, 'b-', lw=2.0, label='Anisotropy Energy')
    ax.plot(thetaM, ZE, 'r-', lw=2.0, label='Zeeman Energy')
    ax.plot(thetaM, TE, 'g-', lw=2.0, label='Total Energy')

    ax.grid(which="major")
    ax.set_xlabel(r"$\theta_{M} \ /\ \degree$")
    ax.set_ylabel(r"$E_{v}\ /\ erg \cdot cm^{-3}$")
    ax.legend(loc="upper right")
    ax.set_xlim(xmin=0.0, xmax=90.0)
    ax.yaxis.set_major_formatter(StrMethodFormatter('{x: ,.2e}'))


    
interactive_plot = interactive(plot_func, thetaH=[90,0], 
                               H=widgets.IntSlider(value=0.0, min=0.0, max=28000, 
                                                   step=1000, description="H / Oe"),
                               M0=fixed(M0), K1=fixed(K1), K2=fixed(K2))

output = interactive_plot.children[-1]
output.layout.height = '360px'
interactive_plot<jupyter_output><empty_output><jupyter_text>Let's redefine our fuction so that it only returns the total energy (we don't need the other terms individually now):<jupyter_code>def SW_E(thetaM, thetaH, H, M0, K1, K2):
    """
    Calculate the energy per unit volume resulting from the
    balance of the anisotropy energy and Zeeman energy, given
    the arguments.
    
    Args:
        thetaM = angle in degrees between the magnetisation and the applied field
        thetaH = angle in degrees between the easy axis and the applied field
        H = applied field in Oe
        M0 = spontaneous magnetisation in erg/Oe.cm^3 (== kA/m)
        K1 = second order anisotropy constant in erg/cm^3
        K2 = fourth order anisotropy constant in erg/cm^3
    
    Returns: 
        E = float, the energy per unit volume for the conditions given, unit: erg/cm^3.
    """
    t1 = K1*numpy.sin(numpy.deg2rad(thetaM-thetaH))**2
    t2 = K2*numpy.sin(numpy.deg2rad(thetaM-thetaH))**4
    t3 = M0*H*numpy.cos(numpy.deg2rad(thetaM))
    return t1+t2-t3<jupyter_output><empty_output><jupyter_text>We can use `scipy.optimize.minimize` to find the minimum of $E_{v}$ for a particular applied field and angle $\theta_{H}$ like this:<jupyter_code>q = minimize(SW_E, 90, args=(90, 8000, M0, K1, K2), method="Nelder-Mead")
q
q.x[0]<jupyter_output><empty_output><jupyter_text>We can use `scipy.optimize.minimize` in a **loop** to find the minimum of $E_{v}$ for a several applied field values like this:<jupyter_code>def getData_single(FR, M0, K1, K2, thetaH):
    """
    For each value of the field in the field range (FR), find the minimum
    of the absolute value of the derivative of the energy equation with
    respect to the angle theta.
    
    Args:
        FR : array of applied field values in Oe
        thetaH = angle in degrees between the easy axis and the applied field
        M0 = spontaneous magnetisation in erg/Oe.cm^3 (== kA/m)
        K1 = second order anisotropy constant in erg/cm^3
        K2 = fourth order anisotropy constant in erg/cm^3
        
    Returns:
        D : a list of angles thetaM in degrees at which the derivative of the energy
            equation is zero for the field values defined in FR
    """
    D = []
    IG = thetaH
    for H in FR:
        q = minimize(SW_E, IG, args=(thetaH, H, M0, K1, K2), method="Nelder-Mead")
        D.append(q.x[0])
        IG = q.x[0]

        
    return D
# field range (start, stop, step) in Oe
FR = numpy.arange(0, 6.6e4, 1.0e3)
D = getData_single(FR, M0, K1, K2, 90.0)<jupyter_output><empty_output><jupyter_text>Now,  $\cos\theta_{M} = \frac{M}{M_{0}}$, so to plot M(H) we need to multiply the cosine of the angles, $\theta_{M}$ in D by $M_{0}$<jupyter_code>plt.plot(FR, numpy.cos(numpy.deg2rad(D))*M0, 'r-')
plt.xlabel("H / Oe")
plt.ylabel("M / kA/m")
plt.xlim(xmin=0)
plt.ylim(ymin=0)
plt.axvline(x=25800, ls=':')
plt.annotate(r"$\mathregular{H = \frac{2(K_{1} + 2K_{2})}{M_{0}}}}$", (25800, 1053), (32000, 800), 
             arrowprops={"arrowstyle" : "-", "relpos" : (0,1), "ls" : ':'}, fontsize=14)<jupyter_output><empty_output><jupyter_text>Experimental results for a $Y_{2}Fe_{14}B$  single crystal were measured by [Givord *et al.* Solid State Comm. 51 (1984) 857-860](https://doi.org/10.1016/0038-1098(84)91087-1).Two text files containing experimental data points extracted from the above paper are in the folder "resources" and can be loaded using  `numpy.loadtxt`.<jupyter_code># load Givord's measurements of a Y2Fe14B SX at 300K

Ge = numpy.loadtxt("resources/Givord_Y2Fe14B_001_conv.txt")
Gh = numpy.loadtxt("resources/Givord_Y2Fe14B_100_and_110_conv.txt")
plt.plot(Gh[:,0],Gh[:,1],'ko', label="Experiment, [100] or [110]")
plt.plot(FR, numpy.cos(numpy.deg2rad(D))*M0, 'r-', label="SW curve, hard axis")
plt.xlabel("H / Oe")
plt.ylabel("M / kA/m")
plt.xlim(xmin=0)
plt.ylim(ymin=0)
plt.legend()<jupyter_output><empty_output><jupyter_text>We have very good agreement between the experimentally measured data for a single crystal and our calculated Stoner-Wohlfarth curve. Now we can use the same code to calculate an easy axis curve to compare with the experimental data:<jupyter_code>D1 = getData_single(FR, M0, K1, K2, 0.0)
plt.plot(Gh[:,0],Gh[:,1],'ko', label="Experiment, [100] or [110]")
plt.plot(FR, numpy.cos(numpy.deg2rad(D))*M0, 'r-', label="SW curve, hard axis")
plt.plot(Ge[:,0],Ge[:,1],'go', label="Experiment, [001]")
plt.plot(FR, numpy.cos(numpy.deg2rad(D1))*M0, 'b-', label="SW curve, easy axis")

plt.xlabel("H / Oe")
plt.ylabel("M / kA/m")
plt.xlim(xmin=0)
plt.ylim(ymin=0)
plt.legend()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1277054056826341, 'fraction_numerical': 0.03743973740896502, 'mean_word_length': 3.8896690070210633, 'pattern_counts': {'":': 1, '<': 40, '<?xml version=': 0, '>': 40, 'https://': 3, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '9', 'dir': 'jupyter-structured-clean-dedup', 'id': '374330', 'license': 'MIT', 'n_tokens_mistral': 3624, 'n_tokens_neox': 3301, 'n_words': 1045, 'path': 'SW_curves_1_single.ipynb', 'repo_name': 'tgwoodcock/Stoner-Wohlfarth-Tutorial', 'size': '9762'}