Lazzatti 2

import pylab

r0 = int(10e15)
p0 = 1.7e-14
r = r0
T0 = int(10e7)
rfinal = int(r*10000)
def P(r):
    return p0*pow(r/r0, -3)

r = 10**pylab.arange(5,10,001.)

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

pylab.figure(1)
pylab.loglog(r,T)
pylab.show()