generate_rest_eos_table_ext.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Generate a restricted SCvH EOS table for a H-He mixture\n",
    "Generate a new EOS table based on the SCvH equation of state [Saumon et al. 1995](https://ui.adsabs.harvard.edu/abs/1995ApJS...99..713S/abstract) for Hydrogen and Helium with extensions to lower temperatures [Vazan et al. 2013](https://academic.oup.com/mnras/article/434/4/3283/961935) restricted to low densities and temperatures where the isotherms behave sensible. For low densities and temperatures the EOS is in excellent agreement with an ideal gas."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Import modules"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from __future__ import print_function\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib as mpl\n",
    "from scipy import interpolate as interp\n",
    "#from scipy import optimize\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib widget\n",
    "%config InlineBackend.figure_format = 'retina'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mpl.rcParams['figure.dpi'] = 100"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def read_eos_table_dt(filename, delimiter=None):\n",
    "    data = np.loadtxt(filename, delimiter=delimiter) \n",
    "\n",
    "    logT_table   = data[:, 0]\n",
    "    logrho_table = data[:, 1]\n",
    "    logP_table   = data[:, 2]\n",
    "    logu_table   = data[:, 3]\n",
    "    logs_table   = data[:, 4]\n",
    "\n",
    "    # All the SCvH EOS tables are tabulated along isotherms\n",
    "    logT = np.unique(logT_table)\n",
    "    nT = np.size(logT)\n",
    "\n",
    "    # For the mixture tables the number of grid points in rho the same for each isotherm\n",
    "    nRho = int(np.size(logT_table)/nT)\n",
    "    \n",
    "    print(\"nT = {:} nRho = {:}\".format(nT, nRho))\n",
    "    \n",
    "    logrho = logrho_table[0:nRho]\n",
    "    \n",
    "    logrho_min = np.min(logrho)\n",
    "    logrho_max = np.max(logrho)\n",
    "    logT_min   = np.min(logT)\n",
    "    logT_max   = np.max(logT)\n",
    "\n",
    "    print(\"logrho_min = {:}\".format(logrho_min))\n",
    "    print(\"logrho_max = {:}\".format(logrho_max))\n",
    "    print(\"logT_min   = {:}\".format(logT_min))\n",
    "    print(\"logT_max   = {:}\".format(logT_max))\n",
    "    print()\n",
    "\n",
    "    # Split into arrays of constant T\n",
    "    logP_array = np.split(logP_table, nT)\n",
    "    logu_array = np.split(logu_table, nT)\n",
    "    logs_array = np.split(logs_table, nT)\n",
    "\n",
    "    # Generate 2d arrays\n",
    "    logP = np.vstack(logP_array)\n",
    "    logu = np.vstack(logu_array)\n",
    "    logs = np.vstack(logs_array)\n",
    "      \n",
    "    eos_table_dt = {\n",
    "        \"nT\":         nT,\n",
    "        \"nRho\":       nRho,\n",
    "        \"logT\":       logT,\n",
    "        \"logrho\":     logrho,\n",
    "        \"logrho_min\": logrho_min,\n",
    "        \"logrho_max\": logrho_max,\n",
    "        \"logT_min\":   logT_min,\n",
    "        \"logT_max\":   logT_max,\n",
    "        \"logP\":       logP,\n",
    "        \"logu\":       logu,\n",
    "        \"logs\":       logs,\n",
    "    }\n",
    "    \n",
    "    return eos_table_dt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Write an EOS table to an output file.\n",
    "def write_eos_table(output_file, logrho_axis, logT_axis, logP_array, logu_array, logs_array, delimiter=\" \", comments=\"#\", input_file=\"\"):\n",
    "    nRho = np.size(logrho_axis)\n",
    "    nT = np.size(logT_axis)\n",
    "    \n",
    "    # Store the different isotherms in 1d arrays\n",
    "    logrho = np.tile(logrho_axis, nT)\n",
    "    logT = np.repeat(logT_axis, nRho)\n",
    "\n",
    "    logP = logP_array.flatten(order='C')\n",
    "    logu = logu_array.flatten(order='C')\n",
    "    logs = logs_array.flatten(order='C')\n",
    "    \n",
    "    header = \" nT = {:} nRho= {:} (input file: {:})\\n\"\\\n",
    "             \" {:} {:} {:} {:} {:}\".format(nT, nRho, input_file, \"logT [K]\", \"logRho [g/cc]\", \"logP [barye]\", \"logE [erg/g]\", \"logS [erg/g/K]\")\n",
    "    \n",
    "    np.savetxt(output_file, np.column_stack([logT, logrho, logP, logu, logs]), header=header, fmt='%.8e', delimiter=delimiter, comments=comments)\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def press_ideal_gas(rho, T, mu=2.3, gamma=5.0/3.0):\n",
    "    kB = 1.38e-16  # erg/K\n",
    "    m_H = 1.67e-24 # g\n",
    "    return kB/(mu*m_H)*rho*T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def intenergy_ideal_gas(T, mu=2.3, gamma=5.0/3.0):\n",
    "    kB = 1.38e-16  # erg/K\n",
    "    m_H = 1.67e-24 # g\n",
    "    return kB/((gamma-1.0)*mu*m_H)*T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def entropy_ideal_gas(T, mu=2.3, gamma=5.0/3.0):\n",
    "    kB = 1.38e-16  # erg/K\n",
    "    m_H = 1.67e-24 # g\n",
    "    return kB/(mu*m_H)*(np.log(T/T0) - 1.0/(gamma-1.0)*np.log(rho/rho0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The mixture was obtained from the additive volume law (AVL) with Y=0.278 using bilinear spline interpolation and extrapolation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ext_table_file = \"scvh_extended_dt_hydrogen_722_helium_278.data\"\n",
    "scvh_hhe_ext = read_eos_table_dt(ext_table_file)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Test if write_eos_table works\n",
    "#write_eos_table(\"output.txt\", scvh_hhe_ext['logrho'], scvh_hhe_ext['logT'], scvh_hhe_ext['logP'], scvh_hhe_ext['logu'], scvh_hhe_ext['logs'], input_file=ext_table_file)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Skip some curves otherwise the plots are very hard to read\n",
    "nSkipT = 4\n",
    "nSkipRho = 10\n",
    "\n",
    "nSkipT = 1\n",
    "nSkipRho = 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot log(P), log(u) and log(s) of the original table\n",
    "fig, ax = plt.subplots(3, 2)\n",
    "\n",
    "x, y = fig.get_size_inches()\n",
    "\n",
    "fig.set_size_inches(2*x, 3*y)\n",
    "\n",
    "# P(rho, T=const)\n",
    "for i in range(0, scvh_hhe_ext['nT'], nSkipT):\n",
    "    ax[0][0].plot(scvh_hhe_ext['logrho'], scvh_hhe_ext['logP'][i,:], '-')\n",
    "    \n",
    "ax[0][0].set(xlabel=\"log(rho) [g cm^${-3}$]\", ylabel=\"log(P) [erg cm$^{-3}$]\")\n",
    "ax[0][0].set(title=\"SCvH EOS H (T=const.)\")\n",
    "\n",
    "# P(rho=const, T)\n",
    "for i in range(0, scvh_hhe_ext['nRho'], nSkipRho):\n",
    "    ax[0][1].plot(scvh_hhe_ext['logT'], scvh_hhe_ext['logP'][:,i], '-')\n",
    "\n",
    "ax[0][1].set(xlabel=\"log(T) [K]\", ylabel=\"log(P) [erg cm$^{-3}$]\")\n",
    "ax[0][1].set(title=\"SCvH EOS H (rho=const.)\")\n",
    "\n",
    "# u(rho, T=const)\n",
    "for i in range(0, scvh_hhe_ext['nT'], nSkipT):\n",
    "    ax[1][0].plot(scvh_hhe_ext['logrho'], scvh_hhe_ext['logu'][i,:], '-')\n",
    "\n",
    "ax[1][0].set(xlabel=\"log(rho) [g cm^${-3}$]\", ylabel=\"log(u) [erg g$^{-1}$]\")\n",
    "\n",
    "# u(rho=const, T)\n",
    "for i in range(0, scvh_hhe_ext['nRho'], nSkipRho):\n",
    "    ax[1][1].plot(scvh_hhe_ext['logT'], scvh_hhe_ext['logu'][:,i], '-')\n",
    "\n",
    "ax[1][1].set(xlabel=\"log(T) [K]\", ylabel=\"log(u) [erg g$^{-1}$]\")\n",
    "\n",
    "# s(rho, T=const)\n",
    "for i in range(0, scvh_hhe_ext['nT'], nSkipT):\n",
    "    ax[2][0].plot(scvh_hhe_ext['logrho'], scvh_hhe_ext['logs'][i,:], '-')\n",
    "\n",
    "ax[2][0].set(xlabel=\"log(rho) [g cm^${-3}$]\", ylabel=\"log(s) [erg g$^{-1}$ K$^{-1}$]\")\n",
    "\n",
    "# s(rho=const, T)\n",
    "for i in range(0, scvh_hhe_ext['nRho'], nSkipRho):\n",
    "    ax[2][1].plot(scvh_hhe_ext['logT'], scvh_hhe_ext['logs'][:,i], '-')\n",
    "\n",
    "ax[2][1].set(xlabel=\"log(T) [K]\", ylabel=\"log(s) [erg g$^{-1}$ K$^{-1}$]\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def generate_restricted_eos_table(scvh_eos_table, logrho_min, logrho_max, logT_min, logT_max):\n",
    "    # Find all grid point that are in the range of the restricted EOS table\n",
    "    index_logrho_grid = np.where((scvh_eos_table['logrho'] >= logrho_min) & (scvh_eos_table['logrho'] <= logrho_max))[0]\n",
    "    index_logT_grid = np.where((scvh_eos_table['logT'] >= logT_min) & (scvh_eos_table['logT'] < logT_max))[0]\n",
    "\n",
    "    nRho_grid = np.size(index_logrho_grid)\n",
    "    nT_grid = np.size(index_logT_grid)\n",
    "    print(\"nRho= {:} nT= {:}\".format(nRho_grid, nT_grid))\n",
    "\n",
    "    logrho_grid = scvh_eos_table['logrho'][index_logrho_grid]\n",
    "    logT_grid = scvh_eos_table['logT'][index_logT_grid]\n",
    "    \n",
    "    index_logrho_grid_min = np.min(index_logrho_grid)\n",
    "    index_logrho_grid_max = np.max(index_logrho_grid)  \n",
    "    index_logT_grid_min = np.min(index_logT_grid)\n",
    "    index_logT_grid_max = np.max(index_logT_grid)\n",
    "    \n",
    "    # Note that when slicing the last index is not included\n",
    "    logP_grid = scvh_eos_table['logP'][index_logT_grid_min:index_logT_grid_max+1,index_logrho_grid_min:index_logrho_grid_max+1]\n",
    "    logu_grid = scvh_eos_table['logu'][index_logT_grid_min:index_logT_grid_max+1,index_logrho_grid_min:index_logrho_grid_max+1]\n",
    "    logs_grid = scvh_eos_table['logs'][index_logT_grid_min:index_logT_grid_max+1,index_logrho_grid_min:index_logrho_grid_max+1]\n",
    "\n",
    "    scvh_eos_restricted = {\n",
    "        \"nT\":         nT_grid,\n",
    "        \"nRho\":       nRho_grid,\n",
    "        \"logT\":       logT_grid,\n",
    "        \"logrho\":     logrho_grid,\n",
    "        \"logrho_min\": logrho_grid_min,\n",
    "        \"logrho_max\": logrho_grid_max,\n",
    "        \"logT_min\":   logT_grid_min,\n",
    "        \"logT_max\":   logT_grid_max,\n",
    "        \"logP\":       logP_grid,\n",
    "        \"logu\":       logu_grid,\n",
    "        \"logs\":       logs_grid,\n",
    "    }\n",
    "    \n",
    "    return scvh_eos_restricted"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Limit the data to the region in which we are interested\n",
    "logrho_grid_min = -14.0\n",
    "logrho_grid_max = -4.0\n",
    "\n",
    "logT_grid_min = np.min(scvh_hhe_ext['logT'])\n",
    "logT_grid_max = 3.46\n",
    "logT_grid_max = 3.49\n",
    "\n",
    "\"\"\"\n",
    "print(\"logrho_grid_min= {:} rho_grid_min= {:}\".format(logrho_grid_min, 10.0**logrho_grid_min))\n",
    "print(\"logrho_grid_max= {:} rho_grid_max= {:}\".format(logrho_grid_max, 10.0**logrho_grid_max))\n",
    "print(\"logT_grid_min= {:} T_grid_min= {:}\".format(logT_grid_min, 10.0**logT_grid_min))\n",
    "print(\"logT_grid_max= {:} T_grid_max= {:}\".format(logT_grid_max, 10.0**logT_grid_max))\n",
    "print()\n",
    "\"\"\"\n",
    "scvh_hhe_rest = generate_restricted_eos_table(scvh_hhe_ext, logrho_grid_min, logrho_grid_max, logT_grid_min, logT_grid_max)\n",
    "\n",
    "print(\"logrho_grid_min= {:} rho_grid_min= {:}\".format(scvh_hhe_rest['logrho_min'], 10.0**scvh_hhe_rest['logrho_min']))\n",
    "print(\"logrho_grid_max= {:} rho_grid_max= {:}\".format(scvh_hhe_rest['logrho_max'], 10.0**scvh_hhe_rest['logrho_max']))\n",
    "print(\"logT_grid_min= {:} T_grid_min= {:}\".format(scvh_hhe_rest['logT_min'], logT_grid_min, 10.0**scvh_hhe_rest['logT_min']))\n",
    "print(\"logT_grid_max= {:} T_grid_max= {:}\".format(scvh_hhe_rest['logT_max'], 10.0**scvh_hhe_rest['logT_max']))\n",
    "print()\n",
    "\n",
    "# Write a new EOS table\n",
    "write_eos_table(\"scvh_extended_dt_hydrogen_722_helium_278_lowrhot.txt\", scvh_hhe_rest['logrho'], scvh_hhe_rest['logT'], scvh_hhe_rest['logP'], scvh_hhe_rest['logu'], scvh_hhe_rest['logs'], input_file=ext_table_file)\n",
    "\n",
    "# Now load the new EOS table to check if it was properly written\n",
    "scvh_hhe_rest = read_eos_table_dt(\"scvh_extended_dt_hydrogen_722_helium_278_lowrhot.txt\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Skip some curves otherwise the plots are very hard to read\n",
    "nSkipT = 3\n",
    "nSkipRho = 10"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot the restricted grid and check that logP, logu and logs agree with logrho and logT\n",
    "fig, ax = plt.subplots(3, 2)\n",
    "\n",
    "x, y = fig.get_size_inches()\n",
    "\n",
    "fig.set_size_inches(2*x, 3*y)\n",
    "\n",
    "fig.suptitle(r\"{:} $\\leq \\log(\\rho) \\leq$ {:}, {:} $\\leq \\log(T) \\leq$ {:}\".format(scvh_hhe_rest['logrho_min'], scvh_hhe_rest['logrho_max'], scvh_hhe_rest['logT_min'], scvh_hhe_rest['logT_max']), x=0.5, y=0.0, verticalalignment='bottom', horizontalalignment='center')\n",
    "\n",
    "# Generate interpolation functions\n",
    "logP_ext_int = interp.interp2d(scvh_hhe_ext['logrho'], scvh_hhe_ext['logT'], scvh_hhe_ext['logP'], kind='linear')\n",
    "logu_ext_int = interp.interp2d(scvh_hhe_ext['logrho'], scvh_hhe_ext['logT'], scvh_hhe_ext['logu'], kind='linear')\n",
    "logs_ext_int = interp.interp2d(scvh_hhe_ext['logrho'], scvh_hhe_ext['logT'], scvh_hhe_ext['logs'], kind='linear')\n",
    "\n",
    "# P(rho, T=const)\n",
    "for i in range(0, scvh_hhe_rest['nT'], nSkipT):\n",
    "    ax[0][0].plot(scvh_hhe_rest['logrho'], scvh_hhe_rest['logP'][i,:], '-')\n",
    "    ax[0][0].plot(scvh_hhe_rest['logrho'], logP_ext_int(scvh_hhe_rest['logrho'], scvh_hhe_rest['logT'][i]), '--')\n",
    "\n",
    "    \n",
    "ax[0][0].set(xlabel=\"log(rho) [g cm^${-3}$]\", ylabel=\"log(P) [erg cm$^{-3}$]\")\n",
    "ax[0][0].set(title=\"SCvH EOS H (T=const.)\")\n",
    "\n",
    "# P(rho=const, T)\n",
    "for i in range(0, scvh_hhe_rest['nRho'], nSkipRho):\n",
    "    ax[0][1].plot(scvh_hhe_rest['logT'], scvh_hhe_rest['logP'][:,i], '-')\n",
    "    ax[0][1].plot(scvh_hhe_rest['logT'], logP_ext_int(scvh_hhe_rest['logrho'][i], scvh_hhe_rest['logT']), '--')\n",
    "\n",
    "\n",
    "ax[0][1].set(xlabel=\"log(T) [K]\", ylabel=\"log(P) [erg cm$^{-3}$]\")\n",
    "ax[0][1].set(title=\"SCvH EOS H (rho=const.)\")\n",
    "\n",
    "# u(rho, T=const)\n",
    "for i in range(0, scvh_hhe_rest['nT'], nSkipT):\n",
    "    ax[1][0].plot(scvh_hhe_rest['logrho'], scvh_hhe_rest['logu'][i,:], '-')\n",
    "    ax[1][0].plot(scvh_hhe_rest['logrho'], logu_ext_int(scvh_hhe_rest['logrho'], scvh_hhe_rest['logT'][i]), '--')\n",
    "\n",
    "\n",
    "ax[1][0].set(xlabel=\"log(rho) [g cm^${-3}$]\", ylabel=\"log(u) [erg g$^{-1}$]\")\n",
    "\n",
    "# u(rho=const, T)\n",
    "for i in range(0, scvh_hhe_rest['nRho'], nSkipRho):\n",
    "    ax[1][1].plot(scvh_hhe_rest['logT'], scvh_hhe_rest['logu'][:,i], '-')\n",
    "    ax[1][1].plot(scvh_hhe_rest['logT'], logu_ext_int(scvh_hhe_rest['logrho'][i], scvh_hhe_rest['logT']), '--')\n",
    "\n",
    "\n",
    "ax[1][1].set(xlabel=\"log(T) [K]\", ylabel=\"log(u) [erg g$^{-1}$]\")\n",
    "\n",
    "# s(rho, T=const)\n",
    "for i in range(0, scvh_hhe_rest['nT'], nSkipT):\n",
    "    ax[2][0].plot(scvh_hhe_rest['logrho'], scvh_hhe_rest['logs'][i,:], '-')\n",
    "    ax[2][0].plot(scvh_hhe_rest['logrho'], logs_ext_int(scvh_hhe_rest['logrho'], scvh_hhe_rest['logT'][i]), '--')\n",
    "\n",
    "\n",
    "ax[2][0].set(xlabel=\"log(rho) [g cm^${-3}$]\", ylabel=\"log(s) [erg g$^{-1}$ K$^{-1}$]\")\n",
    "\n",
    "# s(rho=const, T)\n",
    "for i in range(0, scvh_hhe_rest['nRho'], nSkipRho):\n",
    "    ax[2][1].plot(scvh_hhe_rest['logT'], scvh_hhe_rest['logs'][:,i], '-')\n",
    "    ax[2][1].plot(scvh_hhe_rest['logT'], logs_ext_int(scvh_hhe_rest['logrho'][i], scvh_hhe_rest['logT']), '--')\n",
    "\n",
    "ax[2][1].set(xlabel=\"log(T) [K]\", ylabel=\"log(s) [erg g$^{-1}$ K$^{-1}$]\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Skip some curves otherwise the plots are very hard to read\n",
    "nSkipT = 3\n",
    "nSkipRho = 10"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compare the reduced grid to an ideal gas\n",
    "fig, ax = plt.subplots(2, 2)\n",
    "\n",
    "x, y = fig.get_size_inches()\n",
    "\n",
    "fig.set_size_inches(2*x, 2*y)\n",
    "\n",
    "# P(rho, T=const)\n",
    "ax[0][0].set_prop_cycle(None)\n",
    "for i in range(0, scvh_hhe_rest['nT'], nSkipT):\n",
    "    ax[0][0].plot(scvh_hhe_rest['logrho'], scvh_hhe_rest['logP'][i,:], '-')\n",
    "\n",
    "ax[0][0].set_prop_cycle(None)\n",
    "for i in range(0, scvh_hhe_rest['nT'], nSkipT):\n",
    "    ax[0][0].plot(scvh_hhe_rest['logrho'], np.log10(press_ideal_gas(10.0**scvh_hhe_rest['logrho'], 10.0**scvh_hhe_rest['logT'][i])), '--')\n",
    "\n",
    "ax[0][0].set(xlabel=\"log(rho) [g cm^${-3}$]\", ylabel=\"log(P) [erg cm$^{-3}$]\")\n",
    "ax[0][0].set(title=\"SCvH EOS H (T=const.)\")\n",
    "\n",
    "# P(rho=const, T)\n",
    "ax[0][1].set_prop_cycle(None)\n",
    "for i in range(0, scvh_hhe_rest['nRho'], nSkipRho):\n",
    "    ax[0][1].plot(scvh_hhe_rest['logT'], scvh_hhe_rest['logP'][:,i], '-')\n",
    "\n",
    "ax[0][1].set_prop_cycle(None)\n",
    "for i in range(0, scvh_hhe_rest['nRho'], nSkipRho):\n",
    "    ax[0][1].plot(scvh_hhe_rest['logT'], np.log10(press_ideal_gas(10.0**scvh_hhe_rest['logrho'][i], 10.0**scvh_hhe_rest['logT'])), '--')\n",
    "\n",
    "    \n",
    "ax[0][1].set(xlabel=\"log(T) [K]\", ylabel=\"log(P) [erg cm$^{-3}$]\")\n",
    "ax[0][1].set(title=\"SCvH EOS H (rho=const.)\")\n",
    "\n",
    "# u(rho, T=const)\n",
    "ax[1][0].set_prop_cycle(None)\n",
    "for i in range(0, scvh_hhe_rest['nT'], nSkipT):\n",
    "    ax[1][0].plot(scvh_hhe_rest['logrho'], scvh_hhe_rest['logu'][i,:], '-')\n",
    "\n",
    "ax[1][0].set_prop_cycle(None)\n",
    "for i in range(0, scvh_hhe_rest['nT'], nSkipT):\n",
    "    ax[1][0].plot(scvh_hhe_rest['logrho'], np.ones(scvh_hhe_rest['nRho'])*np.log10(intenergy_ideal_gas(10.0**scvh_hhe_rest['logT'][i])), '--')\n",
    "\n",
    "        #ax[1][0].plot(logrho_grid, *np.log10(intenergy_ideal_gas(10**logT_grid[i])), '--')\n",
    "\n",
    "    \n",
    "ax[1][0].set(xlabel=\"log(rho) [g cm^${-3}$]\", ylabel=\"log(u) [erg g$^{-1}$]\")\n",
    "ax[1][0].set(title=\"SCvH EOS H (T=const.)\")\n",
    "\n",
    "\n",
    "# u(rho=const, T)\n",
    "ax[1][1].set_prop_cycle(None)\n",
    "for i in range(0, scvh_hhe_rest['nRho'], nSkipRho):\n",
    "    ax[1][1].plot(scvh_hhe_rest['logT'], scvh_hhe_rest['logu'][:,i], '-')\n",
    "\n",
    "ax[1][1].set_prop_cycle(None)\n",
    "for i in range(0, scvh_hhe_rest['nRho'], nSkipRho):\n",
    "    ax[1][1].plot(scvh_hhe_rest['logT'], np.log10(intenergy_ideal_gas(10**scvh_hhe_rest['logT'])), '--')\n",
    "\n",
    "ax[1][1].get_xlim()\n",
    "ax[1][1].get_ylim()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"logrho: nRho= {:}\".format(scvh_hhe_rest['nRho']))\n",
    "print(scvh_hhe_rest['logrho'])\n",
    "print()\n",
    "print(\"logT: nT= {:}\".format(scvh_hhe_rest['nT']))\n",
    "print(scvh_hhe_rest['logT'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "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.10.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}