{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# triangles.ipynb\n",
    "\n",
    "Discussion:\n",
    "\n",
    "  This Jupyter notebook investigates how triangles can be defined and analyzed.\n",
    "  \n",
    "Licensing:\n",
    "\n",
    "  This code is distributed under the GNU LGPL license.\n",
    "\n",
    "Modified:\n",
    "\n",
    "  21 October 2016\n",
    "  \n",
    "Author:\n",
    "\n",
    "  John Burkardt, Lukas Bystricky\n",
    "  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Using matplotlib backend: agg\n"
     ]
    }
   ],
   "source": [
    "# Import libraries and set the plot option.\n",
    "#\n",
    "#  Lukas uses \"matplotlib inline\" but I have to use \"matplotlib\" and I still don't get any graphics.\n",
    "#\n",
    "%matplotlib\n",
    "%config InlineBackend.figure_format = 'png'\n",
    "\n",
    "import math\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Triangles #\n",
    "\n",
    "The **triangle** is the fundamental interesting object in 2D.\n",
    "\n",
    "We will look at:\n",
    "* how a triangle can be defined and plotted.\n",
    "* how properties of a triangle can be computed.\n",
    "* how triangles can be sampled regularly or at random.\n",
    "* how to determine if a point is inside a triangle.\n",
    "* a special coordinate system for a triangle.\n",
    "* estimating an integral over a triangle.\n",
    "\n",
    "Some of these ideas will lead us to understand how to work\n",
    "with polygons (by breaking them into triangles) and to\n",
    "analyze a collection of points by organizing them into a triangulation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0 0]\n",
      " [5 0]\n",
      " [0 3]]\n",
      "[[0 0]\n",
      " [1 0]\n",
      " [0 1]]\n"
     ]
    }
   ],
   "source": [
    "#  Defining a triangle:\n",
    "#\n",
    "#  1) A Pythagorean 3,4,5 triangle can be defined by the points\n",
    "#    (0,0), (5,0), (0,3).\n",
    "#  We will define Python variable T1 to be a triangle by simply\n",
    "#  storing this data as an array of 3 rows and 2 columns.\n",
    "#  Each row contains the coordinates of a point.  \n",
    "#  Column 0 is the X, and column 1 the Y coordinate.\n",
    "#\n",
    "#  Print t1 and make sure you believe what you see.\n",
    "#\n",
    "t1 = np.array ( [ [ 0, 0 ], [ 5, 0 ], [0, 3] ] )\n",
    "print ( t1 )\n",
    "#\n",
    "#  2) Suppose we have already defined A, B, and C, which contain point coordinates.\n",
    "#  Then we could build a triangle by creating an array of these points.\n",
    "#  The reference triangle has coordinates (0,0), (1,0), (0,1).\n",
    "#\n",
    "#  Print t2 and make sure that the information looks right to you.\n",
    "#\n",
    "a = np.array ( [ 0, 0 ] )\n",
    "b = np.array ( [ 1, 0 ] )\n",
    "c = np.array ( [ 0, 1 ] )\n",
    "t2 = np.array ( [ a, b, c ] )\n",
    "print ( t2 )\n",
    "#\n",
    "#  3) Another way to build a triangle is to assume that there is already a list \"xy\"\n",
    "#  of coordinates of many points, in which case we can define the triangle simple\n",
    "#  by selecting the indices of those points.\n",
    "#\n",
    "#  We can either copy those indices into our new object t3, or we can just\n",
    "#  store the indices into an object we'll call it3, in which case we must\n",
    "#  keep \"xy\" around.\n",
    "#\n",
    "#  The equilateral triangle t3 has coordinates (0,0), (1,0), (0.5,sqrt(3)/2).\n",
    "#\n",
    "#  Print xy.\n",
    "#\n",
    "xy = np.array ( [ [ 0, 0 ], [ 5, 0 ], [ 0, 3 ], [ 1, 0 ], [ 0, 1 ], [ 0.5, np.sqrt(3.0)/2.0 ] ] )\n",
    "#\n",
    "#  Print t3, which selects rows of XY.\n",
    "#\n",
    "t3 = np.array ( [ xy[0], xy[3], xy[5] ] )\n",
    "#\n",
    "#  Print xy[it3], which should produce the same results.\n",
    "#\n",
    "it3 = np.array ( [ 0, 3, 5 ] )\n",
    "#\n",
    "#  When we work with triangulations, we will actually find it more convenient\n",
    "#  to use the indexed description of a triangle, as exemplified by it3.\n",
    "#"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The counter clockwise convention ##\n",
    "\n",
    "We describe a triangle as a list of three vertices, which means there are six\n",
    "possible orderings.  However, there is a mathematical convention that triangles\n",
    "(and polygons, in general) should be described by listing the vertices in\n",
    "counter-clockwise order.  Right now, it's not clear why such a rule should be used.\n",
    "However, if we follow this rule, then it makes it easy to write down correct \n",
    "formulas for triangle angles and areas; it makes it possible to determine\n",
    "whether a point is inside or outside a triangle; it makes it possible to determine\n",
    "the outward normal direction along the edges of a triangle.  The three triangles\n",
    "we used above all followed the counter-clockwise convention, and we will continue\n",
    "to do so."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "13.8309518948\n",
      "3.41421356237\n",
      "3.0\n"
     ]
    }
   ],
   "source": [
    "## Perimeters and Edge Lengths ##\n",
    "#\n",
    "#  Let's compute the perimeter of each triangle t1, t2 and t3.\n",
    "#  This is simply the sum of the lengths of the sides.\n",
    "#  Each calculation is a little different.\n",
    "#  \n",
    "#  To compute the Euclidean norm of an edge, we want to use the \"norm\" command\n",
    "#  to measure the distance between each pair of vertices.\n",
    "#\n",
    "#  Print the values of p1, p2 and p3 and see if you find them believable.\n",
    "#\n",
    "from numpy.linalg import norm\n",
    "\n",
    "p1 = norm ( t1[1] - t1[0] ) + norm ( t1[2] - t1[1] ) + norm ( t1[0] - t1[2] )\n",
    "p2 = norm ( b - a ) + norm ( c - b ) + norm ( a - c )\n",
    "p3 = norm ( xy[3] - xy[0] ) + norm ( xy[5] - xy[3] ) + norm ( xy[0] - xy[4] )\n",
    "print ( p1 )\n",
    "print ( p2 )\n",
    "print ( p3 )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The angles of a triangle ##\n",
    "\n",
    "The first thing to remember is that most scientific calculations use radians\n",
    "for angle measure, so we expect the angles of a triangle to sum to pi, not 180.\n",
    "\n",
    "The second thing to point out is the law of cosines.  If the vertices are A, B, C\n",
    "and the corresponding angles are a, b, c, then we have:\n",
    "\n",
    "  cos ( a ) = ( ||B-A||^2 + ||A-C||^2 - ||C-B||^2 ) / ( 2 ||B-A|| || A-C|| )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3.14159265359\n"
     ]
    }
   ],
   "source": [
    "## The angles of a triangle:\n",
    "#\n",
    "#  Here is a calculation for the angle \"a\" in triangle t1:\n",
    "#\n",
    "bma = norm ( t1[1] - t1[0] )\n",
    "cmb = norm ( t1[2] - t1[1] )\n",
    "amc = norm ( t1[0] - t1[2] )\n",
    "cosa = ( bma**2 + amc**2 - cmb**2 ) / ( 2 * bma * amc )\n",
    "a1 = np.arccos ( cosa )\n",
    "#\n",
    "#  Compute the other two angles b1 and c1, and compare their sum to pi.\n",
    "#  (You can use the value np.pi)\n",
    "#"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The centroid of a triangle ##\n",
    "\n",
    "The centroid of a triangle is essentially the center of mass, assuming\n",
    "that the triangle is a physical object of uniform thickness.\n",
    "That means the centroid is always inside the triangle.\n",
    "\n",
    "The centroid can be found by averaging the coordinates of the three vertices."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 2.33333333  1.66666667]\n"
     ]
    }
   ],
   "source": [
    "## The centroid of a triangle##\n",
    "#\n",
    "#  Compute the centroid of triangle t2, call it \"cent\", and print its value here.\n",
    "#\n",
    "cent = ( a + b + c ) / 3.0\n",
    "print ( cent )\n",
    "#\n",
    "#  Then run the following commands to make sure your centroid is inside your triangle.\n",
    "#\n",
    "plt.plot ( [a[0], b[0], c[0], a[0] ], [ a[1], b[1], c[1], a[1] ], 'b-' )\n",
    "plt.plot ( cent[0], cent[1], 'ro' )\n",
    "plt.savefig ( 'triangles_fig1.png' )\n",
    "plt.close ( )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The area of a triangle##\n",
    "\n",
    "You may have learned a formula for the area of a triangle that multiplies\n",
    "the length of the base times 1/2 the height.  \n",
    "\n",
    "A simpler approach just uses the coordinates of the vertices directly, \n",
    "and has the following form:\n",
    "\n",
    "  area = 1/2 ( Ax * ( By - Cy ) + Bx * ( Cy - Ay ) + Cx * ( Ay - By ) )\n",
    "\n",
    "where the vertices are A, B, C and Ax and Ay mean the x and y coordinates of A."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "7.5\n"
     ]
    }
   ],
   "source": [
    "## The area of a triangle##\n",
    "#\n",
    "#  Compute and print the areas of triangles t1, t2 and t3.\n",
    "#\n",
    "#  Because we have defined t1, t2 and t3 in different ways, each formula for\n",
    "#  the area will look somewhat different.\n",
    "#\n",
    "#  Make a triangle t4 which is the same as triangle t2, except that the vertices\n",
    "#  are listed in clockwise order.  Compute and print the area of t4, and notice\n",
    "#  how it is different from the area of t2.\n",
    "#\n",
    "area1 = 0.5 * ( t1[0,0] * ( t1[1,1] - t1[2,1] ) + t1[1,0] * ( t1[2,1] - t1[0,1] ) + t1[2,0] * ( t1[0,1] - t1[1,1] ) )\n",
    "print ( area1 )\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## A regular grid of points in a triangle##\n",
    "\n",
    "Suppose we want to sample points in a triangle using a regular grid.\n",
    "Then the grid can be thought of as consisting of N+1 equally spaced points\n",
    "along one edge, then a layer of N points, then a layer of N-1 points,\n",
    "..., all the way to a layer of just 1 point, the vertex.  As it happens,\n",
    "we should see a stack of N+1 layers no matter which edge of the triangle\n",
    "we start at.  How can we compute such a grid of points?\n",
    "\n",
    "It turns out that there is a very interesting formula for doing this:\n",
    "\n",
    "Let A, B and C be the vertices, and G one of the grid points.\n",
    "\n",
    "  Let I run from 0 up to N\n",
    "    Let J run from 0 up to N - I\n",
    "      K = N - I - J\n",
    "      G = ( I * A + J * B + K * C ) / N"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "## A regular grid of points in a triangle##\n",
    "#\n",
    "#  Define triangle \"t5\", whose vertices are:\n",
    "#    A=(4,0), B=(3,4), C=(0,1)\n",
    "#\n",
    "#  Generate a grid that has N=10 (so 11 layers)\n",
    "#\n",
    "a = np.array ( [4,0] )\n",
    "b = np.array ( [3,4] )\n",
    "c = np.array ( [0,1] )\n",
    "\n",
    "plt.plot ( [ a[0], b[0], c[0], a[0] ], [ a[1], b[1], c[1], a[1] ], 'b-' )\n",
    "        \n",
    "n = 10\n",
    "for i in range ( 0, n + 1 ):\n",
    "    for j in range ( 0, n - i + 1 ):\n",
    "        k = n - i - j\n",
    "        g = ( i * a + j * b + k * c ) / float ( n )\n",
    "        plt.plot ( g[0], g[1], 'ro' )\n",
    "    \n",
    "plt.savefig ( 'triangles_fig2.png' )\n",
    "plt.close ( )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Random (but not uniform!) points in a triangle ##\n",
    "\n",
    "Instead of a regular grid, we might like to choose random points in the triangle.\n",
    "\n",
    "We saw that grid points can be defined by blending the three vertices using\n",
    "positive coefficients that add up to 1.  So could we get random values by\n",
    "picking three random values I, J and K, normalizing them by their sum?  \n",
    "\n",
    "That does seem random, but unfortunately it's not uniform.  That is, some parts\n",
    "of the triangle will be sampled more often than others.  It's easiest to see this\n",
    "by generating and plotting a bunch of points.  \n",
    "\n",
    "Generate and plot 10 points, then 100, then 1,000 and then judge whether this\n",
    "method is uniform.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## \"Random\" (but not uniform!) points in a triangle\n",
    "#\n",
    "a = np.array ( [4,0] )\n",
    "b = np.array ( [3,4] )\n",
    "c = np.array ( [0,1] )\n",
    "\n",
    "plt.plot ( [ a[0], b[0], c[0], a[0] ], [ a[1], b[1], c[1], a[1] ], 'b-' )\n",
    "\n",
    "for l in range ( 0, 1000 ):\n",
    "    i = np.random.random ( )\n",
    "    j = np.random.random ( )\n",
    "    k = np.random.random ( )\n",
    "    s = i + j + k\n",
    "    g = ( i * a + j * b + k * c ) / s\n",
    "    plt.plot ( g[0], g[1], 'ro' )\n",
    "    \n",
    "plt.savefig ( 'triangles_fig3.png' )\n",
    "plt.close ( )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Uniform random points in a triangle ##\n",
    "\n",
    "It turns out that we have to be more careful in our sampling strategy.\n",
    "\n",
    "One way to generate uniform random values in a triangle is as follows:\n",
    "    \n",
    "    r1 = random\n",
    "    r2 = sqrt ( random )\n",
    "    i = 1 - r2\n",
    "    j = ( 1 - r1 ) * r2\n",
    "    k = r1 * r2\n",
    "    g = ( i * a + j * b + k * c )\n",
    "    \n",
    "Generate and plot 10 points, then 100, then 1,000 and then judge whether this\n",
    "method is uniform."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## Random (uniform) points in a triangle\n",
    "#\n",
    "a = np.array ( [4,0] )\n",
    "b = np.array ( [3,4] )\n",
    "c = np.array ( [0,1] )\n",
    "\n",
    "plt.plot ( [ a[0], b[0], c[0], a[0] ], [ a[1], b[1], c[1], a[1] ], 'b-' )\n",
    "\n",
    "for l in range ( 0, 1000 ):\n",
    "    r1 = np.random.random ( )\n",
    "    r2 = np.sqrt ( np.random.random ( ) )\n",
    "    i = 1 - r2\n",
    "    j = ( 1 - r1 ) * r2\n",
    "    k = r1 * r2\n",
    "    s = i + j + k\n",
    "    g = ( i * a + j * b + k * c )\n",
    "    plt.plot ( g[0], g[1], 'ro' )\n",
    "\n",
    "plt.savefig ( 'triangles_fig4.png' )\n",
    "plt.close ( )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Is a point inside a triangle? ##\n",
    "\n",
    "Suppose we have a triangle A, B, C and a point D.  Is D inside the triangle?\n",
    "\n",
    "To answer this question, imagine drawing the following three triangles:\n",
    "  * D,B,C\n",
    "  * A,D,C\n",
    "  * A,B,D\n",
    "\n",
    "If D is inside the triangle, then each of these triangles will actually \n",
    "have positive area.  For instance, D,B,C will be in counterclockwise order,\n",
    "as will A,D,C and A,B,D.\n",
    "\n",
    "But if D is outside the triangle A,B,C, then at least one of the new triangles\n",
    "will have vertices listed in the clockwise direction, so the area formula\n",
    "will be negative.  \n",
    "\n",
    "If we observe that, then the point is not inside the triangle."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3.0\n",
      "1.0\n",
      "3.5\n"
     ]
    }
   ],
   "source": [
    "## Is a point inside a triangle?##\n",
    "#\n",
    "#  Use the same triangle t5 as in the previous exercise.\n",
    "#\n",
    "#  1) Let D be the centroid of A,B,C and compute the three areas and print them.\n",
    "#     These should be positive, because the centroid is always inside the triangle.\n",
    "#\n",
    "#  2) Let E be the point (1,4).  What are the three areas?  Where is E?\n",
    "#\n",
    "#  3) Let F be the point (2,3).  What are the three areas?  Where is F?\n",
    "#\n",
    "a = np.array ( [4,0] )\n",
    "b = np.array ( [3,4] )\n",
    "c = np.array ( [0,1] )\n",
    "\n",
    "d = ( a + b + c ) / 3\n",
    "dbc = 0.5 * ( d[0] * ( b[1] - c[1] ) + b[0] * ( c[1] - d[1] ) + c[0] * ( d[1] - b[1] ) )\n",
    "adc = 0.5 * ( a[0] * ( d[1] - c[1] ) + d[0] * ( c[1] - a[1] ) + c[0] * ( a[1] - d[1] ) )\n",
    "abd = 0.5 * ( a[0] * ( b[1] - d[1] ) + b[0] * ( d[1] - a[1] ) + d[0] * ( a[1] - b[1] ) )\n",
    "print ( dbc )\n",
    "print ( adc )\n",
    "print ( abd )\n",
    "\n",
    "e = np.array ( [ 1.0, 4.0 ] )\n",
    "\n",
    "f = np.array ( [ 2.0, 3.0 ] )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The Barycentric coordinates of a point in a triangle ##\n",
    "\n",
    "  We saw that, in a triangle A, B, C, we could define points inside the triangle\n",
    "  by picking nonnegative coefficients that added up to 1, forming a \"blend\" of\n",
    "  the vertices.  Thus, a point D would be defined by\n",
    "\n",
    "    alpha * Ax + beta * Bx + gamma * Cx = Dx\n",
    "    alpha * Ay + beta * By + gamma * Cy = Dy\n",
    "    alpha      + beta      + gamma      = 1\n",
    "\n",
    "  So far, we have chosen alpha, beta, and gamma, to compute the coordinates Dx, Dy.\n",
    "\n",
    "  For instance, in the regular grid, (alpha,beta,gamma) was (I/N, J/N, K/N).\n",
    "\n",
    "  Interestingly enough, we can also choose points Dx and Dy, and figure out the\n",
    "  corresponding values of alpha, beta, and gamma, essentially by inverting the\n",
    "  linear system above.\n",
    "\n",
    "  However, the easiest way to compute the values is as ratios of subtriangle areas:\n",
    "\n",
    "    alpha = area ( D,B,C) / area (A,B,C)\n",
    "    beta  = area ( A,D,C) / area (A,B,C)\n",
    "    gamma = area ( A,B,D) / area (A,B,C)\n",
    "\n",
    "  Compare this to the previous exercise.  A point is inside the triangle if alpha, \n",
    "  beta and gamma are positive.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.4\n",
      "0.133333333333\n",
      "0.466666666667\n"
     ]
    }
   ],
   "source": [
    "# The Barycentric coordinates of a point in a triangle\n",
    "#\n",
    "#  What are the barycentric coordinates of points D, E and F from the previous exercise?\n",
    "#\n",
    "#  Note that alpha + beta + gamma should equal 1.\n",
    "#\n",
    "a = np.array ( [4,0] )\n",
    "b = np.array ( [3,4] )\n",
    "c = np.array ( [0,1] )\n",
    "\n",
    "d = ( a + b + c ) / 3\n",
    "abc = 0.5 * ( a[0] * ( b[1] - c[1] ) + b[0] * ( c[1] - a[1] ) + c[0] * ( a[1] - b[1] ) )\n",
    "dbc = 0.5 * ( d[0] * ( b[1] - c[1] ) + b[0] * ( c[1] - d[1] ) + c[0] * ( d[1] - b[1] ) )\n",
    "adc = 0.5 * ( a[0] * ( d[1] - c[1] ) + d[0] * ( c[1] - a[1] ) + c[0] * ( a[1] - d[1] ) )\n",
    "abd = 0.5 * ( a[0] * ( b[1] - d[1] ) + b[0] * ( d[1] - a[1] ) + d[0] * ( a[1] - b[1] ) )\n",
    "alpha = dbc / abc\n",
    "beta =  adc / abc\n",
    "gamma = abd / abc\n",
    "print ( alpha )\n",
    "print ( beta )\n",
    "print ( gamma )\n",
    "\n",
    "e = np.array ( [ 1.0, 4.0 ] )\n",
    "\n",
    "f = np.array ( [ 2.0, 3.0 ] )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## Estimating integrals over a triangle##\n",
    "\n",
    "  Suppose we need to estimate the integral of f(x,y) over some triangle A,B,C.\n",
    "\n",
    "  Unless we have a very simple integrand function, a simple approach is to use random\n",
    "  sampling.  That is, we pick N points inside the triangle uniformly at random,\n",
    "  and estimate the integral by multiplying the area of the triangle by the average of\n",
    "  the function values.\n",
    "\n",
    "  For the triangle t5, estimate the integral of f(x,y) = y * exp ( x ).\n",
    "  Compute the estimate using N = 10, 100, and then 1000 points."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "19.7948360457\n"
     ]
    }
   ],
   "source": [
    "## Estimating integrals over a triangle##\n",
    "#\n",
    "#  For the triangle t5, estimate the integral of f(x,y) = y * exp ( x ).\n",
    "#\n",
    "#  Compute the estimate using N = 10, 100, 1000, and 10,000 points.\n",
    "#\n",
    "a = np.array ( [4,0] )\n",
    "b = np.array ( [3,4] )\n",
    "c = np.array ( [0,1] )\n",
    "\n",
    "area = 0.5 * ( a[0] * ( b[1] - c[1] ) + b[0] * ( c[1] - a[1] ) + c[0] * ( a[1] - b[1] ) )\n",
    "\n",
    "n = 10000\n",
    "estimate = 0.0\n",
    "\n",
    "for l in range ( 0, n ):\n",
    "    i = np.random.random ( )\n",
    "    j = np.random.random ( )\n",
    "    k = np.random.random ( )\n",
    "    s = i + j + k\n",
    "    g = ( i * a * j * b + k * c ) / s\n",
    "    estimate = estimate + g[1] * np.exp ( g[0] )\n",
    "\n",
    "estimate = estimate * area / float ( n )\n",
    "print ( estimate )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}