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=sqrt(2*m*Vn)/hbar*1e9*sqrt(1.60e-19)

# 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