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
|
using QuadGK
using FastGaussQuadrature
using SpecialFunctions
using FFTW
# numerical values
hbar=6.58e-16 # eV.s
m=9.11e-31 # kg
Un=9 # eV
En=parse(Float64,ARGS[1])*1e9 # V/m
Kn=4.5 # eV
# dimensionless quantities
U=1
E=En*hbar/(Un^1.5*m^0.5)*sqrt(1.60e-19)
k0=sqrt(2*Kn/Un)
# cutoffs
p_cutoff=20*k0
p_npoints=4096
# airy approximations
airy_threshold=30
airy_order=5
# order for Gauss-Legendre quadrature
order=10
# compute at these points
X=[(2*U-k0*k0)/(2*E),10*(2*U-k0*k0)/(2*E)]
include("FN_base.jl")
# compute the weights and abcissa for gauss-legendre quadratures
gl_data=gausslegendre(order)
ps=Array{Array{Array{Complex{Float64}}}}(undef,length(X))
dps=Array{Array{Array{Complex{Float64}}}}(undef,length(X))
intJ=Array{Array{Complex{Float64}}}(undef,length(X))
for i in 1:length(X)
# wave function
ps[i]=psi(X[i],k0,E,U,p_npoints,p_cutoff)
dps[i]=dpsi(X[i],k0,E,U,p_npoints,p_cutoff)
# integrated current
intJ[i]=zeros(Complex{Float64},p_npoints)
for l in 1:order
eval=current(X[i],k0/2*(gl_data[1][l]+1),E,U,p_npoints,p_cutoff)
for j in 1:length(eval)
intJ[i][j]=intJ[i][j]+k0/2*gl_data[2][l]*eval[j]
end
end
end
for j in 1:p_npoints
for i in 1:length(X)
print(real(ps[i][1][j])*hbar/Un*1e15,' ',abs(ps[i][2][j])^2,' ',J(ps[i][2][j],dps[i][2][j])/(2*k0),' ',real(intJ[i][j]/k0^2),' ')
end
print('\n')
end
|
