plot_test_derivs.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Plot the derivatives"
   ]
  },
  {
   "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",
    "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": [
    "data = np.loadtxt(\"scvh_hhe_y0.275_dt_cgs_lowrhot.csv\", delimiter=\",\", skiprows=0) \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 SCVH EOS tables obtained from Ravit have the same size\n",
    "#nRho = 201\n",
    "#nT   = 100\n",
    "nRho = 80\n",
    "nT   = 48\n",
    "\n",
    "logrho_table = logrho_table[0:nRho]\n",
    "logT_table = logT_table[0:np.size(logT_table):nRho]\n",
    "\n",
    "logrho_min = np.min(logrho_table)\n",
    "logrho_max = np.max(logrho_table)\n",
    "logT_min   = np.min(logT_table)\n",
    "logT_max   = np.max(logT_table)\n",
    "    \n",
    "dlogrho = logrho_table[1:]-logrho_table[:-1]\n",
    "dlogT = logT_table[1:]-logT_table[:-1]\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",
    "#print(\"logrho_table = {:}\".format(logrho_table))\n",
    "#print(\"logT_table   = {:}\".format(logT_table))\n",
    "#print()\n",
    "\n",
    "# Mark where the original SCvH EOS table was extended below T=100K\n",
    "index = np.min(np.where(logT_table>=2.0))\n",
    "print(\"The EOS table was extended until index={:} logT={:}\".format(index, logT_table[index]))\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)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Normal variables\n",
    "T_table = 10.0**logT_table\n",
    "rho_table = 10.0**logrho_table\n",
    "\n",
    "P = 10.0**logP\n",
    "u = 10.0**logu\n",
    "s = 10.0**logs\n",
    "\n",
    "rho_min = 10.0**logrho_min\n",
    "rho_max = 10.0**logrho_max\n",
    "T_min = 10.0**logT_min\n",
    "T_max = 10.0**logT_max\n",
    "\n",
    "print(\"rho_min = {:}\".format(rho_min))\n",
    "print(\"rho_max = {:}\".format(rho_max))\n",
    "print(\"T_min   = {:}\".format(T_min))\n",
    "print(\"T_max   = {:}\".format(T_max))\n",
    "print()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we generate interpolation function from ```scipy.interp``` to obtain an independent calculation of the derivatives. Since the ```scipy.interp.interp2d``` does not provide derivatives we use the $1D$ interpolation functions. These also support *Akima splines* which have many desirable properties."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Generate interpolation functions to calculate the derivatives\n",
    "logP_table_int_rho = list()\n",
    "logP_table_int_T   = list()\n",
    "logu_table_int_rho = list()\n",
    "logu_table_int_T   = list()\n",
    "logs_table_int_rho = list()\n",
    "logs_table_int_T   = list()\n",
    "\n",
    "dlogP_dlogrho_int_rho = list()\n",
    "dlogP_dlogT_int_T = list()\n",
    "\n",
    "dlogu_dlogrho_int_rho = list()\n",
    "dlogu_dlogT_int_T = list()\n",
    "\n",
    "dlogs_dlogrho_int_rho = list()\n",
    "dlogs_dlogT_int_T = list()\n",
    "\n",
    "# Interpolation for T=const.\n",
    "# Akima splines\n",
    "\"\"\"\n",
    "for i in range(0, nT):\n",
    "    logP_table_int_rho.append(interp.Akima1DInterpolator(logrho_table, logP[i,:]))\n",
    "    logu_table_int_rho.append(interp.Akima1DInterpolator(logrho_table, logu[i,:]))\n",
    "    logs_table_int_rho.append(interp.Akima1DInterpolator(logrho_table, logs[i,:]))\n",
    "\"\"\"\n",
    "# Linear interpolation\n",
    "for i in range(0, nT):\n",
    "    logP_table_int_rho.append(interp.InterpolatedUnivariateSpline(logrho_table, logP[i,:], k=1))\n",
    "    logu_table_int_rho.append(interp.InterpolatedUnivariateSpline(logrho_table, logu[i,:], k=1))\n",
    "    logs_table_int_rho.append(interp.InterpolatedUnivariateSpline(logrho_table, logs[i,:], k=1))\n",
    "\n",
    "# Calculate first derivatives\n",
    "for i in range(0, nT):\n",
    "    dlogP_dlogrho_int_rho.append(logP_table_int_rho[i].derivative())\n",
    "    dlogu_dlogrho_int_rho.append(logu_table_int_rho[i].derivative())\n",
    "    dlogs_dlogrho_int_rho.append(logs_table_int_rho[i].derivative())\n",
    "    \n",
    "# Interpolation for rho=const.\n",
    "\"\"\"\n",
    "# Akima splines\n",
    "for i in range(0, nRho):\n",
    "    logP_table_int_T.append(interp.Akima1DInterpolator(logT_table, logP[:,i]))\n",
    "    logu_table_int_T.append(interp.Akima1DInterpolator(logT_table, logu[:,i]))\n",
    "    logs_table_int_T.append(interp.Akima1DInterpolator(logT_table, logs[:,i]))\n",
    "\"\"\"\n",
    "# Linear interpolation\n",
    "for i in range(0, nRho):\n",
    "    logP_table_int_T.append(interp.InterpolatedUnivariateSpline(logT_table, logP[:,i], k=1))\n",
    "    logu_table_int_T.append(interp.InterpolatedUnivariateSpline(logT_table, logu[:,i], k=1))\n",
    "    logs_table_int_T.append(interp.InterpolatedUnivariateSpline(logT_table, logs[:,i], k=1))\n",
    "\n",
    "# Calculate first derivatives\n",
    "for i in range(0, nRho):\n",
    "    dlogP_dlogT_int_T.append(logP_table_int_T[i].derivative())\n",
    "    dlogu_dlogT_int_T.append(logu_table_int_T[i].derivative())\n",
    "    dlogs_dlogT_int_T.append(logs_table_int_T[i].derivative())\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Calculate the normal derivatives, e.g.,\n",
    "$$ \\left. \\frac{\\partial P}{\\partial \\rho} \\right|_{T} = \\frac{P}{\\rho} \\left. \\frac{\\partial \\log P}{\\partial \\log \\rho} \\right|_{\\log T}   $$\n",
    "$$ \\left. \\frac{\\partial P}{\\partial T} \\right|_{\\rho} = \\frac{P}{T} \\left. \\frac{\\partial \\log P}{\\partial \\log T} \\right|_{\\log \\rho}   $$\n",
    "from ```scipy.interp```.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#### This is probably not useful\n",
    "dP_drho_int_rho = list()\n",
    "dP_dT_int_T = list()\n",
    "\n",
    "du_drho_int_rho = list()\n",
    "du_dT_int_T = list()\n",
    "\n",
    "ds_drho_int_rho = list()\n",
    "ds_dT_int_T = list()\n",
    "\n",
    "# Derivatives for T=const.\n",
    "for i in range(0, nT):\n",
    "    dlogP_dlogrho_int_rho.append(logP_table_int_rho[i].derivative())\n",
    "    dlogu_dlogrho_int_rho.append(logu_table_int_rho[i].derivative())\n",
    "    dlogs_dlogrho_int_rho.append(logs_table_int_rho[i].derivative())\n",
    "    \n",
    "# Interpolation for rho=const.\n",
    "\"\"\"\n",
    "# Akima splines\n",
    "for i in range(0, nRho):\n",
    "    logP_table_int_T.append(interp.Akima1DInterpolator(logT_table, logP[:,i]))\n",
    "    logu_table_int_T.append(interp.Akima1DInterpolator(logT_table, logu[:,i]))\n",
    "    logs_table_int_T.append(interp.Akima1DInterpolator(logT_table, logs[:,i]))\n",
    "\"\"\"\n",
    "# Linear interpolation\n",
    "for i in range(0, nRho):\n",
    "    logP_table_int_T.append(interp.InterpolatedUnivariateSpline(logT_table, logP[:,i], k=1))\n",
    "    logu_table_int_T.append(interp.InterpolatedUnivariateSpline(logT_table, logu[:,i], k=1))\n",
    "    logs_table_int_T.append(interp.InterpolatedUnivariateSpline(logT_table, logs[:,i], k=1))\n",
    "\n",
    "# Calculate first derivatives\n",
    "for i in range(0, nRho):\n",
    "    dlogP_dlogT_int_T.append(logP_table_int_T[i].derivative())\n",
    "    dlogu_dlogT_int_T.append(logu_table_int_T[i].derivative())\n",
    "    dlogs_dlogT_int_T.append(logs_table_int_T[i].derivative())\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load derivatives\n",
    "logrho_int = np.loadtxt(\"testscvheosinterp_rhoaxis.txt\")\n",
    "logT_int = np.loadtxt(\"testscvheosinterp_taxis.txt\")\n",
    "\n",
    "dlogPdlogrho_int = np.loadtxt(\"testscvheosinterp_dlogpdlogrho.txt\")\n",
    "dlogPdlogT_int = np.loadtxt(\"testscvheosinterp_dlogpdlogt.txt\")\n",
    "\n",
    "dlogudlogrho_int = np.loadtxt(\"testscvheosinterp_dlogudlogrho.txt\")\n",
    "dlogudlogT_int = np.loadtxt(\"testscvheosinterp_dlogudlogt.txt\")\n",
    "\n",
    "dlogsdlogrho_int = np.loadtxt(\"testscvheosinterp_dlogsdlogrho.txt\")\n",
    "dlogsdlogT_int = np.loadtxt(\"testscvheosinterp_dlogsdlogt.txt\")\n",
    "\n",
    "nRho_int = np.size(logrho_int)\n",
    "nT_int = np.size(logT_int)\n",
    "\n",
    "print(\"nRho_int= {:} nT_int= {:}\".format(nRho_int, nT_int))\n"
   ]
  },
  {
   "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 = 8\n",
    "\n",
    "nSkipT = 1\n",
    "nSkipRho = 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First we verify that the interpolation functions from ```scipy.interp``` produce sensible results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot log(P), log(u) and log(s) to verify that the python interpolation functions work\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, nT, nSkipT):\n",
    "    ax[0][0].plot(logrho_table, logP[i,:], '-')\n",
    "    ax[0][0].plot(logrho_int, logP_table_int_rho[i](logrho_int), '--')\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-He, T=const.)\")\n",
    "\n",
    "# P(rho=const, T)\n",
    "for i in range(0, nRho, nSkipRho):\n",
    "    ax[0][1].plot(logT_table, logP[:,i], '-')\n",
    "    ax[0][1].plot(logT_int, logP_table_int_T[i](logT_int), '--')\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-He (rho=const.)\")\n",
    "\n",
    "# u(rho, T=const)\n",
    "for i in range(0, nT, nSkipT):\n",
    "    ax[1][0].plot(logrho_table, logu[i,:], '-')\n",
    "    ax[1][0].plot(logrho_int, logu_table_int_rho[i](logrho_int), '--')\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, nRho, nSkipRho):\n",
    "    ax[1][1].plot(logT_table, logu[:,i], '-')\n",
    "    ax[1][1].plot(logT_int, logu_table_int_T[i](logT_int), '--')\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, nT, nSkipT):\n",
    "    ax[2][0].plot(logrho_table, logs[i,:], '-')\n",
    "    ax[2][0].plot(logrho_int, logs_table_int_rho[i](logrho_int), '--')\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, nRho, nSkipRho):\n",
    "    ax[2][1].plot(logT_table, logs[:,i], '-')\n",
    "    ax[2][1].plot(logT_int, logs_table_int_T[i](logT_int), '--')\n",
    "\n",
    "ax[2][1].set(xlabel=\"log(T) [K]\", ylabel=\"log(s) [erg g$^{-1}$ K$^{-1}$]\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we plot the derivatives obtained with ```scipy.interp``` and compare to the results from ```testscvheosderivs.c```."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot the derivatives\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",
    "# dlogP/dlogrho(rho, T=const)\n",
    "ax[0][0].set_prop_cycle(None)\n",
    "for i in range(0, nT, nSkipT):\n",
    "    ax[0][0].plot(logrho_table, dlogP_dlogrho_int_rho[i](logrho_table), '-')\n",
    "\n",
    "for i in range(0, nT_int, nSkipT):\n",
    "    ax[0][0].plot(logrho_int, dlogPdlogrho_int[i,:], '--')\n",
    "    \n",
    "ax[0][0].set(xlabel=\"log(rho) [g cm^${-3}$]\", ylabel=\"dlog(P)/dlog(rho)\")\n",
    "ax[0][0].set(title=\"SCvH EOS H-He, T=const.)\")\n",
    "\n",
    "# dlogP/dlogT(rho=const, T)\n",
    "ax[0][1].set_prop_cycle(None)\n",
    "\n",
    "for i in range(0, nRho, nSkipRho):\n",
    "    ax[0][1].plot(logT_table, dlogP_dlogT_int_T[i](logT_table), '-')\n",
    "\n",
    "for i in range(0, nRho_int, nSkipRho):\n",
    "    ax[0][1].plot(logT_int, dlogPdlogT_int[:,i], '--')\n",
    "\n",
    "ax[0][1].set(xlabel=\"log(T) [K]\", ylabel=\"dlog(P)/dlog(T)\")\n",
    "ax[0][1].set(title=\"SCvH EOS H-He (rho=const.)\")\n",
    "\n",
    "# dlogu/dlogrho(rho, T=const)\n",
    "ax[1][0].set_prop_cycle(None)\n",
    "for i in range(0, nT, nSkipT):\n",
    "    ax[1][0].plot(logrho_table, dlogu_dlogrho_int_rho[i](logrho_table), '-')\n",
    "\n",
    "for i in range(0, nT_int, nSkipT):\n",
    "    ax[1][0].plot(logrho_int, dlogudlogrho_int[i,:], '--')\n",
    "    \n",
    "ax[1][0].set(xlabel=\"log(rho) [g cm^${-3}$]\", ylabel=\"dlog(u)/dlog(rho)\")\n",
    "ax[1][0].set(title=\"SCvH EOS H-He, T=const.)\")\n",
    "\n",
    "# dlogu/dlogT(rho=const, T)\n",
    "ax[1][1].set_prop_cycle(None)\n",
    "\n",
    "for i in range(0, nRho, nSkipRho):\n",
    "    ax[1][1].plot(logT_table, dlogu_dlogT_int_T[i](logT_table), '-')\n",
    "\n",
    "for i in range(0, nRho_int, nSkipRho):\n",
    "    ax[1][1].plot(logT_int, dlogudlogT_int[:,i], '--')\n",
    "\n",
    "ax[1][1].set(xlabel=\"log(T) [K]\", ylabel=\"dlog(u)/dlog(T)\")\n",
    "ax[1][1].set(title=\"SCvH EOS H-He (rho=const.)\")\n",
    "\n",
    "# dlogs/dlogrho(rho, T=const)\n",
    "ax[2][0].set_prop_cycle(None)\n",
    "for i in range(0, nT, nSkipT):\n",
    "    ax[2][0].plot(logrho_table, dlogs_dlogrho_int_rho[i](logrho_table), '-')\n",
    "\n",
    "for i in range(0, nT_int, nSkipT):\n",
    "    ax[2][0].plot(logrho_int, dlogsdlogrho_int[i,:], '--')\n",
    "    \n",
    "ax[2][0].set(xlabel=\"log(rho) [g cm^${-3}$]\", ylabel=\"dlog(s)/dlog(rho)\")\n",
    "\n",
    "# dlogs/dlogT(rho=const, T)\n",
    "ax[2][1].set_prop_cycle(None)\n",
    "\n",
    "for i in range(0, nRho, nSkipRho):\n",
    "    ax[2][1].plot(logT_table, dlogs_dlogT_int_T[i](logT_table), '-')\n",
    "\n",
    "for i in range(0, nRho_int, nSkipRho):\n",
    "    ax[2][1].plot(logT_int, dlogsdlogT_int[:,i], '--')\n",
    "\n",
    "ax[2][1].set(xlabel=\"log(T) [K]\", ylabel=\"dlog(s)/dlog(T)\")\n",
    "ax[2][1].set(title=\"SCvH EOS H-He (rho=const.)\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally we plot the derivatives in normal variables with ```scipy.interp``` and compare to the results from ```testscvheosderivs.c```."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot the derivatives\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",
    "# dP/drho(rho, T=const)\n",
    "ax[0][0].set_prop_cycle(None)\n",
    "for i in range(0, nT, nSkipT):\n",
    "    ax[0][0].semilogx(rho_table, P[i,:]/rho_table[i]*dlogP_dlogrho_int_rho[i](np.log10(rho_table)), '-')\n",
    "\n",
    "#for i in range(0, nT_int, nSkipT):\n",
    "#    ax[0][0].plot(logrho_int, dlogPdlogrho_int[i,:], '--')\n",
    "    \n",
    "ax[0][0].set(xlabel=\"rho [g cm^${-3}$]\", ylabel=\"dP/drho\")\n",
    "ax[0][0].set(title=\"SCvH EOS H-He, T=const.)\")\n",
    "\n",
    "# dP/dT(rho=const, T)\n",
    "ax[0][1].set_prop_cycle(None)\n",
    "\n",
    "for i in range(0, nRho, nSkipRho):\n",
    "    ax[0][1].semilogx(T_table, P[:,i]/T_table*dlogP_dlogT_int_T[i](np.log10(T_table)), '-')\n",
    "\n",
    "#for i in range(0, nRho_int, nSkipRho):\n",
    "#    ax[0][1].plot(logT_int, dlogPdlogT_int[:,i], '--')\n",
    "\n",
    "ax[0][1].set(xlabel=\"T [K]\", ylabel=\"dP/dT\")\n",
    "ax[0][1].set(title=\"SCvH EOS H-He (rho=const.)\")\n",
    "\n",
    "# du/drho(rho, T=const)\n",
    "ax[1][0].set_prop_cycle(None)\n",
    "for i in range(0, nT, nSkipT):\n",
    "    #ax[1][0].plot(logrho_table, dlogu_dlogrho_int_rho[i](logrho_table), '-')\n",
    "    ax[1][0].semilogx(rho_table, u[i,:]/rho_table[i]*dlogu_dlogrho_int_rho[i](np.log10(rho_table)), '-')\n",
    "\n",
    "\n",
    "#for i in range(0, nT_int, nSkipT):\n",
    "#    ax[1][0].plot(logrho_int, dlogudlogrho_int[i,:], '--')\n",
    "    \n",
    "ax[1][0].set(xlabel=\"rho [g cm^${-3}$]\", ylabel=\"du/drho\")\n",
    "ax[1][0].set(title=\"SCvH EOS H-He, T=const.)\")\n",
    "\n",
    "# du/dT(rho=const, T)\n",
    "ax[1][1].set_prop_cycle(None)\n",
    "\n",
    "for i in range(0, nRho, nSkipRho):\n",
    "    ax[1][1].semilogx(T_table, u[:,i]/T_table*dlogu_dlogT_int_T[i](np.log10(T_table)), '-')\n",
    "\n",
    "#for i in range(0, nRho_int, nSkipRho):\n",
    "#    ax[1][1].plot(logT_int, dlogudlogT_int[:,i], '--')\n",
    "\n",
    "ax[1][1].set(xlabel=\"T [K]\", ylabel=\"du/dT\")\n",
    "ax[1][1].set(title=\"SCvH EOS H-He (rho=const.)\")\n",
    "\n",
    "# ds/drho(rho, T=const)\n",
    "ax[2][0].set_prop_cycle(None)\n",
    "for i in range(0, nT, nSkipT):\n",
    "    ax[2][0].semilogx(rho_table, s[i,:]/rho_table[i]*dlogs_dlogrho_int_rho[i](np.log10(rho_table)), '-')\n",
    "\n",
    "#for i in range(0, nT_int, nSkipT):\n",
    "#    ax[2][0].plot(logrho_int, dlogsdlogrho_int[i,:], '--')\n",
    "    \n",
    "ax[2][0].set(xlabel=\"rho [g cm^${-3}$]\", ylabel=\"ds/drho\")\n",
    "\n",
    "# ds/dT(rho=const, T)\n",
    "ax[2][1].set_prop_cycle(None)\n",
    "\n",
    "for i in range(0, nRho, nSkipRho):\n",
    "    ax[2][1].semilogx(T_table, s[:,i]/T_table*dlogs_dlogT_int_T[i](np.log10(T_table)), '-')\n",
    "\n",
    "#for i in range(0, nRho_int, nSkipRho):\n",
    "#    ax[2][1].plot(logT_int, dlogsdlogT_int[:,i], '--')\n",
    "\n",
    "ax[2][1].set(xlabel=\"T [K]\", ylabel=\"ds/dT\")\n",
    "ax[2][1].set(title=\"SCvH EOS H-He (rho=const.)\")\n"
   ]
  },
  {
   "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
}