plot_test_soundspeed.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Plot the sound speed of the SCvH equation of state"
   ]
  },
  {
   "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": "markdown",
   "metadata": {},
   "source": [
    "It seems that the sound speed $c_{s} \\left( \\rho, T \\right)$ is not provided for the SCvH EOS. It can be calculated from\n",
    "$$ c_{s} \\left( \\rho, T \\right) = \\left. \\frac{\\partial P}{\\partial \\rho} \\right|_{T} + \\frac{T}{\\rho^2 c_{v}} \\left( \\left. \\frac{\\partial P}{\\partial T} \\right|_{\\rho} \\right)^2 $$\n",
    "but that would require accurate derivatives. For the ideal (perfect) gas with\n",
    "$$ P \\left( \\rho, T \\right) = \\frac{k_B}{\\mu m_{H}} \\rho T $$\n",
    "and\n",
    "$$ u \\left( \\rho, T \\right) = c_V T$$\n",
    "with\n",
    "$$ c_V = \\frac{k_B}{\\left( \\gamma - 1 \\right) \\mu m_{H}} $$\n",
    "this expression results in\n",
    "$$ c_{s} \\left( \\rho, T \\right) = \\gamma \\left( \\gamma - 1 \\right) c_V T $$\n",
    "which is indeed the sound speed of an ideal (perfect gas)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def cs_ideal_gas(T, gamma=5.0/3.0, mu= 2.3):\n",
    "    kB = 1.38e-16  # erg/K\n",
    "    m_H = 1.67e-24 # g\n",
    "    \n",
    "    cv = kB/((gamma-1.0)*mu*m_H)\n",
    "    \n",
    "    return np.sqrt(gamma*(gamma-1.0)*cv*T)"
   ]
  },
  {
   "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)\n",
    "\n",
    "# Load sound speed for the eos table\n",
    "logcs = np.loadtxt(\"testscvheossoundspeed_table.txt\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load sound speed table\n",
    "logrho_int = np.loadtxt(\"testscvheossoundspeed_rhoaxis.txt\")\n",
    "logT_int = np.loadtxt(\"testscvheossoundspeed_taxis.txt\")\n",
    "\n",
    "logcs_int = np.loadtxt(\"testscvheossoundspeed.txt\")\n",
    "\n",
    "rho_int = 10**logrho_int\n",
    "T_int = 10**logT_int\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": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot log(cs) for the eos table and interpolated values\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",
    "# cs(rho, T=const)\n",
    "ax[0][0].set_prop_cycle(None)\n",
    "\n",
    "for i in range(0, nT, nSkipT):\n",
    "    ax[0][0].plot(logrho_table, logcs[i,:], '-')\n",
    "\n",
    "ax[0][0].set(xlabel=\"log(rho) [g cm^${-3}$]\", ylabel=\"log(cs) [cm s$^{-1}$]\")\n",
    "ax[0][0].set(title=\"SCvH EOS (H-He, T=const.)\")\n",
    "\n",
    "# cs(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, logcs[:,i], '-')\n",
    "\n",
    "ax[0][1].set(xlabel=\"log(T) [K]\", ylabel=\"log(cs) [cm s$^{-1}$]\")\n",
    "ax[0][1].set(title=\"SCvH EOS H-He (rho=const.)\")\n",
    "\n",
    "# cs(rho, T=const) for interpolated values\n",
    "ax[1][0].set_prop_cycle(None)\n",
    "\n",
    "for i in range(0, nT, nSkipT):\n",
    "    ax[1][0].plot(logrho_table, logcs[i,:], '-')\n",
    "\n",
    "for i in range(0, nT_int, nSkipT):\n",
    "    ax[1][0].plot(logrho_int, logcs_int[i,:], '--')\n",
    "    \n",
    "ax[1][0].set(xlabel=\"log(rho) [g cm^${-3}$]\", ylabel=\"log(cs) [cm s$^{-1}$]\")\n",
    "ax[1][0].set(title=\"SCvH EOS (H-He, T=const.)\")\n",
    "\n",
    "# cs(rho=const, T) for interpolated values\n",
    "ax[1][1].set_prop_cycle(None)\n",
    "\n",
    "for i in range(0, nRho, nSkipRho):\n",
    "    ax[1][1].plot(logT_table, logcs[:,i], '-')\n",
    "\n",
    "for i in range(0, nRho_int, nSkipRho):\n",
    "    ax[1][1].plot(logT_int, logcs_int[:,i], '--')\n",
    "\n",
    "ax[1][1].set(xlabel=\"log(T) [K]\", ylabel=\"log(cs) [cm s$^{-1}$]\")\n",
    "ax[1][1].set(title=\"SCvH EOS H-He (rho=const.)\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot log(cs) for the eos table\n",
    "fig, ax = plt.subplots(1, 2)\n",
    "\n",
    "x, y = fig.get_size_inches()\n",
    "\n",
    "fig.set_size_inches(2*x, y)\n",
    "\n",
    "# cs(rho, T=const)\n",
    "ax[0].set_prop_cycle(None)\n",
    "\n",
    "for i in range(0, nT, nSkipT):\n",
    "    ax[0].plot(logrho_table, logcs[i,:], '-')\n",
    "\n",
    "for i in range(0, nT, nSkipT):\n",
    "    ax[0].plot(logrho_table, np.log10(np.ones(nRho)*cs_ideal_gas(10**logT_table[i])), '--')\n",
    "    \n",
    "ax[0].set(xlabel=\"log(rho) [g cm^${-3}$]\", ylabel=\"log(cs) [cm s$^{-1}$]\")\n",
    "ax[0].set(title=\"SCvH EOS H-He, T=const.)\")\n",
    "\n",
    "# cs(rho=const., T)\n",
    "ax[1].set_prop_cycle(None)\n",
    "\n",
    "for i in range(0, nRho, nSkipRho):\n",
    "    ax[1].plot(logT_table, logcs[:,i], '-')\n",
    "\n",
    "for i in range(0, nRho, nSkipRho):\n",
    "    ax[1].plot(logT_table, np.log10(cs_ideal_gas(10**logT_table)), '--')\n",
    "\n",
    "ax[1].set(xlabel=\"log(T) [K]\", ylabel=\"log(cs) [cm s$^{-1}$]\")\n",
    "ax[1].set(title=\"SCvH EOS H-He (rho=const.)\")\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}