Exercise 5 - Cosmic rays solutions.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<font color='darkred'>\n",
    "    \n",
    "__Rubric__ (20 points total):\n",
    "\n",
    "- Part (1) __10 points__\n",
    "     - Writing code to generating samples $s$ correctly (with correct inverse CDF) __5 points__\n",
    "     - Converting to $\\theta$ and plotting a reasonable histogram __5 points__\n",
    "      - (Note: $N=10^5$ is OK instead of $N=10^6$ for sample size)\n",
    "\n",
    "- Part (2) __10 points__\n",
    "     - Writing code to generate samples for $\\theta$ __5 points__\n",
    "     - Plotting sensible histograms, plotting everything together, and showing agreement between parts (a) and (b) __5 points__\n",
    "</font>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Sean Tulin__\n",
    "<br>\n",
    "Date: Jan. 21, 2021\n",
    "<br>\n",
    "PHYS 2030 W22"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# <center><font color=#46769B>Exercise 5: Cosmic rays</font></center>\n",
    "\n",
    "## <font color=#46769B>Introduction:</font>\n",
    "\n",
    "Our goals for this notebook are:\n",
    "- Make a simulation for the angle of cosmic rays when hitting the Earth.\n",
    "\n",
    "Required reading:\n",
    "- *Lesson 3: Sampling from arbitrary distributions*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <font color=#46769B>Cosmic rays:</font>\n",
    "\n",
    "Cosmic rays are high-energy protons and heavier elements accelerated in stars and other astronomical environments that ultimately smash into Earth's atmosphere. These collisions produce showers of particles.\n",
    "\n",
    "<div>\n",
    "<img src=\"https://cdn.mos.cms.futurecdn.net/JZG8KC9eZovtTaXkvbbdiX-1024-80.jpg.webp\" width=\"400\">\n",
    "</div>\n",
    "\n",
    "Figure credit: [space.com](https://www.space.com/32644-cosmic-rays.html)\n",
    "\n",
    "## <font color=#46769B>Goal:</font>\n",
    "\n",
    "Suppose cosmic rays are isotropic and homogenous in our local galactic neighborhood. They are likely to hit the Earth from any direction and at any position. When a cosmic ray hits the Earth, it will have an angle $\\theta$ with respect to the atmosphere, shown in the Figure here.\n",
    "\n",
    "<div>\n",
    "<img src=\"https://github.com/PHYS-2030-Computational-Methods/Lecture-notes/raw/main/figures/cosmic_rays.png\" width=\"500\">\n",
    "</div>\n",
    "\n",
    "The goal here is figure out what is the probability distribution function for $\\theta$. (All values of $\\theta$ are *not* equally likely.) We will solve this problem two ways and show that they agree."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# <font color=#46769B>Exercise</font>\n",
    "\n",
    "The Earth is (approximately) a sphere of radius $R$. Since the cross sectional area of a sphere is a circle, to a cosmic ray, the Earth looks like a dart board.\n",
    "\n",
    "<div>\n",
    "<img src=\"https://upload.wikimedia.org/wikipedia/commons/f/fb/Darts_in_a_dartboard.jpg\" width=\"400\">\n",
    "</div>\n",
    "\n",
    "Here is another picture to have in mind. Without loss of generality, we orient our coordinate frame so the cosmic ray is pointing down.\n",
    "\n",
    "<div>\n",
    "<img src=\"https://github.com/PHYS-2030-Computational-Methods/Lecture-notes/raw/main/figures/cosmic-ray-fig-2.png\" width=\"800\">\n",
    "</div>\n",
    "\n",
    "The cosmic ray hits a sphere, the Earth __(a)__. In projection, this looks like a circle __(b)__. Let's suppose the radius $s$ is the location where the cosmic ray hits the circular cross section of the Earth. (Note: $s \\le R$, where $R$ is the radius of the earth.) Once we know $s$, we can figure out the angle the cosmic ray makes with some trigonometry __(c)__:\n",
    "$$\\sin\\theta = s/R \\, .$$\n",
    "\n",
    "What is the probability of landing on a particular value of $s$? We need to know the PDF $P(s)$.\n",
    "\n",
    "Note $P(s) ds$ is the probability of landing between $s$ and $s+ds$. This is\n",
    "$$ P(s) \\, ds = \\frac{2 \\pi s \\, ds}{\\pi R^2}$$\n",
    "where the numerator is the area of a ring of radius $s$ and thickness $ds$, and the denominator is the area $\\pi R^2$ of the entire circle.\n",
    "The PDF is therefore\n",
    "$$P(s) = \\left\\{ \\begin{array}{cc} \\frac{2 s}{R^2} & {\\rm for} \\; 0 \\le s \\le R \\\\ 0 & {\\rm otherwise} \\end{array} \\right. \\, .\\qquad {\\rm (1)}$$ \n",
    "\n",
    "\n",
    "## <font color=#46769B>Part (a)</font> \n",
    "\n",
    "Perform a simulation generating $N=10^6$ random samples for $\\theta$ by sampling $s$ from $P(s)$ in Eq. (1). Here are the steps:\n",
    "- Determine (analytically) the CDF $C(s)$ corresponding to the PDF $P(s)$. \n",
    "- Analytically solve the equation $r = C(s)$ to determine the inverse CDF, $C^{-1}$.\n",
    "- Write a function `sample_s(num)` that will draw `num` random samples for $s$ from $P(s)$. (You compute $N$ random numbers $r$ in the domain $[0,1]$ and set $s=C^{-1}(r)$).\n",
    "- For each $s$ in your sample, compute $\\theta = \\sin^{-1}(s/R)$ to obtain a sample for $\\theta$.\n",
    "\n",
    "Make a histogram plot for your $\\theta$ samples.\n",
    "\n",
    "### <font color=#46769B>Hint:</font> \n",
    "\n",
    "You can work in units of \"Earth radii\" setting $R=1$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <font color='darkred'>Solutions:</font>\n",
    "\n",
    "<font color='darkred'>\n",
    "    \n",
    "Some results (you don't need to show this explicitly, but it is needed for your code).\n",
    "    \n",
    "- Step 1: Compute the CDF $C(s)$ by integrating the PDF $P(z)$.\n",
    "\n",
    "$$C(s) = \\int_{-\\infty}^s dz \\, P(z) = \\int_{-\\infty}^{{\\rm min}(s,R)} dz \\, \\frac{2 z}{R^2} \n",
    "= \\left\\{\\begin{array}{cc} \\frac{s^2}{R^2} & {\\rm for} \\; 0 \\le s \\le R \\\\ 1 & {\\rm for} \\; s > R \\\\ 0 & {\\rm otherwise} \\end{array} \\right.$$\n",
    "\n",
    "- Step 2: Compute the inverse CDF $C^{-1}$ by solving the equation $r = C(s) = s^2/R^2$, where as usual $r$ will be a random number in domain $[0,1]$. Note we only care about the range $0 \\le s \\le R$.\n",
    "    \n",
    "$$s = C^{-1}(r) = R \\sqrt{r} $$\n",
    "\n",
    "</font>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.9725639 , 0.9271037 , 0.7770798 , 0.92365275, 0.30635974,\n",
       "       0.51932186, 0.93716366, 0.52691546, 0.54666664, 0.48104851])"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# Your code here\n",
    "\n",
    "R = 1\n",
    "\n",
    "def sample_s(num):\n",
    "    \n",
    "    r = np.random.rand(num)\n",
    "    s = R*np.sqrt(r)\n",
    "    return s\n",
    "\n",
    "# Try it out\n",
    "sample_s(10)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYMAAAD4CAYAAAAO9oqkAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAASgUlEQVR4nO3df6zddX3H8edLEMUpAlIca3FlsSZWMhTuoNNl/sBhAWNZgrNMpZpmTRg6t5lN3JKx+SOpWzYdCWo6aShmWpk/RqNg1yGMbVLgVpGfc1yRwR3GFgsMJcKK7/1xPnXHcm7v97a95562z0dycr7f9/dzznmfc++5r/v9cb4nVYUk6eD2jLluQJI09wwDSZJhIEkyDCRJGAaSJODQuW5gTx1zzDG1cOHCuW5DkvYbW7Zseaiq5g1att+GwcKFCxkfH5/rNiRpv5Hkv6Za5mYiSZJhIEkyDCRJGAaSJAwDSRKGgSQJw0CShGEgScIwkCSxH38CWRq2hRd9pdO4+1afPcudSPueYSDtY4aG9keGgQ56Xf94Swcy9xlIklwzkOaKm5M0SjqtGSS5L8ntSW5NMt5qRyfZlOSedn1UqyfJJUkmktyW5OS++1nRxt+TZEVf/ZR2/xPtttnXT1SSNLWZrBm8tqoe6pu/CLi2qlYnuajNvw84E1jULqcBnwBOS3I0cDEwBhSwJcmGqnq4jVkFbAauBpYC1+zVM5MOEK5BaBj2ZjPRMuA1bXodcD29MFgGXFFVBWxOcmSS49rYTVW1HSDJJmBpkuuBI6rqxla/AjgHw0B7yR3DUndddyAX8E9JtiRZ1WovrKrvAbTrY1t9PvBA320nW2139ckB9adJsirJeJLxbdu2dWxdkjSdrmsGr6qqB5McC2xK8h+7GTtoe3/tQf3pxao1wBqAsbGxgWMkSTPXac2gqh5s11uBLwGnAt9vm39o11vb8Eng+L6bLwAenKa+YEBdkjQk04ZBkp9L8ryd08AZwB3ABmDnEUErgKva9Abg/HZU0RLg0bYZaSNwRpKj2pFHZwAb27LHkixpRxGd33dfkqQh6LKZ6IXAl9rRnocCn6mqrya5BbgyyUrgfuDNbfzVwFnABPA48E6Aqtqe5IPALW3cB3buTAYuAC4HDqe349idx5I0RNOGQVXdC5w0oP4D4PQB9QIunOK+1gJrB9THgRM79CtJmgV+Aln7FQ8XndpMXhs/k6BdeW4iSZJhIEkyDCRJGAaSJAwDSRKGgSQJDy3ViPCQUWluGQbSQcjvSNCu3EwkSTIMJEmGgSQJw0CShGEgScIwkCRhGEiSMAwkSfihM80yP1ks7R9cM5AkuWYgaWqetuLg4ZqBJMkwkCQZBpIkDANJEoaBJAnDQJKEYSBJwjCQJGEYSJIwDCRJGAaSJGZwbqIkhwDjwH9X1RuTnACsB44GvgG8vaqeTPIs4ArgFOAHwFuq6r52H+8HVgJPAb9XVRtbfSnwt8AhwKeqavU+en6ShsBzGO3/ZnKiuvcAdwNHtPmPAB+tqvVJPknvj/wn2vXDVfXiJMvbuLckWQwsB14G/ALwz0le0u7rUuA3gEngliQbququvXxumkWemlo6sHTaTJRkAXA28Kk2H+B1wOfbkHXAOW16WZunLT+9jV8GrK+qJ6rqu8AEcGq7TFTVvVX1JL21jWV7+8QkSd113WfwMeCPgZ+0+RcAj1TVjjY/Ccxv0/OBBwDa8kfb+J/Wd7nNVPWnSbIqyXiS8W3btnVsXZI0nWnDIMkbga1VtaW/PGBoTbNspvWnF6vWVNVYVY3NmzdvN11Lkmaiyz6DVwFvSnIW8Gx6+ww+BhyZ5ND23/8C4ME2fhI4HphMcijwfGB7X32n/ttMVZckDcG0awZV9f6qWlBVC+ntAP5aVb0VuA44tw1bAVzVpje0edryr1VVtfryJM9qRyItAm4GbgEWJTkhyWHtMTbsk2cnSepkb7728n3A+iQfAr4JXNbqlwGfTjJBb41gOUBV3ZnkSuAuYAdwYVU9BZDkXcBGeoeWrq2qO/eiL0nSDM0oDKrqeuD6Nn0vvSOBdh3zY+DNU9z+w8CHB9SvBq6eSS+SpH3HTyBLkgwDSZJhIEnCMJAkYRhIkjAMJEns3ecMdADybKSaTZ7qenS5ZiBJMgwkSYaBJAnDQJKEYSBJwjCQJGEYSJIwDCRJGAaSJAwDSRKGgSQJw0CShGEgScIwkCThKawljSBPdT18hsFBwu8pkLQ7biaSJBkGkiTDQJKEYSBJwjCQJGEYSJIwDCRJGAaSJDqEQZJnJ7k5ybeS3JnkL1r9hCQ3JbknyeeSHNbqz2rzE235wr77en+rfzvJG/rqS1ttIslF+/5pSpJ2p8uawRPA66rqJODlwNIkS4CPAB+tqkXAw8DKNn4l8HBVvRj4aBtHksXAcuBlwFLg40kOSXIIcClwJrAYOK+NlSQNybRhUD0/bLPPbJcCXgd8vtXXAee06WVtnrb89CRp9fVV9URVfReYAE5tl4mqureqngTWt7GSpCHptM+g/Qd/K7AV2AR8B3ikqna0IZPA/DY9H3gAoC1/FHhBf32X20xVH9THqiTjSca3bdvWpXVJUgedwqCqnqqqlwML6P0n/9JBw9p1plg20/qgPtZU1VhVjc2bN2/6xiVJnczoaKKqegS4HlgCHJlk51lPFwAPtulJ4HiAtvz5wPb++i63maouSRqSLkcTzUtyZJs+HHg9cDdwHXBuG7YCuKpNb2jztOVfq6pq9eXtaKMTgEXAzcAtwKJ2dNJh9HYyb9gXT06S1E2X7zM4DljXjvp5BnBlVX05yV3A+iQfAr4JXNbGXwZ8OskEvTWC5QBVdWeSK4G7gB3AhVX1FECSdwEbgUOAtVV15z57hpKkaaX3T/v+Z2xsrMbHx+e6jTnnl9boYOY3nc1Mki1VNTZomd90Jmm/5ddj7juejkKSZBhIkgwDSRKGgSQJw0CShGEgScIwkCRhGEiSMAwkSRgGkiQMA0kShoEkCcNAkoRhIEnCMJAkYRhIkjAMJEn4TWcjya+ylPatmbynDtZvRXPNQJJkGEiSDANJEoaBJAnDQJKEYSBJwjCQJGEYSJIwDCRJGAaSJAwDSRIdwiDJ8UmuS3J3kjuTvKfVj06yKck97fqoVk+SS5JMJLktycl997Wijb8nyYq++ilJbm+3uSRJZuPJSpIG67JmsAN4b1W9FFgCXJhkMXARcG1VLQKubfMAZwKL2mUV8AnohQdwMXAacCpw8c4AaWNW9d1u6d4/NUlSV9OGQVV9r6q+0aYfA+4G5gPLgHVt2DrgnDa9DLiiejYDRyY5DngDsKmqtlfVw8AmYGlbdkRV3VhVBVzRd1+SpCGY0T6DJAuBVwA3AS+squ9BLzCAY9uw+cADfTebbLXd1ScH1Ac9/qok40nGt23bNpPWJUm70TkMkjwX+ALw+1X1P7sbOqBWe1B/erFqTVWNVdXYvHnzpmtZktRRpzBI8kx6QfD3VfXFVv5+28RDu97a6pPA8X03XwA8OE19wYC6JGlIuhxNFOAy4O6q+pu+RRuAnUcErQCu6quf344qWgI82jYjbQTOSHJU23F8BrCxLXssyZL2WOf33ZckaQi6fO3lq4C3A7cnubXV/gRYDVyZZCVwP/Dmtuxq4CxgAngceCdAVW1P8kHgljbuA1W1vU1fAFwOHA5c0y4HHL/OUtKomjYMqurfGLxdH+D0AeMLuHCK+1oLrB1QHwdOnK4XSdLs6LJmIEkHja5r8PetPnuWOxkuT0chSTIMJEmGgSQJw0CShGEgScIwkCRhGEiSMAwkSRgGkiQMA0kShoEkCcNAkoRhIEnCMJAkYRhIkjAMJEkYBpIkDANJEn7t5T7hF91L2t8ZBpK0Bw6070p2M5EkyTCQJBkGkiQMA0kShoEkCcNAkoRhIEnCMJAkYRhIkjAMJEl0CIMka5NsTXJHX+3oJJuS3NOuj2r1JLkkyUSS25Kc3HebFW38PUlW9NVPSXJ7u80lSbKvn6Qkafe6rBlcDizdpXYRcG1VLQKubfMAZwKL2mUV8AnohQdwMXAacCpw8c4AaWNW9d1u18eSJM2yacOgqm4Atu9SXgasa9PrgHP66ldUz2bgyCTHAW8ANlXV9qp6GNgELG3LjqiqG6uqgCv67kuSNCR7us/ghVX1PYB2fWyrzwce6Bs32Wq7q08OqA+UZFWS8STj27Zt28PWJUm72tc7kAdt7689qA9UVWuqaqyqxubNm7eHLUqSdrWnYfD9tomHdr211SeB4/vGLQAenKa+YEBdkjREexoGG4CdRwStAK7qq5/fjipaAjzaNiNtBM5IclTbcXwGsLEteyzJknYU0fl99yVJGpJpv+ksyWeB1wDHJJmkd1TQauDKJCuB+4E3t+FXA2cBE8DjwDsBqmp7kg8Ct7RxH6iqnTulL6B3xNLhwDXtIkkaomnDoKrOm2LR6QPGFnDhFPezFlg7oD4OnDhdH5Kk2eMnkCVJ068ZHMy6fuG1JE2l69+R+1afPcud7J5rBpIkw0CSZBhIkjAMJEkYBpIkDANJEoaBJAnDQJKEYSBJwjCQJGEYSJIwDCRJGAaSJAwDSRKGgSQJw0CShGEgScIwkCRhGEiSMAwkSRgGkiTg0LluYC4svOgrc92CJP2Mrn+X7lt99qw8vmsGkiTDQJJkGEiSMAwkSRgGkiQMA0kShoEkiREKgyRLk3w7yUSSi+a6H0k6mIxEGCQ5BLgUOBNYDJyXZPHcdiVJB4+RCAPgVGCiqu6tqieB9cCyOe5Jkg4ao3I6ivnAA33zk8Bpuw5KsgpY1WZ/mOTbe/h4xwAP7eFtZ9Oo9gX2tqfsbeZGtS8Ygd7ykSkXdentF6daMCphkAG1elqhag2wZq8fLBmvqrG9vZ99bVT7AnvbU/Y2c6PaFxzYvY3KZqJJ4Pi++QXAg3PUiyQddEYlDG4BFiU5IclhwHJgwxz3JEkHjZHYTFRVO5K8C9gIHAKsrao7Z/Eh93pT0ywZ1b7A3vaUvc3cqPYFB3BvqXrapnlJ0kFmVDYTSZLmkGEgSTpww2C601skeVaSz7XlNyVZOEK9/WGSu5LcluTaJFMeGzzs3vrGnZukkgztMLsuvSX5rfba3ZnkM6PSW5IXJbkuyTfbz/WsIfW1NsnWJHdMsTxJLml935bk5GH01bG3t7aebkvy9SQnjUJffeN+JclTSc4dRl9de0vymiS3tvfAv3S+86o64C70dkJ/B/gl4DDgW8DiXcb8LvDJNr0c+NwI9fZa4Dlt+oJR6q2Nex5wA7AZGBuV3oBFwDeBo9r8sSPU2xrggja9GLhvSL39OnAycMcUy88CrqH3WZ8lwE3D6Ktjb6/s+1meOazepuur72f+NeBq4NwRes2OBO4CXtTmO78HDtQ1gy6nt1gGrGvTnwdOTzLow29D762qrquqx9vsZnqfuxiGrqcF+SDwl8CPh9RX195+B7i0qh4GqKqtI9RbAUe06eczpM/RVNUNwPbdDFkGXFE9m4Ejkxw3Cr1V1dd3/iwZ4vugw2sG8G7gC8CwfseATr39NvDFqrq/je/c34EaBoNObzF/qjFVtQN4FHjBiPTWbyW9/9yGYdrekrwCOL6qvjyknnbq8rq9BHhJkn9PsjnJ0hHq7c+BtyWZpPff5LuH09q0Zvr7OFeG+T7YrSTzgd8EPjnXvQzwEuCoJNcn2ZLk/K43HInPGcyCLqe36HQKjFnQ+XGTvA0YA149qx31PeSA2k97S/IM4KPAO4bUT78ur9uh9DYVvYbef5H/muTEqnpkBHo7D7i8qv46ya8Cn269/WSWe5vOXL0POkvyWnph8Gtz3UvzMeB9VfXUcDYmzMihwCnA6cDhwI1JNlfVf3a54YGoy+ktdo6ZTHIovVX36VYNh9UbSV4P/Cnw6qp6Ygh9dentecCJwPXtTfDzwIYkb6qq8TnubeeYzVX1v8B324kMF9H7hPtc97YSWApQVTcmeTa9E4sNdTPDACN9Kpgkvwx8Cjizqn4w1/00Y8D69h44BjgryY6q+se5bQvo/TwfqqofAT9KcgNwEjBtGAxlp8ewL/RC7l7gBP5/h97LdhlzIT+7A/nKEertFfR2SC4atddtl/HXM7wdyF1et6XAujZ9DL3NHy8Ykd6uAd7Rpl9K7w9uhvTaLWTqHY5n87M7kG8e8u/c7np7ETABvHKYPU3X1y7jLmeIO5A7vGYvBa5tv5PPAe4ATuxyvwfkmkFNcXqLJB8AxqtqA3AZvVX1CXprBMtHqLe/Ap4L/EP77+P+qnrTiPQ2Jzr2thE4I8ldwFPAH9UQ/pvs2Nt7gb9L8gf0NsO8o9q7dzYl+Sy9zWbHtP0VFwPPbH1/kt7+i7Po/dF9HHjnbPc0g97+jN5+vI+398GOGsIZQzv0NWem662q7k7yVeA24CfAp6pqt4fI/vS+h/D7KEkacQfq0USSpBkwDCRJhoEkyTCQJGEYSJIwDCRJGAaSJOD/AIb9ToCl2v1GAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Convert to theta\n",
    "\n",
    "num = 10**6\n",
    "s_arr = sample_s(num)\n",
    "th_arr = np.arcsin(s_arr/R)\n",
    "\n",
    "plt.hist(th_arr,bins=30)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <font color=#46769B>Chain rule for PDFs</font>\n",
    "\n",
    "Next, we will obtain $P(\\theta)$ directly and sample from that. First, let me clarify the notation a bit. Every variable has *its own* PDF. So, I should write $P_s(s)$ for the PDF for $s$ (given above) and $P_\\theta(\\theta)$ will be the PDF for $\\theta$ that I will derive now.\n",
    "\n",
    "This is just to clarify that $P_s$ is a completely different function than $P_\\theta$. They are related by the chain rule.\n",
    "\n",
    "The key idea is that probabilities (not the PDFs) should be invariant no matter what variables I use. If a cosmic ray hits a certain patch of sky 5\\% of the time, it doesn't matter if I describe that patch in terms of $s$ or $\\theta$. This idea is written formally as \n",
    "$$P_s(s) \\, ds = P_\\theta(\\theta) \\, d\\theta \\, ,$$\n",
    "which says that the probability of being within the window $[s,s+ds]$ is equal to the probability of being within the window $[\\theta,\\theta+d\\theta]$.\n",
    "Rearranging this equation, we have\n",
    "$$P_\\theta(\\theta) = P_s(s) \\, \\frac{ds}{d\\theta}\\, .$$\n",
    "Plugging in $P_s(s)$, and using the facts that $s=R \\sin\\theta$ and $ds/d\\theta = R \\cos\\theta$, we have\n",
    "$$P_\\theta(\\theta) = \\left\\{ \\begin{array}{cc} \\sin 2\\theta & {\\rm for} \\; 0 \\le \\theta \\le \\frac{\\pi}{2} \\\\ \n",
    "0 & {\\rm otherwise} \\end{array} \\right. \\, . \\qquad {\\rm (2)}$$\n",
    "\n",
    "## <font color=#46769B>Part (b)</font>\n",
    "\n",
    "Perform a simulation generating $N=10^6$ random samples for $\\theta$ by sampling $\\theta$ from $P_\\theta(\\theta)$ in Eq. (2). Here are the steps:\n",
    "- Determine (analytically) the CDF $C_\\theta(\\theta)$ corresponding to the PDF $P_\\theta(\\theta)$. \n",
    "- Analytically solve the equation $r = C_\\theta(\\theta)$ to determine the inverse CDF, $C_\\theta^{-1}$.\n",
    "- Write a function `sample_theta(num)` that will draw `num` random samples for $\\theta$ from $P_\\theta(\\theta)$.\n",
    "\n",
    "Make a histogram plot for your $\\theta$ samples.\n",
    "\n",
    "Compare your histogram for parts (a) and (b) and a plot of the curve for $P_\\theta(\\theta)$. (Make one figure with everything in it.) Everything should agree!\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <font color='darkred'>Solutions:</font>\n",
    "\n",
    "<font color='darkred'>\n",
    "    \n",
    "Some results (you don't need to show this explicitly, but it is needed for your code).\n",
    "    \n",
    "- Step 1: Calculate PDF for $\\theta$: $P_\\theta(\\theta) = P_s(s) \\, \\frac{ds}{d\\theta}$. Note $s = R \\sin\\theta$ and $\\frac{ds}{d\\theta} = R \\, \\cos\\theta$. Using these relations, we have\n",
    "$$P_\\theta(\\theta) = \\frac{2s}{R^2} \\, R \\cos\\theta = 2 \\sin\\theta \\cos\\theta = \\sin 2\\theta $$\n",
    "Note: this only applies in the domain $0 \\le s \\le R$. The corresponding range for $\\theta$ is $0 \\le \\theta \\le \\frac{\\pi}{2}$. You can also look at the Figures above and convince yourself that $\\theta$ lies in the range of $0^\\circ$ to $90^\\circ$.\n",
    "\n",
    "- Step 2: Calculate the CDF $C_\\theta(\\theta)$:\n",
    "$$C_\\theta(\\theta) = \\int_{-\\infty}^\\theta P_\\theta(z) \\, dz \n",
    "= \\int_0^{{\\rm min}(\\theta,\\tfrac{\\pi}{2})} dz \\, \\sin 2z \n",
    "= \\left\\{ \\begin{array}{cc} \\sin^2 \\theta & {\\rm for} \\; 0 \\le \\theta \\le \\frac{\\pi}{2} \\\\\n",
    "0 & {\\rm for} \\; \\theta < 0 \\\\\n",
    "1 & {\\rm for} \\; \\theta > \\frac{\\pi}{2} \\end{array} \\right.$$    \n",
    "    \n",
    "- Step 3: Compute the inverse CDF $C^{-1}$ by solving the equation $r = C_\\theta(\\theta) = \\sin^2\\theta$. We get\n",
    "    \n",
    "$$\\theta = C^{-1}(r) = \\sin^{-1}\\left(\\sqrt{r}\\right) $$\n",
    "\n",
    "</font>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.73021651, 0.70133472, 0.97062466, 0.47018543, 0.57870708,\n",
       "       0.27796208, 0.72095426, 0.24257719, 1.27502924, 0.38280859])"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def sample_th(num):\n",
    "    \n",
    "    r = np.random.rand(num)\n",
    "    th = np.arcsin(np.sqrt(r))\n",
    "    return th\n",
    "\n",
    "# Try it out\n",
    "sample_th(10)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3dd3jV5fnH8fd9VvYAEmaAgAwhgIARsIqCoIJatRUqqDhbWlvb/qodUi1WrXXUqrVFEcHZKo46cFKQoewle4awwkpIIDs56/n9kYAhJOQEkvM9Oblf18V1nfGc8/3kkNx58h33I8YYlFJKNX02qwMopZRqGFrQlVIqTGhBV0qpMKEFXSmlwoQWdKWUChMOqzaclJRkUlNTrdq8Uko1SatXrz5ijEmu6TnLCnpqaiqrVq2yavNKKdUkicie2p7TXS5KKRUmtKArpVSY0IKulFJhwrJ96Eqp5s3j8ZCVlUVZWZnVUUJSZGQkKSkpOJ3OgF+jBV0pZYmsrCzi4uJITU1FRKyOE1KMMeTm5pKVlUWXLl0Cfp3uclFKWaKsrIxWrVppMa+BiNCqVat6//WiBV0pZRkt5rU7k8+mzoIuIq+ISLaIbKzleRGR50UkQ0TWi8jAeqdQSil11gKZob8GjDrN86OB7pX/JgIvnn0spZRqfHa7nf79+9OnTx/Gjh1LSUnJSY+npaVx3nnn8cwzz+D3+wFYsGABCQkJ9O/fn/79+zNy5Egrv4ST1HlQ1BjztYiknmbIdcAbpmKljGUikigi7YwxBxsoo1INY/7jAOT4StntKSDLW8Q+bxF5vjIK/G7y/W48+CEhBQCXzUV8RDzxrnhaRbWiY1xHUmJTSE1IJSkqycqvRDWQqKgo1q5dC8DNN9/M1KlTuffee096PDs7m5tuuon8/HwefvhhAIYOHcqnn35qWe7aNMRZLh2AfVXuZ1U+dkpBF5GJVMzi6dSpUwNsWqmaTfpgAwB+4+OYfydHvBtxlCxgn62AQpv7xDibEWKNkyicRBsHDmy05AAA5cbHTr+bAr+bPH85fr5b3au1PYo0Vyv6nvtDBrcbTFqrNOw2e3C/SNWghg4dyvr16095vHXr1kybNo0LLriAP//5z8EPVg8NUdBr2nNf47p2xphpwDSA9PR0XftONYoidxFZnkUc8q7kiHcjXkoBoY0tmp6+lqR44mjrj6GViaaFicBefc9jDScW+PCTJ2Xk2kqJb2tjszuPTe485n/7PHwLca44hrQbwhWdr+CSlEuIdkYH5WsNF0+ueJKteVsb9D3PbXkufxj0h4DGer1evvjiC0aNqnnvcteuXfH7/WRnZwPwzTff0L9/fwDGjh3LAw880DChz1JDFPQsoGOV+ylQOcVRKkh8fh+L9i/igx0fsGj/Itx+N5HSkvbOISTZ+5LkSGP04XfP+P3t2Eg20ST7ohkZ3+bE43kX/pTlB5ez9MBSvtn/DXP2zCHSHsklKZfww+4/5ML2F2ITPZksVJWWlp4ozEOHDuWuu+6qdWzV9ZfDeZfLLOAeEZkJDAbydf+5Cpbc0lzenf1L/lu0k8O+EpJskfwopgsts2Pp7E/AhgCbKv81jLlbDn93Z8sjOIFLgIvpxy7bMY4kFzJ730L+t+d/dLDHMCa2G2NHTyEhIqHBMoSbQGfSDa3qvvLTyczMxG6307p1a7Zs2RKEZGemzoIuIm8Dw4AkEckCHgKcAMaYqcDnwFVABlAC3NFYYZU6bn/Rfl7b+BofZnxIua+c70W25f4W53NpVAecYmPuocN1v0kDsyGc42/BT1uey+9aDOSrkn28V5TBP/LXMe39yxnbYywTek+gbUzboGdTZy4nJ4ef/exn3HPPPSF/3nwgZ7mMr+N5A/yiwRIpdRqHiw/z4roX+SjjI/wGUpxDubmwjDbFMZALC8mxOuKJGbyTKG6iL8OkiLfinby5+d+8tfUtxvYYy8R+E/VMmRB2fFeMx+PB4XAwYcIE7r33Xqtj1Ul7uagmocBdwPQN03lry1v4jI8be97I4X0XEmVrRZtjM6yOd1rtTSwDou6ip38sGe6Pmbn1Hd7d+gFdXVdxjuv7OCTyxNjHf9jXwqTNT1FRUY2P+3y+Wl8zbNgwhg0b1kiJzo4WdBW65j+O3xg+Kd7FM8fWctRfxjXRqfwisR8dSmBu9kdWJ6yXaFsy/SJ/TFfX1Wwrf48d7g/Z51lI74ibaecYEvJ/zqvQpwVdhazt7mM8mreCte4jnOdK4sWWw+jtaml1rDMy4lD1vyJi2G1L57+ubawx/6S779+MKe8JvGxFPBUmtKCrkOPxe5ixYQZTD35BFA7Gu3uRXtyOA0c9HCD4BzsbS6o/gd+UXcBSx34+d+3k6agVuDe+xoTeE/QiJXVGtKCrkLL96HYeWPQAW/O2MtDXhh+U9yAWl9WxGo0N4SJvCn19ybzn2srfV/+dOXvm8JeL/0KXhMD7YCsFWtCVhY5fng8VF23s8cxhc/l/cEgU50f+hltyT70MO1zFmwjuLO+HZ/iV/HXFX7nx0xuZNGgS13e7Xvetq4BpQVeWc5tC1pVN47B3Ncn28+gf+TMibAlA8ynoAILwzdqOpDseY633BSYvmczUFV/QN/IunBKtZ8CoOmlBV5Y65stkdelzlJmj9I6YQBfnlUgzv1Q+ytaSIVF/JMM9i+3u98kv3sX5Ub8BwrugV/2LrSE0x1+AWtBV8FW2sY078i1furYRa5z8rHwgnYoPAK9am81iVc+GGQlk2gbwWsQGlhX/kS8/ncWomM4VTw6fZE3AMGO32+nbty9er5devXrx+uuvEx0deGO15557jokTJ9b6mjFjxvDUU0/RtWvXWt9j3LhxPProo3Tv3r3e+atr3lMhZQmv8fPXvFW8HbGFLv4E7isdRCd/vNWxQlJXfyL3lQ2ivT+W3+Uu5tlja/EbbVTaUI73ctm4cSMul4upU6cG/Fqfz8dzzz13YlGM6jZt2oTP5zttMQe4++67eeqpp+qVuzZa0FVQFbmLuCdnIW8XbWeYpxM/Lesf1mexNIQEE8Evygbyo9huvFKwmfuOLKLUW2p1rLAzdOhQMjIyALj++us5//zzSUtLY9q0aSfGxMbGMnnyZAYPHsxjjz3GgQMHGD58OMOHDz/l/f7zn/9w3XXXnbh/9913k56eTlpaGg899NBJ2507dy5er/esvwYt6CpoDhQdYMIXE1hedoiHWg7iOnf3U3uRqxo5sPFgiwv4feJAvirdxx1f3kFOifV9a8LF8X7offtW7Hd/5ZVXWL16NatWreL5558nNzcXgOLiYvr06cPy5cuZPHky7du3Z/78+cyfP/+U91y8eDHnn3/+ifuPPfYYq1atYv369SxcuPDEYho2m41u3bqxbt26s/469KdJBcW2vG3c8vktHC4+zAuthzMmtpvVkZocEWFC/Lk8n3QJmfmZ3Pz5zezK32V1rCbteBOu9PR0OnXqdKIf+vPPP895553HkCFD2LdvHzt27AAq9rnfcMMNAb33wYMHSU5OPnH/3XffZeDAgQwYMIBNmzaxefPmE8+1bt2aAwfOfhkJLeiq0a0+vJo7vrwDEeGN0W9wYaS2jz0bw6JTeG3Ua5T7yrnti9vYdKTher03N8f3oa9du5Z//vOfuFwuFixYwNy5c1m6dCnr1q1jwIABlJVVLGMVGRmJ3R7YVbxRUVEnXrdr1y6efvppvvrqK9avX8/VV1994jmAsrIyoqKizvrr0bNcVKO6893XWF36PFG2JPq6JjFjfikjLOhVHg6qLqzx1VEf/ewPsLzkCW7+7HbSo37D9B/dZGG6sxcqpxnm5+fTokULoqOj2bp1K8uWLat1bFxcHIWFhSQlndoKuVevXmRkZJCamkpBQQExMTEkJCRw+PBhvvjii5M6Nm7fvp20tLSzzq4zdNVovtz9JatKnyXO1pGLoh4iyqb9vxtSrK0dF0X/mWhbMitL/8bXWV9bHSksjBo1Cq/XS79+/fjTn/7EkCFDah07ceJERo8eXeNB0auvvpoFCxYAcN555zFgwADS0tK48847ueiii06MO3z4MFFRUbRr1+6ss4ux6BSo9PR0s2rVKku2rRrfJzs/4cHFD5Jo68GgqN/hkO/+nDy186Cqr6/afrf2pdsUsbzkCYrZy9OXPM2IziMsTBa4LVu20KtXL6tjNJrS0lKGDx/O4sWLT7ub5tlnnyU+Pr7G9Uxr+oxEZLUxJr2m99JdLqrhVF4w9GHRTh7KW86giDb8IK8tEQVvWRws/FT/pTiMjsxskc99C37DE62+pxcghYCoqCgefvhh9u/fT6dOnWodl5iYyIQJExpkm7rLRTWoj4syeShvOd+LbMe/ki8lAm0DGwxROJjW+jLOi0ji/twlzCnZa3UkBVx55ZWnLeYAd9xxBw5Hw8yttaCrBvNl8R4m5y1ncGRb/pF8CZE2/QMwmGJsTl5MHkZfVyt+f2QxC0qyrI6kgkwLumoQX+35ivtzlzAgIpl/JF1ChOjM3ArRNicvtB5GT1cL7j2yiCX7l1gdSQWRTqHUWZn0wQZyvBtYUfoUnXxxjMnrxZK8XKtjNWtxNhcvJQ/nruyv+PX8XzPtimkMaD3A6lgqCLSgq7NyzLeTVaXPEGtrz0+LuxCp31KWqXqeOsDN9OXJ2B3c9cXPuDD6IaaMvcqiZAGqPKjeYJrhAWHd5aLOWGZ+JstLn8Ql8QyOup8onFZHUlXE4WJw9CRsEsHy0ifYX7Tf6kghx263079/f/r06cPYsWNr7ZxYm9N1W4SK9rmZmZlARWOvmvz2t79l3rx59dpubbSgqzOSXZLNT+f8FMHGkOhJRNpaWB1J1SDalszgqPvxGzc/nfNTjpYdtTpSSAmF9rm//OUveeKJJ+qVuzZa0FW9FXuK+cVXv6CgvIDBUX8gxqa9WUJZvL0jF0T/joNFB/nVvF9R5i2r+0XNUGO3zwW47777GDhwICNGjCAnp6JbZufOncnNzeXQoUNn/TVoQVd1m//4iX/eeY9x3wc/YEfeNp5uMYgf5sxjxKEZevVniGtp78HjQx9nXc46Hlj0AH7jtzpSSAlG+9zi4mIGDhzImjVruPTSS3n44YdPPDdw4EAWL1581l+HFnQVMGMMfz26isVlB3mg5QUMjWpvdSRVD1ekXsF96ffxvz3/49nVz1odJyQEs32uzWbjxhtvBOCWW25h0aJFJ55rqPa5ekqCCtibhdt4ryiDO+N7M1b7mTcJVf9ymvsCtMdwsSuF1za9RumaZQzxdgBg5M+fsyqipY7vQ6+qavvc6Ohohg0bdtbtc2siIidua/tcFTRztxxmi/0IL0eso58vmbSDbZl7UFvgNkWCcL27O0ekhPdd20j2R3OOP0QOaIfIaYaN0T4XwO/38/777zNu3DjeeustLr744hNjt2/fztixY886u+5yUXU6JEW8EbGR9v5YbipPw4bU/SIVsuzYuLW8D61MFK9GbiBXdH3SqhqjfS5ATEwMmzZt4vzzz2fevHlMnjwZAI/HQ0ZGBunpNTZQrBdtn6tOK788n+vfupxyfPym7AJamEirI6kGkiMlPBe1kngTwUe3fEWMMyao29f2uRU+/PBD1qxZw6OPPnrKc/VtnxvQDF1ERonINhHJEJH7a3i+k4jMF5FvRWS9iIT4JWkqED6/j99//XuOShl3lvfTYh5mkk00t5X15bAU8+CiB7FqcheuqrbPPR2v18t9993XINuss6CLiB2YAowGegPjRaR3tWEPAu8aYwYA44AXGiSdstQ/v/0nSw4sYYy7J6n+BKvjqEbQw9+Sa93dmbt3LjM2Bv/U03D/JRJI+9yxY8eSmJh4yuNn8tkEMkMfBGQYYzKNMW5gJnBdtTEGiK+8nQCc/fk3ylJz9sxhxsYZjOkx5sSZECo8XertyOjU0Ty/5nkW7V9U9wsaSGRkJLm5uWFf1M+EMYbc3FwiI+v3V3EgZ7l0APZVuZ8FDK425s/A/0Tkl0AMMLKmNxKRicBEoM7fWso6mccyeXDRg/RL7sekQZP4+tvfWx1JNSJB+PP3/szO/J38/uvf884179AxrmOjbzclJYWsrKwTV0yqk0VGRpKSklKv1wRS0Gs6paH6r9TxwGvGmL+LyIXAmyLSx5iTL0czxkwDpkHFQdF6JVVBUeIp4d4F9xLpiOSZS5/BZXdZHUkFQbQzmueGPceNn97Ibxf+ljdGv0GEPaJRt+l0OunSpUujbqO5CaSgZwFVf12ncOoulbuAUQDGmKUiEgkkAdkNEVI1rkkfbAAq/sxbW/Yi+72ZDImaxHOzs4FsmsaSw+psdYzvyF8u/gu/nv9r/rbybzw45EGrI6l6CmQf+kqgu4h0EREXFQc9Z1Ubsxcqfu5FpBcQCejfUU3MXs889nsX0cN1A0mOPlbHURa4rNNl3J52O+9se4fPMj+zOo6qpzoLujHGC9wDzAa2UHE2yyYReURErq0cdh/wExFZB7wN3G70SEeTku/bzabyN0i296O763qr4ygL/WrgrxjYeiAPL32YzPxMq+Ooegjo0n9jzOfA59Uem1zl9mbgooaNpoLFa0pZU/o8LoljQOTPEdELiJubuS/830n3vy9RbIvycvd/b+b/ytIZ3auyEVuIXJ6vaqa9XJqzyiW/jhydR4njED8vG0i3ovcsDqVCQaKJZHx5b6ZHrmOWK4PRaGfNpkCnYs3cJ8W7WOU8xBWeLnQLlSZNKiSk+ZK41NORRc4svirZV/cLlOW0oDdjezwFPJq3kq6+RC73pFodR4Wga9zdSPHFMTlvOYe8xVbHUXXQgt5Mefwe7s9dglNs3FKehl2/FVQNHJWdGT3Gzx9zl+Lz+6yOpE5Df4qbqRfXvshGdx4PtRykTbfUaSWbaCa1OJ+V5dm8vvl1q+Oo09CC3gytObyGGRtncF1MV66I1hYMqm7Xx3RlZFRH/vntP9mSu8XqOKoWWtCbmUJ3IX9c9Efax7RnUovz636BUlQsl/ZQy0G0iGjB/d/cT6lXF8UIRXraYjPz5IonOVh8kNdHvU7Mpi+sjqOaiLlbKpYc7Jr0Y5bnP86P3v0TfSJvA+DxH/a1MpqqQmfozcj8vfP5eOfH3NXnLvq37m91HNUEJTv6kuq8kt2e2RzxbrI6jqpGZ+jNwKQPNuD2F7Cg5EHibZ3Zv2sok3ZvYMQhXehZ1V+viHHkeNezruwlLol5wuo4qgqdoTcDxhg2lL+CxxTTP/JubKK/x9WZs0sE/aN+RqnJZXPZm1bHUVVoQW8GDniXctC7gp6uscTb9awWdfZa2LvTzXUt+7wLWbBvgdVxVCWdqoW5I6VH2Fj2Gom2bnR1XW11HNXEjTj03bqjl+Ln2chYHln4ewa0u5oEW7XFULSRV9BpQQ83lQ23oGJXy2NHFmFMCT8rSqJN4WvW5VJhx4GNce5e/MO+iqeOruaxVhdaHanZ010uYWx2yV7mlu5jlKcLbUyM1XFUGOroj+fO+N7MKt7FN6W6NrzVtKCHqTxfGX89uoo+rpYM8+h+c9V4fpbQh3OcCTyct4JCv9vqOM2aFvQw9cTR1RT6PTzScog23lKNyiV2Hm05mBxfKc8cW2t1nGZNf9LD0MLS/XxRsoeJCWl0dyVaHUc1A30jkpgQ15P3izJYWabXN1hFC3qYKfZ7eDRvJd2cCfw4vrfVcVQz8ouEfqQ4Ynk4bwXlRtvsWkHPcgkTkz7YAMCxvKVkO0oYV5bOwmNHLE6lmoPjfV4ArrF1Z2rUt9y/czmRxyq+J7XXS/DoDD2M5Pm2s9iRxcXeFFL9CVbHUc1QT39LBnnaMc+5hwLfHqvjNDta0MOE33jZUDadRBPJVe5zrI6jmrFr3d2JxsG6spcxxm91nGZFC3qYyHR/RqE/ixvcPYnUPWnKQjE4ub68B/n+THZ75lgdp1nRgh4G9hXsY7v7A9o6LiDNl2R1HKUY6GtDkr0v28rfJbsk2+o4zYYW9CbOGMNjyx/DhoO0iNusjqMUAILQN/IO/Hh5YoW22A0WLehN3Je7v2TxgcX0jBhLlK2l1XGUOiHG1pburh8wZ88cvs762uo4zYIW9Cas0F3IUyufoner3qQ6r7A6jlKnOMd1DV0TuvLX5X/VdUiDQAt6EzZl7RRyS3OZPGQyIvpfqUKPTRw8OORB9hftZ/qG6VbHCXtaBZqorXlbeXvr24ztMZa0pDSr4yhVqwvaXsDVXa/m1Y2vsqdAz01vTGKMsWTD6enpZtWqVZZsu0mq0ufcbwy3Hp7DPm8hs9p/nwSb66Sr9ZQKJSN7teGIr5TvH/iUfhGtmJo8HBH5boAuhFEvIrLaGJNe03M6Q2+CPirOZJ37CL9JHHDqKjFKhaAkexT3JPZjSdkh/le6z+o4YUuvQGkijs/Ai/HwVPQauvgTiN4XzVx0Zq5C2/Hv3STi6RAZx6PZK/GW2omoLD8jh1uZLrwENEMXkVEisk1EMkTk/lrG/EhENovIJhF5q2FjquO+dGVSgocb3D2xIXW/QKkQYcfGDe4e5NvKmevcbXWcsFTnDF1E7MAU4HIgC1gpIrOMMZurjOkOTAIuMsYcFZHWjRW4OdtvK2SxI4uLvCl08MdZHUepeuviTyTd05b5zr0M8rYn2URbHSmsBDJDHwRkGGMyjTFuYCZwXbUxPwGmGGOOAhhj9FrfBmYwfODaTjRORrm7Wh1HqTP2fU83HNj40LXd6ihhJ5CC3gGoehQjq/KxqnoAPURksYgsE5FRNb2RiEwUkVUisionJ+fMEjdTa+yHybQf42r3OcTgtDqOUmcs3kQwyt2VLY5cNtm1Z39DCqSg17Sjtvq5jg6gOzAMGA9MF5FT1j4zxkwzxqQbY9KTk5Prm7XZKvGU8IkrgxRfHIO97a2Oo9RZG+pNoY0/mg9d23H7dGHphhJIQc8COla5nwIcqGHMx8YYjzFmF7CNigKvGsD0DdPJt5XrgVAVNuzY+EF5T3Jtpby5+U2r44SNQAr6SqC7iHQRERcwDphVbcxHwHAAEUmiYhdMZkMGba72Fe7j9U2vc763ra5CpMJKT39L+niTmLZ+Gjklugu2IdRZ0I0xXuAeYDawBXjXGLNJRB4RkWsrh80GckVkMzAf+J0xJrexQjcnz6x6BrvNzjW6CpEKQ9e6u+P2u/nHmn9YHSUsBHQeujHmc2NMD2PMOcaYxyofm2yMmVV52xhj7jXG9DbG9DXGzGzM0M3F8oPLmbt3Lnf1uYtEE2l1HKUaXLKJZkLvCXy882M2HtlodZwmTy/9D1E+v48nVz5J+5j23JamC1eo8DWx70RaRbbiiRVPYFVvqXChBT1EfZDxATuO7uDe9HuJdOjsXIWvWFcsvx74a9blrOPL3V9aHadJ04IegorcRfzr238xsPVAruisC1eo8HftOddybstzeXb1s5R5y6yO02Rpc64QMumDDQBsKX+bPHcevW2/4Y8fVuxXHGFlMKUa0fHv+0TvGLaW/oVb3vsb3SOuB+DxH/a1MlqTozP0EFPiz2aX+ws6OC4m0a5ntqjmI8nRmzaOdDLcsyjzH7M6TpOkM3SrVVm4YsShw7wesQG73c+PCwyJ+TMsDKZU8PWKGM9C7+/Z5n6P8yJ/YnWcJkcLegjZbctnrSObK91d9DRF1WyMOHTyxMW42vO1mc+4/GKY3+27J3RlozrpLpcQYYzhY9cO4v0uhns6Wx1HKctc4e5CFA4+cWVYHaXJ0YIeIuaU7mO3PZ/RnnOIwG51HKUsE42Ty91d2OrIY3Fp9bZR6nS0oIcAj/Hx7LG1tPPHMMjbzuo4SlnuYm8KrfxR/P3YWnzGb3WcJkMLegiYWbiDLG8R33d3126KSgEObFztPocdnmPMKt5ldZwmQwu6xfL9bqYWbOTCyLac62tpdRylQkZ/X2v6uVrxz/z1lPi9VsdpErSgW2xG/iYK/W7uTRyA6OxcqRME4bctBpLjK+XNwq1Wx2kStKBb6GDRQf5TuI1rolM519XC6jhKhZwBEckMj0rh1YLN5JXlWR0n5GlBt9C/1v4LA9yTeJ7VUZQKWf+XeB6lxsdL616yOkrI04JukW152/hk5yfcHNeT9o4Yq+MoFbK6OhP4Yew5vLvtXfYW7LU6TkjTK0WD6HgTIoDlJU/iIIpuh5KZe+iwhamUCn0/T+jLZ2X7ef7b53n60qetjhOytKBb4Ih3Ezm+dfSKGE9McbbVcZQKWXO3fDfZSWlxJbN3f0Rx9neN67Qb48l0l0uQGWPYWj6TSGlJqvNKq+Mo1WR0dV2DS+LYWq4rXNZGC3qQHfKu5Jh/Jz0ixmAXl9VxlGoynBJNN9f1HPFtIse7oe4XNENa0IPIb3xsdb9DrK0DKY6hVsdRqsnp7BxJlCSxtXwmRlsCnEILehBleb6m2H+Qnq4fYRNtwKVUfdnFSY+IMeT7d3HQu8LqOCFHC3qQlHnL2O7+L4m2brR1pFsdR6kmK8VxMXG2jmwrfxeP32N1nJCiBT1IZm6dSZnJo1fEeET0En+lzpSIjXMjfkSxOcRHGR9ZHSekaEEPgiJ3EdM3TifZfh6tHL2sjqNUk9faPpAWth5MXTeVMm+Z1XFChp6H3liqrBX6+rH15Jfnc2+pi44Fuk6oUmei+lJ1qbZYpvi3885nE7ktvspEqRkvVacz9EaW5yvjjcKtXB7VkY7+eKvjKBU2uvlb8L3Itkwv2EyR7ksHtKA3uukFmygzPu5J7Gd1FKXCzq8S+3PMX84bBVusjhIStKA3okPeYt4p3MG1MV3o6kywOo5SYSfN1ZLLozryeuFW8ny6L10LeiOamr8RP3B3gvabUKqx3JPYjzLj45WCzVZHsZwW9Eay11PIR8WZjI3tpu1xlWpEXZ0JXBOTysyiHWR7S6yOYyk9y6UBVW2Pe/DoSmx24dzs1sw9rO1xlWoMx7sx9pV2fBq1m8l7VtKi8uewOXZiDGiGLiKjRGSbiGSIyP2nGTdGRIyINOtLIQt9WayxH2KopyPxJsLqOEqFvSQTxRBve5Y5DlDiz7E6jmXqLOgiYgemAKOB3sB4Eeldw7g44FfA8oYO2dRsd7+PCzuXeTpbHUWpZuNyTyqCsMP9odVRLBPIDH0QkGGMyTTGuGET2iEAABM+SURBVIGZwHU1jHsUeApo1oea8327OehdwTBPJ2JwWh1HqWYj0URykbcDWZ6vKfIftDqOJQIp6B2AfVXuZ1U+doKIDAA6GmM+Pd0bichEEVklIqtycsLzz6Jt5e/hJIZLPZ2sjqJUszPCnYoNJ9vL/2t1FEsEUtBr6iRlTjwpYgOeBe6r642MMdOMMenGmPTk5OTAUzYRR30ZZPu+5RzXNUTp8Walgi4OF6muKzngXUrG0Qyr4wRdIAU9C+hY5X4KcKDK/TigD7BARHYDQ4BZzfHA6Pby93BJHKkuXVpOKauc47oaB5G8sO4Fq6MEXSAFfSXQXUS6iIgLGAfMOv6kMSbfGJNkjEk1xqQCy4BrjTGrGiVxiFpzeA05vg2c4/o+Dom0Oo5SzZZL4ujiGsWcPXPYmrfV6jhBVWdBN8Z4gXuA2cAW4F1jzCYReURErm3sgE3FlLVTiJAEUp2XWx1FqWavi2s0cc44XljbvGbpAe3oNcZ8Dnxe7bHJtYwddvaxmpYVB1ew4tAKekdMwC563rlSVnNJLLem3cqUtVPYdGQTaUlpVkcKCr30/ywZY5iydgqto1rT2TnC6jhKqUq39LqFhIgEpqydYnWUoNGCfpaWHVzGmuw13NX3LuzisjqOUqpSrCuW23rfxjf7v2F9znqr4wSFnlt3FowxvLD2BdpEt+GGHjewYct2qyMp1eydWNlofhtu8nt4wxbBC/N+y9TWw08eGIYrG2lBr68qS8stLT3I2py1PNjiAiK+foYRh7QJl1KhJMbm5Pb4Xjx3bC1ry3PoHxF+179UpbtczpAxhin562lrj+YHsV2tjqOUqsX42O60sEXwYv6Gugc3cTpDP0NLyg6y3p3Ln1pcgEvsVsdRSlVzvLUuwMXOjnziz+Dlbdvo4k8EYOTw2l7ZdOkM/QwYY3ghfwPtdHauVJNwkSeFWOPkS9cuq6M0Kp2h19PcLYfZas9lfWQuY8vPZeHWI1ZHUkrVIQI7l7k7Mysig122Yydm6eFGZ+j1ZDB86cykhT+SQd52VsdRSgXoe96KWfpsZ/jO0rWg19M2ex577AWM9HTGoR+fUk1GBHaGezqzzZHHblu+1XEahVakejDGMNu5i0R/BIO87a2Oo5Sqp4s8KcQYJ7OdmVZHaRRa0Oth6YGl7LbnM9KTqrNzpZqg47P0rY481uWsszpOg9OqFCBjDC+ue5FEfwSDdXauVJN1sacDMcbJi+tetDpKg9OCHqDlh5azNmctI3R2rlSTFoGD4Z5OLN6/mA054XWxkVamAE1dN5XWUa0ZrGe2KNXkXeRJISEigZfWv2R1lAalBT0AKw+tZPXh1dzZ906c6FWhSjV1kTi4tfetLMxayObczVbHaTBa0APw0rqXSIpK4obuN1gdRSnVQMafO544VxwvrQufWboW9DqsObyG5YeWc0faHUQ6dK1QpcJFnCuOCb0mMG/fPLblbbM6ToPQgl6Hl9a/RMvIloztOdbqKEqpBnZTr5uIdcaGzb50LeinsT5nPUsOLOG2tNuIckRZHUcp1cASIhIYf+545uyZQ8bRDKvjnDVtzlXNpA++O41pRcnTOCWWLdv6Mml7xeO6aqhS4eXW3rfy7y3/ZtqGaTx1yVNWxzkrOkOvRb5vF9m+b+nqvAqH6L5zpcJVYmQi484dx+zds9mdv9vqOGdFZ+i12OH+CCfRpLqusDqKUqoRzH3h/07c7oIbezQ88t5EbnL3BmBkrzYVTzahtUe1oB9XuVboiEOHOShFfBq9kivcXRh1+G2LgymlGlscLi70dmCRI4srPV1oZZrmMTPd5VKDOa7dRBg7l3g6Wh1FKRUkl3k6Y0P4yrnb6ihnTAt6NdlSzFr74Yo2mzitjqOUCpIEE8FgbztWOA5yTMqsjnNGtKBX85VzDw5sDPN0sjqKUirILvN0xgDznHusjnJGtKBXkeUtYpXjEBd6OxCHy+o4Sqkga2miuMDblmWOAxzxlVodp960oFfxSsFmBBius3Olmq0RnlS8+HmjYKvVUepNC3qlQ94SPirKZLC3PYlGzztXqrlKNtEM8LVhZtEOjvnKrY5TL1rQK71euAU/hhGezlZHUUpZbKQ7lVLj5d+FTatplxZ0ILc0l/eLMrg6JpWWTfT8U6VUw2lnYhkZ1ZG3CrdR6C60Ok7AAiroIjJKRLaJSIaI3F/D8/eKyGYRWS8iX4lIk5rmvrn5TcqNjx/Hp1kdRSkVIn6ckEah8fDOtnesjhKwOgu6iNiBKcBooDcwXkR6Vxv2LZBujOkHvA80mQ43+eX5zNw2kyuiO9HFGW91HKVUiEhzteSiyHa8uflNSr1N44yXQGbog4AMY0ymMcYNzASuqzrAGDPfGFNSeXcZkNKwMRvP21vfpthTzE90dq6UqmZifBp5ZXn8d/t/rY4SkEAKegdgX5X7WZWP1eYu4IuzCRUsJZ4S/r3l3wxLGUZPVwur4yilQszAyNakt0nn1U2v4va5rY5Tp0Cac0kNj5kaB4rcAqQDl9by/ERgIkCnTtad63285/lO92fkl+fjybuMuTkLLcujlAo9c7ccBsCedDnZpY9z+3sv0tl1GQCP/7CvldFqFcgMPQuo2qUqBThQfZCIjAQeAK41xtR48qYxZpoxJt0Yk56cnHwmeRuMz7jJdH9Gkj2NFvZulmZRSoWuJHsfEmxd2emehd/4rI5zWoEU9JVAdxHpIiIuYBwwq+oAERkAvERFMc9u+JgNb59nIeXmGN1c11sdRSkVwkSE7q7rKTHZHPAutTrOadVZ0I0xXuAeYDawBXjXGLNJRB4RkWsrh/0NiAXeE5G1IjKrlrcLCX7jZaf7UxJt3Whlr37CjlJKnayNYyBxthR2umdhjN/qOLUKaIELY8znwOfVHptc5fbIBs7VqA54l1JqcugTeSsiNR0iUEqp74jY6Oa6jm/LpnDYuwY4z+pINWp2V4r6jZ8M9yzibB1pbR9gdRylVBPRzjGEaGnNDvdHGFPjeSGWa3YFfd7eeRT599PNdR0ize7LV0qdIZvYOcf1ffL9mSw7uMzqODVqVhXNGMPLG14mWtrQ3jHE6jhKqSYmxXkJEdKClze8bHWUGoX/ItGViz8DLC09yObczfyo/FwuLHrVwlBKqaZixKEZJ913OZL4+NBK1n55L/0jqpx+PXxSkJOdqlnN0F8u2ERrexQXeNtZHUUp1URd6G1Poi2CGQWbrY5yimZT0NeW57CqPJvb43rhaD5ftlKqgUXg4Oa4Hiwo3c929zGr45yk2VS26fmbSLRFcEOsXhWqlDo74+N6Ei0OZhRssjrKSZpFQd/mPsrCsgPcHNeDaFv4HzZQSjWuBJuLG2O782XJXvZ5QmcBjGZR0GcUbCZaHIyP62l1FKVUmJgQfy4OhFcKt1gd5YSwL+h7PYXMLtnLjbHdSbC5rI6jlAoTyfYofhB7Dh8XZZLtLan7BUEQ9gX9lYLNOBAmxJ9rdRSlVJi5Pb4XfgxvFG61OgoQpuehH+93XuY/yvyiTAZ727NuewFQYG0wpVRYON4rHaB/RBveLtjBwf8uxSWxgHX90sN6hp7p/hwDDPc0qTWrlVJNyAh3Z9ziY7f7f1ZHCd+C7jZF7PHMZYCvDUkmyuo4Sqkw1c7EkuZNYpf7S7ymzNIsYVvQd7v/h49yRrh1dq6UalwjPal4KGKvZ56lOcKyoHtNGbvcX9LGPpB2JtbqOEqpMJfqT6ClvRc73Z/hMx7LcoRlQd/rmYeHIrpFXFv3YKWUagDdXNdSbo6y37PIsgxhV9DdPjeZ7s9pae9FC3sPq+MopZqJZHs/4m2p7HR/gs9vzWLSYVfQP838lDKTRzeXzs6VUsEjInRzXUuxOcTcvXMtyRBWBd3n9/HKxleIt6WSbO9ndRylVDPTzjGIGGnHjA0zLFmmLqwK+ty9c9lTsIdurmt18WelVNCJ2DjHdQ1b8raw5MCSoG8/bAq6MYYZG2aQGp9KO8cgq+MopZqpFOdQWke3ZvqG6UHfdtgU9CUHlrAlbwt39LlDF39WSlnGJg5u630bqw6vYm322uBuO6hba0QzNs6gdVRrrul6jdVRlFLN3JgeY0iISGDGxhl1D25AYdGca13OOlYeWslv03+Ly64tcpVSwVd1MenoRW24KaIjL+5bwI7Zf6C7K/G7gY24mHTTLOjzHz/p7vSchcTbXIw9tAeyH2fEocO1vFAppYLjprievFa4hVcKNvN40veCss0mv8slw32MBaX7uSm2B9E2p9VxlFIKgER7xRrGX5TsYb+3KCjbbJoz9CpeKdhMlNi5WZeXU0qFiOP90s+RJIjazl/2rOEGd0WNGjm88bbbJAv68Q8rT0r5LGoPF3tTWLX9mMWplFLqZIkmknRvW5Y7DnCFO5U4Ihp1e016l8t8514Ahnk6WZxEKaVqdpmnM178fO3MavRtNdmCXoib5Y4DpHvb0sJEWh1HKaVq1NrE0NeXzCJnFmV4G3VbTbagf+Pchxc/l+nyckqpEDfCk0qZeFns3N+o22mSBb0ML4ucWfT1JdPGxFgdRymlTquTP54evpYsdOyl3FfeaNsJqKCLyCgR2SYiGSJyfw3PR4jIO5XPLxeR1IYOWtUS535KxcsIT6NuRimlGsxId2cKbW4+zvi40bZRZ0EXETswBRgN9AbGi0jvasPuAo4aY7oBzwJPNnTQ48p95Sxw7KW7rwWd/PGNtRmllGpQ3fwt6OSL59WNr+L1N86+9EBm6IOADGNMpjHGDcwErqs25jrg9crb7wMjpJH6136y8xMKbW5GulMb4+2VUqpRCMIIT2eyirKYs2dOo2wjkPPQOwD7qtzPAgbXNsYY4xWRfKAVcKTqIBGZCEysvFskItvOJDSQ9As2Hql7mCWSqPZ1hxDNVn+hmgs025myPNtVXFXTw4HmqvVMkEAKek0z7epLcQQyBmPMNGBaANs8fSCRVcaY9LN9n8ag2c5MqGYL1Vyg2c5UqGZriFyB7HLJAjpWuZ8CHKhtjIg4gAQg72yCKaWUqp9ACvpKoLuIdBERFzAOmFVtzCzgtsrbY4B5xooF9ZRSqhmrc5dL5T7xe4DZgB14xRizSUQeAVYZY2YBM4A3RSSDipn5uMYMTQPstmlEmu3MhGq2UM0Fmu1MhWq2s98drRNppZQKD03ySlGllFKn0oKulFJhIqQLeqi1HKhntntFZLOIrBeRr0QkaF3E6spWZdwYETEiEpRTuALJJSI/qvzcNonIW8HIFUg2EekkIvNF5NvK/9MaTyRuhFyviEi2iGys5XkRkecrc68XkYHByBVgtpsrM60XkSUicl6oZKsy7gIR8YnImFDJJSLDRGRt5c/AwnptwBgTkv+oOAC7E+gKuIB1QO9qY34OTK28PQ54J4SyDQeiK2/fHUrZKsfFAV8Dy4D0UMgFdAe+BVpU3m8dKp8ZFQes7q683RvYHaRslwADgY21PH8V8AUV14IMAZYHI1eA2b5X5f9ydChlq/L/Pg/4HBgTCrmARGAz0Knyfr1+BkJ5hh5SLQfqm80YM98YU1J5dxkV5+8HQyCfG8CjwFNAWQjl+gkwxRhzFMAYkx1C2QxwvHlQAqdei9EojDFfc/prOq4D3jAVlgGJItIuFLIZY5Yc/78kuD8DgXxuAL8E/gsE6/sskFw3AR8YY/ZWjq9XtlAu6DW1HOhQ2xhjjBc43nIgFLJVdRcVs6hgqDObiAwAOhpjPg1SpoByAT2AHiKyWESWicioEMr2Z+AWEcmiYkb3y+BEq1N9vxetEsyfgTqJSAfgB8BUq7NU0wNoISILRGS1iNxanxeH8pqiDdZyoBEEvF0RuQVIBy5t1ERVNlnDYyeyiYiNio6Ytwcpz4lN1/BY9c/MQcVul2FUzOa+EZE+xpjGXjA2kGzjgdeMMX8XkQupuO6ijzHG38jZ6mLVz0DARGQ4FQX9YquzVPEc8AdjjC84f9QHzAGcD4wAooClIrLMGLM90BeHqvq0HMgKcsuBQLIhIiOBB4BLjTGN19W+ftnigD7Agspv5LbALBG51hizysJcx8csM8Z4gF2Vzdu6U3G1cmMKJNtdwCgAY8xSEYmkoplS0P5cr0VA34tWEZF+wHRgtDEm1+o8VaQDMyt/BpKAq0TEa4z5yNpYZAFHjDHFQLGIfA2cBwRU0INygOIMDx44gEygC98dqEqrNuYXnHxQ9N0QyjaAigNt3UPtc6s2fgHBOSgayGc2Cni98nYSFbsSWoVIti+A2ytv96KiaEqQ/k9Tqf0g2tWcfFB0RZC/306XrROQAXwvmJkCyVZt3GsE6aBoAJ9ZL+Cryu/JaGAj0CfQ9w7ZGboJzZYD9cn2NyAWeK9yFrDXGHNtiGQLugBzzQauEJHNgA/4nQnCrC7AbPcBL4vIb6jYpXG7qfwJbEwi8jYVu6CSKvffPwQ4K3NPpWJ//lVUFM4S4I7GzlSPbJOpOKb1QuXPgNcEqcthANksUVcuY8wWEfkSWA/4genGmNOeennS+wfhe1IppVQQhPJZLkoppepBC7pSSoUJLehKKRUmtKArpVSY0IKulFJhQgu6UkqFCS3oSikVJv4fCDQF0I5NAI0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "num = 10**6\n",
    "\n",
    "# Plot histogram from part (a)\n",
    "s_arr = sample_s(num)\n",
    "th_arr = np.arcsin(s_arr/R)\n",
    "plt.hist(th_arr,bins=40,density=True,alpha=0.6,label='Part (a)')\n",
    "\n",
    "# Plot histogram from part (b)\n",
    "th_arr = sample_th(num)\n",
    "plt.hist(th_arr,bins=30,density=True,alpha=0.5,label='Part (b)')\n",
    "\n",
    "# Plot curve\n",
    "x = np.linspace(0,np.pi/2,num=100)\n",
    "y = np.sin(2*x)\n",
    "plt.plot(x,y,label='PDF')\n",
    "\n",
    "plt.legend()\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}