Lazzati 3

import pylab

#number of hydrogen and carbon atoms
nH = 1
nC = 10**(-3.4)*nH

#mass of hydrogen and carbon
mH = 1.67e-24
mC = 12*mH

r0 = int(1e16)
p0 = nH*mH + nC*mC
T0 = 1e8


r = 10**pylab.arange(15,20,0.01)
P = p0*(r/r0)**(-3)
T = pylab.zeros(r.size);
T[0] = T0
for i in range(1,r.size):
    T[i] = T[i-1]*(1 + (2./3)*((P[i]-P[i-1])/(P[i-1])))

pylab.figure(1)
pylab.suptitle('Temperature versus Radius', fontsize=12)
pylab.xlabel('Radius (cm)')
pylab.ylabel('Temperature (K)')
pylab.loglog(r,T)
pylab.show()