test_equivalent_source_SIWave-to-HFSS.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "843a3b85-2c98-4350-be20-b14ee673e14b",
   "metadata": {},
   "outputs": [],
   "source": [
    "## Create equivalent current sources from target field profile\n",
    "\n",
    "import numpy as np\n",
    "import meep as mp\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.animation as animation\n",
    "from scipy.interpolate import griddata\n",
    "from IPython.display import Video"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7952207d-6979-484e-9eee-6576408d0b86",
   "metadata": {},
   "outputs": [],
   "source": [
    "# a = 1 mm\n",
    "\n",
    "wvl  = 300\n",
    "freq = 1 / wvl\n",
    "\n",
    "dpml       = 0.2 * wvl\n",
    "pml_layers = [mp.PML(dpml)]\n",
    "\n",
    "resolution = 0.2\n",
    "print(\"resolution               %.2f   pixels / mm\" % resolution)\n",
    "print(\"step size                %.2f   mm\" % (1/resolution))\n",
    "print(\"Pixels per wavelength    %i\" % (resolution*wvl))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6420875d-e6e2-48e2-a85f-1478aff10f79",
   "metadata": {},
   "outputs": [],
   "source": [
    "# cell size\n",
    "Sx = 1300\n",
    "Sy = 1300\n",
    "Sz = 300\n",
    "cell_size = mp.Vector3(Sx, Sy, Sz)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "86048055-aaa1-4d91-b09e-f5d2e5c94d6e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Function: get equivalent sources\n",
    "def get_equiv_sources(field:np.ndarray, src:mp.SourceTime, center:mp.Vector3, size:mp.Vector3):\n",
    "    \n",
    "    # get fields\n",
    "    Ex, Ey, Ez, Hx, Hy, Hz = field\n",
    "    \n",
    "    # Get normal vector nHat\n",
    "    if size[2]: # 3D sim\n",
    "        if size[0] and size[1]:\n",
    "            raise Exception(\"size should be a plane and thus have 2 non-zero values in 3D\")\n",
    "        elif size[0]:\n",
    "            nHat = mp.Vector3(0, 1, 0)\n",
    "        elif size[1]:\n",
    "            nHat = mp.Vector3(1, 0, 0)\n",
    "        else:\n",
    "            raise Exception(\"size should be a plane and thus have 2 non-zero values in 3D\")\n",
    "    elif size[0] and size[1]: # 3D sim ?\n",
    "        nHat = mp.Vector3(0, 0, 1)\n",
    "    else: # 2D\n",
    "        if size[0]:\n",
    "            nHat = mp.Vector3(0, 1)\n",
    "        elif size[1]:\n",
    "            nHat = mp.Vector3(1, 0)\n",
    "        else:\n",
    "            raise Exception(\"size should be a plane and thus have 1 non-zero value in 2D\")\n",
    "    \n",
    "    print(\"nHat =\", nHat)\n",
    "    \n",
    "    # Electric current K = nHat x H\n",
    "    Kx = nHat[1] * Hz - nHat[2] * Hy\n",
    "    Ky = nHat[2] * Hx - nHat[0] * Hz\n",
    "    Kz = nHat[0] * Hy - nHat[1] * Hx\n",
    "    \n",
    "    # Magnetic current N = - nHat x E\n",
    "    Nx = nHat[2] * Ey - nHat[1] * Ez\n",
    "    Ny = nHat[0] * Ez - nHat[2] * Ex\n",
    "    Nz = nHat[1] * Ex - nHat[0] * Ey\n",
    "    \n",
    "    # Source components\n",
    "    components = {\n",
    "        mp.Ex: Kx,\n",
    "        mp.Ey: Ky,\n",
    "        mp.Ez: Kz,\n",
    "        mp.Hx: Nx,\n",
    "        mp.Hy: Ny,\n",
    "        mp.Hz: Nz\n",
    "    }\n",
    "    \n",
    "    # Make sources\n",
    "    sources = [\n",
    "        mp.Source(src, field_comp, center=center, size=size, amp_data=source_comp)\n",
    "        for field_comp, source_comp in components.items() if np.sum(np.abs(source_comp))\n",
    "    ]\n",
    "    \n",
    "    return sources"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a571d673-9131-4fe6-89fd-f5e7a3c84575",
   "metadata": {},
   "outputs": [],
   "source": [
    "## Process input .nfd\n",
    "\n",
    "# File location\n",
    "file = \"SIWave Procedure Documentation_0.nfd\"\n",
    "\n",
    "# Import field data\n",
    "DAT = np.genfromtxt(file, dtype=None, skip_header=1, usecols=range(1,16), delimiter=\",\")\n",
    "\n",
    "# Parse field data\n",
    "x0 = np.empty(len(DAT))\n",
    "y0 = np.empty(len(DAT))\n",
    "z0 = np.empty(len(DAT))\n",
    "Ex0 = np.zeros(len(DAT), dtype=complex)\n",
    "Ey0 = np.zeros(len(DAT), dtype=complex)\n",
    "Ez0 = np.zeros(len(DAT), dtype=complex)\n",
    "Hx0 = np.zeros(len(DAT), dtype=complex)\n",
    "Hy0 = np.zeros(len(DAT), dtype=complex)\n",
    "Hz0 = np.zeros(len(DAT), dtype=complex)\n",
    "for i in range(len(DAT)):\n",
    "    x0[i] = float(DAT[i][0].decode('UTF-8').replace(\"mm\", \"\"))\n",
    "    y0[i] = float(DAT[i][1].decode('UTF-8').replace(\"mm\", \"\"))\n",
    "    z0[i] = float(DAT[i][2].decode('UTF-8').replace(\"mm\", \"\"))\n",
    "    Ex0[i] = DAT[i][3] + 1j * DAT[i][4]\n",
    "    Ey0[i] = DAT[i][5] + 1j * DAT[i][6]\n",
    "    Ez0[i] = DAT[i][7] + 1j * DAT[i][8]\n",
    "    Hx0[i] = DAT[i][9] + 1j * DAT[i][10]\n",
    "    Hy0[i] = DAT[i][11] + 1j * DAT[i][12]\n",
    "    Hz0[i] = DAT[i][13] + 1j * DAT[i][14]\n",
    "\n",
    "z0 = z0 - 2 # apply offset so that top face is at z = 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4af9ea83-6c40-4745-9c18-705686617df0",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Initialize sources\n",
    "sources = []\n",
    "\n",
    "# static src\n",
    "src = mp.ContinuousSource(frequency=freq)\n",
    "\n",
    "# Grid for interpolation\n",
    "Nx, Ny = 500, 500\n",
    "x      = np.linspace(min(x0), max(x0), Nx)\n",
    "y      = np.linspace(min(y0), max(y0), Ny)\n",
    "x, y   = np.meshgrid(x, y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "68f367f2-afdf-4738-b93d-4fc4a50075ab",
   "metadata": {},
   "outputs": [],
   "source": [
    "# # Process Face 1 of 3\n",
    "# zi     = np.unique(z0)[0]\n",
    "# center = mp.Vector3(z=zi)\n",
    "\n",
    "# # grab fields at current zi face\n",
    "# ind = np.where(z0 == zi)\n",
    "# xi  = x0[ind]\n",
    "# yi  = y0[ind]\n",
    "# Exi = Ex0[ind]\n",
    "# Eyi = Ey0[ind]\n",
    "# Ezi = Ez0[ind]\n",
    "# Hxi = Hx0[ind]\n",
    "# Hyi = Hy0[ind]\n",
    "# Hzi = Hz0[ind]\n",
    "\n",
    "# # define size\n",
    "# Sxi  = max(xi) - min(xi)\n",
    "# Syi  = max(yi) - min(yi)\n",
    "# size = mp.Vector3(Sxi, Syi)\n",
    "\n",
    "# # Plot raw data\n",
    "# fig, ax = plt.subplots(ncols=2, dpi=140, constrained_layout=True)\n",
    "\n",
    "# im = ax[0].scatter(xi, yi, c=Exi.real)\n",
    "# ax[0].set_aspect('equal', 'box')\n",
    "# fig.colorbar(im, ax=ax[0], shrink=0.4)\n",
    "# ax[0].set_xlim([min(xi), max(xi)])\n",
    "# ax[0].set_ylim([min(yi), max(yi)])\n",
    "# ax[0].set_xlabel(\"x (mm)\")\n",
    "# ax[0].set_ylabel(\"y (mm)\")\n",
    "# ax[0].set_title(\"(Raw Data) z = \" + str(zi))\n",
    "\n",
    "# # Interpolate data\n",
    "# Nx, Ny = 256, 256\n",
    "# x = np.linspace(min(xi), max(xi), Nx)\n",
    "# y = np.linspace(min(yi), max(yi), Ny)\n",
    "# x, y = np.meshgrid(x, y)\n",
    "\n",
    "# # Interpolate using delaunay triangularization \n",
    "# Exi = griddata((xi, yi), Exi, (x, y))\n",
    "# Eyi = griddata((xi, yi), Eyi, (x, y))\n",
    "# Ezi = griddata((xi, yi), Ezi, (x, y))\n",
    "# Hxi = griddata((xi, yi), Hxi, (x, y))\n",
    "# Hyi = griddata((xi, yi), Hyi, (x, y))\n",
    "# Hzi = griddata((xi, yi), Hzi, (x, y))\n",
    "\n",
    "# # Plot interpolation\n",
    "# im = ax[1].pcolormesh(x, y, Exi.real)\n",
    "# ax[1].set_aspect('equal', 'box')\n",
    "# fig.colorbar(im, ax=ax[1], shrink=0.4)\n",
    "# ax[1].set_xlim([np.amin(x), np.amax(x)])\n",
    "# ax[1].set_ylim([np.amin(y), np.amax(y)])\n",
    "# ax[1].set_xlabel(\"x (mm)\")\n",
    "# ax[1].set_ylabel(\"y (mm)\")\n",
    "# ax[1].set_title(\"(Interpolated Data)\")\n",
    "# plt.show()\n",
    "\n",
    "# field = np.ndarray(shape=(6,Nx,Ny,1), dtype=complex)\n",
    "# for i in range(Nx):\n",
    "#     for j in range(Ny):\n",
    "#         field[0,i,j,0] = Exi[i,j]\n",
    "#         field[1,i,j,0] = Eyi[i,j]\n",
    "#         field[2,i,j,0] = Ezi[i,j]\n",
    "#         field[3,i,j,0] = Hxi[i,j]\n",
    "#         field[4,i,j,0] = Hyi[i,j]\n",
    "#         field[5,i,j,0] = Hzi[i,j]\n",
    "\n",
    "# # update sources\n",
    "# sources += get_equiv_sources(field, src, center, size)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7d324786-0a05-40dd-acf6-8cf7307bb313",
   "metadata": {},
   "outputs": [],
   "source": [
    "# # Process Face 2 of 3\n",
    "# zi     = np.unique(z0)[1]\n",
    "# center = mp.Vector3(z=zi)\n",
    "\n",
    "# # grab fields at current zi face\n",
    "# ind = np.where(z0 == zi)\n",
    "# xi  = x0[ind]\n",
    "# yi  = y0[ind]\n",
    "# Exi = Ex0[ind]\n",
    "# Eyi = Ey0[ind]\n",
    "# Ezi = Ez0[ind]\n",
    "# Hxi = Hx0[ind]\n",
    "# Hyi = Hy0[ind]\n",
    "# Hzi = Hz0[ind]\n",
    "\n",
    "# # define size\n",
    "# Sxi  = max(xi) - min(xi)\n",
    "# Syi  = max(yi) - min(yi)\n",
    "# size = mp.Vector3(Sxi, Syi)\n",
    "\n",
    "# # Plot raw data\n",
    "# fig, ax = plt.subplots(ncols=2, dpi=140, constrained_layout=True)\n",
    "\n",
    "# im = ax[0].scatter(xi, yi, c=Exi.real)\n",
    "# ax[0].set_aspect('equal', 'box')\n",
    "# fig.colorbar(im, ax=ax[0], shrink=0.4)\n",
    "# ax[0].set_xlim([min(xi), max(xi)])\n",
    "# ax[0].set_ylim([min(yi), max(yi)])\n",
    "# ax[0].set_xlabel(\"x (mm)\")\n",
    "# ax[0].set_ylabel(\"y (mm)\")\n",
    "# ax[0].set_title(\"(Raw Data) z = \" + str(zi))\n",
    "\n",
    "# # Interpolate data\n",
    "# Nx, Ny = 256, 256\n",
    "# x = np.linspace(min(xi), max(xi), Nx)\n",
    "# y = np.linspace(min(yi), max(yi), Ny)\n",
    "# x, y = np.meshgrid(x, y)\n",
    "\n",
    "# # Interpolate using delaunay triangularization \n",
    "# Exi = griddata((xi, yi), Exi, (x, y))\n",
    "# Eyi = griddata((xi, yi), Eyi, (x, y))\n",
    "# Ezi = griddata((xi, yi), Ezi, (x, y))\n",
    "# Hxi = griddata((xi, yi), Hxi, (x, y))\n",
    "# Hyi = griddata((xi, yi), Hyi, (x, y))\n",
    "# Hzi = griddata((xi, yi), Hzi, (x, y))\n",
    "\n",
    "# # Plot interpolation\n",
    "# im = ax[1].pcolormesh(x, y, Exi.real)\n",
    "# ax[1].set_aspect('equal', 'box')\n",
    "# fig.colorbar(im, ax=ax[1], shrink=0.4)\n",
    "# ax[1].set_xlim([np.amin(x), np.amax(x)])\n",
    "# ax[1].set_ylim([np.amin(y), np.amax(y)])\n",
    "# ax[1].set_xlabel(\"x (mm)\")\n",
    "# ax[1].set_ylabel(\"y (mm)\")\n",
    "# ax[1].set_title(\"(Interpolated Data)\")\n",
    "# plt.show()\n",
    "\n",
    "# field = np.ndarray(shape=(6,Nx,Ny,1), dtype=complex)\n",
    "# for i in range(Nx):\n",
    "#     for j in range(Ny):\n",
    "#         field[0,i,j,0] = Exi[i,j]\n",
    "#         field[1,i,j,0] = Eyi[i,j]\n",
    "#         field[2,i,j,0] = Ezi[i,j]\n",
    "#         field[3,i,j,0] = Hxi[i,j]\n",
    "#         field[4,i,j,0] = Hyi[i,j]\n",
    "#         field[5,i,j,0] = Hzi[i,j]\n",
    "\n",
    "# # update sources\n",
    "# sources += get_equiv_sources(field, src, center, size)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "acaac74c-7124-4b12-aee0-fa2924d0467f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Process Face 3 of 3\n",
    "zi     = np.unique(z0)[2]\n",
    "center = mp.Vector3(z=zi)\n",
    "\n",
    "# grab fields at current zi face\n",
    "ind = np.where(z0 == zi)\n",
    "xi  = x0[ind]\n",
    "yi  = y0[ind]\n",
    "Exi = Ex0[ind]\n",
    "Eyi = Ey0[ind]\n",
    "Ezi = Ez0[ind]\n",
    "Hxi = Hx0[ind]\n",
    "Hyi = Hy0[ind]\n",
    "Hzi = Hz0[ind]\n",
    "\n",
    "# define size\n",
    "Sxi  = max(xi) - min(xi)\n",
    "Syi  = max(yi) - min(yi)\n",
    "size = mp.Vector3(Sxi, Syi)\n",
    "\n",
    "# Plot raw data\n",
    "fig, ax = plt.subplots(ncols=2, dpi=140, constrained_layout=True)\n",
    "\n",
    "im = ax[0].scatter(xi, yi, c=Exi.real)\n",
    "ax[0].set_aspect('equal', 'box')\n",
    "fig.colorbar(im, ax=ax[0], shrink=0.4)\n",
    "ax[0].set_xlim([min(xi), max(xi)])\n",
    "ax[0].set_ylim([min(yi), max(yi)])\n",
    "ax[0].set_xlabel(\"x (mm)\")\n",
    "ax[0].set_ylabel(\"y (mm)\")\n",
    "ax[0].set_title(\"(Raw Data) z = \" + str(zi))\n",
    "\n",
    "# Interpolate using delaunay triangularization \n",
    "Exi = griddata((xi, yi), Exi, (x, y), method='cubic')\n",
    "Eyi = griddata((xi, yi), Eyi, (x, y), method='cubic')\n",
    "Ezi = griddata((xi, yi), Ezi, (x, y), method='cubic')\n",
    "Hxi = griddata((xi, yi), Hxi, (x, y), method='cubic')\n",
    "Hyi = griddata((xi, yi), Hyi, (x, y), method='cubic')\n",
    "Hzi = griddata((xi, yi), Hzi, (x, y), method='cubic')\n",
    "\n",
    "# Plot interpolation\n",
    "im = ax[1].pcolormesh(x, y, Exi.real)\n",
    "ax[1].set_aspect('equal', 'box')\n",
    "fig.colorbar(im, ax=ax[1], shrink=0.4)\n",
    "# ax[1].set_xlim([np.amin(x), np.amax(x)])\n",
    "# ax[1].set_ylim([np.amin(y), np.amax(y)])\n",
    "ax[1].set_xlabel(\"x (mm)\")\n",
    "ax[1].set_ylabel(\"y (mm)\")\n",
    "ax[1].set_title(\"(Interpolated Data)\")\n",
    "plt.show()\n",
    "\n",
    "field = np.ndarray(shape=(6,Nx,Ny,1), dtype=complex)\n",
    "field[0,:,:,0] = Exi.T\n",
    "field[1,:,:,0] = Eyi.T\n",
    "field[2,:,:,0] = Ezi.T\n",
    "field[3,:,:,0] = Hxi.T\n",
    "field[4,:,:,0] = Hyi.T\n",
    "field[5,:,:,0] = Hzi.T\n",
    "\n",
    "# update sources\n",
    "sources += get_equiv_sources(field, src, center, size)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "54d1588e-7ece-4eb1-8388-b0b81419d1ed",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Simulation\n",
    "sim = mp.Simulation(cell_size=cell_size,\n",
    "                    boundary_layers=pml_layers,\n",
    "                    sources=sources,\n",
    "                    force_complex_fields=False,\n",
    "                    resolution=resolution)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "51d7c809-ddce-41c2-991f-9003722f961a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Test plot of domain\n",
    "vol = mp.Volume(center=mp.Vector3(), size=mp.Vector3(Sx, Sy, 0))\n",
    "\n",
    "sim.plot2D(output_plane=vol)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9f3f3192-8089-4527-b439-3cb98cf603a6",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Test plot of domain\n",
    "vol = mp.Volume(center=mp.Vector3(), size=mp.Vector3(Sx, 0, Sz))\n",
    "\n",
    "sim.plot2D(output_plane=vol)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4df074a5-fbd4-4254-b2d8-149511517fab",
   "metadata": {},
   "outputs": [],
   "source": [
    "## Run simulation\n",
    "\n",
    "# Field probe locations\n",
    "vol = mp.Volume(center=mp.Vector3(z=25), size=mp.Vector3(500, 500, 0))\n",
    "pt  = mp.Vector3(z=25)\n",
    "\n",
    "# Simulation time properties\n",
    "T  = 1/freq\n",
    "tf = 5 * T\n",
    "dt = T / 20\n",
    "\n",
    "Ex = []\n",
    "Ey = []\n",
    "Ez = []\n",
    "Ext = []\n",
    "Eyt = []\n",
    "Ezt = []\n",
    "def updateFields(sim):\n",
    "    Ex.append(sim.get_array(mp.Ex, vol, cmplx=False))\n",
    "    Ey.append(sim.get_array(mp.Ey, vol, cmplx=False))\n",
    "    Ez.append(sim.get_array(mp.Ez, vol, cmplx=False))\n",
    "    Ext.append(sim.get_array(center=pt, size=mp.Vector3(), component=mp.Ex))\n",
    "    Eyt.append(sim.get_array(center=pt, size=mp.Vector3(), component=mp.Ey))\n",
    "    Ezt.append(sim.get_array(center=pt, size=mp.Vector3(), component=mp.Ez))\n",
    "\n",
    "sim.run(mp.at_every(dt, updateFields), until=tf)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "736f3a97-ca42-487a-9af4-7ad7845f96cc",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Find Peaks (discard first half of data due to artificial peaks at start)\n",
    "ind = int(len(Ex)/2)\n",
    "\n",
    "Ex_peak = np.amax(np.abs(Ex), axis=0) # need abs here?\n",
    "Ey_peak = np.amax(np.abs(Ey), axis=0)\n",
    "Ez_peak = np.amax(np.abs(Ez), axis=0)\n",
    "\n",
    "# magnitude\n",
    "E_peak = np.sqrt(Ex_peak**2 + Ey_peak**2 + Ez_peak**2)\n",
    "\n",
    "# Report peak magnitude\n",
    "print(\"Peak field: %.2f V/m\" % (np.amax(E_peak)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b7a13c73-dac1-436d-a1f4-2861f39ebd82",
   "metadata": {},
   "outputs": [],
   "source": [
    "# define axes\n",
    "x = np.linspace(-vol.size[0]/2, vol.size[0]/2, np.shape(Ex)[1])\n",
    "y = np.linspace(-vol.size[1]/2, vol.size[1]/2, np.shape(Ex)[2])\n",
    "\n",
    "# Plot\n",
    "fig, ax = plt.subplots(dpi=150)\n",
    "im = ax.pcolormesh(x, y, E_peak, cmap='jet', vmin=0.228, vmax=23.2389)\n",
    "fig.colorbar(im, ticks=np.linspace(0.228,23.2389, 16))\n",
    "ax.set_xlabel(\"x (mm)\")\n",
    "ax.set_ylabel(\"y (mm)\")\n",
    "plt.axis([min(x), max(x), min(y), max(y)])\n",
    "ax.set_aspect('equal', 'box')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "69b7bbaa-ebc9-4ffc-8e9f-e178d9eaed4e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot Time-Domain\n",
    "fig, ax = plt.subplots(ncols=2, figsize=(14,5), constrained_layout=True)\n",
    "\n",
    "ax[0].plot(Ext[ind:], '.-', label=\"Ex\")\n",
    "ax[0].plot(Eyt[ind:], '.-', label=\"Ey\")\n",
    "ax[0].plot(Ezt[ind:], '.-', label=\"Ez\")\n",
    "ax[0].legend()\n",
    "ax[0].grid()\n",
    "\n",
    "Ext_mag = np.sqrt(np.asarray(Ext)**2 + np.asarray(Eyt)**2 + np.asarray(Ezt)**2)\n",
    "ax[1].plot(Ext_mag[ind:], '.-')\n",
    "ax[1].grid()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "118b38dd-b228-440f-b8aa-21fdced5a5ff",
   "metadata": {},
   "outputs": [],
   "source": [
    "## Animate fields\n",
    "\n",
    "# calculate E_mag at each time\n",
    "E_mag = -np.inf * np.ones((np.shape(Ex)))\n",
    "for i in range(np.shape(Ex)[0]):\n",
    "    E_mag[i] = np.sqrt(Ex[i]**2 + Ey[i]**2 + Ez[i]**2)\n",
    "\n",
    "# Animate\n",
    "fig, ax = plt.subplots(dpi=150)\n",
    "ims = []\n",
    "for i in range(len(E_mag)):\n",
    "    im = ax.pcolormesh(x, y, E_mag[i], animated=True, cmap='jet',\n",
    "                       vmin=0.0287, vmax=23.2389)\n",
    "                       # norm=matplotlib.colors.LogNorm(vmin=1., vmax=100.))\n",
    "    ax.set(xlim=(x[0], x[-1]), ylim=(y[0], y[-1]))\n",
    "    ax.set_xlabel(\"x (mm)\")\n",
    "    ax.set_ylabel(\"y (mm)\")\n",
    "    ax.set_aspect('equal', 'box')\n",
    "    ims.append([im])\n",
    "\n",
    "# cbar = fig.colorbar(im)\n",
    "cbar = fig.colorbar(im, ticks=np.linspace(0.0287, 23.2389, 16))\n",
    "# cbar.set_clims(1, 20)\n",
    "# cbar.set_ticks(np.arange(1, 100, 20))\n",
    "# cbar.set_label(\"hi\")\n",
    "\n",
    "ani = animation.ArtistAnimation(fig, ims, interval=200, blit=True, repeat_delay=1000)\n",
    "ani.save(\"test.mp4\")\n",
    "Video(\"test.mp4\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b4ed809f-96ed-485c-89cd-30098d20f08c",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.10.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}