{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Power Density - All About\n",
    "\n",
    "This document gives some of the details about power density calculations.\n",
    "\n",
    "At the moment the power density calculations require points in space and the surface normal vector at each point.  There are some utilities to easily create rectangular surfaces (made up of such point and normal vectors) and some parametric surfaces.  It is also possible to calculate the power density at arbitrary points in space.  These arbitrary points in space are difficult to connect as a surface for visualization by most plotting software (typically they must be on some mesh), nevertheless the calculation is possible.\n",
    "\n",
    "For the builtin surface types there are simple ways to rotate and translate them in space, typically by giving the power density calculation function the arguments 'rotations=[rx, ry, rz]' and 'translation=[tx, ty, tz]'.  Rotations are done before translation and in the order of rx (rotation about x-axis), ry, rz ad are given in radian.  The translation is done after and given in meters.\n",
    "\n",
    "Any of these can be run in multi-threaded, GPU, or MPI mode.  Results from running on separate nodes on grid/cloud computing can be combined.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# matplotlib plots inline\n",
    "%matplotlib inline\n",
    "\n",
    "# Import the OSCARS SR module\n",
    "import oscars.sr\n",
    "\n",
    "# Import OSCARS plots (matplotlib)\n",
    "from oscars.plots_mpl import *\n",
    "\n",
    "# Import OSCARS 3D tools (matplotlib)\n",
    "from oscars.plots3d_mpl import *\n",
    "\n",
    "# Import OSCARS parametric surfaces\n",
    "from oscars.parametric_surfaces import *"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Create a new OSCARS object.  Default to 8 threads and always use the GPU if available\n",
    "osr = oscars.sr.sr(nthreads=8, gpu=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# For these examples we will make use of a simple undulator field\n",
    "osr.add_bfield_undulator(bfield=[0, 1, 0], period=[0, 0, 0.042], nperiods=31)\n",
    "\n",
    "# Plot the field\n",
    "plot_bfield(osr)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Beam\n",
    "\n",
    "Add a basic beam somewhat like NSLS2.  Filament beam for simple studies."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Add a basic electron beam with zero emittance\n",
    "osr.set_particle_beam(energy_GeV=3, x0=[0, 0, -1], current=0.500)\n",
    "\n",
    "# You MUST set the start and stop time for the calculation\n",
    "osr.set_ctstartstop(0, 2)\n",
    "\n",
    "# Plot trajectory\n",
    "plot_trajectory_position(osr.calculate_trajectory())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Power density - Rectangular planes\n",
    "\n",
    "Calculate the power density on simple rectangular planes.  The 'plane' argument is the initial plane that the surface is in.  From there it can be rotated and translated.  The order of 'XY' or any combination matters since the normal vector is given by X x Y (in the order written in the input).\n",
    "\n",
    "By default the coordinates in the returned power density are relative to the center of the rectangular grid (not absolute space).  If you want the coordinates in absolute space you can specify dim=3, but I give the warning that this makes plotting rotated surfaces difficult.\n",
    "\n",
    "Note: This is the recommended method for calculating power densities on rectangular surfaces.  One may also do the same with parametric surfaces (shown later in this tutorial)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Calculate power density 30 [m] downstream in the XY plane.\n",
    "power_density = osr.calculate_power_density_rectangle(\n",
    "    plane='XY',\n",
    "    width=[0.030, 0.030],\n",
    "    npoints=[51, 51],\n",
    "    translation=[0, 0, 30]\n",
    ")\n",
    "plot_power_density(power_density)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# We can easily rotate the above, here about the z axis.\n",
    "power_density = osr.calculate_power_density_rectangle(\n",
    "    plane='XY',\n",
    "    width=[0.030, 0.030],\n",
    "    npoints=[51, 51],\n",
    "    rotations=[0, 0, osr.pi()/4],\n",
    "    translation=[0, 0, 30]\n",
    ")\n",
    "\n",
    "plot_power_density(power_density)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Calculate the power density on a surface in the XZ plane.\n",
    "# This can be thought of as the power density on the upper beampipe inner surface\n",
    "power_density = osr.calculate_power_density_rectangle(\n",
    "    plane='XZ',\n",
    "    width=[0.020, 1.000],\n",
    "    npoints=[51, 51],\n",
    "    translation=[0, 0.004, 2]\n",
    ")\n",
    "\n",
    "plot_power_density(power_density)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# We'll now take the above plane and tilt is slightly as if it were a tapered beampipe\n",
    "power_density = osr.calculate_power_density_rectangle(\n",
    "    plane='XZ',\n",
    "    width=[0.020, 1.000],\n",
    "    npoints=[51, 51],\n",
    "    rotations=[0.002, 0, 0],\n",
    "    translation=[0, 0.004, 2]\n",
    ")\n",
    "\n",
    "plot_power_density(power_density)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Power Density - Parametric Surfaces\n",
    "\n",
    "We will now calculate the power density on a few parametric surfaces.  A parametric surface in OSCARS is represented by a class.  There is another tutorial which explains how to build your own.\n",
    "\n",
    "We begin with a simple rectangular surface, then explore some others.  These surfaces come from the oscars.parametric_surfaces module."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# First create the surface of interest\n",
    "rectangle = PSRectangle(L=0.030, W=0.030, nu=51, nv=51)\n",
    "\n",
    "# Run calculation and plotting\n",
    "pd = power_density_3d(osr, rectangle, translation=[0, 0, 30])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# It is easy to rotate the above and see it in 3D\n",
    "pd = power_density_3d(osr, rectangle, rotations=[0, osr.pi()/4, 0], translation=[0, 0, 30])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# First create the surface of interest\n",
    "sphere = PSSphere(R=0.020, nu=51, nv=51)\n",
    "\n",
    "# Run calculation and plotting\n",
    "pd = power_density_3d(osr, sphere, translation=[0, 0, 30])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# First create the surface of interest\n",
    "cylinder = PSCylinder(R=0.020, L=0.010, nu=51, nv=51)\n",
    "\n",
    "# Run calculation and plotting\n",
    "pd = power_density_3d(osr, cylinder, rotations=[osr.pi()/2, 0, 0], translation=[0, 0, 30])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Let's do a cylinder that is around the photon beam, downstream a bit\n",
    "cylinder = PSCylinder(R=0.005, L=1.000, nu=51, nv=51)\n",
    "\n",
    "# Run calculation and plotting.  Here we needed to invert the normal due to the way\n",
    "# it is defined in the PSCylinder class\n",
    "pd = power_density_3d(osr, cylinder, translation=[0, 0, 5], normal=-1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Multi-particle power density\n",
    "\n",
    "It is possible to run the power density calculations in multi-particle mode in the case that you have a very large beam, or multiple beams.  It is show here for completeness."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Add a basic electron beam with zero emittance\n",
    "osr.set_particle_beam(\n",
    "    energy_GeV=3,\n",
    "    x0=[0, 0, -1],\n",
    "    current=0.500,\n",
    "    sigma_energy_GeV=0.001*3,\n",
    "    beta=[1.5, 0.8],\n",
    "    emittance=[0.9e-9, 0.008e-9]\n",
    ")\n",
    "\n",
    "# You MUST set the start and stop time for the calculation\n",
    "osr.set_ctstartstop(0, 2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Calculate power density 30 [m] downstream in the XY plane.\n",
    "power_density = osr.calculate_power_density_rectangle(\n",
    "    plane='XY',\n",
    "    width=[0.030, 0.030],\n",
    "    npoints=[51, 51],\n",
    "    translation=[0, 0, 30],\n",
    "    nparticles=3\n",
    ")\n",
    "\n",
    "plot_power_density(power_density)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## On Precision\n",
    "\n",
    "The default relative precision is 0.01 (1%) and is controlled by the parameter:\n",
    "* precision=0.01  (default)\n",
    "\n",
    "You may retrieve the relative precision for all points in a calculation by including the parameter:\n",
    "* quantity='precision'\n",
    "\n",
    "Should you not reach the desired precision withing max_level you will receive a warning message.  To increase max_level you have two options:\n",
    "* max_level=25\n",
    "* max_level_extended=(some number above max_level)\n",
    "\n",
    "The maximum max_level is 25 due to typical memor restrictions (because it is faster).  The 'extended' version runs in non-memory mode which allows higher precision at the cost of CPU time.  Only in rare instances will you need this.  You can also retrieve the 'level' of convergence for all points (which will show -1 for non-converged points) with the addition of:\n",
    "* quantity='level'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Show the precision reached for each point.\n",
    "power_density = osr.calculate_power_density_rectangle(\n",
    "    plane='XY',\n",
    "    width=[0.030, 0.030],\n",
    "    npoints=[51, 51],\n",
    "    translation=[0, 0, 30],\n",
    "    quantity='precision'\n",
    ")\n",
    "plot_power_density(power_density, title='Precision')\n",
    "\n",
    "# Show the level reached for each point.\n",
    "power_density = osr.calculate_power_density_rectangle(\n",
    "    plane='XY',\n",
    "    width=[0.030, 0.030],\n",
    "    npoints=[51, 51],\n",
    "    translation=[0, 0, 30],\n",
    "    quantity='level'\n",
    ")\n",
    "plot_power_density(power_density, title='Level')"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "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.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}