rhoades_ruffini_1974

def main():

    import sympy
    from scipy.integrate import odeint
    import numpy
    import pylab

    r = sympy.Symbol('r')
    p = sympy.Symbol('p')
    G = sympy.Symbol('G')
    rho = sympy.Symbol('rho')
    c = sympy.Symbol('c')
    m = sympy.Symbol('m')

    tov_equation = (p(r)).diff(r) + (G/r**2)*(rho(r)+p(r)/c**2)*(m(r)+4*sympy.pi*r**3*p(r)/c)/(1-2*G*m(r)/(c**2*r))
    tov_equation = tov_equation.subs(p(r),c**2*rho(r))
    tov_equation = tov_equation.subs(rho(r),(m(r)).diff(r)/(4*sympy.pi*r**2))
    tov_equation = tov_equation.subs(G,1).subs(c,1)

    m_rr = sympy.solve(tov_equation.doit(),m(r).diff(r,2))[0]

    def deriv(y,x):
        return [y[1], m_rr.subs(m(r).diff(r),y[1]).subs(m(r),y[0]).subs(r,x).n()]

    x_range = numpy.logspace(-6,0,1000)
    for slope in numpy.logspace(-3,3,7):
        nsol = odeint(deriv, [0,slope], x_range)
        pylab.loglog(x_range, nsol.T[0], label=str(slope))
    pylab.legend(loc='lower right')
    pylab.xlabel('r')
    pylab.ylabel('m')
    pylab.show()

if __name__ == '__main__':

    main()