Creating tables in the CompOSE format I - Mandatory tables only

This section describes how a contributor can provide the mandatory files only (i.e. the grid parameter files eos.t, eos.nb and eos.yq and the table, eos.thermo, holding the thermodynamic quantities for a cold equation of state (EoS) in the CompOSE format.

Functionality to write general purpose tables will be provided in a future release of ComPyTools.

We consider two examples: in the first example, the baryon number density, mass-energy density and pressure are provided. In the second example, we show how to proceed if only the mass-energy density and the pressure are given.

For each of the tables, ComPyTools provides a set of user-friendly interfaces that allows the contributor to provide the required information in order to create the necessary files. With these interfaces, the user can perform the following steps:

  1. Feed the data into a data structure and provide units. This allows the code to perform sanity checking on the units, perform conversions if necessary and compute automatically any additional quantities expected by CompOSE,

  2. This data structure is then coverted into a tabulated form, allowing the user to view and inspect the data. This tabular form also makes it possible to write the data to file further down the line.

Given the thermodynamical quantities and the grid parameter values, we will then describe how we can automatically calculate the corresponding static neutron star properties.

Finally, once we describe how to create the tables, we create the PDF data sheet which summarises the EoS characteristics. This is a document that lives along side the tables on the CompOSE web site that provides a human-readable summary of the EoS.

Example 1 - Using baryon number density, pressure and energy-density as inputs

In this example, we consider an EoS in the Lorene format. The file describes the APR(APR) EoS and contains the following data

Column number

Quantity

Unit

1

Baryon number density

fm\(^{-3}\)

2

Mass-energy density/\(c^{2}\)

g cm\(^{-3}\)

3

Pressure

dyn cm\(^{-2}\)

Let’s first dowload the data and then load it:

[1]:

import os
import tempfile
import requests
import warnings
warnings.filterwarnings("ignore", "Wswiglal-redir-stdio")

import numpy as np

tmpdir = tempfile.TemporaryDirectory()
eos_file = "eos_akmalpr.d"
eos_url = "https://gitlab.in2p3.fr/lpc-caen/compytools/-/raw/main/tests/inputs/lorene/"
source = os.path.join(eos_url, eos_file)
dest = os.path.join(tmpdir.name, eos_file)

response = requests.get(source)
with open(dest, 'wb') as context:
    context.write(response.content)

For the purposes of this tutorial, to avoid polluting the file system with output files, we download the files into a temporary directory. However, if you are following along and you wish to keep the files, you can replace tmpdir.name with the directory path of your choice.

[2]:

nbgrid, rho, p = np.loadtxt(
    dest,
    # Ignore header info
    skiprows=9,
    usecols=(1, 2, 3),
    unpack=True
)

We first describe how to create the grid parameter tables.

Grid parameter files (eos.[t,nb,yq])

[3]:

from compytools import (
    Grid,
    GridData
)

# Load useful units for our input data
from compytools.constants import (
  DIMENSIONLESS,
  MEV,
  PERCM3,
  PERFM3
)

Store the grid points into a data structure via the GridData class. For cold tables there is only a single value for temperature and charge fraction, which are both zero:

[4]:

grid_data = GridData(
    # Temperature grid points in MeV
    t=np.array([0.0]) << MEV,
    # Baryon number density in 1/fm^3
    nb=nbgrid << PERFM3,
    # Charge fraction (no units)
    yq=np.array([0.0]) << DIMENSIONLESS
)

It is mandatory to provide units even if they are already in the units as required by CompOSE. This not only allows ComPyTools to perform sanity checking on the input data, but it also allows units to be added to the table’s metadata.

Here and elsewhere, we use the syntax << (or alternatively *) to provide units for each provided quantity. Next, create tables for each of the grid parameters using the Grid class.

[5]:

grid_params = Grid.from_dataset(grid_data)

We can access the grid points for the baryon number. Note how the units are attached to the data array.

[6]:

grid_params.nb.grid_points

[6]:
$[7.9240596 \times 10^{-15},~2.5058077 \times 10^{-14},~7.9240596 \times 10^{-14},~\dots,~1.3,~1.32,~1.34] \; \mathrm{\frac{1}{fm^{3}}}$

The syntax to access grid points for temperature t and charge fraction yq is similar.

Thermodynamic quantities (eos.thermo)

A similar procedure as described above for the grid parameter data is applicable for the remaining tables.

Here we provide the data for the eos.thermo table.

[7]:

from compytools import (
    ExtraThermoQuantity,
    Nucleons,
    Thermo,
    ThermoData
)

from compytools.constants import (
    CLIGHT_CGS_SQ,
    DYN_PER_CM2,
    ERG,
    GRAM,
    G_PER_CM3,
    KELVIN,
    MEV,
    MEV_PER_FM3
)

The header information that specifies the neutron and proton masses, and the flag indicating whether leptons are included is provided via the Nucleons class:

[8]:

mn = 1.6749286e-24 << GRAM
mp = 1.6726231e-24 << GRAM

nucleons = Nucleons(mn=mn, mp=mp, has_leptons=1)

print(nucleons)

Nucleons(has_leptons=1, mn=<Quantity 939.56603867 MeV>, mp=<Quantity 938.27274802 MeV>)

Note that the nucleon masses have been automatically converted in the required units of MeV.

Now let’s provide the thermodynamic data. Here, we will also demonstrate how to provide an additional thermodynamic quantity by providing the square of the speed of sound, \(c_{s}^{2}\) which will be in units of the square of the speed of light, \(c\). Since we take \(c=1\), \(c_{s}^{2}\) will be dimensionless.

To calculate \(c_{s}^{2}\), we first need to convert the pressure and the mass-energy density into appropriate units of MeV fm\(^{-3}\). astropy allows us to do this easily via the to() method:

[9]:

import astropy.units as u

total_pressure = (p << DYN_PER_CM2).to(MEV_PER_FM3)
rho = rho << G_PER_CM3

energy_density = rho.to(MEV_PER_FM3, equivalencies=u.mass_energy())
cs2 = np.gradient(total_pressure, energy_density, edge_order=2)

Note that the equivalencies=u.mass_energy() argument allows us to convert units into relativistic units. For further information, the interested user is referred to the astropy documentation.

Additional quantities are provided via a dictionary:

  • The key of the dictionary gives the variable name. You can use \(\LaTeX\) syntax which will allow the quantity name to be correctly displayed within the data sheet,

  • The value of the key is a class ExtraThermoQuantity where you provide the array, units and a description. The latter will be added to create the data sheet.

[10]:

extra_thermo = {
    '$c_{s}^{2}$': ExtraThermoQuantity(
        data=cs2,
        unit=DIMENSIONLESS,
        description='Square of the sound speed in units of $c^2$')
}

Now, let’s create the thermodynamic data structure:

[11]:

mub = (total_pressure + energy_density) / grid_params.nb.grid_points
thermo_data = ThermoData(
    # Provide the grid parameter provided above
    grid_params=grid_params,
    # Provide the nucleon masses provided above.
    nucleons=nucleons,
    pressure=total_pressure,
    entropy_density=np.zeros(nbgrid.size) << PERFM3,
    baryon_chemical_potential=mub,
    charge_chemical_potential=np.zeros(nbgrid.size) << MEV,
    lepton_chemical_potential=np.zeros(nbgrid.size) << MEV,
    free_energy_density=energy_density,
    internal_energy_density=energy_density,
    # Additional quantities provided via the `Qextra` argument.
    # This can be omitted otherwise.
    Qextra=extra_thermo
)

# Now create the thermo table
thermo = Thermo.from_dataset(thermo_data)

Also note that if the baryon chemical potential, \(\mu_{\mathrm{B}}\), is not known, the baryon_chemical_potential argument above can be omitted, in which case it will be automatically calculated via

\[\mu_{\mathrm{B}}=\frac{p+e}{n_{\mathrm{B}}}\]

for the case of a cold EoS. Here, \(p\) is the pressure, \(e\) is the internal energy density and \(n_{\mathrm{B}}\) is the baryon number density.

The data, as well as the associated metadata, can be accessed with the respective methods:

[12]:

# max_width specifies the maximum number of characters to use in the displayed table. If a value of -1 is used,
# then no maximum value is imposed.
thermo.pprint(max_width=-1)
thermo.info

 iT jnB kYq     Q1         Q2          Q3         Q4         Q5          Q6          Q7     Nextra $c_{s}^{2}$
               MeV
--- --- --- ---------- ---------- ----------- ---------- ---------- ----------- ----------- ------ -----------
  0   1   1 2.3807e-05 0.0000e+00 -8.9130e-03 0.0000e+00 0.0000e+00 -8.9130e-03 -8.9130e-03      1  4.8046e-08
  0   2   1 4.8003e-05 0.0000e+00 -8.9129e-03 0.0000e+00 0.0000e+00 -8.9130e-03 -8.9130e-03      1  7.9090e-08
  0   3   1 9.6792e-05 0.0000e+00 -8.9128e-03 0.0000e+00 0.0000e+00 -8.9129e-03 -8.9129e-03      1  1.5947e-07
  0   4   1 1.9517e-04 0.0000e+00 -8.9125e-03 0.0000e+00 0.0000e+00 -8.9127e-03 -8.9127e-03      1  3.2156e-07
  0   5   1 3.9353e-04 0.0000e+00 -8.9119e-03 0.0000e+00 0.0000e+00 -8.9124e-03 -8.9124e-03      1  6.4838e-07
  0   6   1 7.9350e-04 0.0000e+00 -8.9108e-03 0.0000e+00 0.0000e+00 -8.9117e-03 -8.9117e-03      1  1.3074e-06
  0   7   1 1.6000e-03 0.0000e+00 -8.9085e-03 0.0000e+00 0.0000e+00 -8.9102e-03 -8.9102e-03      1  2.6361e-06
  0   8   1 3.2262e-03 0.0000e+00 -8.9040e-03 0.0000e+00 0.0000e+00 -8.9074e-03 -8.9074e-03      1  5.3153e-06
  0   9   1 6.5051e-03 0.0000e+00 -8.8948e-03 0.0000e+00 0.0000e+00 -8.9017e-03 -8.9017e-03      1  1.0717e-05
... ... ...        ...        ...         ...        ...        ...         ...         ...    ...         ...
  0 163   1 9.5570e+02 0.0000e+00  1.4679e+00 0.0000e+00 0.0000e+00  4.5071e-01  4.5071e-01      1  1.3312e+00
  0 164   1 9.9224e+02 0.0000e+00  1.5245e+00 0.0000e+00 0.0000e+00  4.6843e-01  4.6843e-01      1  1.3434e+00
  0 165   1 1.0301e+03 0.0000e+00  1.5829e+00 0.0000e+00 0.0000e+00  4.8651e-01  4.8651e-01      1  1.3643e+00
  0 166   1 1.0683e+03 0.0000e+00  1.6420e+00 0.0000e+00 0.0000e+00  5.0497e-01  5.0497e-01      1  1.3824e+00
  0 167   1 1.1076e+03 0.0000e+00  1.7027e+00 0.0000e+00 0.0000e+00  5.2379e-01  5.2379e-01      1  1.3908e+00
  0 168   1 1.1465e+03 0.0000e+00  1.7632e+00 0.0000e+00 0.0000e+00  5.4298e-01  5.4298e-01      1  1.4080e+00
  0 169   1 1.1873e+03 0.0000e+00  1.8262e+00 0.0000e+00 0.0000e+00  5.6254e-01  5.6254e-01      1  1.4326e+00
  0 170   1 1.2282e+03 0.0000e+00  1.8897e+00 0.0000e+00 0.0000e+00  5.8247e-01  5.8247e-01      1  1.4372e+00
  0 171   1 1.2696e+03 0.0000e+00  1.9540e+00 0.0000e+00 0.0000e+00  6.0276e-01  6.0276e-01      1  1.4513e+00
  0 172   1 1.3118e+03 0.0000e+00  2.0196e+00 0.0000e+00 0.0000e+00  6.2342e-01  6.2342e-01      1  1.4713e+00
Length = 172 rows

[12]:

<Thermo length=172>
    name     dtype  unit format                    description                      class
----------- ------- ---- ------ -------------------------------------------------- --------
         iT   int32          %d                             Temperature grid point   Column
        jnB   int32          %d                     Baryon num. density grid point   Column
        kYq   int32          %d                    Charge num. fraction grid point   Column
         Q1 float32  MeV   %.4e       pressure divided by baryon num. density p/nb Quantity
         Q2 float32        %.4e                            entropy per baryon s/nb Quantity
         Q3 float32        %.4e scaled and shifted baryon chem. potential mub/mn-1 Quantity
         Q4 float32        %.4e            scaled charge chemical potential muq/mn Quantity
         Q5 float32        %.4e  scaled effective lepton chemical potential mul/mn Quantity
         Q6 float32        %.4e      scaled free energy per baryon f/(nb x mn) - 1 Quantity
         Q7 float32        %.4e  scaled internal energy per baryon e/(nb x mn) - 1 Quantity
     Nextra   int32          %d                num. extra thermodynamic quantities   Column
$c_{s}^{2}$ float32        %.4e        Square of the sound speed in units of $c^2$ Quantity

If you are using this functionality in your scripts you will need to use Python’s print() function, i.e. print(thermo.pprint(max_width=-1)) and print(thermo.info).

Data can be accessed (for example for plotting) as described in the Section Reading CompOSE tables.

[13]:

import matplotlib.pyplot as plt

plt.figure()

nb = grid_params.nb.grid_points
pressure = nb * thermo['Q1']
plt.loglog(nb, pressure, '-k', linewidth=2)

plt.xlabel(r'$n_{B}$ [fm$^{-3}$]')
plt.ylabel(r'$p$ [MeV fm$^{-3}$]')
plt.grid()
_images/write-eos-mandatory-tables_33_0.png

Neutron star properties

ComPyTools includes the lalsimulation package that is part of the lalsuite library. The lalsimulation package includes a solver for the Tolman-Oppenheimer-Volkoff equations that allows the user to automatically calculate several static neutron star properties based on the input EoS. The quantities are:

  • The radius in km,

  • The gravitational mass in solar masses,

  • The quadrupole tidal deformability (adim, multipole order \(\ell=2\)),

  • The central baryon number density in fm\(^{-3}\),

  • The baryonic mass in solar masses,

  • The octupole tidal deformability (adim, \(\ell=3\)),

  • The hexadecapole tidal deformability (adim, \(\ell=4\)).

These quantities are calculated via the NeutronStarData class by supplying the grid parameters and the thermodynamical table. We then create the corresponding mr data table.

[14]:

from compytools import (
    NeutronStar,
    NeutronStarData,
)

ns_data = NeutronStarData(
    grid_params=grid_params,
    thermo=thermo
)

mr = NeutronStar.from_dataset(ns_data)
mr.pprint(max_width=-1)
mr.info



[11:35:48] INFO     Max. mass: 2.189 Msun                                                                 mr.py:501
                    Radius at max. mass: 9.944 Msun                                                                
                    Radius at 1.4 Msun: 11.370 km                                                                  
                    Lambda_2 at 1.4 Msun: 251.211                                                                  
                    Lambda GW170817: 294.821                                                                       

  radius   mass_grav   Lambda2   central_baryon_num_density mass_baryonic  Lambda3    Lambda4
    km      solMass                       1 / fm3              solMass
---------- ---------- ---------- -------------------------- ------------- ---------- ----------
1.5525e+01 2.1532e-01 3.3481e+06                 2.0573e-01    2.1597e-01 2.3027e+08 1.6920e+10
1.5437e+01 2.1729e-01 3.2046e+06                 2.0678e-01    2.1799e-01 2.1609e+08 1.5484e+10
1.5349e+01 2.1927e-01 3.0644e+06                 2.0782e-01    2.2000e-01 2.0245e+08 1.4137e+10
1.5262e+01 2.2125e-01 2.9298e+06                 2.0885e-01    2.2202e-01 1.8959e+08 1.2906e+10
1.5178e+01 2.2322e-01 2.8024e+06                 2.0987e-01    2.2404e-01 1.7769e+08 1.1803e+10
1.5097e+01 2.2520e-01 2.6832e+06                 2.1089e-01    2.2606e-01 1.6681e+08 1.0827e+10
1.5020e+01 2.2717e-01 2.5706e+06                 2.1190e-01    2.2808e-01 1.5676e+08 9.9474e+09
1.4945e+01 2.2915e-01 2.4640e+06                 2.1290e-01    2.3011e-01 1.4742e+08 9.1535e+09
1.4872e+01 2.3113e-01 2.3628e+06                 2.1390e-01    2.3213e-01 1.3872e+08 8.4366e+09
       ...        ...        ...                        ...           ...        ...        ...
1.0279e+01 2.1715e+00 3.9732e+00                 1.0197e+00    2.6271e+00 1.5454e+00 5.5601e-01
1.0260e+01 2.1734e+00 3.8588e+00                 1.0265e+00    2.6303e+00 1.4903e+00 5.3191e-01
1.0241e+01 2.1754e+00 3.7427e+00                 1.0336e+00    2.6335e+00 1.4349e+00 5.0787e-01
1.0219e+01 2.1774e+00 3.6217e+00                 1.0415e+00    2.6368e+00 1.3777e+00 4.8330e-01
1.0196e+01 2.1794e+00 3.4962e+00                 1.0501e+00    2.6400e+00 1.3190e+00 4.5834e-01
1.0170e+01 2.1814e+00 3.3637e+00                 1.0598e+00    2.6433e+00 1.2576e+00 4.3253e-01
1.0140e+01 2.1833e+00 3.2211e+00                 1.0709e+00    2.6466e+00 1.1924e+00 4.0545e-01
1.0103e+01 2.1853e+00 3.0598e+00                 1.0847e+00    2.6498e+00 1.1197e+00 3.7566e-01
1.0056e+01 2.1873e+00 2.8672e+00                 1.1026e+00    2.6531e+00 1.0344e+00 3.4129e-01
9.9438e+00 2.1893e+00 2.4803e+00                 1.1463e+00    2.6564e+00 8.6817e-01 2.7627e-01
Length = 1000 rows

[14]:

<NeutronStar length=1000>
           name             dtype    unit  format           description             class
-------------------------- ------- ------- ------ -------------------------------- --------
                    radius float64      km   %.4e              Neutron star radius Quantity
                 mass_grav float64 solMass   %.4e  Neutron star gravitational mass Quantity
                   Lambda2 float64           %.4e   Quadrupole tidal deformability Quantity
central_baryon_num_density float64 1 / fm3   %.4e    Central baryon number density Quantity
             mass_baryonic float64 solMass   %.4e       Neutron star baryonic mass Quantity
                   Lambda3 float64           %.4e     Octupole tidal deformability Quantity
                   Lambda4 float64           %.4e Hexadecapole tidal deformability Quantity

[15]:

plt.figure(figsize=(10, 3))
plt.subplot(121)
plt.plot(mr['radius'], mr['mass_grav'], label='Gravitational mass')
plt.plot(mr['radius'], mr['mass_baryonic'], label='Baryonic mass')
plt.legend()
plt.grid()
plt.xlabel(r'$R$ [km]')
plt.ylabel(r'$M$ [M$_{\odot}$]')

plt.subplot(122)
plt.semilogy(mr['mass_grav'], mr['Lambda2'], label=r'$\Lambda_{2}$')
plt.semilogy(mr['mass_grav'], mr['Lambda3'], label=r'$\Lambda_{3}$')
plt.semilogy(mr['mass_grav'], mr['Lambda4'], label=r'$\Lambda_{4}$')
plt.legend()
plt.grid()
plt.xlabel(r'$M$ [M$_{\odot}]$')
plt.ylabel(r'$\Lambda$')

[15]:

Text(0, 0.5, '$\\Lambda$')
_images/write-eos-mandatory-tables_36_1.png

Writing the CompOSE tables

Once we have created a data table associated to each of the CompOSE tables, we now combine the data together using the EoS class. This allows ComPyTools to perform quality checks on the data:

  • Whether the data contains invalid values such as NaNs, infinities or values that are not floats,

  • Sound speed checks (acausality and negative values),

  • Monotonicity checks on pressure and enthalpy, as well as for the grid point values,

  • Whether the entropy and lepton chemical potentials are indeed zero for cold tables.

[16]:

import tempfile

from compytools import EoS

# For the purposes of this tutorial create a temporary directory into which
# to write the tables.
tmpdir = tempfile.TemporaryDirectory()

eos = EoS(
    params=grid_params,
    thermo=thermo,
    mr=mr
)




[11:35:49] INFO     Checking monotonicity in pressure                                                    eos.py:413

           INFO     Checking monotonicity in enthalpy_per_baryon                                         eos.py:413

           INFO     Checking monotonicity in baryon_num_density                                          eos.py:413

           INFO     Checking sound speed causality                                                       eos.py:384

           WARNING  Non-causality found from grid point: 148                                           utils.py:149

           INFO     Checking for negative sound speed.                                                   eos.py:386

           INFO     Checking quantity entropy_per_baryon.                                                eos.py:390

           INFO     Checking quantity lepton_chem_potential.                                             eos.py:390

Next we use the to_compose method to write all of the tables:

[17]:

eos.to_compose(tmpdir.name)

# List the contents of the directory
os.listdir(tmpdir.name)

           INFO     Writing CompOSE tables to /tmp/tmpuu4v1wro                                           eos.py:533

[17]:

['eos.t', 'eos.thermo', 'eos.yq', 'eos.init', 'eos.nb', 'eos.mr']

Creating the data sheet

In addition to the CompOSE tables, the contributor should also provide a data sheet. Python interfaces have been developed that will allow the contributor to provide the required information, from which the PDF document will be automatically created from a \(\LaTeX\) template. Each of the required information will be described in turn.

A brief description of the EoS, such as key physical inputs along with any methods that were used to create it, can be provided in the abstract. Note that you can provide \(\LaTeX\) symbols by using $ symbol, while \(\LaTeX\) commands must be proceeded with an additional backslash \ so that Python does not try and interpret these as Python syntax, for example \\textbf instead of \textbf. Supplementary references can be provided, but this is not mandatory.

[18]:

from compytools import (
    DataSheet,
    SubmissionDetails,
    NuclearMatterProperties,
    NeutronStarProperties
)

abstract = """
This table represents the zero temperature and beta--equilibrium EoS by Akmal,
Pandharipande and Ravenhall using variational techniques [1],
interaction A18 + delta v + UIX*. The inner crust is calculated with SLy4 [3],
the outer crust from Baym, Pethick, Sutherland [2]. No compositional information is available.
"""

references = {
    1: "A. Akmal, V. R. Pandharipande and D. G. Ravenhall, Phys. Rev. C 58, 1804 (1998)",
    2: "G. Baym, C. Pethick and P. Sutherland, Astrophys. J. 170, 299 (1971)",
    3: "F. Douchin, P. Haensel, Astronomy and Astrophysics 380, 151 (2001)"
}

further_refs = {
    1: "Shapiro, Teukolsky, The Physics of Neutron Stars, 1983"
}

The contributor can provide details of the submitting author as well as the name of the EoS as follows:

[19]:

submission_details = SubmissionDetails(
    eos_name='APR(APR)',
    category='Nucleonic',
    submitted_by='Albert Einstein',
    affiliation='University of Princeton',
    email='einstein@email.com',
    abstract=abstract,
)

while the nuclear empirical parameters associated with the EoS should be provided via the NuclearMatterProperties class:

[20]:

nuclear_params = NuclearMatterProperties(
    saturation_density_ns=0.16,
    binding_energy_per_baryon_at_sat_E0=16.0,
    isoscalar_incompressibility_modulus_K=266.0,
    # Indicate unavailable values with None
    minus_isoscalar_skewness_Kprime=None,
    symmetry_energy_J=-32.6,
    symmetry_energy_slope_param_L=57.6,
    symmetry_incompressibility_Ksym=None
)

The arguments to NuclearMatterProperties correspond to the following mathematical notation:

Argument

Notation

saturation_density_ns

\(n_{s}\)

binding_energy_per_baryon_at_sat_E0

\(E_{0}\)

isoscalar_incompressibility_modulus_K

\(K\)

minus_isoscalar_skewness_Kprime

\(K^{\prime}\)

symmetry_energy_J

\(J\)

symmetry_energy_slope_param_L

\(L\)

symmetry_incompressibility_Ksym

\(K_{sym}\)

CompOSE uses the convention that \(E_{0}\) should be positive. If possible, the user should provide the isoscalar parameters to at least second order. If either of these criteria are not met, then ComPyTools will show a warning message.

Finally, the information is combined using the DataSheet class and the PDF is generated using the .write() method. Note that the neutron star static properties can be provided via ns_data that was generated previously via the NeutronStarData class above.

[21]:

sheet = DataSheet(
    eos=eos,
    submission_details=submission_details,
    nuclear_matter=nuclear_params,
    neutron_star_data=ns_data,
    references=references,
    further_references=further_refs
)

content = sheet.write(tmpdir.name)

[11:35:51] INFO     Successfully wrote data sheet to: /tmp/tmpuu4v1wro.                            datasheet.py:340

The further_references argument is not mandatory and can be omitted if it is not used.

[22]:

from pdf2image import convert_from_path
from IPython.display import display

# This allows us to display the PDF data sheet within the notebook documentation.
datasheet_path = os.path.join(tmpdir.name, 'datasheet.pdf')
pages = convert_from_path(datasheet_path, dpi=150)
for page in pages:
    display(page)
_images/write-eos-mandatory-tables_52_0.png
_images/write-eos-mandatory-tables_52_1.png
_images/write-eos-mandatory-tables_52_2.png
_images/write-eos-mandatory-tables_52_3.png

Once all of the CompOSE files and the data sheet have been created, you can arrange to have the data added to the CompOSE web site by contacting the CompOSE core team at the following e-mail address: develop.compose@obspm.fr.

Example 2 - Pressure and energy density as input

In this example, we consider an alternative scenario where only the pressure and the energy density are provided (for example as is done for the RNS format). Indeed, we will use an RNS file as an example which contains the following quantities:

Column number

Quantity

Unit

1

Energy density

g cm\(^{-3}\)

2

Pressure

dyn cm\(^{-2}\)

Let’s download the data and then load it. For the purposes of the tutorial we will download the data into a temporary directory. However, if you wish to follow along, you can download the data into a folder of your choice on your computer.

[23]:

import os
import tempfile
import requests

import numpy as np

tmpdir = tempfile.TemporaryDirectory()
eos_file = "eosAPR.txt"
eos_url = "https://gitlab.in2p3.fr/lpc-caen/compytools/-/raw/main/tests/inputs/RNS"
source = os.path.join(eos_url, eos_file)
dest = os.path.join(tmpdir.name, eos_file)

response = requests.get(source)
with open(dest, 'wb') as context:
    context.write(response.content)

[24]:

rho, p = np.loadtxt(
    dest,
    # Ignore header info
    skiprows=1,
    usecols=(0, 1),
    unpack=True
)

Grid parameter files (eos.[t,nb,yq])

Since we don’t already have the baryon number density \(n_{\mathrm{B}}\), this can be recomputed from the pressure, \(p\) and the energy density \(e\) via:

\[n_{\mathrm{B}}=\frac{(e+p)}{m_{u}}\exp(-h)\]

where \(m_u\) is the atomic mass unit and \(h\) is the pseudo-enthalpy, which can be calculated via

\[h=\int_{0}^{p}\frac{1}{e+p^{\prime}}\mathrm{d}p^{\prime}.\]
[25]:

import astropy.units as u

import compytools.utils

from compytools.constants import (
    CLIGHT_CGS_SQ,
    DYN_PER_CM2,
    G_PER_CM3,
    MEV_PER_FM3
)

# Provide the units to the tabulated quantities so that they can be converted
# to the appropriate units. This will be needed for the CompOSE eos.thermo file
# described below.
# Units are provided via the << syntax.
p = p << DYN_PER_CM2
rho = rho << G_PER_CM3

# Convert to the required units via astropy's .to() method
p = p.to(MEV_PER_FM3)
e = rho.to(MEV_PER_FM3, equivalencies=u.mass_energy())

Note that the equivalencies=u.mass_energy() argument allows us to convert units into relativistic units. For further information, the interested user is referred to the astropy documentation.

Now calculate the baryon number density via the calculate_baryon_num_density utility function:

[26]:

nbgrid = compytools.utils.calculate_baryon_num_density(e, p)

# Note that the resulting points in baryon number density automatically have
# the appropriate units of 1/fm^3
print(nbgrid.unit)

           INFO     Checking non-monotonicity for log-enthalpy.                                        utils.py:305

           WARNING  Re-computing log-enthalpy via trapezoidal method.                                  utils.py:307

1 / fm3

The integration is initially performed using Simpson’s method. However, if the resulting log enthalpy contains non-monotonic points, then we re-perform the integration with the trapezoidal method.

Now we can store the grid points into a data structure via the GridData class. For cold tables there is only a single value for temperature and charge fraction, which are both zero:

[27]:

from compytools import (
    GridData,
    Grid
)

from compytools.constants import (
  DIMENSIONLESS,
  MEV,
  PERFM3
)

grid_data = GridData(
    # Temperature grid points in MeV
    t=np.array([0.0]) << MEV,
    # Units are already provided automatically, so no need to
    # re-supply here.
    nb=nbgrid,
    # Charge fraction (no units)
    yq=np.array([0.0]) << DIMENSIONLESS
)

Here and elsewhere, it is mandatory to provide units even if they are already in the units as required by CompOSE. This not only allows ComPyTools to perform sanity checking on the input data, but it also allows units to be added to the table’s metadata.

Now create tables for each of the grid parameters using the Grid class:

[28]:

grid_params = Grid.from_dataset(grid_data)

Grid points for a given parameter, for example the baryon number density, can be accessed with the following syntax:

[29]:

grid_params.nb.grid_points

[29]:
$[4.7397319 \times 10^{-15},~4.7616103 \times 10^{-15},~4.9176079 \times 10^{-15},~\dots,~0.52686042,~0.82088697,~1.1074668] \; \mathrm{\frac{1}{fm^{3}}}$

Grid points for temperature and charge number fraction can be accessed similarly with grid_params.t.grid_points and grid_params.yq.grid_points.

Thermodynamic quantities (eos.thermo)

Now we provide the data needed to create the eos.thermo table:

[30]:

from compytools import (
    Nucleons,
    Thermo,
    ThermoData
)

from compytools.constants import (
    MNEUTRON,
    MPROTON
)

[31]:

# Provide the header information giving the
# proton and neutron masses and whether
# leptons are considered.

# For the masses, we use those provided by CODATA. These already have astropy units attached to them.
nucleons = Nucleons(
    mn=MNEUTRON,
    mp=MPROTON,
    has_leptons=0
)

# Note that if the units for quantities have been specified earlier
# (e.g. pressure) they do not need to be given again in this interface.
thermo_data = ThermoData(
    nucleons=nucleons,
    grid_params=grid_params,
    pressure=p,
    # Zero for cold tables
    entropy_density=np.zeros((nbgrid.size, )) << PERFM3,
    # Zero for cold tables
    charge_chemical_potential=np.zeros((nbgrid.size, )) << MEV,
    # Zero for cold tables
    lepton_chemical_potential=np.zeros((nbgrid.size, )) << MEV,
    # Free and internal energies are equal for cold matter
    free_energy_density=e,
    internal_energy_density=e
)

Now create the table:

[32]:

thermo = Thermo.from_dataset(thermo_data)

The table and a description of the columns can be shown with the respective commands:

[33]:

thermo.pprint(max_width=-1)
thermo.info

 iT jnB kYq     Q1         Q2          Q3         Q4         Q5          Q6          Q7     Nextra
               MeV
--- --- --- ---------- ---------- ----------- ---------- ---------- ----------- ----------- ------
  0   1   1 1.3312e-10 0.0000e+00 -8.5905e-03 0.0000e+00 0.0000e+00 -8.5905e-03 -8.5905e-03      0
  0   2   1 1.3251e-09 0.0000e+00 -8.5905e-03 0.0000e+00 0.0000e+00 -8.5905e-03 -8.5905e-03      0
  0   3   1 1.2831e-08 0.0000e+00 -8.5905e-03 0.0000e+00 0.0000e+00 -8.5905e-03 -8.5905e-03      0
  0   4   1 1.0796e-07 0.0000e+00 -8.5905e-03 0.0000e+00 0.0000e+00 -8.5905e-03 -8.5905e-03      0
  0   5   1 8.8358e-07 0.0000e+00 -8.5905e-03 0.0000e+00 0.0000e+00 -8.5905e-03 -8.5905e-03      0
  0   6   1 3.9016e-06 0.0000e+00 -8.5905e-03 0.0000e+00 0.0000e+00 -8.5905e-03 -8.5905e-03      0
  0   7   1 2.8395e-05 0.0000e+00 -8.5904e-03 0.0000e+00 0.0000e+00 -8.5904e-03 -8.5904e-03      0
  0   8   1 1.7110e-04 0.0000e+00 -8.5898e-03 0.0000e+00 0.0000e+00 -8.5900e-03 -8.5900e-03      0
  0   9   1 9.6660e-04 0.0000e+00 -8.5847e-03 0.0000e+00 0.0000e+00 -8.5858e-03 -8.5858e-03      0
... ... ...        ...        ...         ...        ...        ...         ...         ...    ...
  0  92   1 5.3225e+00 0.0000e+00  1.8585e-02 0.0000e+00 0.0000e+00  1.2920e-02  1.2920e-02      0
  0  93   1 6.1254e+00 0.0000e+00  1.9850e-02 0.0000e+00 0.0000e+00  1.3331e-02  1.3331e-02      0
  0  94   1 6.7856e+00 0.0000e+00  2.0883e-02 0.0000e+00 0.0000e+00  1.3661e-02  1.3661e-02      0
  0  95   1 8.7036e+00 0.0000e+00  2.3882e-02 0.0000e+00 0.0000e+00  1.4619e-02  1.4619e-02      0
  0  96   1 9.9081e+00 0.0000e+00  2.5761e-02 0.0000e+00 0.0000e+00  1.5215e-02  1.5215e-02      0
  0  97   1 1.1143e+01 0.0000e+00  2.7694e-02 0.0000e+00 0.0000e+00  1.5834e-02  1.5834e-02      0
  0  98   1 3.8326e+01 0.0000e+00  7.6409e-02 0.0000e+00 0.0000e+00  3.5619e-02  3.5619e-02      0
  0  99   1 1.8007e+02 0.0000e+00  3.2486e-01 0.0000e+00 0.0000e+00  1.3320e-01  1.3320e-01      0
  0 100   1 6.9951e+02 0.0000e+00  1.1991e+00 0.0000e+00 0.0000e+00  4.5462e-01  4.5462e-01      0
  0 101   1 2.7447e+03 0.0000e+00  4.6167e+00 0.0000e+00 0.0000e+00  1.6955e+00  1.6955e+00      0
Length = 101 rows

[33]:

<Thermo length=101>
 name   dtype  unit format                    description                      class
------ ------- ---- ------ -------------------------------------------------- --------
    iT   int32          %d                             Temperature grid point   Column
   jnB   int32          %d                     Baryon num. density grid point   Column
   kYq   int32          %d                    Charge num. fraction grid point   Column
    Q1 float32  MeV   %.4e       pressure divided by baryon num. density p/nb Quantity
    Q2 float32        %.4e                            entropy per baryon s/nb Quantity
    Q3 float32        %.4e scaled and shifted baryon chem. potential mub/mn-1 Quantity
    Q4 float32        %.4e            scaled charge chemical potential muq/mn Quantity
    Q5 float32        %.4e  scaled effective lepton chemical potential mul/mn Quantity
    Q6 float32        %.4e      scaled free energy per baryon f/(nb x mn) - 1 Quantity
    Q7 float32        %.4e  scaled internal energy per baryon e/(nb x mn) - 1 Quantity
Nextra   int32          %d                num. extra thermodynamic quantities   Column

Accessing the required data can be done as already described in the Section on writing the eos.thermo file in Example 1. Writing out all of the CompOSE and creating the data sheet are respectively described in Section Writing the CompOSE tables and Creating the data sheet.

Once all of the CompOSE files and the data sheet have been created, you can arrange to have the data added to the CompOSE web site by contacting the CompOSE core team at the following e-mail address: develop.compose@obspm.fr.