17 files changed, 885 insertions, 0 deletions
diff --git a/figs/Millikan-Lauritsen_current.png b/figs/Millikan-Lauritsen_current.png
new file mode 100644
index 0000000..5645be8
--- /dev/null
+++ b/figs/Millikan-Lauritsen_current.png
diff --git a/figs/animation/FN_base.jl b/figs/animation/FN_base.jl
new file mode 100644
index 0000000..af2a1ee
--- /dev/null
+++ b/figs/animation/FN_base.jl
@@ -0,0 +1,170 @@
+# fractional power with an arbitrary branch cut
+function pow(x,a,cut)
+  if(angle(x)/cut<=1)
+    return(abs(x)^a*exp(1im*angle(x)*a))
+  else
+    return(abs(x)^a*exp(1im*(angle(x)-sign(cut)*2*pi)*a))
+  end
+end
+
+# asymptotic airy functions
+# specify a branch cut for the fractional power
+function airyai_asym(x,cut)
+  if(abs(real(pow(x,3/2,cut)))<airy_threshold)
+    return(exp(2/3*pow(x,3/2,cut))*airyai(x))
+  else
+    ret=0
+    for n in 0:airy_order
+      ret+=gamma(n+5/6)*gamma(n+1/6)*(-3/4)^n/(4*pi^(3/2)*factorial(n)*pow(x,3*n/2+1/4,cut))
+    end
+    return ret
+  end
+end
+function airyaiprime_asym(x,cut)
+  if(abs(real(pow(x,3/2,cut)))<airy_threshold)
+    return(exp(2/3*pow(x,3/2,cut))*airyaiprime(x))
+  else
+    ret=0
+    for n in 0:airy_order
+      ret+=gamma(n+5/6)*gamma(n+1/6)*(-3/4)^n/(4*pi^(3/2)*factorial(n))*(-1/pow(x,3*n/2-1/4,cut)-(3/2*n+1/4)/pow(x,3*n/2+5/4,cut))
+    end
+    return ret
+  end
+end
+
+# solutions of (-\Delta+V-ip)phi=0
+# assume that p has an infinitesimal real part (and adjust the branch cuts appropriately)
+function phi(p,x,E,V)
+  return(airyai_asym(exp(-1im*pi/3)*(E^(1/3)*x-E^(-2/3)*(V-1im*p)),pi))
+end
+function dphi(p,x,E,V)
+  return(exp(-1im*pi/3)*E^(1/3)*airyaiprime_asym(exp(-1im*pi/3)*(E^(1/3)*x-E^(-2/3)*(V-1im*p)),pi))
+end
+function eta(p,x,E,V)
+  return(exp(-1im*pi/3)*airyai_asym(-(E^(1/3)*x-E^(-2/3)*(V-1im*p)),pi/2))
+end
+function deta(p,x,E,V)
+  return(-exp(-1im*pi/3)*E^(1/3)*airyaiprime_asym(-(E^(1/3)*x-E^(-2/3)*(V-1im*p)),pi/2))
+end
+
+# Laplace transform of psi
+# assume that p has an infinitesimal real part (and adjust the branch cuts appropriately)
+# for example, (1im*p-V)^(3/2) becomes pow(1im*p-V,3/2,-pi/2) because when 1im*p is real negative, its square root should be imaginary positive
+function f(p,x,k0,E,V)
+  T=2im*k0/(1im*k0-sqrt(V-k0*k0))
+  R=T-1
+
+  if x>=0
+    C2=-1im*T/(pow(-1im*p,1/2,pi/2)*phi(p,0,E,V)-dphi(p,0,E,V))*((sqrt(V-k0*k0)+pow(-1im*p,1/2,pi/2))/(-1im*p+k0*k0)-2im*E^(-1/3)*pi*quadgk(y -> (pow(-1im*p,1/2,pi/2)*eta(p,0,E,V)-deta(p,0,E,V))*phi(p,y,E,V)*exp(-sqrt(V-k0*k0)*y)*exp(2im/3*(pow(E^(1/3)*y+E^(-2/3)*(1im*p-V),3/2,-pi/2)-E^(-1)*pow(1im*p-V,3/2,-pi/2))),0,Inf)[1])
+    FT=2*E^(-1/3)*pi*(quadgk(y -> phi(p,x,E,V)*eta(p,y,E,V)*exp(-sqrt(V-k0*k0)*y)*exp(2im/3*(pow(E^(1/3)*x+E^(-2/3)*(1im*p-V),3/2,-pi/2)-pow(E^(1/3)*y+E^(-2/3)*(1im*p-V),3/2,-pi/2))),0,x)[1]+quadgk(y -> eta(p,x,E,V)*phi(p,y,E,V)*exp(-sqrt(V-k0*k0)*y)*exp(2im/3*(pow(E^(1/3)*y+E^(-2/3)*(1im*p-V),3/2,-pi/2)-pow(E^(1/3)*x+E^(-2/3)*(1im*p-V),3/2,-pi/2))),x,Inf)[1])
+    main=C2*phi(p,x,E,V)*exp(2im/3*(pow(E^(1/3)*x+E^(-2/3)*(1im*p-V),3/2,-pi/2)-E^(-1)*pow(1im*p-V,3/2,-pi/2)))+T*FT
+
+    # subtract the contribution of the pole, which will be added back in after the integration
+    pole=psi_pole(x,k0,E,V)/(p+1im*k0*k0)
+    return(main-pole)
+  else
+    C1=-1im*T*((sqrt(V-k0*k0)*phi(p,0,E,V)+dphi(p,0,E,V))/(-1im*p+k0*k0)/(pow(-1im*p,1/2,pi/2)*phi(p,0,E,V)-dphi(p,0,E,V))+E^(-1/3)*quadgk(y -> phi(p,y,E,V)/(pow(-1im*p,1/2,pi/2)*phi(p,0,E,V)-dphi(p,0,E,V))*exp(-sqrt(V-k0*k0)*y)*exp(2im/3*(pow(E^(1/3)*y+E^(-2/3)*(1im*p-V),3/2,-pi/2)-E^(-1)*pow(1im*p-V,3/2,-pi/2))),0,Inf)[1])
+    FI=-1im*exp(1im*k0*x)/(-1im*p+k0*k0)
+    FR=-1im*exp(-1im*k0*x)/(-1im*p+k0*k0)
+    main=C1*exp(pow(-1im*p,1/2,pi/2)*x)+FI+R*FR
+
+    # subtract the contribution of the pole, which will be added back in after the integration
+    pole=psi_pole(x,k0,E,V)/(p+1im*k0*k0)
+    return(main-pole)
+  end
+end
+# its derivative
+function df(p,x,k0,E,V)
+  T=2im*k0/(1im*k0-sqrt(V-k0*k0))
+  R=T-1
+
+  if x>=0
+    C2=-1im*T/(pow(-1im*p,1/2,pi/2)*phi(p,0,E,V)-dphi(p,0,E,V))*((sqrt(V-k0*k0)+pow(-1im*p,1/2,pi/2))/(-1im*p+k0*k0)-2im*E^(-1/3)*pi*quadgk(y -> (pow(-1im*p,1/2,pi/2)*eta(p,0,E,V)-deta(p,0,E,V))*phi(p,y,E,V)*exp(-sqrt(V-k0*k0)*y)*exp(2im/3*(pow(E^(1/3)*y+E^(-2/3)*(1im*p-V),3/2,-pi/2)-E^(-1)*pow(1im*p-V,3/2,-pi/2))),0,Inf)[1])
+    dFT=2*E^(-1/3)*pi*(quadgk(y -> dphi(p,x,E,V)*eta(p,y,E,V)*exp(-sqrt(V-k0*k0)*y)*exp(2im/3*(pow(E^(1/3)*x+E^(-2/3)*(1im*p-V),3/2,-pi/2)-pow(E^(1/3)*y+E^(-2/3)*(1im*p-V),3/2,-pi/2))),0,x)[1]+quadgk(y -> deta(p,x,E,V)*phi(p,y,E,V)*exp(-sqrt(V-k0*k0)*y)*exp(2im/3*(pow(E^(1/3)*y+E^(-2/3)*(1im*p-V),3/2,-pi/2)-pow(E^(1/3)*x+E^(-2/3)*(1im*p-V),3/2,-pi/2))),x,Inf)[1])
+    main=C2*dphi(p,x,E,V)*exp(2im/3*(pow(E^(1/3)*x+E^(-2/3)*(1im*p-V),3/2,-pi/2)-E^(-1)*pow(1im*p-V,3/2,-pi/2)))+T*dFT
+
+    # subtract the contribution of the pole, which will be added back in after the integration
+    pole=dpsi_pole(x,k0,E,V)/(p+1im*k0*k0)
+    return(main-pole)
+  else
+    C1=-1im*T*((sqrt(V-k0*k0)*phi(p,0,E,V)+dphi(p,0,E,V))/(-1im*p+k0*k0)/(pow(-1im*p,1/2,pi/2)*phi(p,0,E,V)-dphi(p,0,E,V))+E^(-1/3)*quadgk(y -> phi(p,y,E,V)/(pow(-1im*p,1/2,pi/2)*phi(p,0,E,V)-dphi(p,0,E,V))*exp(-sqrt(V-k0*k0)*y)*exp(2im/3*(pow(E^(1/3)*y+E^(-2/3)*(1im*p-V),3/2,-pi/2)-E^(-1)*pow(1im*p-V,3/2,-pi/2))),0,Inf)[1])
+    dFI=k0*exp(1im*k0*x)/(-1im*p+k0*k0)
+    dFR=-k0*exp(-1im*k0*x)/(-1im*p+k0*k0)
+    main=C1*pow(-1im*p,1/2,pi/2)*exp(pow(-1im*p,1/2,pi/2)*x)+dFI+R*dFR
+
+    # subtract the contribution of the pole, which will be added back in after the integration
+    pole=dpsi_pole(x,k0,E,V)/(p+1im*k0*k0)
+    return(main-pole)
+  end
+end
+
+# psi (returns t,psi(x,t))
+function psi(x,k0,E,V,p_npoints,p_cutoff)
+  fft=fourier_fft(f,x,k0,E,V,p_npoints,p_cutoff)
+  # add the contribution of the pole
+  for i in 1:p_npoints
+      fft[2][i]=fft[2][i]+psi_pole(x,k0,E,V)*exp(-1im*k0*k0*fft[1][i])
+  end
+  return(fft)
+end
+# its derivative
+function dpsi(x,k0,E,V,p_npoints,p_cutoff)
+  fft=fourier_fft(df,x,k0,E,V,p_npoints,p_cutoff)
+  # add the contribution of the pole
+  for i in 1:p_npoints
+    fft[2][i]=fft[2][i]+dpsi_pole(x,k0,E,V)*exp(-1im*k0*k0*fft[1][i])
+  end
+  return(fft)
+end
+
+# compute Fourier transform by sampling and fft
+function fourier_fft(A,x,k0,E,V,p_npoints,p_cutoff)
+  fun=zeros(Complex{Float64},p_npoints)
+  times=zeros(p_npoints)
+
+  # prepare fft
+  for i in 1:p_npoints
+    fun[i]=p_cutoff/pi*A(1im*(-p_cutoff+2*p_cutoff*(i-1)/p_npoints),x,k0,E,V)
+    times[i]=(i-1)*pi/p_cutoff
+  end
+
+  ifft!(fun)
+  
+  # correct the phase
+  for i in 2:2:p_npoints
+    fun[i]=-fun[i]
+  end
+  return([times,fun])
+end
+
+# asymptotic value of psi
+function psi_pole(x,k0,E,V)
+  if x>=0
+    return(1im*phi(-1im*k0*k0,x,E,V)*2*k0/(1im*k0*phi(-1im*k0*k0,0,E,V)+dphi(-1im*k0*k0,0,E,V))*exp(2im/3*(pow(E^(1/3)*x+E^(-2/3)*(k0*k0-V),3/2,-pi/2)-E^(-1)*pow(k0*k0-V,3/2,-pi/2))))
+  else
+    return((1im*k0*phi(-1im*k0*k0,0,E,V)-dphi(-1im*k0*k0,0,E,V))/(1im*k0*phi(-1im*k0*k0,0,E,V)+dphi(-1im*k0*k0,0,E,V))*exp(-1im*k0*x)+exp(1im*k0*x))
+  end
+end
+function dpsi_pole(x,k0,E,V)
+  if x>=0
+    return(1im*dphi(-1im*k0*k0,x,E,V)*2*k0/(1im*k0*phi(-1im*k0*k0,0,E,V)+dphi(-1im*k0*k0,0,E,V))*exp(2im/3*(pow(E^(1/3)*x+E^(-2/3)*(k0*k0-V),3/2,-pi/2)-E^(-1)*pow(k0*k0-V,3/2,-pi/2))))
+  else
+    return(-1im*k0*(1im*k0*phi(-1im*k0*k0,0,E,V)-dphi(-1im*k0*k0,0,E,V))/(1im*k0*phi(-1im*k0*k0,0,E,V)+dphi(-1im*k0*k0,0,E,V))*exp(-1im*k0*x)+1im*k0*exp(1im*k0*x))
+  end
+end
+
+# current
+function J(ps,dps)
+  return(2*imag(conj(ps)*dps))
+end
+
+# complete computation of the current
+function current(x,k0,E,V,p_npoints,p_cutoff)
+  ps=psi(x,k0,E,V,p_npoints,p_cutoff)
+  dps=dpsi(x,k0,E,V,p_npoints,p_cutoff)
+  Js=zeros(Complex{Float64},p_npoints)
+  for i in 1:p_npoints
+    Js[i]=J(ps[2][i],dps[2][i])
+  end
+  return(Js)
+end
diff --git a/figs/animation/Makefile b/figs/animation/Makefile
new file mode 100644
index 0000000..29098b3
--- /dev/null
+++ b/figs/animation/Makefile
@@ -0,0 +1,12 @@
+PROJECTNAME=animate
+
+all: animate.dat
+
+run: animate.dat
+	python3 animate.py animate.dat
+
+animate.dat:
+	julia animate_compute.jl > animate.dat
+
+clean:
+	rm -f animate.dat
diff --git a/figs/animation/animate.py b/figs/animation/animate.py
new file mode 100644
index 0000000..6203fb7
--- /dev/null
+++ b/figs/animation/animate.py
@@ -0,0 +1,128 @@
+from matplotlib import pyplot as pl
+from matplotlib import animation
+import sys
+
+# read data
+# time dependent data
+frames=[]
+# asymptotic data (located in the first block)
+asym=[]
+infile=open(sys.argv[1],'r')
+row=[]
+for line in infile:
+    # read first block
+    if len(asym)==0:
+        if line=='\n':
+            asym=row
+            row=[]
+        else:
+            dat=[]
+            for n in line.split():
+                dat.append(float(n))
+            row.append(dat)
+    # read other blocks
+    else:
+        if line=='\n':
+            frames.append(row)
+            row=[]
+        else:
+            dat=[]
+            for n in line.split():
+                dat.append(float(n))
+            row.append(dat)
+infile.close()
+
+
+# set up plot
+fig = pl.figure()
+pl.subplot(211)
+axr=fig.gca()
+asym_rho, = axr.plot([],[],linewidth=3.5,color='#00FF00')
+rho, = axr.plot([],[],color='red')
+
+pl.subplot(212)
+axJ=fig.gca()
+asym_J, = axJ.plot([],[],linewidth=3.5,color='#00FF00')
+J, = axJ.plot([],[],color='red')
+
+# plot ranges
+xmax=0
+maxyr=0
+maxyJ=0
+for frame in frames:
+    for i in range(len(frame)):
+        if frame[i][1]>xmax:
+            xmax=frame[i][1]
+        if frame[i][2]>maxyr:
+            maxyr=frame[i][2]
+        if frame[i][3]>maxyJ:
+            maxyJ=frame[i][3]
+for i in range(len(asym)):
+    if asym[i][0]>xmax:
+        xmax=asym[i][0]
+    if asym[i][1]>maxyr:
+        maxyr=asym[i][1]
+    if asym[i][2]>maxyJ:
+        maxyJ=asym[i][2]
+xmin=0
+minyr=0
+minyJ=0
+for frame in frames:
+    for i in range(len(frame)):
+        if frame[i][1]<xmin:
+            xmin=frame[i][1]
+        if frame[i][2]<minyr:
+            minyr=frame[i][2]
+        if frame[i][3]<minyJ:
+            minyJ=frame[i][3]
+for i in range(len(asym)):
+    if asym[i][0]<xmin:
+        xmin=asym[i][0]
+    if asym[i][1]<minyr:
+        minyr=asym[i][1]
+    if asym[i][2]<minyJ:
+        minyJ=asym[i][2]
+
+
+# plot asymptotes
+asym_rho_datax=[]
+asym_rho_datay=[]
+for i in range(len(asym)):
+    asym_rho_datax.append(asym[i][0])
+    asym_rho_datay.append(asym[i][1])
+asym_rho.set_data(asym_rho_datax,asym_rho_datay)
+asym_J_datax=[]
+asym_J_datay=[]
+for i in range(len(asym)):
+    asym_J_datax.append(asym[i][0])
+    asym_J_datay.append(asym[i][2])
+asym_J.set_data(asym_J_datax,asym_J_datay)
+
+# animate
+def init_plot():
+    axr.set_ylim(minyr,maxyr)
+    axr.set_xlim(xmin,xmax)
+    axJ.set_ylim(minyJ,maxyJ)
+    axJ.set_xlim(xmin,xmax)
+
+    axr.vlines(0,minyr,maxyr,linestyles="dotted")
+    axJ.vlines(0,minyJ,maxyJ,linestyles="dotted")
+    return rho,J
+def update(frame):
+    axr.set_title("t=% .3f fs" % (frame[0][0]))
+    xdata=[]
+    ydata=[]
+    for i in range(len(frame)):
+        xdata.append(frame[i][1])
+        ydata.append(frame[i][2])
+    rho.set_data(xdata,ydata)
+
+    xdata=[]
+    ydata=[]
+    for i in range(len(frame)):
+        xdata.append(frame[i][1])
+        ydata.append(frame[i][3])
+    J.set_data(xdata,ydata)
+    return rho,J
+anim = animation.FuncAnimation(fig, update, frames=frames, blit=False, interval=100, repeat=True, init_func=init_plot)
+pl.show()
diff --git a/figs/animation/animate_compute.jl b/figs/animation/animate_compute.jl
new file mode 100644
index 0000000..b9c2d02
--- /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=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
diff --git a/figs/contour.fig/Makefile b/figs/contour.fig/Makefile
new file mode 100644
index 0000000..bade7d0
--- /dev/null
+++ b/figs/contour.fig/Makefile
@@ -0,0 +1,30 @@
+PROJECTNAME=contour
+
+PDFS=$(addsuffix .pdf, $(PROJECTNAME))
+
+all: $(PDFS)
+
+$(PDFS): poles.tikz.tex
+	pdflatex -jobname $(basename $@) -file-line-error $(patsubst %.pdf, %.tikz.tex, $@)
+
+poles.tikz.tex:
+	python3 contour.py > poles.tikz.tex
+
+install: $(PDFS)
+	cp $^ $(INSTALLDIR)/
+
+$(LIBS):
+	ln -fs libs/$@ ./
+
+clean-libs:
+	rm -f $(LIBS)
+
+clean-aux:
+	rm -f $(addsuffix .aux, $(PROJECTNAME))
+	rm -f $(addsuffix .log, $(PROJECTNAME))
+	rm -f poles.tikz.tex
+
+clean-tex:
+	rm -f $(PDFS)
+
+clean: clean-libs clean-aux clean-tex
diff --git a/figs/contour.fig/contour.py b/figs/contour.fig/contour.py
new file mode 100644
index 0000000..102df86
--- /dev/null
+++ b/figs/contour.fig/contour.py
@@ -0,0 +1,82 @@
+#!/usr/bin/env python3
+
+import cmath
+import math
+import scipy.special as sp
+from scipy import optimize as op
+import random
+import sys
+
+# number of roots
+nr_roots=4
+
+# size of plot
+plotsize_x=3
+plotsize_y=3
+# rescale plot (so roots are not too close together)
+plotsize_scale_x=1
+plotsize_scale_y=6
+
+# numerical values
+hbar=6.58e-16 # eV.s
+m=9.11e-31 # kg
+Vn=9 # eV
+En=20e9 # V/m
+
+V=1
+E=En*hbar/(Vn**1.5*m**0.5)*math.sqrt(1.60e-19)
+
+# sqrt with branch cut along iR_+
+def sqrt_p(x):
+    r,p=cmath.polar(x)
+    if(p<cmath.pi/2):
+        return(cmath.rect(math.sqrt(r),p/2))
+    else:
+        return(cmath.rect(math.sqrt(r),(p-2*math.pi)/2))
+
+# solution of (-\Delta+V-ip)phi=0
+def phi(p,x,E,V):
+    return(sp.airy(cmath.exp(-1j*math.pi/3)*(E**(1/3)*x-E**(-2/3)*(V-1j*p)))[0])
+# its derivative
+def dphi(p,x,E,V):
+    return(cmath.exp(-1j*math.pi/3)*E**(1/3)*sp.airy(cmath.exp(-1j*math.pi/3)*(E**(1/3)*x-E**(-2/3)*(V-1j*p)))[1])
+
+
+# the function whose roots are to be computed and its derivative
+def f(p):
+    return(sqrt_p(-1j*p)*phi(p,0,E,V)-dphi(p,0,E,V))
+def df(p):
+    return(-1j/(2*sqrt_p(-1j*p))*phi(p,0,E,V)+sqrt_p(-1j*p)*dp_phi(p,0,E,V)-dp_dphi(p,0,E,V))
+
+# derivatives of phi with respect to p
+def dp_phi(p,x,E,V):
+    return(1j*cmath.exp(-1j*math.pi/3)*E**(-2/3)*sp.airy(cmath.exp(-1j*math.pi/3)*(E**(1/3)*x-E**(-2/3)*(V-1j*p)))[1])
+def dp_dphi(p,x,E,V):
+    return(-1j*(x-(V-1j*p)/E)*phi(p,x,E,V))
+
+# check whether the root was already found
+def check(root,roots):
+    # reject if the root doesn't fit on the plot
+    if(plotsize_scale_x*root.real<-plotsize_x or abs(plotsize_scale_y*root.imag)>plotsize_y):
+        return(False)
+    for x in roots:
+        if(abs(root-x)<1e-12):
+            return(False)
+    return(True)
+
+# find roots
+roots=[]
+while len(roots)<nr_roots:
+    try:
+        root=op.newton(f, -random.random()*plotsize_x/plotsize_scale_x+(2*random.random()-1)*plotsize_y/plotsize_scale_y*1j, fprime=df)
+        if(check(root,roots)):
+            roots.append(root)
+            print("found "+str(len(roots))+" roots of "+str(nr_roots), file=sys.stderr)
+    except RuntimeError:
+        root=0
+    except:
+        break
+
+# print result
+for root in roots:
+    print("\\pole{(% .3f,% .3f)}" % (plotsize_scale_x*root.real,plotsize_scale_y*root.imag))
diff --git a/figs/contour.fig/contour.tikz.tex b/figs/contour.fig/contour.tikz.tex
new file mode 100644
index 0000000..953d8f0
--- /dev/null
+++ b/figs/contour.fig/contour.tikz.tex
@@ -0,0 +1,34 @@
+\documentclass{standalone}
+
+\usepackage{tikz}
+
+\def\pole#1{
+  \draw#1++(0.1,0.1)--++(-0.2,-0.2);
+  \draw#1++(-0.1,0.1)--++(0.2,-0.2);
+  \draw[color=red, line width=1pt]#1circle(0.3);
+}
+
+\begin{document}
+\begin{tikzpicture}
+
+% branch cut
+\fill[color=gray](-3.3,-0.1)--++(3.3,0)--++(0,0.2)--++(-3.3,0)--cycle;
+
+% axes
+\draw(-3.3,0)--++(6.6,0);
+\draw(0,-3.3)--++(0,6.6);
+
+% -ik0^2
+\pole{(0,-2)}
+
+% other poles
+\input{poles.tikz}
+
+% contour
+\draw[color=red, line width=1pt](-3.3,-0.3)--++(3.3,0);
+\draw[color=red, line width=1pt](0,-0.3)..controls(0.3,-0.3)and(0.3,0.3)..(0,0.3);
+\draw[color=red, line width=1pt](-3.3,0.3)--++(3.3,0);
+
+
+\end{tikzpicture}
+\end{document}
diff --git a/figs/emitter.jpg b/figs/emitter.jpg
new file mode 100644
index 0000000..d84ad1f
--- /dev/null
+++ b/figs/emitter.jpg
diff --git a/figs/fowler-nordheim.fig/FN_base.jl b/figs/fowler-nordheim.fig/FN_base.jl
new file mode 100644
index 0000000..af2a1ee
--- /dev/null
+++ b/figs/fowler-nordheim.fig/FN_base.jl
@@ -0,0 +1,170 @@
+# fractional power with an arbitrary branch cut
+function pow(x,a,cut)
+  if(angle(x)/cut<=1)
+    return(abs(x)^a*exp(1im*angle(x)*a))
+  else
+    return(abs(x)^a*exp(1im*(angle(x)-sign(cut)*2*pi)*a))
+  end
+end
+
+# asymptotic airy functions
+# specify a branch cut for the fractional power
+function airyai_asym(x,cut)
+  if(abs(real(pow(x,3/2,cut)))<airy_threshold)
+    return(exp(2/3*pow(x,3/2,cut))*airyai(x))
+  else
+    ret=0
+    for n in 0:airy_order
+      ret+=gamma(n+5/6)*gamma(n+1/6)*(-3/4)^n/(4*pi^(3/2)*factorial(n)*pow(x,3*n/2+1/4,cut))
+    end
+    return ret
+  end
+end
+function airyaiprime_asym(x,cut)
+  if(abs(real(pow(x,3/2,cut)))<airy_threshold)
+    return(exp(2/3*pow(x,3/2,cut))*airyaiprime(x))
+  else
+    ret=0
+    for n in 0:airy_order
+      ret+=gamma(n+5/6)*gamma(n+1/6)*(-3/4)^n/(4*pi^(3/2)*factorial(n))*(-1/pow(x,3*n/2-1/4,cut)-(3/2*n+1/4)/pow(x,3*n/2+5/4,cut))
+    end
+    return ret
+  end
+end
+
+# solutions of (-\Delta+V-ip)phi=0
+# assume that p has an infinitesimal real part (and adjust the branch cuts appropriately)
+function phi(p,x,E,V)
+  return(airyai_asym(exp(-1im*pi/3)*(E^(1/3)*x-E^(-2/3)*(V-1im*p)),pi))
+end
+function dphi(p,x,E,V)
+  return(exp(-1im*pi/3)*E^(1/3)*airyaiprime_asym(exp(-1im*pi/3)*(E^(1/3)*x-E^(-2/3)*(V-1im*p)),pi))
+end
+function eta(p,x,E,V)
+  return(exp(-1im*pi/3)*airyai_asym(-(E^(1/3)*x-E^(-2/3)*(V-1im*p)),pi/2))
+end
+function deta(p,x,E,V)
+  return(-exp(-1im*pi/3)*E^(1/3)*airyaiprime_asym(-(E^(1/3)*x-E^(-2/3)*(V-1im*p)),pi/2))
+end
+
+# Laplace transform of psi
+# assume that p has an infinitesimal real part (and adjust the branch cuts appropriately)
+# for example, (1im*p-V)^(3/2) becomes pow(1im*p-V,3/2,-pi/2) because when 1im*p is real negative, its square root should be imaginary positive
+function f(p,x,k0,E,V)
+  T=2im*k0/(1im*k0-sqrt(V-k0*k0))
+  R=T-1
+
+  if x>=0
+    C2=-1im*T/(pow(-1im*p,1/2,pi/2)*phi(p,0,E,V)-dphi(p,0,E,V))*((sqrt(V-k0*k0)+pow(-1im*p,1/2,pi/2))/(-1im*p+k0*k0)-2im*E^(-1/3)*pi*quadgk(y -> (pow(-1im*p,1/2,pi/2)*eta(p,0,E,V)-deta(p,0,E,V))*phi(p,y,E,V)*exp(-sqrt(V-k0*k0)*y)*exp(2im/3*(pow(E^(1/3)*y+E^(-2/3)*(1im*p-V),3/2,-pi/2)-E^(-1)*pow(1im*p-V,3/2,-pi/2))),0,Inf)[1])
+    FT=2*E^(-1/3)*pi*(quadgk(y -> phi(p,x,E,V)*eta(p,y,E,V)*exp(-sqrt(V-k0*k0)*y)*exp(2im/3*(pow(E^(1/3)*x+E^(-2/3)*(1im*p-V),3/2,-pi/2)-pow(E^(1/3)*y+E^(-2/3)*(1im*p-V),3/2,-pi/2))),0,x)[1]+quadgk(y -> eta(p,x,E,V)*phi(p,y,E,V)*exp(-sqrt(V-k0*k0)*y)*exp(2im/3*(pow(E^(1/3)*y+E^(-2/3)*(1im*p-V),3/2,-pi/2)-pow(E^(1/3)*x+E^(-2/3)*(1im*p-V),3/2,-pi/2))),x,Inf)[1])
+    main=C2*phi(p,x,E,V)*exp(2im/3*(pow(E^(1/3)*x+E^(-2/3)*(1im*p-V),3/2,-pi/2)-E^(-1)*pow(1im*p-V,3/2,-pi/2)))+T*FT
+
+    # subtract the contribution of the pole, which will be added back in after the integration
+    pole=psi_pole(x,k0,E,V)/(p+1im*k0*k0)
+    return(main-pole)
+  else
+    C1=-1im*T*((sqrt(V-k0*k0)*phi(p,0,E,V)+dphi(p,0,E,V))/(-1im*p+k0*k0)/(pow(-1im*p,1/2,pi/2)*phi(p,0,E,V)-dphi(p,0,E,V))+E^(-1/3)*quadgk(y -> phi(p,y,E,V)/(pow(-1im*p,1/2,pi/2)*phi(p,0,E,V)-dphi(p,0,E,V))*exp(-sqrt(V-k0*k0)*y)*exp(2im/3*(pow(E^(1/3)*y+E^(-2/3)*(1im*p-V),3/2,-pi/2)-E^(-1)*pow(1im*p-V,3/2,-pi/2))),0,Inf)[1])
+    FI=-1im*exp(1im*k0*x)/(-1im*p+k0*k0)
+    FR=-1im*exp(-1im*k0*x)/(-1im*p+k0*k0)
+    main=C1*exp(pow(-1im*p,1/2,pi/2)*x)+FI+R*FR
+
+    # subtract the contribution of the pole, which will be added back in after the integration
+    pole=psi_pole(x,k0,E,V)/(p+1im*k0*k0)
+    return(main-pole)
+  end
+end
+# its derivative
+function df(p,x,k0,E,V)
+  T=2im*k0/(1im*k0-sqrt(V-k0*k0))
+  R=T-1
+
+  if x>=0
+    C2=-1im*T/(pow(-1im*p,1/2,pi/2)*phi(p,0,E,V)-dphi(p,0,E,V))*((sqrt(V-k0*k0)+pow(-1im*p,1/2,pi/2))/(-1im*p+k0*k0)-2im*E^(-1/3)*pi*quadgk(y -> (pow(-1im*p,1/2,pi/2)*eta(p,0,E,V)-deta(p,0,E,V))*phi(p,y,E,V)*exp(-sqrt(V-k0*k0)*y)*exp(2im/3*(pow(E^(1/3)*y+E^(-2/3)*(1im*p-V),3/2,-pi/2)-E^(-1)*pow(1im*p-V,3/2,-pi/2))),0,Inf)[1])
+    dFT=2*E^(-1/3)*pi*(quadgk(y -> dphi(p,x,E,V)*eta(p,y,E,V)*exp(-sqrt(V-k0*k0)*y)*exp(2im/3*(pow(E^(1/3)*x+E^(-2/3)*(1im*p-V),3/2,-pi/2)-pow(E^(1/3)*y+E^(-2/3)*(1im*p-V),3/2,-pi/2))),0,x)[1]+quadgk(y -> deta(p,x,E,V)*phi(p,y,E,V)*exp(-sqrt(V-k0*k0)*y)*exp(2im/3*(pow(E^(1/3)*y+E^(-2/3)*(1im*p-V),3/2,-pi/2)-pow(E^(1/3)*x+E^(-2/3)*(1im*p-V),3/2,-pi/2))),x,Inf)[1])
+    main=C2*dphi(p,x,E,V)*exp(2im/3*(pow(E^(1/3)*x+E^(-2/3)*(1im*p-V),3/2,-pi/2)-E^(-1)*pow(1im*p-V,3/2,-pi/2)))+T*dFT
+
+    # subtract the contribution of the pole, which will be added back in after the integration
+    pole=dpsi_pole(x,k0,E,V)/(p+1im*k0*k0)
+    return(main-pole)
+  else
+    C1=-1im*T*((sqrt(V-k0*k0)*phi(p,0,E,V)+dphi(p,0,E,V))/(-1im*p+k0*k0)/(pow(-1im*p,1/2,pi/2)*phi(p,0,E,V)-dphi(p,0,E,V))+E^(-1/3)*quadgk(y -> phi(p,y,E,V)/(pow(-1im*p,1/2,pi/2)*phi(p,0,E,V)-dphi(p,0,E,V))*exp(-sqrt(V-k0*k0)*y)*exp(2im/3*(pow(E^(1/3)*y+E^(-2/3)*(1im*p-V),3/2,-pi/2)-E^(-1)*pow(1im*p-V,3/2,-pi/2))),0,Inf)[1])
+    dFI=k0*exp(1im*k0*x)/(-1im*p+k0*k0)
+    dFR=-k0*exp(-1im*k0*x)/(-1im*p+k0*k0)
+    main=C1*pow(-1im*p,1/2,pi/2)*exp(pow(-1im*p,1/2,pi/2)*x)+dFI+R*dFR
+
+    # subtract the contribution of the pole, which will be added back in after the integration
+    pole=dpsi_pole(x,k0,E,V)/(p+1im*k0*k0)
+    return(main-pole)
+  end
+end
+
+# psi (returns t,psi(x,t))
+function psi(x,k0,E,V,p_npoints,p_cutoff)
+  fft=fourier_fft(f,x,k0,E,V,p_npoints,p_cutoff)
+  # add the contribution of the pole
+  for i in 1:p_npoints
+      fft[2][i]=fft[2][i]+psi_pole(x,k0,E,V)*exp(-1im*k0*k0*fft[1][i])
+  end
+  return(fft)
+end
+# its derivative
+function dpsi(x,k0,E,V,p_npoints,p_cutoff)
+  fft=fourier_fft(df,x,k0,E,V,p_npoints,p_cutoff)
+  # add the contribution of the pole
+  for i in 1:p_npoints
+    fft[2][i]=fft[2][i]+dpsi_pole(x,k0,E,V)*exp(-1im*k0*k0*fft[1][i])
+  end
+  return(fft)
+end
+
+# compute Fourier transform by sampling and fft
+function fourier_fft(A,x,k0,E,V,p_npoints,p_cutoff)
+  fun=zeros(Complex{Float64},p_npoints)
+  times=zeros(p_npoints)
+
+  # prepare fft
+  for i in 1:p_npoints
+    fun[i]=p_cutoff/pi*A(1im*(-p_cutoff+2*p_cutoff*(i-1)/p_npoints),x,k0,E,V)
+    times[i]=(i-1)*pi/p_cutoff
+  end
+
+  ifft!(fun)
+  
+  # correct the phase
+  for i in 2:2:p_npoints
+    fun[i]=-fun[i]
+  end
+  return([times,fun])
+end
+
+# asymptotic value of psi
+function psi_pole(x,k0,E,V)
+  if x>=0
+    return(1im*phi(-1im*k0*k0,x,E,V)*2*k0/(1im*k0*phi(-1im*k0*k0,0,E,V)+dphi(-1im*k0*k0,0,E,V))*exp(2im/3*(pow(E^(1/3)*x+E^(-2/3)*(k0*k0-V),3/2,-pi/2)-E^(-1)*pow(k0*k0-V,3/2,-pi/2))))
+  else
+    return((1im*k0*phi(-1im*k0*k0,0,E,V)-dphi(-1im*k0*k0,0,E,V))/(1im*k0*phi(-1im*k0*k0,0,E,V)+dphi(-1im*k0*k0,0,E,V))*exp(-1im*k0*x)+exp(1im*k0*x))
+  end
+end
+function dpsi_pole(x,k0,E,V)
+  if x>=0
+    return(1im*dphi(-1im*k0*k0,x,E,V)*2*k0/(1im*k0*phi(-1im*k0*k0,0,E,V)+dphi(-1im*k0*k0,0,E,V))*exp(2im/3*(pow(E^(1/3)*x+E^(-2/3)*(k0*k0-V),3/2,-pi/2)-E^(-1)*pow(k0*k0-V,3/2,-pi/2))))
+  else
+    return(-1im*k0*(1im*k0*phi(-1im*k0*k0,0,E,V)-dphi(-1im*k0*k0,0,E,V))/(1im*k0*phi(-1im*k0*k0,0,E,V)+dphi(-1im*k0*k0,0,E,V))*exp(-1im*k0*x)+1im*k0*exp(1im*k0*x))
+  end
+end
+
+# current
+function J(ps,dps)
+  return(2*imag(conj(ps)*dps))
+end
+
+# complete computation of the current
+function current(x,k0,E,V,p_npoints,p_cutoff)
+  ps=psi(x,k0,E,V,p_npoints,p_cutoff)
+  dps=dpsi(x,k0,E,V,p_npoints,p_cutoff)
+  Js=zeros(Complex{Float64},p_npoints)
+  for i in 1:p_npoints
+    Js[i]=J(ps[2][i],dps[2][i])
+  end
+  return(Js)
+end
diff --git a/figs/fowler-nordheim.fig/Makefile b/figs/fowler-nordheim.fig/Makefile
new file mode 100644
index 0000000..a6e5e53
--- /dev/null
+++ b/figs/fowler-nordheim.fig/Makefile
@@ -0,0 +1,29 @@
+PROJECTNAME=asymptotic
+
+PDFS=$(addsuffix .pdf, $(PROJECTNAME))
+TEXS=$(addsuffix .tikz.tex, $(PROJECTNAME))
+
+all: $(PDFS)
+
+$(PDFS): $(addsuffix .dat, $(PROJECTNAME))
+	gnuplot $(patsubst %.pdf, %.gnuplot, $@) > $(patsubst %.pdf, %.tikz.tex, $@)
+	pdflatex -jobname $(basename $@) -file-line-error $(patsubst %.pdf, %.tikz.tex, $@)
+
+asymptotic.dat:
+	julia asymptotic.jl > asymptotic.dat
+
+install: $(PDFS)
+	cp $^ $(INSTALLDIR)/
+
+clean-aux:
+	rm -f $(addsuffix .aux, $(PROJECTNAME))
+	rm -f $(addsuffix .log, $(PROJECTNAME))
+
+clean-dat:
+	rm -f $(addsuffix .tikz.tex, $(PROJECTNAME))
+	rm -f short-time.dat
+
+clean-tex:
+	rm -f $(PDFS)
+
+clean: clean-dat clean-aux clean-tex
diff --git a/figs/fowler-nordheim.fig/asymptotic.gnuplot b/figs/fowler-nordheim.fig/asymptotic.gnuplot
new file mode 100644
index 0000000..3296380
--- /dev/null
+++ b/figs/fowler-nordheim.fig/asymptotic.gnuplot
@@ -0,0 +1,59 @@
+datafile="asymptotic.dat"
+
+## can also set the following options
+#set title ""
+set ylabel "$|\\psi_{\\mathrm{FN}}|^2$" tc ls 1 #norotate
+set y2label "$J_{\\mathrm{FN}}$" tc ls 2 #norotate
+set xlabel "$x$"
+#
+#set xrange[:]
+#set yrange [:]
+set y2range [0:0.004]
+#
+## start ticks at 0, then every x
+#set xtics 0,x
+#set ytics 0,x
+## puts 4 minor tics between tics (5 intervals, i.e. every 0.01)
+set mxtics 5
+set mytics 5
+set my2tics 5
+
+# default output canvas size: 12.5cm x 8.75cm
+set term lua tikz size 12.5,8.75 standalone
+# run
+## gnuplot gnuplots && gnuplot_tikz out/latext/minimizer.tex
+
+set key off
+
+# 3=1+2 draw bottom and left sides of the box
+#set border 3
+# don't show tics on opposite sides
+set xtics nomirror
+set ytics nomirror tc ls 1
+set y2tics nomirror tc ls 2
+
+# Mathematica colors:
+## 3f3d99 (dark blue)
+## 9c4275 (dark pink)
+## 9a8d3f (dark yellow)
+## 3d9956 (dark green)
+# My colors
+## 4169E1 (pastel blue)
+## DC143C (bright red)
+## 32CD32 (bright green)
+## 4B0082 (deep purple)
+## DAA520 (ochre)
+
+# set linestyle
+set style line 1 linetype rgbcolor "#4169E1" linewidth 3
+set style line 2 linetype rgbcolor "#DC143C" linewidth 3
+set style line 3 linetype rgbcolor "#32CD32" linewidth 3
+set style line 4 linetype rgbcolor "#4B0082" linewidth 3
+set style line 5 linetype rgbcolor "#DAA520" linewidth 3
+
+set pointsize 0.6
+
+set arrow to 0, graph 1 nohead lt 0
+
+plot datafile using 1:2 with lines linestyle 1 ,\
+ datafile using 1:3 with lines linestyle 2 axes x1y2
diff --git a/figs/fowler-nordheim.fig/asymptotic.jl b/figs/fowler-nordheim.fig/asymptotic.jl
new file mode 100644
index 0000000..fd1d492
--- /dev/null
+++ b/figs/fowler-nordheim.fig/asymptotic.jl
@@ -0,0 +1,46 @@
+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
+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 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]),'\n')
+end
diff --git a/figs/libs/Makefile b/figs/libs/Makefile
new file mode 100644
index 0000000..33b81e2
--- /dev/null
+++ b/figs/libs/Makefile
@@ -0,0 +1,28 @@
+PROJECTNAME=$(basename $(basename $(wildcard *.tikz.tex)))
+LIBS=$(notdir $(wildcard libs/*))
+
+PDFS=$(addsuffix .pdf, $(PROJECTNAME))
+
+all: $(PDFS)
+
+$(PDFS): $(LIBS)
+	echo $(LIBS)
+	pdflatex -jobname $(basename $@) -file-line-error $(patsubst %.pdf, %.tikz.tex, $@)
+
+install: $(PDFS)
+	cp $^ $(INSTALLDIR)/
+
+$(LIBS):
+	ln -fs libs/$@ ./
+
+clean-libs:
+	rm -f $(LIBS)
+
+clean-aux:
+	rm -f $(addsuffix .aux, $(PROJECTNAME))
+	rm -f $(addsuffix .log, $(PROJECTNAME))
+
+clean-tex:
+	rm -f $(PDFS)
+
+clean: clean-libs clean-aux clean-tex
diff --git a/figs/potential.fig/Makefile b/figs/potential.fig/Makefile
new file mode 120000
index 0000000..704310e
--- /dev/null
+++ b/figs/potential.fig/Makefile
@@ -0,0 +1 @@
+../libs/Makefile
+\ No newline at end of file
diff --git a/figs/potential.fig/potential.tikz.tex b/figs/potential.fig/potential.tikz.tex
new file mode 100644
index 0000000..8254e5d
--- /dev/null
+++ b/figs/potential.fig/potential.tikz.tex
@@ -0,0 +1,15 @@
+\documentclass{standalone}
+
+\usepackage{tikz}
+
+\begin{document}
+\begin{tikzpicture}
+
+\draw(-3,0)--(3,0);
+\draw(0,-1.5)--(0,3);
+
+\draw[line width=1.5pt](-2.5,0)--(0,0)--(0,2.5)--(1.975,-1.25);
+
+\end{tikzpicture}
+\end{document}
+
diff --git a/figs/potential.fig/potential_square.tikz.tex b/figs/potential.fig/potential_square.tikz.tex
new file mode 100644
index 0000000..0fc8875
--- /dev/null
+++ b/figs/potential.fig/potential_square.tikz.tex
@@ -0,0 +1,15 @@
+\documentclass{standalone}
+
+\usepackage{tikz}
+
+\begin{document}
+\begin{tikzpicture}
+
+\draw(-3,0)--(3,0);
+\draw(0,-1.5)--(0,3);
+
+\draw[line width=1.5pt](-2.5,0)--(0,0)--(0,2.5)--(2.5,2.5);
+
+\end{tikzpicture}
+\end{document}
+