{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "22e5e8ed-fa9f-44d9-a69c-9f6c88e456f7",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy.optimize import fsolve\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "03789e89-3351-42f0-bdcd-91c3a351630d",
   "metadata": {},
   "outputs": [],
   "source": [
    "alo = {\n",
    "    'co' : np.array([3.2624510E+00,  1.5119409E-03, -3.8817550E-06,  5.5819440E-09, -2.4749510E-12, -1.4310539E+04,  4.8488970E+00]),\n",
    "    'co2': np.array([2.2757240E+00,  9.9220720E-03, -1.0409113E-05,  6.8666860E-09, -2.1172800E-12, -4.8373140E+04,  1.0188488E+01]),\n",
    "    'h2' : np.array([3.2981240E+00,  8.2494410E-04, -8.1430150E-07, -9.4754340E-11,  4.1348720E-13, -1.0125209E+03, -3.2940940E+00]),\n",
    "    'h'  : np.array([2.5000000E+00,  0.0000000E+00,  0.0000000E+00,  0.0000000E+00,  0.0000000E+00,  2.5471620E+04, -4.6011760E-01]),\n",
    "    'oh' : np.array([3.6372660E+00,  1.8509100E-04, -1.6761646E-06,  2.3872020E-09, -8.4314420E-13,  3.6067810E+03,  1.3588605E+00]),\n",
    "    'h2o': np.array([3.3868420E+00,  3.4749820E-03, -6.3546960E-06,  6.9685810E-09, -2.5065880E-12, -3.0208110E+04,  2.5902320E+00]),\n",
    "    'n2' : np.array([3.2986770E+00,  1.4082404E-03, -3.9632220E-06,  5.6415150E-09, -2.4448540E-12, -1.0208999E+03,  3.9503720E+00]),\n",
    "    'n'  : np.array([2.5030710E+00, -2.1800180E-05,  5.4205290E-08, -5.6475600E-11,  2.0999040E-14,  5.6098900E+04,  4.1675660E+00]),\n",
    "    'no' : np.array([3.3765410E+00,  1.2530634E-03, -3.3027500E-06,  5.2178100E-09, -2.4462620E-12,  9.8179610E+03,  5.8295900E+00]),\n",
    "    'o2' : np.array([3.2129360E+00,  1.1274864E-03, -5.7561500E-07,  1.3138773E-09, -8.7685540E-13, -1.0052490E+03,  6.0347370E+00]),\n",
    "    'o'  : np.array([2.9464280E+00, -1.6381665E-03,  2.4210310E-06, -1.6028431E-09,  3.8906960E-13,  2.9147640E+04,  2.9639950E+00]),\n",
    "    'ch4': np.array([7.7874150E-01,  1.7476680E-02, -2.7834090E-05,  3.0497080E-08, -1.2239307E-11, -9.8252290E+03,  1.3722195E+01]),\n",
    "}\n",
    "ahi = {\n",
    "    'co' : np.array([3.0250780E+00,  1.4426885E-03, -5.6308270E-07,  1.0185813E-10, -6.9109510E-15, -1.4268350E+04,  6.1082170E+00]),\n",
    "    'co2': np.array([4.4536230E+00,  3.1401680E-03, -1.2784105E-06,  2.3939960E-10, -1.6690333E-14, -4.8966960E+04, -9.5539590E-01]),\n",
    "    'h2' : np.array([2.9914230E+00,  7.0006440E-04, -5.6338280E-08, -9.2315780E-12,  1.5827519E-15, -8.3503400E+02, -1.3551101E+00]),\n",
    "    'h'  : np.array([2.5000000E+00,  0.0000000E+00,  0.0000000E+00,  0.0000000E+00,  0.0000000E+00,  2.5471620E+04, -4.6011760E-01]),\n",
    "    'oh' : np.array([2.8827300E+00,  1.0139743E-03, -2.2768770E-07,  2.1746830E-11, -5.1263050E-16,  3.8868880E+03,  5.5957120E+00]),\n",
    "    'h2o': np.array([2.6721450E+00,  3.0562930E-03, -8.7302600E-07,  1.2009964E-10, -6.3916180E-15, -2.9899210E+04,  6.8628170E+00]),\n",
    "    'n2' : np.array([2.9266400E+00,  1.4879768E-03, -5.6847600E-07,  1.0097038E-10, -6.7533510E-15, -9.2279770E+02,  5.9805280E+00]),\n",
    "    'n'  : np.array([2.4502680E+00,  1.0661458E-04, -7.4653370E-08,  1.8796520E-11, -1.0259839E-15,  5.6116040E+04,  4.4487580E+00]),\n",
    "    'no' : np.array([3.2454350E+00,  1.2691383E-03, -5.0158900E-07,  9.1692830E-11, -6.2754190E-15,  9.8008400E+03,  6.4172930E+00]),\n",
    "    'o2' : np.array([3.6975780E+00,  6.1351970E-04, -1.2588420E-07,  1.7752810E-11, -1.1364354E-15, -1.2339301E+03,  3.1891650E+00]),\n",
    "    'o'  : np.array([2.5420590E+00, -2.7550610E-05, -3.1028030E-09,  4.5510670E-12, -4.3680510E-16,  2.9230800E+04,  4.9203080E+00]),\n",
    "    'ch4': np.array([1.6834780E+00,  1.0237236E-02, -3.8751280E-06,  6.7855850E-10, -4.5034230E-14, -1.0080787E+04,  9.6233950E+00]),    \n",
    "}\n",
    "\n",
    "Rg = 8314.46      # J/kmol*K\n",
    "\n",
    "def cp(sp, T):\n",
    "    a = alo[sp] if T <= 1000 else ahi[sp]\n",
    "    return Rg*( a[0] + T*(a[1] + T*(a[2] + T*(a[3] + T*(a[4])))) )\n",
    "\n",
    "def h(sp,T):\n",
    "    a = alo[sp] if T <= 1000 else ahi[sp]\n",
    "    return Rg*( a[5] + T*(a[0] + T*(a[1]/2 + T*(a[2]/3 + T*(a[3]/4 + T*(a[4]/5))))) )\n",
    "\n",
    "def s(sp,T):\n",
    "    a = alo[sp] if T <= 1000 else ahi[sp]\n",
    "    return Rg*( a[0]*np.log(T) + a[6] + T*(a[1] + T*(a[2]/2 + T*(a[3]/3 + T*(a[4]/4)))) )\n",
    "    \n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8caef12d-4f40-495b-9c5d-1e7cf2f20bd7",
   "metadata": {},
   "source": [
    "### Compute $T_{ad}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "991161fc-03ac-44a2-938e-5b5e2de68fb7",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Tp = 2328.51 K\n"
     ]
    }
   ],
   "source": [
    "Tr = 300.0\n",
    "hr = h('ch4', Tr) + 2*h('o2', Tr) + 7.52*h('n2', Tr)\n",
    "Tp = fsolve(lambda T: h('co2',T) + 2*h('h2o',T) + 7.52*h('n2',T) - hr, Tr)[0]\n",
    "print(f\"Tp = {Tp:.2f} K\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "9e7f589e-be48-42b0-b7ae-8e7c1447e0a6",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(2328.5096476447256, -5.960464477539063e-08)"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def Ftozero(T):\n",
    "    Tr = 300.0\n",
    "    hr = h('ch4', Tr) + 2*h('o2', Tr) + 7.52*h('n2', Tr)\n",
    "    hp = h('co2',T) + 2*h('h2o',T) + 7.52*h('n2',T)\n",
    "    return hr - hp\n",
    "    \n",
    "Tp = fsolve(Ftozero, Tr)[0]\n",
    "Tp, Ftozero(Tp)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "28c02b13-e583-4839-ba52-a0d463e1563c",
   "metadata": {},
   "outputs": [],
   "source": [
    "T = np.linspace(300,5000,1000)\n",
    "cc = np.zeros(len(T))\n",
    "for i in range(len(T)):\n",
    "    cc[i] = cp('co2',T[i])\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "af93d30c-3763-450d-9e56-b5974d7556a3",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x122095520>]"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYMAAAD6CAYAAABDPiuvAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAj8UlEQVR4nO3deXxV9Z3/8dcne8KWAAFJAiQIolgRMEW01nWquFSt48zg1IrWKe2onc7WVjvdpsv8pu1MbTszdRm1aqetWqdWxrFFqmIdq2DYN4GwJmFJQhYCScj2+f1xv8RL2AIkuTe57+fjcR/33M895/I9R3Pe537P95xr7o6IiCS2pFg3QEREYk9hICIiCgMREVEYiIgICgMREUFhICIidDMMzCzbzJ43s/fMbL2ZXWRmXzezCjNbER7XRc3/gJmVmtkGM7smqj471ErN7P6oepGZLQ71Z80srWdXU0REjse6c52BmT0FvOnuj4UddRbw18B+d/+XLvNOAX4BzATygN8BZ4W3NwIfAcqBd4Hb3H2dmT0H/MrdnzGzh4GV7v7Q8do0cuRILyws7PaKiogILF26tNrdc7vWU060oJkNAy4F7gRw9xagxcyOtchNwDPufhDYamalRIIBoNTdt4TPfQa4yczWA1cCfx7meQr4OnDcMCgsLKSkpOREzRcRkShmtv1o9e50ExUBVcBPzGy5mT1mZoPCe/eZ2Soze8LMckItHyiLWr481I5VHwHUuXtbl/rRVmKemZWYWUlVVVU3mi4iIt3RnTBIAWYAD7n7dOAAcD+RI/czgWnALuBfe6mNndz9UXcvdvfi3NwjvuWIiMgp6k4YlAPl7r44vH4emOHue9y93d07gP/k/a6gCmBs1PIFoXas+l4g28xSutRFRKSPnDAM3H03UGZmk0PpKmCdmY2Jmu1jwJowPR+YY2bpZlYETAKWEDlhPCmMHEoD5gDzPXIG+3Xg1rD8XODF01wvERE5CSc8gRx8FvhZ2IlvAe4CfmRm0wAHtgGfBnD3tWF00DqgDbjX3dsBzOw+YAGQDDzh7mvD538ReMbMvgUsBx4//VUTEZHu6tbQ0nhUXFzsGk0kInJyzGypuxd3resKZBER6XY3kYiInAR3p7GlnX3NrTS2tNPU0k5Ta3uYbqOxJTLd0tZBhzvtHU67Ox0dTodDe0ek1yY12UhOSiI12UhJMlKSk/izD44lNblnj+UVBiIi3dDc2k5Vw0Gq9h+MPDccpOZAC3WNrdQ1tVDf2EpdUyt1jS3UN7VR39RCa3vvdMPfekEBqck9+5kKAxFJaO0dzp59zVTUNVFR28Tufc1UNRyksuEgVQ3NnTv+fc1tR11+cHoKwzJTGZaZSnZWKpPPGMKwzDSys1LJzkxlSEYqWWnJZKYlkxUemakpna/TU5JISjKSzUhOMpLMSDJITjLcod2dtnantaODtnanraOD9JSe7+FXGIjIgNbc2s7OuqbOnf1hz3VN7K5vpq3j8CP4zNRkRg1NZ9SQdCafMYRLJo4kd0g6uUPSGTUko3M6JyuNtF7YMR9iBkkYqcmQSQ9/FehCYSAi/Vp9U2tkZx+1g6+obaI8PFfvP3jY/EkGZwzNID8nkwvG55CfnUl+Tib52ZkU5GRyxrBMBqcn3q4x8dZYRPqNjg6n+sDBI4/oo54bDh7efZOWkhTZwWdnctXZozp39IeezxiW0eMnXwcChYGIxExrewe765spP2wn38jOuubOo/yWto7DlhmSkdJ5FH9h0fCwk88iPyeTvOwMRg5KJynpmHdVlmNQGIhIr2lsaTui+yb6ec++Zrp015M7JJ387Eym5A3l6imjycuOOrLPyWRoRmpsVmaAUxiIyClxd+oaW6moazriyP7QdG1j62HLpCQZZwzLID87k4vOHEFBZ/dN5Mh+zLAMMnp6zKR0i8JARI6qvcOpbGjuPIqP3uEfGp3T2NJ+2DKZqcmdffNTC7I7u3MOHdmPGpJBsrpw4pLCQCRBHWxrZ2ddc+dInPIuR/a76o4ccpmTlUp+TiYTcgfx4Um55GVnhJ195Mg+JyuV4/wKosQxhYHIANXQ3HpEP3151OuqhsOHXFoYcpmXncn0sTncMPX9I/qC7EzysjMZlIBDLhOF/suK9EPuzt4DLZGum0NH87VNVBwahVPbeMQVs2nJSeRlR8bXXzE5t/No/v3x9RpymcgUBiJxqr6xlbLaRsprGymraQrTTZTVRJ6bWg/vrx+SntK5c/9gYeRiqryoI/uRgzXkUo5NYSASIwfb2imraWT73kbKahopq206bMff0OXIfkh6CgXDsygaOYhLz8qlICeTgpyszq6cYZkacimnTmEg0ota2zsor21iW/UBtlQfYFv1AbbtPcDW6gPsrGs6bIx9RmoSBTlZjM3JpLgwh7E5WRTkZDJ2eBZjc7IYlqWdvfQehYFID9jX3Epp5X427Wlg0579lFbtZ1v1AcpqmzrvSw+Rq2eLRg5ixrgcbplRQNHILMaPGMTYnCxGDk7TSByJGYWByEnYf7CNDbsbKK1sYOOe/WwKAbCrvrlznozUJCaMHMy5ecO4fuoYCkcMYkLuIApHDGL4IO3wJT4pDESOYe/+g6zduY+1O/exZmc963buY2v1gc73M1KTmDhqMLMmjGDS6MGcNWoIk0YPpiAnSxdWSb+jMBAhMiZ/ZVk9y3bUsqq8jrU79x12tF+Qk8m5eUO5ZXo+Z48Zylna6csAozCQhNPR4Wyp3s+y7XUsL6tl2fY6NlY24B658GrCyEFcWDScc/OGcW7eUKbkDSU7Ky3WzRbpVQoDGfDcnY179vP25mre3rKXxVtrqAs3UBuakcL0cTlcd94YZozP5vyx2borpiQkhYEMSNuqD/BmaTXvbN7LO1v2svdACwBjh2dy9ZTRFBcOZ8a4HCaMHKQLsURQGMgA0dzazpKtNby+oZJFG6o6T/SOGZbBZWflctGZIyK3TM7JinFLReKTwkD6rer9B3ll7R5ee28Pb5Xupam1nfSUJC46cwR3XlzIpWflUjgiS0M5RbpBYSD9yp59zSxYu5uXV+9iydYaOjwy0udPigu4YvIoZk0YQWaafhxF5GQpDCTu1Rxo4aVVO5m/YidLd9TiDhNHDea+KyYy+wNjOGfMEB39i5ymboWBmWUDjwEfABz4JLABeBYoBLYBf+rutRb5q/whcB3QCNzp7svC58wFvhw+9lvu/lSoXwA8CWQCLwOfc/cuv4wqieRgWzuvv1fJfy+r4PX3KmnrcM4+Ywh/80dnce0HzmDS6CGxbqLIgNLdbwY/BH7r7reaWRqQBXwJeNXd/9nM7gfuB74IXAtMCo8LgYeAC81sOPA1oJhIoCw1s/nuXhvm+RSwmEgYzAZ+00PrKP3I+l37+PniHcxfuZP6plZyh6Rz14cK+dj0AqbkDY1180QGrBOGgZkNAy4F7gRw9xagxcxuAi4Psz0FLCISBjcBT4cj+3fMLNvMxoR5F7p7TfjchcBsM1sEDHX3d0L9aeBmFAYJo7m1nd+s2cV/vbODpdtrSUtJYva5Z3DLjHwumTiSFP3gikiv6843gyKgCviJmZ0PLAU+B4x2911hnt3A6DCdD5RFLV8easerlx+lfgQzmwfMAxg3blw3mi7xrLy2kZ++vZ3nSsqobWylcEQW/3DdOdx6QQE5g3TFr0hf6k4YpAAzgM+6+2Iz+yGRLqFO7u5m1ut9/O7+KPAoQHFxsc4p9FNrKup55PdbeHl15FjiI+eM5vZZ47n4zBG6AEwkRroTBuVAubsvDq+fJxIGe8xsjLvvCt1AleH9CmBs1PIFoVbB+91Kh+qLQr3gKPPLAOLuvLGxikd/v4U/bN7L4PQU7r6kiDsvLiQvOzPWzRNJeCcMA3ffbWZlZjbZ3TcAVwHrwmMu8M/h+cWwyHzgPjN7hsgJ5PoQGAuAfzKznDDf1cAD7l5jZvvMbBaRE8h3AP/Wg+soMeTuLFy3hx/8bhPrdu1j9NB0Hrj2bG67cJzuASQSR7o7muizwM/CSKItwF1AEvCcmd0NbAf+NMz7MpFhpaVEhpbeBRB2+t8E3g3zfePQyWTgHt4fWvobdPK433N3Xnuvkh/8bhOrK+oZPyKL7946lZun5ZOWohPCIvHG+utw/uLiYi8pKYl1M+Qo3txUxb8s2MDK8nrGDs/kr66cxMem52tUkEgcMLOl7l7cta4rkKXHbNjdwD+9vJ43NlaRn53Jd/74PG6ZUUCqQkAk7ikM5LRVNjTz4MKNPPtuGYPTU/jy9efwiYvGk56iewSJ9BcKAzllLW0dPP5/W/m31zbR2t7BnRcX8VdXTdSvgon0QwoDOSVvb97LV15cQ2nlfq6eMpovXXcOhSMHxbpZInKKFAZyUiobmvmn/13Pr1fsZNzwLH5y5we54uxRsW6WiJwmhYF0i7vz/NJyvvHSOg62dvBXV07knismkpGq8wIiA4HCQE5oV30TD/xqNYs2VDGzcDj//MfnMSF3cKybJSI9SGEgx+Tu/LKknG++tI62DufrH53CHRcV6v5BIgOQwkCOau/+g3zh+VW8+l4lM4uG871bpzJ+hE4QiwxUCgM5wlul1fzNsyuoa2zlqzdM4c6L9W1AZKBTGEin1vYO/vWVjTzy+81MGDmIJ++aqV8XE0kQCgMBIj80c+/Pl7OyrI7bZo7jqzdMITNNI4VEEoXCQPi/TdV89hfLaGt3fvzxGVx33phYN0lE+pjCIIG5Ow+9sZl/WbCBiaMG88gniinSVcQiCUlhkKAamlv5/C9X8du1u7lh6hi+88dTGZSu/x1EEpX++hNQWU0jdz/1LpurDvDl68/h7kuKMNNoIZFEpjBIMMt31PKpp0toaevg6U/O5EMTR8a6SSISBxQGCeSlVTv5u+dWMnpoBs/M+yATR+mWEiISoTBIAO7Ojxdt5nsLNlA8PodH7yhm+CD95oCIvE9hMMB1dDjfeGkdT/5hGzeen8d3b52qO42KyBEUBgNYa3sHf//Llby4Yid/cUkRX7ruHN1WQkSOSmEwQDW2tHHPz5axaEMVX5g9mb+87EyNGBKRY1IYDED1Ta188sl3Wb6jlv93y3ncNnNcrJskInFOYTDA1De28oknFrN+1z7+/c91awkR6R6FwQBS19jC7Y8vZuPu/Tx8+wVcdc7oWDdJRPoJhcEAUdfYwscfW8ymPft55BMX6EfqReSkKAwGgNoDkSAordrPo3dcwOWTFQQicnIUBv1cQ3MrdzyxhNKq/fznHcVcdlZurJskIv1QUndmMrNtZrbazFaYWUmofd3MKkJthZldFzX/A2ZWamYbzOyaqPrsUCs1s/uj6kVmtjjUnzUzXR7bDc2t7dz9VAnrd+3j4dtnKAhE5JR1KwyCK9x9mrsXR9UeDLVp7v4ygJlNAeYA5wKzgR+bWbKZJQP/AVwLTAFuC/MCfCd81kSgFrj79FZr4Gtp6+Av/2sp726r4ft/No0rz9bJYhE5dScTBt11E/CMux90961AKTAzPErdfYu7twDPADdZ5EqoK4Hnw/JPATf3QrsGjPYO52+fW8HrG6r49s3nceP5ebFukoj0c90NAwdeMbOlZjYvqn6fma0ysyfMLCfU8oGyqHnKQ+1Y9RFAnbu3dakfwczmmVmJmZVUVVV1s+kDi7vz1RfX8NKqXdx/7dn8+YW6oExETl93w+ASd59BpIvnXjO7FHgIOBOYBuwC/rVXWhjF3R9192J3L87NTcz+8R8v2szPFu/g05dN4DOXnRnr5ojIANGtMHD3ivBcCbwAzHT3Pe7e7u4dwH8S6QYCqADGRi1eEGrHqu8Fss0spUtduvj18gq+t2ADN03L44vXnB3r5ojIAHLCMDCzQWY25NA0cDWwxsyi73PwMWBNmJ4PzDGzdDMrAiYBS4B3gUlh5FAakZPM893dgdeBW8Pyc4EXT3/VBpY/bK7m88+vZNaE4Xz31qm6+6iI9KjuXGcwGngh3PEyBfi5u//WzH5qZtOInE/YBnwawN3XmtlzwDqgDbjX3dsBzOw+YAGQDDzh7mvDv/FF4Bkz+xawHHi8Z1ZvYNi4p4FP/3QphSMG8cjtxaSn6PcIRKRnWeTAvP8pLi72kpKSWDej19U3tnLDv79Jc2sHL9xzMQU5WbFukoj0Y2a2tMslAkDvDC2VHtLR4fz1s8vZXd/MI5+4QEEgIr1GYRDH/u21Ul7fUMVXP3ouM8blnHgBEZFTpDCIU3/YXM0PXt3ILdPzuV3XEohIL1MYxKF9za18/perKBwxiG997AP6uUoR6XW6a2kc+sf569i9r5nnP3MRWWn6TyQivU/fDOLMgrW7+e9l5dx7+ZlM13kCEekjCoM4sq+5la/8eg3n5g3ls1dNinVzRCSBqA8ijnz/lY1U7T/IY3OLSU1WTotI39EeJ06sKq/jqbe3cces8UwtyI51c0QkwSgM4kB7h/OlF1aTOzidv7tmcqybIyIJSGEQB36+ZAdrKvbxlRumMDQjNdbNEZEEpDCIsX3NrTy4cCMXFg3nhqljTryAiEgvUBjE2I9f30zNgRa+fP0UXVwmIjGjMIihsppGnnhrK7dMz+e8gmGxbo6IJDCFQQx9b8EGDPh7nTQWkRhTGMTImop65q/cyac+PIG87MxYN0dEEpzCIEa+v3AjwzJTmXfZhFg3RUREYRALy3bU8tp7lcy7dIKGkopIXFAYxMCDCzcyYlAad15cGOumiIgACoM+t2RrDW9uquYzl53JoHTdGkpE4oPCoI99f+EGcoekc/us8bFuiohIJ4VBH1q6vZZ3ttTw6UsnkJmWHOvmiIh0Uhj0oYff2Ex2Viq3zdRvGotIfFEY9JGNexpYuG4Pcy8q1LkCEYk7CoM+8vAbm8lMTWauRhCJSBxSGPSB8tpG5q/YyZyZYxk+KC3WzREROYLCoA889uZWAP7iw7raWETik8Kgl9U3tvLsu2XcOC2PfN2DSETiVLfCwMy2mdlqM1thZiWhNtzMFprZpvCcE+pmZj8ys1IzW2VmM6I+Z26Yf5OZzY2qXxA+vzQsO2Bu7P9cSRlNre3cfUlRrJsiInJMJ/PN4Ap3n+buxeH1/cCr7j4JeDW8BrgWmBQe84CHIBIewNeAC4GZwNcOBUiY51NRy80+5TWKI+0dzlNvb2Nm0XDOzdPvFYhI/DqdbqKbgKfC9FPAzVH1pz3iHSDbzMYA1wAL3b3G3WuBhcDs8N5Qd3/H3R14Ouqz+rXfrd9DeW0Td2kEkYjEue6GgQOvmNlSM5sXaqPdfVeY3g2MDtP5QFnUsuWhdrx6+VHq/d6Tb20jPzuTj0wZfeKZRURiqLtXP13i7hVmNgpYaGbvRb/p7m5m3vPNO1wIonkA48bF91W863ft4+0te7n/2rNJSdZ5ehGJb93aS7l7RXiuBF4g0ue/J3TxEJ4rw+wVwNioxQtC7Xj1gqPUj9aOR9292N2Lc3Nzu9P0mHnqD9vISE1izgfHnnhmEZEYO2EYmNkgMxtyaBq4GlgDzAcOjQiaC7wYpucDd4RRRbOA+tCdtAC42sxywonjq4EF4b19ZjYrjCK6I+qz+qX6plZ+vaKCj03PJztLF5mJSPzrTjfRaOCFMNozBfi5u//WzN4FnjOzu4HtwJ+G+V8GrgNKgUbgLgB3rzGzbwLvhvm+4e41Yfoe4EkgE/hNePRbv15eQXNrBx+/ULepFpH+wSIDePqf4uJiLykpiXUzjuDuXPvDN0lNTuJ/PntJrJsjInIYM1sadYlAJ53Z7GEryup4b3cDc2bqXIGI9B8Kgx72iyU7yEpL5sbz82LdFBGRblMY9KCG5lb+Z+Uubjw/jyEZqbFujohItykMetCLK3bS1NrOHP2SmYj0MwqDHvSLJTs4Z8xQzi/QfYhEpH9RGPSQ9bv2sXbnPv6suIABdNNVEUkQCoMe8sLyClKSjBunDYjbKolIglEY9IC29g5eWF7BFWeP0s9aiki/pDDoAW9t3ktVw0Fuma5vBSLSPykMesALy8oZmpHCleeMinVTREROicLgNO0/2MZv1+7mo+fnkZ6SHOvmiIicEoXBafrtmt00t3Zwy4yCE88sIhKnFAan6VfLyikckcWMcdmxboqIyClTGJyG3fXNvL1lLzdPz9e1BSLSrykMTsPLq3fhDh/VTelEpJ9TGJyG/129i3PGDOXM3MGxboqIyGlRGJyinXVNLN1eyw1Tx8S6KSIip01hcIpeXr0LgOvPUxiISP+nMDhFL63axQfyh1I4clCsmyIictoUBqegrKaRFWV13DBVJ45FZGBQGJwCdRGJyECjMDgFL63axfkFwxg7PCvWTRER6REKg5NUVtPI6op6rtcoIhEZQBQGJ+mVdXsAuObcM2LcEhGRnqMwOEkL1+1m8ughjB+hUUQiMnAoDE5C7YEW3t1Wy0emjI51U0REepTC4CS89l4l7R3O1ecqDERkYFEYnISF6/ZwxtAMzssfFuumiIj0KIVBNzW3tvPGxir+aMoo3a5aRAacboeBmSWb2XIzeym8ftLMtprZivCYFupmZj8ys1IzW2VmM6I+Y66ZbQqPuVH1C8xsdVjmRxaHe9u3Sqtpam3n6ikaRSQiA8/JfDP4HLC+S+3z7j4tPFaE2rXApPCYBzwEYGbDga8BFwIzga+ZWU5Y5iHgU1HLzT75Veldr6zdw5D0FGZNGBHrpoiI9LhuhYGZFQDXA491Y/abgKc94h0g28zGANcAC929xt1rgYXA7PDeUHd/x90deBq4+RTWpde0dzivvreHyybnkpainjURGXi6u2f7AfAFoKNL/duhK+hBM0sPtXygLGqe8lA7Xr38KPUjmNk8Mysxs5KqqqpuNv30rSyvo3p/i4aUisiAdcIwMLMbgEp3X9rlrQeAs4EPAsOBL/Z88w7n7o+6e7G7F+fm5vb2P9dp0YYqkgwuO6vv/k0Rkb7UnW8GHwJuNLNtwDPAlWb2X+6+K3QFHQR+QuQ8AEAFMDZq+YJQO1694Cj1uPHGhkqmjc0mOyst1k0REekVJwwDd3/A3QvcvRCYA7zm7reHvn7CyJ+bgTVhkfnAHWFU0Syg3t13AQuAq80sJ5w4vhpYEN7bZ2azwmfdAbzYs6t56qr3H2RleT1XTB4V66aIiPSalNNY9mdmlgsYsAL4TKi/DFwHlAKNwF0A7l5jZt8E3g3zfcPda8L0PcCTQCbwm/CIC7/fGDk3cbnCQEQGsJMKA3dfBCwK01ceYx4H7j3Ge08ATxylXgJ84GTa0lcWbahi5OA0zs0bGuumiIj0Go2TPI72Duf3m6q47KxRJCXF3XVwIiI9RmFwHCvK6qhrbOXyyRpFJCIDm8LgON7YUEmSwYcnjYx1U0REepXC4DgWbaxi+rgcDSkVkQFPYXAMtQdaWF1Rz6WT1EUkIgOfwuAY3t6yF3e4ZJJuTCciA5/C4BjeKq1mcHoKUwuyY90UEZFepzA4hrdKq5k1YTipydpEIjLwaU93FGU1jWzb28iHJmoUkYgkBoXBUfxhczWAwkBEEobC4CjeKt1L7pB0Jo0aHOumiIj0CYVBFx0dzlul1VwycaR++F5EEobCoIsNexrYe6CFi8/UkFIRSRwKgy7eKtX5AhFJPAqDLt7evJeikYPIy86MdVNERPqMwiBKe4ezZFsNsyYMj3VTRET6lMIgynu799HQ3MbMIoWBiCQWhUGUxVsiv8J5YZFOHotIYlEYRFmytYaCnEydLxCRhKMwCNwj5wv0rUBEEpHCICit3E/NgRYu1PkCEUlACoPgna3hfIFGEolIAlIYBEu21jB6aDrjhmfFuikiIn1OYUDkfMHiLXu5sGiE7kckIglJYQBs39tIZcNBXV8gIglLYQAs2Xbo+gKFgYgkJoUBsGx7LcMyUzkzV79fICKJSWEALN1ey/Rx2SQl6XyBiCSmboeBmSWb2XIzeym8LjKzxWZWambPmllaqKeH16Xh/cKoz3gg1DeY2TVR9dmhVmpm9/fg+p1QfVMrmyr3c8G4nL78Z0VE4srJfDP4HLA+6vV3gAfdfSJQC9wd6ncDtaH+YJgPM5sCzAHOBWYDPw4Bkwz8B3AtMAW4LczbJ5bvqAVgxniFgYgkrm6FgZkVANcDj4XXBlwJPB9meQq4OUzfFF4T3r8qzH8T8Iy7H3T3rUApMDM8St19i7u3AM+EefvEsh11JBmcPza7r/5JEZG4091vBj8AvgB0hNcjgDp3bwuvy4H8MJ0PlAGE9+vD/J31Lsscq34EM5tnZiVmVlJVVdXNph/fsu21nH3GUAanp/TI54mI9EcnDAMzuwGodPelfdCe43L3R9292N2Lc3NzT/vz2juc5TtqmTE++/QbJyLSj3XncPhDwI1mdh2QAQwFfghkm1lKOPovACrC/BXAWKDczFKAYcDeqPoh0cscq96rNu5p4EBLOxfofIGIJLgTfjNw9wfcvcDdC4mcAH7N3T8OvA7cGmabC7wYpueH14T3X3N3D/U5YbRRETAJWAK8C0wKo5PSwr8xv0fW7gSWbg8njzWSSEQS3Ol0lH8ReMbMvgUsBx4P9ceBn5pZKVBDZOeOu681s+eAdUAbcK+7twOY2X3AAiAZeMLd155Gu7pt2Y5aRg5O083pRCThnVQYuPsiYFGY3kJkJFDXeZqBPznG8t8Gvn2U+svAyyfTlp6wbHst08fl6OZ0IpLwEvYK5L37D7Jtb6O6iERESOAwWFVRD8A0XV8gIpLAYVBWjxmcVzAs1k0REYm5xA2D8jrOzB2si81EREjQMHB3VpbXM1XfCkREgAQNg131zVTvP8j5BdmxboqISFxIyDBYVV4HoG8GIiJBQobByvJ6UpKMc8YMjXVTRETiQkKGwaryOiafMYSM1ORYN0VEJC4kXBh0dDiryuuZqvMFIiKdEi4Mttc00tDcxvk6XyAi0inhwuD9k8fZMW2HiEg8SbgwWFlWT0ZqEmeNHhzrpoiIxI2EC4NV5XWcmzeMlOSEW3URkWNKuHsxnD82mzHDMmLdDBGRuJJwYfCVG6bEugkiInFHfSUiIqIwEBERhYGIiKAwEBERFAYiIoLCQEREUBiIiAgKAxERAczdY92GU2JmVcD2WLejD40EqmPdiBjTNtA2AG2D013/8e6e27XYb8Mg0ZhZibsXx7odsaRtoG0A2ga9tf7qJhIREYWBiIgoDPqTR2PdgDigbaBtANoGvbL+OmcgIiL6ZiAiIgoDERFBYRBTZvaEmVWa2Zqo2nAzW2hmm8JzTqibmf3IzErNbJWZzYhaZm6Yf5OZzY3FupwKMxtrZq+b2TozW2tmnwv1RNoGGWa2xMxWhm3wj6FeZGaLw7o+a2ZpoZ4eXpeG9wujPuuBUN9gZtfEaJVOiZklm9lyM3spvE6o9Qcws21mttrMVphZSaj13d+Cu+sRowdwKTADWBNV+y5wf5i+H/hOmL4O+A1gwCxgcagPB7aE55wwnRPrdevm+o8BZoTpIcBGYEqCbQMDBofpVGBxWLfngDmh/jDwl2H6HuDhMD0HeDZMTwFWAulAEbAZSI71+p3Edvhb4OfAS+F1Qq1/WIdtwMgutT77W4j5Bkj0B1DYJQw2AGPC9BhgQ5h+BLit63zAbcAjUfXD5utPD+BF4COJug2ALGAZcCGRK0xTQv0iYEGYXgBcFKZTwnwGPAA8EPVZnfPF+wMoAF4FrgReCuuTMOsf1eajhUGf/S2omyj+jHb3XWF6NzA6TOcDZVHzlYfaser9Svi6P53IkXFCbYPQRbICqAQWEjmqrXP3tjBL9Pp0rmt4vx4YQf/eBj8AvgB0hNcjSKz1P8SBV8xsqZnNC7U++1tIOdVWS+9zdzezAT/218wGA/8N/LW77zOzzvcSYRu4ezswzcyygReAs2Pbor5jZjcAle6+1Mwuj3FzYu0Sd68ws1HAQjN7L/rN3v5b0DeD+LPHzMYAhOfKUK8AxkbNVxBqx6r3C2aWSiQIfubuvwrlhNoGh7h7HfA6kW6RbDM7dLAWvT6d6xreHwbspf9ugw8BN5rZNuAZIl1FPyRx1r+Tu1eE50oiBwUz6cO/BYVB/JkPHBoBMJdIP/qh+h1hFMEsoD58fVwAXG1mOWGkwdWhFvcs8hXgcWC9u38/6q1E2ga54RsBZpZJ5JzJeiKhcGuYres2OLRtbgVe80jn8HxgThhtUwRMApb0yUqcBnd/wN0L3L2QyAnh19z94yTI+h9iZoPMbMihaSL/D6+hL/8WYn3SJJEfwC+AXUArkb69u4n0f74KbAJ+BwwP8xrwH0T6k1cDxVGf80mgNDzuivV6ncT6X0Kkn3QVsCI8rkuwbTAVWB62wRrgq6E+gcjOrBT4JZAe6hnhdWl4f0LUZ/1D2DYbgGtjvW6nsC0u5/3RRAm1/mF9V4bHWuAfQr3P/hZ0OwoREVE3kYiIKAxERASFgYiIoDAQEREUBiIigsJARERQGIiICPD/AVu7PsMlvMsvAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(T,cc)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5883d1b0-b45d-478b-b869-31977feb35ec",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9a295310-4ace-408c-9890-c2b2f2ff8b5d",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}