import pylab
from scipy.integrate import ode


Gamma0 = 
K0 = 
rho1 = 
P1 = 
K1 =

R_km = 
Msun = 

def eos(rho):
    if rho < rho1: return K0*rho**Gamma0
    else: return K1*rho**Gamma1

def inveos(P):
    if P < P1: return (P/K0)**(1.0/Gamma0)
    else: return (P/K1)**(1.0/Gamma1)

def tov(r, y, arg1):
    P, m = y[0], y[1]
    rho = inveos(P)
    dPdr = -G*(rho + P/c**2)*(m + 4.0*pi*r**3*P/c**2)
    dPdr = dPdr/(r*(r - 2.0*G*m/c**2))
    dmdr = 4.0*pi*r**2*rho
    return [dPdr, dmdr]

def tovsolve(rhoc):
    r = pylab.arange (10.0, 20000.0, 10)
    m = pylab.zeros_like(r)
    P = pylab.zeros_like(r)
    m[0] = 4.0*pi*r[0]**3*rhoc
    P[0] = eos(rhoc)
    sol = ode(tov).set_integrator('vode', method='adams', rtol=1.0e-10, 
        atol=1.0e-10)
    sol.set_initial_value([P[0], m[0]],r[0])
    sol.set_f_params(1.0) 
    i = 0
    while P[i] > 0.0 and i < len(r) - 1 and sol.successful():
        sol.integrate(r[i+1])
        P[i+1] = sol.y[0]
        m[i+1] = sol.y[1]
        i = i + 1
    return m[i-1]/Msun, r[i-1]/R_km