Lesson5-2_sedov_exercise.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.2, 4)"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import numpy as np, matplotlib.pyplot as plt\n",
    "d_air = 1.25 # Kg/m3 air density\n",
    "ton = 4.184e9 # Joule/ton\n",
    "\n",
    "0.2, 4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.linspace(0,10)\n",
    "fit, cov = np.polyfit(x, x, deg=2, cov=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "time, radius = np.loadtxt(\"course_materials/sedov.txt\").transpose()\n",
    "time=time/1000\n",
    "tl, rl = np.log10(time), np.log10(radius)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 0, 'log t')"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEKCAYAAAD9xUlFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xl8VOW9x/HPL2xJANnCGghBdlARibutC2oVEWv11oW6VZvrUpdaqr1SxQ1bl7qVKuJyES/a3qpVQKtiLaJVVPCKGvYtENn3JWxJfvePGWKIM0mIOTOZme/79cormTnPyfwOE/LNeZ5znsfcHREREYC0eBcgIiL1h0JBRETKKRRERKScQkFERMopFEREpJxCQUREyikURESknEJBRETKKRRERKRcw3gXcKCysrI8Nzc33mWIiCSUWbNmrXf3ttW1S7hQyM3NZebMmfEuQ0QkoZhZYU3aqftIRETKKRRERKScQkFERMopFEREpJxCQUREyikURESknEJBRETKKRREROq53SWlPDltMZ8v3xT4aykURETqsWnz1/LgpaM4+5xjGZjbBnJzYeLEwF4v4e5oFhFJBSs2FnP3lDlk/O0vPPD2GNL37g5tKCyE/PzQ18OH1/nrKhREROqRXXtDXUVj319MgzTjk89e+jYQ9ikuhpEjFQoiIsnK3XlnzhrumTKHok07OXtAJ24b0ofm966KvMPy5YHUoTEFEZE4WzP2OTa07cRph3Tilfsv5p32RfzpooF0bJEBOTmRd4r2/PekUBARiZPtu0t4/ZYHaX79tWRtWE0aTvtNa+j1u5u/HUwePRoyM/ffMTMz9HwAFAoiIjHm7rz+xTcM/uM0Bj31IJklUcYMIDRuMG4cdO0KZqHP48YFMp4AGlMQEYmpuau2MmpSAZ8u3chhnVuQvW195IYVxwyGDw8sBCpTKIiIxMCWnXt5ZOoCJny8jBYZjfj9Tw7lgrwu2EM5octMKwtozKA6CgURkQCVlTkvzyri/rfmsal4Dz87pis3n9aLlpmNQw1Gjw7dd1Bc/O1OAY4ZVCewUDCzLsAEoANQBoxz98cqtRkO3Bp+uB24xt1nB1WTiEgszV6xmTsmFTB7xWbyurZiwjlH0b9Ti/0b7esWGjky1GWUkxMKhBh1F1UW5JlCCfBrd//czJoDs8xsqrvPqdBmKXCiu28yszOBccDRAdYkIhK4Ddt38+Db8/nrzBVkNWvCIxcM4MeHZ2NmkXeI4ZhBdQILBXdfBawKf73NzOYC2cCcCm0+qrDLDKBzUPWIiAStpLSMFz9dzkNvz6d4TylXndCNGwb3pHl6o3iXVmMxGVMws1xgIPBJFc2uBP4Ri3pEROrap0s3csfrXzNv9TZO6JHFncP60aNd83iXdcACDwUzawa8Atzk7lujtDmZUCicEGV7PpAPkBOnEXkRkUjWbN3F79+cy2tfrKRTi3SeHH4EZxzSIXpXUT0XaCiYWSNCgTDR3V+N0uYw4BngTHffEKmNu48jNN5AXl6eB1SuiEiN7SkpY/xHS3ns3YXsLXOuP6UH157Ug4zGDeJd2vcS5NVHBjwLzHX3h6O0yQFeBS5x9wVB1SIiUpc+WLiOOycVsHjdDk7t247bh/aja5um8S6rTgR5pnA8cAnwlZl9EX7uNiAHwN3HAncAbYAnwqdaJe6eF2BNIiK1VrSpmLdve5gfvfg4U7euZ3enbDIe+AO0OTLepdWZIK8++hCoslPN3a8CrgqqBhGRurBrbylPvb+E5WOe5t43HycjvL5BxsqiQBe8iQfd0SwiEoW78+7ctdw9pYAVG3cy698vlAdCuQAXvIkHzZIqIhLB0vU7uGL8Z/xiwkzSGzbgxauOps2GNZEbB7TgTTzoTEFEpIIdu0sY869FPPvBUpo0TOP2of249NiuNGqQFpqCoh5NXhcEhYKICKGuoilfrmL0G3NZvXUX5x3RmVvP7E275unfNqpnk9cFQaEgIilv/uptjJr0NTOWbKR/p4P48/CBDOra+rsN69nkdUFQKIhIytqycy+PvruACR8X0jy9IaPPPYQLj8yhQVoVF07Wo8nrgqCBZhFJOWVlzsf3PE5xx87cPuxQZj5zFR92XsXwo7tWHQgpQGcKIpJSvizazNTfPcI1E+8vXxu51bpVcP210KRhUp8F1ITOFEQkJWzcsYf/evUrzvnzv7n4tbHlgVBu3/0GKU5nCiKS1ErLnBc/KeShdxawfXcJVx7fjQ73r4vcOInuN6gthYKIJK2ZyzZyx+sFzFm1leO6t+HOYf3p1b55StxvUFsKBRFJOmu37uIP/5jHq//3DR1bpPPni49gyKEV1jhIgfsNakuhICJJY29pGeP/vYzH/rmQPSVl/PLkHlx7cncyG1f6VZcC9xvUlkJBRJLChwvXc+fkAhat3c4pfdpxx9B+5GZVscZBkt9vUFsKBRFJaN9s3sm9U+bwj69X07VNJs9elsfgvu3jXVbCUiiISELatbeUp6cv4c/TFgEw4vReXPWDg0lvlNjLYcabQkFEEs4/567hrslzWL6xmCGHdmDkWf3IbpkR77KSgkJBRBLG0vU7uHtyAf+av44e7Zox8aqjOb5HVrzLSioKBRGp94r3lPDnfy3i6elLadwwjd+d1ZfLjssNrXEgdUqhICL1lrvzxlehNQ5WbdnFT47I5rdn9KHdQenV7yy1opgVkXpp1RPPsj6rE0MGdOa1B4fzXodvePinhysQAqZQEJF6Zeuuvbw64gFa3HgdbTeuJg2n/aY1HDzyVzBxYrzLS3oKBRGpF8rKnJdnFXHKQ+9z1NMPaRbTOAksFMysi5n9y8zmmlmBmd0YoY2Z2eNmtsjMvjSzI4KqR0Tqr6+/2cL5Yz9ixN9m06V1Btnb1kduqFlMAxfkQHMJ8Gt3/9zMmgOzzGyqu8+p0OZMoGf442jgyfBnEUkBm3bs4aF35vPip8tp07QxD55/GOcd0Rl7QLOYxktgoeDuq4BV4a+3mdlcIBuoGArnABPc3YEZZtbSzDqG9xWRJFVa5rz06XIeemc+23aVcMVx3bjptJ4clN4o1ECzmMZNTC5JNbNcYCDwSaVN2cCKCo+Lws8pFESS1KzC0BoHBSu3cszBrblr2CH07tB8/0aaxTRuAg8FM2sGvALc5O5bK2+OsItH+B75QD5Ajk4fRRLS2m3hNQ4+/4YOB6Xzp4sGMvSwjt+ucVCZZjGNi0BDwcwaEQqEie7+aoQmRUCXCo87AysrN3L3ccA4gLy8vO+EhojUX3tLy3j+o2U8+u5CdpeUcu1J3bnu5B40baJ7Z+ujwN4VC8X/s8Bcd384SrNJwC/N7C+EBpi3aDxBJHl8tGg9oyYVsHDtdk7q3ZY7hvbj4LbN4l2WVCHIqD4euAT4ysy+CD93G5AD4O5jgTeBIcAioBi4IsB6RCRGVm7eyeg35vLGV6vo0jqDZy7NY3DfdtG7iqTeCPLqow+JPGZQsY0D1wVVg4jE1u6SUp75YClj3ltEmTs3n9aL/B9qjYNEok49EakT781bw92T57BsQzFn9O/AyLP60qV1ZrzLkgOkUBCR76Vwww7unjyHf85bS/e2TXnhyqP4Qc+28S5LakmhICK1snNPKU9MW8RT05fQKM24bUgfLj+uG40bakq1RKZQEJED4u784+vVjH5jLt9s3sm5A7P57Zl9aK8prZOCQkFEamzhmm3cObmAfy/aQJ8Ozfnf/zyWo7q1jndZUocUCiJSrW279vL4Pxfy3/9eRmbjBtx9Tn8uPiqHhloOM+koFEQkKnfntS++4b4357F++24uPLILI07vTZtmTeJdmgREoSAiERWs3MKo1wuYWbiJAV1a8syleQzo0jLeZUnAFAoisp/NxXv44zsLmPhJIa0yG/PAeYdx/qDOpKXpbuRUoFAQESC0xsFfP1vBg2/PY8vOvVx6bC6/Oq0XLTIaxbs0iSGFgojw+fJNjHq9gK++2cJR3Vpz17D+9O14ULzLkjhQKIiksHXbdvPAW/P426wi2h/UhMcuPJxhAzpp4roUplAQSUElpWVM+LiQR6YuYFdJKf954sFcf0pPmmmNg5SnnwCRFPPx4g3cOamA+Wu28YOeWdw5rD/dtcaBhCkURFLEqi07ue/NeUyevZLOrTJ46pJBnN6vvbqKZD8KBZEkt7uklGc/DK1xUFrm3HRqT64+sbvWOJCIFAoiSWza/LXcNXkOS9fv4PR+7bl9aD+tcSBVUiiIJKEVG4u5e8ocps5Zw8FZTXn+50dxYi+tcSDVUyiIJJGde0p58v3FjH1/MQ3TjN+e2YefH681DqTmFAoiiWriRBg5EpYvx7t04ctrb+Va78M3m3cybEAnbhvSlw4ttMaBHBiFgkgimjgR8vOhuBgAW76cnr+7mbMv+A0n3X0TxxzcJs4FSqLSOaVIIpk4EXJz4Wc/Kw+EfTJLdnPrBxMUCPK96ExBpL7b101UWAhm4B61qa1YEcPCJBkFdqZgZs+Z2Voz+zrK9hZmNtnMZptZgZldEVQtIglrXzdRYWHocRWBAEBOTvA1SVILsvtoPHBGFduvA+a4+wDgJOCPZtY4wHpEEkcV3URRZWbC6NGBliXJL7BQcPfpwMaqmgDNLXSPfbNw25Kg6hFJGJXPDmqia1cYNw6GDw+uLkkJ8RxTGANMAlYCzYEL3L0sjvWIxE+Fy0tJS4PS0prtl5mpMJA6Fc+rj34EfAF0Ag4HxphZxFU9zCzfzGaa2cx169bFskaR4FU8M3CvPhD2TWCnswMJQDxD4QrgVQ9ZBCwF+kRq6O7j3D3P3fPattWt+pIE9o0ZpKXBZZfVfNyga1d44YVQeCxbpkCQOhfP7qPlwGDgAzNrD/QGlsSxHpHYqHTjWY26itRNJDESWCiY2UuErirKMrMiYBTQCMDdxwL3AOPN7CvAgFvdfX1Q9YjUGyNH1uzMoEEDKCsLXWY6erQCQWIisFBw94uq2b4SOD2o1xept5Yvr76NzgwkTjTNhUiM7eqYHXlDgwahQWQNIEscKRREYmTFxmLyJ8zkliMuYFejJvtvzMyE558PdRdpAFniSKEgErBde0t59N0FnPrw+3ywcD29b76aBs8+HToj0JmB1DOaEE8kIO7O1DlruHvKHIo27WToYR25bUhfOrXMAHrAJZfEu0SR71AoiARg8brt3DV5DtMXrKNX+2a8+IujOa57VrzLEqnWAYeCmZ0G3OLupwVQj0hC27G7hD+9t4hnP1xCesMG3D60H5ce25VGDdRTK4khaiiY2SnAWELTULwG3AdMIHRPgaZiFKnA3Zn85SpGvzGHNVt3c/6gztx6Rh/aNm9S/c4i9UhVZwp/BPKBj4EzgRnA7e7+WCwKE0kU81ZvZdTrBXyydCOHZB/EE8MHMahrq3iXJVIrVYWCu/u08Nevmdk6BYLIt7bs3MsjUxfwwoxCmqc35L5zD+WCI7vQIM3iXZpIrVUVCi3N7CcVHlvFx+7+anBlidRfZWXOy7OKuP+teWwq3sPFR+fw69N606qp1oiSxFdVKLwPnB3lsQMKBUleFdc3qDD30JdFm7n99QJmr9jMoK6teH7YURyS3SLe1YrUmaih4O5aM1lSU+VZTAsL8V/k89fPlvNfGYfRpmkTHv7pAM4dmI2Zuookueg+BRGoduUz21nMD/77Ea6cOI0bT+1J8/RGcSpUJFgKBZEarm/Qadt6fje0XwwLE4k93VEjqakWK59ZTk7wdYnEWbVnCpWuQNpnC/CVu6+t+5JEAlbblc9G655NSX416T66EjgW+Ff48UmEbmTrZWZ3u/sLAdUmUneqGTOISCufSQqqSSiUAX3dfQ1AeD3lJ4GjgemAQkHqN62JLFJjNRlTyN0XCGFrgV7uvhHYG0xZInWohmsiu1Y+E6lRKHxgZlPM7DIzuwyYBEw3s6bA5mDLE6mFioPIublQWFj9PpmZmFY+E6lR99F1wE+AEwjNkPo88Iq7O3BygLWJHLgIN55hBu7faVrWoAFpGjMQ2U+1oeDubmYfAnsITW/xaTgQROqfSF1F7jiG8e2PrWdmkqYuIpHvqLb7yMx+CnwKnA/8FPjEzM4PujCRGqlhV5HjbMzqiIfHDEyBIBJRTbqPRgJH7rsnwczaAu8CL1e1k5k9BwwF1rr7IVHanAQ8CjQC1rv7iTUvXVLeAXQV7c3uTOuiFTEuUCTx1GSgOa3STWobarjfeOCMaBvNrCXwBDDM3fsD/1GD7ynyrahdRZWeysykyf1/iFlZIomsJr/c3zKzt83scjO7HHgDeLO6ndx9OrCxiiYXA6+6+/Jwe90dLdFV7iaaODF0I1oEDmxqq64ikdqoyUDzb8zsPOB4QlcfjXP3v9fBa/cCGpnZNKA58Ji7T6iD7yvJJlI3UX4+tG4NGzZ8p3lJdhdaFUUODBGpWo1mSXX3V4BXAnjtQcBgIAP42MxmuPuCyg3NLJ/QetHkaFKy1BOpm6i4mLL0dEqapNN4967ypz0zk8b3/z7GBYokj6jdR2a2zcy2RvjYZmZb6+C1i4C33H2Hu68nNGXGgEgN3X2cu+e5e17btm3r4KUloUTpJmLjJm750S/ZrK4ikTpT1cprzQN+7deBMWbWEGhMaC6lRwJ+TUlEOTkRLzXd0Lo91zx1By07PBiHokSSU2DrKZjZS8DHQG8zKzKzK83sajO7GsDd5wJvAV8Sug/iGXf/Oqh6JAFEGkwGGD0az8jcr2lJegZZjz9E7w5B/+0ikloCW3nN3S+qQZsHAf2ZJ1EHk0vKnPG5x7HgjF9y4z/H02nbOrxLFxred5+6iUQCoJXXpH6IMpi84cYR3PvGXNaefR57Fi/ByspIKyxUIIgERGs0S/0QZTC57aa1PHNpHoP7tsPMYlyUSOrRmYLEVrRxg2iXGud04dR+7RUIIjGiUJDY2TduUFgYmp9o301oEydScN0t7GrUZP/2mZmk3XdffGoVSVEKBYmdKOMG628cwVkbuvLH80ewq1NnrX4mEkcaU5DYiTJu0HrDGm4b0ofLjzuTxg3vjXFRIlKRzhQkdqKMG3iXzuT/sDuNG+rHUSTe9L9Q6k60QeSwVbfczu7G6fvvk5lJg99rriKR+kKhIHWjikHkbbv2cu+UOfygqCN3DL2B7e2zy+cq0riBSP1iibbccl5ens+cOTPeZUhlUZbC3NEhmxOvG8+GHbu58MgujDi9N22aNfnu/iISKDOb5e551bXTQLPUjSiDyBmrV5LdKoNnL8tjQJeWMS5KRA6UQkHqRpSZTHd26MTfrzmOtDTdfCaSCDSmIHWi7N7RlKRn7PecZ2TS9KH7FQgiCUShIN/brMJNDNuUw82nXsu61u2/XfDmaQ0iiyQahYJEV80lpuu27WbE32Zz3pMfsW7bbgbf8yuy1q/Cyspg2TIFgkgC0piCRBZlfQOAvRdexISPC3l06gJ2lZRy9Ynduf6UHjRtoh8nkUSn/8USWZR5inbd8luGre3MgjXb+UHPLO4c1p/ubZvFp0YRqXMKBYksyiWmjVd+Q/GeUp66ZBCna0prkaSjUJDIolxiur1dR969+UTSGzWIQ1EiEjQNNEtko0dTWukS07KMDA56+AEFgkgSUyjIdyzfUMxVJb341anXsrrlt5eYpj39tK4oEkly6j6Scjv3lPLktEWMnb6EhmnGDSOuofXfH8A0pbVIylAoCO7O2wWruWfKXL7ZvJNhAzpx25C+dGiRXv3OIpJUAvsT0MyeM7O1ZvZ1Ne2ONLNSMzs/qFpSXhU3oS1au51Ln/uUq//nc5qnN+Qv+cfw+EUDFQgiKSrIM4XxwBhgQrQGZtYAuB94O8A6UluUm9B2lZTxcNs8nvtwKRmNG3Dn2f342TFdadhAXUUiqSywUHD36WaWW02z64FXgCODqiPlRbkJbdONIxh39XNckNeF35zRmyytcSAixHFMwcyygXOBU6gmFMwsH8gHyImyzq9EEeUmtPZb1vHadcdzuNY4EJEK4tlX8Chwq7uXVtfQ3ce5e56757Vt2zYGpSWRKCFqOV0UCCLyHfEMhTzgL2a2DDgfeMLMfhzHepJOWZnz0c9vZmejSl1DmZnYfffFpygRqdfiFgru3s3dc909F3gZuNbdX4tXPcnmixWbOfeJf3NxcXfG/ey37MnuAuGb0BindQ5EJLLAxhTM7CXgJCDLzIqAUUAjAHcfG9Trprr123fzwFvz+N+ZRbRr3oTHLjycYQOGYHZnvEsTkQQQ5NVHFx1A28uDqiNVlJSW8T8zCvnj1AXs3FPKf/7wYK4f3JNmWuNARA6AfmMkgU+WbGDUpALmrd7GD3pmMers/vRopzUOROTAKRQS2Ootu7jvzblMmr2S7JYZjP3ZIH7UX2sciEjtKRQS0J6SMp79cCl/em8hJWXODYN7cs2J3clorCmtReT7USgkmPcXrOOuSQUsWb+D0/q15/az+pHTJjPeZYlIklAoJIgVG4u5Z8oc3pmzhm5ZTfnvK47k5N7t4l2WiCQZhUI9t2tvKWPfX8yT0xaTZsYtZ/TmyhO60aShuopEpO4pFOopd+edOWu4Z8ocijbt5OwBnbhtSB86tsiofmcRkVpSKNRDi9dt567Jc5i+YB292zfnpV8cw7Hd28S7LBFJAZo8vx7Z9fwLbGmfTbd2B/GHm89mfJOFvHHDCQoEEYkZhUI94O7MvG8M/otf0GLtStJwOm1Zy0kPjaThX16Kd3kikkIUCnE2d9VWLhg3gw73303G3t37bywuDi2SIyISIxpTiJMtO/fyyNQFvDCjkIPSG5K9bX3khlEWyRERCYJCIcbKypyXZxVx/1vz2FS8h+FHd+XXp/fCnssJrZ9cmVaaE5EYUijE0OwVm7ljUgGzV2xmUNdWPD/sKA7JbhHaOHo05Ofvv55yZmboeRGRGFEoxMCG7bt58O35/HXmCrKaNeHhnw7g3IHZ+09ct2/Rm5EjQ11GOTmhQNBiOCISQwqFAJWUlvHip8t56O35FO8p5aoTunHD4J40T28UeYfhwxUCIhJXCoWAfLZsI3e8XsDcVVs5vkcb7jy7Pz3bN493WSIiVVIo1LE1W3fx+zfn8toXK+nUIp0nhx/BGYd00BoHIpIQFAp1ZE9JGeM/Wspj7y5kb6lz/Sk9uOak7mQ21j+xiCQO/caqAx8sXMedkwpYvG4Hg/u04/ah/cjNahrvskREDphC4Xso2lTMvVPm8lbBarq2yeS5y/M4pU/7eJclIlJrCoVa2LW3lHHTl/DEtEUYxm9+FFrjIL2R1jgQkcSmUDgA7s67c9dyz5Q5LN9YzFmHdWTkkL50aqk1DkQkOQQWCmb2HDAUWOvuh0TYPhy4NfxwO3CNu88Oqp7va+n6Hdw1uYBp89fRs10zXrzqaI7rkRXvskRE6lSQZwrjgTHAhCjblwInuvsmMzsTGAccHWA9tVK8p4Qx7y3imQ+W0rhhGr87qy+XHZdLowaaYFZEkk9goeDu080st4rtH1V4OAPoHFQtteHuTPlyFfe9OZdVW3Zx3hGdufXM3rRrnh7v0kREAlNfxhSuBP4RbaOZ5QP5ADkxmDV0/uptjJr0NTOWbKR/p4MYc/FABnVtHfjriojEW9xDwcxOJhQKJ0Rr4+7jCHUvkZeX50HVsnXXXh6dupDnP15G8/SG3PvjQ7joqBwapOluZBFJDXENBTM7DHgGONPdN8SrjrIy55XPQ2scbNixh4uPymHE6b1p1bRxvEoSEYmLuIWCmeUArwKXuPuCeNXxVdEW7pj0Nf+3fDNH5LRk/BUV1jgQEUkxQV6S+hJwEpBlZkXAKKARgLuPBe4A2gBPhCeLK3H3vKDqqWzTjj08+M58Xvp0OW2aNuah/xjATwZmk6auIhFJYUFefXRRNduvAq4K6vWjKS3z8jUOtu8u4efHd+PGU3tyULQ1DkREUkjcB5pjaWZ4jYM5q7Zy7MFtuOuc/vTSGgciIuVSJhRe/byIm/93Nh1bpDPm4oGcdWhHrXEgIlJJyoTC4D7tuenUnuT/8GCtcSAiEkXK/HZskdmIm07tFe8yRETqNU3gIyIi5RQKIiJSTqEgIiLlFAoiIlJOoSAiIuUUCiIiUk6hICIi5RQKIiJSTqEgIiLlFAoiIlJOoSAiIuUUCiIiUi41QmHiRMjNhbS00OeJE+NdkYhIvZT8s6ROnAj5+VBcHHpcWBh6DDB8ePzqEhGph5L/TGHkyG8DYZ/i4tDzIiKyn+QPheXLD+x5EZEUlvyhkJNzYM+LiKSw5A+F0aMhM3P/5zIzQ8+LiMh+AgsFM3vOzNaa2ddRtpuZPW5mi8zsSzM7IpBChg+HceOga1cwC30eN06DzCIiEQR5pjAeOKOK7WcCPcMf+cCTgVUyfDgsWwZlZaHPCgQRkYgCCwV3nw5srKLJOcAED5kBtDSzjkHVIyIi1YvnmEI2sKLC46LwcyIiEifxDAWL8JxHbGiWb2YzzWzmunXrAi5LRCR1xTMUioAuFR53BlZGauju49w9z93z2rZtG5PiRERSUTxDYRJwafgqpGOALe6+Ko71iIikPHOP2GPz/b+x2UvASUAWsAYYBTQCcPexZmbAGEJXKBUDV7j7zBp833VAYS3LygLW13Lf+i5Zjy1ZjwuS99iS9bggsY+tq7tX29USWCjUR2Y2093z4l1HEJL12JL1uCB5jy1ZjwuS+9j2Sf47mkVEpMYUCiIiUi7VQmFcvAsIULIeW7IeFyTvsSXrcUFyHxuQYmMKIiJStVQ7UxARkSqkRCiY2QgzczPLirL9MjNbGP64LNb11YaZ3ROeXfYLM3vHzDpFaVcabvOFmU2KdZ0H6gCOKxHfswfNbF74+P5uZi2jtFtmZl+F/w2qvUw73g7guM4ws/nhmZF/G+s6a8PM/sPMCsyszMyiXnWUaO9ZVZK++8jMugDPAH2AQe6+vtL21sBMII/QNBuzwu02xbrWA2FmB7n71vDXNwD93P3qCO22u3uzmBdYSzU5rgR+z04H3nP3EjO7H8Ddb43QbhmQV/lntb6qyXGZWQNgAXAaodkMPgMucvc5sa73QJhZX6AMeAoYEe1eqkR7z6qSCmcKjwC3EGVeJeBHwFR33xj+pTKVqqf8rhf2/eIMa0qobRmUAAADvUlEQVT040soNTyuRH3P3nH3kvDDGYSmdkl4NTyuo4BF7r7E3fcAfyE0U3K95u5z3X1+vOuIpaQOBTMbBnzj7rOraJaws7Wa2WgzWwEMB+6I0iw9PJngDDP7cQzLq7UaHFfCvmcV/Bz4R5RtDrxjZrPMLD+GNdWFaMeVDO9ZVRL5PdtPw3gX8H2Z2btAhwibRgK3AadX9y0iPFcv/uqu6tjc/XV3HwmMNLP/An5JaCqRynLcfaWZHQy8Z2ZfufviAMuuVh0cV8K+Z+E2I4ESYGKUb3N8+D1rB0w1s3nh9Unipg6OK6Hfsxqod+9ZbSV8KLj7qZGeN7NDgW7A7NA0S3QGPjezo9x9dYWmRYTmaNqnMzAtkGIPULRji+BF4A0ihIK7rwx/XmJm04CBQFxDoQ6OK2Hfs/Cg+FBgsEcZ0Kvwnq01s78T6nqJ6y+YOjiuGs+KHGsH8PNY1feod+9ZbSVt95G7f+Xu7dw9191zCf1QHlEpEADeBk43s1Zm1orQmcXbMS73gJlZzwoPhwHzIrRpZWZNwl9nAccD9X1gr9rjInHfszOAW4Fh7l4cpU1TM2u+72tCxxZxnfP6oibHRWhguaeZdTOzxsCFhGZKTniJ+J5Vyd1T4gNYBmSFv84Dnqmw7efAovDHFfGutYbH8wqhH7wvgclAduVjA44DvgJmhz9fGe+66+K4Evg9W0SoX/2L8MfY8POdgDfDXx8cfr9mAwWEujDiXvv3Pa7w4yGErkBanAjHFa75XEJ/UO4mNNvz28nwnlX1kfSXpIqISM0lbfeRiIgcOIWCiIiUUyiIiEg5hYKIiJRTKIiISDmFgkg1zGx7HX2fy6PN+ipSXygURGLnckLXt4vUWwoFkRqykAfN7Ovw3PkXhJ9PM7MnwvPuTzGzN83s/Er7nk/oBryJ4Tn3M+JxDCLVSfi5j0Ri6CfA4cAAIAv4zMymE5o+JBc4FGgHzAWeq7iju79sZr+kijn5ReoDnSmI1NwJwEvuXurua4D3gSPDz//N3cs8NLfWv+JZpMj3oVAQqblI0z9X9bxIwlEoiNTcdOACM2tgZm2BHwKfAh8C54XHFtqz/7TeFW0DmsekUpFa0piCSM39HTiW0GyYDtzi7qvN7BVgMKHZXRcAnwBbIuw/HhhrZjuBY919Z0yqFjkAmiVVpA6YWTN3325mbQidPRzv3127Q6Te05mCSN2YYmYtgcbAPQoESVQ6UxARkXIaaBYRkXIKBRERKadQEBGRcgoFEREpp1AQEZFyCgURESn3/7ZLRi0VBW73AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "cc, cov = np.polyfit(tl, rl, deg=1, cov=True)\n",
    "plt.plot(tl, np.polyval(cc, tl))\n",
    "plt.plot(tl, rl, 'or')\n",
    "plt.ylabel(\"log R\")\n",
    "plt.xlabel(\"log t\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.sqrt(cov[0,0])\n",
    "a = cc[1]\n",
    "E = d_air * 10**(5*a)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2.286e+01\n"
     ]
    }
   ],
   "source": [
    "print (\"%0.3e\"%(E/(ton*1000)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[5.70694937e-05, 1.48995193e-04],\n",
       "       [1.48995193e-04, 4.16689828e-04]])"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cov"
   ]
  },
  {
   "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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}