{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ca606099",
   "metadata": {
    "tags": [
     "remove-input"
    ]
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "62388ef2",
   "metadata": {},
   "source": [
    "# Dirac Strings in Kagome ASI\n",
    "Here we demonstrate how flatspin was to produce Dirac Strings following the experimental setup given by {cite}`mengotti2011`.\n",
    "\n",
    "[“Real-space observation of emergent magnetic monopoles and associated Dirac strings in artificial kagome spin ice” Nature Phys 7, 68–74 (2011)](https://doi.org/10.1038/nphys1794)\n",
    "\n",
    "The flatspin results are discussed in more detail in our paper {cite}`Flatspin2022`."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5ad2a3f5",
   "metadata": {},
   "source": [
    "## Creating the Dataset\n",
    "We set our parameters to closely match the experimental setup in {cite}`mengotti2011`.\n",
    "\n",
    "Below is the `flatspin-run` command used to generate the dataset for this section, followed by an explaination of the parameters used.\n",
    "\n",
    "```flatspin-run -m KagomeSpinIceRotated -e Triangle -p phase=0 -p phi=-3.6 -p \"size=(29,29)\" -p temperature=300 -p H=0.1 -p periods=1  -p sw_beta=2.5 -p sw_b=0.212 -p sw_gamma=3 -p sw_c=1 -p hc=0.216 -p alpha=0.00103 -p use_opencl=1 -p neighbor_distance=10 -p disorder=0.05  -p \"m_therm=0.05*1.29344e-15\" -p timesteps=2000 -p spp=2000 -o diracStringsPublish2```\n",
    "\n",
    "`-m KagomeSpinIceRotated` and `-p size=(29,29)` define the geometry and the number of magnets similar to the experimental setup.\n",
    "\n",
    "`-e Triangle`, `-p phi=-3.6` and `-p phase=0` uses a Triangle encoder to set up a reversal field at an angle of -3.6 degrees.  As `H=0.1`, the field starts at -0.1 T ramps up linearly to 0.1 T then back down to -0.1 T\n",
    "\n",
    "A `temperature` of 300 K is used to simulate room temperature.\n",
    "\n",
    "Using a magnetization saturation ($M_S$) of 860 kA/m, and the volume of the magnets given in the experimental setup (~1.5e-21 m^3), the `m_therm` parameter is taken to be 5% of volume * msat.\n",
    "\n",
    "The value of `alpha` is calculated from $\\alpha = \\frac{\\mu_0 M}{4\\pi a^3}$ (see [](theory)), with $M$ = 860e3 * 1.504e-21 Am^2 and $a$ = 500 nm\n",
    "\n",
    "Micromagnetic simulations of magnets with this msat and the given dimensions (470 nm * 160 nm * 20 nm) were used to obtain the switching parameters `sw_beta=2.5`, `sw_b=0.212`, `sw_gamma=3`, `sw_c=1` and `hc=0.216`.\n",
    "\n",
    "5% disorder was used to account for variations in the magnets used in the experiments."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "24e03245",
   "metadata": {},
   "source": [
    "## Calculate $H_c$\n",
    "First we analyze the dataset created by the above `flatspin-run` command to calculate the switching field of the full lattice, $H_c$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0df5ab79",
   "metadata": {},
   "outputs": [],
   "source": [
    "from flatspin import data\n",
    "\n",
    "ds = data.Dataset.read(\"/data/flatspin/diracStringsPublish2\")\n",
    "\n",
    "# for now we're only interested in the fist half of the run\n",
    "t = slice(None, 1000)\n",
    "\n",
    "# using grid_size=(1,1) returns the average magnetization over the whole lattice\n",
    "mag = data.load_output(ds, \"mag\", grid_size=(1, 1), t=t)\n",
    "\n",
    "# get timestep where array switches (the time where the absolute magnetization in the x-direction is minimized)\n",
    "t_min = np.argmin(abs(mag[:, 0]))\n",
    "print(f\"Timestep where array switches, t_min = {t_min}\")\n",
    "\n",
    "# load the data for the external field and use t_min to get find HC\n",
    "h_ext = data.load_output(ds, \"h_ext\", t=t)\n",
    "Hc = h_ext[t_min][0]\n",
    "print(f\"H_c = {Hc}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5d6e4421",
   "metadata": {},
   "source": [
    "## Find the field values of interest\n",
    "In {cite}`mengotti2011` they show the state of ASI at field values: `[0.8HC, 0.85HC, 0.92HC, 0.95HC, 0.99HC, 1.06HC]`.\n",
    "To allow us to compare the results of flatspin to the experimental setup, we will find the timesteps in our dataset where the field is closest to these values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3c784b23",
   "metadata": {},
   "outputs": [],
   "source": [
    "# calculate the fields of interest in terms of our HC\n",
    "foi = [0.8, 0.85, 0.92, 0.99, 1.06]\n",
    "foiHC = Hc * np.array(foi)\n",
    "print(f\"foiHC = {foiHC}\")\n",
    "\n",
    "# find the nearest times by minimizing the absolute difference between the field and the HC\n",
    "nearest_time = [np.argmin(np.abs(field - h_ext[:, 0])) for field in foiHC]\n",
    "print(f\"nearest_time = {nearest_time}\")\n",
    "print(f\"nearest fields = {[str(round(h_ext[t, 0]/Hc,2))+'HC' for t in nearest_time]}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "84d330c0",
   "metadata": {},
   "source": [
    "Below we animate the state of the ASI, as it evolves from the first to the last field value of interest (`0.8HC` to `1.06HC`)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "56b0a4fc",
   "metadata": {},
   "outputs": [],
   "source": [
    "from IPython.display import HTML\n",
    "from matplotlib.animation import FuncAnimation\n",
    "from flatspin.plotting import plot_vectors\n",
    "\n",
    "def animate_dirac_strings(ds, times):\n",
    "    fig, ax = plt.subplots(figsize=(10, 10), facecolor=(0.4, 0.4, 0.4))\n",
    "    fig.subplots_adjust(left=0, right=1, bottom=0, top=0.95, wspace=0, hspace=0)\n",
    "    ax.set_axis_off()\n",
    "\n",
    "    _, UV = data.read_vectors(ds.tablefiles(), \"mag\", times)\n",
    "    positions, _ = data.read_geometry(ds.tablefile('geometry'))\n",
    "\n",
    "    def animate(i):\n",
    "        plot_vectors(positions, UV[i], arrows=True, replace=True, ax=ax, cmap=\"peem180\")\n",
    "        ax.set_title(f\"{round(h_ext[times[i],0]/Hc,2)}$H_c$\", fontsize=20, color=\"white\")\n",
    "\n",
    "\n",
    "    anim = FuncAnimation(\n",
    "        fig, animate, init_func=lambda: None,\n",
    "        frames=len(times), interval=100, blit=False\n",
    "    )\n",
    "    plt.close() # Only show the animation\n",
    "    return HTML(anim.to_jshtml(fps=8))\n",
    "#animate_dirac_strings(ds, times=nearest_time)\n",
    "#animate_dirac_strings(ds, times=list(range(710,750,1)))\n",
    "animate_dirac_strings(ds, times=list(range(nearest_time[0], nearest_time[-1]+1)))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fdb9d1c2",
   "metadata": {},
   "source": [
    "## Hysteresis\n",
    "Now we plot the hysterisis of our dataset, as well as a sketch of the hysteresis shown in {cite}`mengotti2011`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0333f8bf",
   "metadata": {},
   "outputs": [],
   "source": [
    "mag = data.load_output(ds, \"mag\", grid_size=(1, 1))  # now we want full run length\n",
    "h_ext = data.load_output(ds, \"h_ext\")\n",
    "plt.figure(figsize=(8, 5))\n",
    "plt.plot(h_ext[:, 0] / Hc, mag[:, 0], label=\"flatspin\")\n",
    "\n",
    "meng = np.loadtxt(\"scaledMeng.txt\").round(2)  # our rough estimate of the Mengotti et al. dataset inferred from their graph\n",
    "plt.plot(meng[:, 0], meng[:, 1], label=\"Mengotti\")\n",
    "\n",
    "#make text marker at the fields of interest\n",
    "rom_num = [\"I\", \"II\", \"III\", \"IV\", \"V\"]\n",
    "offsetx = [-.1,0,.1,-.3,-.3]\n",
    "offsety = [-.15,-.15,-.15,0,0]\n",
    "for i, field in enumerate(foi):\n",
    "    plt.plot(field, mag[nearest_time[i], 0], \".\", label=f\"{rom_num[i]}\", color=\"black\")\n",
    "    plt.text(field+offsetx[i], mag[nearest_time[i], 0] + offsety[i], rom_num[i], color=\"black\")\n",
    "    \n",
    "plt.xlabel(\"$H/H_c$\")\n",
    "plt.ylabel(\"$M/M_S$\")\n",
    "plt.ylim(-1.2, 1.2)\n",
    "plt.legend([\"flatspin\", \"Mengotti et al.\"]);"
   ]
  }
 ],
 "metadata": {
  "jupytext": {
   "formats": "ipynb,md:myst",
   "text_representation": {
    "extension": ".md",
    "format_name": "myst",
    "format_version": 0.13,
    "jupytext_version": "1.13.8"
   }
  },
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.10"
  },
  "source_map": [
   15,
   20,
   29,
   52,
   57,
   76,
   82,
   92,
   96,
   123,
   128
  ]
 },
 "nbformat": 4,
 "nbformat_minor": 5
}