I've been working on a code to solve the Tolman-Oppenheimer-Volkoff (TOV) equations for a while and recently I've got it right but only for one specific equation of state, the bag model, which is not what I need since it should run for any hadronic equation of state. The code has given only wrong results for any EoS other than the bag model and I've tried everything I could think of but nothing has made it so that the final diagrams for other EoS are correct. I'll put the code and some of the wrong results for some of the EoS I've been studying below.

import math
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
from scipy.interpolate import UnivariateSpline
from scipy.interpolate import interp1d

pi = math.pi
km2_gevfm3 = (1/(1.3234*10**(-6)))*0.001
gevfm3_km = 1/(km2_gevfm3)
B = 7.424*10**(-5) 
cm = 1/1.4760

eos = np.genfromtxt("C:\\Users\\felip\\Desktop\\TOV\\Equações de estado\\PtableTM1.dat") 
eden, p = (gevfm3_km*eos[:,0])/1000, (gevfm3_km*eos[:,1])/1000 
press = interp1d(eden, p, kind = 'cubic', bounds_error = False, fill_value = 'extrapolate') 
# press = UnivariateSpline(eden, p, k = 5, s = 5)  
dpde = press._spline.derivative(nu = 1) 
# dpde = press.derivative(n = 1)  

def TOV(y, r):
     eden, M = float(y[0]), float(y[1])
     dedr = -(1*(dpde(eden))**(-1))*((M*eden)/(r**2))*(1+press(eden)/eden)*(1+(4*pi*press(eden)*r**(3))/(M))*((1-((2*M)/r))**(-1)) 
     dMdr = 4*pi*eden*(r**2) 
     return [dedr, dMdr] 

r0 = 10**(-10)
rm = 5000  #km
N = 1000
rn = np.linspace(r0, rm, N) 

emin, emax =   7.7*10e-15, 0.0003
h_e = 1000
e_range = np.linspace(emin, emax, h_e) 
R,M,E =[],[],[]

for e0 in e_range:
    M0 = (4*pi*e0*r0**(3))/3 
    contorno = [e0, M0]
    soltov = odeint(TOV, contorno, rn, atol = 10e-10, rtol = 10e-14, mxstep=5000)
    rmax = r0 
    mmax = M0
    esol = soltov[:,0]
    psol = press(esol)  
    for i in psol[:-1]:
        h=rn[k+1]-rn[k] + r0
        rmax = rmax + h
        mmax = soltov[k,1] 
        if i <= 0:

immax = np.argmax(M)
Mmax = M[immax]
Rmax = R[immax]

f1 = plt.figure()
ax1 = f1.add_subplot(111)
ax1.title.set_text('M x E')
ax1.set_xlabel('$(\epsilon)_c [GeV/fm^3]$')
ax1.set_ylabel('$M [M_\odot]$')    

f2 = plt.figure()
ax2 = f2.add_subplot(111)
ax2.title.set_text('M x R')
ax2.set_xlabel('$R [Km]$')
ax2.set_ylabel('$M [M_\odot]$')    

I couldn't get the ident right but it is correct on my code. I have no idea what could be changed here and I've gone through every line in this code and can't figure out why it shouldn't work for other EoS.

For the BPS EoS, the results are : enter image description hereenter image description here.

  • 4
    $\begingroup$ You may be more likely to get answers if you write down the equations you're solving and the discretization you're using. $\endgroup$ Jun 25 at 4:57


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.