MagDipole2LoopsCoupling.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# MagDipoleLoops3D\n",
    "### An interactive notebook application for learning the magnetic flux  linkage between two loops\n",
    "\n",
    "* Plot loops and the dipole directions in 3D\n",
    "* Plot magnetic field lines from a magnetic dipole (a small loop) in 3D\n",
    "* Explore how two loops interact through magnetic flux linkage in 3D\n",
    "* \"Run All\" to get the default result of Loop 1 transmitting and Loop 2 receiving\n",
    "* Plot many loops and lines the way you like by changing the plotting codes\n",
    "\n",
    "**Bugs and improvement: yangdikun@gmail.com**\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Import packages\n",
    "* Note the 3D plotting requires **plotly**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from ipywidgets import interact, interactive, fixed, Layout\n",
    "import ipywidgets as widgets"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3.6.1\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<script type=\"text/javascript\">window.PlotlyConfig = {MathJaxConfig: 'local'};</script><script type=\"text/javascript\">if (window.MathJax) {MathJax.Hub.Config({SVG: {font: \"STIX-Web\"}});}</script><script>requirejs.config({paths: { 'plotly': ['https://cdn.plot.ly/plotly-latest.min']},});if(!window._Plotly) {require(['plotly'],function(plotly) {window._Plotly=plotly;});}</script>"
      ],
      "text/vnd.plotly.v1+html": [
       "<script type=\"text/javascript\">window.PlotlyConfig = {MathJaxConfig: 'local'};</script><script type=\"text/javascript\">if (window.MathJax) {MathJax.Hub.Config({SVG: {font: \"STIX-Web\"}});}</script><script>requirejs.config({paths: { 'plotly': ['https://cdn.plot.ly/plotly-latest.min']},});if(!window._Plotly) {require(['plotly'],function(plotly) {window._Plotly=plotly;});}</script>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import plotly.plotly as py\n",
    "import plotly.graph_objs as go\n",
    "from plotly import __version__\n",
    "from plotly.offline import download_plotlyjs, init_notebook_mode, plot, iplot\n",
    "print(__version__) # requires version >= 1.9.0\n",
    "init_notebook_mode(connected=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Function definitions\n",
    "* VerticalMagneticDipoleField: Calculate the magnetic flux intensity (B field) from a vertical magnetic dipole (a horizontal loop) at the origin\n",
    "* VerticalMagneticDipoleLine: Calculate the paths of B field lines from a vertical magnetic dipole at the origin\n",
    "* MagneticDipoleLine: Calculate the paths of B field lines from a magnetic dipole with an arbitrary orientation and location.\n",
    "* MagneticDipoleLoop: Calculate the paths of an arbitrarily oriented and located loop and its normal direction."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# x, y, z: location of the observing point\n",
    "# Bx, By, Bz: vector B field at the observing point in vaccum\n",
    "def VerticalMagneticDipoleField(x, y, z):\n",
    "    r = (x**2 + y**2 + z**2)**(0.5)\n",
    "    Bx = 1e-7/r**3 * 3*x*z/r**2\n",
    "    By = 1e-7/r**3 * 3*y*z/r**2\n",
    "    Bz = 1e-7/r**3 * (3*z**2/r**2 - 1)\n",
    "    return Bx, By, Bz\n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# radius: a float defining how far the field lines expand\n",
    "# stepSize: length of the segments that make the field lines\n",
    "# yloc, zloc: a single field line path in the y-z plane\n",
    "def VerticalMagneticDipoleLine(radius,stepSize=1.):\n",
    "    x, y, z = 0., [radius], [-stepSize/4.]\n",
    "    while (y[-1]>=0.) & (z[-1] <= 0.):\n",
    "        _, By, Bz = VerticalMagneticDipoleField(x, y[-1], z[-1])\n",
    "        B = (By**2 + Bz**2)**(0.5)\n",
    "        By, Bz = By/B*stepSize, Bz/B*stepSize\n",
    "        y.append(y[-1]+By)\n",
    "        z.append(z[-1]+Bz)\n",
    "    y.pop()\n",
    "    z.pop()\n",
    "    yloc = np.r_[y[-1:0:-1], y]\n",
    "    zloc = np.r_[-np.array(z[-1:0:-1]),z]\n",
    "    return yloc, zloc\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# dipoleLoc: (x,y,z) of a magnetic dipole\n",
    "# dipoleDec, dipoleInc: declination and inclination angles of the dipole direction in degree\n",
    "# radii: a float vector defining the radii of the field line shells\n",
    "# naz: number of azimuths for the plotting of field lines\n",
    "# stepSize: length of the segments that make the field lines\n",
    "# xloc, yloc, zloc: a collection of all field line paths \n",
    "def MagneticDipoleLine(dipoleLoc,dipoleDec,dipoleInc,radii,naz=10,stepSize=1.):\n",
    "    x0, y0, z0 = dipoleLoc[0], dipoleLoc[1], dipoleLoc[2]\n",
    "    \n",
    "    # rotation matrix\n",
    "    theta, alpha = -np.pi*(dipoleInc+90.)/180., -np.pi*dipoleDec/180.\n",
    "    Rx = np.array([[1.,0.,0.],[0.,np.cos(theta),-np.sin(theta)],[0.,np.sin(theta),np.cos(theta)]])\n",
    "    Rz = np.array([[np.cos(alpha),-np.sin(alpha),0.],[np.sin(alpha),np.cos(alpha),0.],[0.,0.,1.]])\n",
    "    R = np.dot(Rz,Rx) # Rz @ Rx\n",
    "\n",
    "    azimuth = np.linspace(0.,2*np.pi,num=naz,endpoint=False)\n",
    "    xloc, yloc, zloc = [], [], []\n",
    "    for r in radii:\n",
    "        hloc, vloc = VerticalMagneticDipoleLine(r,stepSize)\n",
    "        for a in azimuth:\n",
    "            x, y, z = np.sin(a)*hloc, np.cos(a)*hloc, vloc\n",
    "            xyz = np.dot(R, np.vstack((x,y,z)))\n",
    "            xloc.append(xyz[0]+x0)\n",
    "            yloc.append(xyz[1]+y0)\n",
    "            zloc.append(xyz[2]+z0)\n",
    "    return xloc, yloc, zloc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# dipoleLoc: (x,y,z) of a magnetic dipole\n",
    "# dipoleDec, dipoleInc: declination and inclination angles of the dipole direction in degree\n",
    "# radius: radius of loop surface\n",
    "# nseg: number of line segments that make a loop\n",
    "# xloc, yloc, zloc: a collection of line paths for a loop and its normal direction\n",
    "def MagneticDipoleLoop(dipoleLoc,dipoleDec,dipoleInc,radius,nseg=20):\n",
    "    x0, y0, z0 = dipoleLoc[0], dipoleLoc[1], dipoleLoc[2]\n",
    "    \n",
    "    # rotation matrix\n",
    "    theta, alpha = -np.pi*(dipoleInc+90.)/180., -np.pi*dipoleDec/180.\n",
    "    Rx = np.array([[1.,0.,0.],[0.,np.cos(theta),-np.sin(theta)],[0.,np.sin(theta),np.cos(theta)]])\n",
    "    Rz = np.array([[np.cos(alpha),-np.sin(alpha),0.],[np.sin(alpha),np.cos(alpha),0.],[0.,0.,1.]])\n",
    "    R = np.dot(Rz,Rx) # Rz @ Rx\n",
    "    # loop path\n",
    "    azimuth = np.linspace(0.,2*np.pi,num=nseg+1,endpoint=True)\n",
    "    x, y, z = np.sin(azimuth)*radius, np.cos(azimuth)*radius, azimuth*0.\n",
    "    xyz = np.dot(R, np.vstack((x,y,z)))\n",
    "    xloc1, yloc1, zloc1 = xyz[0]+x0, xyz[1]+y0, xyz[2]+z0\n",
    "    # normal direction pointer\n",
    "    x, y, z = [0.,0.],[0.,0.],[0.,radius]\n",
    "    xyz = np.dot(R, np.vstack((x,y,z)))\n",
    "    xloc2, yloc2, zloc2 = xyz[0]+x0, xyz[1]+y0, xyz[2]+z0\n",
    "    return [xloc1,xloc2], [yloc1,yloc2], [zloc1,zloc2]\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plotting code\n",
    "* Loop 1 is fixed to the origin\n",
    "* 8 azimuths (field lines) around the dipole\n",
    "* Other parameters adjustable through sliders\n",
    "* Choose what to plot: Loop 1 and/or field lines of Loop 1 and/or Loop 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "db171d9d5aa744f2b8882a81da3b6316",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(IntSlider(value=0, description='Loop 1 dec', max=360), IntSlider(value=-90, description=…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<function __main__.plot_magDipoleEM(x1, y1, z1, dec1, inc1, nazimuth, nshell, rmax, x2, y2, z2, dec2, inc2, loop1flag, line1flag, loop2flag)>"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def plot_magDipoleEM(x1,y1,z1,dec1,inc1,nazimuth,nshell,rmax,x2,y2,z2,dec2,inc2,loop1flag,line1flag,loop2flag):\n",
    "    # internal parameters\n",
    "    loopRadius1, loopRadius2 = 0.5, 0.5\n",
    "    stepSize = 0.2 # length of field line segment    \n",
    "    # get line path for the 1st loop\n",
    "    loopx1, loopy1, loopz1 = MagneticDipoleLoop((x1,y1,z1),dec1,inc1,loopRadius1)\n",
    "    # get line path for the field lines of the 1st loop\n",
    "    radii = np.linspace(rmax,0.,nshell,endpoint=False)\n",
    "    linex1, liney1, linez1 = MagneticDipoleLine((x1,y1,z1),dec1,inc1,radii,nazimuth,stepSize)\n",
    "    # get line path for the 2nd loop\n",
    "    loopx2, loopy2, loopz2 = MagneticDipoleLoop((x2,y2,z2),dec2,inc2,loopRadius2)\n",
    "\n",
    "    # create plot\n",
    "    data = []\n",
    "    if loop1flag:\n",
    "        for xx,yy,zz in zip(loopx1,loopy1,loopz1):\n",
    "            data.append(go.Scatter3d(x=xx, y=yy, z=zz, mode='lines', line=dict(color='red',width=3)))\n",
    "    if line1flag:\n",
    "        for xx,yy,zz in zip(linex1,liney1,linez1):\n",
    "            data.append(go.Scatter3d(x=xx, y=yy, z=zz, mode='lines', line=dict(color='blue',width=3)))\n",
    "    if loop2flag:\n",
    "        for xx,yy,zz in zip(loopx2,loopy2,loopz2):\n",
    "            data.append(go.Scatter3d(x=xx, y=yy, z=zz, mode='lines', line=dict(color='green',width=3)))\n",
    "    layout = dict(width=1000, height=500, autosize=False, showlegend=False, \n",
    "                  margin=go.layout.Margin(l=5,r=5,b=5,t=5,pad=4),paper_bgcolor='#efefef',\n",
    "                 scene = dict(aspectmode = 'data', aspectratio = dict(x=1, y=1, z=1)))\n",
    "    fig = dict(data=data, layout=layout)\n",
    "    iplot(fig)\n",
    "    return\n",
    "\n",
    "interact(plot_magDipoleEM, \n",
    "         x1=fixed(0.), y1=fixed(0.), z1=fixed(0.),\n",
    "         dec1=widgets.IntSlider(min=0,max=360,step=1,value=0,description='Loop 1 dec'),\n",
    "         inc1=widgets.IntSlider(min=-90,max=90,step=1,value=-90,description='Loop 1 inc'),\n",
    "         nazimuth=fixed(8),\n",
    "         nshell=widgets.IntSlider(min=1,max=5,step=1,value=2,description='Number of line layers',style = {'description_width': 'initial'}),\n",
    "         rmax=widgets.FloatSlider(min=1.,max=20.,step=1.,value=5.,description='Max. radius of lines',style = {'description_width': 'initial'}),\n",
    "         x2=widgets.FloatSlider(min=-20.,max=20.,step=1.,value=10.,description='Loop 2 X'),\n",
    "         y2=widgets.FloatSlider(min=-20.,max=20.,step=1.,value=0.,description='Loop 2 Y'),\n",
    "         z2=widgets.FloatSlider(min=-20.,max=20.,step=1.,value=0.,description='Loop 2 Z'),\n",
    "         dec2=widgets.IntSlider(min=0,max=360,step=1,value=0,description='Loop 2 dec'),\n",
    "         inc2=widgets.IntSlider(min=-90,max=90,step=1,value=-90,description='Loop 2 inc'),\n",
    "         loop1flag=widgets.Checkbox(value=True,disabled=False,description='Show Loop 1'),\n",
    "         line1flag=widgets.Checkbox(value=False,disabled=False,description='Show field lines of Loop 1'),\n",
    "         loop2flag=widgets.Checkbox(value=False,disabled=False,description='Show Loop 2')\n",
    "         )\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "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.7.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}