#
# 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