{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this post, we will develop a simple ray tracing model, which can be a very useful tool in many domains. For instance in acoustics, optics, fluids dynamics... Why so many different domains? Well, this is related to the fact that under certain hypotheses (high frequency and paraxial approximations) waves (light, sound, water...) behave like rays."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Modelling our ray tracer"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There's many good introductions to ray tracing, but here we will define things as we go. The way we want to trace rays is very basic. Essentially, we want to be able to:\n",
    "\n",
    "- draw any geometry \n",
    "- account for different wave speeds in different matters\n",
    "\n",
    "To do this, we will model the geometry as a sequence of oriented segments, with wave speed properties on each side of the segment. Any type of geometry will thus be represented as a list of segments. Also, we will only work in 2d.\n",
    "\n",
    "A single segment might be modelled as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "class Segment:\n",
    "    \"A class representing segments.\"\n",
    "    def __init__(self, A, B, normal_vector, celerity_front, celerity_back):\n",
    "        \"Inits a segment from A to B, with a normal vector and front and back celerities.\"\n",
    "        self.A = A\n",
    "        self.B = B\n",
    "        self.normal_vector = normal_vector\n",
    "        self.celerity_front = celerity_front\n",
    "        self.celerity_back = celerity_back"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A sphere would thus be represented by a sequence of segments. Any geometry could thus just be a list of segments.\n",
    "\n",
    "So how does one trace a ray? There's three things here:\n",
    "\n",
    "- intersection of a ray with a segment\n",
    "- refraction\n",
    "- reflection\n",
    "\n",
    "Let's tackle these. But first, we need to define a ray."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Defining a ray"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A ray needs an origin and an end."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "class Ray:\n",
    "    \"A class representing a ray.\"\n",
    "    def __init__(self, origin, end):\n",
    "        self.origin = origin\n",
    "        self.end = end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we also need sources of rays, which are defined by a point and a direction, which is an unit vector."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "class Source:\n",
    "    \"A ray source.\"\n",
    "    def __init__(self, origin, direction):\n",
    "        self.origin = origin\n",
    "        self.direction = direction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we also need points:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from collections import namedtuple\n",
    "\n",
    "Point = namedtuple('Point', 'x, y')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Point(x=1, y=2)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Point(1, 2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And vectors:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "Vector = namedtuple('Vector', 'u_x, u_y')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Vector(u_x=3.0, u_y=2.4)"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Vector(3., 2.4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Intersection of a source with a segment to create a ray"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Once we have a source, we can trace it until a segment with which it intersects. The intersection point allows us to return a ray."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def trace_ray(source, geometric_object):\n",
    "    \"Traces a ray from a source to a geometric object.\"\n",
    "    intersected_segment = min(geometric_object.segments, key=lambda seg: tof(source, seg))\n",
    "    intersection_point = intersection(source, intersected_segment)\n",
    "    return Ray(source.origin, intersection_point)\n",
    "\n",
    "# TODO: tof function, intersection function"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's now turn to the intersection function. Our rays are always straight lines. So it is easy to come up with a way for intersecting a line."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "ename": "SyntaxError",
     "evalue": "unexpected EOF while parsing (<ipython-input-9-30eed575786f>, line 6)",
     "output_type": "error",
     "traceback": [
      "\u001b[0;36m  File \u001b[0;32m\"<ipython-input-9-30eed575786f>\"\u001b[0;36m, line \u001b[0;32m6\u001b[0m\n\u001b[0;31m    else:\u001b[0m\n\u001b[0m         ^\u001b[0m\n\u001b[0;31mSyntaxError\u001b[0m\u001b[0;31m:\u001b[0m unexpected EOF while parsing\n"
     ]
    }
   ],
   "source": [
    "def intersection(source, segment):\n",
    "    \"Computes the intersection point between a source and a segment.\"\n",
    "    # they're parallel: no intersection or infinity\n",
    "    if scalar(source.direction, segment.normal_vector) == 0.:\n",
    "        return None\n",
    "    else:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "A = Point(0, 0)\n",
    "B = Point(5, 0)\n",
    "S = Point(1, -4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "plt.style.use('bmh')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def draw_points(points, labels=None):\n",
    "    \"Draws a list of points with labels.\"\n",
    "    if labels is None:\n",
    "        labels = ['Point {}'.format(ind) for ind, val in enumerate(points)]\n",
    "    for p, label in zip(points, labels):\n",
    "        plt.plot(p.x, p.y, 'o', label=label)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-5, 1)"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAD5CAYAAAA6JL6mAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAGY5JREFUeJzt3X1sVYd5x/Hfgy82DPMSSGjApCQVJR1iBSqXNDLNmrg0\nWdI2CVetUimblqZqWq1TsxWlY0hb12hL0qohi7Yum9KyTX2JljiEvDECLBZQURIwhDokdaGrU+yB\nk9r3OqZcG9vP/vC93rF9je/xeXzPedjvI0WNfa7v/QrTx5dzHk5EVUFERH5NizuAiIii4SAnInKO\ng5yIyDkOciIi5zjIiYic4yAnInLOZJCLyPdFpENEmi2ej4iISmf1jvxfAdxk9FxERBSCySBX1b0A\nOi2ei4iIwuE5ciIi51LleJEdO3bo6dOnhz++5JJLMG/ePFRWVkJE0N/fj4GBgTFfV+7jqgoRie31\nSz1e6ExqX0FFRQVSqVRi+wrH+/r6UOxWFUnp4+9P++OqiqqqqsT2FY6fO3funfr6+svGHBylLIO8\nuroat912WzleKpJMJoN58+bFnTEhdtpipy0PnR4aAaCpqam1lMeV5dRKKlWWnxeRefjGAuy0xk5b\nHjo9NIZhtX74YwAHAFwtIqdE5O7gcS93WOzv7487oSTstMVOWx46PTSGYbW18jlVXaSq01V1iap+\nL3i8r6/P4mWmXFtbW9wJJWGnLXba8tDpoTEMH+c8iIgmQVXR09Mz5qxAZWUluru7Y6oaS0RQXV09\n4mJ2GBzkRHTR6unpQVVVFSorK0d8vqqqClVVVTFVjdXX14eenh7Mnj17Ul/PPXIiumip6pghnkSV\nlZWRriVykBMROVeWQV5RUVGOl4ls7ty5cSeUhJ222GnLQ2e5Z9ILL7yA+fPno6WlZUqen3vkAV52\nS9lpi522PHSON5P2nOjEnU8048bHj+DOJ5qx54TNLaQaGhrwkY98BA0NDSbPN1pZBjn3yG2x0xY7\nbXnoLDaT9pzoxCP73kJHz3kogI6e83hk31uRh3lPTw8OHjyIRx99FNu2bYv0XOMpyyDnHrktdtpi\npy0PncVm0tZD7egdGDngewcUWw+1R3qtHTt24IYbbsCyZctwySWX4OjRo5Gerxhe7CQiAvB2z/lQ\nny9VQ0MDNmzYAADYsGHDlJxe8XHymohoil1WPR0dRYb2ZdXTJ/2cXV1d2LdvH44fPw4RwcDAAEQE\n3/zmNyf9l3+K4TtyIiIAd9UuRlXFyOFaVSG4q3bxpJ9z+/bt+OxnP4tjx47htddeQ3NzM5YuXYoD\nBw5EzR2Bg5yICED9svm496PvxcLq6RAAC6un496Pvhf1y+ZP+jmffvpp3HLLLSM+96lPfcr89EpZ\nTq1wj9wWO22x05aHzvFmUv2y+ZEG92jPPvvsmM/dc889Zs9fwD3yAA/7rwA7rbHTlodOLzOpVNwj\nD/Cw/wqw0xo7bXno9DKTSsU98gAP+68AO62x05aHTi8zqVS82ElE5BwHORGRcxzkRETOXVyXbomI\nEubSSy/FihUroKqoqKjAQw89hGuuucb0NbhHHuBh/xVgpzV22vLQOd5Mam/YiZYHHkOurQMzahZi\n+aYvYXH6xkivNXPmTOzduxcAsGfPHtx///14/vnnIz3naGUZ5F52Nj3svwLstMZOWx46i82k9oad\naN74IAbP9QIAcqfOoHnjgwAQeZgXvPvuu1Py61OWCetlZ7O/v9/FDx122mKnLQ+dqjrmplUtDzw2\nPMQLBs/1ouWBxyIN8nPnzuG6665Db28vTp8+je3bt0/6ucbDPfIAD/uvADutsdOWh85iMynX1lH0\nseN9vlSFUysHDx7Ek08+iS9/+cvmb265tUJEBGBGzcJQn5+MtWvXorOzE++8847ZcwIc5EREAIDl\nm76EaTOrRnxu2swqLN/0JbPXaGlpwcDAAObPt7sxF8D1QyIiAP93QdN6a6VwjhwYOjf/3e9+13yT\nz2SQi8hNAP4eQAWAx1X1QYvnJSIqp8XpG802VAqsT6MUE3mQi0gFgH8EsB7AKQCvisizqnq88Jik\n75HvOdGJrYfa8XbPeVxW/S7uql1sek9iax72dAF2Wkt6Z3AH+7+N3s1OlaTPpLAs3pGvBXBCVX8J\nACLyBIBbAQwP8iSvIu050YlH9r01/F/P7ug5j0f2vQUAiR3mHvZ0AXZaS3JnOXawLSV5Jk2GxcXO\nGgC/Dnx8Kv+5YUneI996qH14iBf0Dii2HmqPqWhiHu73DLDTWpI7L7SDnURJnkmTUZYfS2fOnMEX\nvvCF4Y/T6TTS6TRqamqQSqWQyWSQzWbHfF05jr9d5L+aDQy9My/8xYY4+4odz2azmDt3biJ+/S50\nfHBwEFdddVVi+wrHT548iRkzZiS2r3C88H1PYl+u7cyY48DQDnacfcDQD8CBgYERxwYGBjBz5kyI\nSNHjAFBZWVnW42fPnkVXV9eY/lJI1J9MInItgG+o6o35jzcBgKo+UHhMY2Ojrlq1KtLrTJU7n2hG\nR5FhvrB6On5wx8oYiibW2tqKpUuXxp0xIXbaSnJnY+3tyJ0aO8xnLHkPPnZoWwxFQ7q7uzFnzpwx\nn+/t7UVVVVWRr4hPsdampqbD9fX1tRN9rcWplVcBvF9ErhKRSgB3ABj7XxxNqLtqF6OqYuRf1a2q\nENxVuzimIiJ/yrGDTeOLfGpFVftF5CsAdmJo/fD7qvp65LIyKVzQ3HqoHR0957Gwenrit1aIkmaq\ndrAvFt/5znfw1FNPoaKiAtOmTcPDDz+M2toJ32iXzOQcuaq+COBFi+eKQ/2y+ahfNj/Rf3QlSrrC\nDrbn/x8dP9KO/bta0J3JYc68GVi3fjlWrIn2p/NXXnkFO3fuRGNjI6qqqvCb3/zG/P5TvB95QNL3\ndAvYaYudtjx0FptJx4+046VnmtF/fhAA0J3J4aVnmgEg0jA/c+YMFixYMHxOfsGCBZN+rvGU5V4r\nXnY2k7ynG8ROW+y05aGz2Ezav6tleIgX9J8fxP5dLZFe6/rrr0dbWxs+/OEPY+PGjfjJT34S6fmK\nKcsg97KzmeQ93SB22mKnLQ+dxWZSdyZX9LHjfb5U1dXVePnll7FlyxYsWLAAd999N370ox9Fes7R\neD/yAA/3UQbYaY2dtjx0FptJc+aN/bsEF/p8GBUVFVi3bh02bdqEhx56CM8991zk5wzibWyJiACs\nW78cqekjR2Jq+jSsW7880vP+4he/wMmTJ4c/bm5uxhVXXBHpOUfzcfKaiGiKFS5oWm+tnD17Fl//\n+teRzWaRSqXwvve9D1u2bLFIHsZBTkSUt2LN4siDe7TVq1dj586dps85Gk+tEBE5V5ZBzj1yW+y0\nxU5bHjq9zKRScY88wMP+K8BOa+y05aHTy0wqFffIAzzsvwLstMZOWx46vcykUnGPPMDD/ivATmvs\ntJWkThEpOn+SNpP6+vogIhM/cBwX158viIgCqqur0dPTg1xu5N/OPHv2LGbNmhVT1Vgigurq6kl/\nPQc5EV20RASzZ88e8/muri4sWrQohqKpwfVDIiLnOMiJiJzjHnmAh/1XgJ3W2GnLQ6eHxjC4Rx7g\nYf8VYKc1dtry0OmhMQzukQd42H8F2GmNnbY8dHpoDIN75AFJ2n+9EHbaYqctD50eGsPgxU4iIuc4\nyImInOMgJyJyjoOciMg57pEHeNktZactdtry0OmhMQzukQd42S1lpy122vLQ6aExDO6RB3jZLWWn\nLXba8tDpoTEM7pEHeNktZactdtry0OmhMYxIg1xEPiMir4vIoIjUWkUREVHpor4jbwawAcBegxYi\nIpqESFchVfUNAJH+E0VERBRNWdZJurq6UFdXN/xxOp1GOp1GTU0NUqkUMpkMstnsmK8r9/FsNovW\n1tbYXr/U44XOpPYVDA4OAkBi+wrHc7nciO970vr4+9P+eDabRX9/f2L7CsdLJRNtlIjIbgCXFzm0\nWVW35x/TCGCjqh4q9hz79u3TlStXlhwVl0wm42ItiZ222GnLQ6eHRgBoamo6XF9fP+H1xwnfkavq\nx6PGcI/cFjttsdOWh04PjWFwjzzAy24pO22x05aHTg+NYURdP7xdRE4BuBbACyKys9jjuEdui522\n2GnLQ6eHxjCibq1sA7DNqIWIiCaBdz8kInKOg5yIyDkOciIi53g/8gAv9yhmpy122vLQ6aExDN6P\nPMDLbik7bbHTlodOD41hcI88wMtuKTttsdOWh04PjWHwfuQBXnZL2WmLnbY8dHpoDIMXO4mInOMg\nJyJyjoOciMg5DnIiIue4Rx7gZbeUnbbYactDp4fGMLhHHuBlt5Sdtthpy0Onh8YwuEce4GW3lJ22\n2GnLQ6eHxjC4Rx7gZbeUnbbYactDp4fGMHixk4jIOQ5yIiLnOMiJiJzjICcico575AFedkvZaYud\ntjx0emgMg3vkAV52S9lpi522PHR6aAyDe+QBXnZL2WmLnbY8dHpoDIN75AFedkvZaYudtjx0emgM\ngxc7iYic4yAnInKOg5yIyDkOciIi5yINchH5toi8KSLHRGSbiBTd6eEeuS122mKnLQ+dHhrDiPqO\nfBeAlar6QQAtADYVexD3yG2x0xY7bXno9NAYRqRBrqovqWphIfOnAJaM87goL1M2XnZL2WmLnbY8\ndHpoDMPyHPnnAewodoB75LbYaYudtjx0emgMY8JzHiKyG8DlRQ5tVtXt+cdsBtAP4IfFnqOrqwt1\ndXXDH6fTaaTTadTU1CCVSiGTySCbzY75unIfz2azaG1tje31Sz1e6ExqX8Hg4CAAJLavcDyXy434\nvietj78/7Y9ns1n09/cntq9wvFQS9bSHiPwxgHsA1Kvqb4s9prGxUVetWhXpdcqhtbUVS5cujTtj\nQuy0xU5bHjo9NAJAU1PT4fr6+tqJHhfpKqSI3ATgPgC/P94QJyKiqRX1HPk/AJgNYJeIHBWRxwya\niIgohEjvyFV1WSmP4x65LXbaYqctD50eGsPg/cgDvOyWstMWO2156PTQGAbvRx7gZbeUnbbYactD\np4fGMHg/8gAvu6XstMVOWx46PTSGwZtmERE5x0FOROQcBzkRkXMc5EREzpVlkHOP3BY7bbHTlodO\nD41hcI88wMtuKTttsdOWh04PjWFwjzzAy24pO22x05aHTg+NYXCPPMDLbik7bbHTlodOD41h8GIn\nEZFzHORERM5xkBMROcdBTkTkHPfIA7zslrLTFjtteej00BgG98gDvOyWstMWO2156PTQGAb3yAO8\n7Jay0xY7bXno9NAYBvfIA7zslrLTFjtteej00BgGL3YSETnHQU5E5BwHORGRcxzkRETOcY88wMtu\nKTttsdOWh04PjWFwjzzAy24pO22x05aHTg+NYXCPPMDLbik7bbHTlodOD41hcI88wMtuKTttsdOW\nh04PjWFEGuQicr+IHBORoyLykogstgojIqLSRH1H/m1V/aCqrgbwPIC/MmgiIqIQIg1yVe0OfDgL\ngI+T4UREF5HI6yQi8rcA/ghAFsD1kYuIiCiUCQe5iOwGcHmRQ5tVdbuqbgawWUQ2AfgKgL8e/cBM\nJoO6urrhj9PpNNLpNGpqapBKpZDJZJDNZse8QLmP53I5tLa2xvb6pR4vdCa1r6C6uhoAEttXOF5R\nUTHi+560Pv7+tD+ey+XQ39+f2L7C8VKJ1WqgiLwXwIuqunL0sQMHDugHPvABk9chIvr/oqmp6XB9\nfX3tRI+LurXy/sCHtwJ4s9jjuEdui5222GnLQ6eHxjCibq08KCLNInIMwCcAfLXYg7hHboudtthp\ny0Onh8YwIl3sVNW0VQgREU0O735IROQcBzkRkXMc5EREzvF+5AFe7lHMTlvstOWh00NjGLwfeYCX\nexSz0xY7bXno9NAYBu9HHuBlt5Sdtthpy0Onh8YweD/yAC+7pey0xU5bHjo9NIbBi51ERM5xkBMR\nOcdBTkTkHAc5EZFz3CMP8LJbyk5b7LTlodNDYxjcIw/wslvKTlvstOWh00NjGNwjD/CyW8pOW+y0\n5aHTQ2MY3CMP8LJbyk5b7LTlodNDYxi82ElE5BwHORGRcxzkRETOcZATETnHPfIAL7ul7LTFTlse\nOj00hsE98gAvu6XstMVOWx46PTSGwT3yAC+7pey0xU5bHjo9NIbBPfIAL7ul7LTFTlseOj00hsGL\nnUREznGQExE5x0FOROQcBzkRkXMmg1xEviYiKiKXFjvOPXJb7LTFTlseOj00hhF5wVtErgDwCQBv\njfsi3CM3cfxIO/bvakF3Joc582Zg3frlWLFmcdxZ40r6r2cBO2156PTQGIbFO/ItAO4DMO6yOPfI\nozt+pB0vPdOM7kwOANCdyeGlZ5px/Eh7zGXjS/KvZxA7bXno9NAYRqRBLiK3AmhT1dcu9DjukUe3\nf1cL+s8Pjvhc//lB7N/VElPRxJL86xnETlseOj00hjHhOQ8R2Q3g8iKHNgP4SwydVrmgrq4u1NXV\nDX+cTqeRTqdRU1ODVCqFTCaDbDY75uvKfTybzaK1tTW217/Q8cI78dG6M7lE9BU7Pjg49IMnqX2F\n47lcbsT3PWl9Hn5/Bo8XOpPaV2js7+9PbF/heKlksqc9ROT3AOwB8Nv8p5YAaAewVlVPBx/b2Nio\nq1atmtTrlFNrayuWLl0ad0ZR//KtxqLDfM68GfjifR8rf1AJkvzrGcROWx46PTQCQFNT0+H6+vra\niR436VMrqvozVV2oqleq6pUATgH40OghTjbWrV+O1PSR367U9GlYt355TEVElBQ+1kloeDvF09YK\nEZWH2SDPvysvinvkNlasWYwVaxYjk8m4WJ9K+q9nATtteej00BgG70ce4GE4Auy0xk5bHjo9NIbB\n+5EHeNktZactdtry0OmhMQzejzzAy24pO22x05aHTg+NYfCmWUREznGQExE5x0FOROQcBzkRkXNl\nGeTcI7fFTlvstOWh00NjGNwjD/CyW8pOW+y05aHTQ2MY3CMP8LJbyk5b7LTlodNDYxjcIw/wslvK\nTlvstOWh00NjGLzYSUTkHAc5EZFzHORERM5xkBMROcc98gAvu6XstMVOWx46PTSGwT3yAC+7pey0\nxU5bHjo9NIbBPfIAL7ul7LTFTlseOj00hsE98gAvu6XstMVOWx46PTSGwYudRETOcZATETnHQU5E\n5BwHORGRc9wjD/CyW8pOW+y05aHTQ2MY3CMP8LJbyk5b7LTlodNDYxjcIw/wslvKTlvstOWh00Nj\nGNwjD/CyW8pOW+y05aHTQ2MYvNhJRORcpEEuIt8QkTYROZr/5+Zij8tkMlFepmwaGhriTigJO22x\n05aHTg+NANDZ2XlpKY+zeEe+RVVX5/95sdgDurq6DF5m6nn55rLTFjtteej00AgA3d3dl5XyOJ5a\nISJyTqJslIjINwDcBSAL4BCAr6nqmLffzz33XK6jo2Og8PGcOXPenj9//juTfuEp0tnZeWkSu0Zj\npy122vLQ6aERAHp7e6+++eabZ0/0uAkHuYjsBnB5kUObAfwUwDsAFMD9ABap6ufD5xIR0WRFekc+\n4olErgTwvKquNHlCIiIqSdStlUWBD28H0Bwth4iIwop6sfNbIvIzETkG4HoAfzbeA0XkMyLyuogM\nikhtxNc1JyI3icjPReSEiPxF3D3FiMj3RaRDRBL9A1NErhCRl0XkeP57/tW4m0YTkRki8oqIvJZv\n/Ju4my5ERCpE5IiIPB93y3hE5Ff5eXBURA7F3TMeEZknIk+JyJsi8oaIXBt302gicnVgrfuoiHSL\nyL3jPr5cf31eRH4XwCCAfwawUVUT840WkQoALQDWAzgF4FUAn1PV47GGjSIi1wHoAfDvST6Flf+T\n2iJVbRKR2QAOA7gtSb+eIiIAZqlqj4hMB7AfwFdV9acxpxUlIn8OoBbAHFX9ZNw9xYjIrwDUqmqi\nLyKKyL8B2Keqj4tIJYDfUdXE/mWX/HxqA3CNqrYWe0zZ1g9V9Q1V/Xm5Xi+ktQBOqOovVbUPwBMA\nbo25aQxV3QugM+6Oiajq/6hqU/7f3wXwBoCaeKtG0iE9+Q+n5/9J5E2BRGQJgFsAPB53i3ciMhfA\ndQC+BwCq2pfkIZ5XD+DkeEMc4B55QQ2AXwc+PoWEDR6v8hfB1wA4GG/JWPnTFUcBdADYpaqJa8x7\nBMB9GPoTbZIpgN0iclhEvhh3zDiuAvA2gK35U1WPi8isuKMmcAeAH1/oAaaDXER2i0hzkX8S9+6W\npp6IVANoAHCvqnbH3TOaqg6o6moASwCsFZHEna4SkU8C6FDVw3G3lGBd/tfzDwD8Sf5UYNKkAHwI\nwD+p6hoAZwEk8poYAORP/XwawJMXepzpjcJV9eOWz1dGbQCuCHy8JP85mqT8eecGAD9U1afj7rkQ\nVc2IyMsAbkLyNq/qAHw6fx+jGQDmiMgPVPXOmLvGUNW2/P92iMg2DJ2y3Btv1RinAJwK/OnrKSR4\nkGPoh2KTqp650IN4amXIqwDeLyJX5X8C3gHg2Zib3MpfSPwegDdU9eG4e4oRkctEZF7+32di6EL3\nm/FWjaWqm1R1iapeiaHfl/+VxCEuIrPyF7aRP1XxCSTvhyJU9TSAX4vI1flP1QNIzEX4Ij6HCU6r\nAGUc5CJyu4icAnAtgBdEZGe5XnsiqtoP4CsAdmLowtx/qOrr8VaNJSI/BnAAwNUickpE7o67aRx1\nAP4QwA0T3RkzRosAvJxfnX0VQ+fIE7va58B7AOwXkdcAvALgBVX9z5ibxvOnAH6Y/96vBvB3MfcU\nlf+BuB7AhH+iLdv6IRERTQ2eWiEico6DnIjIOQ5yIiLnOMiJiJzjICcico6DnIjIOQ5yIiLnOMiJ\niJz7X6+cvkn+TWd1AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x8dbb860>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "draw_points([A, B, S], ['A', 'B', 'S'])\n",
    "plt.legend()\n",
    "plt.xlim(-1, 7)\n",
    "plt.ylim(-5, 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "ename": "SyntaxError",
     "evalue": "invalid syntax (<ipython-input-14-318c53a6be6f>, line 1)",
     "output_type": "error",
     "traceback": [
      "\u001b[0;36m  File \u001b[0;32m\"<ipython-input-14-318c53a6be6f>\"\u001b[0;36m, line \u001b[0;32m1\u001b[0m\n\u001b[0;31m    AB =\u001b[0m\n\u001b[0m         ^\u001b[0m\n\u001b[0;31mSyntaxError\u001b[0m\u001b[0;31m:\u001b[0m invalid syntax\n"
     ]
    }
   ],
   "source": [
    "AB = "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Information about line intersections can be found here:\n",
    "\n",
    "<http://geomalgorithms.com/a05-_intersect-1.html>\n",
    "\n",
    "I thing this is the way to go."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Creating a nice visualization "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Drawing rays that intersect a straight wall."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "s = Point(0, 0)\n",
    "east = Vector(1, 0)\n",
    "source = Source(s, direction=east, aperture=(-30, 30), n=11)\n",
    "top = Point(1, 1)\n",
    "bottom = Point(1, -1)\n",
    "wall = Segment(top, bottom)\n",
    "rays = Intersect(source, wall)\n",
    "plot(rays)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Drawing rays that get focused by spherical shell."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "s = Point(0, 0)\n",
    "east = Vector(1, 0)\n",
    "source = Source(s, direction=east, aperture=(-30, 30), n=11)\n",
    "top = Point(1, 1)\n",
    "bottom = Point(1, -1)\n",
    "wall = Segment(top, bottom)\n",
    "rays = Intersect(source, wall)\n",
    "plot(rays)"
   ]
  }
 ],
 "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.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}