phys512/HW1_solutions.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "51aef09e",
   "metadata": {},
   "source": [
    "# Homework 1 solutions"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "134cf3e5",
   "metadata": {},
   "source": [
    "1. **Long term planetary orbits**"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "19b91861",
   "metadata": {},
   "source": [
    "First write down the equations of the orbit that we need to solve.\n",
    "\n",
    "The acceleration is\n",
    "\n",
    "$$a = {GM\\over r^2}$$\n",
    "\n",
    "towards the Sun, where $M$ is the mass of the Sun and $r$ is the Earth-Sun distance.\n",
    "\n",
    "A circular orbit has $a = v^2/r$, so the velocity is given by $v^2 = GM/r$. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fb8a3016",
   "metadata": {},
   "source": [
    "First compute the orbit and make sure it looks like a circle:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "4caabbe4",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "nsteps = 1000\n",
      "Delta E/E =  -1.9724690657386198e-05\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjgAAAGdCAYAAAAfTAk2AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABr5klEQVR4nO3dd3xT5f4H8E9Gk8500E0LbaHQllVmKYKoVIao4LiCogwRrgMXOOD+BFT04rpeF4oi84qiqAwFy5IllFUos4zSlu7dJp1Jk5zfH2mjEeiApidNP+/X67xoT56cfk9Dej455znPIxEEQQARERGRHZGKXQARERFRS2PAISIiIrvDgENERER2hwGHiIiI7A4DDhEREdkdBhwiIiKyOww4REREZHcYcIiIiMjuyMUuQAxGoxE5OTlwc3ODRCIRuxwiIiJqAkEQUF5ejsDAQEilDZ+jaZcBJycnB8HBwWKXQURERDcgMzMTQUFBDbZplwHHzc0NgOkXpFKpRK6GiIiImkKj0SA4ONh8HG9Iuww49ZelVCoVAw4REVEb05TuJexkTERERHaHAYeIiIjsDgMOERER2R0GHCIiIrI7DDhERERkdxhwiIiIyO4w4BAREZHdYcAhIiIiu8OAQ0RERHbHqgFn3759uOeeexAYGAiJRIKNGzc2+pw9e/agX79+UCqV6Nq1K1atWnVVmyVLliAkJASOjo6IiYnBkSNHWr54IiIiarOsGnAqKyvRp08fLFmypEnt09LSMHbsWNx+++1ISkrCCy+8gCeeeALbtm0zt/n+++8xe/ZsLFy4EMePH0efPn0watQoFBQUWGs3iIiIqI2RCIIgtMoPkkiwYcMGjB8//rptXn31VWzZsgVnzpwxr5s4cSLKysoQHx8PAIiJicHAgQPx2WefAQCMRiOCg4Px7LPPYu7cuU2qRaPRwN3dHWq1mnNRERERtRHNOX7b1GSbCQkJiIuLs1g3atQovPDCCwAAnU6HxMREzJs3z/y4VCpFXFwcEhISrrtdrVYLrVZr/l6j0bRs4UR004xGAdW1BlTq9KjSmv6t1hlQqTOgSqtHpc6Aap0etQYBBqOAWqMReoMAvVGA3mA0rTOYPq9JJYBMKoFUKjF9Lan/WgKZVAJHBxmcHGRwUkjh5CAzf+/oIIOTQgYXpRzuTg5wUciaNKkfEdkemwo4eXl58PPzs1jn5+cHjUaD6upqlJaWwmAwXLPN+fPnr7vdxYsX44033rBKzUR0bUajgOJKHfI1NcjX1KCgXIuSSh3KqnQorao1/1tapUNZ3ffGVjmf3HRyqQQqJwe4OzmY/3V3coCHkwN83JTwcVPCt+5fHzclvF2VcJDx3g0iW2BTAcda5s2bh9mzZ5u/12g0CA4OFrEiorZPXVWLjJIqZJRUIbO0Cjll1XVhRot8TQ0Ky7XQ30BikUgAF4UcznVnUpwcZHBRyuCsMH3tIJdCLpWYFpkEcqkUMqkEDjIJZFIpJBJTuDIKAgxGwCjUfy3AKAAGoxE1tUZU1xpQU7dU1xpQrTOY11fU6KEzGKE3Ciip1KGkUtfk+r1cFPB1U8Lf3REdPZwQ5OmMIE+nusUZ3q4KnhUiagU2FXD8/f2Rn59vsS4/Px8qlQpOTk6QyWSQyWTXbOPv73/d7SqVSiiVSqvUTGSvBMF0BialoAKXCytMQaYu0GQUV0FTo290GxIJ4O2qhL/KEb5uSni6KODp7AAPZwU8nRXwcvnza09nB7g5OsDRQSp6ABAEATW1Rqira6GuNp1dqv/a9H0tiiq0KCzXoqDc9G9RhdYiEJ3PK7/mtpVyKTp6OiHY0xmh3i7o4uOCMB9XdPFxhZ9KKfq+E9kLmwo4sbGx2Lp1q8W6HTt2IDY2FgCgUCjQv39/7Nq1y9xZ2Wg0YteuXZg1a1Zrl0tkFwRBQHZZNVIKKsxhJqWgApcKKlBWVdvgc33clAj2dEInL2cEejjB390RfirT4q9yhLerAvI2eMlGIpHASWHqj+Pv7tik5xiNAkqrdCis0KJAo0WuuhpZpaYlu7QaWaVVyNXUQKs3IrWwEqmFldh7sdBiGy4KGUJ9XBDmbQo83f1dEeGvQicvZ0ilDD5EzWHVgFNRUYGUlBTz92lpaUhKSoKXlxc6deqEefPmITs7G2vWrAEAPPnkk/jss8/wyiuv4PHHH8fvv/+OH374AVu2bDFvY/bs2ZgyZQoGDBiAQYMG4aOPPkJlZSWmTZtmzV0hsgu1BiMuF1bgbLYGZ3M0OJujxrlcDcqvczZGIgGCPZ3RxccFId4u6OTljGBPZ3TqYLrs4qywqc9IopJKJejgqkQHVyUirnNCWac3Ik9dg6wy01mw1KJKpBZW4HJhJTJKqlCpM+BMtgZnsi1vhHBRyNDd3w2RASrzEuHvBhclf/9E12PVd8exY8dw++23m7+v7wczZcoUrFq1Crm5ucjIyDA/Hhoaii1btuDFF1/Exx9/jKCgIHz99dcYNWqUuc2ECRNQWFiIBQsWIC8vD9HR0YiPj7+q4zFRe2cwCriYX47jGaU4nWUKMufzyqHTG69qK5dKEOLtgnBfV3T9yxLm7QonhUyE6u2TQi5Fpw6mgDiki+VjOr0RGSVVSC2sQGpRJVIKKnA+T4OL+RWo1BlwPKMMxzPKzO0lEiDU2wXRQR6I7uSBPkEeiAxQQSFve2fMiKyh1cbBsSUcB4fsUUmlDicySnE8oxQnMspwMrMMlTrDVe1clXJEBagQFahCj0DTv+G+bjww2ii9wYi0okqcy9UgObccybkaJOdqUFCuvaqtQiZFVKAK0cEeiA72QN9OHujk5cx+PWQ3mnP8ZsBhwKE2KrOkCodSi3EotQSJV0qQXlx1VRtXpRzRwR7oE+yOHoHu6BGoQrAn+3PYg+IKLU5lq5GUUYaTWaZAW3qNPlN+KiViQjtgUKgXYkK90NXXlYGH2iwGnEYw4FBblF1WjUOXi5GQWoxDqcXIKq2+qk0XHxf06+SJfp090beTB8J93SBjmGkXBEFARkkVkjLLzMuZbLV58MN6Xi4KDArxwqBQLwwO64AIfzcGXmozGHAawYBDbYG6uhYHU4qw92IhDl4uRkaJ5RkauVSC3kHuGBxm+nTeN9gT7s4OIlVLtqim1oATGWU4nFaMI2klOJ5Rippayz5Y3q5KDAv3rlt84OPGITXIdjHgNIIBh2yRIAg4m6PB3ouF2HuhEIkZpTD8ZaA8mVSCXh1NgSa2SwcM6OzJu2ioWXR6I05nl+FwWgkOp5bgaHoJqv7WTysyQIVbu3nj1nAfDAjxhFLOTuZkOxhwGsGAQ7aiQqvHngsF2H2+EPsuFaLwbx1Hu/i4YHg3XwwL98aAEE+4OfIMDbUcnd6IxCul2H/J9P/v77enOznIcGs3b9wZ5Y87Inzh5aIQqVIiEwacRjDgkJiKKrTYeS4f287m4UBKMXSGPy8ZOCtkGNLFG8O7++C2bj4I9nIWsVJqb4ortPgjpQj7LhZh/6VCizu1pBJgQGcv3Bnlhzuj/BDi7SJipdReMeA0ggGHWltmSRW2n8vHtjN5OHalxGJSyVBvF8RF+uK27r68JEA2o/6S6fZz+dhxLh/JuZZnd7r6umJUDz+M7RWIyAA33plFrYIBpxEMONQa8tQ1+OVkDjafzMHpbLXFYz07qjAqyh+je/rztl1qEzJLqrAzOR87k/NxOLXEYiLVLj4uuLt3IO7pE4Cuvm4iVkn2jgGnEQw4ZC1lVTpsPZ2HzSezcTitBPXvLqkEGBjihVE9/DGyhx+CPHnpidoudVUt9lwswJZTudhzsdBidOwIfzfc0ycQd/cOQOcOvIxFLYsBpxEMONSSqnUG7EjOx+akbOy9WGgx7sjAEE/cG90Rd/X0RwdX3n5L9kdTU4sdZ/Px66kc7L9UZHFmJzrYAw/2D8I9fQLh7sQO8nTzGHAawYBDN0sQBJzMUuOHY5n4JSkH5do/J6uMDFBhXHQg7ukTiI4eTiJWSdS6yqp02HY2D7+czMXBy0XmvmYKuRSjevjjwf5BGNrVm4NP0g1jwGkEAw7dqJJKHTacyMYPRzNxIb/cvD7
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkAAAAGdCAYAAAD60sxaAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABZiUlEQVR4nO3deVxU5eIG8GdmYGbYZhBZBhQE3HAhcUVcUpOisrqWlZq5XdOuaWV2Ky2vdqub3W7br7JMy7RSK8vMzChzVxAUwVxwB0FwWERmWGTW8/sDmaJAQRnOLM/385mPn4b3DM9Mwjy+c973SARBEEBERETkRqRiByAiIiJqbSxARERE5HZYgIiIiMjtsAARERGR22EBIiIiIrfDAkRERERuhwWIiIiI3A4LEBEREbkdD7EDOCqr1YrCwkL4+flBIpGIHYeIiIiaQBAEVFRUICwsDFJp4/M8LECNKCwsRHh4uNgxiIiI6Drk5+ejffv2jX6dBagRfn5+AGpfQJVKJXIaIiIiagq9Xo/w8HDb+3hjWIAaUfexl0qlYgEiIiJyMtc6fYUnQRMREZHbYQEiIiIit8MCRERERG6HBYiIiIjcDgsQERERuR0WICIiInI7LEBERETkdliAiIiIyO2wABEREZHbYQEiIiIit8MCRERERG6HBYiIiIjcDi+GSkREDktXbcLJ4gpodTUo0tegvNoEs1WAxWqF0lMGtZcn2njLEdHWG52CfNHGRy52ZHISLEBEROQwzl+qxq6TpdhzugS/ndfh/KXLzTo+0FeB/pFtEB8VgCGdA9Ep2M9OScnZsQAREZGotLoabMgqwIbMAhzXVvzl6+38vdCujRdCVEoEeHvCUyaFTCpBjcmC8ssmXKw0Iqe0CgXll1FaacBPR7T46YgWANAlxBd3xobi3t7t0KGtT2s/NXJgEkEQBLFDOCK9Xg+1Wg2dTgeVSiV2HCIilyIIAtJzyrB8dw62Hi9C3TuRTCpBnwh/3Nw5CP0iA9A9TAW1l2eTHrPKYEb2BT3Scsqw7+xF7Dt7ESZL7QNLJMCwLkGYnBCJYV2CIJVK7PXUSGRNff9mAWoECxARUcsTBAE7T5bg7V9P4VB+ue3+/pFtcG/v9hgVGwq1d9MKz7XoLpuw5VgRvs8qwO5Tpbb7YzR+eOrWLritewgkEhYhV8MCdINYgIiIWtah/HK89tNxpJ69CACQe0gxpk97TBsShU7Bvnb93jmlVfhi3zl8tT8flQYzAKBnOxUW3d0D/SMD7Pq9qXU19f27VZbBL1myBJGRkVAqlYiPj0d6evpVx69btw4xMTFQKpWIjY3F5s2b631dEAQsXLgQoaGh8PLyQmJiIk6dOlVvzD333IOIiAgolUqEhoZi4sSJKCwsbPHnRkREV6e7bMKCDYcx+oO9SD17EXKZFI8MiULKvFuw+L5Yu5cfAIgK9MG/7uqO3c+OwKwRHeEjl+FIgR4PLE3FnC8zUaSvsXsGcix2L0BfffUV5s6di0WLFuHgwYPo1asXkpKSUFxc3OD4lJQUjB8/HtOmTUNmZiZGjx6N0aNH48iRI7Yxr7/+Ot59910sXboUaWlp8PHxQVJSEmpqfv8LPGLECHz99dc4ceIEvv32W5w5cwb333+/vZ8uERH9webDF5D41k58sS8PggDc17sdtj8zHAvu6o5AX0Wr52njI8czSTHY9ewIjB8QDokE2JBViMS3duLbjPPghyLuw+4fgcXHx6N///54//33AQBWqxXh4eF4/PHHMW/evL+MHzt2LKqqqrBp0ybbfQMHDkRcXByWLl0KQRAQFhaGp59+Gv/85z8BADqdDiEhIVi5ciXGjRvXYI6NGzdi9OjRMBgM8PS89ufL/AiMiOj6VRrMWLjhCNZnFgAAogN98J97Y5HQsa3Iyer77Xw5/rXhCA6d1wEAErsF49X7YhHspxQ5GV0vh/gIzGg0IiMjA4mJib9/Q6kUiYmJSE1NbfCY1NTUeuMBICkpyTY+JycHWq223hi1Wo34+PhGH7OsrAyrV6/GoEGDGi0/BoMBer2+3o2IiJrvUH45Rr27G+szCyCVAE/c0gk/zRnqcOUHAG5q749vZw7CM0ld4SmT4NfsYox6dw/SrpynRK7LrgWotLQUFosFISEh9e4PCQmBVqtt8BitVnvV8XV/NuUxn3vuOfj4+KBt27bIy8vD999/32jWxYsXQ61W227h4eFNe5JERGTz+b5zGPNhCs5drEY7fy98/WgC5t7WFQoPmdjRGuUhk2LWiE744fEh6Bzsi5IKAx76OA1Ld56B1cqPxFyVS18L7JlnnkFmZiZ++eUXyGQyTJo0qdHPd+fPnw+dTme75efnt3JaIiLnZTRb8cJ3h/GvDUdgtgq4M1aDzU8MRT8nWmEVo1Hh+9mDcW/vdrBYBbz203E8tvogLhstYkcjO7DrTtCBgYGQyWQoKiqqd39RURE0Gk2Dx2g0mquOr/uzqKgIoaGh9cbExcX95fsHBgaiS5cu6NatG8LDw7Fv3z4kJCT85fsqFAooFK1/Qh4RkbMrqzJi5hcZSMspg0QCPJsUg38Mi3bKPXa85R5468Fe6B8ZgBc3HkXyUS0uLEvF8sn9eF6Qi7HrDJBcLkffvn2xdetW231WqxVbt25tsIQAQEJCQr3xALBlyxbb+KioKGg0mnpj9Ho90tLSGn3Muu8L1J7rQ0RELaOg/DLuX5qCtJwy+Co88PGkfpg5vKNTlp86EokED8VHYPX0eLTx9sSh8zrcuyQFJ4v+epkOcl52/whs7ty5WL58OVatWoXs7GzMnDkTVVVVmDp1KgBg0qRJmD9/vm38k08+ieTkZLz55ps4fvw4XnzxRRw4cACzZ88GUPsXc86cOXjllVewceNGHD58GJMmTUJYWBhGjx4NAEhLS8P777+PrKwsnDt3Dtu2bcP48ePRsWPHq5YkIiJqulNFFRjzQQrOllQhTK3E+scGYWS3kGsf6CT6RwZg/WODERXoU1v0PkxBZt4lsWNRC7F7ARo7dizeeOMNLFy4EHFxccjKykJycrLtJOa8vDxcuHDBNn7QoEFYs2YNli1bhl69euGbb77Bhg0b0LNnT9uYZ599Fo8//jhmzJiB/v37o7KyEsnJyVAqa6cnvb29sX79eowcORJdu3bFtGnTcNNNN2Hnzp38mIuIqAVk5l3CAx+lQquvQadgX3z72CB0CXG9K69HBfpg/cxB6BPhD32NGQ9/nMYVYi6Cl8JoBPcBIiJqWGbeJUz8JB2VBjPiwv3x6ZT+aOMjFzuWXVUZzJj+2QGknLkIpacUH03sh2FdgsSORQ1wiH2AiIjItWTll2PSlfIzMDoAqx+Jd/nyAwA+Cg+smNIft8QEo8ZkxfRVB7D7VInYsegGsAAREVGT/Ha+HBM/SUOFwYwBUQFYMaU/fBR2XUzsUJSeMix9uC9u76GB0WLF9M8OID2nTOxYdJ1YgIiI6JpOF1dg0op0VNSYMSAyAJ9O6Q9vufuUnzpyDyneHd8bI7oGocZkxd9X7kdWfrnYseg6sAAREdFVXdBdxqRP0lFebUJcuD9WTHWvmZ8/k3tI8eHDfZEQ3RaVBjMmr0jHKS6RdzosQERE1ChdtQmTV6SjUFeD6CAfrJjSH75uXH7qKD1l+HhyP/SJ8IfusglTPt2PYn2N2LGoGViAiIioQTUmCx75bD9OFlUiRKXAZ38fgAA3OOG5qXwUHvhkcn/bPkF/X7UfVQaz2LGoiViAiIjoLwRBwD/XHcL+3EvwU3pg1d8HoH0bb7FjOZw2PnKsnNofbX3kOFKgx+w1B2G2WMWORU3AAkRERH/x/rbT2PTbBXhIJVg2sR9iNNwPrTEd2vrg48n9oPSUYvuJErz4w1GxI1ETsAAREVE9yUcu4M0tJwEAr4zuiYSObUVO5Ph6R7TB/43rDYkE+GJfHlannRM7El0DCxAREdkcLdThqa8OAQCmDo7EuAERIidyHkk9NPjnbV0BAC9uPIoDudwjyJGxABEREQCgrMqIGZ9l4LLJgqGdA/HCnd3EjuR0HhveEXfGamCyCPjHFweh1XFlmKNiASIiIli
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkAAAAGdCAYAAAD60sxaAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABaJElEQVR4nO3de1yT9eIH8M822AbIhlwHCiLe0LyLIl7SksKyi2Wp5fH2M7WOdio9lZZpp1PHsuupLMtKPeWl48nKzCzCuyAoSt4AL4AgOC4iGxfZ2Pb8/gBmFCgo49nl83699vKc7XnGZ0vdx++e7/crEQRBABEREZELkYodgIiIiKitsQARERGRy2EBIiIiIpfDAkREREQuhwWIiIiIXA4LEBEREbkcFiAiIiJyOSxARERE5HLcxA5grywWCwoKCuDt7Q2JRCJ2HCIiImoGQRBQXl6OkJAQSKVNj/OwADWhoKAAoaGhYscgIiKiG5CXl4eOHTs2+TgLUBO8vb0B1L6BKpVK5DRERETUHHq9HqGhodbP8aawADWh/msvlUrFAkRERORgrnf5Ci+CJiIiIpfDAkREREQuhwWIiIiIXA4LEBEREbkcFiAiIiJyOSxARERE5HJYgIiIiMjlsAARERGRy2EBIiIiIpfDAkREREQuhwWIiIiIXA4LEBEREbkcboZK1AST2YJzxZXIK62CVl+NkgoDDCYLzBYBEgBqT3f4eMgRrFaia2A7dPDxgFR67c33iIjIPrAAEdWprjEjObsUezKLkZp7GRkX9TCYLM0+38Ndhr4d1Rga4YeYLn6I6tQebjIOshIR2SOJIAiC2CHskV6vh1qthk6ng0qlEjsO2YjJbMG+MyX435ELSEgvRHVNw8LTTuGGcH9PaFRKBHgroXSXwl0mhdkiQHelBpcrjci7XIXskkrUmBv+UfLzkmNsbw3u6xeCIZ19IZFwdIiIyNaa+/nNAtQEFiDnpquqwVfJ57EuMQdF5Qbr/RqVEqO6B2B4N3/07aBGmK9ns77WMpktyLlUiZTsyziYdQn7zhTjclWN9fFIjTemxYTjgQEd4CGX2eQ1ERFR8z+/22R8fuXKlQgPD4dSqUR0dDRSUlKuefzmzZsRGRkJpVKJPn36YPv27Q0eFwQBS5cuRXBwMDw8PBAbG4szZ840OOa+++5DWFgYlEolgoODMXXqVBQUFLT6ayPHcrnSiNd+PIWY1xPw5s+ZKCo3wNdLjhnDwvHD/BFIWnw73nioL+7rF4Jwf69mX9PjJpOia6A3Ho0Ow/uPDEDKi7H4ctYQTIzqCA93GTK05Xjh2+MYuWInPtuXheoas41fKRERXYvNR4C+/vprTJs2DatWrUJ0dDTee+89bN68GZmZmQgMDPzT8YmJibj11luxfPly3HPPPdiwYQPeeOMNHDlyBL179wYAvPHGG1i+fDnWrVuHzp0746WXXsLx48dx6tQpKJVKAMC7776LmJgYBAcHIz8/H3//+9+tz98cHAFyLtU1Zqw5kIOPdp9FebUJQO2ozNxRERjXJwRyN9v9W0B3pQabD+dhXVIO8kqvAAACvRV4Nq4HJgzsyAuniYhakd18BRYdHY3Bgwfjww8/BABYLBaEhobiySefxKJFi/50/KRJk1BZWYlt27ZZ7xs6dCj69++PVatWQRAEhISEYOHChdZSo9PpEBQUhLVr12Ly5MmN5ti6dSvGjx8Pg8EAd3f36+ZmAXIe+84U48VvTyC3tAoA0DNYhefG9sDo7gFtel1OjdmCb1Iv4IOdZ5FfVluEBob54JX7e6N3B3Wb5SAicmZ28RWY0WhEamoqYmNjr/5AqRSxsbFISkpq9JykpKQGxwNAXFyc9fjs7GxotdoGx6jVakRHRzf5nKWlpVi/fj2GDRvWrPJDzqGsyogFX6dh6ucpyC2tgkalxDsT++HHJ0fgth6BbX5RsrtMislDwrDz76Ow+K5IeMplOJJbhvs+3I+3fs6EsQUzzoiI6ObYtACVlJTAbDYjKCiowf1BQUHQarWNnqPVaq95fP2vzXnO559/Hl5eXvDz80Nubi6+//77JrMaDAbo9foGN3JcB7Mu4a5/78OWo/mQSIAZw8Lx68JReNAOvnJSuMkwd1QX7Fw4GuP6BsMiAB/uOov7Vx5Ahpa/74iI2oJTL1Ly7LPP4ujRo/jll18gk8kwbdo0NPWN3/Lly6FWq6230NDQNk5LraHGbMHbv2TikdUHcVFXjc7+XtjyxDC8fN8taKewr2WvNGolVj46EB9NGYj2nu5Iv6jH/R8ewObDeWJHIyJyejYtQP7+/pDJZCgsLGxwf2FhITQaTaPnaDSaax5f/2tzntPf3x/du3fHHXfcgU2bNmH79u04ePBgoz938eLF0Ol01lteHj+EHE1ppRFTP0/GBzvPQhCAiVEdse3JERgQ1l7saNd0d59g/PLMKIzuEQCDyYJn/3cMi7cc40wxIiIbsmkBksvlGDRoEBISEqz3WSwWJCQkICYmptFzYmJiGhwPAPHx8dbjO3fuDI1G0+AYvV6P5OTkJp+z/ucCtV91NUahUEClUjW4kePI1Jbj/pX7cTCrFF5yGT54ZABWPNQPXnY26tOUAG8Fvpg+GAvv6A6JBNiYkoeJnyShSF8tdjQiIqdk80+HBQsWYPr06YiKisKQIUPw3nvvobKyEjNnzgQATJs2DR06dMDy5csBAE899RRGjRqFt99+G+PGjcOmTZtw+PBhfPrppwAAiUSCp59+Gq+++iq6detmnQYfEhKC8ePHAwCSk5Nx6NAhjBgxAu3bt8e5c+fw0ksvoUuXLtcsSeSYEtIL8beNR1FpNCPM1xOfTY9C9yBvsWO1mFQqwZNjuqF/mA/+tvEojl3Q4YGPErFm5mCHfD1ERPbM5gVo0qRJKC4uxtKlS6HVatG/f3/s2LHDehFzbm4upNKrA1HDhg3Dhg0bsGTJErzwwgvo1q0bvvvuO+saQADw3HPPobKyEnPmzEFZWRlGjBiBHTt2WNcA8vT0xJYtW7Bs2TJUVlYiODgYY8eOxZIlS6BQKGz9kqkNbT6ch0VbjsNsERAT4Vd7PY2XXOxYN2VktwB8N284Zq45hKySSkz4OBGf/GUQhnX1FzsaEZHT4FYYTeA6QPbvkz3nsPynDADAQ4M6YvmDfeDuRJuPXq40Ys6Xh3Eo5zLcZRK8P3kA7uoTLHYsIiK7ZhfrABHZgiAIeP2nDGv5mXNrBN58qK9TlR8AaO8lx5ezojGuTzBqzALmbzyK79PyxY5FROQUHOMKUaI69eXnk71ZAIDFd0Vi7qguIqeyHaW7DO8/MgAechn+l3oBT3+dBkONBRMHc5kGIqKbwQJEDkMQBKz4OdNafv45vjemDu0kcirbk0klWDGhLxRuUqxPzsVz3xwDAJYgIqKb4FzfGZBTeyf+ND7efQ4A8Mr9t7hE+aknlUrw6vjemDk8HACwaMsx/HjsorihiIgcGAsQOYTP9mXhg51nAQDL7u2FaTHh4gYSgUQiwdJ7euGRIWGwCMDTXx/FrswisWMRETkkFiCye9+n5ePVH9MBAM+PjcTM4Z1FTiQeiaR2JOjefiGoMQt4/MtUpGSXih2LiMjhsACRXdt3phh/3/wbAGDm8HA8PipC5ETik0kleGdiP8T2DITBZMHs/xzGueIKsWMRETkUFiCyWyfydXj8y1TUmAXc0zcYL43rBYlE3J3c7YW7TIoPHx2IAWE+0F2pwcw1h3CpovFtXoiI6M9YgMguFeqr8X9rD6HSaMawLn54e2I/SKUsP7+ndJdh9bQohPp6ILe0Co/95zA3UCUiaiYWILI71TVmzPnPYRSVG9AtsB1WTR0EhZtM7Fh2yb+dAmtmDIHawx1Hc8uw4L9psFi4uDsR0fWwAJFdEQQBz39zDL9d0MHH0x2fTx8MldJd7Fh2rWtgO3wydRDcZRJsP67Fh7vOih2JiMjusQCRXfl4zzl8n1YAN6kEH00ZiDA/T7EjOYShEX547YE+AIB3fz2NnRmFIiciIrJvLEBkN3ZlFuHNnzMBAMvuuwXDunD385aYGBWKvwwNgyAAT21KQ3ZJpdi
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAi8AAAGsCAYAAAAVGEevAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABXVElEQVR4nO3deVxU5eIG8OfMMAz7oOwIyKKIua8I7mlqWWZRlllqkbZo93fVW0l1s26L1W2717Q9t0yzUlMzc99xQ3EHBWURZBMZ9mGW8/tjkC6JCMjMmeX5fj7z+dyGc5yHc5V5eOe87yuIoiiCiIiIyErIpA5ARERE1BwsL0RERGRVWF6IiIjIqrC8EBERkVVheSEiIiKrwvJCREREVoXlhYiIiKwKywsRERFZFZYXIiIisiosL0RERGRVbLq87NmzB/fddx8CAwMhCALWrVtn0td74403IAhCvUdUVJRJX5OIiMje2HR5qaioQI8ePbBw4UKzvWaXLl1w5cqVuse+ffvM9tpERET2wEHqAKZ099134+67777p1zUaDV599VWsXLkSJSUl6Nq1K95//30MGzasxa/p4OAAf3//Fp9PREREjbPpkZdbmTlzJhITE7Fq1SqcPHkSDz/8MMaMGYMLFy60+M+8cOECAgMDER4ejkmTJiErK6sVExMREZEgiqIodQhzEAQBa9euxfjx4wEAWVlZCA8PR1ZWFgIDA+uOGzlyJPr3749333232a/x+++/o7y8HJ06dcKVK1fw5ptvIicnB6dPn4a7u3trfStERER2zaY/NmrMqVOnoNfrERkZWe95jUYDLy8vAEBKSgo6d+7c6J/z8ssv47333gOAeh9Rde/eHdHR0Wjfvj1Wr16N+Pj4Vv4OiIiI7JPdlpfy8nLI5XIkJSVBLpfX+5qbmxsAIDw8HOfOnWv0z7ledBri6emJyMhIpKWl3X5gIiIiAmDH5aVXr17Q6/UoKCjA4MGDGzzG0dHxtqY6l5eXIz09HU888USL/wwiIiKqz6bLS3l5eb1Rj0uXLiE5ORlt27ZFZGQkJk2ahMmTJ+Ojjz5Cr169UFhYiO3bt6N79+4YO3Zss1/vH//4B+677z60b98eubm5mDdvHuRyOSZOnNia3xYREZFds+kbdnft2oXhw4ff8PyUKVOwZMkSaLVavP3221i2bBlycnLg7e2NAQMG4M0330S3bt2a/XqPPvoo9uzZg6tXr8LHxweDBg3CO++8g4iIiNb4doiIiAg2Xl6IiIjI9tj1Oi9ERERkfVheiIiIyKrY3A27BoMBubm5cHd3hyAIUschIiKiJhBFEWVlZQgMDIRM1vjYis2Vl9zcXAQHB0sdg4iIiFogOzsbQUFBjR5jc+Xl+jL82dnZ8PDwkDgNERERNUVpaSmCg4ObtJ2OzZWX6x8VeXh4sLwQERFZmabc8sEbdomIiMiqsLwQERGRVWF5ISIiIqvC8kJERERWheWFiIiIrArLCxEREVkVlhciIiKyKiwvREREZFVYXoiIiMiqmKW8LFy4EKGhoXByckJ0dDQOHz7c6PE//fQToqKi4OTkhG7dumHTpk3miElERERWwOTl5ccff8Ts2bMxb948HDt2DD169MDo0aNRUFDQ4PEHDhzAxIkTER8fj+PHj2P8+PEYP348Tp8+beqoREREZAUEURRFU75AdHQ0+vXrh88++wwAYDAYEBwcjBdeeAFz58694fhHHnkEFRUV2LhxY91zAwYMQM+ePfHFF1/c8vVKS0uhUqmgVqu5txEREZGVaM77t0k3ZqypqUFSUhISEhLqnpPJZBg5ciQSExMbPCcxMRGzZ8+u99zo0aOxbt26Bo/XaDTQaDR1/11aWnr7wRugrtTigz9S4O6kgLuTAzycFQhq44xwb1e083SGg5y3DxERkW2p0OiQml+G3JIqFJRqcLVCg6oaA9q6KjDzzo6S5TJpeSkqKoJer4efn1+95/38/JCSktLgOXl5eQ0en5eX1+Dx8+fPx5tvvtk6gRtRWK7BikNZDX7N0UGGnkGe6BvaBrER3ogObwsFywwREVmZgtJq7L1QhL0XCnHishoZVyvQ0Ocz4T6utltezCEhIaHeSE1paSmCg4Nb/XU8nB3wfyM6oqxah7JqLUqqtMgursSlogpodAYczijG4YxiLNqVDpWzAiM7++GBXu0QG+EFmezW23sTERFJobiiBr+dzMWa4zk4nlVyw9d93ZVo7+UCX3cneLs5wkXpAB83pfmD/g+Tlhdvb2/I5XLk5+fXez4/Px/+/v4NnuPv79+s45VKJZRK019EX3cnzLor8obnDQYRl65WICnjGg5dKsbu8wUoKq/BL8cu45djlxHm7YpJ0SF4pF8w3J0UJs9JRETUFGdy1fh27yVsOJkLrd44vCIIQLd2Kgzu6I3oMC90CfSAl8RFpSFmuWG3f//+WLBgAQDjDbshISGYOXPmTW/YraysxIYNG+qei42NRffu3a3ihl29QcTRjGJsOJmLdcdzUa7RAQA8XRR4elAYpsSGssQQEZFkkrNL8NGWVOy9UFT3XJdADzzQqx3G9QiEr4eTJLma8/5t8vLy448/YsqUKfjyyy/Rv39/fPrpp1i9ejVSUlLg5+eHyZMno127dpg/fz4A41TpoUOH4r333sPYsWOxatUqvPvuuzh27Bi6du16y9eTurz8rwqNDr8m5+KbfRdxsbACgLHEzLkrEhP7h/AmXyIiMpv0wnJ8sDkFf5wxfrohlwkY2y0A8YPC0CPYU9pwsKDZRoBxJKWwsBCvv/468vLy0LNnT2zevLnuptysrCzIZH++icfGxuKHH37Aa6+9hldeeQUdO3bEunXrmlRcLI2r0gGP1X5ktPFkLv67/QLSCyvwz1/PYMWhLPzr/q7oH9ZW6phERGTDqrV6fLYjDV/uSYdWL0ImAHG9g/C3ER0R3NZF6ngtYvKRF3OzpJGXv9LpDfjhcBY+2nIe6iotAGBqbCheGtMJLo5Wf+80ERFZmAPpRZj7yylkFVcCAIZ18sEr93RGpJ+7xMluZFEfG5mbJZeX665V1OD9zSlYdSQbANDeywWfPNITvUPaSJyMiIhsgUanx8dbzuOrvRchioC/hxPeGHcHRnfxhyBY5gxYlhcLLy/X7T5fiLm/nMQVdTUcZAIS7umMpwaGWuxfLCIisnzpheX428rjOJNrXLT10X7BeO3eO+CmtOwRfpYXKykvAFBarUXCmlP47eQVAMCYLv744OHu8OCMJCIiaqatZ/Mx68dklGt0aOOiwHtx3TG6S8NLjVia5rx/c7qLxDycFPhsYi+8Oa4LFHIBm8/k4YGF+5F1tVLqaEREZCUMBhGfbD2PacuOolyjQ/+wtvjj70Osprg0F8uLBRAEAVNiQ/Hzs7EIUDkhvbACDyzaj6TMa1JHIyIiC1et1eP5Fcfwn+0XABgngqx4Olqy9VrMgeXFgvQI9sS6GQPRJdADVytqMPHrg3UfJxEREf1VSWUNHv/mEDafyYOjXIYPHuqON8Z1sfn99Wz7u7NCfh5OWP1MDEZ29kWNzoCZK4/hxyMNbwhJRET2K7ekCg9/kYijmdfg7uSAZfH9MaFv6+/tZ4lYXiyQq9IBXz7RFxP7h0AUgZd/OYUl+y9JHYuIiCxERlEF4j4/gAsF5fD3cMJPz8ZgQLiX1LHMhuXFQsllAt59oCueHhQGAHhjw1ks2pUmcSoiIpLapaIKPPJVIq6oqxHh44pfno9FlL/lz65tTSwvFkwQBLw6tjP+NqIjAOCDzan4es9FiVMREZFU0gvL8ciXicgv1SDSzw2rpsegnaez1LHMjuXFwgmCgNl3ReIfoyIBAO9sOofvD2ZKnIqIiMwt82oFHv3qIArKNOjk544fpg2Aj7tS6liSYHmxEjPv7Ijnh0UAAP7562msOXZZ4kRERGQuBaXVePzbQygs0yDK3x0/TIuGt5t9FheA5cWqvDi6E6bGhkIUgX/8dAI7UvKljkRERCamrtRi8neHkV1chfZeLlgW3x9edlxcAJYXqyIIAl6/9w481CcIBhGY+cNxnLqsljoWERGZSFWNHvFLjyA
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# Constants\n",
    "GM = 1.3271e20   # Gravitational parameter for the Sun (SI units)\n",
    "rAU = 1.496e11    # AU = astronomical unit (Earth-Sun distance) in m\n",
    "secperyear = 3600*24*365\n",
    "\n",
    "# Compute the velocity and energy for a circular orbit\n",
    "vEarth = np.sqrt(GM / rAU)\n",
    "E0 = - 0.5 * GM / rAU\n",
    "\n",
    "# Initial conditions\n",
    "y = 0\n",
    "x = rAU\n",
    "vx = 0\n",
    "vy = vEarth\n",
    "\n",
    "# Timestep \n",
    "dt = 1e-3\n",
    "nsteps = int(1.0 / dt)\n",
    "dt *= secperyear\n",
    "t=0.0\n",
    "\n",
    "print('nsteps = %d' % (nsteps,))\n",
    "\n",
    "t_vec = np.zeros(nsteps)\n",
    "x_vec = np.zeros(nsteps)\n",
    "y_vec = np.zeros(nsteps)\n",
    "vx_vec = np.zeros(nsteps)\n",
    "vy_vec = np.zeros(nsteps)\n",
    "\n",
    "for i in range(nsteps):\n",
    "    \n",
    "    # compute the acceleration\n",
    "    rr = np.sqrt(x**2 + y**2)\n",
    "    ax = - GM * x / rr**3\n",
    "    ay = - GM * y / rr**3\n",
    "    \n",
    "    # take the timestep\n",
    "    # for semi-implicit Euler (second order), do the velocity update first\n",
    "    # for explicit Euler (first order), do the position updates before the velocity\n",
    "    vx = vx + ax * dt\n",
    "    vy = vy + ay * dt\n",
    "    x = x + vx * dt\n",
    "    y = y + vy * dt\n",
    "\n",
    "    t = t + dt\n",
    "\n",
    "    # store the results\n",
    "    t_vec[i] = t\n",
    "    x_vec[i] = x\n",
    "    y_vec[i] = y\n",
    "    vx_vec[i] = vx\n",
    "    vy_vec[i] = vy\n",
    "\n",
    "    \n",
    "rr = np.sqrt(x_vec**2 + y_vec**2)\n",
    "energy = - GM / rr  + 0.5 * (vx_vec**2 + vy_vec**2)\n",
    "print('Delta E/E = ', np.mean((energy-E0)/E0))\n",
    "\n",
    "plt.plot(x_vec/rAU, y_vec/rAU)\n",
    "plt.show()\n",
    "\n",
    "plt.clf()\n",
    "plt.plot(t_vec/secperyear, (np.sqrt(x_vec**2 + y_vec**2)-rAU)/rAU)\n",
    "plt.show()\n",
    "\n",
    "plt.clf()\n",
    "plt.plot(t_vec/secperyear, (np.sqrt(vx_vec**2 + vy_vec**2)-vEarth)/vEarth)\n",
    "plt.show()\n",
    "\n",
    "plt.clf()\n",
    "plt.plot(t_vec/secperyear, (energy-E0)/E0)\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f2bcee48",
   "metadata": {},
   "source": [
    "Now do this for different timesteps"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "fb481046",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[-1.  -1.5 -2.  -2.5 -3.  -3.5 -4.  -4.5 -5.  -5.5 -6.  -6.5 -7.  -7.5\n",
      " -8.  -8.5 -9. ]\n",
      "-1.0 10 -0.07183938624734791 5.125999450683594e-05\n",
      "-1.5 31 -0.0005620920258684807 5.2928924560546875e-05\n",
      "-2.0 100 -2.5277116776878066e-07 0.0001468658447265625\n",
      "-2.5 316 -1.0393759622775783e-08 0.0004467964172363281\n",
      "-3.0 1000 -2.1997291674745818e-10 0.0010802745819091797\n",
      "-3.5 3162 -2.6377522768445045e-11 0.002825021743774414\n",
      "-4.0 10000 -2.093522445015285e-12 0.008417129516601562\n",
      "-4.5 31622 -2.2576017123215092e-13 0.021621227264404297\n",
      "-5.0 100000 -3.8836124694102156e-14 0.057350873947143555\n",
      "-5.5 316227 -6.745928926103558e-14 0.1693429946899414\n",
      "-6.0 1000000 1.5937593040555416e-13 0.5223250389099121\n",
      "-6.5 3162277 2.394670387712458e-13 1.6560380458831787\n",
      "-7.0 10000000 -9.946885639645126e-13 5.324723243713379\n",
      "-7.5 31622776 -1.84908330723476e-13 16.520429134368896\n",
      "-8.0 100000000 6.651862188090161e-13 51.796481132507324\n",
      "-8.5 316227766 1.1378043867991797e-12 162.00765204429626\n",
      "-9.0 1000000000 -1.0711513838639733e-12 515.1064240932465\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkkAAAG0CAYAAAAmZLNuAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABPJklEQVR4nO3deVhU9eI/8PcM+74KyCauuIOCoKYmihqVuWTZcg3XW92bVtwW7ZaWN29+s8xMi37lWma2XJeb5VVRck1kMzdQFBRBRhAYYJAZZub8/kCnkFFZZjizvF/Pw+MzZ4Yz7wPmvDufzzkfiSAIAoiIiIioEanYAYiIiIhMEUsSERERkR4sSURERER6sCQRERER6cGSRERERKQHSxIRERGRHixJRERERHrYih3AnGm1WhQXF8PNzQ0SiUTsOERERNQMgiCguroagYGBkErvfL6IJakNiouLERISInYMIiIiaoXCwkIEBwff8XmWpDZwc3MD0PBDdnd3FzkNERERNUdVVRVCQkJ0n+N3wpLUBreG2Nzd3VmSiIiIzMy9pspw4jYRERGRHixJRERERHqwJBERERHpwZJEREREpAdLEhEREZEeVl+SJk2aBC8vL0yZMkXsKERERGRCrL4kvfjii9i4caPYMYiIiMjEWH1JGjly5D1vJkVERETWx6RL0oEDBzB+/HgEBgZCIpFg27ZtTV6zevVqhIWFwdHREbGxsUhLS2v/oERERGRxTPqO2wqFAhEREZg5cyYmT57c5PktW7YgKSkJycnJiI2NxYoVKzBu3Djk5ubCz88PABAZGQm1Wt3ke3fv3o3AwMAW5VEqlVAqlbrHVVVVLTwiIiIiMhcmXZISEhKQkJBwx+eXL1+OOXPmYMaMGQCA5ORk7Ny5E2vXrsX8+fMBANnZ2QbL89577+Gdd94x2P6IiIjIdJn0cNvdqFQqZGRkID4+XrdNKpUiPj4eR48eNcp7LliwAHK5XPdVWFholPchIiIi8Zn0maS7KSsrg0ajgb+/f6Pt/v7+yMnJafZ+4uPjceLECSgUCgQHB+P777/HkCFD9L7WwcEBDg4ObcpNRERE93ZeVg1HOxuEeDuLlsFsS5Kh7N27V+wIREREdJtFO04jLb8cyx7rj0kDgkXJYLbDbb6+vrCxsYFMJmu0XSaTISAgQKRURERE1FYZl8px5MJ1AMCgMG/RcphtSbK3t0dUVBRSUlJ027RaLVJSUu44XEZERESmb2VKHgDg0YHBCPbicJteNTU1yMvL0z3Oz89HdnY2vL29ERoaiqSkJCQmJiI6OhoxMTFYsWIFFAqF7mo3IiIiMi8nCivx67lS2Egl+FtcV1GzmHRJSk9PR1xcnO5xUlISACAxMRHr16/H1KlTUVpaioULF6KkpASRkZHYtWtXk8ncREREZB4+2ddwcmRCRCA6+biImkUiCIIgagIzVlVVBQ8PD8jlcri7u4sdh4iIyKydKa7CgysPQiIB9rx8P7r5uRrlfZr7+W22c5KIiIjIsqzafx4A8FC/jkYrSC3BkkRERESiOy+rxi+nSgAAL4zqJnKaBixJREREJLpV+/MgCMC4Pv7oGWAaU1hYkoiIiEhU+WUK/PdEMQBg7qjuIqf5A0sSERERiWr1/jxoBWBUTz/0DfIQO44OSxIRERGJprC8FluzigAAc01kLtItLElEREQkmk9TL0CjFTC8uy8GhHqJHacRliQiIiISRXHlDfyQUQjAtOYi3cKSRERERKL4/NcLqNcIiO3sjZjO4i1keycsSURERNTurlXVYfPxhrNI80ab3lkkgCWJiIiIRPD/DlyESq3FwFBPDO3qI3YcvViSiIiIqF1dr1Fi07HLAIC5o7tDIpGInEg/liQiIiJqV18eyseNeg36B3tgZI8OYse5I5YkIiIiajeVtSpsPFIAoOGKNlM9iwSwJBEREVE7Wnu4AAqVBr06uiO+l5/Yce6KJYmIiIjaRVVdPdYdzgfQcHdtUz6LBLAkERERUTvZcLgA1XVqdPdzxQN9AsSOc08sSURERGR0NUo11tw8i/TCqG6QSk37LBLAkkRERETt4OvfLqGyth6dfV3wcP9AseM0C0sSERERGdUNlQZfHrwIAPjbyK6wMYOzSABLEhERERnZN2mXUVajQrCXEyYOCBI7TrOxJBEREZHR1NVr8PmvFwAAfxvZDXY25lM9zCcpERERmZ3v0wtxrVqJjh6OeDTKfM4iASxJREREZCQqtRafpTacRXru/q5wsLUROVHLsCQRERGRUfwn8wqK5XXo4OaAqYNCxI7TYixJREREZHBqjRaf3jyL9OyILnC0M6+zSABLEhERERnB9uxiXC6vhbeLPZ6KDRU7TquwJBEREZFBabQCVu/PAwDMHt4Zzva2IidqHZYkIiIiMqidJ6/iYpkCHk52eGZImNhxWo0liYiIiAxGqxWwat95AMDM+zrD1cE8zyIBLElERERkQLvPlOCcrAZuDraYfl+Y2HHahCWJiIiIDEIQBHyyr2EuUuLQMHg42YmcqG1YkoiIiMgg9uVcw+niKjjb22DmsM5ix2kzliQiIiJqM0EQsPLmWaRpgzvB28Ve5ERtx5JEREREbXbwfBlOFFbC0U6K2cO7iB3HIFiSiIiIqE0a5iI1XNH2ZEwoOrg5iJzIMFiSiIiIqE1+u1iO4wUVsLeR4tkRXcWOYzAsSURERNQmt84iPT4oGAEejiKnMRyWJCIiImq1jEvlOHLhOmylEjx3v+WcRQJYkoiIiKgNVqY0XNE2JSoYwV7OIqcxLKsuSZWVlYiOjkZkZCT69u2LL774QuxIREREZuNEYSV+PVcKG6kEfxvZTew4Bme+C6oYgJubGw4cOABnZ2coFAr07dsXkydPho+Pj9jRiIiITN6tuUgTIgMR6mNZZ5EAKz+TZGNjA2fnhl+qUqmEIAgQBEHkVERERKbvdLEce89eg0QC/D3O8s4iASZekg4cOIDx48cjMDAQEokE27Zta/Ka1atXIywsDI6OjoiNjUVaWlqL3qOyshIREREIDg7Gq6++Cl9fXwOlJyIislyrbt5d++H+gejawVXkNMZh0iVJoVAgIiICq1ev1vv8li1bkJSUhEWLFiEzMxMREREYN24crl27pnvNrflGt38VFxcDADw9PXHixAnk5+fjm2++gUwma5djIyIiMlfnZNX45VQJAOAFCz2LBJj4nKSEhAQkJCTc8fnly5djzpw5mDFjBgAgOTkZO3fuxNq1azF//nwAQHZ2drPey9/fHxERETh48CCmTJmi9zVKpRJKpVL3uKqqqplHQkREZDlunUV6oE8AwgPcRE5jPCZ9JuluVCoVMjIyEB8fr9smlUoRHx+Po0ePNmsfMpkM1dXVAAC5XI4DBw4gPDz8jq9/77334OHhofsKCQlp20EQERGZmYulNfjp94bRmBdGWe5ZJMCMS1JZWRk0Gg38/f0bbff390dJSUmz9nHp0iUMHz4cERERGD58OObOnYt+/frd8fULFiyAXC7XfRUWFrbpGIiIiMzN6v0XoBWA0T390DfIQ+w4RmXSw23GFhMT0+zhOABwcHCAg4NlLNpHRETUUpev12JbdhEAYO7o7iKnMT6zPZPk6+sLGxubJhOtZTIZAgICREpFRERkuT77NQ8arYDh3X0RGeIpdhyjM9uSZG9vj6ioKKSkpOi2abVapKSkYMiQISImIyIisjxFlTfwQ8YVAMA8KziLBJj4cFtNTQ3y8vJ0j/Pz85GdnQ1vb2+EhoYiKSkJiYmJiI6ORkxMDFasWAGFQqG72o2IiIgM4/NfL6BeI2BwF28MCvMWO067MOmSlJ6ejri4ON3jpKQkAEBiYiLWr1+PqVOnorS0FAsXLkRJSQkiIyOxa9euJpO5iYiIqPUUSjW2HG+4WGneKOs4iwQAEoHrcLRaVVUVPDw8IJfL4e7uLnYcIiIio/jl5FU8vykTnXyckfrKSEgkErEjtUlzP7/Ndk4SERERtY89ZxoukhrTy9/sC1JLsCQRERHRHak1WuzLbVjua0xv65rOwpJEREREd5R+qQKVtfXwdLZDVCcvseO0K5Y
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkkAAAG1CAYAAADtOGDLAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABrSUlEQVR4nO3dd3zT1f7H8VfSvUtbaCktlL1pmWUIAlaWggsVXIBeVNziuOK9PxCviooDFbxOLrhxgbgQRBFkUygKMgoUKKMtBbp3kt8fgZRKmR3ftH0/H488bE6++eYTIuTdc873HJPNZrMhIiIiImWYjS5ARERExBkpJImIiIiUQyFJREREpBwKSSIiIiLlUEgSERERKYdCkoiIiEg5FJJEREREyqGQJCIiIlIOV6MLqMmsViuHDh3Cz88Pk8lkdDkiIiJyHmw2G9nZ2YSHh2M2n7m/SCGpAg4dOkRkZKTRZYiIiMhFSE5OJiIi4oyPKyRVgJ+fH2D/Q/b39ze4GhERETkfWVlZREZGOr7Hz0QhqQJODrH5+/srJImIiNQw55oqo4nbIiIiIuVQSBIREREph0KSiIiISDkUkkRERETKoZAkIiIiUo46H5KuueYa6tWrx8iRI40uRURERJxInQ9JDz74IB988IHRZYiIiIiTqfMhqX///udcTEpERETqHqcOScuXL2f48OGEh4djMplYsGDBacfMmjWLqKgoPD09iY2NZd26ddVfqIiIiNQ6Tr3idm5uLtHR0dx+++1ce+21pz0+b948Jk6cyFtvvUVsbCwzZsxg8ODB7NixgwYNGgAQExNDSUnJac9dvHgx4eHhF1RPYWEhhYWFjvtZWVkX+I5ERESkpnDqkDR06FCGDh16xsdfeeUVxo8fz7hx4wB46623+P7775k9ezZPPPEEAAkJCZVWz7Rp05g6dWqlnU9EREScl1MPt51NUVER8fHxxMXFOdrMZjNxcXGsXr26Sl5z0qRJZGZmOm7JyclV8joAxRZrlZ1bREREzq3GhqT09HQsFguhoaFl2kNDQ0lJSTnv88TFxXH99dfzww8/EBERcdaA5eHh4djMtio3tf1pawr9py9jZ2p2lZxfREREzs2ph9uqw88//2x0CWXYbDY+X5/MwYx8HvtiM19N6I2rS43NsiIiIjVWjf32DQkJwcXFhdTU1DLtqamphIWFGVRVxZlMJp69piN+nq5sPpDJe78nGV2SiIhInVRjQ5K7uztdu3Zl6dKljjar1crSpUvp1auXgZVVXFiAJ/93ZTsAXlmyk91HcgyuSEREpO5x6pCUk5NDQkKC4wq1pKQkEhIS2L9/PwATJ07k3XffZe7cuWzbto0JEyaQm5vruNqtJru+awT9WtWnqMTK41/+gcVqM7okERGROsWpQ9KGDRvo3LkznTt3BuyhqHPnzkyePBmAG2+8kZdeeonJkycTExNDQkICixYtOm0yd01kMpmYdm1HfD1cid93nDmr9hpdkoiISJ1istls6qK4SFlZWQQEBJCZmVllV7p9vHYf/5q/BU83M4se7EdUiE+VvI6IiEhdcb7f307dkyRwU4/G9G4eTEGxlce/+gOrht1ERESqhUKSkzOZTLxwXSe83V1Yl3SMj9buM7okERGROkEhqQaIDPLmn0PaAPD8j9tJPpZncEUiIiK1n0JSDXFrzyb0iAoir8jCpK//RFPJREREqpZCUg1hNpt4YWQnPFzN/L4rnc/WV92+cSIiIqKQVKM0DfHhscGtAXj2+20cysg3uCIREZHaSyGphhnXpyldGgeSU1iiYTcREZEqpJBUw7iYTbw4Mhp3VzO/7TzCl/EHjC5JRESkVlJIqoFaNPDl4bhWAPznu79IzSowuCIREZHaRyGphhrftynREQFkFZTwr/kadhMREalsCklOKrsom13Hd53xcVcXMy+OjMbNxcTP29JYuPlQNVYnIiJS+ykkOaln1z7LqO9H8e3ub894TOswP+4f2BKAKQu3ciS7sLrKExERqfUUkpxQoaWQzMJMSqwlNPZvfNZjJ/RvTruG/mTkFTP5my3VVKGIiEjtp5DkhDxcPHjzsjf5aNhHRNePdrQfLzh+2rFuLmamX98JV7OJH7ek8P0fh6uzVBERkVpLIclJmUwmOoR0cNxPz0/n6m+u5unVT5NfUnYRyfbhAdzTvzkAk7/ZwtEcDbuJiIhUlEJSDbHy4EqOFxxn85HNuJhcTnv8voEtaR3qx9HcIp769i8DKhQREaldXI0uQM7PVS2uIswnjGDPYNxd3AGw2WzYsGE2mXF3tQ+7XfPmKr7dfIgrOzVkcPswg6sWERGpudSTVIPENoylRb0Wjvvf7fmOcYvGcSjHfvl/p4hA7uzXDIB/L9hCRl6RIXWKiIjUBgpJNVSRpYgZG2ewMW0jPyT94Gh/8LKWNK/vw5HsQp7+TsNuIiIiF0shqYZyd3FnzpA53NTmJsa1H+do93Rz4cWR0ZhM8PXGg/yyPdXAKkVERGouhaQaLNIvkkmxk3Ax2ydyW21Wpqyags1jN3f0aQrAk19vIaug2MgyRUREaiSFpFpkfuJ8vk78mnuW3sMdl4YSFexNSlYBz363zejSREREahyFpFpkSNMhXN3iah7u+jAN/YIdw27zNiSzfOcRo8sTERGpURSSahEfNx/+0+c/jGo9CoAeTYO4rocnboFreeLrP8gpLDG4QhERkZpDIakWMplMAJRYSzjkPhvPhvNJd13I8z9q2E1EROR8KSTVYmaTmcFRg/BzC6Q4ozsfrdnPqt3pRpclIiJSIygk1WJmk5kx7cfw8/U/MbqrfaPcJ776k+XJqykoKTC4OhEREeemkFQHeLt5M2loG8IDPDmQu5v7f7mHUd+N4njBcaNLExERcVoKSXWEn6cb067rhMkln5IST/xcwwj0CDS6LBEREaelkFSHXNqqPte160fenofYv+NKCoqtABRbi0nLSzO4OhEREeeikFTH/OuKdjTwCWbfEROv/rwTgLc3v80131zD0v1LDa5ORETEeSgk1TEBXm5Mu7YjAO+t2MOGvUdYfXg1WUVZFFu0fYmIiMhJCkl10MA2oVzbuRFWGzzx9V+8fdn7vHTpSwxpOsRxjAKTiIjUdQpJddTk4e0I8fVgV1oOb/66l8FRgx2PFZQUMOr7UbyZ8CYlVq3SLSIidZNCUh0V6O3OM1d3AODt5Xv480Cm47HF+xaz8/hOvtj5BTlFOUaVKCIiYihXowsQ4wzpEMaVnRry3R+HeezLzXx7/yW4uZgZ0XwEriZXAj0DCfQMNLpMERERQ6gnqY6bOqI99bzd2J6SzY9bUhztw5oNo3d4b8f9+NR4Hv3tUTILM8s7jYiISK2jkFTHBft6MLZ3UwDmrEwq95gSawmTV07mp70/8fYfb1dneSIiIoZRSBJuim2Mm4uJjfsz2JyccdrjrmZXXrz0RS5pdAn3xtxb/QWKiIgYQCFJqO/nwfBO4QDMXbW33GPaB7fnv3H/xcfNx9H2vy3/Y9fxXdVRooiISLVTSBIAxvSOAuDbPw6Rll1wzuOXH1jOK/GvMPr70drSREREaiWFJAEgOjKQLo0DKbbY+HRt8jmPbxfcjksaXcINrW+ggXeDaqhQRESkeikkicPYPvYJ3B+t3UdRifWsx4Z4hfDmZW/yUNeHHG2ZhZksP7C8KksUERGpNgpJ4jC0Qxih/h4cyS7kxy2Hz3m8yWTCzewGgM1m45k1z3Dv0nt5a/NbVV2qiIhIlVNIEgc3FzO3xDYB4H8r917Qc602K2E+YbiaXekT3qcKqhMREaleCklSxujYxri7mElIzmDT/uPn/TwXswuPdHuEH6/9kY71OzrakzKTsFgtVVGqiIhIlVJIAvLy8mjSpAmPPvqo0aUYLsTXg+HRZ18O4GzCfMIcPx/JO8JtP97G7T/dztH8o5VVooiISLVQSAKeffZZevbsaXQZTmPsieUAvv/zMGlZ514O4EwSjydSZCkiryQPP3e/SqpORESketT5kJSYmMj27dsZOnSo0aU4jY4RAXRrUo9ii42P1+6/6PP0btSbL0d8yYv9XsTdxR2wT/D
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import time\n",
    "\n",
    "def evolve(dt):\n",
    "    \n",
    "    # Initial conditions\n",
    "    y = 0.0\n",
    "    x = rAU\n",
    "    vx = 0.0\n",
    "    vy = vEarth\n",
    "\n",
    "    nsteps = int(10**(0-dt))\n",
    "    dt = 10.0**dt * secperyear\n",
    "    \n",
    "    for i in range(nsteps):\n",
    "            \n",
    "        # compute the acceleration\n",
    "        rr = (x**2 + y**2)**(3/2)\n",
    "        ax = - GM * x / rr\n",
    "        ay = - GM * y / rr\n",
    "        \n",
    "        # take the timestep\n",
    "        # for semi-implicit Euler (second order), do the velocity update first\n",
    "        # for explicit Euler (first order), do the position updates before the velocity\n",
    "        vx = vx + ax * dt\n",
    "        vy = vy + ay * dt\n",
    "\n",
    "        x = x + vx * dt\n",
    "        y = y + vy * dt\n",
    "\n",
    "    energy = - GM / (x**2 + y**2)**0.5  + 0.5 * (vx**2 + vy**2)\n",
    "     \n",
    "    return energy, nsteps\n",
    "\n",
    "# Constants\n",
    "GM = 1.3271e20   # Gravitational parameter for the Sun (SI units)\n",
    "rAU = 1.496e11    # AU = astronomical unit (Earth-Sun distance) in m\n",
    "secperyear = 3600*24*365\n",
    "\n",
    "# Compute the velocity and energy for a circular orbit\n",
    "vEarth = np.sqrt(GM / rAU)\n",
    "E0 = - 0.5 * GM / rAU\n",
    "\n",
    "# The different values of log10 timestep to try\n",
    "dt_vals = np.arange(-1,-9.1,-0.5)\n",
    "print(dt_vals)\n",
    "\n",
    "E_vals = np.array([])\n",
    "n_vals = np.array([], dtype =np.int_)\n",
    "for dt in dt_vals:\n",
    "    t0 = time.time()\n",
    "    E, nsteps = evolve(dt)\n",
    "    t1 = time.time()\n",
    "    print(dt, nsteps, (E-E0)/E0, t1-t0)\n",
    "    E_vals = np.append(E_vals, np.abs((E-E0)/E0))\n",
    "    n_vals = np.append(n_vals, nsteps)\n",
    "\n",
    "plt.plot(dt_vals, E_vals)\n",
    "plt.yscale('log')\n",
    "plt.xlabel(r'$\\log_{10}\\Delta t\\ (\\mathrm{yr})$')\n",
    "plt.ylabel(r'$\\Delta E/E$')\n",
    "plt.show()\n",
    "\n",
    "plt.clf()\n",
    "plt.plot(n_vals, E_vals)\n",
    "plt.plot(n_vals, 1e-16 * np.sqrt(n_vals), \":\")\n",
    "plt.plot(n_vals, 1/n_vals**2, \":\")\n",
    "plt.yscale('log')\n",
    "plt.xscale('log')\n",
    "plt.xlabel(r'Number of steps')\n",
    "plt.ylabel(r'$\\Delta E/E$')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cc939434",
   "metadata": {},
   "source": [
    "2. **Interpolation and thermodynamics**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "1170d357",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgoAAAG1CAYAAACYtdxoAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAy1klEQVR4nO3deXQUdbr/8U8nkA6QdAQmgQBhEdwRF0REkIkOghwPY+aMqOjIchnvoB1HZBwVl2EZoXGd8bpEdJT4u04OjChwLyIKSEBGcGG5igsDCiaiAZRDAgGydNXvD6RnepIK6eqlKvT7dU4d7Up9q55uSHjyPN/6lsc0TVMAAACNSHE6AAAA4F4kCgAAwBKJAgAAsESiAAAALJEoAAAASyQKAADAEokCAACwRKIAAAAskSgAAABLJAoAAMASiQIAALBEooCTTnFxsTwej3bt2uV0KADQ4pEoAC5zPNE5vqWnp+v0009XYWGh9uzZ43R4AJJMK6cDANC4mTNnqlevXjp69KjWrVunoqIiLVu2TFu3blXbtm2dDg9AkiBRAFxq5MiRuuiiiyRJv/71r9WxY0c98cQTWrJkicaMGdPomOrqarVr1y6RYUatJcYMJBNaD0gKmzdv1siRI+Xz+ZSRkaGf/exn2rBhQ4PjSktLddFFFyk9PV29e/fW3LlzNX36dHk8HgeiDnfFFVdIknbu3ClJobg+++wz3XjjjWrfvr2GDBkSOn737t36j//4D3Xq1Eler1fnnHOOXnrppbBzHjx4UJMnT1bPnj3l9XqVk5OjK6+8Ups2bYromPHjx6tnz54NYv73zy4WMQNILCoKOOl9+umnuuyyy+Tz+XT33XerdevWmjt3rvLz87VmzRoNHDhQ0rFk4qqrrlJubq5mzJihYDComTNnKjs72+F3cMyXX34pSerYsWPY/tGjR+u0007T7NmzZZqmJGnPnj265JJL5PF4VFhYqOzsbL355puaOHGiqqqqNHnyZEnSpEmTtHDhQhUWFurss8/WDz/8oHXr1unzzz/XhRde2OxjIhVNzAASzAROMvPmzTMlmTt37jRN0zQLCgrMtLQ088svvwwd8+2335qZmZnm0KFDQ/tGjRpltm3b1ty9e3do3/bt281WrVqZifxWOR7/ypUrzX379pnl5eXm/PnzzY4dO5pt2rQxv/nmG9M0TXPatGmmJHPMmDENzjFx4kQzNzfX/P7778P233DDDWZWVpZ5+PBh0zRNMysry/T7/U3G05xjxo0bZ/bo0aPB/uMx/vvraGIGkFi0HnBSCwaDevvtt1VQUKBTTz01tD83N1c33nij1q1bp6qqKgWDQa1cuVIFBQXq0qVL6Lg+ffpo5MiRToSuYcOGKTs7W3l5ebrhhhuUkZGhRYsWqWvXrmHHTZo0Key1aZp67bXXNGrUKJmmqe+//z60jRgxQpWVlaG2wSmnnKL3339f3377rWUczTkmUtHEDCCxSBRwUtu3b58OHz6sM844o8HXzjrrLBmGofLycu3du1dHjhxRnz59GhzX2L6ioiJdeOGFat26taZPn97gmldffbXatWunM844Q6tWrbIV+zPPPKMVK1Zo9erV+uyzz/TVV19pxIgRDY7r1atXg+sfOHBAzz//vLKzs8O2CRMmSJL27t0rSXrkkUe0detW5eXl6eKLL9b06dP11VdfhZ2vOcdEKpqYASQWcxQAG3JzczV9+nSVlJQ0+Jrf71fnzp21b98+rVy5Utddd522b9+uDh06RHSNiy++OHTXQ1PatGkT9towDEnSr371K40bN67RMf369ZMkXXfddbrsssu0aNEivf3223r00Uf18MMP6/XXXw9VUppzjNVkz2AwGPOYASQWiQJOatnZ2Wrbtq22bdvW4GtffPGFUlJSlJeXp3bt2ik9PV07duxocFxj+woKCiRJy5YtC9t/6NAhLV68WF999ZXatm2rn//85zr33HO1ZMmS0G/G8Zadna3MzEwFg0ENGzbshMfn5ubqtttu02233aa9e/fqwgsv1KxZs8JaLic6pn379jpw4ECDc3/99ddxiRlA4tB6wEktNTVVw4cP15IlS8KWdN6zZ49KSko0ZMgQ+Xw+paamatiwYVq8eHFYL37Hjh168803m3297du3KyMjQ926dQvtO/fcc/Xpp5/G5P00R2pqqn75y1/qtdde09atWxt8fd++fZKO/bZfWVkZ9rWcnBx16dJFNTU1zT5Gknr37q3Kykp9/PHHoX3fffedFi1aFNOYASQeFQWc9B566CGtWLFCQ4YM0W233aZWrVpp7ty5qqmp0SOPPBI6bvr06Xr77bc1ePBg3XrrrQoGg3r66afVt29fbdmypVnXOnTokHw+X9g+n8+nH374IZZv6YTmzJmj1atXa+DAgbrlllt09tlna//+/dq0aZNWrlyp/fv36+DBg+rWrZuuvfZanXfeecrIyNDKlSv14Ycf6vHHH5ekZh0jSTfccIPuuece/eIXv9Bvf/tbHT58WEVFRTr99NObPQmxOTEDSDwSBZz0zjnnHL377ruaOnWqAoGADMPQwIED9corr4TWUJCk/v37680339Rdd92lBx98UHl5eZo5c6Y+//xzffHFF826VkZGhqqqqsL2VVVVKSMjI6bv6UQ6deqkDz74QDNnztTrr7+uZ599Vh07dtQ555yjhx9+WJLUtm1b3XbbbXr77bf1+uuvyzAM9enTR88++6xuvfXWZh8jHVvbYdGiRZoyZYruvvtu9erVS4FAQNu3b292otCcmAEknsc0f1ztBECjCgoK9Omnn2r79u0NvjZp0iR17tw5dOfDoUOH1KFDB+3cuTN0G+Pll1+usWPHJmyOAgDEEnMUgH9x5MiRsNfbt2/XsmXLlJ+fH7a/vr5eR48eVTAYDPv/jIwMXXPNNZo2bZqOHDmipUuX6uOPP9Y111yTwHcBALFDRQH4F7m5uRo/frxOPfVUff311yoqKlJNTY02b96s0047LXTc9OnTNWPGjLCx8+bN0/jx47Vv3z6NGzdOpaWl6tatm5599llm8gNosUgUgH8xYcIErV69WhUVFfJ6vRo0aJBmz55t+5kGANDSRdR6KCoqUr9+/eTz+eTz+TRo0KAT3jr26quv6swzz1R6errOPffcBvedA24yb9487dq1S0ePHlVlZaWWL19OkgAgqUWUKHTr1k1z5szRxo0b9dFHH+mKK67QNddcY3mP+HvvvacxY8Zo4sSJ2rx5swoKClRQUNDofdIAAMB9om49dOjQQY8++qgmTpzY4GvXX3+9qqurtXTp0tC+Sy65ROeff76ee+65aC4LAAASwPY6CsFgUK+++qqqq6s1aNCgRo9Zv369pkyZErZvxIgRWrx4cZPnrqmpCVv1zTAM7d+/Xx07drRcUx4AAOnY00gPHjyoLl26KCUlPjf3HT16VLW1tTE5V1pamtLT02NyrniIOFH45JNPNGjQIB09ejT02Nuzzz670WMrKirUqVOnsH2dOnVSRUVFk9cIBAINZpQDABCJ8vLysOXUY+Xo0aPq1SNDFXsbf+hZpDp37qydO3e6NlmIOFE444wztGXLFlVWVmrhwoUaN26c1qxZY5ks2DF16tSwSkRlZaW6d++uoedMVqtUb8yuE0uHu7RzOgRL1bmpTofQpCM57q4S1eTE5odBPKRlH3Y6hCb1+om7l10e0H6X0yE0aXC7hot8ucVgd/6bJkmqOmSox4W7lJmZGZfz19bWqmJvUDs39pAvM7qKRdVBQ736f63a2tqTJ1FIS0tTnz59JB1b8vbDDz/Uk08+qblz5zY4tnPnztqzZ0/Yvj179qhz585NXsPr9crrbZgQtEr1ujZRaNXanX/AkpSa5u5EIdXr7kQhpY17E4XUtobTITSpdbs0p0NoUnpGa6dDaFK7DPeuiedz74+8kHi3qn2ZKVEnCi1B1O/QMIyw+QT/atCgQVq1alXYvhUrVljOaQAAoKUImkZMNreLqKIwdepUjRw5Ut27d9fBgwdVUlKi0tJSvfXWW5KksWPHqmvXrgoEApKkO+64Qz/96U/1+OOP6+qrr9b8+fP
      "text/plain": [
       "<Figure size 640x480 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgcAAAG1CAYAAABtS1fYAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA+dUlEQVR4nO3dfVxUZf7/8feAcqMyo6iIKKhZq7belZl5k2HratbDohtttfWm3MrdwVL3a2Y3P7U7ym6sLVerLdlNWe0Odc2lKAVzlVpRS9qiNF3cDK1MUBQYmPP7w5j1LILMmYEZ8PV8PM6jx7nmXNf5zITwmc91nXNshmEYAgAA+ElIoAMAAADBheQAAACYkBwAAAATkgMAAGBCcgAAAExIDgAAgAnJAQAAMCE5AAAAJiQHAADAhOQAAACYkBwAAAATkgPgHLVo0SL17NlTbrfbb2MeP35cISEhWrx4seUxli1bpoSEBJWVlfktLgDeITlAk5WamiqbzVbjlpOT4/WYW7du1YIFC3T06FH/B9yAiouL9cQTT2ju3LkKCTH/Gvj00091yy23qFOnTgoLC1NsbKyGDRumhQsXVhvH7Xarffv2WrRokSQpLy9PhmGod+/elmObOnWqysvL9eKLL1oeA4BvmgU6AKC+PfTQQ+rWrVu19vPPP9/rsbZu3aqFCxdq6tSpat26tR+iC4xXX31VFRUVmjBhgqn97bff1oQJE5SQkKA777xTcXFxKigo0Pvvv68VK1Zo/vz5puM//vhjff/997rmmmskSbt375Yk9enTx3JsERERmjJlip555hnNmDFDNpvN8lgArCE5QJM3ZswYXXLJJQ1+3pKSErVs2bLBz1sXy5cv17XXXquIiAhP248//qjbbrtNAwcO1MaNGxUWFuZ57aGHHtLBgwerjbNhwwZ16dJFP//5zyWdSg7atWun2NhYn+IbP368Fi1apE2bNunKK6/0aSwA3mNaAee8BQsWyGazac+ePZ6KgMPh0K233qoTJ06YjpszZ44kqVu3bp7pif3793vG+Ne//qWJEyeqTZs2GjZsmKfvzp07NWbMGNntdrVq1Uq/+MUvqk1rVI3xxRdfaPz48bLb7Wrbtq3uvvtulZaWeo7btGmTbDab0tPTq72XtLQ02Ww2bdu2rcb3u2/fPn366acaOXKkqf39999XUVGRpk2bZkoMqsTFxVVre+eddzxVA+lUclCVKFR5+eWXFRYWppkzZ6qysrLGuE43YMAARUdHa+3atXU6HoB/UTlAk1dUVKTvv//e1Gaz2dS2bVtT2/jx49WtWzelpKRox44d+tOf/qSYmBg98cQTkqQbbrhBX375pf76179q8eLFateunSSpffv2njHGjRunCy64QI899pgMw5AkffbZZ7r88stlt9t1zz33qHnz5nrxxReVmJio7OxsDRo0qFocXbt2VUpKinJycvSHP/xBP/74o/7yl79IkhITExUfH6+VK1fq+uuvN/VduXKlunfvrsGDB9f4eWzdulWSdPHFF5vaS0pKJJ1ac1AXhYWF2rlzpx566CFP2+7duz1TFRUVFZo5c6ZeeuklLVmyRLfffnudxq1y8cUX6x//+IdXfQD4iQE0UcuXLzcknXELDw/3HDd//nxDknHbbbeZ+l9//fVG27ZtTW1PPvmkIcnYt2+fqb1qjAkTJlSLIykpyQgLCzP27t3raTt48KARFRVlDB8+vNoY1157ran/7373O0OS8cknn3ja5s2bZ4SHhxtHjx71tB0+fNho1qyZMX/+/Fo/lwceeMCQZBw7dszUvn//fqNFixaGJOOCCy4w7rnnHmPjxo1GRUXFGcd55ZVXjMjISOPEiROe9yTJWLZsmfHDDz8YV155pREdHW1s2rSp1nhqcscddxiRkZGW+gLwDdMKaPKWLFmizMxM0/b3v/+92nHTp0837V9++eX64YcfVFxcXOdz/e8YlZWVeu+995SUlKTzzjvP096xY0dNnDhRW7ZsqTa+0+k07c+YMUPSqfn9KpMnT1ZZWZnefPNNT9vq1atVUVGhX//617XG+MMPP6hZs2Zq1aqVqb1Lly7atm2bxo8fr4MHD2rRokW68sordd555+n999+vNs6GDRs0YsQIRUZGSvpvxcFms2ngwIE6ePCgPvroIyUmJtYaT03atGmjkydPmqZ2ADQMkgM0eZdeeqlGjhxp2kaMGFHtuISEBNN+mzZtJJ1aqFdX/3tVxHfffacTJ06oR48e1Y7t1auX3G63Dhw4YGq/4IILTPvdu3dXSEiI9u/f72nr2bOnBg4cqJUrV3raVq5cqcsuu8zSVRhV+vbtq9WrV+vIkSPatGmTbr31Vh04cEC/+tWvTH+kXS6XMjMzq603kKTk5GR16NBB27ZtO2Ms+/fvl81mU6tWrdSiRQu1a9dO9913X7XjjJ+mZbhaAWh4JAfAT0JDQ8/YXvVHqi6qvkX7U01/HCdPnqzs7Gz95z//0d69e5WTk3PWqoEktW3bVhUVFTp27FiNx4SFhSkxMVGvvvqqrrnmGv3www/Kz8/3vF5V8bj66qs9bbt371aXLl10+eWXa+/evTp+/PgZx/7kk0904YUX6vjx4zpx4oTee+89paSkKC8vz3Tcjz/+qBYtWtTLZwqgdiQHgBe8/Rbbvn17tWjRwvSHtcoXX3yhkJAQxcfHm9q/+uor0/6ePXvkdrvVtWtXU/uvfvUrhYaG6q9//atWrlyp5s2b6+abbz5rTD179pR06qqFuggPD5ckORwOT9s777yjCy+80BTT7t271b9/f61atUoRERG6/vrrTVdZVKlKDqr069dPzZs3l8vlMh23b98+9erVq04xAvAvkgPAC1X3LajrHRJDQ0M1atQorV271jQtcOjQIaWlpWnYsGGy2+2mPkuWLDHtP//885JO3a/hdO3atdOYMWO0YsUKrVy5UldddZXnCoraVF3JsH37dk/bli1bdPLkyWrHfvrpp8rIyNBFF11kWjOxYcMG05RCZWWlPv/8c/Xp00ft27fX22+/rby8PP32t7+tNubpyUFpaakeeeQRDRgwQP379zcdt2PHDg0ZMuSs7weA/3EpI5q8v//97/riiy+qtQ8ZMsT0B68uBgwYIEm6//779atf/UrNmzfX2LFja+3zyCOPKDMzU8OGDdPvfvc7NWvWTC+++KLKyso8tx0+3b59+3Tttdfqqquu0rZt27RixQpNnDhR/fr1q3bs5MmTddNNN0mSHn744Tq9h/POO0+9e/fW+++/r9tuu02SdO+99+rLL7/UuHHj1K9fP1VUVGjXrl167bXX5HA49Nprr5ni+/zzz7V06VJP21dffaXS0lLPnREHDBigpUuX6tZbb9WAAQOUnJzsOfaTTz5RRkaGnnvuOR0/flydOnVSdna2qSqTm5urI0eO6LrrrqvTewLgZ4G+XAKoL7VdyijJWL58uWEY/72E8Lvvvjtj//+9bPHhhx82OnXqZISEhHher2mMKjt27DBGjx5ttGrVymjRooUxYsQIY+vWraZjqsb417/+Zdx0001GVFSU0aZNGyM5Odk4efLkGcctKysz2rRpYzgcjhqPOZNnnnnGaNWqlecyxLffftuYMGGCcf755xstW7Y0IiIijF69ehlz5swxDh8+bOr7wgsvGA6Hw3C5XJ62119/3ZBkfPbZZ6Zjf/e73xnNmzc3srOzDcMwjGPHjhk2m8346quvDMMwjNLSUmPq1KnGuHHjTP3mzp1rJCQkGG63u87vCYD/2AzDi9VWAOrNggULtHDhQn333Xd1mh6QTt1oKC4uTmPHjtUrr7xS53MVFRXpvPPO06JFizRt2jSv4rz66qvVqlUrvf766171k07dgOnqq6/Wjz/+6KkUvPTSS3rttdf04YcfSpLKysrUtWtX3Xvvvbr77ru9PgcA37HmAGjE1qxZo++++06TJ0/2qp/D4dA999yjJ5980utHNicmJmrWrFle9anyySef6KKLLvIkBl999ZWWLFmipKQkzzHLly9X8+bNq90zAkDDoXIABAlvKgcfffSRPv30Uz388MNq166dduzY0UBR+mb
      "text/plain": [
       "<Figure size 640x480 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAhcAAAGzCAYAAACcvDUtAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA/oElEQVR4nO3de1RVZf7H8c8BuSqgKAgqeKESJ3VqWRrqGKYpTlPDL7OLJtowli1sMm1SSnPUjCzLGqex62AXHcvSLGtszAt28VImkzrJpMlgKKYygmICcfbvD/NMZ+QgZ58NZ0vv11p7Lc5mP/v5nu0RvnyfZz/bYRiGIQAAAIsE+DsAAADQtJBcAAAAS5FcAAAAS5FcAAAAS5FcAAAAS5FcAAAAS5FcAAAAS5FcAAAAS5FcAAAAS5FcAAAAS5Fc4Cdt0aJFcjgctW5Tp071d3gAcF5q5u8AADuYNWuWOnfu7Lave/fufooGAM5vJBeApGHDhumyyy6r17GnTp1ScHCwAgJ+eoU/p9OpqqoqhYaG+jsUADb20/vpCHhhw4YNcjgcWrp0qaZNm6b27dsrPDxc5eXlkqQtW7YoLS1NUVFRCg8P15VXXqmPP/74rPMUFxfrN7/5jdq2bauQkBBdfPHF+stf/lLvOF599VX16tVLYWFhio6O1s0336z9+/e7HZOamqru3bvrn//8pwYOHKjw8HC1b99ejz766Fnnq6ys1IwZM3TBBRcoJCRECQkJuu+++1RZWel2nMPh0IQJE7R48WJdfPHFCgkJ0erVqyVJX3zxha688kqFhYWpQ4cOeuihh5SbmyuHw6HCwkJJ0pgxY9SmTRtVV1efFcOQIUPUtWvXel8DAOcPKheApLKyMh05csRtX5s2bVxfz549W8HBwbr33ntVWVmp4OBgrVu3TsOGDVOvXr00Y8YMBQQEKDc3V1dddZU+/PBD9e7dW5J06NAhXXHFFa5f1DExMfrb3/6mzMxMlZeXa+LEiXXGNmfOHE2fPl033nijfvvb3+rw4cNasGCBBgwYoO3bt6tly5auY//zn/8oLS1N119/vW688Ua98cYbmjJlinr06KFhw4ZJOl19uO666/TRRx/p9ttvV7du3bRjxw7Nnz9f//rXv/TWW2+59b9u3Tq9/vrrmjBhgtq0aaNOnTqpuLhYAwcOlMPhUHZ2tpo3b64XXnhBISEhbm1Hjx6tl19+We+//75+9atfufaXlJRo3bp1mjFjRn3/iQCcTwzgJyw3N9eQVOtmGIaxfv16Q5LRpUsX4+TJk652TqfTuPDCC42hQ4caTqfTtf/kyZNG586djauvvtq1LzMz04iPjzeOHDni1vfNN99sREVFuZ33fxUWFhqBgYHGnDlz3Pbv2LHDaNasmdv+K6+80pBkvPzyy659lZWVRlxcnDF8+HDXvldeecUICAgwPvzwQ7dzPvPMM4Yk4+OPP3btk2QEBAQYu3btcjv2rrvuMhwOh7F9+3bXvqNHjxrR0dGGJGPfvn2GYRhGTU2N0aFDB+Omm25ya//EE08YDofD+Prrrz2+dwDnL4ZFAElPP/201qxZ47b92JgxYxQWFuZ6nZ+fr6+++kojR47U0aNHdeTIER05ckQVFRUaNGiQNm7cKKfTKcMw9Oabb+raa6+VYRiu444cOaKhQ4eqrKxMn3/+uce4li9fLqfTqRtvvNGtbVxcnC688EKtX7/e7fgWLVro1ltvdb0ODg5W79699fXXX7v2LVu2TN26dVNycrLbOa+66ipJOuucV155pX72s5+57Vu9erVSUlJ0ySWXuPZFR0dr1KhRbscFBARo1KhRevvtt3X8+HHX/sWLF6tv375nTaIF0DSQXACSevfurcGDB7ttP/a/vwS/+uorSaeTjpiYGLfthRdeUGVlpcrKynT48GEdO3ZMzz333FnH3XbbbZKkb7/91mNcX331lQzD0IUXXnhW+y+//PKsth06dJDD4XDb16pVK/3nP/9xO+euXbvOOt9FF11Uazy1JQD//ve/dcEFF5y1v7Z9GRkZ+u6777RixQpJUkFBgbZt26bRo0d7fN+AWRs3btS1116rdu3ayeFwnDXM54/+Dh06pLFjx6pdu3YKDw9XWlqa62dIU8WcC6Aefly1kE7PW5Ckxx57zO2v9x9r0aKFjh49Kkm69dZbNWbMmFqP69mzp8d+nU6nHA6H/va3vykwMLDWPn6stmMkyTAMt3P26NFDTzzxRK3HJiQkuL3+3/furZ/97Gfq1auXXn31VWVkZOjVV19VcHCwbrzxRp/OC9SmoqJCP//5z/Wb3/xG119/vd/7MwxD6enpCgoK0sqVKxUZGaknnnhCgwcP1j//+U81b968wWP0B5ILwISkpCRJUmRk5FlVjh+LiYlRRESEampq6jyurn4Mw1Dnzp1dlQVfJSUl6R//+IcGDRp0VpWjvjp27Kg9e/actb+2fdLp6sWkSZN08OBBLVmyRNdcc41atWplqm+gLsOGDXNNXq5NZWWlHnjgAf31r3/VsWPH1L17d82dO1epqakN0t9XX32lzZs3a+fOnbr44oslSQsXLlRcXJz++te/6re//a2pfu2OYRHAhF69eikpKUnz5s3TiRMnzvr+4cOHJZ2uJAwfPlxvvvmmdu7c6fE4T66//noFBgZq5syZbtUH6fRfRGcqI9648cYbVVxcrOeff/6s73333XeqqKg45zmGDh2qTZs2KT8/37WvtLRUixcvrvX4W265RQ6HQ3fffbe+/vprt3khQGOaMGGCNm3apKVLl+qLL77QiBEjGnSY4szt3T9eGyYgIEAhISH66KOPGqRPO6ByAZgQEBCgF154QcOGDdPFF1+s2267Te3bt1dxcbHWr1+vyMhIvfPOO5KkRx55ROvXr1efPn00btw4/exnP1Npaak+//xzffDBByotLfXYT1JSkh566CFlZ2ersLBQ6enpioiI0L59+7RixQrdfvvtuvfee72KffTo0Xr99dc1fvx4rV+/Xv369VNNTY12796t119/Xe+///45FxS777779Oqrr+rqq6/WXXfd5boVNTExUaWlpWdVRGJiYpSWlqZly5apZcuWuuaaa7yKGbBCUVGRcnNzVVRUpHbt2kmS7r33Xq1evVq5ubl6+OGHLe8zOTlZiYmJys7O1rPPPqvmzZtr/vz5+uabb3Tw4EHL+7MLkgvApNTUVG3atEmzZ8/Wn/70J504cUJxcXHq06eP7rjjDtdxbdu21datWzVr1iwtX75cf/7zn9W6dWtdfPHFmjt37jn7mTp1qi666CLNnz9fM2fOlHR6XsSQIUN03XXXeR13QECA3nrrLc2fP18vv/yyVqxYofDwcHXp0kV33313vYZfEhIStH79ev3ud7/Tww8/rJiYGGVlZal58+b63e9+V+sKnhkZGVq1apVuvPHGs9bDABrDjh07VFNTc9ZnvLKyUq1bt5Yk7d69W926davzPFOmTNEjjzxSrz6DgoK0fPlyZWZmKjo6WoGBgRo8eLCGDRt2VjWyKXEYTfndAWhUEydO1LPPPqsTJ06cNbl05cqVSk9P18aNG/WLX/zCTxHip8ThcGjFihVKT0+XJL322msaNWqUdu3addbns0WLFoqLi1NVVZXbrdu1ad26tWJiYs7Z3/8qKytTVVWVYmJi1KdPH1122WV6+umnTb03u6NyAcCU7777zu1OkqNHj+qVV15R//79a71r5fnnn1eXLl3Uv3//xgwTcLn00ktVU1Ojb7/91mOCGxwcrOTk5AbpPyoqStLpSZ6fffaZZs+e3SD92AHJBQBTUlJSlJqaqm7duunQoUN68cUXVV5erunTp7sdd2bi3LvvvqunnnrK9B0qQH2cOHHC7a6lffv2KT8/X9HR0brooos0atQoZWRk6PHHH9ell16qw4cPa+3aterZs6epuUB19ZeYmCjp9MJ1MTExSkxM1I4dO3T33XcrPT1dQ4YM8f0N25WfVgYFcJ7Lzs42LrzwQiMsLMwIDw83+vfvb6xZs+as4yQZLVq0MDIzM43q6mo/RIqfkjNL9v/vNmbMGMMwDKOqqsp48MEHjU6
      "text/plain": [
       "<Figure size 640x480 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# First calculate P, S, and F on a coarse grid and make some color\n",
    "# maps to show what they look like\n",
    "\n",
    "import numpy as np\n",
    "import scipy.interpolate\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "def P(rho, T):\n",
    "    n = rho / (28*mu)\n",
    "    return n * kB * T\n",
    "\n",
    "def S(rho, T):\n",
    "    n = rho / (28*mu)\n",
    "    nQ = (28*mu*kB*T/(2*np.pi*hbar**2))**1.5\n",
    "    return kB * (2.5 - np.log(n/nQ))\n",
    "\n",
    "def F(rho, T):\n",
    "    n = rho / (28*mu)\n",
    "    nQ = (28*mu*kB*T/(2*np.pi*hbar**2))**1.5\n",
    "    return kB * T * (np.log(n/nQ) - 1.0)\n",
    "\n",
    "# constants\n",
    "mu = 1.66e-27 # kg\n",
    "kB = 1.381e-23 # SI\n",
    "hbar = 6.626e-34 # SI\n",
    "\n",
    "# grid in log T and log rho\n",
    "Tp = np.linspace(2,3,10)\n",
    "rhop = np.linspace(-6,0.0,10)\n",
    "Tgrid, rhogrid = np.meshgrid(Tp, rhop, indexing='ij')\n",
    "Pgrid = np.log10(P(10**rhogrid, 10**Tgrid))\n",
    "Sgrid = S(10**rhogrid, 10**Tgrid)/kB\n",
    "Fgrid = F(10**rhogrid, 10**Tgrid)\n",
    "\n",
    "plt.imshow(Pgrid, aspect='auto', origin='lower', extent=(rhop[0],rhop[-1],Tp[0],Tp[-1]), \n",
    "           interpolation='none')\n",
    "plt.title(r'$\\log_{10}$ Pressure')\n",
    "plt.colorbar()\n",
    "plt.show()\n",
    "\n",
    "plt.clf()\n",
    "plt.imshow(Sgrid, aspect='auto', origin='lower', extent=(rhop[0],rhop[-1],Tp[0],Tp[-1]), \n",
    "           interpolation='none')\n",
    "plt.title(r'Entropy ($S/k_B$)')\n",
    "plt.colorbar()\n",
    "plt.show()\n",
    "\n",
    "plt.clf()\n",
    "plt.imshow(Fgrid, aspect='auto', origin='lower', extent=(rhop[0],rhop[-1],Tp[0],Tp[-1]), \n",
    "           interpolation='none')\n",
    "plt.title(r'Free energy')\n",
    "plt.colorbar()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "5988543f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgcAAAG1CAYAAABtS1fYAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA+cElEQVR4nO3de3wU1f3/8fcmkAuQXQgQQiABRAsoGCxG5CIGS0HkgcYqWLBclKq0GxToFxEvP8DbKt6rFNQqaYUUFA1QpNEoJEghKkHUeEFBaFIxeEESCOS68/sDs2UMCdnZTXYJr2cf8+hjzs6Z+ewSk89+zpkzNsMwDAEAAPwkJNABAACA4EJyAAAATEgOAACACckBAAAwITkAAAAmJAcAAMCE5AAAAJiQHAAAABOSAwAAYEJyAAAATEgOAACACckBcIZatGiRevfuLbfb7bdzHjlyRCEhIXriiScsn2Pp0qVKSEhQeXm53+IC4B2SAzRbaWlpstlsdW65ublen3Pr1q1asGCBDh065P+Am1BJSYkefvhhzZ07VyEh5l8DH330ka6//np16dJFYWFhio2N1dChQ7Vw4cJa53G73erYsaMWLVokScrPz5dhGOrbt6/l2KZOnaqKigo9++yzls8BwDctAh0A0Njuvfde9ejRo1b72Wef7fW5tm7dqoULF2rq1Klq27atH6ILjBdffFFVVVWaMGGCqf21117ThAkTlJCQoFtuuUVxcXEqKCjQW2+9peXLl2v+/Pmm49977z19//33GjNmjCTp448/liT169fPcmwRERGaMmWKHn/8cc2YMUM2m83yuQBYQ3KAZm/06NG68MILm/y6paWlat26dZNftyGWLVumK6+8UhEREZ62H3/8UTfeeKOSkpK0ceNGhYWFeV679957tX///lrn2bBhg7p166bzzjtP0vHkoEOHDoqNjfUpvvHjx2vRokXatGmTLrvsMp/OBcB7DCvgjLdgwQLZbDbt3r3bUxFwOBy64YYbdPToUdNxc+bMkST16NHDMzyxb98+zzk+/fRTTZw4Ue3atdPQoUM9fT/44AONHj1adrtdbdq00a9+9atawxo15/j88881fvx42e12tW/fXrfddpvKyso8x23atEk2m00ZGRm13kt6erpsNpu2bdtW5/vdu3evPvroI40YMcLU/tZbb6m4uFjTpk0zJQY14uLiarW9/vrrnqqBdDw5qEkUajz//PMKCwvTzJkzVV1dXWdcJxowYICio6O1du3aBh0PwL+oHKDZKy4u1vfff29qs9lsat++valt/Pjx6tGjh1wul3bs2KG//vWviomJ0cMPPyxJ+s1vfqMvvvhC//jHP/TEE0+oQ4cOkqSOHTt6zjFu3Didc845evDBB2UYhiTpk08+0SWXXCK73a7bb79dLVu21LPPPqvk5GTl5ORo4MCBteLo3r27XC6XcnNz9ec//1k//vij/v73v0uSkpOTFR8frxUrVujqq6829V2xYoV69uypQYMG1fl5bN26VZL0y1/+0tReWloq6ficg4YoKirSBx98oHvvvdfT9vHHH3uGKqqqqjRz5kw999xzWrx4sW666aYGnbfGL3/5S/373//2qg8APzGAZmrZsmWGpJNu4eHhnuPmz59vSDJuvPFGU/+rr77aaN++vantkUceMSQZe/fuNbXXnGPChAm14khJSTHCwsKMPXv2eNr2799vREVFGcOGDat1jiuvvNLU/49//KMhyfjwww89bfPmzTPCw8ONQ4cOedq+/fZbo0WLFsb8+fPr/VzuvvtuQ5Jx+PBhU/u+ffuMVq1aGZKMc845x7j99tuNjRs3GlVVVSc9zwsvvGBERkYaR48e9bwnScbSpUuNH374wbjsssuM6OhoY9OmTfXGU5ebb77ZiIyMtNQXgG8YVkCzt3jxYmVlZZm2f/3rX7WOmz59umn/kksu0Q8//KCSkpIGX+vn56iurtabb76plJQUnXXWWZ72zp07a+LEidqyZUut8zudTtP+jBkzJB0f368xefJklZeXa/Xq1Z62VatWqaqqSr/73e/qjfGHH35QixYt1KZNG1N7t27dtG3bNo0fP1779+/XokWLdNlll+mss87SW2+9Ves8GzZs0PDhwxUZGSnpfxUHm82mpKQk7d+/X++++66Sk5Prjacu7dq107Fjx0xDOwCaBskBmr2LLrpII0aMMG3Dhw+vdVxCQoJpv127dpKOT9RrqJ/fFfHdd9/p6NGj6tWrV61j+/TpI7fbrcLCQlP7OeecY9rv2bOnQkJCtG/fPk9b7969lZSUpBUrVnjaVqxYoYsvvtjSXRg1zj//fK1atUoHDx7Upk2bdMMNN6iwsFC//e1vTX+kKysrlZWVVWu+gSSlpqaqU6dO2rZt20lj2bdvn2w2m9q0aaNWrVqpQ4cOuvPOO2sdZ/w0LMPdCkDTIzkAfhIaGnrS9po/Ug1R8y3an+r64zh58mTl5OTov//9r/bs2aPc3NxTVg0kqX379qqqqtLhw4frPCYsLEzJycl68cUXNWbMGP3www/atWuX5/WaiscVV1zhafv444/VrVs3XXLJJdqzZ4+OHDly0nN/+OGHOvfcc3XkyBEdPXpUb775plwul/Lz803H/fjjj2rVqlWjfKYA6kdyAHjB22+xHTt2VKtWrUx/WGt8/vnnCgkJUXx8vKn9yy+/NO3v3r1bbrdb3bt3N7X/9re/VWhoqP7xj39oxYoVatmypa677rpTxtS7d29Jx+9aaIjw8HBJksPh8LS9/vrrOvfcc00xffzxx+rfv79WrlypiIgIXX311aa7LGrUJAc1EhMT1bJlS1VWVpqO27t3r/r06dOgGAH4F8kB4IWadQsaukJiaGioRo4cqbVr15qGBQ4cOKD09HQNHTpUdrvd1Gfx4sWm/aefflrS8fUaTtShQweNHj1ay5cv14oVK3T55Zd77qCoT82dDNu3b/e0bdmyRceOHat17EcffaTMzExdcMEFpjkTGzZsMA0pVFdX67PPPlO/fv3UsWNHvfbaa8rPz9cf/vCHWuc8MTkoKyvT/fffrwEDBqh///6m43bs2KHBgwef8v0A8D9uZUSz969//Uuff/55rfbBgweb/uA1xIABAyRJd911l37729+qZcuWGjt2bL197r//fmVlZWno0KH64x//qBYtWujZZ59VeXm5Z9nhE+3du1dXXnmlLr/8cm3btk3Lly/XxIkTlZiYWOvYyZMn69prr5Uk3XfffQ16D2eddZb69u2rt956SzfeeKMk6Y477tAXX3yhcePGKTExUVVVVdq5c6deeuklORwOvfTSS6b4PvvsMy1ZssTT9uWXX6qsrMyzMuKAAQO0ZMkS3XDDDRowYIBSU1M9x3744YfKzMzUU089pSNHjqhLly7KyckxVWXy8vJ08OBBXXXVVQ16TwD8LNC3SwCNpb5bGSUZy5YtMwzjf7cQfvfddyft//PbFu+77z6jS5cuRkhIiOf1us5RY8eOHcaoUaOMNm3aGK1atTKGDx9ubN261XRMzTk+/fRT49prrzWioqKMdu3aGampqcaxY8dOet7y8nKjXbt2hsPhqPOYk3n88ceNNm3aeG5DfO2114wJEyYYZ599ttG6dWsjIiLC6NOnjzFnzhzj22+/NfV95plnDIfDYVRWVnraXn75ZUOS8cknn5iO/eMf/2i0bNnSyMnJMQzDMA4fPmzYbDbjyy+/NAzDMMrKyoypU6ca48aNM/WbO3eukZCQYLjd7ga/JwD+YzMML2ZbAWg0CxYs0MKFC/Xdd981aHhAOr7QUFxcnMaOHasXXnihwdcqLi7WWWedpUWLFmnatGlexXnFFVeoTZs2evnll73qJx1fgOmKK67Qjz/+6KkUPPfcc3rppZf0zjvvSJLKy8vVvXt33XHHHbrtttu8vgYA3zHnADiNrVmzRt99950mT57sVT+Hw6Hbb79djzzyiNePbE5OTtasWbO86lPjww8/1AUXXOBJDL788kstXrxYKSkpnmOWLVumli1b1lozAkDToXIABAlvKgfvvvuuPvroI913333q0KGDduzY0URR+mb
      "text/plain": [
       "<Figure size 640x480 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Here we calculate the entropy by differentiating the free energy\n",
    "\n",
    "Finterp = scipy.interpolate.RectBivariateSpline(Tp,rhop,Fgrid)\n",
    "SS = Finterp.partial_derivative(1,0)\n",
    "#dSdrho = Sinterp.partial_derivative(0,1)\n",
    "\n",
    "plt.clf()\n",
    "plt.imshow(SS(Tp,rhop)/ (10.0**Tgrid * np.log(10) * kB) * -1.0, aspect='auto', origin='lower', extent=(rhop[0],rhop[-1],Tp[0],Tp[-1]),\n",
    "           interpolation='none')\n",
    "#plt.imshow(Sgrid, aspect='auto', origin='lower', extent=(rhop[0],rhop[-1],Tp[0],Tp[-1]), \n",
    "#         interpolation='none')\n",
    "plt.title(r'Entropy ($S/k_B$)')\n",
    "plt.colorbar()\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "cc9dc872",
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgkAAAGzCAYAAACl24R2AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABP8klEQVR4nO3de3yUxb0/8M8mkIuQhGsIlwApFm9I5cdNwAtYBKly5JyWWorcRA9qUBCrgoqIFgKtWjy05dY2eE5NoWARiwqHcgxUgSoorVhFUZAUDaBIwFiSsPv8/qCseWY3O5mdeS67+3m/Xvt6sbvPMzP77JMwme/MdwKWZVkgIiIiEqR53QAiIiLyJ3YSiIiIKCp2EoiIiCgqdhKIiIgoKnYSiIiIKCp2EoiIiCgqdhKIiIgoKnYSiIiIKCp2EoiIiCgqdhKIiIgoKnYSiOKwcuVKBAKB8CMrKwvdu3fH1KlTceTIEa+bR0RkRBOvG0CUyB577DEUFRXh9OnTePXVV7FkyRK89NJL2Lt3L8477zyvm0dEpIWdBCINI0aMQJ8+fQAAt956K1q3bo2nnnoK69evx5gxYyKOr66uRrNmzdxuZtwsy8Lp06eRnZ3tdVOIyAMMNxAZdM011wAADhw4gIkTJ6J58+b48MMP8Z3vfAc5OTkYO3YsACAUCmHRokW45JJLkJWVhXbt2mHKlCn44osvbOXt2rULw4cPR5s2bZCdnY2ioiLccssttmNWrVqF3r17IycnB7m5ubj00kvx9NNPh99/9NFHEQgEItp6LmRy8ODB8Gtdu3bFDTfcgE2bNqFPnz7Izs7GsmXLAAAnTpzA9OnTUVhYiMzMTJx//vlYuHAhQqGQkWtHRP7DkQQigz788EMAQOvWrQEAZ86cwfDhw3HFFVfgiSeeCIcgpkyZgpUrV2LSpEm4++67ceDAAfz85z/HW2+9hddeew1NmzbF0aNHMWzYMLRt2xYzZ85EixYtcPDgQfzhD38I17d582aMGTMG3/72t7Fw4UIAwLvvvovXXnsN06ZNi+sz7Nu3D2PGjMGUKVNw22234YILLsBXX32Fq6++GocPH8aUKVPQuXNnbN++HbNmzcKnn36KRYsWaVw1IvIti4iUlZaWWgCsP/3pT9axY8esiooKa9WqVVbr1q2t7Oxs6x//+Ic1YcIEC4A1c+ZM27l//vOfLQDWs88+a3t948aNttfXrVtnAbDeeOONBtsxbdo0Kzc31zpz5kyDx8yZM8eK9qN+7jMcOHAg/FqXLl0sANbGjRttxz7++ONWs2bNrPfff9/2+syZM6309HTr0KFDDdZPRImL4QYiDUOHDkXbtm1RWFiIH/zgB2jevDnWrVuHjh07ho+54447bOesWbMGeXl5uPbaa/HZZ5+FH71790bz5s3xyiuvAABatGgBANiwYQPq6uqi1t+iRQtUV1dj8+bNxj5TUVERhg8fHtHmK6+8Ei1btrS1eejQoQgGg9i2bZux+onIP9hJINLwi1/8Aps3b8Yrr7yCv//97/joo49s/8E2adIEnTp1sp3zwQcfoKqqCvn5+Wjbtq3t8eWXX+Lo0aMAgKuvvhrf/e53MXfuXLRp0wY33ngjSktLUVNTEy7rzjvvRPfu3TFixAh06tQJt9xyCzZu3Kj1mYqKiiJe++CDD7Bx48aI9g4dOhQAwm0m8tq2bdswcuRIdOjQAYFAAM8//7yv6rv99tsRCAQSJkTHOQlEGvr16xde3RBNZmYm0tLsffFQKIT8/Hw8++yzUc9p27YtACAQCGDt2rXYuXMn/vjHP2LTpk245ZZb8OSTT2Lnzp1o3rw58vPzsWfPHmzatAkvv/wyXn75ZZSWlmL8+PF45plnwuVEEwwGo74ebSVDKBTCtddei/vvvz/qOd27d49+AYhcVl1djW9961u45ZZb8B//8R++qm/dunXYuXMnOnTo4Hi7TGEngchl3bp1w5/+9CcMGjSoUUsLL7/8clx++eWYN28eysrKMHbsWKxatQq33norACAjIwMjR47EyJEjEQqFcOedd2LZsmWYPXs2zj//fLRs2RLA2dUJ50IYAPDxxx8rtfnLL78MjxwQ+dWIESMwYsSIBt+vqanBQw89hN/97nc4ceIEevTogYULF2Lw4MGO1HfO4cOHcdddd2HTpk24/vrr46rLCww3ELns+9//PoLBIB5//PGI986cOYMTJ04AAL744gtYlmV7/7LLLgOAcMjh888/t72flpaGnj172o7p1q0bANjmDVRXV4dHGhrb5h07dmDTpk0R7504cQJnzpxpdFlEXpo6dSp27NiBVatW4W9/+xtGjx6N6667Dh988IFjdYZCIYwbNw733XcfLrnkEsfqcQJHEohcdvXVV2PKlCkoKSnBnj17MGzYMDRt2hQffPAB1qxZg6effhrf+9738Mwzz+CXv/wl/v3f/x3dunXDqVOnsGLFCuTm5uI73/kOgLMJnI4fP45rrrkGnTp1wscff4zFixfjsssuw0UXXQQAGDZsGDp37ozJkyfjvvvuQ3p6On7zm9+gbdu2OHToUKPafN999+GFF17ADTfcgIkTJ6J3796orq7G22+/jbVr1+LgwYNo06aNY9eMyIRDhw6htLQUhw4dCg/5/+hHP8LGjRtRWlqK+fPnO1LvwoUL0aRJE9x9992OlO8kdhKIPLB06VL07t0by5Ytw4MPPogmTZqga9euuPnmmzFo0CAAZzsTr7/+OlatWoUjR44gLy8P/fr1w7PPPhueXHjzzTdj+fLl+OUvf4kTJ06goKAAN910Ex599NHwXIimTZti3bp1uPPOOzF79mwUFBRg+vTpaNmyJSZNmtSo9p533nnYunUr5s+fjzVr1uC///u/kZubi+7du2Pu3LnIy8tz5kIRGfT2228jGAxGzKGpqakJ5zZ57733wh3shjzwwANYsGBBo+rcvXs3nn76abz55psNzg/ys4AljmcSERElgUAggHXr1mHUqFEAgNWrV2Ps2LF45513kJ6ebju2efPmKCgoQG1tLT766KOY5bZu3To8wThWfQCwaNEizJgxwzaBORgMIi0tDYWFhbaMp37EkQQiIkoJvXr1QjAYxNGjR3HllVdGPSYjIwMXXnihsTrHjRsXMeF3+PDhGDduXKNH8rzETgIRESWNL7/8Evv37w8/P3DgAPbs2YNWrVqhe/fuGDt2LMaPH48nn3wSvXr1wrFjx7Blyxb07NkzrlUHserr3LkzWrduHQ5lnNO0aVMUFBTgggsuiP+DuoSdBCIiShq7du3CkCFDws9nzJgBAJgwYQJWrlyJ0tJS/PjHP8a9996Lw4cPo02bNrj88stxww03OFJfolOak7BkyRIsWbIkHEO55JJL8Mgjj8RcI7pmzRrMnj0bBw8exDe/+U0sXLgwPDObiIiI/EspT0KnTp2wYMEC7N69G7t27cI111yDG2+8Ee+8807U47dv344xY8Zg8uTJeOuttzBq1CiMGjUKe/fuNdJ4IiIico726oZWrVrhpz/9KSZPnhzx3k033YTq6mps2LAh/Nrll1+Oyy67DEuXLtWploiIiBwW95yEYDCINWvWoLq6GgMGDIh6zI4dO8LxmXOGDx8u3QCjpqbGtolNKBTC8ePH0bp164RcZ0pERO6xLAunTp1Chw4dIvZOMeX06dOora01UlZGRgaysrKMlGWacifh7bffxoABA3D69OnwtrgXX3xx1GMrKyvRrl0722vt2rVDZWVlzDpKSkowd+5c1aYRERGFVVRUROzCasLp06dR1KU5Ko9G3yRNVUFBAQ4cOODLjoJyJ+GCCy7Anj17UFVVhbVr12LChAnYunVrgx2FeMyaNcs2AlFVVYXOnTuj93UPIb1pAxdRHGBQDaIonp92xn5AMNNeQMDMvVOvQIVjxbarDr54fb5u+arvJxPVayGSXZtUupaNYfp6mmT6u0mgn7Ng3Wns3jgPOTk5jpRfW1uLyqNBHNjdBbk5eiMVJ0+FUNT7Y9TW1iZHJyEjIwPnn38+AKB3795444038PTTT2PZsmURxxYUFODIkSO2144cOYKCgoKYdWR
      "text/plain": [
       "<Figure size 640x480 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAhcAAAGzCAYAAACcvDUtAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAACqJElEQVR4nO2deXiU1dn/70kgySRkJvtKSAIEAoR9K7hhXXCpSmttX21F/Vmr74tWa9/WYrVYbUWr1rbW2trXQt/XUq2tS6tWRSoqAsouARJIyL5vM5NlsjAzvz9Sk7m/T+YcJjMkmfT+XFeui5NznmWeeZ6Hk3N/7+9t8ng8HhIEQRAEQQgSYaN9AoIgCIIgjC9kciEIgiAIQlCRyYUgCIIgCEFFJheCIAiCIAQVmVwIgiAIghBUZHIhCIIgCEJQkcmFIAiCIAhBRSYXgiAIgiAEFZlcCIIgCIIQVGRyIQiCIAhCUJHJhSAEgc2bN5PJZPL5s3v3br/29+abb9IDDzxwZk5WEAThDDNhtE9AEMYTDz74IOXm5hp+P336dL/28+abb9LTTz8tEwxBEEISmVwIQhC59NJLacmSJSN6zFOnTpHb7aaIiIgRPa4gCIIvJCwiCCNEeXk5mUwmevzxx+nZZ5+ladOmUWRkJC1dupT27NkzMO7GG2+kp59+moiIhVZwHz//+c8H9nH06FEiIvrnP/9J55xzDsXExFBcXBxdddVVdOzYMXYeDzzwAJlMJioqKqKvfOUrZLFYKDExke68807q7u4eGHfeeefR/Pnzh/wsM2fOpNWrVwf1+giCMH6QlQtBCCJ2u52am5vZ70wmEyUmJg60t2zZQu3t7XTrrbeSyWSin/70p/SlL32JTp48SRMnTqRbb72VamtraevWrfR///d/Qx5n06ZN1N3dTd/85jcpMjKSEhIS6N1336VLL72Upk6dSg888AA5nU566qmn6KyzzqL9+/dTTk4O28dXvvIVysnJoY0bN9Lu3bvpl7/8JbW1tdH//u//EhHR9ddfT7fccgsVFhZSQUHBwHZ79uyh48eP03333RekqyYIwrjDIwhCwGzatMlDREP+REZGejwej6esrMxDRJ7ExERPa2vrwLavvfaah4g8f//73wd+t27dOs9Qj+dn+7BYLJ7GxkbWt2DBAk9KSoqnpaVl4HeHDh3yhIWFedauXTvwuw0bNniIyHPllVey7f/rv/7LQ0SeQ4cOeTwej8dms3mioqI899xzDxv3rW99yxMTE+Pp6Ojw9zIJgvBvgqxcCEIQefrpp2nGjBnsd+Hh4az91a9+leLj4wfa55xzDhERnTx58rSPc/XVV1NycvJAu66ujg4ePEjf+973KCEhYeD38+bNo4suuojefPNNwz7WrVvH2nfccQf9+te/pjfffJPmzZtHVquVrrrqKvrTn/5EGzduJJPJRC6Xi1588UVas2YNxcTEnPb5CoLw74VoLgQhiCxbtowuvPBC9nP++eezMVOmTGHtzyYabW1tp30czEipqKggon4tBDJr1ixqbm6mzs5O9vu8vDzWnjZtGoWFhVF5efnA79auXUuVlZX04YcfEhHRu+++Sw0NDXT99def9rkKwkjywQcf0BVXXEEZGRlkMpno1VdfHRPHO3bsGF155ZVktVopJiaGli5dSpWVlWf03EYTmVwIwgiDKxmf4fF4TnsfZrM5WKczwGeiUW9Wr15Nqamp9PzzzxMR0fPPP09paWl04YUXBv34ghAMOjs7af78+QOi6LFwvNLSUjr77LMpPz+ftm/fTp9++indf//9FBUVNSLnOBpIWEQQxiBD/UevIjs7m4iIiouLDX1FRUWUlJRkCGOcOHGCrYCUlJSQ2+1mws/w8HC67rrraPPmzfToo4/Sq6++SrfccovPCZIgjDaXXnopXXrppT77e3p66Ac/+AH96U9/IpvNRgUFBfToo4/SqlWrzsjxiIh+8IMf0GWXXUY//elPB343bdq0YR0vVJCVC0EYg3w2EbDZbKc1Pj09nRYsWEB/+MMf2DaFhYX0zjvv0GWXXWbYBv/Seuqpp4iIDC/K66+/ntra2ujWW2+ljo4O+vrXv+7HJxGEscXtt99Ou3btohdeeIE+/fRTuuaaa+iSSy6hEydOnJHjud1ueuONN2jGjBm0evVqSklJoeXLl5/xcM1oIysXghBE/vGPf1BRUZHh9ytXrqSwsNOfyy9evJiIiL71rW/R6tWrKTw8nP7jP/5Duc1jjz1Gl156Ka1YsYJuvvnmgVRUq9U6pNNnWVkZXXnllXTJJZfQrl276Pnnn6frrrvO4G2xcOFCKigooJdeeolmzZpFixYtOu3PIQhjicrKStq0aRNVVlZSRkYGERH993//N7311lu0adMmevjhh4N+zMbGRuro6KBHHnmEfvzjH9Ojjz5Kb731Fn3pS1+i9957j84777ygH3NMMNrpKoIwHlClohKRZ9OmTQNppI899phheyLybNiwYaB96tQpzx133OFJTk72mEymgbRU1T48Ho/n3Xff9Zx11lkes9nssVgsniuuuMJz9OhRNuazVNSjR496vvzlL3tiY2M98fHxnttvv93jdDqH3O9Pf/pTDxF5Hn744WFeIUEYeYjI88orrwy0X3/9dQ8ReWJiYtjPhAkTPF/5ylc8Ho/Hc+zYMeWzTESG9Gxfx/N4PJ6amhoPEXmuvfZa9vsrrrjC8x//8R9B/bxjCVm5EIQgcOONN9KNN96oHefxIdrE34eHh9Mvf/lL+uUvf8l+n5OToxR+XnDBBXTBBRfoT5iIkpOT6aWXXjqtsREREWQymehrX/vaaY0XhLFIR0cHhYeH0759+wy6oUmTJhER0dSpUw2utoi3KZ6OpKQkmjBhAs2ePZv9ftasWbRjx47T3k+oIZMLQRCUeDweeu655+i8884zpNEKQiixcOFCcrlc1NjYOOAvg0RERFB+fn7QjhkREUFLly41iK2PHz8+IMQej8jkQhCEIens7KS//e1v9N5779Hhw4fptddeG+1TEgQtHR0dVFJSMtAuKyujgwcPUkJCAs2YMYO+9rWv0dq1a+mJJ56ghQsXUlNTE23bto3mzZtHl19+eVCP99lk/Lvf/S599atfpXPPPZfOP/98euutt+jvf/87bd++PeDPO2YZ3aiMIAgjzWeai6amJuW4z/QdcXFxnnvvvXeEzk4QAuO9994bUidxww03eDwej6e3t9fzwx/+0JOTk+OZOHGiJz093fPFL37R8+mnn56R433Gc88955k+fbonKirKM3/+fM+rr74a4Ccd25g8ntN37nnmmWfomWeeGXDwmzNnDv3whz9U5vi+9NJLdP/991N5eTnl5eXRo48+OmRanCAIgiAI4wO/fC4mT55MjzzyCO3bt4/27t1Ln//85+mqq66iI0eODDl+586ddO2119LNN99MBw4coDVr1tCaNWuosLAwKCcvCIIgCMLYw6+Vi6FISEigxx57jG6++WZD31e/+lXq7Oyk119/feB3n/vc52jBggX0m9/8JpDDCoIgCIIwRhm2oNPlctFLL71EnZ2dtGLFiiHH7Nq1i+6++272u9WrV2udyXp6eqinp2eg7Xa7qbW1lRITE/22RRYEQRD+vfB4PNTe3k4ZGRl+mdf5Q3d3N/X29gZlXxEREeOuzojfk4vDhw/TihUrqLu7myZNmkSvvPKKIX/3M+rr6yk1NZX9LjU1lerr65XH2LhxI/3oRz/y99QEQRAEYYCqqiqaPHly0Pfb3d1NudmTqL7RFZT9paWlUVlZ2biaYPg9uZg5cyYdPHiQ7HY7/eUvf6EbbriB3n//fZ8TjOGwfv16tuJht9tpypQplPHY9ynM3H/xo6zdyn102/mXFG7nH9VlPcU3gAURM+zfaVd/6eE22H/cKR8j+4mO4/vPsNhZu9Zh5adnGoxeddrUFTFj4pzKfu99DYXzpIW1XfH8s0yK71KOj2jjfyn0xrv5/hL6WDsWzvfUvjjl+RmAj9Ob4B563GliyeOlzx0n4lnbncjPHwlrmcjHJwX2101Yc4TP/WFf4swW1u7YmczaPYnqaxP
      "text/plain": [
       "<Figure size 640x480 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAf4AAAGzCAYAAADOqXWXAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAuXElEQVR4nO3df3hU1Z3H8c8EySRCEgETEkggKSqoCCIiAkKDIorIGrelfRAN+LCtuAMF8Wf8hVhg8AePuNRGRJt0t6ahYAMtIhipCVJAIcBK2Brll4lIAEUDREnCzN0/kNFJSDJzZzI/nPfrec6ze+/cc+6ZKe0353vOPddiGIYhAAAQEaKC3QEAABA4BH4AACIIgR8AgAhC4AcAIIIQ+AEAiCAEfgAAIgiBHwCACELgBwAgghD4AQCIIAR+IEQ89dRTslgsfmsvPT1dkydP9lt7AH4cCPxAiEpPT5fFYnGVpKQkDR8+XEVFRT631VzJz8/3/xcBEFLOC3YHADTvyiuv1P333y9J+vzzz7VkyRL9+7//u3JzczV16lSP21m0aJFOnjzpOl6zZo3+/Oc/64UXXtCFF17oOj906FD/dR5ASCLwAyGse/fuuvPOO13H2dnZuuiii/TCCy94FfizsrLcjqurq/XnP/9ZWVlZSk9P91NvAYQDUv1AEGzcuFGDBg1STEyMevXqpSVLlnhULzk5WZdeeqn279/vOmcYhubOnavU1FSdf/75GjlypHbv3t1WXQcQ5hjxAwG2a9cujR49WomJiXrqqad0+vRpzZ49W127dm21bkNDg6qqqtSlSxfXuSeffFJz587VLbfcoltuuUXbt2/X6NGjVV9f35ZfA0CYIvADAfbkk0/KMAy999576tGjhyTpZz/7ma644oom1zY0NOiLL76QdGaO32636/Dhw5o+fbok6ejRo3r22Wc1duxY/f3vf3c9FfDYY49p/vz5AfpGAMIJqX4ggBwOh9atW6esrCxX0JekSy+9VDfddFOT699++20lJiYqMTFR/fv31/Lly3XXXXfpmWeekSS98847qq+v1/Tp090eBZw5c2abfxcA4YnADwTQ0aNH9e233+riiy9u8lnv3r2bnBs8eLCKi4v1zjvvaNOmTfriiy/03//934qNjZUkffrpp5LUpL3ExER16tSpDb4BENo2bNigcePGqVu3brJYLFq5cmXA7r1gwQJZLJYmf3jfc8896tWrl2JjY5WYmKjbbrtNH330UcD61RiBHwhhF154oUaNGqUbbrhBQ4YM0QUXXBDsLgEhrba2Vv3799dLL70U0Ptu3bpVS5YsUb9+/Zp8NnDgQOXl5elf//qX1q1bJ8MwNHr0aDkcjoD28SwCPxBAiYmJio2N1SeffNLks4qKCq/b69mzpyQ1ae/o0aP66quvzHUSCGNjxozR3Llzdfvtt5/z87q6Oj3wwAPq3r27OnTooMGDB6ukpMSne548eVITJ07U0qVLz5lp+/Wvf60RI0YoPT1dV111lebOnauqqiodOHDAp/uaReAHAqhdu3a66aabtHLlSlVWVrrOnx0JeGvUqFFq3769Fi9eLMMwXOcXLVrkj+4CPzrTpk3T5s2bVVhYqA8//FDjx4/XzTfffM4/xj1ls9k0duxYjRo1qtVra2trlZeXp4yMDKWlpZm+py9Y1Q8E2Jw5c7R27VoNHz5c//mf/6nTp09r8eLFuvzyy/Xhhx961VZiYqIeeOAB2e123Xrrrbrlllu0Y8cOvfXWW2478gGQKisrlZeXp8rKSnXr1k2S9MADD2jt2rXKy8sz9SRMYWGhtm/frq1bt7Z43e9//3s99NBDqq2tVe/evVVcXKzo6GhT38NXjPiBAOvXr5/WrVunxMREPfnkk/rDH/6gOXPmNJuabM3cuXM1Z84c7dixQw8++KD27t2rt99+Wx06dPBzz4HwtmvXLjkcDl1yySXq2LGjq5SWlmrv3r2SpI8++qjVd1o88sgjkqSqqirNmDFDr7/+umJiYlq898SJE7Vjxw6Vlpbqkksu0S9+8QudOnWqzb/zuViMH+YHAQD4kbBYLCoqKnJtWb1s2TJNnDhRu3fvVrt27dyu7dixo5KTk1VfX699+/a12G6XLl2UmJiolStX6vbbb3dry+FwyGKxKCoqSnV1dU3uI0n19fXq1KmTXn31VU2YMMH3L+olUv0AgIgwYMAAORwOHTlyRMOHDz/nNdHR0erTp49H7d1www3atWuX27m7775bffr00cMPP3zOoC+d2WbbMAzV1dV59wX8hMAPAPjROHnypPbs2eM63r9/v3bu3KnOnTvrkksu0cSJE5Wdna2FCxdqwIABOnr0qNavX69+/fpp7NixXt0rLi5Offv2dTvXoUMHdenSxXV+3759WrZsmWub7s8++0wLFixQbGysbrnlFt+/sAkEfgDAj8a2bds0cuRI1/GsWbMkSZMmTVJ+fr7y8vI0d+5c3X///Tp48KAuvPBCXXvttbr11lvbpD8xMTF67733tGjRIn311Vfq2rWrRowYoU2bNikpKalN7tkar+b4c3NzlZub63r28PLLL9eTTz6pMWPGNFtn+fLleuKJJ3TgwAFdfPHFeuaZZ4L2Vw4AAJHOq1X9qampWrBggcrKyrRt2zZdf/31uu2225p9BeimTZs0YcIETZkyRTt27FBWVpaysrJUXl7ul84DAADv+Lyqv3Pnznruuec0ZcqUJp/98pe/VG1trVavXu06d+211+rKK6/Uyy+/7MttAQCACabn+B0Oh5YvX67a2loNGTLknNds3rzZNb9y1tldy1pSV1fnttrR6XTq2LFj6tKli9sbyAAAaMwwDJ04cULdunVTVFTbbFdz6tQp1dfX+6Wt6OjoVvcB8CevA/+uXbs0ZMgQnTp1Sh07dlRRUZEuu+yyc15bXV2trl27up3r2rWrqqurW7yH3W7XnDlzvO0aAAAuVVVVSk1N9Xu7p06dUkbPjqo+4p+X7CQnJ2v//v0BC/5eB/7evXtr586dqqmp0YoVKzRp0iSVlpY2G/zNyMnJccsU1NTUqEePHvp0e7riO5756+2vJ+Pd6qz/yv3+u4+luB0f+7Kj27Hlq/Zux+1r3P8qjD7h3qf2J91nRNp/43583rdOt+N2p9yPz/vG/R9IVN1p9+MG9+st9e6fq+H7Y8vpRp+dbvSPz+HeltH4c2drx+7fzWj8BqnGs0NOZ4sfN9GofRnOc1/n+tjL2ahW2vMae1wBYeW0GrRRaxQXF9cm7dfX16v6iEP7y3oqPs63jMLxE05lDPxU9fX1oRv4o6OjddFFF0k686rBrVu36sUXX9SSJUuaXJucnKzDhw+7nTt8+LCSk5NbvIfVapXVam1yPr5jlOtHPt/ivjFCdIP7nsftTrnXj/rG/Qe1fOse+Nudcv8Pr12jDE67evf/8W/X0Cjwn24U+Bsdn3deo8DvaBT4GwVfS+ONH5zfH1ucjT6LahyY3Y+Nxp+r5etlaRT4LY3rNwqElkaBv/HnjVkaf95K4G9yfWv8HPhb+z4AQst3/5Vt66nh+LgonwN/MPjcY6fT2ezuQ0OGDNH69evdzhUXFze7JgAAgHDhMJx+KYHm1Yg/JydHY8aMUY8ePXTixAkVFBSopKTE9TrR7Oxsde/eXXa7XZI0Y8YM/fSnP9XChQs1duxYFRYWatu2bXrllVf8/00AAAggpww5fcwI+lrfDK8C/5EjR5Sdna1Dhw4pISHB9ZaxG2+8UdKZVx7+cAXl0KFDVVBQoMcff1yPPvqoLr74Yq1cubLJFocAAIQbp5w+Tyz63oL3vAr8r732Woufl5SUNDk3fvx4jR8/3qtOAQCAtsFe/QAAmOAwDDl8fOrH1/pmEPgBADAhXOf4w+85BAAAYBojfgAATHDKkCMMR/wEfgAATCDVDwAAQh4jfgAATGBVPwAAEcQp398MEvjte0j1AwAQURjxAwBggsMPq/p9rW8GgR8AABMcxpniaxuBRuAHAMAE5vgBAEDII/ADAGCCUxY5fCxOWby+78GDB3XnnXeqS5cuio2N1RVXXKFt27Z5XJ9UPwAAJjiNM8X
      "text/plain": [
       "<Figure size 640x480 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgoAAAGzCAYAAABO7D91AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAxrUlEQVR4nO3de3QUVbr38V8nJJ0ASQzYgQAhXFRwRBQRIyA3h+EyHjQuZTyIIC4Oim/H4aZCFEQcoHHAGRQURZ3EdY4MKAxwTlSQQQIioIigwkAcQEwMF3GQwETJpbveP5Ce6UCFdHUn3aG/n7X2kqquvWt3CeTheXZV2QzDMAQAAHABUaGeAAAACF8ECgAAwBSBAgAAMEWgAAAATBEoAAAAUwQKAADAFIECAAAwRaAAAABMESgAAABTBApAHXn66adls9lCMl5ubq5sNps+/fTToJ0fQGQgUABCpLy8XM8//7y6dOmixMREXXbZZbrmmmv04IMPat++faGeHgBIkhqEegJApLrrrrv03nvvadiwYRozZowqKiq0b98+5eXlqUePHurYsWOopwgABApAKGzfvl15eXmaNWuWnnjiCZ/PFi5cqJMnTwblPB6PR+Xl5UEZC0BkovQA1ILNmzerW7duiouLU/v27fXKK6/4fH7gwAFJUs+ePc/rGx0draZNm/o13jk2m01ZWVl68803dc0118hut2vNmjXez8vKyjRx4kQ5HA41atRId955p44fP37eOC+99JK3f4sWLeR0OoMWvACoX2y8ZhoIri+//FIZGRlyOBx6+OGHVVlZqYULF6pZs2b64osvZBiGtm7dqh49emjMmDF66aWX1KCBeXKvJuOdY7PZdPXVV+v7779XVlaWLr/8cvXo0UO7du3SAw88oC5duig5OVl33nmnDh06pPnz5+uuu+7SsmXLvGM8/fTTmjFjhvr376877rhDBQUFWrRokW644QZ99NFHiomJqdXrByDMGACCKjMz04iLizO++eYb776//e1vRnR0tHHuj5zH4zH69OljSDKaNWtmDBs2zHjxxRd9+vgz3jmSjKioKGPPnj0++3NycgxJRv/+/Q2Px+PdP2HCBCM6Oto4efKkYRiG8d133xmxsbHGgAEDDLfb7T1u4cKFhiTjT3/6UwBXBkB9ROkBCCK32621a9cqMzNTrVu39u6/+uqrNXDgQO+2zWbT2rVrNXPmTCUnJ+vPf/6znE6n0tPTdc8993jT/DUd79/16dNHv/jFLy742YMPPuhzS2WvXr3kdrv1zTffSJL++te/qry8XOPHj1dU1L/+ehgzZowSExP1zjvv+H9RANRrBApAEB0/flw//fSTrrzyyvM+69Chg8+23W7Xk08+qb179+rw4cP685//rJtvvllvvfWWsrKy/B7vnLZt25rO79+DDUlKTk6WJP3www+S5A0Yqo4dGxurdu3aeT8HwtWmTZs0ZMgQtWjRQjabTatWraqzc8+ZM0c2m03jx4/32f/QQw+pffv2io+Pl8Ph0B133FGvboEmUADCQGpqqv7zP/9TmzZt0pVXXqm33npLlZWVlsaKj483/Sw6OvqC+w2WKuESUVpaquuuu04vvvhinZ53+/bteuWVV9S5c+fzPuvatatycnK0d+9erV27VoZhaMCAAXK73XU6R6sIFIAgcjgcio+P19///vfzPisoKLho/5iYGHXu3FkVFRX6/vvvAx7PX+np6Rccu7y8XF9//bX3cyBcDR48WDNnztSdd955wc/Lysr06KOPqmXLlmrUqJEyMjKUn58f0Dn/+c9/avjw4Xr11Ve9Wbp/9+CDD6p3795q06aNbrjhBs2cOVNFRUU6dOhQQOetKwQKQBBFR0dr4MCBWrVqlQoLC737z/1L4py///3vPp+fc/LkSW3dulXJyclyOBw1Hi9Y+vfvr9jYWL3wwgs+WYbXX39dJSUluu2224J+TqAuZWVlaevWrVq6dKm++OILDR06VIMGDbpgMF5TTqdTt912m/r373/RY0tLS5WTk6O2bdsqLS3N8jnrEg9cAoJsxowZWrNmjXr16qX/9//+nyorK7VgwQJdc801+uKLLyRJn3/+ue69914NHjxYvXr1UpMmTVRcXKw33nhDhw8f1vz5871lgpqMFywOh0PZ2dmaMWOGBg0apNtvv10FBQV66aWX1K1bN913331BPR9QlwoLC5WTk6PCwkK1aNFCkvToo49qzZo1ysnJ0ezZs/0ec+nSpfrss8+0ffv2ao976aWX9Pjjj6u0tFQdOnTQunXrFBsba+l71LkQ33UBXJI2btxodO3a1YiNjTXatWtnvPzyy8b06dO9tzMeO3bMmDNnjtGnTx8jNTXVaNCggZGcnGzceuutxvLly/0e7xxJhtPpPK//udsjt2/f7rN/w4YNhiRjw4YNPvsXLlxodOzY0YiJiTGaNWtmPPzww8YPP/wQ2EUB6pgkY+XKld7tvLw8Q5LRqFEjn9agQQPjN7/5jWEYhrF3715DUrVt8uTJhmEYRmFhoZGSkmJ8/vnn3nP06dPHGDdu3HlzOXnypPHVV18ZGzduNIYMGWLccMMNxk8//VSr3z9YeOASAOCSZLPZtHLlSmVmZkqSli1bpuHDh2vPnj3nLext3LixmjdvrvLych08eLDacZs2bSqHw6FVq1bpzjvv9BnL7XbLZrMpKipKZWVlF1xAXF5eruTkZL322msaNmxY4F+0llF6AABEhC5dusjtduu7775Tr169LnhMbGxsjV/I9stf/lJffvmlz74HHnhAHTt21OTJk6u9y8gwDJWVlfn3BUKEQAEAcMn45z//qf3793u3v/76a+3atUtNmjTRVVddpeHDh2vkyJF67rnn1KVLFx0/flzr169X586d/V6sm5CQoE6dOvnsa9SokZo2berdf/DgQS1btkwDBgyQw+HQt99+qzlz5ig+Pl6//vWvA//CdYBAAQBwyfj000/Vr18/7/bEiRMlSffff79yc3OVk5OjmTNnatKkSSouLtbll1+um2++Wf/xH/9RK/OJi4vThx9+qPnz5+uHH35Qs2bN1Lt3b23ZskUpKSm1cs5g82uNwqJFi7Ro0SLvvZ/XXHONnnrqKQ0ePNi0z9tvv61p06bp0KFDuvLKK/Xss8/WmygKAIBI59dzFFq1aqU5c+Zox44d+vTTT3Xrrbfqjjvu0J49ey54/JYtWzRs2DCNHj1aO3fuVGZmpjIzM7V79+6gTB4AANSugO96aNKkiebOnavRo0ef99k999yj0tJS5eXleffdfPPNuv766/Xyyy8HcloAAFAHLK9RcLvdevvtt1VaWqru3btf8JitW7d660PnnHvKXHXKysp8VoN6PB6dOHFCTZs29XnzHQAAVRmGodOnT6tFixY+b0ENpjNnzqi8vDwoY8XGxiouLi4oY9UGvwOFL7/8Ut27d9eZM2fUuHFjrVy50vSVtkePHlWzZs189jVr1kxHjx6t9hwul0szZszwd2oAAHgVFRWpVatWQR/3zJkzapveWEe/C85LnZo3b66vv/46bIMFvwOFDh06aNeuXSopKdHy5ct1//33a+PGjabBghXZ2dk+mYiSkhK1bt1at+jXaqAYSVLUL67y6VPaNsF3u7nv/atnLvfNRJQne3y2PUkVPttxjX3vb72s0U8+203iq2zHlvpsJ8f86Pt5jO/niQ18+ydEVd32PX+87V/bjaJ8o1i7zfctgzE23+8WZ/P9zdxAvtWmmCpJmqp3/sZUyeJEy3e7apan6udVRVVZGhN9kSxR1eMvJuoi5/dXtI1XogD1yal/epR+wyElJCRc/GALysvLdfQ7t77eka7EhMD+fjh12qO2Xb9ReXn5pRMoxMbG6oorrpB09tWZ27dv1/PPP69XXnnlvGObN2+uY8eO+ew7duyYmjdvXu057Ha77Hb7BSYbowa2nwOFaN/PG8T4XuDoWN8fd9F23x8eUXG+P0wVX+X4hr7HN2jke3xMvO8P31i7b6Bhj6kSeMTE+J6uge/n8VG+/ysaRvv+8G9o+9f8GlVJpcVV+UFW9Qd/nK1KYBDkQCGKQAFAGKr
      "text/plain": [
       "<Figure size 640x480 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgsAAAGzCAYAAAChLlRLAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAC1BklEQVR4nO29eZgfVZX/f6rqs/benaS7sycQ9ggCAgaUZQZBZJyJ8x101GH7iTPMExwQRx3cWFSig4h8GQdUlMyofMFlkBnwwWHAwDCACsooiGHLnnSSTu/LZ6uq3x8J3Z/zqnR90p1Ouj/kvJ4nz5PbVXXr1q17q2/fc877OGEYhmIYhmEYhjEG7lQ3wDAMwzCM6Y0tFgzDMAzDiMUWC4ZhGIZhxGKLBcMwDMMwYrHFgmEYhmEYsdhiwTAMwzCMWGyxYBiGYRhGLLZYMAzDMAwjFlssGIZhGIYRiy0WDGMac91114njOFPdDMMwDnJssWAYIvLqq6/K3/zN38ghhxwimUxGGhoa5LTTTpNbb71VhoeHp7p542LLli1y3XXXyXPPPTfVTTEM4w2CY7khjIOdBx98UC644AJJp9Ny0UUXydKlS6VQKMgTTzwhP/7xj+WSSy6Rb37zm1PStlKpJKVSSTKZzF5f88wzz8hJJ50kd911l1xyySX7r3GGYRw0JKa6AYYxlaxdu1b+8i//UhYuXCiPPvqozJ49e+TYihUr5JVXXpEHH3xwytqXSCQkkbBpahjG1GJmCOOg5h//8R9lYGBAvv3tb6uFwussWbJErrzyShHZ9Vf+5z//eTn00EMlnU7LokWL5FOf+pTk83l1zaJFi+RP/uRP5IknnpCTTz5ZMpmMHHLIIfKv//qv6rxisSjXX3+9HHbYYZLJZGTGjBnytre9TR5++OGRc/bks/Dwww/L2972NmlqapK6ujo54ogj5FOf+pSIiKxevVpOOukkERG59NJLxXEccRxHVq1aNXL9L37xC3nnO98pjY2NUlNTI2eccYb8z//8j7rH6/d95ZVX5JJLLpGmpiZpbGyUSy+9VIaGhiL99L3vfU9OPvlkqampkebmZjn99NPlP//zP0VE5OKLL5aZM2dKsViMXHfOOefIEUccEfm5YRjTC1ssGAc1//Ef/yGHHHKInHrqqRXPveyyy+Rzn/ucnHDCCXLLLbfIGWecIStXrpS//Mu/jJz7yiuvyF/8xV/IO97xDrn55pulublZLrnkEnnhhRdGzrnuuuvk+uuvl7POOkv+6Z/+ST796U/LggUL5Ne//vWYbXjhhRfkT/7kTySfz8sNN9wgN998s/zpn/7pyC/7o446Sm644QYREfnrv/5r+e53vyvf/e535fTTTxcRkUcffVROP/106evrk2uvvVZuvPFG6enpkT/6oz+SX/7yl5H7vfe975X+/n5ZuXKlvPe975VVq1bJ9ddfr865/vrr5cILL5RkMik33HCDXH/99TJ//nx59NFHRUTkwgsvlJ07d8rPfvYzdV1HR4c8+uij8ld/9VcV+94wjCkmNIyDlN7e3lBEwj/7sz+reO5zzz0Xikh42WWXqZ///d//fSgi4aOPPjrys4ULF4YiEj7++OMjP9u+fXuYTqfDj33sYyM/O+6448Lzzz8/9r7XXnttWD5Nb7nlllBEwh07dox5za9+9atQRMK77rpL/TwIgvCwww4Lzz333DAIgpGfDw0NhYsXLw7f8Y53RO77//1//5+q4z3veU84Y8aMkfLLL78cuq4bvuc97wl934/cLwzD0Pf9cN68eeH73vc+dfyrX/1q6DhO+Nprr8X2gWEYU4/tLBgHLX19fSIiUl9fX/Hcn/70pyIicvXVV6uff+xjHxMRifg1HH300fL2t799pDxr1iw54ogj5LXXXhv5WVNTk7zwwgvy8ssv73Wbm5qaRETk/vvvlyAI9vo6EZHnnntOXn75ZfnABz4gO3fulM7OTuns7JTBwUH54z/+Y3n88ccjdV5++eWq/Pa3v1127tw50nc/+clPJAgC+dznPieuqz8nr5tPXNeVD37wg/Lv//7v0t/fP3L8+9//vpx66qmyePHicT2HYRgHHlssGActDQ0NIiLqF9hYrF+/XlzXlSVLlqift7e3S1NTk6xfv179fMGCBZE6mpubpbu7e6R8ww03SE9Pjxx++OHypje9ST7+8Y/Lb3/729h2vO9975PTTjtNLrvsMmlra5O//Mu/lB/84Ad7tXB4fVFy8cUXy6xZs9S/O++8U/L5vPT29sY+R3Nzs4jIyHO8+uqr4rquHH300bH3vuiii2R4eFjuu+8+ERFZs2aNPPvss3LhhRdWbLdhjIfHH39c3v3ud8ucOXPEcRz5yU9+Mq3ud/nll4vjOPK1r31tv7ZrsrHFgnHQ0tDQIHPmzJHnn39+r6/ZW4Ekz/P2+POwLFL59NNPl1dffVW+853vyNKlS+XOO++UE044Qe68884x681ms/L444/Lf/3Xf8mFF14ov/3tb+V973ufvOMd7xDf92Pb9PqC4qabbpKHH354j//q6urG/Rx7w9FHHy0nnniifO973xORXQ6RqVRK3vve946rHsOoxODgoBx33HHy9a9/fdrd77777pOnn35a5syZcwBaNrnYYsE4qPmTP/kTefXVV+Wpp56KPW/hwoUSBEHEZLBt2zbp6emRhQsXTuj+LS0tcumll8r/+3//TzZu3CjHHnusXHfddbHXuK4rf/zHfyxf/epX5fe//7188YtflEcffVR+/vOfi8jYC5pDDz1URHYtks4+++w9/ksmk+Nq/6GHHipBEMjvf//7iudedNFF8uijj8rWrVvl7rvvlvPPP39kp8IwJovzzjtPvvCFL8h73vOePR7P5/Py93//9zJ37lypra2VU045RVavXr3f7vc6mzdvlo985CPy/e9/f9zzbDpgiwXjoOYTn/iE1NbWymWXXSbbtm2LHH/11Vfl1ltvlXe9610iIpGtw69+9asiInL++eeP+947d+5U5bq6OlmyZEkkFLOcrq6uyM/e/OY3i4iMXFdbWysiIj09Peq8E088UQ499FD5yle+IgMDA5F6duzYMZ7mi4jI8uXLxXVdueGGGyKmEO4+vP/97xfHceTKK6+U1157zaIgjCnhiiuukKeeekruuece+e1vfysXXHCBvPOd7xyX79B4CYJALrzwQvn4xz8uxxxzzH67z/7E1F6Mg5pDDz1U7r77bnnf+94nRx11lFJwfPLJJ+WHP/yhXHLJJXLllVfKxRdfLN/85jelp6dHzjjjDPnlL38p//Iv/yLLly+Xs846a9z3Pvroo+XMM8+UE088UVpaWuSZZ56RH/3oR3LFFVeMec0NN9wgjz/+uJx//vmycOFC2b59u/zzP/+zzJs3T972treNPFNTU5PccccdUl9fP/LX0+LFi+XOO++U8847T4455hi59NJLZe7cubJ582b5+c9/Lg0NDfIf//Ef43qGJUuWyKc//Wn5/Oc/L29/+9vlz//8zyWdTsuvfvUrmTNnjqxcuXLk3FmzZsk73/lO+eEPfyhNTU0TWmAZxr6wYcMGueuuu2TDhg0jpoC///u/l4ceekjuuusuufHGG/fLfb/85S9LIpGQv/u7v9sv9R8QpjgawzCmBS+99FL44Q9/OFy0aFGYSqXC+vr68LTTTgtvu+22MJfLhWEYhsViMbz++uvDxYsXh8lkMpw/f354zTXXjBx/nYULF+4xJPKMM84IzzjjjJHyF77whfDkk08Om5qawmw2Gx555JHhF7/4xbBQKIycw9DJRx55JPyzP/uzcM6cOWEqlQrnzJkTvv/97w9feuklda/7778/PProo8NEIhEJo/zNb34T/vmf/3k4Y8aMMJ1OhwsXLgzf+973ho888kjkvgzRvOuuu0IRCdeuXat+/p3vfCc8/vjjw3Q6HTY3N4dnnHFG+PDDD0f64Ac/+EEoIuFf//VfR44ZxmQjIuF99903Un7ggQdCEQlra2vVv0QiEb73ve8NwzAMX3zxxVBEYv998pOf3Kv7hWEYPvPMM2FbW1u4efPmkZ8tXLgwvOWWWyb7cfc
      "text/plain": [
       "<Figure size 640x480 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-1.390557629454562e-14\n"
     ]
    }
   ],
   "source": [
    "# Now interpolate P and S, plot the fractional errors in the interpolation\n",
    "# and check for thermodynamic consistency using the derivatives\n",
    "\n",
    "#Pinterp = scipy.interpolate.RegularGridInterpolator((Tp,rhop),Pgrid, bounds_error=False, fill_value=None)\n",
    "#Sinterp = scipy.interpolate.RegularGridInterpolator((Tp,rhop),Sgrid, bounds_error=False, fill_value=None)\n",
    "\n",
    "Pinterp = scipy.interpolate.RectBivariateSpline(Tp,rhop,Pgrid)\n",
    "Sinterp = scipy.interpolate.RectBivariateSpline(Tp,rhop,Sgrid)\n",
    "\n",
    "dPdT = Pinterp.partial_derivative(1,0)\n",
    "dSdrho = Sinterp.partial_derivative(0,1)\n",
    "\n",
    "# now compute the function on a finer grid\n",
    "Tp2 = np.linspace(2,3,100)\n",
    "rhop2 = np.linspace(-6,0.0,100)\n",
    "Tgrid2, rhogrid2 = np.meshgrid(Tp2, rhop2, indexing='ij')\n",
    "ngrid2 = 10.0**rhogrid2 / (28*mu)\n",
    "Pgrid2 = np.log10(P(10**rhogrid2, 10**Tgrid2))\n",
    "Sgrid2 = S(10**rhogrid2, 10**Tgrid2)/kB\n",
    "\n",
    "# plot the fractional error\n",
    "plt.imshow((Pinterp(Tp2,rhop2)-Pgrid2)/Pgrid2, origin='lower',\n",
    "          extent=(rhop2[0],rhop2[-1],Tp2[0],Tp2[-1]), aspect='auto')\n",
    "plt.title('Pressure')\n",
    "plt.colorbar()\n",
    "plt.show()\n",
    "\n",
    "plt.clf()\n",
    "plt.imshow((Sinterp(Tp2,rhop2)-Sgrid2)/Sgrid2, origin='lower',\n",
    "           extent=(rhop2[0],rhop2[-1],Tp2[0],Tp2[-1]), aspect='auto')\n",
    "plt.title('Entropy')\n",
    "plt.colorbar()\n",
    "plt.show()\n",
    "\n",
    "plt.clf()\n",
    "f1 = 10**Pinterp(Tp2,rhop2) * dPdT(Tp2,rhop2)/ 10**Tgrid2 / ngrid2**2\n",
    "plt.imshow(f1, origin='lower',\n",
    "           extent=(rhop2[0],rhop2[-1],Tp2[0],Tp2[-1]), aspect='auto')\n",
    "plt.title('dPdT')\n",
    "plt.colorbar()\n",
    "plt.show()\n",
    "\n",
    "plt.clf()\n",
    "f2 = kB * dSdrho(Tp2,rhop2)/(np.log(10)*ngrid2)\n",
    "plt.imshow(f2, origin='lower',\n",
    "           extent=(rhop2[0],rhop2[-1],Tp2[0],Tp2[-1]), aspect='auto')\n",
    "plt.title('dSdrho')\n",
    "plt.colorbar()\n",
    "plt.show()\n",
    "\n",
    "plt.clf()\n",
    "plt.imshow(np.abs((f1+f2)/f1), origin='lower',\n",
    "           extent=(rhop2[0],rhop2[-1],Tp2[0],Tp2[-1]), aspect='auto')\n",
    "plt.title('Consistency')\n",
    "plt.colorbar()\n",
    "plt.show()\n",
    "\n",
    "print(np.min((f1+f2)/(f1-f2)))\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.11.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}