{ "cells": [ { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "# Import the numpy module to provide numerical functionality\n", "import numpy as np\n", "\n", "# Import the matplotlib.pyplot module to provide plotting functionality\n", "import matplotlib.pyplot as plt\n", "\n", "# Tell matplotlib.pyplot to do inline plots\n", "%matplotlib inline\n", "\n", "# Import the mesa-web module to simplify reading MESA-Web files\n", "\n", "import mesa_web as mw" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", "dict_keys(['model_number', 'num_zones', 'initial_mass', 'initial_z', 'star_age', 'time_step', 'Teff', 'photosphere_L', 'photosphere_r', 'center_eta', 'center_h1', 'center_he3', 'center_he4', 'center_c12', 'center_n14', 'center_o16', 'center_ne20', 'star_mass', 'star_mdot', 'star_mass_h1', 'star_mass_he3', 'star_mass_he4', 'star_mass_c12', 'star_mass_n14', 'star_mass_o16', 'star_mass_ne20', 'he_core_mass', 'c_core_mass', 'o_core_mass', 'si_core_mass', 'fe_core_mass', 'neutron_rich_core_mass', 'tau10_mass', 'tau10_radius', 'tau100_mass', 'tau100_radius', 'dynamic_time', 'kh_timescale', 'nuc_timescale', 'power_nuc_burn', 'power_h_burn', 'power_he_burn', 'power_neu', 'burn_min1', 'burn_min2', 'time_seconds', 'version_number', 'compiler', 'build', 'MESA_SDK_version', 'date', 'mass', 'radius', 'luminosity', 'pressure', 'logRho', 'logT', 'energy', 'entropy', 'cp', 'gamma1', 'grada', 'mu', 'free_e', 'ye', 'pgas', 'prad', 'gradr', 'gradT', 'velocity', 'conv_vel', 'opacity', 'eps_nuc', 'pp', 'cno', 'tri_alfa', 'eps_nuc_neu_total', 'non_nuc_neu', 'eps_grav', 'h1', 'he3', 'he4', 'c12', 'n14', 'o16', 'ne20', 'mg24', 'si28', 's32', 'ar36', 'ca40', 'ti44', 'cr48', 'fe52', 'fe54', 'fe56', 'ni56', 'eta', 'log_omega', 'v_rot', 'j_rot', 'dynamo_log_B_r', 'dynamo_log_B_phi', 'log_D_conv', 'log_D_semi', 'log_D_ovr', 'log_D_thrm'])\n" ] } ], "source": [ "# Read history data. Be sure to replace the MMDDNNNNNN with the\n", "# specific digits of your folder, and P with the profile index number\n", "\n", "prof_data = mw.read_profile('MESA-Web_Job_1008205501/profile8.data')\n", "\n", "# Inspect the prof_data variable\n", "\n", "print(type(prof_data))\n", "print(prof_data.keys())" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "scrolled": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Help on function read_profile in module mesa_web:\n", "\n", "read_profile(filename)\n", " Read data from a MESA-Web profile file\n", " \n", " Parameters\n", " ----------\n", " \n", " filename : string giving name of profile file\n", " \n", " Returns\n", " -------\n", " \n", " prof_data: dict containing header and profile data (see below for\n", " details)\n", " \n", " Header Data\n", " -----------\n", " \n", " The following keys/value pairs in the data dict contain header\n", " data -- i.e., scalars describing position-independent properties\n", " of the star. Where applicable, units are given in square brackets\n", " [].\n", " \n", " star_mdot -- mass-loss rate [Msun/year]\n", " model_number -- model number\n", " num_zones -- number of zones\n", " initial_mass -- initial mass [Msun]\n", " initial_z -- initial metal mass fraction\n", " star_age -- stellar age [years]\n", " time_step -- current time-step [s]\n", " Teff -- effective temperature [K]\n", " photosphere_L -- photospheric luminosity [Lsun]\n", " photosphere_r -- photospheric radius [Rsun]\n", " center_eta -- center electron chemical potential [kB*T]\n", " center_h1 -- center 1H mass fraction\n", " center_he3 -- center 3He mass fraction\n", " center_he4 -- center 4He mass fraction\n", " center_c12 -- center 12C mass fraction\n", " center_n14 -- center 14N mass fraction\n", " center_o16 -- center 16O mass fraction\n", " center_ne20 -- center 20Ne mass fraction\n", " star_age -- stellar age [years]\n", " star_mass -- stellar mass [Msun]\n", " star_mdot -- mass-loss rate [Msun/year]\n", " star_mass_h1 -- total mass in 1H [Msun]\n", " star_mass_he3 -- total mass in 3He [Msun]\n", " star_mass_he4 -- total mass in 4He [Msun]\n", " star_mass_c12 -- total mass in 12C [Msun]\n", " star_mass_n14 -- total mass in 14N [Msun]\n", " star_mass_o16 -- total mass in 16O [Msun]\n", " star_mass_ne20 -- total mass in 20Ne [Msun]\n", " he_core_mass -- mass of helium core [Msun]\n", " c_core_mass -- mass of carbon core [Msun]\n", " o_core_mass -- mass of oxygen core [Msun]\n", " si_core_mass -- mass of silicon core [Msun]\n", " fe_core_mass -- mass of iron core [Msun]\n", " neutron_rich_core_mass -- mass of neutron-rich core [Msun]\n", " tau10_mass -- mass coordinate of optical depth 10 [Msun]\n", " tau10_radius -- radius coordinate of optical depth 10 [Rsun]\n", " tau100_mass -- mass coordinate of optical depth 100 [Msun]\n", " tau100_radius -- radius coordinate of optical depth 100 [Rsun]\n", " dynamic_timescale -- dynamical timescale [s]\n", " kh_timescale -- Kelvin-Helmholtz timescale [s]\n", " nuc_timescale -- nuclear timescale [s]\n", " log_LH -- log10(total H-burning luminosity, excluding neutrinos [Lsun])\n", " log_LHe -- log10(total He-burning luminosity, excluding neutrinos [Lsun])\n", " log_LZ -- log10(total metal-burning luminosity, excluding neutrinos [Lsun])\n", " power_nuc_burn -- total nuclear burning luminosity, excluding photodisintegrations [Lsun]\n", " power_h_burn -- total H-burning luminosity, excluding neutrinos [Lsun]\n", " power_he_burn -- total He-burning luminosity, excluding neutrinos [Lsun]\n", " power_neu -- total neutrino luminosity [Lsun]\n", " burn_min1 -- 1st limit for reported burning [erg/g/s]\n", " burn_min2 -- 2nd limit for reported burning [erg/g/s]\n", " \n", " Profile Data\n", " ------------\n", " \n", " The following keys/value pairs in the data dict contain profile\n", " data -- i.e., describing local properties of the star over a\n", " sequence of spatial zones. Where applicable, units are given in\n", " square brackets [].\n", " \n", " mass -- mass coordinate [Msun]\n", " radius -- radius coordinate [Rsun]\n", " luminosity -- luminosity [Lsun]\n", " pressure -- pressure [dyn/cm^2]\n", " logRho -- log10(density [g/cm^3])\n", " logT -- log10(temperature [K])\n", " energy -- log10(specific internal energy [erg/g])\n", " entropy -- log10(specific entropy [kB*N_avo])\n", " cp -- specific heat at constant pressure [erg/K/g]\n", " gamma1 -- first adiabatic exponent\n", " grada -- adiabatic temperature gradient\n", " mu -- mean molecular weight\n", " free_e -- free electrons per nucleon\n", " ye -- average charge per baryon [e]\n", " pgas -- gas pressure [dyn/cm^2]\n", " prad -- radiation pressure [dyn/cm^2]\n", " gradr -- radiative temperature gradient\n", " gradT -- physical temperature gradient\n", " velocity -- velocity [km/s] \n", " conv_vel -- convective velocity [km/s]\n", " opacity -- opacity [cm^2/g]\n", " eps_nuc -- nuclear energy release rate, excluding neutrinos [erg/s/g]\n", " pp -- pp energy release rate [erg/s/g]\n", " cno -- CNO energy release rate [erg/s/g]\n", " tri_alfa -- triple-alpha energy release rate [erg/s/g]\n", " eps_nuc_neu_total -- energy loss rate as nuclear neutrinos [erg/s/g]\n", " non_nuc_neu -- energy loss rate as non-nuclear neutrinos [erg/s/g]\n", " eps_grav -- thermal energy release rate [erg/s/g]\n", " h1 -- 1H mass fraction\n", " he3 -- 3He mass fraction\n", " he4 -- 4He mass fraction\n", " c12 -- 12C mass fraction\n", " n14 -- 14N mass fraction\n", " o16 -- 16O mass fraction\n", " ne20 -- 20Ne mass fraction\n", " mg24 -- 24Mg mass fraction\n", " si28 -- 28Si mass fraction\n", " s32 -- 32S mass fraction\n", " ar36 -- 36Ar mass fraction\n", " ca40 -- 40Ca mass fraction\n", " ti44 -- 44Ti mass fraction\n", " cr48 -- 48Cr mass fraction\n", " fe52 -- 52Fe mass fraction\n", " fe54 -- 54Fe mass fraction\n", " fe56 -- 56Fe mass fraction\n", " ni56 -- 56Ni mass fraction\n", " eta -- electron chemical potential [kB*T]\n", " log_omega -- log10(rotation angular velocity [rad/s])\n", " v_rot -- rotation velocity [km/s]\n", " j_rot -- specific angular momentum [cm^2/s]\n", " dynamo_log_B_r -- log10(dynamo-generated radial field strength [gauss])\n", " dynamo_log_B_phi -- log10(dynamo-generated azimuthal field strength [gauss])\n", " log_D_conv -- log10(convective diffusivity [cm^2/s])\n", " log_D_semi -- log10(semiconvective diffusivity [cm^2/s])\n", " log_D_ovr -- log10(overshoot diffusivity [cm^2/s])\n", " log_D_thrm -- log10(thermohaline diffusivity [cm^2/s])\n", "\n" ] } ], "source": [ "# Print out documentation for the read_history function\n", "\n", "help(mw.read_profile)" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAEGCAYAAABmXi5tAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3deXhV1bn48e+beSSBJJCQgQAis4wqdUDUUsGZVutAi1oVa4vaem9H+/z0Vr31er0OrfWqFVEUR7SIFHACFLjIoAFlBhOGhDkESICQaf3+WDlTchIynCEn5/08z37Ons7eaxN9373XXmctMcaglFKq84sIdgGUUkoFhgZ8pZQKExrwlVIqTGjAV0qpMKEBXymlwkRUsAvQnPT0dJOfnx/sYiilVMj46quvDhljMrxt69ABPz8/nzVr1gS7GEopFTJEZGdT27RKRymlwoQGfKWUChMa8JVSKkxowFdKqTChAV8ppcKEBnyllAoTGvCVUipMdOh2+G32t79BaSlERdkpMtI1Hx8PCQmQmOj6dExdu0Jqqt1fKaU6mc4Z8J97DjZvbvv3U1KgWzebALp1g8xMO2VlNZ5PTQUR35VdKaX8pHMG/E2bwBiorYWaGtdnTQ2cPAknTsDx467pxAkoL4cjR+DwYSgrs5+HD9snhe3bYe9eOHWq8bliY23wz82FvDzXp/t8SoomBaVU0HXOgA82wDqqcXzBGDh6FPbts8F/3z7X/J49sHs3LF8OxcU2sbhLTnYF//x86NPHNfXubZ8SlFLKzzpvwPc1ERuYU1NhwICm96uthf37YdcumwR27XLN79wJq1fbpwZ3Xbs2TgKO+V69fJe0lFJhTSOJr0VGQs+edhozxvs+R49CUREUFrqmoiJYtw4++ACqqlz7RkdD376u6YwzoF8/6N/fPjHoC2alVAtpwA+GlBQYPtxODdXW2ioiRyLYuhW2bIHvvoMlS+w7B4fYWDjzTBv8G04pKQG7HKVUaNCA39FERtr6/txcuOgiz23G2OoiRxJwTOvWwT//aZOFQ48e3hNB795aRaRUmNL/80OJiKtZ6NixntuqquwTgXsi2LwZ3n/f851BdLQN/IMGeU79+kFMTGCvRykVUBrwO4uYGPsy2dsL5dJSVxLYtMlOa9bAu+/apwawTxb9+tngP3CgKxH0729/rKaUCnka8MNBWhqcd56d3J04YauHNm50TRs22BfHjuohEdtaqOETwYABkJQU+GtRSrWZBvxwlpDg/eXxqVOwbZtnIti4ERYuhOpq1369ejVOBAMH6gtjpTooDfiqsdhYGDLETu6qq+17goaJYPFiqKx07dezp6tayP0zI0N/caxUEIlx1OF2QKNHjzY6iHkIqK2FHTtsdZDjHcHGjfazosK1X1qa90SQk6OJQCkfEZGvjDGjvW3TO3zVfpGRrh+GXX21a70xtqsJ9wSwcSPMnm37KXJITrbvBBomgt699YdlSvmQ3uGrwDMGDh5snAg2bbI/OnOIjXU1IXVPBNqEVKkm6R2+6lhEoHt3OzX8cdmRI/b3A+6JYOVKePttzyakZ5zhai00cKD97N8funQJ/PUoFSL8FvBFpA/wAJBijLmuft044GFgA/CWMWaJv86vQlRqqu2DqGE/RCdO2N8RNHwi+PBDz95Js7M9k4BjPitL3xOosNeqgC8iLwNXAgeMMUPc1k8AngEigZeMMY8ZYwqB20VkttshDFABxAHF7S28CiMJCTBihJ3cVVfbfoY2b7bTpk3289VX7RgHDo73BA0TQd++9tfHSoWBVtXhi8hYbMCe6Qj4IhIJbAXGY4P4auAmY8zG+u2z3e7wI4wxdSLSA3jSGDO5ufNpHb5qM2PsWAXuScAxX1Li2i8qygZ9RyJwrx7S3xOoEOSzOnxjzBcikt9g9TnA9vo7ekTkLeAaYKOX79fVz5YBsa05t1KtIuLqpvqSSzy3lZe7uplwTwTz5nlWD2VmNu6a+swz7UtjTQYqBPmiDj8b2O22XAycKyJpwKPACBH5gzHmLyLyQ+AyIBV41tvBRGQqMBUgLy/PB8VTqoHkZBg92k7uqqvtuASORODolnrRIpg503PfjAxX8HdMZ55pk0JiYuCuRalWaHWzzPo7/HluVTrXA5cZY+6oX/4pcI4x5p72Fk6rdFSHcfKkDf7bttn+h7Ztc01793ru27OnKwG4J4S+fSEuLjjlV2HD380yi4Fct+UcYE8T+yoVmuLjvXc3AbaKaPt2zySwdasdo+DQIdd+InaUMm/JoHdvfXms/M4XAX810E9EegMlwI3AzT44rlKhITnZewsisL8r8PZUMGuWHerSITLSBv1+/VzDWDrm8/N10BrlE61tlvkmMA5IF5Fi4EFjzHQRmQZ8hG2W+bIxZoPPS6pUKEpNhbPPtpM7Y+zdf8Ongu3bYelSzz6IoqJsMnAkAveEoIPcq1bQrhWU6mgcQ1k6qoncq4u2b/eeDNwTgeNTk0FY0q4VlAol7kNZXnCB57aGycA9IXz+uecg9+7JoGFC0GQQlvQvrlQoaUky8PZU0DAZREd7VhO5J4S8PE0GnZT+VZXqLNyTwYUXem4zBvbt815N1FQy8FZNpMkgpOlfTqlwIGI7kMvKOn0ycE8IS5Y0nQy8VRPp+AUdmgZ8pcJdS5KBt2oib8mgT5+mq4k0GQSdBnylVNPck8HYsZ7b3JNBw4SweLHt0trBkQwaNisdNMj+Mlm7rg4IDfhKqbY5XTLYu9d7NdGiRZ7JoGtXGDrUcxoyRAez8QMN+Eop33PvrbSpZLB1qx34/ttv7TRzpucYBr16uRLAWWfByJH26SAiIrDX0olowFdKBZZ7Mhg3zrXeGNi1y5UAHNPCha5uq7t0sYF/1Cg7jR5tO6XTJNAiGvCVUh2DiL2r79ULrrzStb6qynZZ/dVXrunZZ+HUKbvdkQRGj7ZDY55/vm2aqhrRrhWUUqGnutqOa/zVV7Bmjf1ct86VBPr0sYH/ggvs58CBYfMU0FzXChrwlVKdQ1UVFBTA8uWwbJn9PHDAbuva1Qb/Sy+F73/ftg7qpC2DNOArpcKPMXbQGkcCWLLEthQCW+Vz6aV2Gj8ecnKCWlRf0oCvlFIAO3fCZ5/Z6dNPXU8Aw4fD1VfDVVfZ9wEhXP2jAV8ppRoyBtavt62APvzQPgnU1dnWQ1ddBTfcYJuUhtgvhDXgK6XU6Rw6BPPn2+C/YIHtNiIrC268EW6+2TYDDYF6/+YCfug+tyillC+lp8OUKfDuu7aq55134Nxz4e9/tyOWDR4MzzwDZWXBLmmbacBXSqmGEhLg+uvtQPT798NLL0FKCvzqV7bK57bbbIugEKMBXymlmpOaCrffDitW2CB/660we7Z9uTthAnzxhX0fEAI04CulVEsNHw7/+79QXAyPPWYTwEUX2W6lly0LdulOSwO+Ukq1VkoK/O53sGOH7eZhxw4b9K+/3rb976A04CulVFvFx8MvfwlbtsB//Idt5TNoEPzlL64O3zoQDfhKKdVeiYnw//6f7e//mmvgj3+0Hblt3BjsknnwW8AXkT4iMl1EZje3TimlOo2ePW1zznfftV09n302vP12sEvl1KqALyIvi8gBEVnfYP0EEdkiIttF5PcAxphCY8zt7vt5W6eUUp3OddfZ3jtHjLA/3Prd7zpES57W3uG/AkxwXyEikcDfgYnAIOAmERnkk9IppVSoysqywzn+/Ofw+ONw551QWxvUIrVqABRjzBcikt9g9TnAdmNMIYCIvAVcA7Sp8kpEpgJTAfLy8tpyCKWU6hhiYuC55yAjAx5+2Hbh/MorQeuczRdnzQZ2uy0XA9kikiYizwMjROQPAN7WNWSMedEYM9oYMzojI8MHxVNKqSASgT//2U6vvWZf6AaJL4Y49NabkDHGlAI/b7Cy0TqllAoLf/oTlJTAf/0XDBkCP/lJwIvgizv8YiDXbTkH2OOD4yqlVOchYn+kdcEFcPfdtglngPki4K8G+olIbxGJAW4E5vrguEop1blERcEbb0B0tO2ALcAtd1rbLPNNYAXQX0SKReR2Y0wNMA34CNgEvGOM2eD7oiqlVCeQmwv//d92wJWZMwN6ah0ARSmlAq2uzlbtFBVBYaHtosFHdAAUpZTqSCIibH87+/bZvvYDddqAnUkppZSLo1vlxx6D6uqAnFIDvlJKBctvfgN79sC//hWQ02nAV0qpYJk40XbBMH16QE6nAV8ppYIlKgomT4aFC6G83O+n04CvlFLBNGGCHSzl88/9fioN+EopFUznnw9xcfDpp34/lQZ8pZQKprg421rnk0/8fioN+EopFWzjx9vhEEtK/HoaDfhKKRVsY8faz5Ur/XoaDfhKKRVsZ50FkZFQUODX02jAV0qpYIuPh4EDNeArpVRYGDgQtm716yk04CulVEfQv7/tObOqym+n0ICvlFIdwZlnQm2t7TLZTzTgK6VUR9Cvn/3049CHGvCVUqojyMmxn35si68BXymlOoLMTDswSnGx306hAV8ppTqCqCgb9PUOXymlwkB2tgZ8pZQKCxrwlVIqTPTsCXv3+u3wGvCVUqqjSE+HsjLbHt8PovxyVC9E5EJgcv05BxljzvPXuT79FCoq7AvvyEj72dQUGQkxMaefogL2L6WUCltpaWAMHDli532sXWFMRF4GrgQOGGOGuK2fADwDRAIvGWMeM8YsBZaKyLXA6vac93SmTYMtW3x7zIgIG/jj4uwUH+/6dJ9v+JmQ4H1KTGx+XWSkb8uvlAoB3brZz9LSjhfwgVeAZ4GZjhUiEgn8HRgPFAOrRWSuMWZj/S43A3e087zNmjMHKivtU1FdXdNTba2dqqtt9xXu06lT3tedOgUnT9qpstI1f/IkHD7sua6yEk6csMdvrZiYlieHptY3TDgNl2NjQcT3//5KqTZyBPnDh/1y+HYFfGPMFyKS32D1OcB2Y0whgIi8BVwDbBSRPOCoMeZYU8cUkanAVIC8vLw2lWvAgDZ9zW9qamzgbzgdP97y9e7rSksbrz91qvXlEnElAfdk0HC+PZPjiSguziaYCH1rpFTTHAG/tNQvh/dHzXQ2sNttuRg4t37+dmBGc182xrwIvAgwevRo44fyBVxUFHTpYid/qa21TxXuScCx3HDe27L7upMn7TEOHnQtu0+mHX8V92qxhtVj/p70iUZ1eO5VOn7gj4Dv7X8pA2CMedAP51PYOv+kJDv5kzG2eqthEvA2VVa2fjpypOltbXmKaSg29vRJoSWJo6376LsZ1ayUFPtZXu6Xw/sj4BcDuW7LOcAeP5xHBYGIDWaxsZCaGthz19XZZNNUQmhJkjl50iYO9yTivv3oUdi/v3Giccy3V1RUyxOHP5JPTIw+5XRo8fH28+RJvxzeHwF/NdBPRHoDJcCN2Be1SrVLRIQrcAWD4+mmqaeP5hLN6ba779PcU44vmmf7IrEkJdmbUUdVpWPe8anVZ23kCPgnTvjl8O1tlvkmMA5IF5Fi4EFjzHQRmQZ8hG2W+bIxZkO7S6pUkLk/3TievAOtpsa3ycXbVFbW9PaWVqtFR7uCf2qqrZru2tU1NbfcpUsYJ4vISPsfWEcM+MaYm5pYPx+Y355jK6Uai4qyU2JicM5fV2eDfkUFHDtmp6NHXfMNl48etQmkrMz2+ltWZlscNtdUOSLCezLo1g0yMqB7d/vpPt+tWyd6PxIf3zEDvlIqvEREuJrcZmS07RjG2HjmCP6OhNBw2X2+sNA2XGmqeXpEhG3R6EgEjmTQvTtkZdkuarKz7ZSe3sGbByckhFQdvlJKNUnEPqEkJroGeWqp6mob+A8etNOBA6559+X16+28twQRHW0TgHsScEy9e0OfPjZRBK1aKSFB7/CVUio62o4RkpnZsv2rq22rq5IS79M338CCBfZ3J+4SE23gd5/69oXBgyE318/JQAO+Ukq1XnS0fYpo7knCGNvsffdu2LHDVh8VFsJ338H27fDxx541LMnJNvAPHgxDhsDQoTB6tA9f5GsdvlJK+YeIbRnkCOINGWOfErZtgw0b7LR+PXzwAUyf7jrGgAEwZgycey5ceCEMHNjGJwG9w1dKqeAQcVUjXXih57YDB2DdOli5Er78Ej78EGbUdx6TkwOXXQYTJsDEia1oWZWQAPv2+fQaHDTgK6VUG3XvDuPH2wns08B338HixfDRRzB7tn0KSEyEa6+Fn/7U7ttsKyE/3uF35MZJSikVUkTgjDPgzjttsD90yAb/yZPhX/+yd/uDBsE//tFMVx1+rMPXgK+UUn4SFQXjxsELL9hamlmz7N3+1Km2jv+997z0Pqt3+EopFdpiY+Hmm2HNGtvyJzkZrrvOrquocNvRjz+80oCvlFIBJGLr8b/+Gh59FN55x77cdf4WwFGl056BJ5qgAV8ppYIgKgr++Ed4+23bwuemm2xfRSQkuPoC9zEN+EopFUTXXQdPP22bdD72GPCLX9ifAUdH+/xcGvCVUirIpk2DG26ABx+EgqJU29GPH3p404CvlFJBJgLPPWfHDrjvPr9U3wMa8JVSqkPo1g0efhiWLrVt9v1BA75SSnUQd9xh++9/6SX/HF8DvlJKdRBRUbaJ5qpV/jm+BnyllOpAMjPtIC/+oAFfKaU6kKQk2wTfD83wtbdMpZTqSG64wQ6o4o9xdzXgK6VUB3LGGXbyh4BV6YjIQBF5XkRmi8jdgTqvUkopq10BX0ReFpEDIrK+wfoJIrJFRLaLyO8BjDGbjDE/B34MjG7PeZVSSrVee+/wXwEmuK8QkUjg78BEYBBwk4gMqt92NbAM+Kyd51VKKdVK7Qr4xpgvgMMNVp8DbDfGFBpjqoC3gGvq959rjDkPmNye8yqllGo9f7y0zQZ2uy0XA+eKyDjgh0AsML+pL4vIVGAqQF5enh+Kp5RS4ckfAV+8rDPGmCXAktN92RjzIvAiwOjRo/3UhZBSSoUff7TSKQZy3ZZzgD1+OI9SSqlW8EfAXw30E5HeIhID3AjM9cN5lFJKtUJ7m2W+CawA+otIsYjcboypAaYBHwGbgHeMMRvaX1SllFLt0a46fGPMTU2sn08zL2b9bc7mOVRUVRAVEUV0RDRREVHOKTEmkaSYJJJikkiOSSYpJomE6AREvL16UEqpzqNTdq3w+09/z5bSLS3eXxCSY5NJi08jLSGN9IR0O8Wnu+YT0slMyiQrOYvMpEwSohP8eAVKKeV7nTLgfzrlUyprKqmuraamroaauhqq66qprq3mRPUJKqoqnFN5VTkVVRUcO3WM0pOlHDpxiEMnDrH50GZKT5RSXlXu9RxdYruQlZTlTALZydn0TO5JdnI22V3sfM/knsRFxQX46pVSyrtOGfBzuuT47Finak5RerKUA8cPsK9iH/sq9rG3fK/9rNjL3oq9rC5ZzZzyOVTWVDb6frf4bq4kkNST7C7ZruRQnxi6J3YnQrSnaqWUf3XKgO9LsVGxzrv15hhjOFJ5hD3leygpL6HkWIlz3vG5bt869h/fT52p8/huVEQUWUlZziSQk5xDXkoeOV1yPBJEbFSsPy9VKdXJacD3ERGha3xXusZ3ZXD3wU3uV1NXw/6K/Z7J4FiJc37zoc188t0nXquSeiT2IC8lj9yUXPK65JGX4jl1T+yuL5+VUk3SgB9gURFR9q69SzZnc7bXfRxPC44nBcfn7mO72XV0F5sObuKj7R9xvPq4x/diI2NtMnAkgS55Hsu5XXJJjEkMxGUqpTogDfgdkPvTwpDuQ7zuY4yhrLKM3UdtEnBOx+znp4Wfsqd8T6Pqo7T4NPJS8shPzad3am96d+1N79Te5Kfmk5+arwlBqU5MA36IEhG6xXejW3w3hmUO87pPdW01e8r3eCaEo7vYeXQnmw9tZuH2hZysOenxne6J3V3JoD4hOJbzUvL0PYJSIUwDficWHRlNr9Re9Ert5XW7MYb9x/ez48gOisqKKDpSRFFZETuO7mDNnjW8t+k9aupqnPsLQs/kns6nAseTgWM5u0s2URH6n5RSHZUY03E7pBw9erRZs2ZNsIsRtmrratlTvseZCIqOFNnkUL9cfKwYg+u/n6iIKHK75HokBOd81970SOyhL5WV8jMR+coY43VUQb0dU02KjIgkNyWX3JRcxvYa22h7VW0Vu4/u9poQ5m2dx/7j+z32j4+K93giaJgQUuNSA3VpSoUlDfiqzWIiY+jbrS99u/X1uv1E9YlG1UVFR4ooLCtk2a5lHDt1zGP/1LhUzyTQ4KVyfHR8IC5LqU5LA77ym4ToBAZlDGJQxqBG2xytjBomg6IjRWw4sIF/bf0Xp2pPeXwnMynTIwn06drHuZzTJUffHyh1GlqHrzqkOlPHvop9XhNCUVkRu4/t9mhyGimR5KXk6fsDFfa0Dl+FnAiJcHZpcX7e+Y22V9dWs/vYbq8JoSXvD/JT8+mVYlsw9Urppb9SVmFBA74KSdGR0fTp2oc+Xft43d7U+4OiI0Us37Wco6eOeuwfFxVHXkqeTQJuicDxqU1OVWeg/wWrTqm59wcARyqPsPPITnYe3en6rJ9ft38dB44f8Ng/UiLJ7pLtmQzc5vNS8vSlsurwNOCrsJQal0pqZmqTv1I+WX3S+avkhgnhi51fUHKshFpT6/GdjIQMr8kgNyWX3C65pCeka7WRCioN+Ep5ER8dT//0/vRP7+91e01dDSXHSjwTQv3ntwe+5V/b/tVofIS4qDhyuuSQ2yXXmQQ85lNySYlN0aSg/EYDvlJtEBUR5eq2wkvPFcYYDhw/wK6ju9h9bDe7j+62n/Xzi4oWee3cLikmqdmEoD2eqvbQgK+UH4gIPZJ60COpB2dne+8Gu6auhr3lez0TglticAyY01DXuK5NJgTH6Go65rLyRgO+UkESFRHl7LqCXO/7nKo5RUl5ideEsPvoblYUr+DwycONvucYWtMxzrL7kJrZydlkJWfRPbG7tjwKM/rXVqoDi42Kbbb5KcDxquMUHyumpLyE4mPFzskxqto3+7/xOrRmhETQI7EHWcl2eM2eST1d88k9ncNudk/sTmREpL8vVQVAwAK+iPQBHgBSjDHXBeq8SnV2iTGJzb5gBtfQmiXlJewt38ue8j3srbCfe8r3UHysmFUlqxo1RwWbGDKTMp0JwD0Z9Ey2SSIzKVOfGEJAi/46IvIycCVwwBgzxG39BOAZIBJ4yRjzWFPHMMYUAreLyOz2FVkp1VruQ2s2p7q2mv3H9zsTgSM5OBLErqO7+LL4Sw6eONjou4KQkZjhTA6ZSZmN5+uTQ3JMsrZGCoKWpuNXgGeBmY4VIhIJ/B0YDxQDq0VkLjb4/6XB939mjGl866CU6lCiI6PJ6ZJDTpecZverqq1if8V+ZyLYV7GPveX2c99xO7/x4Eb2Veyjuq660fcTohOcScA9KWQlZdmmq/UvopNjk/11qWGpRQHfGPOFiOQ3WH0OsL3+zh0ReQu4xhjzF+zTQJuIyFRgKkBeXl5bD6OU8qOYyBjXC+dmGGM4fPKwTQQV+5zJwX1+08FNLC5aTFllWaPvp8SmNNtENadLjv7CuRXaU+GWDex2Wy4Gzm1qZxFJAx4FRojIH+oTQyPGmBeBF8H2ltmO8imlgkxESEtIIy0hjcHdBze776maU+yt2Ntki6TVe1Zz6MShRt/LTMp0vtju27Wvx2dmUqZWHblpT8D39q/YZIA2xpQCP2/H+ZRSnVhsVCz5qfnkp+Y3uc/J6pMUHyt2JoFdR3ex48gOviv7js93fM6sb2Z5DLsZHxVvE0C3vgxIG8DAjIEMTB/IgPQBpMSlBOCqOpb2BPxiPFsP5wB72lccpZRqWnx0PP3S+tEvrZ/X7adqTrHjyA4KywopLCvku7LvKCwrZPvh7SzYtsDjfULP5J4MTLcJYGDGQIb1GMawzGEkxSQF6nICrj0BfzXQT0R6AyXAjcDNPimVUkq1QWxUbJNNVGvqaigsK2TTwU1sOlQ/HdzEK+teoaKqArAtjc5MO5ORWSMZkTmCkVkjGZk1kq7xXQN9KX7RohGvRORNYByQDuwHHjTGTBeRy4GnsS1zXjbGPOrLwnkb8aq6upri4mIqKyub+JYKlri4OHJycoiOjg52UZRqMWMMu4/tZu2+tRTsLeDrfV9TsLeA3cdcrygHZwzmgrwLnFOvlF4d9t1AcyNehdwQh0VFRSQnJ5OWltZh/8HDkTGG0tJSysvL6d27d7CLo1S7HTpxiIK9BawqWcXy3ctZvns5x04dAyA7OZvxfcdz+RmXM77veFLjUoNcWpdONcRhZWUl+fn5Guw7GBEhLS2Ngwcb/yBHqVCUnpDO+L7jGd93PAC1dbVsOLiBZbuWsWTHEuZsnsMra18hUiL5Xu73uLb/tdw45MbT/rgtmELuDn/Tpk0MHDgwSCVSp6N/HxUuaupqWFm8kgXbFzB/23wK9hUgCBflX8TkoZO5achNQenKurk7/IhAF0YppTqDqIgozs87n0cueYSv7/qaLdO28OBFD1JyrIQ7P7yTnKdy+O0nv2XX0V3BLqqTBvw2iIyMZPjw4QwZMoTrr7+eEydONNrHGMMll1zCsWPHGm176KGHeOKJJ9p8/urqakaNGtXm7zd01113kZiYyKJFizzWP/nkkwwaNIizzjqLSy+9lJ07dwJw8OBBJkyY4LPzK9UZnJl2Jg+Oe5At07aw9LalfL/P9/mfFf/DGX89g3sX3Ou1Y7pA04DfBvHx8axdu5b169cTExPD888/32if+fPnM2zYMLp06eLz8y9btozzzjvPJ8d65JFHKCsrY+XKlfzyl7/km2++cW4bMWIEa9as4ZtvvuG6667jt7/9LQAZGRlkZWWxfPlyn5RBqc5ERLgg7wLevf5diu4r4rbht/Hc6ufo+9e+PLf6uUbdVAdSyL209fCrX8Hatb495vDh8PTTLd79wgsv9AiSDrNmzWLq1KnO5UcffZSZM2eSm5tLRkYGo0aN4rvvvuP666/n66+/BmDbtm3ceOONfPXVV+Tn53PLLbfw4YcfUl1dzbvvvsuAAQMAWLhwIRMnTgRg5syZPPHEE4gIZ511Fq+99hq33nor8fHxbN68mZ07dzJjxgxeffVVVqxYwbnnnssrr7wCwKuvvgbI8EAAABJ7SURBVMr69et54403iIqKYu7cufzkJz9h9uzZ5ObmcvHFFzvLP2bMGF5//XXn8rXXXsusWbM4//zzW/5vq1SYyUvJ44WrXuD+793PvQvv5Zfzf8n7m97nzR+9SUZiRsDLo3f47VBTU8OCBQsYOnRoo23Lly93Vrt89dVXvPXWWxQUFPD++++zevVqAPr27UtKSgpr65PWjBkzuPXWW53HSE9P5+uvv+buu+/2qAJavHgx48aNY8OGDTz66KMsWrSIdevW8cwzzzj3KSsrY9GiRTz11FNcddVV/PrXv2bDhg18++23zvPdcsstvPXWW0RF2bzfr18/Vq5cSW5u4w6xpk+f7kwyAKNHj2bp0qVt/adTKqz0T+/PwskLeeHKF1i+eznnvHQOGw9uDHg5QvsOvxV34r508uRJhg8fDtg7/Ntvv73RPocPHyY52XbtunTpUiZNmkRCgh1n9Oqrr3bud8cddzBjxgyefPJJ3n77bVatWuXc9sMf/hCAUaNG8f777wOwZ88eunXrRkJCAosWLeK6664jPT0dgG7dujm/e9VVVyEiDB06lB49ejiT0uDBg9mxY4ez/C3x+uuvs2bNGj7//HPnuu7du7Nnj/akoVRLiQhTR01lROYIrn7rai559RKW3ra0yW4i/EHv8NvAUYe/du1a/va3vxETE9Non6ioKOrqXHV1Tf1u4Ec/+hELFixg3rx5jBo1irS0NOe22NhYwL4krqmpAWDBggVcdtllgH0x3NRxHd+NiIhwzjuWHcdqiU8//ZRHH32UuXPnehynsrKS+Hjtllap1jo7+2wW37KYWlPLD17/gdcxif1FA76f9O/fn8LCQgDGjh3LP//5T06ePEl5eTkffvihc7+4uDguu+wy7r77bm677bbTHte9/v7SSy/lnXfeobS0FLBPFb5UUFDAXXfdxdy5c+nevbvHtq1btzJkyJAmvqmUas6A9AHMu2keJcdKuHXOrQTq91Aa8P3kiiuuYMmSJQCMHDmSG264geHDh/OjH/2ICy+80GPfyZMnIyL84Ac/aPaYtbW1bNu2zfnydvDgwTzwwANcdNFFDBs2jPvvv9+n1/Cb3/yGiooKrr/+eoYPH+5RFbV48WKuuOIKn55PqXBybs65PD7+cT7c+iGzvp0VkHPqL239ZO/evUyZMoVPPvnktPs+8cQTHD16lIcffrjZ/ZYtW8brr7/utRlooI0dO5YPPviArl09exEMlb+PUh1BnaljzEtj2H1sN9vv2e6TX+bqL22DICsrizvvvNPrD6/cTZo0iZkzZ3Lfffed9pgXXHBBhwj2Bw8e5P77728U7JVSrRMhETx12VPsq9jHywUv+/18eoevfEr/Pkq13sgXRhITGcOXd3zZ7mPpHb5SSnVgNwy+gZUlKykqK/LreTTgK6VUkP148I8BmL1xtl/PowFfKaWCrHfX3pzR7QxWFK/w63k04CulVAcwMmskBfsK/HoODfg+kJTUvlHu58yZw5///OcW7TthwgRKSkradT6HN954g5iYGB555BGP9atWrWL48OEMHz6cYcOG8c9//hOAqqoqxo4d26pf6iqlWmZk5kh2HNlB2ckyv51DA347GGM8uk9oq8cff5xf/OIXp93v5MmTHD58mOzs9g+htmjRIh5//HE2btzIJ5984uxBE2DIkCGsWbOGtWvXsnDhQu666y5qamqIiYnh0ksv5e233273+ZVSnkZkjQDw611+SHee9quFv2LtPt92jzw8czhPT2i6U7YdO3YwceJELr74YlasWMGcOXMAeOCBB5g3bx7x8fF88MEH9OjRg507d/Kzn/2MgwcPkpGRwYwZM8jLy/M43tatW4mNjXV2gPbdd98xefJkamtrmThxIk8++SQVFRUALFmyhHHjxgGwevVq7rvvPo4fP05sbCyfffYZ7733HnPmzKG2tpb169fzb//2b1RVVfHaa68RGxvL/Pnz6datG99++y1/+tOf+Oijj+jRowfz589n0qRJZGVlcdlllzk7eQPbZ457fz3XXnstf/jDH5g8ebJP/r2VUtaIzPqAv7eAS3pf4pdz6B1+G2zZsoUpU6ZQUFBAr169OH78OGPGjGHdunWMHTuWf/zjHwBMmzaNKVOm8M033zB58mTuvffeRsdavnw5I0eOdC7fd9993HfffaxevZqePXt67LtgwQImTJhAVVUVN9xwA8888wzr1q3j008/dXZk5ujfftWqVTzwwAMkJCRQUFDA9773PWbOnAnA0KFD+b//+z969OgBQGJiIh9//LGzUzaAlStXMnjwYIYOHcrzzz/v7EJ5yJAhzu6dlVK+k5GYQWZSJhsObvDbOUL6Dr+5O3F/6tWrF2PGjHEux8TEcOWVVwK2K2NHdworVqxwdmv805/+1DlilLu9e/eSkeEaCMH9qeHmm2/m3//9353bli9fzhNPPMGWLVvIysri7LPPBvAYVeviiy8mOTmZ5ORkUlJSuOqqqwAb5L0N1NKUc889lw0bNrBp0yZuueUWJk6cSFxcHJGRkcTExFBeXu7s/lkp5Rt9uvah6Ij/2uIH7A5fRMaJyFIReV5ExgXqvP6QmOjZ30V0dLSz2sO9K+OGvHVlHB8fT2Vl5WnPWVhYSG5uLjExMS3qFhk8u0ZubbfIDgMHDiQxMZH169c71506dYq4uLhWH0sp1bzeqb3ZcWSH347fooAvIi+LyAERWd9g/QQR2SIi20Xk96c5jAEqgDiguG3FDS3nnXceb731FmCHPLzgggsa7TNw4EC2b9/uXB4zZgzvvfcegPO74KrOARgwYAB79uxxVq2Ul5f7tOVMUVGR83g7d+5ky5Yt5OfnA1BaWkpGRgbR0dE+O59SyuqR2IP9Ffv9dvyW3uG/AkxwXyEikcDfgYnAIOAmERkkIkNFZF6DqTuw1BgzEfgd8B++u4SO669//SszZsxwjjXrPgShw9ixYykoKHD2h/3000/z5JNPcs4557B3715SUlIA2w++I+DHxMTw9ttvc8899zBs2DDGjx/foqeEllq2bBnDhg1j+PDhTJo0ieeee875Unnx4sVcfvnlPjuXUsolIzGDkzUnOV513D8nMMa0aALygfVuy98DPnJb/gPwhxYcJwaY3cz2qcAaYE1eXp5paOPGjY3Whbp7773XfPLJJ8YYY44fP27q6uqMMca8+eab5uqrrzaVlZVm1KhRwSyi06RJk8zmzZub3N4Z/z5KBcpLX71keAhTVFbU5mMAa0wT8bU9L22zgd1uy8XAuU3tLCI/BC4DUoFnm9rPGPMi8CLY3jLbUb6Q8cc//pGVK1cCdsDzadOmYYwhNTWVl19+mdjYWBr2GhoMVVVVXHvttfTv3z/YRVGqU8pItA04Dp04RH5qvs+P356A7+2tYZMB2hjzPvB+O87XafXo0cM5mtSFF17IunXrglwi72JiYpgyZUqwi6FUp5WRYAP+weMH/XL89rTSKQZy3ZZzgD3tK45SSoUvxx3+wRMdL+CvBvqJSG8RiQFuBOb6plhKKRV+OsQdvoi8CawA+otIsYjcboypAaYBHwGbgHeMMf77iZhSSnVyXWK7EB0R7bc7/BbV4Rtjbmpi/Xxgvk9LpJRSYUpEyEjM6JB1+GGptLTU2XVwZmYm2dnZzuWqqirnfsYYLrnkEucg5u5dKM+fP59+/fqxa9cunn32WWbMmBHw61BKdUzd4rtxuPKwX44d0n3pBENaWhpr19oeOh966CGSkpI8+rtxmD9/PsOGDfPo5wbgs88+45577uHjjz8mLy+Pn/3sZ5x//vncdtttASm/UqpjS4pJ8tsPr0I64P/qV7DWt70jM3w4PO2DPtlmzZrF1KlTPdYtXbqUO++8k/nz59O3b18AEhISyM/PZ9WqVZxzzjntP7FSKqQlxSRRUVXhl2NrlY6fLF++nFGjRjmXT506xTXXXMOcOXMYMGCAx76jR49m6dKlgS6iUqoD8mfAD+k7fF/cifvL4cOHPboPjo6O5rzzzmP69OmN+tTp3r07mzdvDnQRlVIdUFJMEser/VOlo3f4fhIVFeUx/GFERATvvPMOq1ev5j//8z899q2srHQOYKKUCm9J0VqlE3L69+9PYWGhx7qEhATmzZvHrFmzmD59unP91q1bGTJkSKCLqJTqgBJjEjXgh5orrriCJUuWNFrfrVs3Fi5cyCOPPMIHH3wA2Pr+73//+wEuoVKqI0qKSeJE9Qlq62p9fuyQrsMPtoceeqjJbXfccQdTpkzhjjvuAHAORA6Qm5tLUZEdxqygoIDBgwc7+5tXSoW3YT2GceOQG6mpqyEyItKnx9aA7ydZWVnceeedHDt2rFFbfHeHDh3i4YcfDmDJlFId2aSBk5g0cJJfjq0B349+/OMfn3af8ePHB6AkSikVonX4xoTFuCghR/8uSnVsIRfw4+LiKC0t1eDSwRhjKC0tJS4uLthFUUo1IeSqdHJyciguLubgQf/0JqfaLi4ujpycnGAXQynVhJAL+NHR0fTu3TvYxVBKqZATclU6Siml2kYDvlJKhQkN+EopFSakI7d2EZGDwM42fj0dOOTD4oSCcLvmcLte0GsOF+255l7GmAxvGzp0wG8PEVljjBkd7HIEUrhdc7hdL+g1hwt/XbNW6SilVJjQgK+UUmGiMwf8F4NdgCAIt2sOt+sFveZw4Zdr7rR1+EoppTx15jt8pZRSbjTgK6VUmAjpgC8iE0Rki4hsF5Hfe9keKyJv129fKSL5gS+lb7Xgmu8XkY0i8o2IfCYivYJRTl863TW77XediBgRCfkmfC25ZhH5cf3feoOIvBHoMvpaC/7bzhORxSJSUP/f9+XBKKeviMjLInJARNY3sV1E5K/1/x7fiMjIdp/UGBOSExAJfAf0AWKAdcCgBvv8Ani+fv5G4O1glzsA13wxkFA/f3c4XHP9fsnAF8CXwOhglzsAf+d+QAHQtX65e7DLHYBrfhG4u35+ELAj2OVu5zWPBUYC65vYfjmwABBgDLCyvecM5Tv8c4DtxphCY0wV8BZwTYN9rgFerZ+fDVwqIhLAMvraaa/ZGLPYGHOifvFLINT7K27J3xngYeBxoDKQhfOTllzzncDfjTFlAMaYAwEuo6+15JoN4BgvNAXYE8Dy+Zwx5gvgcDO7XAPMNNaXQKqIZLXnnKEc8LOB3W7LxfXrvO5jjKkBjgJpASmdf7Tkmt3djr1DCGWnvWYRGQHkGmPmBbJgftSSv/OZwJkislxEvhSRCQErnX+05JofAn4iIsXAfOCewBQtaFr7//tphVx/+G683ak3bGPakn1CSYuvR0R+AowGLvJrifyv2WsWkQjgKeDWQBUoAFryd47CVuuMwz7FLRWRIcaYI34um7+05JpvAl4xxvyPiHwPeK3+muv8X7yg8Hn8CuU7/GIg1205h8aPeM59RCQK+xjY3CNUR9eSa0ZEvg88AFxtjDkVoLL5y+muORkYAiwRkR3Yus65If7itqX/bX9gjKk2xhQBW7AJIFS15JpvB94BMMasAOKwnYx1Vi36/701Qjngrwb6iUhvEYnBvpSd22CfucAt9fPXAYtM/duQEHXaa66v3ngBG+xDvV4XTnPNxpijxph0Y0y+MSYf+97iamPMmuAU1yda8t/2HOwLekQkHVvFUxjQUvpWS655F3ApgIgMxAb8zjzW6VxgSn1rnTHAUWPM3vYcMGSrdIwxNSIyDfgI+4b/ZWPMBhH5M7DGGDMXmI597NuOvbO/MXglbr8WXvN/A0nAu/Xvp3cZY64OWqHbqYXX3Km08Jo/An4gIhuBWuA3xpjS4JW6fVp4zf8G/ENEfo2t2rg1lG/gRORNbJVcev17iQeBaABjzPPY9xSXA9uBE8Bt7T5nCP97KaWUaoVQrtJRSinVChrwlVIqTGjAV0qpMKEBXymlwoQGfKWUChMh2yxTqWARkT9g24T3w/ZpcxDb4dfDxpg3g1k2pZqjd/hKtVD9D2AigB8AH9evfsoYMxzb0dULIhIdtAIqdRoa8JVqhojki8gmEXkO+Br7U/cYY4zHLzyNMduwP47pGoRiKtUiGvCVOr3+2G5qRwCjgM8a7lA/OMW2TtKdheqkNOArdXo76/sjB5iAZ5fTvxaRLcBKbPe9SnVYGvCVOr3jbvPnAKvclp8yxvQHbgBmikhcQEumVCtowFeqhURkMLDZGFPbcJsx5n1gDa7eWZXqcDTgK9VyE4GFzWz/M3B/fUsepToc7S1TqRYSkU+AKe3tk1ypYNGAr5RSYUIfPZVSKkxowFdKqTChAV8ppcKEBnyllAoTGvCVUipMaMBXSqkw8f8B5xNnVQTChpwAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# Define constants\n", "\n", "RSUN = 6.957E10 # Solar radius, in cm\n", "MSUN = 1.989E33 # Solar mass, in g\n", "\n", "# Plot pressure, density and temperature versus fractional radius\n", "\n", "P = prof_data['pressure']\n", "rho = 10**prof_data['logRho']\n", "T = 10**prof_data['logT']\n", "\n", "r = prof_data['radius']*RSUN # Convert radii to cm\n", "R = prof_data['photosphere_r']*RSUN # Same\n", "\n", "plt.figure()\n", "\n", "plt.plot(r/R, P, color='r', label='P (dyn/cm^2)')\n", "plt.plot(r/R, rho, color='g', label='rho (g/cm^3)')\n", "plt.plot(r/R, T, color='b', label='T (K)')\n", "\n", "plt.xlabel('r/R')\n", "\n", "plt.yscale('log')\n", "\n", "plt.legend()" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAEJCAYAAACXCJy4AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3deXhU933v8fd3tAICsSM2IQwYMMTERMZO4iR2bGycGJOndn1NkydxTM1N0+TpTbNcN21u2iZP6b1Jb9w89o1LEpuQtGCXUgevmMRxvBQ7krxgViNjAwIMYhNikdBI3/vHkdCClpFmNHNG83k9z3nOmTNnznyPBJ/56XfO/I65OyIiMvBFUl2AiIgkhwJfRCRDKPBFRDKEAl9EJEMo8EVEMoQCX0QkQyjwRUQyhAJfRCRDZCfrjcxsCPD/gPPA8+7+rz29ZvTo0V5SUtLfpYmIDCgVFRVH3X1Mx/VxBb6ZPQTcAhxx97lt1i8C/hnIAn7m7v8I/BGwzt0fN7NHgB4Dv6SkhPLy8nhKFBHJOGa2t7P18XbprAIWdXijLOAB4GbgMmCpmV0GTAL2N2/WGOf7iohIL8UV+O7+AnC8w+oFQKW773H388BaYAlQRRD6cb+viIj0Xn8E70RaW/IQBP1EYD1wm5n9BHi8qxeb2XIzKzez8urq6n4oT0QkM/XHSVvrZJ27+xngiz292N1XAisBSktLNZSniEiC9EcLvwqY3ObxJOBgb3ZgZovNbGVNTU1CCxMRyWT9EfhlwAwzm2pmucCdwIbe7MDdH3f35YWFhf1QnohIZoor8M1sDbAZmGlmVWa2zN2jwFeAjcAO4FF339bL/aqFLyKSYBbmO16Vlpa6rsPPQO7Q1BRMjY3B1Nflrta5t75Py3Jfp3j30dnr2/4sOv5serPcl9cMhPccCO6+G/r4xVMzq3D30o7rk/ZN294ws8XA4unTp/fp9ecbz7PnxB5mjZ6V2MLSmTucOxdMdXXtp/r6i9d1NTU09G6KRnvepiWQW8K4qSnVPy1JV9bZNSNp6oYb+hz4XRmQLfw//vc/5sW9L/LYnY9x9aSr+6GyJKmvh5Mng6mmpnW55XFtLZw5A6dPB/O2yx3XnT0bfwsoLw9ycyEnJ7YpOzv27SIRyMoKpr4ux7ptJBIEQ8s8ninefXT3+hYdQ6yr5xL5mnR5T+lUWrXw4/Wdj3+HW/7tFj788w+z8JKF3DP/HpbMWkJuVm5qCnIPQre6OpiOHGldbrvuxIn2gV5X1/1+zWDIkNapoKB1edy4ztcPGgT5+bFPeXnBPDdX/9FE0lwoW/htunTu2b17d5/2UVtfy49e+RE/e+1n7D+1nzGDx7B07lIWTV/EJ0o+weCcwYkp1j0I6n372k8HDsDBg63zM2c6f/2gQTBmTDCNGgXDh7dOhYVdPy4sDAJcISwiHXTVwg9l4LdIxEnbxqZGnn3nWX762k95uvJp6qJ15GXl8bEpH+OmaTdx07SbmDt2LtZTcNbXw/bt8OabwbRzZ2u4nz7dftu8PJg4ESZMaJ2KimDs2GBqCfixY4PQFhFJoIwN/LbONZzjxX0vsrFyIxvf2ci26uBq0fEF47lx2o3cNO0mrr/kesbmj4KtW+Hll2HzZnjjjSDgo9FgR4MGwaxZwQmVKVOguDiYWpbHjFHLW0RSJq0CPxFdOrGoOlXFs+88y8adT7JpzyZORGsBKD5llFY5Vx6E0nMj+NDEKxkxtxTmzQum6dODE4AiIiGUVoHfol+vw9+3D1avhvXr4c03afQmyifCSwvGUz5jCOXDaqlsOHxh82kjpnHjtBu5csKVzB4zm9mjZ1OYr28Ci0j4KPAhuDRx/XpYtQqeey444XrNNXDddfDRj8LVVwcnQ5udOHeCikMVlB8s55WqV9i0ZxNnG85eeH58wfgL4T979Gxmj5nNrNGzGF8wvudzAiIi/SSzA7+2Fr73PXjwwWB56lS46y74/Od79cWGaFOUd0+8y46jO9hRvYPtR7ez8+hOdlTvoPZ87YXtCvMKmTV6FrPHzGbmqJlMKZzClOFTKC4sZnzBeLIi6g4Skf6TmYHvDmvXwje/GVwe+dnPwvLlQas+krhx49ydg7UH2XF0x4UPgJblQ6cPtds2O5LNpGGTKC4sDqZhxRc+DFqmgtyChNUmIpknrQI/ISdtjxyBO+6A3/8e5s+H+++HD384oXXGora+lv2n9rOvZh97T+5lX80+9p3ad+Fx1akqGr39HR9HDhrZ7gOhuDD4UJgwdALjhoyjqKCIgtwCdRuJSKfSKvBb9LmFX1cXtOK3b4f77oNly0J7VU1jUyOHTh9q/TBonvbWtD6uqb941NBB2YMoKiiiqKCIcQXjKBrSZrmgiHFDxjF68GhGDBrB8PzhZEcG5JeqRaQTGTW0AvfdBxUV8Otfw623prqabmVFspg0bBKThk3io3y0021q6mrYV7OPQ6cPcfj0Yd4//T7vn36fw2eC5crjlby07yWOnj3a5fsU5hUyYtAIRg4ayYj8EYwaPIpRg5qnwaMYnj+cYXnDGJo7lGF5w9pNQ3KHEDHdhlgk3Q3MwG9ogMWLQx/2sSrML+QD+R/gA+M+0O12DY0NVJ+tvvCBcPzccY6fO86JcyeCed2JC+v2v7+fY2ePcaLuBE3e/eiUhjE0r/WDYGjuUIbmDaUgt4CC3AKG5AxhcM5gBucMbr+c2/P6vKw8dU2JJMnADPzvfCfVFaRETlYOE4ZOYMLQCTG/psmbOFl3klP1py6aautrL15/vvW5w6cPU3u+lrMNZznbcJYz58/g9K6LMGIR8rPzycvKIy8778K8s3Xt5lnN23TyfE4kh+xINjlZzfNITq+XW17fcTliEX1ASdoKZeDHOx6+xC5iEUYOGsnIQSPj3pe7U99Y3+4D4MJyw5ku19dF66iP1lPfGEztHkfrORc9x8m6kxce10XrLiy3zHv7QROPnEgOWZEsIhbp9WRY719jsb3GyJwPokz40P3+dd9nztg5Cd1nKAPf3R8HHi8tLb0n1bVI7MyM/Ox88rPzE/IBEit3J9oUvRD+0aYoDU0NNDQ2dLocbYrS0NgQ03JX+2r0RtydJm/qeaL945hf18kUbYp2ur6xqbHnH9QAkcwP91Sqi/YwPHofhDLwRXrDzIIul6wcfYdBpBu69EJEJEMo8EVEMoQCX0QkQyjwRUQyRCgD38wWm9nKmpqLhxQQEZG+CWXgu/vj7r68sFA3GBERSZRQBr6IiCSeAl9EJEMo8EVEMoQCX0QkQyjwRUQyhAJfRCRDKPBFRDJE0gLfzC4xs5+b2bpkvaeIiLSKKfDN7CEzO2JmWzusX2Rmu8ys0szu7W4f7r7H3ZfFU6yIiPRdrOPhrwLuB1a3rDCzLOABYCFQBZSZ2QYgC1jR4fV3u/uRuKsVEZE+iynw3f0FMyvpsHoBUOnuewDMbC2wxN1XALckskgREYlfPH34E4H9bR5XNa/rlJmNMrMHgSvM7K+62W65mZWbWXl1dXUc5YmISFvx3OKws7sId3mzSXc/Bnypp526+0pgJUBpaWlm3LxSRCQJ4mnhVwGT2zyeBByMr5yAhkcWEUm8eAK/DJhhZlPNLBe4E9iQiKI0PLKISOLFelnmGmAzMNPMqsxsmbtHga8AG4EdwKPuvi0RRamFLyKSeOYe3m7y0tJSLy8vT3UZIiJpxcwq3L2043oNrSAikiFCGfjq0hERSbxQBr5O2oqIJF4oA18tfBGRxAtl4KuFLyKSeKEMfBERSTwFvohIhghl4KsPX0Qk8UIZ+OrDFxFJvFAGvoiIJJ4CX0QkQ4Qy8NWHLyKSeKEMfPXhi4gkXigDX0REEk+BLyKSIRT4IiIZIpSBr5O2IiKJF8rA10lbEZHEC2Xgi4hI4inwRUQyhAJfRCRDKPBFRDKEAl9EJEMo8EVEMkQoA1/X4YuIJF4oA1/X4YuIJF4oA19ERBJPgS8ikiEU+CIiGUKBLyKSIRT4IiIZQoEvIpIhFPgiIhkiqYFvZp8xs5+a2a/N7MZkvreISKaLOfDN7CEzO2JmWzusX2Rmu8ys0szu7W4f7v6Yu98D3AX8tz5VLCIifZLdi21XAfcDq1tWmFkW8ACwEKgCysxsA5AFrOjw+rvd/Ujz8t80v05ERJIk5sB39xfMrKTD6gVApbvvATCztcASd18B3NJxH2ZmwD8CT7v7a529j5ktB5YDFBcXx1qeiIj0IN4+/InA/jaPq5rXdeWrwA3A7Wb2pc42cPeV7l7q7qVjxoyJszwREWnRmy6dzlgn67yrjd39x8CPe9yp2WJg8fTp0+MoTURE2oq3hV8FTG7zeBJwMM59arRMEZF+EG/glwEzzGyqmeUCdwIb4i1K4+GLiCReby7LXANsBmaaWZWZLXP3KPAVYCOwA3jU3bfFW5Ra+CIiidebq3SWdrH+KeCphFUkIiL9It6Ttv0i3pO254+eJ5IfIWtIFsGVoPFx7/I8dIcNY91hjNtFSEj9IiIQ0sB398eBx0tLS+/py+vf/fa7HPrpoeBB27yMNWhDJJIfwfKMSH6ESF6ESH6EnFE55Bblkj81n8EzBzP82uEMvnRwqksVkZALZeDHa+zSseRfkk/T2aaLn+yswRxDIzrmlnasDfKetnPwqNNU10RTfVPr/FwTDccaOPv2WY4/c5ymc8Ex5k3KY+iCoYxaPIoxt48hu2BA/mpFJA4Wc3dFErXp0rln9+7dqS4ntLzJqXu3juPPHKfm5RpqXq6hfl89kcERir9VTPG3i4nkaEBUkUxjZhXuXnrR+jAGfovS0lIvLy9PdRlpw905tfkUVf9cRfWj1eSOz2X8n45n0tcmkTMiJ9XliUiSdBX4av4NIGZG4UcKmfPIHC5/5nIKrihg7/f38vo1r3Nuz7lUlyciKRbKwNcXr+I38qaRXP7k5cx7bh71++v5w6w/sPsvdtNwvCHVpYlIioQy8PXFq8QZce0IFuxYQNFdRRx44AAVH6qg9o3aVJclIimgSzkyQN7EPGaunMn4ZePZettWXv/I68xZP4dRi0alurR2olE4dgxqauDUKaitvXh+5gw0NLRO0Wj7x91N0Si4dz01NXX/fHdTf0rGabaBcAwDzX/8B1x9dWL3qcDPIMOuGkZpRSlbFm1hx9IdzH9lPoNnJuf6/YYG2LMnmN59N5gfOABHjsDhw8H86NHYgiEnB7Kzg3msU34+ZGVBJAJmF09dre/N1J+S8f27gXAMA8nIkYnfZyiv0tFlmf3r3LvnqLiyAq93pt83naK7ixL6jd76eti2DSoq4LXXgvmWLcH6Fvn5MHEijBsXTGPHts6HD4dhw2Do0IvngwYpOER6ossypZ26/XXsvGsnJ587yZTvTmHq307t877cYfNm+MUvoKwMtm4NWvQAhYUwf34wXX45TJsGU6dCUVHQqhaRxOsq8NWlk6HyJ+czb9M8dn5hJ3u/t5fCjxQy8sbY/4ZsagpCfv36YHrvvaAVftVV8PWvt4b8JZeoRS4SFgr8DGYRY8b9M6h9rZYtN2+h5H+VMOVvpmBZXSe0O6xbB1/7WtAHn5sLCxfCd74Dd9wBBQVJPAAR6RX9UZ3hsguzmf/qfMb9yTje+9v32L50O00NnYxBBDz/PFxzTRDsRUWwZg1UV8MTT8DddyvsRcIulC183dM2ubILspm1ehZD5g1hzzf34FHnsrWXEckN2gNVVUGLft264ETrT34Cy5YFV7+ISPoIZQtfX7xKPjOj+BvFTL9vOkf/8yjbbt9GY10Tq1fDZZcFrfjvfQ8qK+FLX1LYi6SjULbwJXUm/cUkLMfY/ee7+ddpW/nTg3P4yCeyePjh4OoaEUlfoWzhS2r5rRN5dPKlTDp4nEemb+XZJxoV9iIDgAJf2tm5M7i08lcnJ9DwP2Yy4p0TbF/yFo1nGlNdmojESYEvF+zcCdddF4w589JLcNOPxjNr9SxOPn+SNz75Buerz6e6RBGJgwJfACgvh2uvDa6z/93vgm/FAhR9roi56+dyZssZXv/I65zdfTaldYpI3ynwhV//Gj7+8WB8m+efD67KaWv0ktHMe24e0ZNRKq6s4NhTx1JSp4jEJ5SBrxugJE9lJdx+O8ydC6++CrNmdb5d4YcLmV82n0FTB/HWLW9xaNWh5BYqInELZeDrOvzk+bu/C66p37AhGK2yO4NKBnHFy1cw4oYR7Fq2i8q/rKT+UH33LxKR0Ahl4EtyrFkDv/oVfPWrwVAJscganMWc9XMo+nwRVT+u4pWpr/DOve8QPRXt32JFJG4K/Az17LPwhS/AJz4RtPJ7I7sgm1kPz+KqXVcx9o6x7P/f+3l1+qtUr68mzMNti2Q6jYefgU6cgClTgqGLf//7YMz6eJwqP8Xb//1tTr92mvySfEZ/ZjSjloyi8JpCItlqU4gkm8bDlwt++cvg/rAPPxx/2AMMKx3G/M3zOfzLw1T/ZzUHfnKAqvuqyB6ZzajFoxj16VEMWzCMvOK8hN5ZS0R6Ry38DOMeXJEzZAj84Q/98x7R01FObDzB0ceOcuyJY0RPBv372SOyKZhfwJC5Q8ifkk9+cT55k/PIHZdL9vBssoZmYRF9IIjESy18AeCxx2D7dli1qv/eI7sgmzG3jWHMbWNoamji9GunqX29ltOvn+b0a6c59LNDNJ3pZMx9C8bnzx6eTVZhFlmDs4jkR4jkR7BsC6Ysu7BMFu3XZQUT1rq/tvN2f10kYpu2n036nJIEK/piEYNKBiV0nwr8DNLYGNyZ6tJL4bOfTc57RnIiDLtqGMOuGnZhnbsTPRmlfl89dXvraKhuIFoTJXqyeWpebjrXRNO5JhpqG/Co443eft5mmUYuPA7epPW92j5utxzPNp1tK5JAI64fkb6Bb2azgb8ARgO/dfefJOu9JbBmDWzbBo88Atkp/Kg3M3JG5JAzIoeCebpNlkiyxHQJhZk9ZGZHzGxrh/WLzGyXmVWa2b3d7cPdd7j7l4A7gIv6lqR/1dfDd78LH/xg8M1aEck8sbbzVgH3A6tbVphZFvAAsBCoAsrMbAOQBazo8Pq73f2Imd0K3Nu8L0mif/on2LMHNm6EiK6UFMlIMQW+u79gZiUdVi8AKt19D4CZrQWWuPsK4JYu9rMB2GBmTwL/1teipXfeew++/3247Ta48cZUVyMiqRJPT+5EYH+bx1XAVV1tbGbXAn8E5AFPdbPdcmA5QHFxcRzlSYuvfQ3M4Ec/SnUlIpJK8QR+ZxeidXm9grs/Dzzf007dfaWZHQIW5+bmfqjP1QkA+/cHl2J+97sweXKqqxGRVIqnN7cKaBshk4CD8ZUT0GiZifOb3wTz225LbR0iknrxBH4ZMMPMpppZLnAnsCExZUmibNoUjIQ5d26qKxGRVIv1ssw1wGZgpplVmdkyd48CXwE2AjuAR919WyKK0g1QEiMaDVr4N9wQ9OGLSGaL9SqdpV2sf4puTsD2lbs/DjxeWlp6T6L3nUnWrYPqanXniEgglFdkq4UfP3f4h3+A2bPh1ltTXY2IhEEoA18nbeP35JPw1lvwV3+lL1qJSEBRMECtWBHc5OTOO1NdiYiERSgDX1068Tl4EP7rv+DLXw5uUC4iAiENfHXpxKflxiYf+1hq6xCRcAll4Et8Xn01GP74iitSXYmIhEkoA19dOvH5wx9g3jzIz091JSISJqEMfHXp9F1jI5SVwYIFqa5ERMImlIEvfVdeDrW18PGPp7oSEQkbBf4A88wzwTAKCxemuhIRCZtQBr768Ptu40a48koYNSrVlYhI2IQy8NWH3zcnTgRX6CxalOpKRCSMQhn40je/+Q00NcFNN6W6EhEJIwX+APLMMzB8uK7QEZHOKfAHCPeg//6GG4IvXYmIdKTAHyAqKuDAAfj0p1NdiYiEVSgDX1fp9N66dUHLXmPfi0hXQhn4ukqnd9yDwL/+ehg5MtXViEhYhTLwpXe2bIF33tGtDEWkewr8AWDduuCuVp/5TKorEZEwU+APAOvWwbXXwpgxqa5ERMJMgZ/mjh2DnTvh5ptTXYmIhJ0CP829/XYwnz07tXWISPiFMvB1WWbsdu0K5pdemto6RCT8Qhn4uiwzdm+/HVx/P3VqqisRkbALZeBL7N5+G6ZN03AKItIzBX6a27VL3TkiEhsFfhprbITKSgW+iMRGgZ/Gdu6Eujq4/PJUVyIi6UCBn8YqKoJ5aWlq6xCR9KDAT2MVFTBkCMycmepKRCQdKPDTWEUFfPCDkJWV6kpEJB0kNfDNbIiZVZjZLcl834GosRFefx0+9KFUVyIi6SKmwDezh8zsiJlt7bB+kZntMrNKM7s3hl39T+DRvhQq7e3cCWfPKvBFJHaxfl1nFXA/sLplhZllAQ8AC4EqoMzMNgBZwIoOr78buBzYDuTHV7KATtiKSO/FFPju/oKZlXRYvQCodPc9AGa2Flji7iuAi7pszOw6YAhwGXDOzJ5y96Y4as9oOmErIr0VzxfyJwL72zyuAq7qamN3/2sAM7sLONpV2JvZcmA5QHFxcRzlDWw6YSsivRXPSVvrZJ339CJ3X+XuT3Tz/Ep3L3X30jG6o0endMJWRPoinsCvAia3eTwJOBhfOQENj9y9lhO26r8Xkd6IJ/DLgBlmNtXMcoE7gQ2JKErDI3ev5YStWvgi0huxXpa5BtgMzDSzKjNb5u5R4CvARmAH8Ki7b0tEUWrhd08nbEWkL8y9x273lCktLfXy8vJUlxE611wTzF96KbV1iEg4mVmFu1/U6RvKoRXUwu+aTtiKSF+FMvDVh981nbAVkb4KZeBL13TCVkT6KpSBry6drumErYj0VSgDX106XdM3bEWkr0IZ+NI5nbAVkXgo8NPI9u3BCdsrr0x1JSKSjkIZ+OrD71xZWTBX4ItIX4Qy8NWH37myMhg2DGbMSHUlIpKOQhn40rmysuD6+4h+ayLSB4qONFFXB1u2qDtHRPoulIGvPvyLvfkmNDQo8EWk70IZ+OrDv5hO2IpIvEIZ+HKxsjIYNw4mT+55WxGRzijw00RZWdC6t85uLCkiEgMFfho4ejQYJfOqLm8RLyLSs1AGvk7atvfb34I73HBDqisRkXQWysDXSdv2Nm2CwkKNgS8i8Qll4MfrzTfhF79IdRWJ4R4E/ic/CdnZqa5GRNLZgAz8H/4QvvhF+PM/hx07Ul1NfHbvhn37YOHCVFciIuluQAb+gw/Cl78M//IvcNllMG8erFgBb78dtJjTyaZNwVyBLyLxGpCBP2QI3H8/7N8PP/4xFBTAt78d3CVq8mT4/Ofh4Ydh795UV9qzTZtg6lSYNi3VlYhIujMPcZO3tLTUy8vLE7KvvXvhmWfguefgd7+D6upg/dSpQf/4ddcF04QJCXm7hDh1CiZNgqVLg79WRERiYWYV7n7RZR6hDHwzWwwsnj59+j27d+9O+P7dYdu21vB//nk4eTJ4bvp0WLAg+JLTlVfCFVfA4MEJLyEmP/gBfOtbraNkiojEIq0Cv0UiW/jdaWwMrux57jl4+eUgYA8cCJ7LyoI5c4Lwv/nm4MNg4sT+H6K4ri7462Pu3NZ+fBGRWHQV+LrQjyDU588Ppm98I1h36FAQ/C3T+vXw858HzxUWBt0/H/hAEMhz5sCll0JOTuJqWr0a3n8ffvWrxO1TRDKbWvgxikZh8+agK6isDF56CSoroakpeD43Nwj+efPggx8M5vPmwYgRvX+vysrgA2X8eHj1VY2fIyK9oy6dflBXB7t2wVtvBTcnefPNYDp8uHWbiRODrpkpU6CkpHVeUgLFxZCX17rt+fNBt9LnPhc8fvppDYcsIr2nLp1+kJ/f2pJv6/DhIPjfeAO2bw+uEHr5ZVi7Njhf0Nb48TB8eHAi+b33gg+RSy+FJ58MTiCLiCSKAr8fjBsHN94YTG1Fo3DwYBDse/cG8/feg9raIPA/9Sm4+mpYtAiGDk1B4SIyoCnwkyg7O+jGKS5OdSUikokG5DdtRUTkYkkLfDO71sxeNLMHzezaZL2viIgEYgp8M3vIzI6Y2dYO6xeZ2S4zqzSze3vYjQOngXygqm/liohIX8Xah78KuB9Y3bLCzLKAB4CFBAFeZmYbgCxgRYfX3w286O6/N7NxwP8FPhtf6SIi0hsxBb67v2BmJR1WLwAq3X0PgJmtBZa4+wrglm52dwLI6+Z5ERHpB/FcpTMR2N/mcRXQ5W22zeyPgJuA4QR/LXS13XJgOUCxLmcREUmYeAK/sy/8d/m1XXdfD6zvaafuvhJYCcE3bftcnYiItBPPVTpVwOQ2jycBB+MrJ2Bmi81sZU1NTSJ2JyIi9GIsneY+/CfcfW7z42zgbeB64ABQBvyJu29LWHFm1UBf70s1GjiaqFrShI45M+iYB754j3eKu4/puDKmLh0zWwNcC4w2syrgu+7+czP7CrCR4MqchxIZ9gCdFRwrMyvvbPCggUzHnBl0zANffx1vrFfpLO1i/VPAUwmtSERE+oWGVhARyRADOfBXprqAFNAxZwYd88DXL8cb6hugiIhI4gzkFr6IiLSR9oHf0wBuZpZnZo80P/9qJ0NEpJ0YjvkvzWy7mW0xs9+a2ZRU1JlIsQ7UZ2a3m5mbWVpf0RHL8ZrZHc2/521m9m/JrjHRYvh3XWxmvzOz15v/bX8qFXUmUlcDU7Z53szsx80/ky1mNj+uN3T3tJ0ILgd9B7gEyAXeBC7rsM2XgQebl+8EHkl13Uk45uuAwc3Lf5YJx9y83VDgBeAVoDTVdffz73gG8Dowovnx2FTXnYRjXgn8WfPyZcB7qa47Acf9cWA+sLWL5z8FPE0wssHVwKvxvF+6t/AvDODm7ueBtcCSDtssAX7RvLwOuN7MOhsWIl30eMzu/jt3P9v88BWCb0Gns1h+zwDfA/4PUJfM4vpBLMd7D/CAu58AcPcjSa4x0WI5ZgeGNS8XkqBv9qeSu78AHO9mkyXAag+8Agw3s/F9fb90D/zOBnCb2NU27h4FaoBRSamuf8RyzG0tI2ghpLMej9nMrgAmu/sTySysn8TyO74UuNTMXjazV8xsUdKq6x+xHPPfAp9r/vLnU8BXk1NaSvX2/3u30v2etrEM4NarQd7SQMzHY2afA0qBT/RrRfvguzIAAAL0SURBVP2v22M2swjwI+CuZBXUz2L5HWcTdOtcS/AX3ItmNtfdT/Zzbf0llmNeCqxy938ysw8Dv2w+5qb+Ly9lEppf6d7Cj2UAtwvbNI//U0j3f0KFXUyD1pnZDcBfA7e6e32SausvPR3zUGAu8LyZvUfQ17khjU/cxvrv+tfu3uDu7wK7CD4A0lUsx7wMeBTA3TcT3D1vdFKqS52EDlKZ7oFfBswws6lmlktwUnZDh202AF9oXr4deM6bz4akqR6Publ7418Iwj7d+3ahh2N29xp3H+3uJe5eQnDe4lZ3L09NuXGL5d/1YwQn5zGz0QRdPHuSWmVixXLM+wgGa8TMZhMEfnVSq0y+DcDnm6/WuRqocfdDfd1ZWnfpuHu0swHczOzvgXJ33wD8nOBPv0qClv2dqas4fjEe8w+AAuDfm89P73P3W1NWdJxiPOYBI8bj3QjcaGbbgUbgm+5+LHVVxyfGY/468FMz+xpBt8Zdad5463RgSiAHwN0fJDhX8SmgEjgLfDGu90vzn5eIiMQo3bt0REQkRgp8EZEMocAXEckQCnwRkQyhwBcRyRAKfJE+MLMcM6toXnYz+2Wb57LNrNrMBsIwDzKAKPBF+uYa4L+al88Ac81sUPPjhcCBlFQl0g0FvkgbZlZiZjvN7GdmttXM/tXMbmgepGy3mS1o3nQR7Qelexr4dPPyUmBNMusWiYUCX+Ri04F/Bi4HZgF/QtCi/wbw7eZtrgOeb/OatcCdZpbf/LpXk1WsSKwU+CIXe9fd32oehXEb8Nvmr/C/BZSY2QTgeJt7DuDuW4ASgtb9UymoWaRHaT2Wjkg/aTu6aFObx00E/2duJhjzpaMNwA8JxkZJ53suyAClFr5I73Xsv2/xEPD37v5WkusRiYla+CK9kwXMcPedHZ9w9yqCvn+RUNJomSK9YGbXAJ9z9y+luhaR3lLgi4hkCPXhi4hkCAW+iEiGUOCLiGQIBb6ISIZQ4IuIZAgFvohIhlDgi4hkiP8PoxV7p8xfjjUAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# Plot abundances versus fractional mass\n", "\n", "X = prof_data['h1']\n", "Y = prof_data['he4']\n", "Z_C12 = prof_data['c12']\n", "Z_N14 = prof_data['n14']\n", "\n", "m = prof_data['mass']*MSUN # Convert masses to g\n", "M = prof_data['star_mass']*MSUN # Same\n", "\n", "plt.figure()\n", "\n", "plt.plot(m/M, X, color='r', label='X')\n", "plt.plot(m/M, Y, color='g', label='Y')\n", "plt.plot(m/M, Z_C12, color='b', label='Z_C12')\n", "plt.plot(m/M, Z_N14, color='m', label='Z_N14')\n", "\n", "plt.xlabel('m/M')\n", "\n", "plt.yscale('log')\n" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "# Define constants\n", "\n", "G = 6.674E-8\n", "\n", "# Evaluate lhs and rhs of the hydrostatic equilibrium equation\n", "\n", "lhs = np.gradient(P, r)\n", "\n", "rhs = -G*m/r**2 * rho" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAERCAYAAAB4jRxOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3deXxU9fX/8dfJToAsQBAEAsgiILKJgKIsIouoUEVbtCp+3Wpt3bXqD1vXVuu+V6kLKlpttSjUBVBU3BBRFlllhwCyJQFCErKd3x93JhsTMmEmc2cy5/l4xDv3zidzzwV855PPvfdzRVUxxhjT8MW4XYAxxpjQsMA3xpgoYYFvjDFRwgLfGGOihAW+McZECQt8Y4yJEmEf+CLysojsFJFlfrR9XEQWe75+FpHcUNRojDGRQML9OnwRGQLkAa+pas86fN+1QF9VvazeijPGmAgS9j18VZ0HZFfeJiKdRORjEflBRL4UkW4+vvUC4F8hKdIYYyJAnNsFHKEpwNWqukZEBgLPAad53xSR9kBHYK5L9RljTNiJuMAXkSbAycB/RMS7ObFas4nAO6paGsrajDEmnEVc4OMMQ+Wqap/DtJkI/CFE9RhjTEQI+zH86lR1H7BBRM4HEEdv7/siciyQDnzrUonGGBOWwj7wReRfOOF9rIhkicjlwG+By0VkCbAcGF/pWy4A3tJwv/zIGGNCLOwvyzTGGBMcYd/DN8YYExxhfdK2RYsW2qFDB7fLMMaYiPHDDz/sVtUMX++FdeB36NCBhQsXul2GMcZEDBHZVNN7NqRjjDFRwgLfGGOihAW+McZECQt8Y4yJEhb4xhgTJYIS+CIyRkRWi8haEbndx/uJIvK25/3vRKRDMPZrjDHGfwEHvojEAs8CZwA9gAtEpEe1ZpcDOaraGXgc+Hug+zXGGFM3wbgOfwCwVlXXA4jIWzhz26yo1GY8cLfn9TvAMyIiNt+NiTZappQUFFN8sIzi2CRKdudSsi+fksISDqa2pLiwFNmwnrKiEkoSG1PQriuJ61Yg+3LRkjL2t++JlimNf16ElpZR1CiV3C4nkrJiPnG5u9CSMnYfO5iyUqXFss/R0jIKm7Rge48RHLVkNonZ26G0lC09z0AVMn/6AFWlIKUVm3qdTeaSmSTnbkOBzcefRZkKHZe+j6qSn3o063v9io5LptM4ZysA647/FYrQecm7oEpeahvW9D6PzovfoUnuFgB+Pv48ALou/jeosj+tLat6T6Tb4rdomuNcMr6y10QUocfiN0GVfantWNbnIo5bPI3U7I0A/NT7IhDh+MWvA7A3NZOf+l7C8YtfJzV3E6oVbXotfq28zdK+k+i1+LXyNkt7XwxA7yVOm9y09iztM4lei18l1VPPkt6XANBnyatOm9T2LOl7Kb0Xv0pqzsaKNiL0WTy1vM3ivv9H70VTSct12izuPQlE6Lv4lfI2i/peRp9Fr5S3WdT7UhCh3+KXy9v8fOrl/OlPQfknV5WqBvQFnAe8WGn9YuCZam2WAW0rra8DWtTweVcBC4GFmZmZakwolJWU6oHd+bpli+r2r9epfvONln74seau260z/7FFpw99TD8dcre+O/oFvf561XcGPqRfdfytftvmXL1myE/6uxN/0J+aDNTVSb10esrFmpmp+kHirzRb0nUvTXVEwhd6WsxnqqAK+hlDFVQ/Y2j5tmHM1aFUtPmcIQqqcxlWpc0QPrc2DajNcObqUD7XUkRLEf2Modqq1ZH/WwYWag15HfDkaZ5piker6hWe9YuBAap6baU2yz1tsjzr6zxt9hzus/v37692p605UmWlSvbqXWQv2UL2L0WsTDuJjJkvk7JmIXE5u3ix5WTyc4t4etPZpJftYT6DOJWvmM1IRvIJAKOYRR5N+IbBAHwhQxnX9HOmHryAfqULKIptxF/b/5PS5KbcsPkmSuIbkdWsFzMH3McZPz9Jq7y1aFw8C/pdDYmJnLh8KsTFkdeiA2sG/JauK9+naf4vSEI8O/qMITYpntZr5hETH0txWgZ7ew4m/efvSMjPReJiyevaj9j4GJpuWIrExlDWNJXCLseTvHkVcYV5SFwsJR06I7ExJP2yEYmNgeRkytq1J37XNmJLDiKxMWjLo0CEuJxdzh9WYiJlzTOIydmDFBchApreDESI3ZfjtElIQNPSidmXCyUliAApKSBCzIH9Tpv4eGjaFDmQB6WlTpvGjREBKSwAQOJiITkZOVgIZWVIjEBCgtOmpNhpExvjbPOuxwjExDjL0krPNYqLc9Y9OSaxMSACZWUVbWJjq3yPxHpGsitnX0xM+feIeP9TjUjFfqTa9jAiIj+oan+f7wUh8E8C7lbV0Z71OwBU9YFKbWZ52nwrInHAL0CG1rJzC3xTm5KCYjZ8t5NV+9sQ+5+3SPvhE5r8so6HG9/Nuu3JfFMyAIAl9KIPS3iHCQyTL9gb34Knuv2DgqM7ceHGv1GS3oLCtl3YPmoSqcu+ZvPy/TRpk0JJlx507duY0wblE5PSxAkrY8LY4QI/GGP43wNdRKQjsBXnaVMXVmszA5iEM6/9ecDc2sL+SBUUwObWAzjQbwj95j5SH7swLtm/cQ9r3v6RH3KOYcWaeK7/YCRtDq4nnjaMYyMvMocTZSZbk7vQuUcxx5zblS+yniL+mHY06tGR9UOh1VHv0ChZaA48Uf7Jz1fb02Afe0+tz0MzJiQCDnxVLRGRPwKzgFjgZVVdLiL34owlzQBeAl4XkbVANs4PhXqhCsV78/nlm/X1tQsTCqps+Xwd699awBt6If3fu5Ordv2VfsB/uY/pR9/GuWm9WNfxPBJ7dWP+ZdAp83mat3qJDIGK519eW+2Dw+vXb2NCKawfgHKkQzo/HXU6JfsO0LfAnnIYCX74AVatglED97Lss128t/QYbnihGx2L1wDQo+kWzuq2lpFp39Ps9H4cc/4JpHdMc7lqY8JTfQ/phJ13xr/OP99IZquG3fkUU81NNypznljG49xIGl9QzHD+2Wg2J3S4hA19WtDuoqEsO7MNMbFtgWFul2tMRGuQgd+x8U7Oy/+C/fuuJSXVEj9clJbCk0/CS88W8tFV0yl69Ck27bqFMy4YyEnfbOP7djeRNmE82VdDUtKdbpdrTIPTIAP/+N1zuZSb+HnVxaQMTHe7nKjmvTouJgZeuG45Nz/Xg7/wdzJvv5s1dKb3gET+32ttiYtbwcnulmpMg9cgJ09r1KEVANnLt7tcSXT79lto1w6e/t0yNvc+m2ue68m9p8xhx9jLGMUsjmU1I584i7gG2e0wJvw0yP/VGg86nme5hta5yW6XErWysuBPw79n28H+NH37RZrnz+OBJn/lxn8N4JOFaZzzYTsA+vs8tWSMqQ8NsofffGhP/siz/FzUwe1SotKmD5eztfdYvjw4gH9MnMdt+yeTWbqejRf+P5q0TeOMM2DMGPjoI+fGTGNMaDTIwG8iB9glGRz7wWNulxJVCvYW8Ze/wNyzH6Nr9rd8O+ERxvxlALvJIJvmDHBueiUx0Qn7MWPcrdeYaNMgA5/kZBpzgLid29yupEErLKx4vfi+mWS36MJn933J/HMfpmjZGk5652Y6dG9U3mbCBBeKNMaUa5Bj+IiQndCKxvvtpG19mTcPhg6FDz6AA39/hvPnXcvqhJ48+kwSA/7QrErbL75wfjik2b1SxriqYQY+cN3ABew4mMZpbhfSwJSUwJlnwuzZcDJf86czU8nmXJoO2MXQ2ZNplJpwyPcMGeJCocaYQzTYwO+lS5CNq4Fr3C6lwdi3D95/H2bPVu7mbv7MfXzS6Gzk/fcZOfIet8szxtSiwQb+qXveo9/OaVjgB8crr8BllzmXUbZjC7fwCF93uoSRPz6NpLhdnTHGHw3zpC0Q2641aZrLx9ML3C4lohUUwO7d8LvfQXs2csnCa9lOawbyHW+PfgVJaep2icYYPzXYwO9x1an8jTv4cUGJ26VEtO7dISMDjitexLecxEVMowtrWE5PmqbYPEXGRJIGO6STce6pPNPqFM7cGb7TP4ez4mKnZ7/JeaYzZzOTYuI5nU/YndEDdkGTJu7WaIypmwbbw2f5ctbvSCbzx/fcriTiFBZCQgKccw4M4Qsu5jXu48/0ZRErOI7mzZ12FvjGRJaGG/jNm5Okhcgvdi1+Xb3n+RnZ9Ls5zGYUt/AIcZSQjZP0JZ5RMm/wG2MiQ8MN/IwMyiSGuJxdblcSMQ4ccMJ8/nwA5R7uYjXHMozPKaFi0ptLL3WWg309+tUYE7YabuDHxvLkvfu44+Dd7N/vdjHhr6TEGaK54gpovGM9bdjKWD5kfOIscqh65+yNN8LBg9Cxo0vFGmOOSMMNfGDg9ve4jJfYutXtSsLfF184yw9f3cm1H4xmJmezjxTi2rU+pG1ysjPGb4yJLA32Kh2AY398k2v4haysy+nWze1qwtPBg/Dii7B9O4DyLhNI2b+Vc/gUJYa2bWHtWqftu+9Cdrab1RpjAtGgAz8+szWt5y/ipyy3KwlfU6bAdddBQrxz+eo93EVjDjCfkwBo27ai7bnnulGhMSZYGvSQTtLEc3iCG8iywK+Rd2jm2uJHeSv+ElLPGsIMxpe/7w38Tp1cKM4YE1QNOvATzjmTN1tcz9bNpW6XErZiYmAMH/EQfyI5ppBjj6v6S1+bNs5vAZ995lKBxpigadCBz4wZZO1OIm7lT25XErby8uB0PmEJvflTy6k0z6j6T6JZM7jySudh5MaYyNagx/DJyACgbOt2oI+7tYSj/ftp9tOPXMojNGU/rRs1PuRmqmbNfH+rMSbyBNTDF5FmIjJHRNZ4luk1tCsVkcWerxmB7LNOWremVGIp2ZUTsl1GlJtu4qKpI8hkM/tJISkJWrRw3kpPh2uugREj3C3RGBM8gQ7p3A58qqpdgE89674UqGofz9e4APfpv/bteejeg0zJu5ACmyW5iuxpH8KLL/Jpn1vYTPvy7Y0bO8sePeDZZyE+voYPMMZEnEADfzzwquf1q8CvAvy84BJh1LLH+A1v2c1XlRw4AGf9uS9vtriOt7pXPKmqoABatXJejxrlUnHGmHoTaOAfparOLTvOsmUN7ZJEZKGIzBeRw/5QEJGrPG0X7toV+Dw43b59hfN4xy7NrOTL4X8mf+MObk14ktyCxPLtBQXO/PfLl8Odd7pYoDGmXtR60lZEPgFa+Xhrch32k6mq20TkGGCuiPykqut8NVTVKcAUgP79+wc8mb20bkXrzdvZEMWBv2kTZGaCCCy9fwZjvr+fBfEJPLS3D3l5Fe2Kipxljx7u1GmMqV+19vBV9XRV7enj631gh4i0BvAsd9bwGds8y/XA50DfoB1BLWL+8Hue45qo7eFv2gQdOsA990BBXikp997MqoTjKbrhNg4cgJwcSE112iYluVqqMaaeBTqkMwOY5Hk9CXi/egMRSReRRM/rFsBgYEWA+/Vb0m8nsCB1FFlbovPJV3v3Ostp0+D+v8VwevFH7H32DTLaOLfYbtsGLT0DcY0auVSkMSYkAg38B4GRIrIGGOlZR0T6i8iLnjbdgYUisgT4DHhQVUMW+Dz+OGv2tmTPxuibI/m222D0aOf10evmMeqB4Zx1XiMGXnE8aWnO9m3boFcviI2Fv/3NvVqNMfUvoBuvVHUPcMiV2qq6ELjC8/ob4PhA9hMQz2UnRZu2AymuleGGhx5ylgkcZApXkRRTxJ2POEnvHcYBZ377d95xoUBjTEg17KkVAFq3Jj8hlfztuW5XEjKq8NZbFeun8BWdWcu6G5+lRXvnQvvKgW9PrjImOjT8wB8+nEcm5/Jx9sDyq1AaurfeggsucF4ncJC5jOC6ses47eEzytt4Az/d573RxpiGqGHPpQNQVsa5X93CQoazbds4OnRwu6D6t2NHxeuXuYxEDtLjof8gUrE9xTO6ZVfmGBM9Gn4PPzaWbl+/xAg+ZfNmt4sJrVP4kt/yJivpTnJjqfKet2c/caILhRljXNHwAx8oa9WGNmxl/Xq3KwkN9VyBOpDv2EQmD3DHIZdcNm8OWVnw8MOhr88Y446oCPyYe+7iJbkyqgI/nWwe5RZ6sZQCkn1eY9+mjXM5pjEmOkRF4Med9yv2Hd0tagI/rriAJfTmr/w/9uGcnbWbqowxURH43HMPX2zrzIZ1ZW5XEhK9v36OdmQxi9Hl22yaY2NMdAR+mzbEaQl71wY++2YkaJq9iY8YwzyGul2KMSaMREfgZ2ayN7UdxbtzOXDA7WLq1/IvdvPEMU8xjtA9WMwYExmiI/DHjePjFzbzM8c26HH8vM3ZtBnWmZavP0IJNoZjjKkqOgI/L4/T3ricMXzEOp+z8DcMSy95mBT28TFj3C7FGBOGoiPwk5Jo8b9XGMR81q51u5jgeewx+PRT5/W+fbD52628yYUsp6e7hRljwlLDn1oBIC4OOeooOuVsZd7PbhcTPDff7CxV4eVHsrmx6DWO7VQCDfi3GGPMkYuOHj7A00/zRdcr+bkBBb7XgfU7uPK+TJ7q8TyDh0bHz3BjTN1FT+CfdhoZnVIaROCvXQtvv12xvuLKx0migJPvPI3Gjd2ryxgT3qIn8CdP5s7Zp7J9O1Ue3B2JevasOunZtm828kXG+ZxwQVefgW8zYhpjIJoCv00bkvP3kEgha9a4XUxgDh6seN2IfH5V+Ba89hoAycmHtm/TJkSFGWPCWvQEfufOFHY6jjRyG8SwDkBj8ljPMTzW4SmGj3YeSl69h//HP8KcOS4UZ4wJO9ET+BMnUrZ0GTto1WAC/3e8QCt2cPINA8ofblI98P/6V+jYMfS1GWPCT/QE/o4dJE86n4kZn0b8kI5XL5aysdNpDLx+UPm26kM6TZuGuChjTNiKnsBv1AjeeYdhqYsaSA9fuZSpfH3ju1W2Vu/hS9UHXRljolj0BH5KCjRtStekzRF/t61Qxo/048/ch6SnVXnPG/hdu8JPP7lQnDEmbEVP4AO88QbLBv+O7Gwoi+Cp8c/gI/qymPUcQ2Ji1fdiPH+jbdo4l28aY4xXdAV+v36kJxWgGlnX4r/0EixYULE+itlspD1v85tDAr+42FkmJISuPmNMZAgo8EXkfBFZLiJlItL/MO3GiMhqEVkrIrcHss+A3HMPE14+E4C9e12ros6uuAIGDvSsqHIDTzCI+ZQQf0iwd+3qLM8/P6QlGmMiQKA9/GXAucC8mhqISCzwLHAG0AO4QER6BLjfI9OuHY327ySBgxEV+FWMH8+9/IUdtAI4pIffubMzc+bll7tQmzEmrAUU+Kq6UlVX19JsALBWVderahHwFjA+kP0esR49yOlxMqnsJTfXlQoCUrR0FcycSQEVTySvHvhgl2IaY3wLxRh+G2BLpfUszzafROQqEVkoIgt37QryM2gnTODnl79mFy0jsoe/+O73KCKel6jovvsKfGOM8aXWwBeRT0RkmY8vf3vpvq4E15oaq+oUVe2vqv0zMjL83IWfdu7kuJtHcxYzIybwS0srXl+x5jaOYzk7Oap8m52cNcb4q9bJ01X19AD3kQW0q7TeFtgW4GcemSZNaPL1bI5nKHv3nu1KCf6aPRuys+FXv3LWJ/IvBi2bz+08WKWd9fCNMf4KxdMyvge6iEhHYCswEbgwBPs9VHIy2rw57fZsCfse/ujRztJb5808SiMppFCrznVsgW+M8Vegl2WeIyJZwEnAByIyy7P9aBH5EEBVS4A/ArOAlcC/VXV5YGUH4J13eTLulrAPfK+iIujLj/TnB5ae/Huqj5BZ4Btj/BVQD19VpwPTfWzfBoyttP4h8GEg+woWaduG7smbyM3t5HYpfikqgsX04XTm8MoLA6j+fHILfGOMv6LrTluAxx9nat6EiOnhl63bwMPcymXPnEDbHimHvG8nbY0x/oq+wO/QgdSyXIp3R0biN3n+Ea7hORqVHUDk0NkvrYdvjPFX9AV+795812wMRTkH3K7kEPv3w/r1FesxlJI8dyazGE1p67bOtmp/Y9XXjTGmJtEXF6NG8eCQj1hfeLTblRxiyBDoVOnUQhkxbJ48hSe5nvh4Z5sFvDHmSEVffOTn8/DXJzFm60tuV3KIxYurrl/Im+RldORzhpeP1VvgG2OOVPTFR6NGtM1dxjF5S92u5LBas41XmUSzD6cBFSdn7QlWxpgjFX2BL0JuekfaFG+oMm1BOFGFy3mJOErZctokAOvhG2MCFpXxMeuSN/k9/2D/frcr8a24GOYxhEe4md1pnQELfGNM4KIyPhqlJzGcz9ibW+Mcbq4qWbqC9RzDrTzCAc/FRHa9vTEmUFEZ+Mf+PJNpXEze5my3S/Ep/k83MJfTAA4JfO8jDI0xpq6iMvC1Q0cAin7e6G4hPgxkPvGfzeFVnLH7mgL/xRdh82YXCjTGRKyoDPyYnj14m1+zrzD8xkkSOUj+2b/mCW4ADg38sjJnOWwYtGt36PcbY0xNojLwk3p1ZSJvsyXteLdLqUIoYyXd2fro2xygCVAR+N4br7xSDp1WxxhjDisqAz81Fb7kFI6beqvbpVQxgk/Joi3y1Zfl22o6aWvPrTXG1FXUBn4y+TTdssLtUqq4hUfYSyq7O55Yvq164Hfr5ixt0jRjTF1FZeAnJcEm6UjK7vW1Nw6VvDySyecB7iCvpOKpVtUDf948+OYbu+PWGFN3URn4APekP8Hfxs0H4KuvYORIly95bNSIIczjSa4nP79ic/Ux/IwMOOmk0JdnjIl8URv4jVNi6bP0NfjlF6ZNg08+gV9+qf/9qsIf/uD01KtsPPFE7uR+yogtD3lwAj821vkyxphARG3gd0nczKU/Xgfff8833zjbcnLqf7+lpfDcczB0aMW24nnfwqJFbMOZsrl64NtdtsaYYIjawM9u6Zz9LPhxJcuWOdtCEfi+ho2+uXs2+2nC2/wG4JAhHQt8Y0wwRG3gx2ek8Xqz61kR1wv1TKmTm1v/+60e+GVlMGnD3XRj1SHX3gPk5VngG2OCI2oDPyMDri19ghl5p5VvC0UPv6Sk6vryO6Zx+6aryaZZ+bbKPfyNG6FDh/qvyxjT8EVt4A8ZArfsvZMb/t6KDu2dLr4bQzqNXnyawXxNIYdeiuk1aVL912WMafiiNvBHjoQdHEW65vCbob8gEvrAz/16OZ2zFzAj4wqg4sL6yj18gLFj678uY0zDF7WBn5EBpZ27UUoMI47ZQGpq6MfwX/2+B6fwJenXXVylTfXA79ix/usyxjR8URv4AM0nDCOZfLpeejLp6aHt4SdSSOJ9k2nRuy3HD21WpU31IR1jjAmGgAJfRM4XkeUiUiYi/Q/TbqOI/CQii0VkYSD7DKabb4/nxyufp/2X00Ie+OcwnauzH+CK4esOmReneg/fGGOCIS7A718GnAu84Efb4aq6O8D9BVVaGqQtfgM2pJKWdlFIh3Su4EU2SEeG3DWcDZuqtrEevjGmPgTUw1fVlaq6OljFuKJbN1i1KuQ9/Ie5lQ+GPkxKWswh19lbD98YUx8C7eH7S4HZIqLAC6o6paaGInIVcBVAZmZm/Vd28cVw8smkL1Rycup/CsriYhjOXFbQgwcfd46v+pDO7rD6PcgY01DUGvgi8gnQysdbk1X1fT/3M1hVt4lIS2COiKxS1Xm+Gnp+GEwB6N+/v/r5+Udu5Ejo14+WKw+Qk9Ok3ndXklfINC4iYdAJtOgzEzg08Ddtgk6dYN26ei/HGBNFah3SUdXTVbWnjy9/wx5V3eZZ7gSmAwOOvOQgW70aWrRg4LbpFBZCYWH97q71v5/gaLaz7Tc3lm/z9TCTiROdZVpa/dZjjIke9X5Zpog0FpGm3tfAKJyTveHhmGMgLo42+1cB9XMt/pdfwn//67ze3/pYbuNBDgysmNLBV+CfcQZMnQoLw+aaJmNMpAtoDF9EzgGeBjKAD0RksaqOFpGjgRdVdSxwFDBdnEc0xQFvqurHAdYdPPHx0KsXjWMKAOfEbStfA1gBGDLEWequ3ezpNICHOIfzKz2U3NfkaCeeCIMHB7cOY0x0CyjwVXU6zhBN9e3bgLGe1+uB3oHsp94tXMjGWQIfBd7D370bfvgBRo921jdurHiv8KGnGPH4g7RgK/HxGeXbfQW+zZBpjAm2qL7TttycOQx46DziKA740sxx42DMmIpr6Wc652WJo5iif7zIpq4j2U1G+SMLwZ5Pa4wJDQt8gD17aPbZuxzH8oAD33tlzd69zjIry1m2ZCcL87px3YqrAaoEvjHGhIIFPsDAgQAMYEHAgd/Ec2Wn93Py8yE1FfaRwgjmMiv+bMAC3xgTehb4AB07UvLSq3zEGQGP4VcP/IICOD5pDbtjWjKe92jUyNlugW+MCTULfAAR4s4+g95JP9dLD/+yg/8gMaaEVSkDyctztlvgG2NCzQLf69lnmVE4koKd+wP6mMaNnWV5Dz9f6VayDH7zG3IbtaaszNlugW+MCTULfK+BA4lBab7xh4A+pnoPvzC/jJuOmwVTplS51LKmwBexJ1wZY+pHqCZPC38nnsjSlMEc2F8W0MckJzvLnByguJgpX/bg3cwbIfmaKoEfV8OffE6Oc5LXGGOCzQLfq0UL7hz6FZs3B/Yx6pnuLTcX+N//aFe4lr0pzqyY/vTwvUNCxhgTbDakU8mFWx9m2op+AX1GSYmzzMkB5s5lR9zRrMgcA1QEvgjExvr+/pp6/sYYEygL/EoSm8TTs3gRbNt2xJ/hfcBJTg7w1FOc2+pbkpo4Ke6dJM1O2Bpj3GCBX8mezs4NWKXffX/En+EN/BEzr2f1xLtYW5RZfu29t4dvgW+McYMFfiUFPU6gK6vJPcW5G7aoCD76qG6fUVwMHdjANTzHx//eS0FBxYncwwV+9+4BFG6MMX6wwK8kNSOBBIqIvf1WKC3llVecSyTr8uSp4mKYxKuUEcND/In8fP8C/8cfYX9gtwAYY8xhWeBXkp4O3VhF2suPwaefMs/zEMZdu/z/jOIi5V7+wgAWsI02lJbi15BOUlLFNfzGGFMfLPArSUuDmZxNcZM0eO01vvrK2e73dAt5eTyw7Cxe7vQ3bnmt4hEA3h6+nbQ1xrjJAr+S9HQoIpEfLnuWXyb8ofyafL8nVHvuOQbnfsj2lGNp3bpisx24N+EAABEeSURBVJ20NcaEAwv8StLTneXS4y7kxzVNOZEFQB0Cf+VKViX15rvM86sEvj9j+MYYU98s8CtJS3OWOdnKCQ9M4NGYW511f4Z0tm+Hl1/mt+2/Jj6eKoHvfUauN/Dt5ipjjBss8CtJTnaGX5YtF/6deAmnls2jW+KG2nv4Bw5AZiY88AD7ShsTH1/x2wJA//7O0nr4xhg3WeBXIgLXXAPTpsHDOy5mxXHn0yKlqPbA//hjZ06Ffv0oLnYCvfJzaps1c5YW+MYYN9ngQjX33AP//S9s2JDJjqf/TfxVW8jJVuDQJ43fcIPTsb8pZgv06gXDh5cHPsCFF1ad+dKu0jHGuMl6+NU0bgyvvw7jx8PgA7OZuzaTU5a/4LPtW/9Sfv7Hp3D99fD995CYWCXw33gDnnuuor2N4Rtj3GSB78PgwfDee5Aw9nQWZozh6tU3wOLFVdoc3LKTK3fez/NrT6f0vRnlaV5SUnMP3hv43imUjTEmlCzwDycmhhdPfY1PksdDRgaceio8/TSUlZHz/jzu4y/MYhTrup9d/i2Ve/jVeQO/LLBnrBhjzBEJKPBF5GERWSUiS0Vkuoik1dBujIisFpG1InJ7IPsMtdhWGUxKehvatHEG7K+7DkaNYl27YZzAQs7kA35aXvHHWFxc85CNDekYY9wUaA9/DtBTVXsBPwN3VG8gIrHAs8AZQA/gAhHpEeB+QyY93bnxShXn8p2pU+G449iwJ4UfOYFS4vjpp4r2h+vhe0/aWuAbY9wQUOCr6mxV9TzjiflAWx/NBgBrVXW9qhYBbwHjA9lvKKWlQWkp5OXhXGs5aRI8+SSbtjvd9bZtKQ/8sjLnq7YhHbtKxxjjhmCO4V8G+Jo9vg2wpdJ6lmebTyJylYgsFJGFu+oyTWU98d59W/1a/C1boEULOPFEWLbM2eZ9+EltgW89fGOMG2oNfBH5RESW+fgaX6nNZKAEeMPXR/jYVuN1Kqo6RVX7q2r/jIwMf46hXnnvmK0e+Js3Q7t2cMwxkJXlbKst8L3ba3qerTHG1Kda+5qqevrh3heRScBZwAhVnxccZgHtKq23BY78obEhVj6/TrX5dLZsgU6dnF67N+hrC3zv1TnWwzfGuCHQq3TGALcB41Q1v4Zm3wNdRKSjiCQAE4EZgew3lA43pNOunRPuxcXOSd3aAr+01Fla4Btj3BDoGP4zQFNgjogsFpHnAUTkaBH5EMBzUvePwCxgJfBvVV0e4H5DxteQzr59sHevE/jecfmSktoDv8RzetsC3xjjhoCiR1U717B9GzC20vqHwIeB7MstvoZ0tnhOQbdrV/G6uNgC3xgT3uxO21p4Jz9bsQJOOQUWLaoI+czMih5+UVHtgd+ypbPs0qX+6jXGmJpYX7MWsbGQkgL/+x9s2wYTJsC4cc57mZkVU+z408M/80z44AMYPbr+6zbGmOqsh++HtDQn7MHp3T/5pBP83pO24F8PXwTGjrXLMo0x7rDA94N3HP/oo+GFF2DkSHjpJWebd0jHnx6+Mca4yQLfD94rdbp0gcsug9mzK8b2veFeXFxxUtYC3xgTjizw/eDt4fs62VqXk7bGGOMmC3w/HC7wK/fwLfCNMeHMAt8PlYd0qrPAN8ZECgt8P3h7+J193Gbma0jHbqwyxoQjC3w/9O3rhL318I0xkcwC3w/jxsGaNZCUdOh7dbkO3xhj3GSBHyC7Dt8YEyks8ANkQzrGmEhhgR8guw7fGBMpLPADVLmHv3atc4WO9y5cY4wJJxb4Aap80nbmTDjtNGja1N2ajDHGFwv8AFU+abt7t/OcW2OMCUcW+AGqPKSTlweNG7tbjzHG1MQCP0DewC8shIICC3xjTPiywA+Qd0hn715n2aSJe7UYY8zhWOAHyNvD9z7k3Hr4xphwZYEfIO9EaRb4xphwZ4EfIBGnl+8NfBvSMcaEKwv8IIiPh9xc57X18I0x4coCPwgSEmxIxxgT/izwg8CGdIwxkSCgZzOJyMPA2UARsA74P1XN9dFuI7AfKAVKVLV/IPsNNwkJsGeP89p6+MaEt+LiYrKysigsLHS7lIAkJSXRtm1b4uswW2OgD+ObA9yhqiUi8nfgDuC2GtoOV9XdAe4vLMXHOzdegQW+MeEuKyuLpk2b0qFDB0TE7XKOiKqyZ88esrKy6Nixo9/fF9CQjqrOVtUSz+p8oG0gnxepKv+AtSEdY8JbYWEhzZs3j9iwBxARmjdvXuffUoI5hn8Z8FEN7ykwW0R+EJGrDvchInKViCwUkYW7du0KYnn1x3u3LVgP35hIEMlh73Ukx1DrkI6IfAK08vHWZFV939NmMlACvFHDxwxW1W0i0hKYIyKrVHWer4aqOgWYAtC/f3/14xhc5+3hJyVBbKy7tRhjTE1q7eGr6umq2tPHlzfsJwFnAb9VVZ8BrarbPMudwHRgQPAOwX3ewLfevTHGH008Y7+ff/45Z511Vsj2G9CQjoiMwTlJO05V82to01hEmnpfA6OAZYHsN9x4h3Qs8I0x4SzQq3SeARJxhmkA5qvq1SJyNPCiqo4FjgKme96PA95U1Y8D3G9Y8fbw7YStMZHlhhtg8eLgfmafPvDEE/63z8vL47zzzmPZsmWccMIJTJs2DRHh9ttvZ8aMGcTFxTFq1CgeeeSRgGsLKPBVtXMN27cBYz2v1wO9A9lPuLMevjHmSC1atIjly5dz9NFHM3jwYL7++mt69OjB9OnTWbVqFSJCbu4htzcdkUB7+AYbwzcmUtWlJ15fBgwYQNu2zhXtffr0YePGjQwaNIikpCSuuOIKzjzzzKCN89vUCkFgQzrGmCOVmJhY/jo2NpaSkhLi4uJYsGABEyZM4L333mPMmDFB2Zf18IPAO6TTqJG7dRhjGoa8vDzy8/MZO3YsgwYNonNnn6PndWaBHwSVr8M3xphA7d+/n/Hjx1NYWIiq8vjjjwflcy3wg8Dbw6/0m5kxxtQoLy8PgGHDhjFs2LDy7c8880z56wULFgR9vzaGHwTWwzfGRAIL/CDwBr718I0x4cwCPwi8QzrWwzfGhDML/CCI8fwpWg/fGBPOLPCDoLTUWVrgG2PCmQV+EJSVOcs6PGnMGGNCzgI/CLyBb3PhG2OORJMQ3aZvgR8E3sCPsT9NY0wdqSpl3hCpZxZRQeB97IsFvjERaNiwii+AqVMr1qdOrXsbP2zcuJHu3btzzTXX0K9fPwoKCpg8eTK9e/dm0KBB7NixA4D//Oc/9OzZk969ezNkyJBAjhKwwA8K70lbC3xjjL9Wr17NJZdcwqJFiwAYNGgQS5YsYciQIfzzn/8E4N5772XWrFksWbKEGTNmBLxPm1ohCGxIx5gI9vnnVdcvvdT5CrRNLdq3b8+gQYMASEhIKJ8C+YQTTmDOnDkADB48mEsvvZRf//rXnHvuuXX6fF8sooLAAt8YU1eNKz1AIz4+Hs9TAcunSAZ4/vnnuf/++9myZQt9+vRhz549Ae3TIioILPCNMfVh3bp1DBw4kHvvvZcWLVqwZcuWgD7PhnSCwHvDlXeKBWOMCYZbb72VNWvWoKqMGDGC3r0De1qsqPcSkzDUv39/Xbhwodtl1ConBx58EO6/326+MibcrVy5ku7du7tdRlD4OhYR+UFV+/tqbz38IEhPh7//3e0qjDHm8GzU2RhjooQFvjEm6oTzULa/juQYLPCNMVElKSmJPXv2RHToqyp79uwhqY4P4bAxfGNMVGnbti1ZWVns2rXL7VICkpSURNu2bev0PQEHvojcB4wHyoCdwKWqus1Hu0nAnZ7V+1X11UD3bYwxdRUfH0/Hjh3dLsMVwRjSeVhVe6lqH+B/wF+qNxCRZsBdwEBgAHCXiKQHYd/GGGP8FHDgq+q+SquNAV8DY6OBOaqarao5wBxgTKD7NsYY47+gjOGLyF+BS4C9wHAfTdoAle8JzvJs8/VZVwFXAWRmZgajPGOMMfh5p62IfAK08vHWZFV9v1K7O4AkVb2r2vffCiSq6v2e9T8D+ar6aC373QVsqrVA31oAu4/weyNRtB0v2DFHCzvmummvqhm+3vCrh6+qp/u5ozeBD3DG6yvLAoZVWm8LfO7Hfn0W7Q8RWVjT7cUNUbQdL9gxRws75uAJeAxfRLpUWh0HrPLRbBYwSkTSPSdrR3m2GWOMCZFgjOE/KCLH4lyWuQm4GkBE+gNXq+oVqprtuXzze8/33Kuq2UHYtzHGGD8FHPiqOqGG7QuBKyqtvwy8HOj+6mBKCPcVDqLteMGOOVrYMQdJWE+PbIwxJnhsLh1jjIkSFvjGGBMlIjrwRWSMiKwWkbUicruP9xNF5G3P+9+JSIfQVxlcfhzzTSKyQkSWisinItLejTqDqbZjrtTuPBFRzwUDEc2fYxaRX3v+rpeLyJuhrjHY/Pi3nSkin4nIIs+/77Fu1BksIvKyiOwUkWU1vC8i8pTnz2OpiPQLeKeqGpFfQCywDjgGSACWAD2qtbkGeN7zeiLwttt1h+CYhwPJnte/j4Zj9rRrCswD5gP93a47BH/PXYBFQLpnvaXbdYfgmKcAv/e87gFsdLvuAI95CNAPWFbD+2OBjwABBgHfBbrPSO7hDwDWqup6VS0C3sKZtbOy8YB3Vs53gBEiIiGsMdhqPWZV/UxV8z2r83Fucotk/vw9A9wHPAQUhrK4euLPMV8JPKvO3FSo6s4Q1xhs/hyzAime16nAIbPyRhJVnQcc7vL08cBr6pgPpIlI60D2GcmB78/8POVtVLUEZ66f5iGprn74PSeRx+U4PYRIVusxi0hfoJ2q/i+UhdUjf/6euwJdReRrEZkvIpE+GaE/x3w3cJGIZAEfAteGpjTX1PX/91pF8gNQfPXUq19j6k+bSOL38YjIRUB/YGi9VlT/DnvMIhIDPA5cGqqCQsCfv+c4nGGdYTi/xX0pIj1VNbeea6sv/hzzBcBUVX1URE4CXvccc1n9l+eKoOdXJPfws4B2ldbbcuiveOVtRCQO59fASL7D159jRkROByYD41T1YIhqqy+1HXNToCfwuYhsxBnrnBHhJ279/bf9vqoWq+oGYDXOD4BI5c8xXw78G0BVvwWScCYZa6j8+v+9LiI58L8HuohIRxFJwDkpO6NamxnAJM/r84C56jkbEqFqPWbP8MYLOGEf6eO6UMsxq+peVW2hqh1UtQPOeYtx6tzpHan8+bf9Hp6pyEWkBc4Qz/qQVhlc/hzzZmAEgIh0xwn8yH5O4eHNAC7xXK0zCNirqtsD+cCIHdJR1RIR+SPOJGyxwMuqulxE7gUWquoM4CWcX/vW4vTsJ7pXceD8POaHgSbAfzznpzer6jjXig6Qn8fcoPh5zN4JCVcApcCtqrrHvaoD4+cx3wz8U0RuxBnauDSSO3Ai8i+cIbkWnvMSdwHxAKr6PM55irHAWiAf+L+A9xnBf17GGGPqIJKHdIwxxtSBBb4xxkQJC3xjjIkSFvjGGBMlLPCNMSZKROxlmca4RUTuwLkmvAvOnDa7cCb8uk9V/+VmbcYcjvXwjfGT5waYGGAUMNuz+XFV7YMz0dULIhLvWoHG1MIC35jDEJEOIrJSRJ4DfsS51T1BVavc4amqa3Bujkl3oUxj/GKBb0ztjsWZprYvcALwafUGnodTrGkg01mYBsoC35jabfLMRw4whqpTTt8oIquB73Cm7zUmbFngG1O7A5VeDwAWVFp/XFWPBX4DvCYiSSGtzJg6sMA3xk8ichywSlVLq7+nqv8FFlIxO6sxYccC3xj/nQF8fJj37wVu8lzJY0zYsdkyjfGTiMwBLgl0TnJj3GKBb4wxUcJ+9TTGmChhgW+MMVHCAt8YY6KEBb4xxkQJC3xjjIkSFvjGGBMl/j/XB6nsPXcOJAAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# Plot the lhs and rhs of the hydrostatic equilibrium equation\n", "\n", "plt.figure()\n", "\n", "plt.plot(r/R, lhs, color='b', label='lhs')\n", "plt.plot(r/R, rhs, color='r', label='rhs', dashes=(2,1))\n", "\n", "plt.xlabel('r/R')\n", "\n", "plt.legend()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.7.3" } }, "nbformat": 4, "nbformat_minor": 2 }