from ClimateUtilities import *
import phys,math,planets

#Computes mass-radius relation given an equation
#of state. 


#Define the equation of state here.
gas = phys.H2
def rho(p,T):
    return p/(gas.R*T)

#Integrand function for
#right hand side of hydrostatic relation
def derivs(r,Y):
    p,M = Y
    dpdr = -rho(p,T)*phys.G*M/(r*r)
    dMdr = rho(p,T)*4.*math.pi*r*r
    return numpy.array([dpdr,dMdr])

T = 1000.
MEarth = 5.972e24
rEarth = planets.Earth.a

def M(pstart):
    p = pstart
    m = integrator(derivs,1000.,numpy.array([pstart,0.]),10000.)
    while p>1.e7:
        res = m.next()
        p = res[1][0]

    return res[0]/planets.Earth.a,res[1][1]/MEarth

rList = []
MList = []
pList = numpy.array([10.e9 + i*(100.e10-10.e9)/100. for i in range(101)])
for p in pList:
    rr,mm = M(p)
    print p/1.e9, rr,mm
    rList.append(rr)
    MList.append(mm)
c = Curve()
c.Xlabel = 'radius/Re'
c.Ylabel = 'mass/Me'
c.addCurve(rList,'r/re')
c.addCurve(MList,'m/me')
c.dump('MassRad%f.txt'%T)

cp = Curve()
cp.addCurve(pList/1.e9,'p (GPa)')
cp.addCurve(rList,'r/re')
cp.addCurve(MList,'m/me')
cp.Xlabel = 'Pressure (GPa)'