#
# graph-diffusion.py
#
# graphs Pb diffusion in zircon with gnuplot over a given period of time
#
# based on a paper by James K.W. Lee, Nature 1997
#
#
# Eric Thern eric NOSPAM(at) thern.org
#
#
import math
import cmath
import sys
# get T and time from command line arguments
#
#T = input("Input Temperature: ")
#time = input("Input Time (in millions of years): ")
T = 0
T2 = 1200
time = 0
time2 = 2000
steps = 25
temp = T
temp2 = T2
# convert T to Kelvin
T = T + 273.16
# E - activation energy ( Kcal Mol^-1)
E = 161
# D - diffusion coefficient (initialized to zero)
D = 0
# D0 (D subzero) - pre-exponential coefficient ( cm2 s-1 )
Dsz = 3.9e+9
# e - exponent (exp)
# math.e
# R - gas constant
# R = 1.98722 cal mol-1 K-1
# R = 0.00198722 Kcal mol-1 K-1 (using this)
R = 0.00198722
# T - temperature (In Kelvin)
temptime = time
temptemp = temp
T_backup = T
while time <= time2:
temp = temptemp
T = T_backup
print ""
print ""
print "#",time
for temp in range(temp, temp2, 25):
D = Dsz*(cmath.e)**(-E/(R*T))
Dum = D * 10e+8 * 3.15569259747e+13
diffusion = Dum * time
print diffusion,temp
T = T + 25
time = time + steps