From 0085a724986ed4628d3cfccde44805a2f7db389a Mon Sep 17 00:00:00 2001 From: ChunHuangPhy Date: Mon, 11 Nov 2024 15:49:57 -0600 Subject: [PATCH] add new documentation files --- .github/workflows/documentation.yaml | 2 +- docs/source/MITbag_EOS.py | 21 + docs/source/MITbag_EOS.rst | 9 + docs/source/Polytrope_EOS.py | 88 ++++ docs/source/Polytrope_EOS.rst | 9 + docs/source/RMF_DDH.py | 635 +++++++++++++++++++++++++++ docs/source/RMF_DDH.rst | 9 + docs/source/SpeedofSound_EOS.py | 106 +++++ docs/source/SpeedofSound_EOS.rst | 9 + docs/source/Strangeon_EOS.py | 36 ++ docs/source/Strangeon_EOS.rst | 9 + 11 files changed, 932 insertions(+), 1 deletion(-) create mode 100644 docs/source/MITbag_EOS.py create mode 100644 docs/source/MITbag_EOS.rst create mode 100644 docs/source/Polytrope_EOS.py create mode 100644 docs/source/Polytrope_EOS.rst create mode 100644 docs/source/RMF_DDH.py create mode 100644 docs/source/RMF_DDH.rst create mode 100644 docs/source/SpeedofSound_EOS.py create mode 100644 docs/source/SpeedofSound_EOS.rst create mode 100644 docs/source/Strangeon_EOS.py create mode 100644 docs/source/Strangeon_EOS.rst diff --git a/.github/workflows/documentation.yaml b/.github/workflows/documentation.yaml index ffaa3d5..7830dad 100644 --- a/.github/workflows/documentation.yaml +++ b/.github/workflows/documentation.yaml @@ -37,7 +37,7 @@ jobs: pip install matplotlib pip install ultranest pip install NumbaMinpack - pip install CompactObject-TOV==1.8.1 + pip install CompactObject-TOV - name : install docs dependencies run: | pip install decorator h5py diff --git a/docs/source/MITbag_EOS.py b/docs/source/MITbag_EOS.py new file mode 100644 index 0000000..e6e29f8 --- /dev/null +++ b/docs/source/MITbag_EOS.py @@ -0,0 +1,21 @@ +import numpy as np +from TOVsolver.unit import g_cm_3, dyn_cm_2, MeV, fm + +def MITbag_compute_EOS(B): + """ + Compute the energy density and pressure based on the given parameters. + + Args: + B: Input value of bag constant; MeVfm^-3 + + Returns: + tuple: Arrays of energy densities in units of gcm^3 and pressures in units of dyncm^2. + """ + + B_cgs = B * (MeV / (fm)**3) # converting input to cgs + energy_density = np.linspace(4 * B_cgs, 10 * B_cgs, 1000) # cgs + # epsilon has a minimum value of 4B so that pressure >= 0 + + pressure = ((energy_density / 3) - (4 * B_cgs / 3)) + + return energy_density, pressure \ No newline at end of file diff --git a/docs/source/MITbag_EOS.rst b/docs/source/MITbag_EOS.rst new file mode 100644 index 0000000..25742b1 --- /dev/null +++ b/docs/source/MITbag_EOS.rst @@ -0,0 +1,9 @@ +.. _MITbag_EOS: + +MIT bag EOS solver +===================== + +Functions to compute MIT bag Equation of state from given parameters. + +.. automodule:: EOSgenerators.MITbag_EOS + :members: diff --git a/docs/source/Polytrope_EOS.py b/docs/source/Polytrope_EOS.py new file mode 100644 index 0000000..cac9e21 --- /dev/null +++ b/docs/source/Polytrope_EOS.py @@ -0,0 +1,88 @@ +import numpy as np +from scipy.optimize import root + +from TOVsolver.unit import g_cm_3, dyn_cm_2 + +def compute_EOS(rhos, theta, + rho_t_start=4.3721E11*g_cm_3, P_t_start=7.7582E29*dyn_cm_2): + """ + Calculate the pressure of a neutron star based on density using a piecewise polytropic equation of state (EOS). + + This function computes the pressure (`pres`) as a function of density (`rho`) by applying different polytropic indices + (`gamma_1`, `gamma_2`, ..., `gamma_n`) within specified density thresholds (`rho_t_start`, `rho_t_1`, `rho_t_2`, ..., `rho_t_n-1`). + The EOS is defined in n distinct regions: + + - **Low-density region:** `rho_t_start < rho <= rho_t_1` + - **Intermediate-density region:** `rho_t_1 < rho <= rho_t_2`, ..., `rho_t_k < rho <= rho_t_k+1`, ..., `rho_t_n-2 < rho <= rho_t_n-1` + - **High-density region:** `rho > rho_t_n-1` + + Parameters + ---------- + rho : array-like + An array of density values (in cgs units) at which to calculate the pressure. + + theta : array-like, length `2 * n - 1` + A list or tuple containing the EOS parameters in the following order: + + - `gamma_1` (float): Polytropic index for the low-density region. + - `gamma_2` to `gamma_2_n-1` (float): Polytropic index for the intermediate-density region. + - `gamma_n` (float): Polytropic index for the high-density region. + - `rho_t_1` (float): Density threshold between the low and intermediate-density regions (in cgs units). + - `rho_t_2` to `rho_t_n-2` (float): Density for the intermediate regions (in cgs units). + - `rho_t_n-1` (float): Density threshold between the intermediate and high-density regions (in cgs units). + + rho_t_start : float + Start point of the density for polytropic EOS + + P_t_start : float + Start point of the pressure for polytropic EOS + Returns + ------- + pres : ndarray + An array of pressure values (in cgs units) corresponding to the input density values. + """ + + n = len(theta) // 2 + 1 + gammas = theta[:n] + rho_ts = theta[n:] + + P_ts, ks = np.zeros((2, n)) + rho_ts = np.insert(rho_ts, 0, rho_t_start) + + # Calculate the values of P_t and k, based on the parameter gamma and rho_t + for i in range(len(ks)): + if i == 0: + P_ts[i] = P_t_start + ks[i] = P_ts[i] / (rho_ts[i]**gammas[i]) + else: + P_ts[i] = ks[i-1] * (rho_ts[i]**gammas[i-1]) + ks[i] = P_ts[i] / (rho_ts[i]**gammas[i]) + + # Construct judgement criteria for each section + conds = [np.array(rhos > rho_ts[i]) & np.array(rhos <= rho_ts[i+1]) for i in range(n-1)] + conds.append(rhos > rho_ts[-1]) + + # Build polytrope functions to compute the pressure for each section + functions = [lambda rho, k=ks[i], g=gammas[i]: k * (rho ** g) for i in range(n)] + + # Calculate pressure by using np.piecewise, based on the conds and functions we defined above. + pres = np.piecewise(rhos, conds, functions) + + return pres + +def fun_gamma_max(rho2, rho1, p1): + """Outputs the maximum gamma for given densities and pressure at both ends of a section of polytropic function. + Args: + rho2 (float): density at the end point. + rho1 (float): density at the start point. + p1 (float): pressure at the start point. + Returns: + gamma_max (float): the maximum gamma. + """ + def fun(x): + a1 = np.log(rho2/p1/x) + a2 = np.log(rho2/rho1) + return a1 / a2 - x + + gamma_max = root(fun, [3/4]).x[0] + return gamma_max diff --git a/docs/source/Polytrope_EOS.rst b/docs/source/Polytrope_EOS.rst new file mode 100644 index 0000000..bd6a9ba --- /dev/null +++ b/docs/source/Polytrope_EOS.rst @@ -0,0 +1,9 @@ +.. _Polytrope_EOS: + +polytrope EOS solver +===================== + +Functions to compute polytrope Equation of state from given parameters. + +.. automodule:: EOSgenerators.Polytrope_EOS + :members: diff --git a/docs/source/RMF_DDH.py b/docs/source/RMF_DDH.py new file mode 100644 index 0000000..1dfb7bc --- /dev/null +++ b/docs/source/RMF_DDH.py @@ -0,0 +1,635 @@ +########################### Imports ############################ +from TOVsolver.unit import g_cm_3, dyn_cm_2, km, Msun, MeV +from scipy import optimize +import pandas as pd +import sympy as sp +import numpy as np +import math + + +c = 3e10 +G = 6.67428e-8 +Msun = 1.989e33 + +dyncm2_to_MeVfm3 = 1.0 / (1.6022e33) +gcm3_to_MeVfm3 = 1.0 / (1.7827e12) +oneoverfm_MeV = 197.33 + +m_e = 2.5896 * 10**-3 +m_mu = 0.53544 +m_n = 4.7583690772 +m_p = 4.7583690772 + +J_B = 1 / 2.0 +b_B = 1 + +m_l = np.array([m_e, m_mu]) +m_b = np.array([m_p, m_n]) + +Matrix_b = np.array( + [[1.0, 1.0, 1 / 2.0, 1.0, 1.0, 1.0], [1.0, 0.0, -1 / 2.0, 1.0, 1.0, 1.0]] +) + +Matrix_l = np.array([[0.0, -1.0, 1 / 2.0], [0.0, -1.0, 1 / 2.0]]) + + + +def Function(type='Typel99', couplings="Default"): + """ + Defines the density-dependent couplings for the sigma, omega, and rho mesons + (g_sigma(rho), g_omega(rho), g_rho(rho)) based on various models. + + The function returns symbolic expressions, allowing users to evaluate the couplings + and their derivatives at any given density rho. Coupling types and default values are + based on literature references, with options for user-defined functions. + + Parameters: + type (str): Specifies the model for density dependence. Options include: + - "Typel99" (default): Uses the DD-MEX model from Typel99 (2008.04491). + - "Malik22": Uses the model from Malik22 (2201.12552). + - "Char23": Uses the model from Char23 (2307.12364). + - "UserDefined": Allows users to specify a custom LaTeX expression for each function. + + couplings (str or list): Sets the model-specific coupling constants. + - If "Default", the default parameters are used: + - "Typel99": DD-MEX model parameters. + - "Malik22": DDBm model parameters. + - "Char23": Model I parameters. + - If not "Default", expect a 1D array with specific coupling values. + + Returns: + tuple: Contains six symbolic expressions to be evaluated at density rho: + - gs (Sympy Expression): Sigma meson coupling. + - gw (Sympy Expression): Omega meson coupling. + - gr (Sympy Expression): Rho meson coupling. + - dgs (Sympy Expression): Derivative of sigma coupling with respect to rho. + - dgw (Sympy Expression): Derivative of omega coupling with respect to rho. + - dgr (Sympy Expression): Derivative of rho coupling with respect to rho. + """ + + if type == 'Typel99': + """ + Ref. 2008.04491v1 + """ + if couplings == "Default": + # DD-MEX model + as_, av, ar, bs, bv, cs, cv, ds, dv, gs0, gv0, gr0, rho0 = [1.3970 , 1.3936 , 0.6202, + 1.3350 , 1.0191 , + 2.0671 , 1.6060 , + 0.4016 , 0.4556 , + 10.7067, 13.3388, 7.2380, 0.153] + else: + as_, av, ar, bs, bv, br, cs, cv, cr, ds, dv, dr, gs0, gv0, gr0, rho0 = couplings + + x = sp.symbols("x") + + gs = sp.lambdify(x, sp.simplify(f"{gs0} * {as_} * (1+{bs}*(x/{rho0} + {ds})**(2))/(1+{cs}*(x/{rho0} + {ds})**(2))")) + gw = sp.lambdify(x, sp.simplify(f"{gv0} * {av} * (1+{bv}*(x/{rho0} + {dv})**(2))/(1+{cv}*(x/{rho0} + {dv})**(2))")) + gr = sp.lambdify(x, sp.simplify(f"{gr0} * exp(-{ar}*(x/{rho0} - 1))")) + + dgs = sp.lambdify(x,sp.simplify(f"2*{gs0}*{as_}*(x/{rho0} + {ds})*({bs}-{cs})/(1+{cs}*(x/{rho0} + {ds})**(2))**(2) / {rho0}")) + dgw = sp.lambdify(x,sp.simplify(f"2*{gv0}*{av} *(x/{rho0} + {dv})*({bv}-{cv})/(1+{cv}*(x/{rho0} + {dv})**(2))**(2) / {rho0}")) + dgr = sp.lambdify(x,sp.simplify(f"- {ar} * ({gr0} * exp(-{ar}*(x/{rho0} - 1))) / {rho0} ")) + + return gs, gw, gr, dgs, dgw, dgr + + elif type == 'Malik22': + """ + https://doi.org/10.3847/1538-4357/ac5d3c + """ + if couplings == "Default": + # DDBm model + as_, av, ar, gs0, gv0, grho0, rho0 = [0.086372, 0.054065, 0.509147, 9.180364, 10.981329, 3.826364*2, 0.150] + else: + as_, av, ar, gs0, gv0, grho0, rho0 = couplings + + x0 = sp.symbols("x0") + + gs = sp.simplify(f"{gs0} * exp(-((x0/{rho0})**{as_} - 1.0))") + gw = sp.simplify(f"{gv0} * exp(-((x0/{rho0})**{av} - 1.0))") + gr = sp.simplify(f"{grho0} * exp(-{ar} * ((x0 / {rho0}) - 1.0))") + + dgs = sp.lambdify(x0,sp.simplify(f"-({gs} * {as_} / {rho0}) * (x0 / {rho0})**({as_} - 1.0)")) + dgw = sp.lambdify(x0,sp.simplify(f"-({gw} * {av} / {rho0}) * (x0 / {rho0})**({av} - 1.0)")) + dgr = sp.lambdify(x0,sp.simplify(f"-({gr} * {ar} / {rho0})")) + + gs = sp.lambdify(x0, gs) + gw = sp.lambdify(x0, gw) + gr = sp.lambdify(x0, gr) + + return gs, gw, gr, dgs, dgw, dgr + + elif type == 'Char23': + """ + Functions used in: PhysRevD.108.103045 e-Print: 2307.12364 [nucl-th] + """ + if couplings == "Default": + # Model I + as_, av, ar, bs, bv, br, cs, cv, cr, ds, dv, dr, rho0 = [8.225494 , 10.426752, 0.64584657, + 2.7079569, 1.6468675, 5.2033131 , + 2.4776689, 6.8349408, 0.4262597 , + 3.8630221, 1.4458185, -0.1824181, 0.16194209] + else: + as_, av, ar, bs, bv, br, cs, cv, cr, ds, dv, dr, rho0 = couplings + + x = sp.symbols("x") + n0 = 0.16 ## In calculation of Char23 the n0 plays as normalization factor and fixed to 0.16 fm-3 + + gs = sp.lambdify(x, sp.simplify(f"{as_} + ({bs} + {ds}*(x/{n0})**(3))*exp(-{cs}*x/{n0})")) + gw = sp.lambdify(x, sp.simplify(f"{av} + ({bv} + {dv}*(x/{n0})**(3))*exp(-{cv}*x/{n0})")) + gr = sp.lambdify(x, sp.simplify(f"({ar} + ({br} + {dr}*(x/{n0})**(3))*exp(-{cr}*x/{n0}))*2")) # Multiplied by 2 ( Lagrange density defintion is different ) + + dgs = sp.lambdify(x,sp.simplify(f"(3*{ds}*(x/{n0})**(2) - {cs}*({bs} + {ds}*(x/{n0})**(3)) )*exp(-{cs}*(x/{n0}))")) + dgw = sp.lambdify(x,sp.simplify(f"(3*{dv}*(x/{n0})**(2) - {cv}*({bv} + {dv}*(x/{n0})**(3)) )*exp(-{cv}*(x/{n0}))")) + dgr = sp.lambdify(x,sp.simplify(f"((3*{dr}*(x/{n0})**(2) - {cr}*({br} + {dr}*(x/{n0})**(3)) )*exp(-{cr}*(x/{n0})))*2")) # Multiplied by 2 + + return gs, gw, gr, dgs, dgw, dgr + + elif type == "UserDefined": + if couplings[-1] == "latex": + from sympy.parsing.latex import parse_latex + couplings = [parse_latex(couplings[i]) for i in range(3)] + + x = sp.symbols("x") + + gs = sp.simplify(couplings[0]) + dgs = sp.lambdify(x,sp.diff(gs, x)) + gs = sp.lambdify(x, gs) + + gw = sp.simplify(couplings[1]) + dgw = sp.lambdify(x,sp.diff(gw, x)) + gw = sp.lambdify(x, gw) + + gr = sp.simplify(couplings[2]) + dgr = sp.lambdify(x,sp.diff(gr, x)) + gr = sp.lambdify(x, gr) + + return gs, gw, gr, dgs, dgw, dgr + + + + +############################### Beta Equilibrium ############################### +def initial_guess(rho, theta): + """ + Computes initial values for the sigma, omega, and rho meson fields, + as well as the chemical potentials for neutrons and electrons at a given density. + + This method provides an initial approximation of these fields and potentials + based on the specified nuclear density and parameters of the chosen RMF model. + + Parameters: + rho (float): Nuclear density at which the initial guess is calculated. + theta (array): Array of parameters specific to the chosen RMF model. + The number and meaning of these parameters vary by model. + + Returns: + tuple: + sigma (float): Initial value of the sigma meson field. + omega (float): Initial value of the omega meson field. + rho_03 (float): Initial value of the rho meson field. + mu_n (float): Initial neutron chemical potential. + mu_e (float): Initial electron chemical potential. + """ + + m_sig, m_w, m_rho, gsf, gwf, grf, dgsf, dgwf, dgrf, rho0 = theta + + g_sigma = gsf(rho) + g_omega = gwf(rho) + g_rho = grf(rho) + + sigma = g_sigma*rho/(m_sig**2) + omega = rho*((m_w**2)/g_omega) + rho_03 = -g_rho*rho/(2.*(m_rho**2)) + + m_eff_n = m_b[1]-(g_sigma*sigma) # Thesis Tiago 2.28 m_eff_Neutron + mu_n = m_eff_n + g_omega*omega + g_rho*rho_03*Matrix_b[1, 2] + mu_e = 0.12*m_l[0]*(rho/rho0)**(2/3.) + + return math.sqrt(sigma), math.sqrt(omega), rho_03, math.sqrt(mu_n), math.sqrt(mu_e) + + + +def beta_equilibrium_function(x, args): + """ + Iteratively adjusts the sigma, omega, and rho meson fields, as well as the chemical potentials + of neutrons and electrons, to achieve beta equilibrium at a specified density. + + Parameters: + x (array): Initial values for sigma, omega, rho meson fields, and chemical potentials + obtained from an initial guess function. + args (array): Model parameters defining the RMF Lagrangian for the chosen RMF model. + + Returns: + tuple: + sigma (float): Adjusted value of the sigma meson field. + omega (float): The adjusted value of the omega meson field. + rho_03 (float): Adjusted value of the rho meson field. + mu_n (float): Adjusted neutron chemical potential. + mu_e (float): Adjusted electron chemical potential. + """ + + m_sig, m_w, m_rho, g_sigma, g_omega, g_rho, dg_sigma, dg_omega, dg_rho, rho0, rho = args + + sigma_sqrt, omega_sqrt, rho_03, mu_n_sqrt, mu_e_sqrt = x + + sigma = sigma_sqrt**2 + omega = omega_sqrt**2 + mu_n = mu_n_sqrt**2 + mu_e = mu_e_sqrt**2 + + rho_B_list = [] + rho_l_list = [] + rho_SB_list = [] + q_list = [] + + m_eff = [m - (g_sigma*sigma) for m in m_b] + + Sigma_0R = dg_omega*omega*rho - dg_sigma*(sigma**2)*(m_sig**2)/g_sigma + dg_rho*(rho_03**2)*(m_rho**2)/g_rho + + for i in range(len(Matrix_b)): + + mu_b = Matrix_b[i,0]*mu_n - Matrix_b[i, 1]*mu_e + + E_fb = mu_b - g_omega*omega - g_rho*rho_03*Matrix_b[i,2] - Sigma_0R + + k_fb_sq = E_fb**2 - m_eff[i]**2 + if k_fb_sq <= 0: + k_fb_sq = 0 + E_fb = m_eff[i] + + k_fb = math.sqrt(k_fb_sq) + + rho_B = k_fb**3 / (3.*math.pi**2) + rho_SB = (m_eff[i]/(2.*math.pi**2))*(E_fb*k_fb - (m_eff[i]**2)*np.log((E_fb + k_fb )/m_eff[i])) + + rho_B_list.append(rho_B) + rho_SB_list.append(rho_SB) + + Q_B = Matrix_b[i,1]*rho_B + q_list.append(Q_B) + + for j in range(len(Matrix_l)): + + mu_l = Matrix_l[j,0]*mu_n - Matrix_l[j,1]*mu_e + + k_fl_sq = mu_l**2 - m_l[j]**2 + if k_fl_sq < 0.0: + k_fl_sq = 0.0 + mu_l = m_l[j] + k_fl = math.sqrt(k_fl_sq) + + rho_l = k_fl**3 / (3.*math.pi**2) + Q_L = Matrix_l[j,1]*rho_l + q_list.append(Q_L) + rho_l_list.append(rho_B) + + f = [(sigma*(m_sig**2)/g_sigma - sum(rho_SB_list)), + ( omega*(m_w**2)/g_omega - sum(rho_B_list) ), + ( rho_03*(m_rho**2)/g_rho - sum([rho_B_list[k]*Matrix_b[k][2] for k in range(len(Matrix_b))])), + ( rho - sum(rho_B_list)), + ( sum(q_list) )] + + return f + +def get_energy_pressure(x, rho, theta): + """ + Computes the energy density and pressure at a specified density, using given parameters + for an RMF model. + + Parameters: + x (array): Array containing initial values of the sigma, omega, and rho meson fields, + as well as the chemical potentials, obtained from an initial guess function. + rho (float): Nuclear density at which to compute energy density and pressure. + theta (array): Array of parameters defining the RMF model in the Lagrangian. + + Returns: + tuple: + energy_density (float): The energy density in natural units, essential for the EOS. + pressure (float): The pressure in natural units, also essential for the EOS. + """ + + sigma_sqrt, omega_sqrt, rho_03, mu_n_sqrt, mu_e_sqrt = x + sigma = sigma_sqrt**2 + omega = omega_sqrt**2 + mu_n = mu_n_sqrt**2 + mu_e = mu_e_sqrt**2 + + m_sig, m_w, m_rho, g_sigma, g_omega, g_rho, dg_sigma, dg_omega, dg_rho, rho0 = theta + + energy_b = 0 + energy_l = 0 + multi = 0 + m_eff = [m - (g_sigma*sigma) for m in m_b] + + rho_S = sigma * m_sig**2 / g_sigma + Sigma_0R = dg_omega*omega*rho - dg_sigma*sigma*rho_S + dg_rho*(rho_03**2)*(m_rho**2)/g_rho + Pressure_bl = 0 + q_list = [] + for i in range(len(Matrix_b)): + + mu_b = Matrix_b[i,0]*mu_n - Matrix_b[i, 1]*mu_e + + E_fb = mu_b - g_omega*omega - g_rho*rho_03*Matrix_b[i,2] - Sigma_0R + + k_fb_sq = E_fb**2 - m_eff[i]**2 + + if k_fb_sq <= 0: + k_fb_sq = 0.0 + E_fb = m_eff[i] + + k_fb = math.sqrt(k_fb_sq) + + rho_B = (k_fb**3) / (3.*math.pi**2) + + energy_baryon = (1/(8.*(math.pi**2)))*(k_fb*E_fb*(2*k_fb**2+m_eff[i]**2) - np.log((k_fb + E_fb)/m_eff[i])*m_eff[i]**4) + + energy_b += energy_baryon + + integral_b = 1/4 * ( 1.5 * m_eff[i]**4*np.arctanh(k_fb/E_fb) - 1.5*k_fb*m_eff[i]**2*E_fb + k_fb**3*E_fb ) + Pressure_bl += 1/3*(1/math.pi**2 * integral_b) + + if i == 0: + alpha = rho_B/rho + + for j in range(len(Matrix_l)): + + mu_l = Matrix_l[j, 0]*mu_n - Matrix_l[j, 1]*mu_e + + k_fl_sq = mu_l**2 - m_l[j]**2 + if k_fl_sq < 0.0: + k_fl_sq = 0.0 + mu_l = m_l[j] + k_fl = math.sqrt(k_fl_sq) + + energy_lepton = (1/(8.*(math.pi**2)))*(k_fl*mu_l*(2*k_fl**2+m_l[j]**2)-(m_l[j]**4)*np.log((k_fl+mu_l)/m_l[j])) + energy_l += energy_lepton + integral_l = 1/4 * ( 1.5 * m_l[j]**4*np.arctanh(k_fl/mu_l) - 1.5*k_fl*m_l[j]**2*mu_l + k_fl**3*mu_l ) + Pressure_bl += 1/3*(1/math.pi**2) * integral_l + + sigma_terms = 0.5*((sigma*m_sig)**2) + omega_terms = 0.5*((omega*m_w)**2) + rho_terms = 0.5*((rho_03*m_rho)**2) + + energy_density = energy_b + energy_l + sigma_terms + omega_terms + rho_terms + + Pressure = Pressure_bl - sigma_terms + omega_terms + rho_terms + Sigma_0R*rho + return energy_density, Pressure, alpha + +def compute_eos(eps_crust, pres_crust, theta): + """ + Computes the core part of the equation of state (EOS) using a chosen RMF model, + complementing it with crust EOS data. + + Parameters: + eps_crust (array): Energy density values for the crust EOS + pres_crust (array): Pressure values for the crust EOS + theta (array): Array of parameters defining the RMF model in the Lagrangian + for the core EOS. + + Returns: + tuple: + energy_density (float): Energy density of the core in natural units, an essential EOS ingredient. + pressure (float): Pressure of the core in natural units, is also essential for the EOS. + """ + + rho0 = theta[9] + x_init = np.array(initial_guess(0.04, theta)) + dt = 0.006 + N_points = 125 + + Density = np.empty(N_points, dtype=float) + Energy = np.empty(N_points, dtype=float) + Pressure = np.empty(N_points, dtype=float) + Alpha = np.empty(N_points, dtype=float) + for i in range(N_points): + rho = 0.04 + i*dt + + theta_in = theta.copy() + theta_in[3] = theta_in[3](rho) + theta_in[4] = theta_in[4](rho) + theta_in[5] = theta_in[5](rho) + theta_in[6] = theta_in[6](rho) + theta_in[7] = theta_in[7](rho) + theta_in[8] = theta_in[8](rho) + theta_in = np.array([float(t) for t in theta_in], dtype=np.float64) + + arg = np.append(theta_in, float(rho)) + sol = optimize.root(beta_equilibrium_function, x_init ,method='lm', args = arg) + if not sol.success: + raise ValueError("Did not converge") + x_init = sol.x + Re = get_energy_pressure(x_init, rho, theta_in) + + Energy[i] = Re[0] + Pressure[i] = Re[1] + Alpha[i] = Re[2] + Density[i] = rho + + end = 0 + for i in range(0, len(Energy) - 1): + if Energy[i] > max(eps_crust) and i > 18: + end = i + 2 + break + end += 1 + + Density = Density[end::] + Energy = Energy[end::] + Pressure = Pressure[end::] + Alpha = Alpha[end::] + + return Density, Energy, Pressure, Alpha + + + + + +###################################### Alpha Depedent #######################3########## +def initial_guess_alpha(rho, theta): + """ + Provides initial estimates for the sigma, omega, and rho meson fields at a specified + nuclear density, based on the parameters of a selected RMF model. + + Parameters: + rho (float): Nuclear density at which to calculate initial field values. + theta (array): Array of parameters defining the RMF model in the Lagrangian. + + Returns: + tuple: + sigma (float): Initial value of the sigma meson field. + omega (float): Initial value of the omega meson field. + rho_03 (float): Initial value of the rho meson field. + """ + m_sig, m_w, m_rho, gsf, gwf, grf, dgsf, dgwf, dgrf, rho0, alpha = theta + g_sigma = gsf(rho) + g_omega = gwf(rho) + g_rho = grf(rho) + + sigma = g_sigma*rho/(m_sig**2) + omega = rho*((m_w**2)/g_omega) + rho_03 = -g_rho*rho/(2.*(m_rho**2)) + + return math.sqrt(sigma), math.sqrt(omega), rho_03 + +def fields_with_alpha(x, args): + """ + Iteratively adjusts the sigma, omega, and rho meson fields at a given nuclear density, + based on a specified proton fraction, using the parameters of a selected RMF model. + + Parameters: + x (array): Initial values of sigma, omega, rho meson fields, and chemical potentials + from an initial guess function. + args (array): Parameters defining the RMF model in the Lagrangian. The last element + in this array specifies the proton fraction (alpha): + - alpha = 0 for pure neutron matter + - alpha = 0.5 for symmetric nuclear matter + - alpha = rho_p / (rho_p + rho_n) + + Returns: + tuple: + sigma (float): Adjusted value of the sigma meson field. + omega (float): Adjusted value of the omega meson field. + rho_03 (float): Adjusted value of the rho meson field. + """ + + m_sig, m_w, m_rho, g_sigma, g_omega, g_rho, dg_sigma, dg_omega, dg_rho, rho0, alpha, rho = args + + sigma_sqrt, omega_sqrt, rho_03 = x + + sigma = sigma_sqrt**2 + omega = omega_sqrt**2 + + m_eff = m_b[1] - (g_sigma*sigma) + + # Proton + rho_p = alpha*rho + kf_p = (rho_p*(3*math.pi**2))**(1/3) + E_fp = (kf_p**2 + m_eff**2)**(1/2) + rho_SB_p = (m_eff/(2.*math.pi**2))*(E_fp*kf_p - (m_eff**(2))*np.arctanh(kf_p/E_fp)) + + #Neutron + rho_n = (1-alpha)*rho + kf_n = (rho_n*(3*math.pi**2))**(1/3) + E_fn = (kf_n**2 + m_eff**2)**(1/2) + rho_SB_n = (m_eff/(2.*math.pi**2))*(E_fn*kf_n - (m_eff**(2))*np.arctanh(kf_n/E_fn)) + + rho_b = rho_p + rho_n + rho_SB = rho_SB_p + rho_SB_n + + f =[( rho - rho_b ), + ( sigma *(m_sig**2)/g_sigma - rho_SB ), + ( omega *(m_w**2) /g_omega - rho ), + ( rho_03*(m_rho**2)/g_rho - (alpha-0.5)*rho )] + return f + +def get_energy_pressure_alpha(x, rho, theta): + """ + Computes the energy density and pressure for a specified nuclear density and proton fraction, + using parameters from a chosen RMF model. + + Parameters: + x (array): Array containing initial values of sigma, omega, rho meson fields, and + chemical potentials from an initial guess function. + rho (float): Central density at which the equation of state (EOS) calculation begins. + theta (array): Parameters defining the RMF model in the Lagrangian. The last element in + this array specifies the proton fraction (alpha): + - alpha = 0 for pure neutron matter + - alpha = 0.5 for symmetric nuclear matter + - alpha = rho_p / (rho_p + rho_n) + + Returns: + tuple: + energy_density (float): Energy density in natural units, essential for the EOS. + pressure (float): Pressure in natural units, is essential for the EOS. + """ + + sigma_sqrt, omega_sqrt, rho_03 = x + + sigma = sigma_sqrt**2 + omega = omega_sqrt**2 + + m_sig, m_w, m_rho, g_sigma, g_omega, g_rho, dg_sigma, dg_omega, dg_rho, rho0, alpha = theta + + m_eff = m_b[1] - (g_sigma*sigma) + + rho_S = sigma * m_sig**2 / g_sigma + Sigma_0R = dg_omega*omega*rho - dg_sigma*sigma*rho_S + dg_rho*(rho_03**2)*(m_rho**2)/g_rho + + #Proton + rho_p = alpha*rho + kf_p = (rho_p*(3*math.pi**2))**(1/3) + E_fp = (kf_p**2 + m_eff**2)**(1/2) + energy_p = (1/(8.*(math.pi**2)))*(kf_p*E_fp*(2*kf_p**2+m_eff**2) - np.log((kf_p + E_fp)/m_eff)*m_eff**4) + integral_p = 1/4 * ( 1.5 * m_eff**4*np.arctanh(kf_p/E_fp) - 1.5*kf_p*m_eff**2*E_fp + kf_p**3*E_fp ) + Pressure_p = 1/3*(1/math.pi**2 * integral_p) + Sigma_0R*rho_p + + #Neutron + rho_n = (1-alpha)*rho + kf_n = (rho_n*(3*math.pi**2))**(1/3) + E_fn = (kf_n**2 + m_eff**2)**(1/2) + energy_n = (1/(8.*(math.pi**2)))*(kf_n*E_fn*(2*kf_n**2+m_eff**2) - np.log((kf_n + E_fn)/m_eff)*m_eff**4) + integral_n = 1/4 * ( 1.5 * m_eff**4*np.arctanh(kf_n/E_fn) - 1.5*kf_n*m_eff**2*E_fn + kf_n**3*E_fn ) + Pressure_n = 1/3*(1/math.pi**2 * integral_n) + Sigma_0R*rho_n + + #Total + energy_b = energy_p + energy_n + + sigma_terms = 0.5*((sigma*m_sig)**2) + omega_terms = 0.5*((omega*m_w)**2) + rho_terms = 0.5*((rho_03*m_rho)**2) + + energy_density = energy_b + sigma_terms + omega_terms + rho_terms + Pressure = Pressure_p + Pressure_n - sigma_terms + omega_terms + rho_terms + + return energy_density, Pressure + +def compute_eos_alpha(theta): + """ + Generates the equation of state (EOS) table for a sequence of densities, considering a specified + proton fraction, using parameters from the chosen RMF model. + + Parameters: + eps_crust (array): Energy density values for the crust EOS, in MeV/fm³ (including a G/c² factor). + pres_crust (array): Pressure values for the crust EOS, in MeV/fm³ (including a G/c⁴ factor). + theta (array): Parameters defining the RMF model in the Lagrangian. The last element specifies + the proton fraction (alpha): + - alpha = 0 for pure neutron matter + - alpha = 0.5 for symmetric nuclear matter + - alpha = rho_p / (rho_p + rho_n) + + Returns: + tuple: + energy_density (float): Energy density in natural units, an essential EOS ingredient. + pressure (float): Pressure in natural units, also essential for the EOS. + """ + rho0 = theta[9] + x_init = np.array(initial_guess_alpha(0.04, theta)) + dt = 0.006 + N_points = 125 + + Density = np.empty(N_points, dtype=float) + Energy = np.empty(N_points, dtype=float) + Pressure = np.empty(N_points, dtype=float) + + for i in range(N_points): + rho = 0.04 + i*dt + + theta_in = theta.copy() + theta_in[3] = theta_in[3](rho) + theta_in[4] = theta_in[4](rho) + theta_in[5] = theta_in[5](rho) + theta_in[6] = theta_in[6](rho) + theta_in[7] = theta_in[7](rho) + theta_in[8] = theta_in[8](rho) + theta_in = np.array([float(t) for t in theta_in], dtype=np.float64) + + arg = np.append(theta_in, float(rho)) + sol = optimize.root(fields_with_alpha, x_init ,method='lm', args = arg) + x_init = sol.x + Re = get_energy_pressure_alpha(x_init, rho, theta_in) + + Density[i] = rho + Energy[i] = Re[0] + Pressure[i] = Re[1] + + return Density, Energy, Pressure diff --git a/docs/source/RMF_DDH.rst b/docs/source/RMF_DDH.rst new file mode 100644 index 0000000..95f80c9 --- /dev/null +++ b/docs/source/RMF_DDH.rst @@ -0,0 +1,9 @@ +.. _RMF_DDH_EOS: + +DDH RMF EOS solver +===================== + +Functions to compute DDH RMF Equation of state from given parameters. + +.. automodule:: EOSgenerators.RMF_DDH + :members: diff --git a/docs/source/SpeedofSound_EOS.py b/docs/source/SpeedofSound_EOS.py new file mode 100644 index 0000000..2ced456 --- /dev/null +++ b/docs/source/SpeedofSound_EOS.py @@ -0,0 +1,106 @@ +import numpy as np +import math +from scipy import optimize +from TOVsolver.unit import g_cm_3, dyn_cm_2, e0 +from scipy.integrate import cumulative_simpson + + +##################### SPEED-OF-SOUND-EOS BEGIN ############################## + + +class compute_EOS: + def __init__(self, x_last, y_last, dydx_last, enablePTcheck=False) -> None: + """ + Speed of sound EOS class. + See https://arxiv.org/abs/1812.08188 for details. + Since the speed of sound core eos need to fit the last point and derivative of outer eos, + the initial values are provided here. + + Args: + x_last (float): Last energy density value. + y_last (float): Last pressure value. + dydx_last (float): Last derivative value. + enablePTcheck (bool, optional): Enable phase transition check. Defaults to False. + """ + self.x_last = x_last + self.y_last = y_last + self.dydx_last = dydx_last + self.enablePTcheck = enablePTcheck + + def cs2(self, x, a): + a1, a2, a3, a4, a5, a6 = a + ret = ( + a1 * np.exp(-0.5 * (((x - a2) / a3) ** 2)) + + a6 + + (1 / 3 - a6) / (1 + np.exp(-a5 * (x - a4))) + ) + return np.clip(ret, 0, 1) # requirement 2 and 4 + + def cal_a6(self, a1, a2, a3, a4, a5): + A = a1 * np.exp(-0.5 * (((self.x_last - a2) / a3) ** 2)) + (1 / 3) / ( + 1 + np.exp(-a5 * (self.x_last - a4)) + ) + B = 1 - 1 / (1 + np.exp(-a5 * (self.x_last - a4))) + # cs2 = A + B * a6 = dydx + return (self.dydx_last - A) / B + + def uniform_prior(self, a, b, r): + return a + (b - a) * r + + def gen_a(self, cube5_): + """ + Generate the EOS parameters from the given cube. + + Args: + cube5_ (array): The input cube. The shape is (5,). The values are in the range of [0, 1]. + + Returns: + tuple: The EOS parameters. + """ + a1 = self.uniform_prior(0.1, 1.5, cube5_[0]) + a2 = self.uniform_prior(1.5, 12, cube5_[1]) + a3 = self.uniform_prior(0.05 * a2, 2 * a2, cube5_[2]) + a4 = self.uniform_prior(1.5, 37, cube5_[3]) + a5 = self.uniform_prior(0.1, 1, cube5_[4]) + a6 = self.cal_a6(a1, a2, a3, a4, a5) + return (a1, a2, a3, a4, a5, a6) + + def check_a(self, a): + a1, a2, a3, a4, a5, a6 = a + + # PT requirement + if self.enablePTcheck and a6 > 0: + return False + + # requirement 5 + if np.any(self.cs2(np.linspace(0.5 * e0, 1.5 * e0, 16), a) > 0.163): + return False + + # requirement 3 + if a6 > 1 / 3: + return False + if np.any(self.cs2(np.linspace(49 * e0, 50 * e0, 16), a) > 1 / 3): + return False + + return True + + def cal_core_p(self, core_e, a): + """ + Calculate the core pressure. + + Args: + core_e (array): The core energy density. + a (array): The EOS parameters. + + Returns: + array: The core pressure. + """ + core_p = cumulative_simpson(self.cs2(core_e, a), x=core_e, initial=self.y_last) + return core_p + # fix check dpde < 0 issue + # dif_p = np.diff(core_p, prepend=self.y_last) + # dif_p = np.maximum(dif_p, 0) + # return np.cumsum(dif_p) + + +##################### SPEED-OF-SOUND-EOS END ########################### ### diff --git a/docs/source/SpeedofSound_EOS.rst b/docs/source/SpeedofSound_EOS.rst new file mode 100644 index 0000000..488b574 --- /dev/null +++ b/docs/source/SpeedofSound_EOS.rst @@ -0,0 +1,9 @@ +.. _SpeedofSound_EOS: + +Speed of Sound EOS solver +===================== + +Functions to compute Speed of Sound Equation of state from given parameters. + +.. automodule:: EOSgenerators.SpeedofSound_EOS + :members: diff --git a/docs/source/Strangeon_EOS.py b/docs/source/Strangeon_EOS.py new file mode 100644 index 0000000..9922c31 --- /dev/null +++ b/docs/source/Strangeon_EOS.py @@ -0,0 +1,36 @@ +import numpy as np +import math +from scipy import optimize +from TOVsolver.unit import g_cm_3, dyn_cm_2 + +def compute_EOS(n, theta): + """ + Compute the energy density and pressure based on the given parameters. + + Args: + n (array): An array of n values. Input values of baryon number densities. + theta (array): An array representing the parameters [epsilon, Nq, ns]. + epsilon: the depth of the potential well; MeV; + Nq: the number of quarks in a strangeon; + ns: the number density of baryons at the surface of the star; fm^-3 + + Returns: + tuple: Arrays of energy densities in units of gcm3 and pressures in units of dyncm2. + """ + + Nq, epsilon, ns = theta + + A12 = 6.2 + A6 = 8.4 + mq = 300 + """ + mq: the mass of the quark in this EOS. + A12 and A6 are fixed throughout the calculation. + """ + + sigma = np.sqrt(A6 / (2 * A12)) * (Nq / (3 * ns)) + + energy_density = 2 * epsilon * (A12 * sigma**4 * n**5 - A6 * sigma**2 * n**3) + n * Nq * mq + pressure = 4 * epsilon * (2 * A12 * sigma**4 * n**5 - A6 * sigma**2 * n**3) + + return energy_density , pressure \ No newline at end of file diff --git a/docs/source/Strangeon_EOS.rst b/docs/source/Strangeon_EOS.rst new file mode 100644 index 0000000..cd94842 --- /dev/null +++ b/docs/source/Strangeon_EOS.rst @@ -0,0 +1,9 @@ +.. _Strangeon_EOS: + +Strangeon EOS solver +===================== + +Functions to compute Strangeon Equation of state from given parameters. + +.. automodule:: EOSgenerators.Strangeon_EOS + :members: