1 files changed, 60 insertions, 0 deletions
diff --git a/figs/plots.fig/current_density.jl b/figs/plots.fig/current_density.jl
new file mode 100644
index 0000000..9ec85f6
--- /dev/null
+++ b/figs/plots.fig/current_density.jl
@@ -0,0 +1,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