{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Coordinate systems"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Introduction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebooks provides examples of curvilinear coordinates. Particularly,\n",
    "it presents the coordinate surfaces for three common coordinate systems."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Use of PyVista in this notebook"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We use `PyVista` to represent the coordinate surfaces and `ìtkwidgets`\n",
    "to render them interactively in the notebook.\n",
    "\n",
    "If installing things manually, the following would be required\n",
    "\n",
    "\n",
    "    conda install -c conda-forge pyvista\n",
    "    conda install -c conda-forge itkwidgets\n",
    "\n",
    "\n",
    "**More information:** https://docs.pyvista.org/examples/index.html"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pyvista as pv\n",
    "from itkwidgets import view"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "red = (0.9, 0.1, 0.1)\n",
    "blue = (0.2, 0.5, 0.7)\n",
    "green = (0.3, 0.7, 0.3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Cylindrical coordinates\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The [ISO standard 80000-2](https://en.wikipedia.org/wiki/ISO_80000-2) recommends the use of $(ρ, φ, z)$, where $ρ$ is the radial coordinate, $\\varphi$ the azimuth, and $z$ the height.\n",
    "\n",
    "For the conversion between cylindrical and Cartesian coordinates, it is convenient to assume that the reference plane of the former is the Cartesian $xy$-plane (with equation $z=0$, and the cylindrical axis is the Cartesian $z$-axis. Then the $z$-coordinate is the same in both systems, and the correspondence between cylindrical $(\\rho, \\varphi)$ and Cartesian $(x, y)$ are the same as for polar coordinates, namely\n",
    "\n",
    "\\begin{align}\n",
    "x &= \\rho \\cos \\varphi \\\\\n",
    "y &= \\rho \\sin \\varphi\n",
    "\\end{align}\n",
    "\n",
    "in one direction, and\n",
    "\n",
    "\\begin{align}\n",
    "    \\rho &= \\sqrt{x^2+y^2} \\\\\n",
    "    \\varphi &= \\begin{cases}\n",
    "        0 & \\mbox{if } x = 0 \\mbox{ and } y = 0\\\\\n",
    "        \\arcsin\\left(\\frac{y}{\\rho}\\right) & \\mbox{if } x \\geq 0 \\\\\t\n",
    "        \\arctan\\left(\\frac{y}{x}\\right) & \\mbox{if } x > 0 \\\\\t\n",
    "        -\\arcsin\\left(\\frac{y}{\\rho}\\right) + \\pi & \\mbox{if } x < 0\n",
    "    \\end{cases}\n",
    "\\end{align}\n",
    "\n",
    "in the other. It is also common to use [$\\varphi = \\mathrm{atan2}(y, x)$](https://en.wikipedia.org/wiki/Atan2)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Cylinder\n",
    "phi, z = np.mgrid[0:2*np.pi:31j, -2:2*np.pi:31j]\n",
    "x = 2*np.cos(phi)\n",
    "y = 2*np.sin(phi)\n",
    "cylinder = pv.StructuredGrid(x, y, z)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Vertical plane\n",
    "rho, z = np.mgrid[0:3:31j, -2:2*np.pi:31j]\n",
    "x = rho*np.cos(np.pi/4)\n",
    "y = rho*np.sin(np.pi/4)\n",
    "vert_plane = pv.StructuredGrid(x, y, z)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Horizontal plane\n",
    "rho, phi = np.mgrid[0:3:31j, 0:2*np.pi:31j]\n",
    "x = rho*np.cos(phi)\n",
    "y = rho*np.sin(phi)\n",
    "z = np.ones_like(x)\n",
    "hor_plane = pv.StructuredGrid(x, y, z)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "29a27643a8434035b33c11db8da76bec",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Viewer(geometries=[{'vtkClass': 'vtkPolyData', 'points': {'vtkClass': 'vtkPoints', 'name': '_points', 'numberO…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "view(geometries=[cylinder, vert_plane, hor_plane],\n",
    "     geometry_colors=[blue, red, green])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Spherical coordinates"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The use of $(r, θ, φ)$ to denote radial distance, inclination (or elevation), and azimuth, respectively, is common practice in physics, and is specified by [ISO standard 80000-2](https://en.wikipedia.org/wiki/ISO_80000-2).\n",
    "\n",
    "\n",
    "The spherical coordinates of a point can be obtained from its [Cartesian coordinate system](https://en.wikipedia.org/wiki/ISO_80000-2) $(x, y, z)$\n",
    "\n",
    "\\begin{align}\n",
    "r&=\\sqrt{x^2 + y^2 + z^2} \\\\\n",
    "\\theta &= \\arccos\\frac{z}{\\sqrt{x^2 + y^2 + z^2}} = \\arccos\\frac{z}{r} \\\\\n",
    "\\varphi &= \\arctan \\frac{y}{x}\n",
    "\\end{align}\n",
    "\n",
    "The inverse tangent denoted in $\\varphi = \\arctan\\left(\\frac{y}{x}\\right)$ must be suitably defined , taking into account the correct quadrant of $(x, y)$ (using the function ``atan2``).\n",
    "\n",
    "Conversely, the Cartesian coordinates may be retrieved from the spherical coordinates (_radius_ $r$, _inclination_ $\\theta$, _azimuth_ $\\varphi$), where $r \\in [0, \\infty)$, $\\theta \\in [0, \\pi]$ and $\\varphi \\in [0, 2\\pi)$, by:\n",
    "\n",
    "\\begin{align}\n",
    "x&=r \\, \\sin\\theta \\, \\cos\\varphi \\\\\n",
    "y&=r \\, \\sin\\theta \\, \\sin\\varphi \\\\\n",
    "z&=r \\, \\cos\\theta\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "theta, phi = np.mgrid[0:np.pi:21j, 0:np.pi:21j]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Sphere\n",
    "x = np.sin(phi) * np.cos(theta)\n",
    "y = np.sin(phi) * np.sin(theta)\n",
    "z = np.cos(phi)\n",
    "sphere = pv.StructuredGrid(x, y, z)\n",
    "sphere2 =  pv.StructuredGrid(-x, -y, z)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Cone\n",
    "x = theta/3 * np.cos(phi)\n",
    "y = theta/3 * np.sin(phi)\n",
    "z = theta/3\n",
    "cone = pv.StructuredGrid(x, y, z)\n",
    "cone2 = pv.StructuredGrid(-x, -y, z)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plane\n",
    "x = theta/np.pi\n",
    "y = theta/np.pi\n",
    "z = phi - np.pi/2\n",
    "plane = pv.StructuredGrid(x, y, z)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "99f523032243438d93c4431e154dc150",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Viewer(geometries=[{'vtkClass': 'vtkPolyData', 'points': {'vtkClass': 'vtkPoints', 'name': '_points', 'numberO…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "view(geometries=[sphere, sphere2, cone, cone2, plane],\n",
    "     geometry_colors=[blue, blue, red, red, green],\n",
    "     geometry_opacities=[1.0, 0.5, 1.0, 0.5, 1.0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Ellipsoidal coordinates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [],
   "source": [
    "v, theta = np.mgrid[0:np.pi/2:21j, 0:np.pi:21j]\n",
    "a = 3\n",
    "b = 2\n",
    "c = 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Ellipsoid\n",
    "lam = 3\n",
    "x =  np.sqrt(a**2  + lam) * np.cos(v) * np.cos(theta)\n",
    "y = np.sqrt(b**2  + lam)* np.cos(v) * np.sin(theta)\n",
    "z = np.sqrt(c**2 + lam) * np.sin(v)\n",
    "ellipsoid = pv.StructuredGrid(x, y, z)\n",
    "ellipsoid2 = pv.StructuredGrid(x, y, -z)\n",
    "ellipsoid3 = pv.StructuredGrid(x, -y, z)\n",
    "ellipsoid4 = pv.StructuredGrid(x, -y, -z)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Hyperboloid of one sheet\n",
    "mu = 2\n",
    "x = np.sqrt(a**2 + mu) * np.cosh(v) * np.cos(theta)\n",
    "y = np.sqrt(b**2 + mu) * np.cosh(v) * np.sin(theta)\n",
    "z = np.sqrt(c**2 + mu) * np.sinh(v)\n",
    "z = np.sqrt(c**2 + mu) * np.sinh(v)\n",
    "hyper = pv.StructuredGrid(x, y, z)\n",
    "hyper2 = pv.StructuredGrid(x, y, -z)\n",
    "hyper3 = pv.StructuredGrid(x, -y, z)\n",
    "hyper4 = pv.StructuredGrid(x, -y, -z)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Hyperboloid of two sheets\n",
    "nu = 1\n",
    "x = np.sqrt(a**2 + nu) * np.cosh(v)\n",
    "y = np.sqrt(c**2 + nu) * np.sinh(v) * np.sin(theta)\n",
    "z = np.sqrt(b**2 + nu) * np.sinh(v) * np.cos(theta)\n",
    "hyper_up = pv.StructuredGrid(x, y, z)\n",
    "hyper_down = pv.StructuredGrid(-x, y, z)\n",
    "hyper_up2 = pv.StructuredGrid(x, -y, z)\n",
    "hyper_down2 = pv.StructuredGrid(-x, -y, z)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "c7d007e315d6443eae22c9b9655b22e1",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Viewer(geometries=[{'vtkClass': 'vtkPolyData', 'points': {'vtkClass': 'vtkPoints', 'name': '_points', 'numberO…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "view(geometries=[ellipsoid, ellipsoid2, ellipsoid3, ellipsoid4,\n",
    "                 hyper, hyper2, hyper3, hyper4,\n",
    "                 hyper_up, hyper_down, hyper_up2, hyper_down2],\n",
    "     geometry_colors=[blue]*4 + [red]*4 + [green]*4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## References"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1. Wikipedia contributors. [\"Cylindrical coordinate system.\"](https://en.wikipedia.org/wiki/Cylindrical_coordinate_system) Wikipedia, The Free Encyclopedia. Wikipedia, The Free Encyclopedia, 12 Dec. 2016. Web. 9 Feb. 2017.\n",
    "\n",
    "2. Wikipedia contributors. [\"Spherical coordinate system.\"](https://en.wikipedia.org/wiki/Spherical_coordinate_system) Wikipedia, The Free Encyclopedia. Wikipedia, The Free Encyclopedia, 29 Jan. 2017. Web. 9 Feb. 2017.\n",
    "\n",
    "3. Sullivan et al., (2019). [PyVista: 3D plotting and mesh analysis through a streamlined interface for the Visualization Toolkit (VTK)](https://joss.theoj.org/papers/10.21105/joss.01450). Journal of Open Source Software, 4(37), 1450, https://doi.org/10.21105/joss.01450"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "\n",
       "<link href='http://fonts.googleapis.com/css?family=Fenix' rel='stylesheet' type='text/css'>\n",
       "<link href='http://fonts.googleapis.com/css?family=Alegreya+Sans:100,300,400,500,700,800,900,100italic,300italic,400italic,500italic,700italic,800italic,900italic' rel='stylesheet' type='text/css'>\n",
       "<link href='http://fonts.googleapis.com/css?family=Source+Code+Pro:300,400' rel='stylesheet' type='text/css'>\n",
       "<style>\n",
       "/* Based on Lorena Barba template available at: https://github.com/barbagroup/AeroPython/blob/master/styles/custom.css*/\n",
       "@font-face {\n",
       "font-family: \"Computer Modern\";\n",
       "src: url('http://mirrors.ctan.org/fonts/cm-unicode/fonts/otf/cmunss.otf');\n",
       "}\n",
       "div.cell{\n",
       "width:800px;\n",
       "margin-left:16% !important;\n",
       "margin-right:auto;\n",
       "}\n",
       "h1 {\n",
       "font-family: 'Alegreya Sans', sans-serif;\n",
       "}\n",
       "h2 {\n",
       "font-family: 'Fenix', serif;\n",
       "}\n",
       "h3{\n",
       "font-family: 'Fenix', serif;\n",
       "margin-top:12px;\n",
       "margin-bottom: 3px;\n",
       "}\n",
       "h4{\n",
       "font-family: 'Fenix', serif;\n",
       "}\n",
       "h5 {\n",
       "font-family: 'Alegreya Sans', sans-serif;\n",
       "}\t\n",
       "div.text_cell_render{\n",
       "font-family: 'Alegreya Sans',Computer Modern, \"Helvetica Neue\", Arial, Helvetica, Geneva, sans-serif;\n",
       "line-height: 135%;\n",
       "font-size: 120%;\n",
       "width:600px;\n",
       "margin-left:auto;\n",
       "margin-right:auto;\n",
       "}\n",
       ".CodeMirror{\n",
       "font-family: \"Source Code Pro\";\n",
       "font-size: 90%;\n",
       "}\n",
       "/* .prompt{\n",
       "display: None;\n",
       "}*/\n",
       ".text_cell_render h1 {\n",
       "font-weight: 200;\n",
       "font-size: 50pt;\n",
       "line-height: 100%;\n",
       "color:#CD2305;\n",
       "margin-bottom: 0.5em;\n",
       "margin-top: 0.5em;\n",
       "display: block;\n",
       "}\t\n",
       ".text_cell_render h5 {\n",
       "font-weight: 300;\n",
       "font-size: 16pt;\n",
       "color: #CD2305;\n",
       "font-style: italic;\n",
       "margin-bottom: .5em;\n",
       "margin-top: 0.5em;\n",
       "display: block;\n",
       "}\n",
       ".warning{\n",
       "color: rgb( 240, 20, 20 )\n",
       "}\n",
       "</style>\n",
       "<script>\n",
       "MathJax.Hub.Config({\n",
       "TeX: {\n",
       "extensions: [\"AMSmath.js\"]\n",
       "},\n",
       "tex2jax: {\n",
       "inlineMath: [ ['$','$'], [\"\\\\(\",\"\\\\)\"] ],\n",
       "displayMath: [ ['$$','$$'], [\"\\\\[\",\"\\\\]\"] ]\n",
       "},\n",
       "displayAlign: 'center', // Change this to 'center' to center equations.\n",
       "\"HTML-CSS\": {\n",
       "styles: {'.MathJax_Display': {\"margin\": 4}}\n",
       "}\n",
       "});\n",
       "</script>\n",
       "\n",
       "\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from IPython.core.display import HTML\n",
    "def css_styling():\n",
    "    styles = open('./styles/custom_barba.css', 'r').read()\n",
    "    return HTML(styles)\n",
    "css_styling()"
   ]
  },
  {
   "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.6.7"
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}