Ian Jauslin
summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorIan Jauslin <jauslin@ias.edu>2019-05-22 22:59:03 -0400
committerIan Jauslin <jauslin@ias.edu>2019-05-22 22:59:03 -0400
commitdf7449e4a29ec8d3182cf7b2aebcb86f7ac596c2 (patch)
tree182c5c85f57264b38c60a02ec268e5c00d937b34 /figs/animation/animate_compute.jl
As presented at VirginiaTech on 2019-05-24HEADv1.0master
Diffstat (limited to 'figs/animation/animate_compute.jl')
-rw-r--r--figs/animation/animate_compute.jl66
1 files changed, 66 insertions, 0 deletions
diff --git a/figs/animation/animate_compute.jl b/figs/animation/animate_compute.jl
new file mode 100644
index 0000000..452320e
--- /dev/null
+++ b/figs/animation/animate_compute.jl
@@ -0,0 +1,66 @@
+using QuadGK
+using SpecialFunctions
+using FFTW
+
+# numerical values
+hbar=6.58e-16 # eV.s
+m=9.11e-31 # kg
+Vn=9 # eV
+#En=14e9 # V/m
+En=10e9 # V/m
+Kn=4.5 # eV
+
+V=1
+E=En*hbar/(2*Vn^1.5*m^0.5)*sqrt(1.60e-19)
+k0=sqrt(Kn/Vn)
+
+# rescale x to nm
+nm_scale=hbar*sqrt(1.6e-19)/sqrt(2*m*Vn)*1e9
+
+# cutoffs
+p_cutoff=20*k0
+p_npoints=256
+
+# airy approximations
+airy_threshold=30
+airy_order=5
+
+# xbounds
+xmax=10
+xmin=-10
+x_npoints=200
+
+include("FN_base.jl")
+
+# compute wave function
+ps=Array{Array{Array{Complex{Float64},1},1}}(undef,x_npoints)
+dps=Array{Array{Array{Complex{Float64},1},1}}(undef,x_npoints)
+for i in 1:x_npoints
+ print(stderr,i,'/',x_npoints,'\n')
+ x=xmin+(xmax-xmin)*i/x_npoints
+ ps[i]=psi(x,k0,E,V,p_npoints,p_cutoff)
+ dps[i]=dpsi(x,k0,E,V,p_npoints,p_cutoff)
+end
+
+# compute asymptotic values
+ps_asym=Array{Complex{Float64}}(undef,x_npoints)
+dps_asym=Array{Complex{Float64}}(undef,x_npoints)
+for i in 1:x_npoints
+ x=xmin+(xmax-xmin)*i/x_npoints
+ ps_asym[i]=psi_pole(x,k0,E,V)
+ dps_asym[i]=dpsi_pole(x,k0,E,V)
+end
+
+# print asymptotic values
+for i in 1:x_npoints
+ print((xmin+(xmax-xmin)*i/x_npoints)*nm_scale,' ',abs(ps_asym[i])^2,' ',J(ps_asym[i],dps_asym[i])/(2*k0),'\n')
+end
+print('\n')
+
+# print values at each time
+for j in 1:p_npoints
+ for i in 1:x_npoints
+ print(real(ps[i][1][j])*hbar/Vn*1e15,' ',(xmin+(xmax-xmin)*i/x_npoints)*nm_scale,' ',abs(ps[i][2][j])^2,' ',J(ps[i][2][j],dps[i][2][j])/(2*k0),'\n')
+ end
+ print('\n')
+end