Chapter 8 - Alice and Bob.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Chapter 8 - Local Realism\n",
    "## with Alice and Bob"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from numpy import sin,cos,pi,sqrt,angle,exp,deg2rad,arange,rad2deg\n",
    "import matplotlib.pyplot as plt\n",
    "from qutip import *\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "H = Qobj([[1],[0]])\n",
    "V = Qobj([[0],[1]])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First define the projection operator for a state at angle $\\theta$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def P(theta):\n",
    "    \"\"\"The projection operator for a state at angle theta\"\"\"\n",
    "    theta_ket = cos(theta)*H + sin(theta)*V\n",
    "    return theta_ket*theta_ket.dag()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create the projection operators for each of the angles, two for Alice, two for Bob"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "Pa1 = P(deg2rad(19))\n",
    "Pa2 = P(deg2rad(-35))\n",
    "Pb1 = P(deg2rad(-19))\n",
    "Pb2 = P(deg2rad(35))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create the state $\\big|\\psi\\big\\rangle = \\sqrt{0.2} \\big|H,H\\big\\rangle + \\sqrt{0.8} \\big|V,V\\big\\rangle$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "psi=sqrt(0.2)*tensor(H,H) + sqrt(0.8)*tensor(V,V)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, find the joint probability that Alice measures A1 and Bob measures B1. We do this by finding the expectation value of the projection operator for the joint state $\\big|\\theta_{A1},\\theta_{B1}\\big\\rangle$. This is formed as the tensor product of the two appropriate projection operators. **In these tensor products, be sure to put Alice's operator first, then Bob's (just like we did for the signal and idler photons).** Each operator acts on the photon corresponding to the order in the `tensor()` function.\n",
    "\n",
    "Notice we'll be using a new function `expect()`. This is equivalent to putting the operator in between the state bra and ket:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "P1 = expect(tensor(Pa1,Pb1),psi)  # joint for A1, B1 (expect 0.09)\n",
    "P2 = psi.dag()*tensor(Pa1,Pb1)*psi\n",
    "P1 == P2.data[0,0]  # The only difference is that we have to pull out the value \n",
    "                    # from the Qobj using the .data[0,0] method so we can compare it to result from `expect`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.09302920042790863"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "P1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Find the conditional probability $P(\\theta_{B2}|\\theta_{A1}) = \\frac{P(\\theta_{B2},\\theta_{A1})}{P(\\theta_{A1})}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.999939256912056"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# B2 conditioned on A1 (expect 1)\n",
    "Prob_b2_a1 = expect(tensor(Pa1,Pb2),psi)\n",
    "#(psi.dag()*tensor(Pa1,Pb2)*psi).data[0,0]   # the joint probability\n",
    "Prob_a1 = expect(tensor(Pa1,qeye(2)),psi)\n",
    "#(psi.dag()*tensor(Pa1,qeye(2))*psi).data[0,0]  # the singular probability\n",
    "Prob_b2a1 = Prob_b2_a1 / Prob_a1  # the conditional probability\n",
    "Prob_b2a1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Find the conditional probability $P(\\theta_{A2}|\\theta_{B1}) = \\frac{P(\\theta_{A2},\\theta_{B1})}{P(\\theta_{B1})}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.999939256912056"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# A2 conditioned on B1 (expect 1)\n",
    "# can do it all on one line:\n",
    "expect(tensor(Pa2,Pb1),psi) / expect(tensor(qeye(2),Pb1),psi)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3.39573005293161e-05"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "expect(tensor(Pa2,Pb2),psi)  # joint for A2, B2 (classically expect 0.09, QM says 0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is what we described in class."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## What if the state was just $|H,H\\rangle$?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "psi2=tensor(H,H)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.7992456137533195"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "expect(tensor(Pa1,Pb1),psi2)  # joint for A1, B1 (expect 0.09)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.6710100716628343"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# B2 conditioned on A1:\n",
    "expect(tensor(Pa1,Pb2),psi2) / expect(tensor(Pa1,qeye(2)),psi2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.6710100716628343"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# A2 conditioned on B1 \n",
    "expect(tensor(Pa2,Pb1),psi2) / expect(tensor(qeye(2),Pb1),psi2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.45025451627296204"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# joint for A2, B2\n",
    "expect(tensor(Pa2,Pb2),psi2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is harder to interpret, but we clearly have different probabilities. Finally, check if we had used a mixed state:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## A mixed state instead of the pure (entangled state).\n",
    "Here we have to use the density matrix (since a `ket` cannot describe a mixed state). First some background:\n",
    "QuTiP has a function that gives the density matrix from a `ket` state: `ket2dm`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "Quantum object: dims = [[2, 2], [2, 2]], shape = [4, 4], type = oper, isherm = True\\begin{equation*}\\left(\\begin{array}{*{11}c}0.200 & 0.0 & 0.0 & 0.0\\\\0.0 & 0.0 & 0.0 & 0.0\\\\0.0 & 0.0 & 0.0 & 0.0\\\\0.0 & 0.0 & 0.0 & 0.800\\\\\\end{array}\\right)\\end{equation*}"
      ],
      "text/plain": [
       "Quantum object: dims = [[2, 2], [2, 2]], shape = [4, 4], type = oper, isherm = True\n",
       "Qobj data =\n",
       "[[ 0.2  0.   0.   0. ]\n",
       " [ 0.   0.   0.   0. ]\n",
       " [ 0.   0.   0.   0. ]\n",
       " [ 0.   0.   0.   0.8]]"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rho_mix = 0.2 * ket2dm(tensor(H,H)) + 0.8 * ket2dm(tensor(V,V))\n",
    "rho_mix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.1688370108679419"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# joint for A1, B1\n",
    "expect(tensor(Pa1,Pb1),rho_mix)  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.5609865818344059"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# B2 conditioned on A1\n",
    "expect(tensor(Pa1,Pb2),rho_mix) / expect(tensor(Pa1,qeye(2)),rho_mix)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.5609865818344059"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# A2 conditioned on B1\n",
    "expect(tensor(Pa2,Pb1),rho_mix) / expect(tensor(Pb1,qeye(2)),rho_mix)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.1766384016124271"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# joint for A2, B2:\n",
    "expect(tensor(Pa2,Pb2),rho_mix)  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see that  $P(\\theta_{B2},\\theta_{A2}) >  P(\\theta_{B1},\\theta_{A1})$ as we said in class for a state that obeys realism."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Now repeat with the pure state but using density matrix techniques.\n",
    "This isn't going to tell us anything new, but it shows how to work with the density matrix if you already know the `ket` state."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "Quantum object: dims = [[2, 2], [2, 2]], shape = [4, 4], type = oper, isherm = True\\begin{equation*}\\left(\\begin{array}{*{11}c}0.200 & 0.0 & 0.0 & 0.400\\\\0.0 & 0.0 & 0.0 & 0.0\\\\0.0 & 0.0 & 0.0 & 0.0\\\\0.400 & 0.0 & 0.0 & 0.800\\\\\\end{array}\\right)\\end{equation*}"
      ],
      "text/plain": [
       "Quantum object: dims = [[2, 2], [2, 2]], shape = [4, 4], type = oper, isherm = True\n",
       "Qobj data =\n",
       "[[ 0.2  0.   0.   0.4]\n",
       " [ 0.   0.   0.   0. ]\n",
       " [ 0.   0.   0.   0. ]\n",
       " [ 0.4  0.   0.   0.8]]"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rho_pure = ket2dm(psi)  # convert from a ket to a density matrix (dm)\n",
    "rho_pure"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The calculations are actually the same in QuTiP, the `expect` function takes either a `ket` state or a density matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.0930292004279086"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# joint for A1, B1\n",
    "expect(tensor(Pa1,Pb1),rho_pure)  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.999939256912056"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# B2 conditioned on A1\n",
    "expect(tensor(Pa1,Pb2),rho_pure) / expect(tensor(Pa1,qeye(2)),rho_pure)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.999939256912056"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# A2 conditioned on B1\n",
    "expect(tensor(Pa2,Pb1),rho_pure) / expect(tensor(Pb1,qeye(2)),rho_pure)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3.3957300529310896e-05"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# joint for A2, B2:\n",
    "expect(tensor(Pa2,Pb2),rho_pure)  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These all agree (as they should)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Explore the angles in more detail:\n",
    "Why these angles, 19 and 35?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "psi=sqrt(0.2)*tensor(H,H) + sqrt(0.8)*tensor(V,V)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "angles = arange(1,90,1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "rads = deg2rad(angles)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Make a list of the probability of joint measurements for a pair of angles:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x110757710>]"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAEACAYAAABI5zaHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAGTJJREFUeJzt3X+wXHdZx/H3JyQVUy4FRcs0MalS29gUKxVur1Ca1QK9\nKUjwx2iig9LRMeNQBK2aoMz0OsMMdNTBOnXUYMUfowShCnGaSOjAcuMQSKAtaX7cJFKaH00pYgea\nVtQkffzj7OaenO69e/buj3P27Oc1c+fu2T337JObvc9+9/k+53sUEZiZWbUsKjoAMzPrPSd3M7MK\ncnI3M6sgJ3czswpycjczqyAndzOzCsqV3CVNSpqRdETSphaPv1DSNkkPSXpY0tt6HqmZmeWmdn3u\nkhYBR4CbgFPAXmB9RMyk9nk38MKIeLeklwCHgUsj4mzfIjczsznlGbmPA0cj4lhEnAG2Ausy+wQw\n1rg9BvyXE7uZWXHyJPdlwInU9snGfWl3A1dLOgV8GXhnb8IzM7OF6NWE6s3AgxFxGfAK4M8kvaBH\nxzYzsw4tzrHPY8CK1Pbyxn1ptwLvA4iIr0j6KrAK+GJ6J0leyMbMbAEiQp3sn2fkvhe4QtJKSRcB\n64FtmX2OAa8DkHQpcCXwyBwBlurrjjvuKDyGYYnLMTmmUYirHzE99VTwuc/Nfl+8OIBgyZLg/vuD\na69Nbl9zTfK1ZEly31NPJT+/EG1H7hFxTtJtwE6SN4N7IuKQpI3Jw7EFeC/wN5L2NX7sdyPiyQVF\nZGZWIadPw2tfCwcOwOrVsH178v3gQbj6ahgfh127Zh+H2dtjY/Mfez55yjJExL8BV2Xu+8vU7cdJ\n6u5mZiPv9GnYvx+uuSb5fuAAnD2bJPTjxy9M5s0EPjEx+/Pp2ws18meo1mq1okNoqYxxOaZ8HFN+\nZYyr25iaI/Ubb0y+r1yZJPElS5KRejOhT0x0NzJvp+1JTD19MikG+XxmZoPSHK0/8wysXZuM1Jcs\ngenpJKF3U2qRRHQ4oerkbmbWpXRdfdWq5L7Dh5OR+q5d3Y/QndzNzAYkW1e/8cbZ0fqOHXDxxd1P\nijYtJLnnmlA1M7NZeTpg+llPz8PJ3cysQ3k7YIrksoyZWQ7pMgwkI/fmSL0XdfX5uOZuZtYH2TLM\nrl3J/YMaqTu5m5n1we7dF06YTk/35kSjvBaS3Ef+JCYzs3auuea5JyKVnUfuZmZzyNbZi5owdVnG\nzKxHWtXZi+qCcVnGzKxHsu2OBw4UHVFnnNzNzBpOn04mT0+fHs46e5rLMmZmFN/uOB/X3M3MFqjo\ndsf5uOZuZrZAw16GyfLI3cxGVrrVcWws2S5DGSbLZRkzs5zK1OrYTt/KMpImJc1IOiJpU4vHf1vS\ng5IekPSwpLOSXtRJIGZmgzTsrY7ttB25S1oEHAFuAk4Be4H1ETEzx/5vAt4VEa9r8ZhH7mZWCs2R\n+6BWduxGvy7WMQ4cjYhjjSfZCqwDWiZ3YAPw4U6CMDMbtLGx8q3B3kt5yjLLgBOp7ZON+55D0ncC\nk8C93YdmZtZb6ZOUIEnoExPVS+zQ+1bInwT+PSK+2ePjmpl1pVmGufHG5HszwVdVnrLMY8CK1Pby\nxn2trKdNSWZqaur87VqtRq1WyxGCmVl3Wk2gluUkpax6vU69Xu/qGHkmVJ8HHCaZUH0c2ANsiIhD\nmf0uAR4BlkfEt+c4lidUzawQwzSBmtWXCdWIOCfpNmAnSRnnnog4JGlj8nBsaez6FuCTcyV2M7Mi\npE9UqvIEapZPYjKzyhqmE5Xm47VlzMxSqn6i0nyc3M2ssqq2GFgnXJYxs0or62JgnfDCYWY28rIr\nPVaBa+5mNtJG7USl+Ti5m1lljPIEapaTu5lVxihPoGa55m5mlVKFCdQsT6ia2cip4gRqlidUzWyk\neAJ1bk7uZja0PIE6Nyd3MxtankCdm2vuZjbUqjiBmuUJVTOzCvKEqplVXvY6qNaak7uZDQ13x+Tn\n5G5mQ8PdMfk5uZvZ0HB3TH6eUDWzoTIK3TFZfZtQlTQpaUbSEUmb5tinJulBSfslfaaTIMzM5pOe\nRB0bg4mJ0UnsC9V25C5pEXAEuAk4BewF1kfETGqfS4DPAW+IiMckvSQivtHiWB65m1lHqnKR6270\na+Q+DhyNiGMRcQbYCqzL7PMLwL0R8RhAq8RuZrYQnkRdmDzJfRlwIrV9snFf2pXAd0n6jKS9kt7a\nqwDNbLR5EnVhFvfwONcBPwFcDOyWtDsi/iO749TU1PnbtVqNWq3WoxDMrIrGxpJSzChNotbrder1\nelfHyFNznwCmImKysb0ZiIi4M7XPJuD5EfEHje2/AnZExL2ZY7nmbmbWoX7V3PcCV0haKekiYD2w\nLbPPJ4AbJD1P0lLgeuBQJ4GYmTV5iYHutS3LRMQ5SbcBO0neDO6JiEOSNiYPx5aImJH0SWAfcA7Y\nEhEH+xq5mVWSu2N6wycxmVmp7N6drB1z9mwyiTo9nfS1jzKvCmlmQ8/dMb3hkbuZlc4oLjEwH1+s\nw8ysglyWMbOh5O6Y3nNyN7NC+QIc/eHkbmaF8tox/eHkbmaFcndMf3hC1cwK5+6Y+blbxsysgtwt\nY2ZDwx0y/eXkbmYD5w6Z/nNyN7OBc4dM/zm5m9nAuUOm/zyhamaFcIdMfu6WMTOrIHfLmFlpuTtm\nsJzczazv3B0zeE7uZtZ37o4ZPCd3M+s7d8cMXq7kLmlS0oykI5I2tXh8jaRvSnqg8fWe3odqZsNq\nbCy50PX0tC94PShtu2UkLQKOADcBp4C9wPqImEntswa4PSLe3OZY7pYxM+tQv7plxoGjEXEsIs4A\nW4F1rZ6/kyc2M7P+yZPclwEnUtsnG/dl/ZikhyTdJ+nqnkRnZkPLrY/FWtyj43wJWBER/y1pLfBx\n4MpWO05NTZ2/XavVqNVqPQrBzMqi2frYPAPVdfbO1Ot16vV6V8fIU3OfAKYiYrKxvRmIiLhznp/5\nKvCjEfFk5n7X3M1GwO7dSU/72bNJh8z0NExMFB3V8OpXzX0vcIWklZIuAtYD2zJPfGnq9jjJm8aT\nmNlIcutj8dqWZSLinKTbgJ0kbwb3RMQhSRuTh2ML8LOSfh04A3wb+Pl+Bm1m5dZsffTCYMXxwmFm\nZiXnhcPMrDDujikXJ3cz65oXBisfJ3cz65oXBisfJ3cz65q7Y8rHE6pm1hO+bF7/+DJ7ZmYV5G4Z\nMzMDnNzNrAtufywvJ3czWxC3P5abk7uZLYjbH8vNyd3MFsTtj+XmbhkzWzC3Pw6GWyHNzCrIrZBm\n1lfujhkeTu5mlou7Y4aLk7uZ5eLumOHi5G5mubg7Zrh4QtXMcnN3TDHcLWNmVkF965aRNClpRtIR\nSZvm2e9Vks5I+ulOgjAzs95qm9wlLQLuBm4GVgMbJK2aY7/3A5/sdZBmVgy3Pg6vPCP3ceBoRByL\niDPAVmBdi/3eAXwM+HoP4zOzgrj1cbjlSe7LgBOp7ZON+86TdBnwloj4c6CjupCZlZNbH4fb4h4d\n50+AdC1+zgQ/NTV1/natVqNWq/UoBDPrpWbr48GDbn0ctHq9Tr1e7+oYbbtlJE0AUxEx2djeDERE\n3Jna55HmTeAlwDPAr0XEtsyx3C1jNkTc+lgOfWmFlPQ84DBwE/A4sAfYEBGH5tj/Q8C/RsQ/t3jM\nyd3MrEN9aYWMiHPAbcBO4ACwNSIOSdoo6dda/UgnAZhZubhDphp8EpOZndfskGmWYnbtcjmmDLzk\nr5l1xR0y1eHkbmbneXGw6nBZxswu4A6Z8vHCYWZmFeSau5mZAU7uZiPPrY/V5ORuNsK8OFh1Obmb\njTC3PlaXk7vZCHPrY3W5W8ZsxLn1sfzcCmlmVkFuhTSzttwdMxqc3M1GiLtjRoeTu9kIcXfM6HBy\nNxsh7o4ZHZ5QNRsx7o4ZPu6WMTOrIHfLmJkZ4ORuVnlufRxNuZK7pElJM5KOSNrU4vE3S/qypAcl\n7ZH0mt6Hamadcuvj6Gpbc5e0CDgC3AScAvYC6yNiJrXP0oj478btlwP/FBE/1OJYrrmbDdDu3Uli\nP3s26ZCZnoaJiaKjsk71q+Y+DhyNiGMRcQbYCqxL79BM7A0vAJ7tJAgz6w+3Po6uxTn2WQacSG2f\nJEn4F5D0FuB9wPcAb+xJdGbWlbEx2LXLrY+jKE9yzyUiPg58XNINwHuB17fab2pq6vztWq1GrVbr\nVQhm1sLYmEsxw6Zer1Ov17s6Rp6a+wQwFRGTje3NQETEnfP8zFeAV0XEk5n7XXM3G4DTp5OlBq65\nxqP1KuhXzX0vcIWklZIuAtYD2zJP/LLU7euAi7KJ3cwGwx0yBjnKMhFxTtJtwE6SN4N7IuKQpI3J\nw7EF+BlJvwT8H/Bt4Of6GbSZza3V4mAuy4weLz9gVjHNkfvBg0mHzK5dLs0MO68tY2aAFwerGid3\nM7MK8sJhZiPK68dYlpO72ZBzd4y14uRuNuR86TxrxcndbMh5/RhrxROqZhXg7phqc7eMmVkFuVvG\nbES4O8bacXI3GzLujrE8nNzNhoy7YywPJ3ezIePuGMvDE6pmQ8jdMaPF3TJmZhXkbhmzCnOHjHXC\nyd1sCLhDxjrl5G42BNwhY51ycjcbAu6QsU55QtVsSLhDZnT1bUJV0qSkGUlHJG1q8fgvSPpy4+vf\nJb28kyDMrL2xseRC107slkfb5C5pEXA3cDOwGtggaVVmt0eAGyPiWuC9wAd7HajZqHF3jHUjz8h9\nHDgaEcci4gywFViX3iEiPh8R32psfh5Y1u6g2ReuX8hms9wdY93Kk9yXASdS2yeZP3n/KrBjvgNm\nX7inTj1324neRpm7Y6xbi3t5MEk/DtwK3DDXPlNTU5w4Afv2QUSNgwdr3Hff7Av5wAFYswYefTSZ\nONq1K/m5/fuTjgHXG20UNLtjDh50d8woqtfr1Ov1ro7RtltG0gQwFRGTje3NQETEnZn9fhi4F5iM\niK/McayIiPMj9+YLd/t2uOWWZHvlyiSxnz2btH3t2AG33z7bJeBkb6PC3THW1K9umb3AFZJWSroI\nWA9syzzxCpLE/ta5Enva2FiSpKenk++XXTa7/dnPXtjPG3Hhx9M9e55bi3S93qog+zp2d4x1I1ef\nu6RJ4C6SN4N7IuL9kjaSjOC3SPog8NPAMUDAmYgYb3GcXH3u6RELXDjK/6M/grVr5x7Zb98Ox455\nVG/DpflpNv0J1a9fa6rsqpB5k/3ixXD55a7X2/DZvTv5NNoctExPJ6N2M6hwcs+aK9nnqdc7wVsZ\nZeeh/Fq1tJFJ7lnNZL9ixezEbKsSzvR0kuQ9krcy8gSqzWVkk3vafCWcZleOR/JWBqdPe6Bh+Ti5\nt5BO9vv3X1jX3LEDli71H5cNnidQrRNO7m2k65pXXZXcNzPjPy4bPE+gWiec3HNojuSfftr1eCuO\nJ1CtE07uHZjrLFl/TLZB8QSq5eXk3qH56vEeyVs/eBLVFqJvF+uoqvTp3dnLmK1Y4SVXrbe8jK8N\n0kgn97TsejfHjj13yVWvYWPd8DK+NkhO7ikeyVs/+SLXNkgjXXNvxz3y1mueRLWF8IRqH7lH3hbC\nE6jWC07ufeYeeeuEz0K1XnG3TJ81a/LXX+96vLXnCVQrkpP7AuTprDHzBKoVyWWZHmh1Kjm4TGOe\nQLXecM29QK2WGnatdfR4AtX6oW81d0mTkmYkHZG0qcXjV0n6nKT/kfRbnQRQFeke+Wytdc8en/w0\nCnwGqpVJ2+QuaRFwN3AzsBrYIGlVZrf/At4B/GHPIxxC6VrrVVfBu97lP/hR4AlUK5M8I/dx4GhE\nHIuIM8BWYF16h4j4RkR8CTjbhxiHTnrC9QMfSPrhvYxB9XkC1cokT3JfBpxIbZ9s3GfzcNvkaEi/\nUWe7qFxztyK5FbLP3DZZXa1q7Om5F7MiLc6xz2PAitT28sZ9CzI1NXX+dq1Wo1arLfRQQ6P5Bw+z\nH92bbZOrV7vDYli1qrH7UnnWC/V6nXq93tUx2rZCSnoecBi4CXgc2ANsiIhDLfa9A3g6Iv54jmNV\nthWyE26brAZfKs8GpW997pImgbtIyjj3RMT7JW0EIiK2SLoU+CIwBjwLPA1cHRFPZ47j5J7R6kLJ\nXqem3NKftMAnKVn/+SSmIeRruQ4XLwZmRfDCYUPIE67Dxb3sNiyc3Eug3RWg3BNfrHS7o3vZbVi4\nLFNCzQnXFStcoilaqzIMuM5ug+WyTEU0R/K+SHfxWpVh3Mtuw8DJvcR8ke7iuQxjw8plmZKb7yLd\nbpvsj+xJZV6T3YrmVsiKc9tk/7nV0crINfeKc9tk/7nV0arCyX3IuG2yP5oT1StXusZu1eCyzJBz\n22T3sqWY7dvh+HHX2K08XJYZQW6b7F62FHP8uFsdbfg5uVeE2yY747NOrepclqkQt03m47NObdi4\nLDPi2k22jvJIPj1S91mnNgqc3CsqT9vkqNTks5fDc0eMjQIn9wrrZCR/6lS1Ev18I/Xjx30ha6s+\n19xHyFw1+cWL4fLL4dFHq9FG2aq18ZZbfDk8G15efsBySy9lsHJlktiHffK1OUp/5hlYu/a5/x5P\nmNqwcnK3jmRPgJpvzRood7JPj9ZXrUruO3zYI3Wrhn5fIPtPmL1A9p0t9vlTYC3wDPC2iHioxT5O\n7iU1Xxvljh1w++3lS/bp1RtbxXzxxR6pWzX0pRVS0iLgbuBmYDWwQdKqzD5rgZdFxA8CG4G/6CSI\nItXr9aJDaGnQcc03+RrRnJCsc/Ag7NlTzGRsepK0OVK/4YZ6yw6Y8fHiWhvL+JoqY0xQzrjKGNNC\n5OmWGQeORsSxiDgDbAXWZfZZB/wdQER8AbhE0qU9jbRPyvofWWRc2TbK669PEueiRfVMsk++r1lz\nYf98Ny2W2Z9tbp86deEbyhe+kDz3s8/WS9cBU8bXVBljgnLGVcaYFiJPcl8GnEhtn2zcN98+j7XY\nx4ZIeiTfTPa33nphsl+yZLbLptlmmB3VZ5P9XMk7PRpPfyJobq9Zc2E7o9R8w5ntVfeJSGaz3Odu\nuYyNwfLlFyb76Wn47GfnKuE8N9m/+tXJV6vknR6NN3/2vvtmtx99NHkjSZdd0m84TuhmF2o7oSpp\nApiKiMnG9mYg0pOqkv4C+ExEfKSxPQOsiYgnMsfybKqZ2QJ0OqG6OMc+e4ErJK0EHgfWAxsy+2wD\n3g58pPFm8M1sYl9IcGZmtjBtk3tEnJN0G7CT2VbIQ5I2Jg/HlojYLukWSf9B0gp5a3/DNjOz+Qz0\nJCYzMxuMgU2oSpqUNCPpiKRNg3reTAz3SHpC0r7UfS+WtFPSYUmflHTJgGNaLunTkg5IeljSbxQd\nl6TvkPQFSQ82Yrqj6JhSsS2S9ICkbSWK6VFJX278vvaUIS5Jl0j6qKRDjdfW9QW/pq5s/H4eaHz/\nlqTfKMHv6Tcl7Ze0T9I/SLqoBDG9s/F311U+GEhyz3Mi1IB8qBFD2mbg/oi4Cvg08O4Bx3QW+K2I\nWA38GPD2xu+msLgi4n+BH4+IVwA/AqyVNF5kTCnvBA6mtssQ07NALSJeERHjJYnrLmB7RPwQcC0w\nU2RMEXGk8fu5DvhRkvLtvxQZk6TLgHcA10XED5OUqTcUHNNq4FeAV5L87b1J0ssWFFNE9P0LmAB2\npLY3A5sG8dwtYlkJ7EttzwCXNm6/FJgpIq5UPB8HXleWuIClwBeBVxUdE7Ac+BRQA7aV5f8P+Crw\n3Zn7CosLeCHwlRb3F/67ajz3G4BdRccEXAYcA15Mkti3Ff23B/ws8MHU9nuA3wEOdRrToMoyeU6E\nKsr3RqOzJyK+BnxvUYFIupzk3frzJP+RhcXVKH88CHwN+FRE7C06JuADJC/09ERR0THRiOdTkvZK\n+tUSxPX9wDckfahRBtkiaWnBMaX9PPCPjduFxRQRp4A/Bo6TnHj5rYi4v8iYgP3AaxtlmKXALcD3\nLSQmn8T0XIXMMEt6AfAx4J0R8XSLOAYaV0Q8G0lZZjkw3vi4WFhMkt4IPBHJgnTztdQW8f/3mkjK\nDbeQlNVe2yKOQca1GLgO+LNGXM+QfFou9DUFIGkJ8Gbgo3PEMMjX1ItIlk5ZSTKKv1jSLxYZU0TM\nAHeSfELdDjwInGu1a7tjDSq5PwasSG0vb9xXBk8018GR9FLg64MOQNJiksT+9xHxibLEBRARTwF1\nYLLgmF4DvFnSI8CHgZ+Q9PfA14r+PUXE443v/0lSVhun2N/VSeBERHyxsX0vSbIvw2tqLfCliPhG\nY7vImF4HPBIRT0bEOZI5gFcXHBMR8aGIeGVE1IBvAocXEtOgkvv5E6EkXURyItS2AT13lrhw5LcN\neFvj9i8Dn8j+wAD8NXAwIu5K3VdYXJJe0pyNl/SdwOtJan6FxRQRvxcRKyLiB0heP5+OiLcC/1pU\nTACSljY+dSHpYpJ68sMU+7t6Ajgh6crGXTcBB4qMKWUDyZtzU5ExHQcmJD1fkkh+TwcLjglJ39P4\nvgL4KZISVucxDXCiYJLkHegosHlQz5uJ4R+BU8D/kvzH3koymXJ/I7adwIsGHNNrSD52PUTyEeyB\nxu/qu4qKC3h5I46HgH3A7zfuLyymTHxrmJ1QLTQmkvp28//u4eZruwRxXUsyqHoI+GfgkhLEtBT4\nT2AsdV/RMd1BMnDZB/wtsKQEMU2T1N4fJOnCWtDvyScxmZlVkCdUzcwqyMndzKyCnNzNzCrIyd3M\nrIKc3M3MKsjJ3cysgpzczcwqyMndzKyC/h8vJqc3jleqVwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10e2a0518>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "out = []\n",
    "for r in rads:\n",
    "    out.append(expect(tensor(P(-r),P(r)),psi))\n",
    "\n",
    "plt.plot(angles,out,\".\") # plot in units of pi"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see that the joint probabilities have a zero at 35˚. Now plug that in to one of the conditional probabilities and see what angle for the conditional probability gives 1:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x110c4ef60>]"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAEACAYAAABI5zaHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFytJREFUeJzt3X+QnPdd2PH3R5FcULg4kIDBFlIgxg6ywWCIEE2VLHHA\nsgtxSpki0wFKaeNpa0jSHyOFMuPrTGeI/2BCmHTKGEygtImmkDRRW2Viu2Fr0zpEwXJs/ThLxMlZ\nthODmwbLIZPKyqd/PLvWo9Webm9vd58f+37N3Nw+u889+7m9vc9+9/N8vt+NzESS1C4bqg5AkjR5\nJndJaiGTuyS1kMldklrI5C5JLWRyl6QWWjW5R8TdEfFMRDxykX1+MyJORsTDEfF9kw1RkrRWo4zc\n3wfcuNKNEXET8OrM/C7gNuC3JhSbJGlMqyb3zPwT4P9eZJdbgP/Q2/dPgUsj4rLJhCdJGsckau5X\nAKdK20/1rpMkVcQTqpLUQhsncIyngG8vbW/pXXeBiHAhG0kaQ2bGWvYfdeQeva9hDgA/BxARO4Ev\nZeYzFwmwVl933HHHxI/53HPJddclGzcmV15ZfIdk06bkvvuK2zZtKr4/9dTw7Q0b7uDaa5Nrrz13\n2733Dj/Wxo3F9+eea95jZUzzG1Nd46pjTONYdeQeEe8HOsArIuIJ4A7gkiJP512ZeTAibo6IPwe+\nDPzCWJE03OnTcOQIXHtt8f3oUXjhBfjc5+BVr4LlZdi+HXbsgAceKG6/5hpYWBi+/Y53wLvfXRy7\nfxsU348dK46Vee5+jh2DT34SNm8uYlhYqOqRkFQHqyb3zPyZEfa5fTLhNNPp07Br17kkfPDg+Un4\n4EF44olzyRtg585zP7+wcOH2li3D9y2/EMC5+7n6anj722FpqbjugQdM8NI8m0TNvdE6nc66j1Ee\nqR87ViTywdH45ZdPJq7BF4L+/Tz/PNx007kY+vfdfzcxiUQ/icdq0oxpNHWMCeoZVx1jGkeMW88Z\n684icpb3N03lMgwUI/f+SL2KUXP/3UP53cLNN59L8o7kpeaKCHKNJ1RN7mMYLMM88EBxfXmkXlVc\n5RH7619fjOQ3bYL775/8SF7SbIyT3O1zH8NgGebo0XPlkiqTZjmGa68tkvmmTcVIfuvW4gXp9a8v\nvp8+XV2ckqbP5L4Gp0/Dgw/Ctm3nJ87+yc066Xfd3H9/8X15+cIXpP7vY6KX2seyzIiGdcQMdsDU\nmTV5qbmsuU/Rgw9eWMMud600wWo1+ab9PtK8sOY+YeWyxWANu46lmNWsVpO3RCO1hyP3FdS1I2aS\n+iP5rVst0Uh15sh9guraETNJ/d9n2MlWSc1mcl9BG8owoxr2u9pJIzWbZZmLKJ+AbMtofSXl3xUu\nLEm1/feX6sxumXUqLykwz8msDZ1BUptYc1+H/glUZ3BappHawOTeM+wE6rwanN0KvvBJTWNy75mn\nE6ijKHcG+cInNY8195J5OoG6FoNLF/RH856fkGbDE6pr5AnU0dlNI1XHE6pr4AnUtbFMIzXL3CZ3\nE9T4XJdGqr+5LcsMqyNbWhid69JIs2PNfY08gbp+TniSps+a+wjKk3HathBYFZzwJNXTSMk9InZH\nxFJEnIiIvUNuf3lEfCgiPh0Rn4iI7ZMPdf08iTp5TniS6mnV5B4RG4D3AjcC1wC3RsRrBnb7FeBw\nZl4H/Dzwm5MOdBI8iToddtJI9TPKyH0HcDIzlzPzDLAfuGVgn+3AxwEy8zHgVRHxzRONdAKchTp9\nPsZSPYyS3K8ATpW2n+xdV/Zp4CcBImIHsBXYMokAJ2mwhGCtffKGPcbW4KXZ2zih47wLeE9EPAQ8\nChwGzg7bcXFx8cXLnU6HTqczoRBG0y8haHrKj/Gwjyv0RVW6uG63S7fbXdcxVm2FjIidwGJm7u5t\n7wMyM++8yM98FviezHx+4PqZt0K6xEC1bJWU1m9arZCHgCsjYltEXALsAQ4M3PGlEbGpd/kfA/9z\nMLFXwe6Y6jmbVarGqmWZzDwbEbcD91C8GNydmccj4rbi5rwL+G7g9yPia8BR4BenGfSohnVuOGqc\nrX4N3tms0my1eoaqSwzUiyUaaTwuPzCESwzUhy+20nhM7qq9wRdbT3hLqzO5q1Fsk5RG48JhOGGm\nSVyqQJqeViV3Wx+bxRUlpelpVXJ3JNgsrigpTU+rkruLVjWPK0pK09G6E6q2PjaXrZLScHbLqPFs\nlZQuZHJXq9gqKRVshVSrWIOXxteK5G77XDt5glwaX+PLMr51b7dyDR6sv2s+zWVZxrfu7Vb+VCd7\n4KXRNT65+9Z9PvgiLq1N48syYG/7PBjWAw+WaTQfbIVUqw3W3z3Xonlhctfc8FOdNE/m8oSq5pPn\nWqSLa+TI3SnpApcq0PyYi7KMfe0axueF2mwuyjK2xGkYnxfS+RqX3K21ahifF9L5RirLRMRu4Dco\nXgzuzsw7B25/GfAfga3AS4Bfz8zfG3KcidXc7WvXIJcqUFtNpeYeERuAE8ANwNPAIWBPZi6V9nkn\n8LLMfGdEvBJ4DLgsM18YOJatkJo66+9qm2nV3HcAJzNzOTPPAPuBWwb2SaD/77MA/J/BxC7NivV3\nabTkfgVwqrT9ZO+6svcC2yPiaeDTwNsmE560dsPq7y4LrXmzcULHuRE4nJlvjIhXA/dGxPdm5vOD\nOy4uLr54udPp0Ol0Vj24/ctai4WFohTjUgVqqm63S7fbXdcxRqm57wQWM3N3b3sfkOWTqhHx34Bf\ny8z/1dv+H8DezPzUwLHWXHO3fqr1cqkCNd20au6HgCsjYltEXALsAQ4M7LMMvKkXxGXAVcDjawlk\nJdZPtV62SWoeraUV8j2ca4V8V0TcRjGCvysivg34PeDbej/ya5n5gSHHGXvkXl7q1ZG71sqlCtRk\nrV1+wL52TZKlPjVNa5O7NEnW4NU0c7G2jLRe1uA1Dxy5ay5Zg1eTWJaRxmANXnVnWUYag+22aqNa\nJnenimuWrMGrjWpXlvEtsqrgcsGqs1aUZXyLrCosLJxrh9y1q2iV3LXLd49qrtold98iq0oOLtQW\ntSvLgDNSVZ1hy12AZRpVy1ZIaQIG6++eA1LVTO7ShLlUgeqgFSdUpTrxHJCaypG7tAqXKlDVGl2W\n8R9GTeA8DFWhsWWZ/j+MvcWqO1sl1RS1SO7+w6gprMGrKWpRlvGj9NQkLlWgWWt8zd2JS2oS6++a\nlUYnd6lp7IHXrDT2hKrURMPq7y5Xrbpw5C6tg0sVaBYsy0gVskyjaZlaWSYidkfEUkSciIi9Q27/\nlxFxOCIeiohHI+KFiHj5WgKRms42SdXJqiP3iNgAnABuAJ4GDgF7MnNphf1/HHh7Zr5pyG0vjtyd\nkao2cqkCTcO0Ru47gJOZuZyZZ4D9wC0X2f9W4AMXO6AzUtVW/U906id2n+eqyijJ/QrgVGn7yd51\nF4iIrwd2Ax+82AGdkap54PNcVdo44eP9BPAnmfmllXZYXFzkq1+FV7wCnn22w/btHWuTaqV+Db4/\n89rnuUbV7XbpdrvrOsYoNfedwGJm7u5t7wMyM+8csu+HgP+cmftXONZ5NXdnpKrtrMFrEqbSChkR\nLwEeozih+nngk8CtmXl8YL9LgceBLZn5lRWOZSuk5pbLFWhcUzmhmplngduBe4CjwP7MPB4Rt0XE\nW0u7vgX42EqJXZp31uA1S05ikmZkcPXTgwdhedkSjVbnDFWp5vo1+K1b4eabLdFoNC4cJtVcvw9+\nedkSjabL5C5VwKUKNG2WZaSK2CapUVlzlxrKNkldjDV3qaFsk9SkmdylGvBTnTRplmWkmvBTnbQS\na+5SS/ipTiqz5i61hK2SWi9H7lJN2SqpPssyUkvZKjnfLMtILWWrpNbK5C41wGANfutW2yR1cZZl\npIZwRcn5ZVlGajFXlNRamNylhrFNUqOwLCM1kG2S88VWSGkO2SbZftbcpTlkm6SGMblLDeeKkhrG\nsozUAq4o2W7W3CW5omQLTa3mHhG7I2IpIk5ExN4V9ulExOGIOBIRf7yWICRNjrNZBSOM3CNiA3AC\nuAF4GjgE7MnMpdI+lwL/G/ixzHwqIl6Zmc8OOZYjd2kGnM3aLtMaue8ATmbmcmaeAfYDtwzs8zPA\nBzPzKYBhiV3S7DibVaMk9yuAU6XtJ3vXlV0FfFNE/HFEHIqIn51UgJLG52zW+bVxgse5Hngj8FLg\nwYh4MDP/fHDHxcXFFy93Oh06nc6EQpA0aGGhKMU4m7VZut0u3W53XccYpea+E1jMzN297X1AZuad\npX32Al+Xmf+mt/07wEcz84MDx7LmLlXI2azNNK2a+yHgyojYFhGXAHuAAwP7fAT4WxHxkojYDPwQ\ncHwtgUiaPmezzo9Vk3tmngVuB+4BjgL7M/N4RNwWEW/t7bMEfAx4BPgEcFdmHpte2JLG4WzW+eEk\nJmnOOJu1eZyhKmlNnM3aDK4KKWlNLNO0lyN3ac5Zpqk/yzKS1sUyTT1ZlpG0Li461h6O3CWdx0XH\n6seRu6R1c9GxdjC5SxrKRceazbKMpBWVO2lcdKw6dstImhoXHauONXdJU+OiY81icpc0EmezNotl\nGUkjczZrNay5S5oZZ7POjjV3STPjbNZ6c+QuaWzOZp0NR+6SZsrZrPVlcpe0bnbS1I9lGUkTYSfN\n9NgtI6kW7KSZLGvukmrBMk31HLlLmgrLNJNjWUZSLVmmWZ+plWUiYndELEXEiYjYO+T2N0TElyLi\nod7Xr64lCEntZplm9lYduUfEBuAEcAPwNHAI2JOZS6V93gD8i8x88yrHcuQuzSnLNOOb1sh9B3Ay\nM5cz8wywH7hl2P2v5Y4lzZf+hKeFBZcPnoVRkvsVwKnS9pO96wb9cEQ8HBH/PSK2TyQ6Sa3kujTT\nt3FCx/kzYGtm/nVE3AR8GLhq2I6Li4svXu50OnQ6nQmFIKkpFhaKUozr0gzX7XbpdrvrOsYoNfed\nwGJm7u5t7wMyM++8yM98FviBzPziwPXW3CWdx06a1U2r5n4IuDIitkXEJcAe4MDAHV9WuryD4kXj\ni0jSKuykmY5VyzKZeTYibgfuoXgxuDszj0fEbcXNeRfwUxHxT4AzwFeAn55m0JLao1yisZNmcpzE\nJKlWLNNcyLVlJDWeZZrJcOQuqXac8HQ+15aR1DqWaSzLSGohJzyNx5G7pNqb9w/iduQuqZX8IO61\nM7lLagw7aUZnWUZSo8xjJ43dMpLmyrx00lhzlzRXLNOszJG7pEabhzKNZRlJc62tZRrLMpLmmmWa\ncxy5S2qVNpZpLMtIUklbyjSWZSSpZJ7XpXHkLqnV2rAujSN3SRowr+vSmNwlzYV566SxLCNpbjS1\nk8ZuGUkaUZM6aay5S9KI2l6mceQuaW41pUwztZF7ROyOiKWIOBERey+y32sj4kxE/ORagpCkKvQ7\naRYW4MiRC7tpmjySXzW5R8QG4L3AjcA1wK0R8ZoV9nsX8LFJBylJ0zZswtOuXUVdfteu5iX4UUbu\nO4CTmbmcmWeA/cAtQ/b7JeCPgL+YYHySNBMLC0Up5v77i+9N74sfJblfAZwqbT/Zu+5FEXE58JbM\n/PfAmupCklQX5TJN05cu2Dih4/wGUK7Fr5jgFxcXX7zc6XTodDoTCkGSJqc/kq9i6YJut0u3213X\nMVbtlomIncBiZu7ube8DMjPvLO3zeP8i8Ergy8BbM/PAwLHslpHUOFX3xE+rW+YQcGVEbIuIS4A9\nwHlJOzO/s/f1HRR19386mNglqama2BO/alkmM89GxO3APRQvBndn5vGIuK24Oe8a/JEpxClJlSmX\naOreE9/nJCZJWqNZl2lcfkCSZqAJZRpH7pI0hlkuXeCqkJJUgcEyzUc/Cps3FyP8SSR5k7skVeD0\n6WLkfuwYXH11cd3S0uRG8SZ3SapIv0zz/PNw002TPdlqcpekipVH8du3FyN3KFadHLdMY3KXpBqY\n9MlWk7sk1cwkeuLtc5ekmqmqJ96RuyRN2XrLNI7cJamGqvg4P5O7JM3QrD7Oz7KMJM1YuUxz5Mjq\nJ1wty0hSA8zi4/wcuUtSxfoj+ZU+zs+RuyQ1UH8kv7x84cnWcZncJakmhvXEj8uyjCTVSPlka7//\n3eUHJKmFrLlLkgCTuyS1ksldklpopOQeEbsjYikiTkTE3iG3vzkiPh0RhyPikxHxusmHKkka1arJ\nPSI2AO8FbgSuAW6NiNcM7HZfZl6Xmd8P/CLwOxOPdEq63W7VIQxVx7iMaTTGNLo6xlXHmMYxysh9\nB3AyM5cz8wywH7ilvENm/nVp8xuAr00uxOmq6x+yjnEZ02iMaXR1jKuOMY1jlOR+BXCqtP1k77rz\nRMRbIuI48F+BfziZ8CRJ45jYCdXM/HBmfjfwFuDfTuq4kqS1W3USU0TsBBYzc3dvex+QmXnnRX7m\nM8BrM/OLA9c7g0mSxrDWSUwbR9jnEHBlRGwDPg/sAW4t7xARr87Mz/QuXw9cMpjYxwlOkjSeVZN7\nZp6NiNuBeyjKOHdn5vGIuK24Oe8C/m5E/Bzw/4CvAH9vmkFLki5upmvLSJJmY2YzVFebCDWjGO6O\niGci4pHSdd8YEfdExGMR8bGIuHTGMW2JiI9HxNGIeDQifrnquCLib0TEn/YmpT0aEXdUHVMptg0R\n8VBEHKhRTJ8rT+KrQ1wRcWlE/GFEHO89t36o4ufUVb3H56He97+KiF+uweP0jog4EhGPRMR/iohL\nahDT23r/d+vKBzNJ7iNOhJqF9/ViKNtHMQnrauDjwDtnHNMLwD/PzGuAHwb+We+xqSyuzPwq8CO9\nSWnfB9wUETuqjKnkbcCx0nYdYvoa0MnM78/MHTWJ6z3AwV4H23XAUpUxZeaJ3uNzPfADwJeB/1Jl\nTBFxOfBLwPWZ+b0UZepbK47pGoqJoD9I8b/34xHx6rFiysypfwE7gY+WtvcBe2dx30Ni2QY8Utpe\nAi7rXf5WYKmKuErxfBh4U13iAjYDnwJeW3VMwBbgXqADHKjL3w/4LPCKgesqiwt4GfCZIddX/lj1\n7vvHgAeqjgm4HFgGvpEisR+o+n8P+Cngt0vbvwr8K+D4WmOaVVlmpIlQFfmWzHwGIDO/AHxLVYFE\nxKsoXq0/QfGHrCyuXvnjMPAF4N7MPFR1TMC7KZ7o5RNFVcdEL557I+JQRPyjGsT1HcCzEfG+Xhnk\nrojYXHFMZT8NvL93ubKYMvNp4NeBJ4CngL/KzPuqjAk4AuzqlWE2AzcD3z5OTK4KeaFKzjBHxDcA\nfwS8LTOfHxLHTOPKzK9lUZbZAuzovV2sLKaI+NvAM5n5MHCxltoq/n6vy6LccDNFWW3XkDhmGddG\n4Hrg3/Xi+jLFu+VKn1MAEbEJeDPwhyvEMMvn1MspllLZRjGKf2lE/P0qY8rMJeBOineoB4HDwNlh\nu652rFkl96eAraXtLb3r6uCZiLgMICK+FfiLWQcQERspEvsfZOZH6hIXQGY+B3SB3RXH9DrgzRHx\nOPAB4I0R8QfAF6p+nDLz873vf0lRVttBtY/Vk8CpzPxUb/uDFMm+Ds+pm4A/y8xne9tVxvQm4PHM\n/GJmnqU4B/A3K46JzHxfZv5gZnaALwGPjRPTrJL7ixOhIuISiolQB2Z034OC80d+B4B/0Lv888BH\nBn9gBn4XOJaZ7yldV1lcEfHK/tn4iPh64Ecpan6VxZSZv5KZWzPzOymePx/PzJ+lWMuokpgAImJz\n710XEfFSinryo1T7WD0DnIqIq3pX3QAcrTKmklspXpz7qozpCWBnRHxdRATF43Ss4piIiG/ufd8K\n/B2KEtbaY5rhiYLdFK9AJ4F9s7rfgRjeDzwNfJXiD/sLFCdT7uvFdg/w8hnH9DqKt10PU7wFe6j3\nWH1TVXEB39OL42HgEeBf966vLKaB+N7AuROqlcZEUd/u/+0e7T+3axDXdRSDqoeBDwGX1iCmzcBf\nAgul66qO6Q6KgcsjwO8Dm2oQ0/0UtffDFF1YYz1OTmKSpBbyhKoktZDJXZJayOQuSS1kcpekFjK5\nS1ILmdwlqYVM7pLUQiZ3SWqh/w9wfwEzlTBJCwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1107a3f60>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "out = []\n",
    "for r in rads:\n",
    "    out.append(expect(tensor(P(r),P(deg2rad(35))),psi) / expect(tensor(P(r),qeye(2)),psi))\n",
    "\n",
    "plt.plot(angles,out,\".\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So only 19 and 35 work. Now, can you derive 19 and 35 given only the state $|\\psi\\rangle$? Try the first plot, i.e. calculate the joint probability $P(\\theta_A,\\theta_B)$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Solution\n",
    "Using the state, write the projection operators for a two photon state with angles $\\theta_A$ and $\\theta_B$. First, recall $$\\big|\\theta_i\\big\\rangle = \\cos\\theta_i\\big|H\\big\\rangle + \\sin\\theta_i\\big|V\\big\\rangle.$$ Next, form the two-photon state: $$\\big|\\theta_A,\\theta_B\\big\\rangle = \\big|\\theta_A\\big\\rangle \\otimes \\big|\\theta_B\\big\\rangle = \\left(\\cos\\theta_A\\big|H\\big\\rangle + \\sin\\theta_A\\big|V\\big\\rangle\\right) \\otimes \\left(\\cos\\theta_B\\big|H\\big\\rangle + \\sin\\theta_B\\big|V\\big\\rangle\\right)$$\n",
    "which we can reduce to:\n",
    "$$=\\cos\\theta_A\\cos\\theta_B\\big|H,H\\big\\rangle + \\cos\\theta_A\\sin\\theta_B\\big|H,V\\big\\rangle + \\sin\\theta_A\\cos\\theta_B\\big|V,H\\big\\rangle + \\sin\\theta_A\\sin\\theta_B\\big|V,V\\big\\rangle.$$\n",
    "Find the probability of a joint measurement of polarizations $\\theta_A$ and $\\theta_B$:\n",
    "$$P(\\theta_A,\\theta_B) = \\big|\\big\\langle\\psi\\big|\\theta_A,\\theta_B\\big\\rangle\\big|^2$$\n",
    "Since $\\big|\\psi\\big\\rangle$ only has $\\big|H,H\\big\\rangle$ and $\\big|V,V\\big\\rangle$ terms, this probability only has two terms:\n",
    "$$P(\\theta_A,\\theta_B) = \\left|\\sqrt{0.2}\\cos\\theta_A\\cos\\theta_B + \\sqrt{0.8}\\sin\\theta_A\\sin\\theta_B\\right|^2$$\n",
    "Plot is shown below for $\\theta_A = -\\theta_B$ and it agrees perfectly with our model above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x11147ed30>]"
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAEACAYAAABI5zaHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X2c1WP+x/HXJ7S7ITcbWaF1syikLG0Uzqpl3G1uFmXd\nxmqtrBa72v2hWatoabGEjcjNknVfbsN2tBQNlUJTIW0prfst3U31+f1xnewxpubMzDnnOjfv5+Mx\nj+ac+c45bzPjc77nc13f6zJ3R0RESkuz2AFERCT7VNxFREqQiruISAlScRcRKUEq7iIiJUjFXUSk\nBGVU3M2swsyqzWyWmV1Sx9dbmtloM5tqZtPN7IysJxURkYxZffPczawZMAvoDiwAqoBe7l6ddszv\ngZbu/nszawXMBFq7+6qcJRcRkXXK5My9MzDb3ee6ew0wCuhZ6xgHNk19vinwiQq7iEg8mRT3NsC8\ntNvzU/eluwlob2YLgDeAC7ITT0REGiNbA6qHAVPcfVugEzDMzDbJ0mOLiEgDbZjBMR8AO6Td3i51\nX7ozgasA3P1dM5sD7A68ln6QmWkhGxGRRnB3a8jxmZy5VwG7mFlbM2sO9AJG1zpmLtADwMxaA7sC\n760jYEF9DBw4MHqGYsmlTMpUDrlylWnNGmfiROecc5wttnC6dnX++MdwX03N+r+3Meo9c3f31WbW\nDxhLeDEY4e4zzKxv+LIPB64ERprZtNS3/c7dP21UIhGRErJ8Odx5J9x4I9TUQJ8+MH06tKk9cpll\nmbRlcPdngN1q3fe3tM8XEvruIiICLF0Kw4fDNddAp05w661w4IFgDWquNF5Gxb2UJRKJ2BHqVIi5\nlCkzypS5QszV1EyrV8OIEXD55dC1KzzxRCju+VbvRUxZfTIzz+fziYjk08SJcP758O1vhzZMtoq6\nmeENHFAt+zN3EZGm+uILuPBCeOYZGDIEfv7z/LVf1kULh4mINMG4cdChAzRvDjNmwCmnxC/soDN3\nEZFGWb4c/vAHeOABuP12OPzw2Im+TsVdRKSB5syB446DnXeGadPgu9+Nneib1JYREWmAsWOhSxc4\n/XR48MHCLOygM3cRkYy4h8HSG24IrZgCnMX5NSruIiL1WLkSzj47DJhOmgTbbx87Uf1U3EVE1mPx\nYjj++DB3/cUXoUWL2Ikyo567iMg6LFwIBx8MO+0EjzxSPIUdVNxFROr03nth+YDjj4dbboENi6zP\nUWRxRURyb/Zs6N49zGP/5S9jp2kcFXcRkTQzZ0KPHjBwYBhELVYq7iIiKW+/DT/5CQweHOaxFzMV\ndxER4J13QmEfMiSsD1PsNKAqImVv/vxQ2CsrS6Owg4q7iJS5jz+GQw+Fc8+FX/widprs0WYdIlK2\n/vtfOOSQcNZ+1VWx06xbYzbryOjM3cwqzKzazGaZ2SV1fP1iM5tiZpPNbLqZrTKzzRsSREQkn2pq\nwhz2ffcNA6ilpt4zdzNrBswCugMLgCqgl7tXr+P4o4D+7t6jjq/pzF1EonMP0xwXLYLHHiv8C5Ry\ntc1eZ2C2u89NPckooCdQZ3EHegP3NySEiEg+DR4MU6bA+PGFX9gbK5O2TBtgXtrt+an7vsHMvgNU\nAA83PZqISPbddx8MHw5PPAGbbBI7Te5k+zXraOAld/88y48rItJkEyZA//7wwguw7bax0+RWJsX9\nA2CHtNvbpe6rSy/qaclUVlZ+9XkikSBR6Cvei0hJmD8fTjgBRo6EvfaKnWb9kskkyWSySY+RyYDq\nBsBMwoDqQmAS0NvdZ9Q6bjPgPWA7d1+2jsfSgKqI5N2yZXDQQWF2zIABsdM0XE4GVN19tZn1A8YS\nevQj3H2GmfUNX/bhqUOPAZ5dV2EXEYnBHfr2DZtZX/KNidylSxcxiUhJu+46uPtuePnl4tpsI11j\nztxV3EWkZI0fDyeeCK++Cm3bxk7TeDm7QlVEpNgsWgQnnxwGUIu5sDeWiruIlJzVq6F3b+jTByoq\nYqeJQ8VdRErOwIHQrFn4t1yV6IW3IlKunn4a7roLXn8dNtggdpp4VNxFpGQsXBhaMQ88AFtvHTtN\nXGrLiEhJWLMGTjstzGk/6KDYaeJTcReRknDttbB8OVx6aewkhUHz3EWk6E2aBEcdBa+9BjvsUP/x\nxUbz3EWk7CxeHOaz33xzaRb2xtKZu4gUtT59wqyY226LnSR3crUTk4hIQXr8cXjxRXjjjdhJCo/O\n3EWkKC1aBB07wkMPQdeusdPklhYOE5Gy4A49e8Kee4b9UEud2jIiUhbuuAPmzQtn7VI3nbmLSFF5\n/33Ybz8YNy6cuZcDTYUUkZLmDmedBRdfXD6FvbFU3EWkaPztb7BkCVx0UewkhU9tGREpCu+/D/vu\nG3ZXat8+dpr8yllbxswqzKzazGaZWZ1bzJpZwsymmNmbZjauISFERNbHHc4+O7Rjyq2wN1a9Z+5m\n1gyYBXQHFgBVQC93r047ZjNgAnCou39gZq3c/eM6Hktn7iLSYMOHhytQJ06EDctwjl+upkJ2Bma7\n+9zUk4wCegLVacecDDzs7h8A1FXYRUQaY/58+MMfIJksz8LeWJm0ZdoA89Juz0/dl25XYEszG2dm\nVWZ2arYCikj5codf/Qr69dPsmIbK1uvghsA+wCHAxsBEM5vo7u/UPrCysvKrzxOJBIlEIksRRKTU\nPPggvPtu+LecJJNJkslkkx4jk557F6DS3StStwcA7u5D0o65BPi2u/8xdft24Gl3f7jWY6nnLiIZ\n+fRT2GMPeOQR2H//2GniytVsmSpgFzNra2bNgV7A6FrHPA50M7MNzKwF8CNgRkOCiIiku/hiOOEE\nFfbGqrct4+6rzawfMJbwYjDC3WeYWd/wZR/u7tVm9iwwDVgNDHf3t3OaXERK1gsvhI8334ydpHjp\nIiYRKSjLlkGHDnD99XDkkbHTFAatLSMiRW/w4LBOuwp70+jMXUQKxowZcNBBYWelbbeNnaZw6Mxd\nRIqWO/zyl3D55Srs2aDiLiIFYeRIWLo0XLQkTae2jIhE9/HHYU7700/DPvvETlN4tIeqiBSls86C\nTTcNM2Tkm7SHqogUnZdfhmeeCYOpkj3quYtINKtWhR770KHQsmXsNKVFxV1EornpJthqKzjppNhJ\nSo967iISxYIF4UrUl1+G3XaLnaawaUBVRIpG796w004waFDsJIVPA6oiUhTGjQtb5o0YETtJ6VLP\nXUTyqqYGzj8f/vIXaNEidprSpeIuInk1bFhYXuDYY2MnKW3quYtI3nz4YdgL9aWXYPfdY6cpHhpQ\nFZGCdsYZsPXW8Oc/x05SXDSgKiIFa8IEeP55XYmaL+q5i0jOrV4N/frBkCFhDRnJPRV3Ecm5O+4I\nM2NOPjl2kvKRUXE3swozqzazWWZ2SR1fP9jMPjezyamPS7MfVUSK0WefwWWXwY03gjWoayxNUW/P\n3cyaATcB3YEFQJWZPe7u1bUOHe/uP81BRhEpYpWVcMwx0KlT7CTlJZMB1c7AbHefC2Bmo4CeQO3i\nrtdkEfmat96C++7TIGoMmbRl2gDz0m7PT91X2/5mNtXMnjSz9llJJyJFyx0uuCDsidqqVew05Sdb\nUyFfB3Zw96VmdjjwGLBrXQdWVlZ+9XkikSCRSGQpgogUksceCxctnXtu7CTFJ5lMkkwmm/QY9V7E\nZGZdgEp3r0jdHgC4uw9Zz/fMAX7o7p/Wul8XMYmUgeXLoX17uO026N49dpri15iLmDJpy1QBu5hZ\nWzNrDvQCRtd64tZpn3cmvGh8ioiUpeuug733VmGPqd62jLuvNrN+wFjCi8EId59hZn3Dl3048DMz\nOxeoAZYB2ldFpEwtWBC2zXv11dhJypvWlhGRrDrjDNhmG7j66thJSofWlhGRqCZNgrFjYebM2ElE\nyw+ISFa4Q//+MHiw1o8pBCruIpIV998PK1fCaafFTiKgnruIZMGXX0K7duFq1G7dYqcpPbmaCiki\nsl7XXgsHHKDCXkh05i4iTTJvHnTsCJMnQ9u2sdOUJm2zJyJ5d8op8P3vw5VXxk5SujQVUkTy6pVX\nYNw4TX0sROq5i0ijrJ36OGgQbLJJ7DRSm4q7iDTK/ffDqlWa+lio1HMXkQZbuhR23x3+/nc48MDY\naUqfpkKKSF4MHQpduqiwFzKduYtIg3zwAXToAK+9BjvuGDtNedBUSBHJuTPOgO99D666KnaS8qGp\nkCKSU6+9FlZ9rK6OnUTqo567iGTEHX7zG7jiCmjZMnYaqY+Ku4hk5KGHYPFiOPPM2EkkE+q5i0i9\nli8Pqz7ecQf8+Mex05QfTYUUkZy44Yaw4bUKe/HIqLibWYWZVZvZLDO7ZD3H7WdmNWZ2XPYiikhM\nixbBNdeEDyke9bZlzKwZMAvoDiwAqoBe7l5dx3HPAcuAO9z9kToeS20ZkSJzzjlh27yhQ2MnKV+5\nmgrZGZjt7nNTTzIK6AnUngx1PvAQsF9DAohI4Zo2DR5/XKs+FqNM2jJtgHlpt+en7vuKmW0LHOPu\ntwANenURkcK0durj5ZfD5pvHTiMNla2LmK4H0nvx6yzwlZWVX32eSCRIJBJZiiAi2TRmDHz4IfTt\nGztJ+UkmkySTySY9RiY99y5ApbtXpG4PANzdh6Qd897aT4FWwJfAOe4+utZjqecuUgRWroQ99oAb\nb4SKithpJFc99ypgFzNrCywEegG90w9w953SQtwJjKld2EWkeAwbBj/4gQp7Mau3uLv7ajPrB4wl\n9OhHuPsMM+sbvuzDa39LDnKKSJ58/DEMHgzjx8dOIk2hK1RF5Gv69Qv/3nRT3BzyP1oVUkSa5K23\n4IEHYMaM2EmkqbT8gIgAYerjhRfC//0ftGoVO400lYq7iADw9NMwdy6cd17sJJINKu4iQk1NOGsf\nOhQ22ih2GskGFXcR4ZZboG1bOOKI2EkkWzRbRqTMffJJWKt93Lhw4ZIUHm2QLSINdv75YTBVUx8L\nl6ZCikiDaOpj6VLPXaRMrV318dJL4bvfjZ1Gsk3FXaRMjRkD8+fDuefGTiK5oLaMSBlasQIuuij0\n2TX1sTTpzF2kDP31r7D77nDYYbGTSK5otoxImVm0KEx5nDABdt01dhrJhKZCiki9zjoLttgCrr02\ndhLJlKZCish6VVWFNWSqa29vLyVHPXeRMrFmDfz61zBoELRsGTuN5JqKu0iZuPdeWL0aTj89dhLJ\nB/XcRcrA4sVhdszDD0OXLrHTSEM1pueuM3eRMjBoEPToocJeTjIq7mZWYWbVZjbLzC6p4+s/NbM3\nzGyKmU0ys67ZjyoijTF7Ntx+O1x9dewkkk/1tmXMrBkwC+gOLACqgF7uXp12TAt3X5r6fC/gH+7e\nro7HUltGJI/cwxrt3bvDxRfHTiONlau2TGdgtrvPdfcaYBTQM/2AtYU9ZRNgTUNCiEhujBkDc+aE\nWTJSXjIp7m2AeWm356fu+xozO8bMZgBjgD7ZiScijbVsGfTvH5YaaN48dhrJt6xdxOTujwGPmVk3\n4ErgJ3UdV1lZ+dXniUSCRCKRrQgikuaaa6BTJzj00NhJpKGSySTJZLJJj5FJz70LUOnuFanbAwB3\n9yHr+Z53gf3c/dNa96vnLpIH778PP/whTJ4c9kaV4parnnsVsIuZtTWz5kAvYHStJ9457fN9gOa1\nC7uI5M+FF4aWjAp7+aq3LePuq82sHzCW8GIwwt1nmFnf8GUfDhxvZqcBK4FlwIm5DC0i6/b00zBt\nGtx3X+wkEpOuUBUpIcuXw557wo03wuGHx04j2aIrVEXK3J//DB06qLCLztxFSsZ778F++2kQtRTp\nzF2kjF1wQbgKVYVdQJt1iJSE0aPDGjIPPxw7iRQKFXeRIvfll2F5gREjdCWq/I967iJF7ne/g4UL\n4Z57YieRXNEeqiJlZto0GDkSpk+PnUQKjQZURYrUmjXQty9ceSW0bh07jRQaFXeRInX77WAGZ58d\nO4kUIvXcRYrQokWw117w/PPhoiUpbY3puau4ixSh3r1h++3DFalS+jSgKlIGnnoKJk0KUx9F1kXF\nXaSILFkC554b+u0tWsROI4VMbRmRItK/P3z2Gdx1V+wkkk9qy4iUsEmTYNQoePPN2EmkGGgqpEgR\nqKmBX/wChg6FVq1ip5FioOIuUgSuvhratIGTT46dRIqFeu4iBW76dDjkEJgyBbbbLnYaiUHruYuU\nmFWroE8fGDxYhV0aJqPibmYVZlZtZrPM7JI6vn6ymb2R+njJzPbKflSR8jN0KGy2mZYYkIarty1j\nZs2AWUB3YAFQBfRy9+q0Y7oAM9z9CzOrACrdvUsdj6W2jEiGqquhWzeoqoIdd4ydRmLKVVumMzDb\n3ee6ew0wCuiZfoC7v+LuX6RuvgK0aUgIEfm61atDO6ayUoVdGieTee5tgHlpt+cTCv66nA08Xd+D\nvvQSzJ0LW20VPrbeGrbZBjbYIINEIiXu2mvh29+GX/0qdhIpVlm9iMnMfgycCXRb1zGVlZVAmAHw\n+ecJmjVL8NFHYZW7zz4Lm/vuvDP84Aew997QsSPssQd861vZTCpSuKZPD8W9qgqaacpDWUomkyST\nySY9RiY99y6EHnpF6vYAwN19SK3jOgAPAxXu/u46Hmu9Pfdly2DOHHj3XZg1C6ZODR/vvgvt28NB\nB8GBB4YPXcghpWjlSvjRj6BfPzjrrNhppFDkZMlfM9sAmEkYUF0ITAJ6u/uMtGN2AF4ATnX3V9bz\nWI0aUF22DCZPhvHj4V//ggkTYKed4Igj4MgjoXNntXOkNFx2WZjPPmZM2IhDBHK4nntqBswNhAHY\nEe5+tZn1JZzBDzez24DjgLmAATXu/o2+fLZmy6xaBRMnhqVPn3wybA587LFw0kmQSKjQS3GaNAmO\nPjq8W/3e92KnkUJStpt1vP8+PPggPPAAzJ8PJ5wQZhp06pT1pxLJiSVLYJ99wn6oJ54YO40UmrIt\n7uneeQfuvTfsCL/llqHIn3IKbL55Tp9WpEnOOSf020eOjJ1ECpGKe5o1a+CFF8JuNc8+C716wfnn\nh4FZkULy6KNw8cWh196yZew0Uoi0tkyaZs3gJz8J61+//Ta0bg3du4f7nn0WdKGsFIIFC8LOSvfe\nq8Iu2VWyZ+51WbEC/vEPGDIENtoIBgyAn/1MA7ASx5o1UFEBBxwQrkQVWRe1ZTK0Zk2YZXPVVfDR\nR2H62cknw4bal0ry6JprQktm/Hj97cn6qbg3kDskkzBwYLhCduDAMJ1SZ/KSaxMnwjHHhOmPbdvG\nTiOFTsW9kdzD4OvAgfDFF+GM/qijdBGJ5MZnn4VpujfcAD171n+8iIp7E7mHds2AAWEa5ZAhsP/+\nsVNJKXGH444LZ+vXXx87jRQLzZZpIrNwxv7GG3DmmeFikhNPDOvdiGTDTTfBvHnhxEEkl1Tc67DB\nBqG4z5wJHTrAvvvC738PixfHTibF7JVX4E9/CldSa5VTyTUV9/Vo0QIuvRSmTQvr1+y2G9x1V5ht\nI9IQ//lPeBd4++1hSWuRXFPPvQEmTQpLsW60EQwbFtaaF6nPqlVw6KFh/GbQoNhppBip555jnTuH\nt9ZnnAGHHRYK/Rdf1PttUuYuvTS0+q64InYSKScq7g3UrBn84hdhSYOaGmjXLixxUMRvSCSHHn00\n/H3cf7+un5D8UlumiSZMgF/+Mqy/PWwY7LJL7ERSKKZPD+sZPfVUGJQXaSy1ZSI44AB4/XXo0QO6\ndIGrrw5n9FLePvoIfvrTMJddhV1i0Jl7Fs2ZE1b4+/BDuO022G+/2IkkhpUrw4t9t24weHDsNFIK\ndIVqAXCH++6Diy4Ki5H96U+w8caxU0m+uEPfvmGtokcfDWM0Ik2Vs7aMmVWYWbWZzTKzS+r4+m5m\nNsHMlpvZhQ0JUGrM4Oc/hzffDHObO3QI69ZIebj++rAo2L33qrBLXPWeuZtZM2AW0B1YAFQBvdy9\nOu2YVkBb4BjgM3f/yzoeq+TP3Gt76qkw4HrooXDttdrur5Q99BD07x8G2XfYIXYaKSW5OnPvDMx2\n97nuXgOMAr62lp27f+zurwOrGvLk5eCII8JZfPPmsOee8MQTsRNJLrz8chhvGTNGhV0KQybFvQ0w\nL+32/NR9kqGWLeHmm8Nb9f79w4bdn3wSO5Vky8yZcPzx4ffbqVPsNCKBuoJ5lEiEFSe32iqcxT/0\nUOxE0lQLF4Z3Z4MHh6uWRQpFJpt7fQCkv9HcLnVfo1SmbRaZSCRIJBKNfaiitPHGcN11cMIJ0KdP\nWCFw2DDYeuvYyaShPvkkbLjep0/4EMmWZDJJMpls0mNkMqC6ATCTMKC6EJgE9Hb3GXUcOxBY4u5D\n1/FYZTeguj7LloWNkUeODAW/d2/t/lQsFi8Oc9kPPjisza7fm+RSzua5m1kFcAOhjTPC3a82s76A\nu/twM2sNvAZsCqwBlgDt3X1JrcdRca9DVVVYP37nneHWW8NSBlK4li0LrZhddw2/LxV2yTVdxFTE\nVqwIy8Heeiv8+c9w+ukqGoVoxYoweLrppmEAVYuBST6ouJeAqVPDWXzr1vC3v4W9NqUwLF8e9j9t\n0SKs8rjRRrETSbnQwmEloGPHsCnIwQeHBaeGDdPOT4Vg2TI45phwxq7CLsVAZ+4FrLoazjorXMZ+\n222w++6xE5WnpUtDYW/VCu6+GzbMZI6ZSBbpzL3E7L47/OtfcNJJcOCBYRGylStjpyovn38OFRWh\nTabCLsVExb3ANWsWtvObPBlefRX22Sds9Se5t2ABHHRQuOr0rrtU2KW4qLgXie23D+uWXHZZGNQ7\n99xwVim5MWsWdO0alm2+/nqt8CjFR3+yRcQstGjefjvcbt8+DO5pGCO7JkwIA9qXXQYDBmhKqhQn\nDagWsYkTw8YQrVvDjTdqwDUbRo6E3/0utGEOPzx2GpFAA6plZv/9w/6tRx4ZBlwHDIAlS+r/Pvmm\n1avh4ovDhWQvvqjCLsVPxb3IbbRRWEZ42jT44ANo1y5s86c3SJn75BM4+miYMiUMWrdrFzuRSNOp\nuJeI730P7rkn9OD/8hc44IBwMZSs38svh9kw7dvDM8/AllvGTiSSHSruJaZbt1DUzzknXHhz6qkw\nd27sVIVnzRq46qqwTszNN4ctEHXVqZQSFfcS1KxZWJ9m5kzYcccwN/6ii7T701rvvx/2tH3yybAi\n51FHxU4kkn0q7iVs003hiivCHq5Ll4bZNIMGhbXIy9GaNXDLLWHNnh49IJkM1w+IlCJNhSwjM2eG\nYv/cc3DhheHK1002iZ0qP2bPDq2qZcvgzjs1aCrFRVMhZb122w3+/vdwxjp1atgc5Mor4dNPYyfL\nnS++gN/+NkwbPfroMICqwi7lQMW9DLVvD6NGwbhx8O67sMsu4Uz+3/+OnSx7Vq2C228PrahPPw2t\nqQsv1OYaUj5U3MtY+/ahRTFtWhiE7dgxzB755z+Ld558TU24yrRduzA1dMwYGDECttkmdjKR/FLP\nXb6yeHEoiMOGheJ+zjlh4aytt46drH5Ll4Zt74YMCbtXXX45JBKxU4lkR8567mZWYWbVZjbLzC5Z\nxzF/NbPZZjbVzDo2JIQUhk03hV/9KrQwbr45LDO8666hV/3gg2EwstBUV4crdNeumjlyZHjnocIu\n5a7e4m5mzYCbgMOAPYDeZrZ7rWMOB3Z29x8AfYFbc5A1J5LJZOwIdYqZyywUx7vvhvnz4YQTYPhw\naNUqybHHhkW1Ys6ZnzcvLMPbtSsccECSjTcOL0RjxoQ1dmIrxL+pQswEhZmrEDM1RiZn7p2B2e4+\n191rgFFAz1rH9ATuBnD3V4HNzKx1VpPmSKH+Igsl1yabwGmnhemT/folOe44GD06XBy1775hJsqT\nT4ZZKbmyfHmY4VNZGWa9dOwYxgkuvRTOOy/JoEGFtZF4ofzu0hViJijMXIWYqTEy2VumDTAv7fZ8\nQsFf3zEfpO5b1KR0UlC+852wnMGpp8KKFeHqznHjYOhQOPHEsL5Np07ho107+P73w8fmm2e2Jvrq\n1fDRR2FO+ptvwltvwRtvhAW99twzvJu44orw79qlAl59NXf/vSLFTBuHSaN861thHZtu3cKmFqtW\nhd2Lpk4NxXjkyHCZ/5w54fgtt4SWLWGzzaBFi3C1aE1N+L4lS+DDD0OrZ/PNw9TMPfYIBf2YY+BH\nPwrjASKSuXpny5hZF6DS3StStwcA7u5D0o65FRjn7g+kblcDB7v7olqPpakyIiKN0NDZMpmcuVcB\nu5hZW2Ah0AvoXeuY0cB5wAOpF4PPaxf2xoQTEZHGqbe4u/tqM+sHjCUMwI5w9xlm1jd82Ye7+1Nm\ndoSZvQN8CZyZ29giIrI+eb2ISURE8iNvyw9kciFUHjKMMLNFZjYt7b4tzGysmc00s2fNbLM8Z9rO\nzP5pZm+Z2XQz+3XsXGb2LTN71cympDINjJ0pLVszM5tsZqMLKNP7ZvZG6uc1qRBymdlmZvagmc1I\n/W39KPLf1K6pn8/k1L9fmNmvC+Dn9Bsze9PMppnZ382seQFkuiD1/12T6kFeinsmF0LlyZ2pDOkG\nAM+7+27AP4Hf5znTKuBCd98D2B84L/WziZbL3VcAP3b3TkBH4HAz6xwzU5oLgLfTbhdCpjVAwt07\nufvaacKxc90APOXu7YC9geqYmdx9Vurnsw/wQ0L79tGYmcxsW+B8YB9370BoU/eOnGkP4CxgX8L/\ne0eZ2c6NyuTuOf8AugBPp90eAFySj+euI0tbYFra7WqgderzbYDqGLnS8jwG9CiUXEAL4DVgv9iZ\ngO2A54AEMLpQfn/AHOC7te6LlgtoCbxbx/3Rf1ap5z4U+FfsTMC2wFxgC0JhHx37/z3gZ8Btabcv\nBX4LzGhopny1Zeq6EKpNnp67Plt7amaPu38IRFsmy8y+T3i1foXwi4yWK9X+mAJ8CDzn7lWxMwHX\nEf7Q0weKYmcilec5M6sys7MLINeOwMdmdmeqDTLczFpEzpTuJOC+1OfRMrn7AmAo8G/ChZdfuPvz\nMTMBbwIHptowLYAjgO0bk0lL/n5TlBFmM9sEeAi4wN2X1JEjr7ncfY2Htsx2QOfU28VomczsSGCR\nu08F1jefhGrSAAACJElEQVSlNsbvr6uHdsMRhLbagXXkyGeuDYF9gGGpXF8S3i1H/ZsCMLONgJ8C\nD64jQz7/pjYnLJ3SlnAWv7GZ/TxmJnevBoYQ3qE+BUwBVtd1aH2Pla/i/gGwQ9rt7VL3FYJFa9fB\nMbNtgP/kO4CZbUgo7Pe4++OFkgvA3f8LJIGKyJm6Aj81s/eA+4FDzOwe4MPYPyd3X5j69yNCW60z\ncX9W84F57v5a6vbDhGJfCH9ThwOvu/vHqdsxM/UA3nP3T919NWEM4IDImXD3O919X3dPAJ8DMxuT\nKV/F/asLocysOeFCqNF5eu7ajK+f+Y0Gzkh9fjrweO1vyIM7gLfd/Ya0+6LlMrNWa0fjzew7wE8I\nPb9omdz9D+6+g7vvRPj7+ae7nwqMiZUJwMxapN51YWYbE/rJ04n7s1oEzDOzXVN3dQfeipkpTW/C\ni/NaMTP9G+hiZt82MyP8nN6OnAkz2yr17w7AsYQWVsMz5XGgoILwCjQbGJCv562V4T5gAbCC8Is9\nkzCY8nwq21hg8zxn6kp42zWV8BZscupntWWsXMBeqRxTgWnA/6Xuj5apVr6D+d+AatRMhP722t/d\n9LV/2wWQa2/CSdVU4BFgswLI1AL4CNg07b7YmQYSTlymAXcBGxVApvGE3vsUwiysRv2cdBGTiEgJ\n0oCqiEgJUnEXESlBKu4iIiVIxV1EpASpuIuIlCAVdxGREqTiLiJSglTcRURK0P8D3O3xHU6qrqgA\nAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x111338160>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Solution:\n",
    "# For the first plot, we can show the joint probability for two angles is given by:\n",
    "\n",
    "plt.plot(rad2deg(rads),(sqrt(0.2)*cos(-rads)*cos(rads) + sqrt(0.8)*sin(-rads)*sin(rads))**2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Challenge:\n",
    "If we change the state to $\\big|\\psi\\big\\rangle = \\sqrt{0.8} \\big|H,H\\big\\rangle + \\sqrt{0.2} \\big|V,V\\big\\rangle$, the two angles that work for this state."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Solution\n",
    "\n",
    "psi3=sqrt(0.8)*tensor(H,H) + sqrt(0.2)*tensor(V,V)\n",
    "\n",
    "out = []\n",
    "for r in rads:\n",
    "    out.append(expect(tensor(P(-r),P(r)),psi3))\n",
    "\n",
    "plt.plot(angles,out,\".\") # plot in units of pi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Solution\n",
    "out = []\n",
    "for r in rads:\n",
    "    out.append(expect(tensor(P(r),P(deg2rad(55))),psi3) / expect(tensor(P(r),qeye(2)),psi3))\n",
    "\n",
    "plt.plot(angles,out,\".\")"
   ]
  }
 ],
 "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.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}