mp573_hw4_q7_MonteCarlo.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Homework 4 \n",
    "By: Lawrence Lechuga\n",
    "Date: 11/04/2020\n",
    "\n",
    "## Problem 7\n",
    "### Method 1:\n",
    "\"Throws darts\" at the upper quadrants and takes the ratio of successes falling within the circle relative to the total number of trials. This is a purely $\\text{binomial}(n,p)$ process where the probability of success is the ratio of the areas and $n$ is the number of trials. \n",
    "### Method 2:\n",
    "Takes the functional form of the circle and calculates the average dimension in y. The area is then calculated by taking the mean value of $f(x)$ times the diameter to estimate the area. This is a Gaussian process that will converge to the correct answer based on the CLT."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "import random\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "<ipython-input-24-38525ae8e93f>:11: RuntimeWarning: invalid value encountered in sqrt\n",
      "  G = np.longdouble(np.sqrt(radius**2-X**2));\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD6CAYAAACxrrxPAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3dfXSb5Zkm8OuW/BHixITEApuExGlrphZlSMFro/Q0zZ5Qm9DpBPrBgXEmlMMgZ1u2UAo2ZXdKuzNMsdltw+wAkUozQ8YGDl1ayMwGbJY9gT2xSAjTFJpkg1NqE0OycWLy6WDF9r1/SHEt6ZUs25Je6dH1O8cnfj8k3ULkyqvnfT5EVUFERLnPYXcBRESUGgx0IiJDMNCJiAzBQCciMgQDnYjIEAx0IiJDMNDJaCLyoIg8Nc3H9orIdamuiShdCuwugCgVROQvANwL4LMATgHYDeBhVf07WwsjyiBeoVPOE5F7AWwA8HcALgGwGMATANZM8jhe0JBRGOiU00TkQgD/BcB3VPVXqnpGVc+p6r+o6v0i8iMRaQ+fWykiKiJ3iMgHAP53eP+dIrJPRE6JyF4RudridRwi8oCI/F5EjonI8yIyP6NvlmgSDHTKdR4AswD8egqP+RKAagANIvJNAD8CsA5AKYA/B3DM4jHfBXBj+LGXAvgYwOPTrpooDfiVk3LdAgBHVXVkCo/5kaqeAQAR+SsAbar6VvjYgTiPaQJwl6r2hx/3IwAfiMhfTvG1idKGgU657hiAMhEpmEKwHpzw+2UAfp/EY5YA+LWIjE3YN4pQm/2HSb4uUVqxyYVyXQDAJwg1hyRr4hSjBwF8OonHHASwWlXnTfiZpaoMc8oaDHTKaap6AsAPATwuIjeKyGwRKRSR1SLSlsRTPAXgPhG5RkI+IyJLLM7bCODh88dExCUiCXvREGUaA51ynqr+FKE+6P8ZwABCV9N3AXgxicf+EsDDAJ5BqP/6iwCseq88BmALgC4ROQXgTQB1qaifKFWEC1wQEZmBV+hERIZgoBMRGYKBTkRkCAY6EZEhGOhERIZgoBMRGYKBTkRkCAY6EZEhGOhERIZgoBMRGYKBTkRkCAY6EZEhGOhERIZgoBMRGYKBTkRkCAY6EZEhGOhERIZgoBMRGYKBTkRkCAY6EZEhGOhERIZgoBMRGYKBTkRkCAY6EZEhGOhERIZgoBMRGYKBTkRkCAY6EZEhGOhERIZgoBMRGYKBTkRkCAY6EZEhGOhERIZgoBMRGYKBTkRkCAY6EZEhGOhERIZgoBMRGYKBTkRkCAY6EZEhGOhERIZgoBMRGYKBTkRkCAY6EZEhGOhERIZgoBMRGYKBTkRkCAY6EZEhGOhERIZgoBMRGYKBTkRkCAY6EZEhGOhERIZgoBMRGYKBTkRkCAY6EZEhGOhERIZgoBMRGYKBTkRkiAK7XrisrEwrKyvtenkiopz09ttvH1VVl9Ux2wK9srISu3btsuvliYhykoj0xTvGJhciIkNMGugisklEjojI7+IcFxH5exE5ICLviMjVqS+TiIgmk8wV+j8BuD7B8dUAqsI/XgBPzrwsoswKBAK46aabUFFRgblz56KkpASFhYVwOBwQEbjdbrtLJJrUpG3oqvqGiFQmOGUNgM2qqgDeFJF5IlKhqodSVCPRlLW0tGDDhg0IBoNwOELXLWNjYxARAEBJSQkKnE6cOn0ao6Ojkz7fvn37xh8LAAUFBXA6nSgqKsLIyAjGxsYwOjqK0dFROBwOrFq1Cp2dnel5c0RxpOKm6EIAByds94f3xQS6iHgRuorH4sWLU/DSlM/8fj8eeughHDlyBKo6HtwCYGRCSI+NjY3/HrruAE6fPj2j1x4ZGcHIyAiGh4djjo2OjqKrqyviH4CJqqursXfv3hm9PpGVVNwUtfq/Vq1OVFW/qtaoao3LZdnrhshSQ0PDePNHcXExnE4nmpqacPjwYYyNjUFVx6+QR5K44rbT+at9h8MBp9MJh8OBuro6u8siA6Qi0PsBXDZhexGAj1LwvJTH/H4/FixYMB7iXV1d41fXwWAw4qo71YqKijB79mwUFKS3V6+qjv9jtHPnTogIRAQlJSVoaWlJ62uTmVIR6FsArAv3drkWwAm2n9NUBQIBfP7znx8P8KamJgwODo6H+Ew4HI4/NseEQ3POnDmYN28enE4nAKCwsBButxs+nw/Dw8M4c+YMzp07B1VFc3NzTLgXFBSguLgYc+fOxQUXXIDi4uKU/QMwNDSEtra2iID3+/0peW4ynKom/AHwLELt4ecQuhq/A8B6AOvDxwXA4wB+D+BdADWTPaeq4pprrlGi5uZmdblcilAz3Yx+nE6nOp1OLSkp0erqavX5fLa8J5/Pp7Nnz07Je4r+mT17tjY3N9vyvig7ANilcXJVNAVXQNNRU1OjHCmanwKBAG644QYcP3582s9RWFiIc+fOweFw4LrrrsuZHiULFizA4OAggNC3hZn8/SspKcFPf/pTeL3eVJVHOUBE3lbVGqtjHClKGeH3+1FSUgIRwfLly6cd5tXV1VBVBIPB8RuhuRLmAHDs2LHxq6nz7eeNjY3Taq45c+YMmpqaUFJSgosvvpjt7sRAp/RpaWlBcXHxeJv40NDQlB5fVVWF7u7uiK+UJnb3a29vH2+vV1XU19fH7fJoZWhoCAMDA2hra0NBQQHWrl2bxmopmzHQKeXq6uogImhra0MwGEzqMUVFRSgtLUVjY+N4sL333nvweDxprjb7dHZ2jl+9+3w+lJeXj9/Unczo6Cg6OjrgdDpRUVHBm6l5hoFOKeP3+1FQUICdO3cm/RgRQWNjI4aHh3HixAm0t7enscLc4/V6cejQIYyOjo43zyQT7mNjYzh8+DCampqwYMGCDFRK2YCBTjMSCARw+eWXjw/0SWYYvcvlGm9KGRsbY4hPQXt7+3i4+3w+zJ8/f9LHDA4Oori4GIFAIAMVkp0Y6DRtLS0tWL58OXp6epIa6FNbWwtVxZEjR/KyKSXVvF4vjh07hubm5og+9VaCwSCWL18+3neeN1DNxECnKfH7/SgrKxtvI0/E6XTC5/ONt4nv2LEjQ1Xml9bWVnz88ccYGRlBd3c3PvOZz8Q9NxgM4vTp02hra8OsWbN4A9UwDHRKWkNDA5qamnDs2LGE5zmdTtTX12NkZIR9pDPM4/Ggp6cH3d3dKCoqSnju8PAwOjo62DPGIAx0mtTatWvhcDjQ1dWV8LzS0lL4fD6MjIzkVN9wE3k8HgwPD6O5uRmLFi1CWVlZ3HPP94y56KKL2M6e4xjoFFcgEIDL5UJHR0fCEY3n+4ufOHGCV+RZprW1FQcPHsTAwAAaGxsxa9asuH3cjx8/juXLl/NqPYcx0CmG3+/H0qVLsXz5chw9ejTueQUFBWhubs7b/uK5pr29HWfPnsX27dsTDlzq6OhAcXExb5zmIAY6RXC73WhqakJvb2/cc2bNmgWfz4dz586htbU1c8VRSng8Hmzfvh0rVqyIe04wGBwfecpgzx0MdAIQaicvKCjAvn37Ep5XX1+Ps2fPsmklx3k8Hrz++uvo7u7GvHnz4p43OjqKtrY2NDQ0ZLA6mi4GOsHtdqOjoyPhoKD6+nqoKm92Gsbj8eDjjz9GY2Njwn7sXV1dcLlcvGma5RjoeaylpQUikvCq/HzPFQa52drb2zEyMoLm5ua4wX706FHeNM1ySQW6iFwvIvtF5ICIPGBx/EIR+RcR+a2I7BGR21NfKqVSS0vLpAOD6uvr2XMlz7S2tmJkZARLliyJe05HRwdmz57Nq/UsNGmgi4gToRWJVgNwA7hVRNxRp30HwF5VvQrASgD/TUQSj2og2zQ0NODRBGG+YsUKdHd386o8j/X29qK+vj7u8bNnz2L58uVsW88yyVyh1wI4oKrvq2oQwHMA1kSdowDmSqgv1BwAgwBGUlopzVggEMDsCy4ILbhscdzhcMDn8+H1119nN0RCZ2cnfD4fCgsL457T1dXFUM8iyQT6QgAHJ2z3h/dN9A8AqgF8hNC6oneravqWZacpOz+R1tlPPrE8Xl5ejtHRUTavUASv14tgMJjwar2rqwtud/SXdrJDMoFuNQIh+gKvAcBuAJcCWAbgH0SkNOaJRLwisktEdg0MDEy5WJqehoaGuO3l53uvHDp0KMNVUS7p7OxEd3c35syZY3l83759qKioyHBVFC2ZQO8HcNmE7UUIXYlPdDuAX4UXpT4A4A8APhv9RKrqV9UaVa1xuVzTrZmmwO12x52DZf78+Wwnp6R5PB6cOnUKjY2NlscPHz6M0tJSrpJko2QC/S0AVSKyNHyj8xYAW6LO+QDAKgAQkUsA/AmA91NZKE2N3+9HcVFR3C6J5eXlk86aSGSlvb0dPp8PJSUlMcdOnTqFpqYmdm20iSSadGn8JJEbAGwA4ASwSVUfFpH1AKCqG0XkUgD/BKACoSaaR1Q14TI0NTU1umvXrhmWT1bcbnfCvuXV1dVGLrZMmVdRUYHDhw9bHluyZEnCKSRoekTkbVWtsTyWTKCnAwM9PRKFeUFBAe69917Ov0IpVVdXF3cdWYZ66iUKdI4UNUhdXV3CJhZOpkXpsGPHDtTW1loe6+vrY7fGDGKgG8Ltdse9SqqurmYvFkqrHTt2xJ02oKurC5WVlZkvKg8x0A1QWVkZ98q8traW7eWUEYmmDejr60NJSQl7wKQZAz3Hud1u9PX1WR6rra3lwsyUcfGmDRgaGkJTUxPnV08jBnoOW7t2bdwr8+bmZoY52aazszPuBF9tbW3s1pgmBXYXQNNTWVkZ98rc5/NxCD/Zrre3N26vq46ODixcuJA36VOMV+g5qLS01DLMXS4Xuru7GeaUNfbu3Qufz2d5s/Sxxx6zoSKzMdBzjNvtxqlTp2L2NzY24siRI5wlkbKO1+vF97///Zj9w8PDuPDCC3mjNIUY6DmkoaHB8uvrkiVL0N6ecGAuka1aW1stb5SePHkSTU1NDPUUYaDniIaGBstJtubOncuReJQTOjs74w5AevDBBzNcjZkY6Dlg7dq1lmHucrlw8uRJGyoimp54o0qPHTvGEaUpwEDPcn6/Hx0dHZbHXnrppQxXQzRzO3bssGx+4epHM8dAz2J+vx93f/e7lsd8Ph9vgFLO6uzsjBvqHHg0fQz0LOX3+9HU1IRPhodjjjU3N7NrIuW8eKH+6KOP8ibpNHH63Cw1Z84cnDlzJmZ/Y2Mje7SQUaxu+IsItm/fzm+hFmY8fa6IXC8i+0XkgIg8EOeclSKyW0T2iMjrMyk439XV1VmGeXNzM8OcjNPZ2YkVK1ZE7FNV3HzzzTZVlLsmDXQRcQJ4HMBqAG4At4qIO+qceQCeAPDnqnoFgG+moda80NLSYjkNbmNjI4dJk7EeeeQRiESuR9/f3w+32x3nEWQlmSv0WgAHVPV9VQ0CeA7Amqhz/gKhRaI/AABVPZLaMvOD3+9HW1tbzP7q6mpemZPRPB4P7r///pj9+/btQ2lpqQ0V5aZkAn0hgIMTtvvD+ya6HMBFIrJNRN4WkXWpKjBfBAIBNDU1xexfuHAh5zOnvNDa2mrZR/3UqVOoqKiwoaLck0ygi8W+6DupBQCuAfAVAA0A/lpELo95IhGviOwSkV0DAwNTLtZkt912m+X+X/7ylxmuhMg+O3bsgMvlitl/+PBhBAIBGyrKLckEej+AyyZsLwLwkcU5r6jqGVU9CuANAFdFP5Gq+lW1RlVrrD60fLV27Vr09PTE7K+vr+ddfso78QbMxbvooT9KJtDfAlAlIktFpAjALQC2RJ3zEoAvikiBiMwGUAfAeuUFihBvJKjL5UJnZ6cNFRHZy+PxoLu7Gw5HZDz19PRw0NEkJg10VR0BcBeAToRC+nlV3SMi60VkfficfQBeAfAOgJ0AnlLV36WvbHNs2LDBcj+H9VM+83g8ePLJJ2P2c8BRYhxYZCO/34/1TU0xNyS44hBRSFVVFQ4cOBCxr76+Pq+/vc54YBGlXktLC5qiwlxEGOZEE2zevDmmf3pXVxfXJI2DgW4Dq/7mIoKNGzcyzIkmiNc/vaOjg71eLLDJxQZXXHFFTN9yl8uFI0c4HovIitV8L1VVVXjvvfdsqsg+bHLJIoFAwHKg0O23325DNUS5obOzE+Xl5RH7enp6OH96FAZ6hln1pV22bBnnaSGaxI9//OOYfV1dXez5MgEDPYP8fr/lAKInnnjChmqIcovX60VVVVXM/u9973s2VJOdGOgZ9NBDD8XsW7ZsGUeDEiXp6aefjtk3NDTEXi9hDPQMWbt2LQ4fPhyzn1fnRMnzeDxobGyM2f/ss8/aUE32YaBnQCAQsBzez3VBiaauvb0dhYWFEfvGxsZQV1dnU0XZg4GeAZs3b47ZV15ezj7nRNNk1W6+c+fOvL9BykDPgNdefjlmn9UdeyJKTmtrK6qrq2P2W92nyicM9DRraWlBT19fxL7KykpenRPN0N69e1FUVBSxL9/nTWegp1EgELBcUu4HP/iBDdUQmeeee+6J2Wf1dy5fMNDTyGrV8sbGRl6dE6VIa2trzAjS3bt321SN/RjoaeL3+9Hf3x+xr6CggIs9E6VY9P2o3t7evL05ykBPkxdeeCFmX1lZmQ2VEJnN6/XG3CDN12bNpAJdRK4Xkf0ickBEHkhw3r8TkVER+UbqSsxNxy0GEbFnC1F6RLelDw4OorKy0p5ibDRpoIuIE8DjAFYDcAO4VUTccc5rRWipurxWV1eHne+8M75dWlrKhSuI0sjr9aKgoCBiX19fX96tQZrMFXotgAOq+r6qBgE8B2CNxXn/EcALAPJ6Um+/34+dO3fG7GeYE6XX0qVLY/ZZrUtqsmQCfSGAgxO2+8P7xonIQgA3AdiYutJyk9Wiz5/61KdsqIQov1hN3BUMBm2oxD7JBLpY7Ite5mgDgBZVHU34RCJeEdklIrsGBgaSrTGnfBA1iAjgBFxEmeDxeFBfXx+xL7pLo+mSCfR+AJdN2F4E4KOoc2oAPCcivQC+AeAJEbkx+olU1a+qNapa43K5plly9mppacGZoaGIfZWVlZyAiyhDOjs7sWTJkvHtvr6+vFrVKJlAfwtAlYgsFZEiALcA2DLxBFVdqqqVqloJ4H8A+LaqvpjyarPcM888E7MvX7tPEdnl1KlTEdtdXV15Mx3ApIGuqiMA7kKo98o+AM+r6h4RWS8i69NdYK4IBAL4MGogUW1tLW+GEmXY6tWrY/bly3QABZOfAqjqVgBbo/ZZ3gBV1W/NvKzc09bWFnNj4cYbY1qdiCjN2tvbsWXLlogr9TfeeMPGijKHI0VT5KOPIm8rOBwOrFy50p5iiPLcqlWrIrYHBwfzok86Az1F7vjCFyK277vvPt4MJbJJc3NzzL6f//znNlSSWUk1uVBigUAAv9m6FSscDnxy9dW448472XZOZCOPx4M5c+bg9OnT4/s++eQTGyvKDF6hz1AgEMAXv/hFbNy/H2+MjeG3776LK6+80u6yiPLet7/97Yjts2fPGj8LIwN9htra2jA6+sfxVMPDw9i2bZt9BRERgNBc6ZUT+qQD1rOgmoSBPkPRk+mLCG+GEmWJHzz4YMS2iQMaJ2Kgz1B0u9zixYt5M5QoS3i9XjR+9avj2x0dHUY3uzDQZ8Dv9+Nw1LznD0ZdERCRvf7twIGI7Z/85Cc2VZJ+DPQZiG6Pc7vd7N1ClGVEIucX7O3tNXYqAAb6DBw/fjxie9GiRTZVQkTx3H333TH7TJ0KgIE+TYFAAG+99VbEvvfff9+maogoHq/XG7McXfTIblMw0Kdp27ZtUI2cveVrX/uaTdUQUSJfiBrJXVVVZVMl6cVAn6YFCxZEbNfX16O1tdWmaogokZ6enojt7du321RJejHQp+nYoUPj//E4ERdRdrv00ksjtnt7e43svshAn6aVLheKATgdDhQXFzPQibKY1WRdv/jFL2yoJL2SCnQRuV5E9ovIARF5wOJ4o4i8E/7pFpGrUl9qdnn31VdxJYCvXn89XnvtNQ4mIspiHo8HCxdGrG2PwcFBm6pJn0kDXUScAB4HsBqAG8CtIuKOOu0PAL6kqn8K4G8AmPddZgK/34+mF1/ETgAvbt2Kd9991+6SiGgSZ86cidj+8MMPbaokfZK5Qq8FcEBV31fVIIDnAKyZeIKqdqvqx+HNNxFaSNpY0QOKTJ/wh8gE8+bNi9g2cfbFZAJ9IYCDE7b7w/viuQPAyzMpKtstW7w4YvvrX/+6TZUQUbKsFmw37WIsmUAXi33Ry2eGThT59wgFuuVaTyLiFZFdIrJrYGAg+SqzSCAQwM/+8R8BhIYUNzc3c7g/UQ7wer1obGyM2GfaxVgygd4P4LIJ24sAxAyzEpE/BfAUgDWqeszqiVTVr6o1qlqTq9NYtrW14Vx4/nNVxXvvvWdzRUSUrBUrVthdQlolE+hvAagSkaUiUgTgFgBbJp4gIosB/ArAX6qq0Qm3f//+hNtElL2im1hM67o4aaCr6giAuwB0AtgH4HlV3SMi60Vkffi0HwJYAOAJEdktIrvSVrHNor9Z5Oo3DaJ8FN3E8pvf/MaomReT6oeuqltV9XJV/bSqPhzet1FVN4Z//ytVvUhVl4V/atJZtJ3mz50bse12R/fgJKJs5fV6ceOf/dn49rlz54xaMpIjRacgEAjgX195ZXy7sLAQ69ats7EiIpqqkgsvjNjes2ePTZWkHgN9CjZv3oyRCQtCf+UrX+EIUaIc8/rrryfczmUM9BkoLy+3uwQimqKysrKE27mMgT4FpaWlCbeJKPtde+21CbdzGQN9Cnbv3p1wm4iy37p161DkdEIAFBUVGXUfjIE+BcuWLUu4TUTZz+PxYNvf/i3WAFhWVWXU5HoFdheQS06ePDn+u4jETPZDRLnh3aNH8SIA7NmDnU1NAGDEFB68Qk9SIBDApk2bxreLioq4qAVRjnohajCRKZN0MdCTtG3bNpw7d258+/bbb2eXRaIcZWrzKQM9ScePH4fqHyeZZA8Xotw175JLxqeRdTgcxjSfMtCTxB4uROZY+eUvYxYAp4hRawIz0JMUPamPafMoE+UTj8eDDRdeiFUVFdiwYYMxzafs5ZIkr9cL7N2LFx57DF+/7z4j7ogT5atAIIB7Tp5E8MQJ/J977sGVV15pRKjzCn0KvB4POgF4v/Utu0shohnYtm0bgqoYBRAMBo2ZcZGBPhWDg6E/58+3tw4impGVK1eiyOmEE2Z1QWagT0Fg9278BECAy84R5TSPx4PXbrgBfzN3Ll577TUjmluAJANdRK4Xkf0ickBEHrA4LiLy9+Hj74jI1akv1V6BQACrNm3CXwNYtXq1UaucEOUjT3U1fhAMwpNPk3OJiBPA4wBWA3ADuFVEopfpWQ2gKvzjBfBkiuu03bZt2xAcGTGuzY0oby1YAAwPA0NDdleSMsn0cqkFcEBV3wcAEXkOwBoAeyecswbAZg2NvHlTROaJSIWqHkp5xcGgLR/AymuuQZHDgeDYmFFtbkR56/y9sN5eYOHCzL52cTFwwQUpf9pkAn0hgIMTtvsB1CVxzkIAqQ/0l14Cbr455U87GQ+A1wBs+/SnsfKf/9mYNjeivHV+gffPfS7zr93SAjzySMqfNplAF4t9Oo1zICJehJpksHjx4iRe2sJVVwE/+9n0HjtDHgCeVauAK6+05fWJKIUaGgCfz54ml5qatDxtMoHeD+CyCduLAHw0jXOgqn4AfgCoqamJCfykXH556IeIaCZmzQIMGyCYTC+XtwBUichSESkCcAuALVHnbAGwLtzb5VoAJ9LSfk5ERHFNeoWuqiMicheATgBOAJtUdY+IrA8f3whgK4AbABwAMATg9vSVTEREVpKay0VVtyIU2hP3bZzwuwL4TmpLIyKiqeBIUSIiQ8jERRsy+sIiAwD6bHnxmSkDcNTuIjKM7zk/8D3nhiWq6rI6YFug5yoR2aWq6elzlKX4nvMD33PuY5MLEZEhGOhERIZgoE+d3+4CbMD3nB/4nnMc29CJiAzBK3QiIkMw0CchIt8UkT0iMiYice+GT7YISC4Rkfki8qqI9IT/vCjOeb0i8q6I7BaRXZmuMxXycfGWJN7zShE5Ef5cd4vID+2oM5VEZJOIHBGR38U5bsbnrKr8SfADoBrAnwDYBqAmzjlOAL8H8CkARQB+C8Btd+0zeM9tAB4I//4AgNY45/UCKLO73hm8z0k/N4SmtHgZoRlFrwWww+66M/CeVwL4V7trTfH7XgHgagC/i3PciM+ZV+iTUNV9qrp/ktPGFwFR1SCA84uA5Ko1AJ4O//40gBttrCWdkvncxhdvUdU3AcwTkYpMF5pCpv2/mhRVfQPAYIJTjPicGeipEW+Bj1x1iYZnywz/eXGc8xRAl4i8HZ7rPtck87mZ9tkm+348IvJbEXlZRK7ITGm2MuJzTmpyLtOJyP8CUG5x6D+p6kvJPIXFvqzuPpToPU/hab6gqh+JyMUAXhWR/xu+EsoVKVu8JYck837+DaHh5adF5AYALyK0XrDJjPicGegAVPW6GT5FUgt8ZJNE71lE/t/5NWHDXzuPxHmOj8J/HhGRXyP0dT6XAj1li7fkkEnfj6qenPD7VhF5QkTKVDXX5jyZCiM+Zza5pEYyi4Dkki0Abgv/fhuAmG8pIlIiInPP/w6gHoBlD4Islo+Lt0z6nkWkXEQk/HstQjlxLOOVZpYRnzOv0CchIjcB+O8AXAD+p4jsVtUGEbkUwFOqeoPGWQTExrJn6hEAz4vIHQA+APBNAJj4ngFcAuDX4b/3BQCeUdVXbKp3WuJ9biYv3pLke/4GgP8gIiMAzgK4RcNdQXKViDyLUO+dMhHpB/AQgELArM+ZI0WJiAzBJhciIkMw0ImIDMFAJyIyBAOdiMgQDHQiIkMw0IkAiMhlIvIHEZkf3r4ovL3E7tqIksVAJwKgqgcBPIlQH3yE//Srap99VRFNDfuhE4WJSCGAtwFsAnAngM+HZyQkygkcKUoUpqrnROR+AK8AqGeYU65hkwtRpNUADgH4nN2FEE0VA50oTESWAfgyQivWfC8XFzig/MZAJ0JoTUmEboreo6ofAHgUwH+1tyqiqWGgE4XcCeADVX01vP0EgM+KyJdsrFqP9UAAAAA6SURBVIloStjLhYjIELxCJyIyBAOdiMgQDHQiIkMw0ImIDMFAJyIyBAOdiMgQDHQiIkMw0ImIDPH/ATPcs0w0OL7hAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The calculated area is A =  3.1401597874737017264\n",
      "The true area is A =  3.141592653589793116\n"
     ]
    }
   ],
   "source": [
    "# -------------- Edit these --------------------------------------#\n",
    "trials = 10 # Modify the number of trials to n = 10, 100, 10000.\n",
    "radius = np.longdouble( 1 ) # Radius of the circle\n",
    "flag = 2 # experiment 1 or experiment 2\n",
    "#-----------------------------------------------------------------#\n",
    "c = np.longdouble(radius+radius/10);\n",
    "\n",
    "lims = np.longdouble(2*radius+radius/2);\n",
    "X = np.longdouble(np.arange(-lims/2,lims/2,0.01))\n",
    "\n",
    "G = np.longdouble(np.sqrt(radius**2-X**2));\n",
    "G[np.isnan(G)] = 0\n",
    "\n",
    "fig, axs = plt.subplots(1)\n",
    "fig.suptitle('Circle')\n",
    "axs.plot(X, G,'r-')\n",
    "# initialize some arrays and a list\n",
    "x = np.empty(trials)\n",
    "u = np.empty(trials)\n",
    "#g = []\n",
    "g = 0\n",
    "\n",
    "if (flag ==1):\n",
    "    f_x = 0\n",
    "    for ii in range(0,trials):\n",
    "        x[ii] = np.longdouble(lims*(random.uniform(0, 1)-0.5))\n",
    "        f_x = np.longdouble(np.sqrt(radius**2-x[ii]**2)) # This can be any function/pdf\n",
    "        u[ii] = np.longdouble(c*random.uniform(0, 1))\n",
    "        if (f_x>u[ii]):\n",
    "            #g.append(x[ii])\n",
    "            g = g+1\n",
    "            axs.plot(x[ii],u[ii],'g.')\n",
    "        else:\n",
    "            axs.plot(x[ii],u[ii],'k.')\n",
    "    \n",
    "    counts = len(g)\n",
    "    calc_area = 2*(counts/trials)*(c*lims)\n",
    "    true_area = np.pi*radius**2\n",
    "else:\n",
    "    f_x = np.empty(trials)\n",
    "    for ii in range(0,trials):\n",
    "        x[ii] = 2*radius*(random.uniform(0, 1)-0.5)\n",
    "        f_x[ii] = np.sqrt(radius**2-x[ii]**2) # This can be any function/pdf\n",
    "        axs.plot(x[ii],f_x[ii],'k.')\n",
    "        \n",
    "    calc_area = 2*(np.sum(f_x)*2*radius/trials)\n",
    "    true_area = np.pi*radius**2\n",
    "\n",
    "axs.set_aspect('equal')\n",
    "axs.set_xlabel('X')\n",
    "plt.show()\n",
    "\n",
    "# Calculations\n",
    "print('The calculated area is A = ',calc_area)\n",
    "print('The true area is A = ',true_area)\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.8.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}