Lazzati 4

import pylab
import math

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

#constants
r0 = int(1e16)
p0 = 1.7e-14
T0 = 1e8
k =  1.38e-16
r = 10**pylab.arange(15,20,0.01)
P = (p0*(r/r0)**(-3))
nH = P/(mH*(1 + (10**(-3.4))*mC))
nC = 10**(-3.4)*nH

#define temperature
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])))

#define equilbrium density
Ns = pylab.zeros(r.size);
for i in range(1,r.size):
        Ns[i] = ((6.9e13)*math.exp(-844282/T[i]))/(k*T[i])

#define saturation
S = nC/Ns

#plotting
pylab.figure(1)
pylab.title('Saturation versus Radius', fontsize=12)
pylab.xlabel('Radius (cm)')
pylab.ylabel('Saturation (cm^-3)')
pylab.loglog(r,S)
pylab.xlim((1e16,1e17))
pylab.ylim((0,10))
pylab.show()