"""Barrier and lifetime arithmetic for the section on one atmosphere: per-site, one-in-N and diffusion-limited requirements, relaxation times, effective barriers read from decomposition temperatures, the AlH3 extrapolation, and the molecular-dynamics exclusion bound.

Run with the project environment: .venv/bin/python numerics/derived/kinetics.py
The printed output is stored beside this file as kinetics.out.txt.
"""
import pathlib
import math
kB=8.617333e-5
yr=3.15576e7; wk=7*86400; day=86400; hr=3600
def E(T,t,nu=1e13): return kB*T*math.log(nu*t)
print('per-site: yr@300', E(300,yr), 'nu1e10', E(300,yr,1e10), 'wk@200', E(200,wk), 'wk@200 nu1e10', E(200,wk,1e10))
print('kT ln10 @300', kB*300*math.log(10), '@200', kB*200*math.log(10), '@270', kB*270*math.log(10))
for N in (1e6,1e9,1e12):
    print('one in N', N, E(300,yr)+kB*300*math.log(N))
# N in a grain
rho=0.085e30 # per m3
for N in (1e12,):
    V=N/rho; print('cube side um', V**(1/3)*1e6, 'sphere diam um', (6*V/math.pi)**(1/3)*1e6)
# diffusion-limited: plate thickness L, tau = L^2/(pi^2 D), D = D0 exp(-E/kT), D0=1e-3 cm2/s
D0=1e-3*1e-4 # m2/s
def Ediff(L,t,T=300,geom=math.pi**2): 
    D=L*L/(geom*t); return kB*T*math.log(D0/D)
for L in (10e-6,1e-3,1e-2):
    print('diff L',L, 'plate', Ediff(L,yr), 'L2/D', Ediff(L,yr,geom=1), 'L2/6D', Ediff(L,yr,geom=6))
def tau(L,Eb,T=300,geom=math.pi**2):
    D=D0*math.exp(-Eb/(kB*T)); return L*L/(geom*D)
for Eb in (0.23,0.46):
    for L in (10e-6,1e-3,1e-2):
        t=tau(L,Eb); print('tau Eb',Eb,'L',L, t,'s', t/hr,'h', t/yr,'yr')
# effective barriers
print('ZrH3', E(200,60), E(270,600))
print('Y3Fe4H20', E(300,30*hr), E(300,30*day), E(300,90*day))
print('AlH3', E(333.15,day), E(373.15,hr), E(413.15,60))
print('anneal 200K 1min..1wk', E(200,60), E(200,wk))
# ratio to 300 K requirement for prefactors
for nu in (1e6,1e8,1e10,1e13):
    print('nu',nu, E(200,60,nu)/E(300,yr,nu), E(200,wk,nu)/E(300,yr,nu), E(200,60,nu),E(200,wk,nu),E(300,yr,nu))
# AlH3 measured: Ea 1.06 eV; one-third decomposed at ~1.5e4..1.8e4 s at 99 C
for Ea in (0.97,1.06,1.10):
  for t13 in (1.5e4,1.8e4):
    k=math.log(1.5)/t13; nu=k*math.exp(Ea/(kB*372.15)); 
    k300=nu*math.exp(-Ea/(kB*300)); print('AlH3 Ea',Ea,'t13',t13,'nu',f'{nu:.2e}','t(1/3) at 300K yr', math.log(1.5)/k300/yr)
# per-site model with 1.06 eV and 1e13: lifetime at 99C and 300K
for T in (372.15,300):
    print('1.06eV,1e13 at',T, math.exp(1.06/(kB*T))/1e13,'s')
# MD exclusion
def Eex(T,t,N,nu=1e13): return kB*T*math.log(nu*t*N)
print('MD 0.5ns 700K 750H', Eex(700,0.5e-9,750), 'lifetime at 300', math.exp(Eex(700,0.5e-9,750)/(kB*300))/1e13)
# PV
print('100 GPa * 1 A3 =', 100*6.2415,'meV', 100*6.2415/67)