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