idl_scripts/lamberts_tunnel.pro at master · cnshingledecker/idl_scripts · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
; The barrier function from Lamberts and Kaestner 2017
; returns the Langmuir-Hinshelwood rate coefficient (s-1)
function lamberts_tunnel, T
  fac1 = 1.51e10*((T/300.0)^0.86)
  fac2 = exp(-1750.0*(T+180.0)/((T^2)+(180^2)))
  return, fac1*fac2
end

; The rate of hopping over the barrier
function thermal_barrier, EA, T
  return, 1.0e12*exp(-EA/T)
end

function kappa, EACT, EA, EB, T
  k1 = lamberts_tunnel(T)
  k2 = thermal_barrier(EACT,T)
  k3 = thermal_barrier(EA,T)
  k4 = thermal_barrier(EB,T)
  return, (k1+k2)/(k1+k2+k3+k4)
end

FUNCTION Exponent, axis, index, number

  ; A special case.
  IF number EQ 0 THEN RETURN, '0'

  ; Assuming multiples of 10 with format.
  ex = String(number, Format='(e8.0)')
  pt = StrPos(ex, '.')

  first = StrMid(ex, 0, pt)
  sign = StrMid(ex, pt+2, 1)
  thisExponent = StrMid(ex, pt+3)

  ; Shave off leading zero in exponent
  WHILE StrMid(thisExponent, 0, 1) EQ '0' DO thisExponent = StrMid(thisExponent, 1)

  ; Fix for sign and missing zero problem.
  IF (Long(thisExponent) EQ 0) THEN BEGIN
    sign = ''
    thisExponent = '0'
  ENDIF

  ; Make the exponent a superscript.
  IF sign EQ '-' THEN BEGIN
    RETURN, first + 'x10!U' + sign + thisExponent + '!N'
  ENDIF ELSE BEGIN
    RETURN, first + 'x10!U' + thisExponent + '!N'
  ENDELSE

END

; Calculate tunneling rate
TSTART = 10 ;K
TINCR = 0.1
NTEMPS = 10000 ; number of temperature values
TEMPAR = FINDGEN(NTEMPS,INCREMENT=TINCR,START=TSTART)
KLHVAL = lamberts_tunnel(TEMPAR)

klhplot = PLOT($
  TEMPAR, $
  KLHVAL,$
  /DEVICE,$
  DIMENSIONS=[420,400],$
  MARGIN=[80,50,25,25],$
  FONT_SIZE=20,$
  XTICKFONT_SIZE=12,$
  YTICKFONT_SIZE=11,$
  YTITLE="k(T) (sec!E-1!N)",$
  XTITLE="Temperature (K)",$
  FONT_NAME="Hershey 3",$
  THICK=2,$
  XLOG=1,$
  YLOG=1,$
  XSTYLE=2,$
  YSTYLE=2$
;  XTICKFORMAT='exponent'$
  )

klhplot.SAVE, "/home/cns/Pictures/klhplot.eps",BORDER=kj0, RESOLUTION=1000

; Calculate rate of hopping over the barrier
EACT = 1.40e3 ; K for H + HOOH -> H2O + OH

HOPVAL = thermal_barrier(EACT,TEMPAR)

hopplot = PLOT($
  TEMPAR, $
  HOPVAL,$
  /DEVICE,$
  DIMENSIONS=[420,400],$
  MARGIN=[80,50,25,25],$
  FONT_SIZE=20,$
  XTICKFONT_SIZE=12,$
  YTICKFONT_SIZE=11,$
  YTITLE="k(T) (sec!E-1!N)",$
  XTITLE="Temperature (K)",$
  FONT_NAME="Hershey 3",$
  THICK=2,$
  XLOG=1,$
  YLOG=1,$
  XSTYLE=2,$
  YSTYLE=2$
  ;  XTICKFORMAT='exponent'$
  )

hopplot.SAVE, "/home/cns/Pictures/hopplot.eps",BORDER=kj0, RESOLUTION=1000

; Calculate kappa, the reaction/diffusion prabability
EH = 230 ; K
EHOOH = 5700/2 ; K

KAPVAL = kappa(EACT,EH,EHOOH,TEMPAR)

kapplot = PLOT($
  TEMPAR, $
  KAPVAL,$
  /DEVICE,$
  DIMENSIONS=[420,400],$
  MARGIN=[100,50,25,25],$
  FONT_SIZE=20,$
  XTICKFONT_SIZE=12,$
  YTICKFONT_SIZE=11,$
  YTITLE="Kappa",$
  XTITLE="Temperature (K)",$
  FONT_NAME="Hershey 3",$
  THICK=2,$
  XLOG=1,$
  YLOG=1,$
  XSTYLE=2,$
  YSTYLE=2,$
  YTICKFORMAT='exponent'$
  )

kapplot.SAVE, "/home/cns/Pictures/kappplot.eps",BORDER=kj0, RESOLUTION=1000

PRINT, "Ending script"

END