9 files changed, 526 insertions, 0 deletions
diff --git a/figs/plots.fig/FN_base.jl b/figs/plots.fig/FN_base.jl
new file mode 100644
index 0000000..cff6d8e
--- /dev/null
+++ b/figs/plots.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+U-ip)phi=0
+# assume that p has an infinitesimal real part (and adjust the branch cuts appropriately)
+function phi(p,x,E,U)
+  return(airyai_asym(2^(1/3)*exp(-1im*pi/3)*(E^(1/3)*x-E^(-2/3)*(U-1im*p)),pi))
+end
+function dphi(p,x,E,U)
+  return(2^(1/3)*exp(-1im*pi/3)*E^(1/3)*airyaiprime_asym(2^(1/3)*exp(-1im*pi/3)*(E^(1/3)*x-E^(-2/3)*(U-1im*p)),pi))
+end
+function eta(p,x,E,U)
+  return(exp(-1im*pi/3)*airyai_asym(-2^(1/3)*(E^(1/3)*x-E^(-2/3)*(U-1im*p)),pi/2))
+end
+function deta(p,x,E,U)
+  return(-2^(1/3)*exp(-1im*pi/3)*E^(1/3)*airyaiprime_asym(-2^(1/3)*(E^(1/3)*x-E^(-2/3)*(U-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-U)^(3/2) becomes pow(1im*p-U,3/2,-pi/2) because when 1im*p is real negative, its square root should be imaginary positive
+function f(p,x,k0,E,U)
+  T=2im*k0/(1im*k0-sqrt(2*U-k0*k0))
+  R=T-1
+
+  if x>=0
+    C2=-2im*T/(pow(-2im*p,1/2,pi/2)*phi(p,0,E,U)-dphi(p,0,E,U))*((sqrt(2*U-k0*k0)+pow(-2im*p,1/2,pi/2))/(-2im*p+k0*k0)-2im*(2*E)^(-1/3)*pi*quadgk(y -> (pow(-2im*p,1/2,pi/2)*eta(p,0,E,U)-deta(p,0,E,U))*phi(p,y,E,U)*exp(-sqrt(2*U-k0*k0)*y)*exp(sqrt(2)*2im/3*(pow(E^(1/3)*y+E^(-2/3)*(1im*p-U),3/2,-pi/2)-E^(-1)*pow(1im*p-U,3/2,-pi/2))),0,Inf)[1])
+    FT=4*(2*E)^(-1/3)*pi*(quadgk(y -> phi(p,x,E,U)*eta(p,y,E,U)*exp(-sqrt(2*U-k0*k0)*y)*exp(sqrt(2)*2im/3*(pow(E^(1/3)*x+E^(-2/3)*(1im*p-U),3/2,-pi/2)-pow(E^(1/3)*y+E^(-2/3)*(1im*p-U),3/2,-pi/2))),0,x)[1]+quadgk(y -> eta(p,x,E,U)*phi(p,y,E,U)*exp(-sqrt(2*U-k0*k0)*y)*exp(sqrt(2)*2im/3*(pow(E^(1/3)*y+E^(-2/3)*(1im*p-U),3/2,-pi/2)-pow(E^(1/3)*x+E^(-2/3)*(1im*p-U),3/2,-pi/2))),x,Inf)[1])
+    main=C2*phi(p,x,E,U)*exp(sqrt(2)*2im/3*(pow(E^(1/3)*x+E^(-2/3)*(1im*p-U),3/2,-pi/2)-E^(-1)*pow(1im*p-U,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,U)/(p+1im*k0*k0/2)
+    return(main-pole)
+  else
+    C1=-2im*T*((sqrt(2*U-k0*k0)*phi(p,0,E,U)+dphi(p,0,E,U))/(-2im*p+k0*k0)/(pow(-2im*p,1/2,pi/2)*phi(p,0,E,U)-dphi(p,0,E,U))+quadgk(y -> phi(p,y,E,U)/(pow(-2im*p,1/2,pi/2)*phi(p,0,E,U)-dphi(p,0,E,U))*exp(-sqrt(2*U-k0*k0)*y)*exp(sqrt(2)*2im/3*(pow(E^(1/3)*y+E^(-2/3)*(1im*p-U),3/2,-pi/2)-E^(-1)*pow(1im*p-U,3/2,-pi/2))),0,Inf)[1])
+    FI=-2im*exp(1im*k0*x)/(-2im*p+k0*k0)
+    FR=-2im*exp(-1im*k0*x)/(-2im*p+k0*k0)
+    main=C1*exp(pow(-2im*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,U)/(p+1im*k0*k0/2)
+    return(main-pole)
+  end
+end
+# its derivative
+function df(p,x,k0,E,U)
+  T=2im*k0/(1im*k0-sqrt(2*U-k0*k0))
+  R=T-1
+
+  if x>=0
+    C2=-2im*T/(pow(-2im*p,1/2,pi/2)*phi(p,0,E,U)-dphi(p,0,E,U))*((sqrt(2*U-k0*k0)+pow(-2im*p,1/2,pi/2))/(-2im*p+k0*k0)-2im*(2*E)^(-1/3)*pi*quadgk(y -> (pow(-2im*p,1/2,pi/2)*eta(p,0,E,U)-deta(p,0,E,U))*phi(p,y,E,U)*exp(-sqrt(2*U-k0*k0)*y)*exp(sqrt(2)*2im/3*(pow(E^(1/3)*y+E^(-2/3)*(1im*p-U),3/2,-pi/2)-E^(-1)*pow(1im*p-U,3/2,-pi/2))),0,Inf)[1])
+    dFT=4*(2*E)^(-1/3)*pi*(quadgk(y -> dphi(p,x,E,U)*eta(p,y,E,U)*exp(-sqrt(2*U-k0*k0)*y)*exp(sqrt(2)*2im/3*(pow(E^(1/3)*x+E^(-2/3)*(1im*p-U),3/2,-pi/2)-pow(E^(1/3)*y+E^(-2/3)*(1im*p-U),3/2,-pi/2))),0,x)[1]+quadgk(y -> deta(p,x,E,U)*phi(p,y,E,U)*exp(-sqrt(2*U-k0*k0)*y)*exp(sqrt(2)*2im/3*(pow(E^(1/3)*y+E^(-2/3)*(1im*p-U),3/2,-pi/2)-pow(E^(1/3)*x+E^(-2/3)*(1im*p-U),3/2,-pi/2))),x,Inf)[1])
+    main=C2*dphi(p,x,E,U)*exp(sqrt(2)*2im/3*(pow(E^(1/3)*x+E^(-2/3)*(1im*p-U),3/2,-pi/2)-E^(-1)*pow(1im*p-U,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,U)/(p+1im*k0*k0/2)
+    return(main-pole)
+  else
+    C1=-2im*T*((sqrt(2*U-k0*k0)*phi(p,0,E,U)+dphi(p,0,E,U))/(-2im*p+k0*k0)/(pow(-2im*p,1/2,pi/2)*phi(p,0,E,U)-dphi(p,0,E,U))+quadgk(y -> phi(p,y,E,U)/(pow(-2im*p,1/2,pi/2)*phi(p,0,E,U)-dphi(p,0,E,U))*exp(-sqrt(2*U-k0*k0)*y)*exp(sqrt(2)*2im/3*(pow(E^(1/3)*y+E^(-2/3)*(1im*p-U),3/2,-pi/2)-E^(-1)*pow(1im*p-U,3/2,-pi/2))),0,Inf)[1])
+    dFI=2*k0*exp(1im*k0*x)/(-2im*p+k0*k0)
+    dFR=-2*k0*exp(-1im*k0*x)/(-2im*p+k0*k0)
+    main=C1*pow(-2im*p,1/2,pi/2)*exp(pow(-2im*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,U)/(p+1im*k0*k0/2)
+    return(main-pole)
+  end
+end
+
+# psi (returns t,psi(x,t))
+function psi(x,k0,E,U,p_npoints,p_cutoff)
+  fft=fourier_fft(f,x,k0,E,U,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,U)*exp(-1im*k0*k0/2*fft[1][i])
+  end
+  return(fft)
+end
+# its derivative
+function dpsi(x,k0,E,U,p_npoints,p_cutoff)
+  fft=fourier_fft(df,x,k0,E,U,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,U)*exp(-1im*k0*k0/2*fft[1][i])
+  end
+  return(fft)
+end
+
+# compute Fourier transform by sampling and fft
+function fourier_fft(A,x,k0,E,U,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,U)
+    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,U)
+  if x>=0
+    return(1im*phi(-1im*k0*k0/2,x,E,U)*2*k0/(1im*k0*phi(-1im*k0*k0/2,0,E,U)+dphi(-1im*k0*k0/2,0,E,U))*exp(sqrt(2)*2im/3*(pow(E^(1/3)*x+E^(-2/3)*(k0*k0/2-U),3/2,-pi/2)-E^(-1)*pow(k0*k0/2-U,3/2,-pi/2))))
+  else
+    return((1im*k0*phi(-1im*k0*k0/2,0,E,U)-dphi(-1im*k0*k0/2,0,E,U))/(1im*k0*phi(-1im*k0*k0/2,0,E,U)+dphi(-1im*k0*k0/2,0,E,U))*exp(-1im*k0*x)+exp(1im*k0*x))
+  end
+end
+function dpsi_pole(x,k0,E,U)
+  if x>=0
+    return(1im*dphi(-1im*k0*k0/2,x,E,U)*2*k0/(1im*k0*phi(-1im*k0*k0/2,0,E,U)+dphi(-1im*k0*k0/2,0,E,U))*exp(sqrt(2)*2im/3*(pow(E^(1/3)*x+E^(-2/3)*(k0*k0/2-U),3/2,-pi/2)-E^(-1)*pow(k0*k0/2-U,3/2,-pi/2))))
+  else
+    return(-1im*k0*(1im*k0*phi(-1im*k0*k0/2,0,E,U)-dphi(-1im*k0*k0/2,0,E,U))/(1im*k0*phi(-1im*k0*k0/2,0,E,U)+dphi(-1im*k0*k0/2,0,E,U))*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,U,p_npoints,p_cutoff)
+  ps=psi(x,k0,E,U,p_npoints,p_cutoff)
+  dps=dpsi(x,k0,E,U,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/plots.fig/Makefile b/figs/plots.fig/Makefile
new file mode 100644
index 0000000..d9fdf82
--- /dev/null
+++ b/figs/plots.fig/Makefile
@@ -0,0 +1,36 @@
+PROJECTNAME=currents-4 currents-8 density-4 density-8 integrated_current-4 integrated_current-8
+
+PDFS=$(addsuffix .pdf, $(PROJECTNAME))
+TEXS=$(addsuffix .tikz.tex, $(PROJECTNAME))
+
+all: $(PDFS)
+
+currents-4.pdf density-4.pdf integrated_current-4.pdf: current_density-4.dat
+	gnuplot $(patsubst %.pdf, %.gnuplot, $@) > $(patsubst %.pdf, %.tikz.tex, $@)
+	pdflatex -jobname $(basename $@) -file-line-error $(patsubst %.pdf, %.tikz.tex, $@)
+
+currents-8.pdf density-8.pdf integrated_current-8.pdf: current_density-8.dat
+	gnuplot $(patsubst %.pdf, %.gnuplot, $@) > $(patsubst %.pdf, %.tikz.tex, $@)
+	pdflatex -jobname $(basename $@) -file-line-error $(patsubst %.pdf, %.tikz.tex, $@)
+
+current_density-4.dat:
+	julia current_density.jl 4 > $@
+current_density-8.dat:
+	julia current_density.jl 8 > $@
+
+install: $(PDFS)
+	cp $^ $(INSTALLDIR)/
+
+clean-aux:
+	rm -f $(addsuffix .tikz.tex, $(PROJECTNAME))
+	rm -f $(addsuffix .aux, $(PROJECTNAME))
+	rm -f $(addsuffix .log, $(PROJECTNAME))
+
+clean-dat:
+	rm -f current_density-4.dat
+	rm -f current_density-8.dat
+
+clean-tex:
+	rm -f $(PDFS)
+
+clean: clean-dat clean-aux clean-tex
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
diff --git a/figs/plots.fig/currents-4.gnuplot b/figs/plots.fig/currents-4.gnuplot
new file mode 100644
index 0000000..1e9dba9
--- /dev/null
+++ b/figs/plots.fig/currents-4.gnuplot
@@ -0,0 +1,43 @@
+## can also set the following options
+set title "$E=4\\ \\mathrm{V}\\cdot\\mathrm{nm}^{-1}$"
+set ylabel "$\\displaystyle\\frac j{2k}$" norotate
+set xlabel "$t$ (fs)"
+#
+## 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
+
+# default output canvas size: 12.5cm x 8.75cm
+set term lua tikz size 12.5,8.75 standalone
+
+set key  spacing 1.5
+
+# 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
+
+# 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 xrange [:12]
+
+plot "current_density-4.dat" using 1:3 every ::1 with lines linestyle 1 title "$x=x_0$",\
+     "current_density-4.dat" using 5:7 every ::1 with lines linestyle 2 title "$x=10x_0$"
diff --git a/figs/plots.fig/currents-8.gnuplot b/figs/plots.fig/currents-8.gnuplot
new file mode 100644
index 0000000..62448f7
--- /dev/null
+++ b/figs/plots.fig/currents-8.gnuplot
@@ -0,0 +1,43 @@
+## can also set the following options
+set title "$E=8\\ \\mathrm{V}\\cdot\\mathrm{nm}^{-1}$"
+set ylabel "$\\displaystyle\\frac j{2k}$" norotate
+set xlabel "$t$ (fs)"
+#
+## 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
+
+# default output canvas size: 12.5cm x 8.75cm
+set term lua tikz size 12.5,8.75 standalone
+
+set key spacing 1.5
+
+# 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
+
+# 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 xrange [:12]
+
+plot "current_density-8.dat" using 1:3 every ::1 with lines linestyle 1 title "$x=x_0$",\
+     "current_density-8.dat" using 5:7 every ::1 with lines linestyle 2 title "$x=10x_0$"
diff --git a/figs/plots.fig/density-4.gnuplot b/figs/plots.fig/density-4.gnuplot
new file mode 100644
index 0000000..8139a7a
--- /dev/null
+++ b/figs/plots.fig/density-4.gnuplot
@@ -0,0 +1,44 @@
+## can also set the following options
+set title "$E=4\\ \\mathrm{V}\\cdot\\mathrm{nm}^{-1}$"
+set ylabel "$|\\psi|^2$" norotate
+set xlabel "$t$ (fs)"
+#
+## 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
+
+# default output canvas size: 12.5cm x 8.75cm
+set term lua tikz size 12.5,8.75 standalone
+
+set key spacing 1.5
+
+# 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
+
+# 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 xrange [:12]
+
+plot \
+     "current_density-4.dat" using 1:2 every ::1 with lines linestyle 1 title "$x=x_0$",\
+     "current_density-4.dat" using 5:6 every ::1 with lines linestyle 2 title "$x=10x_0$"
diff --git a/figs/plots.fig/density-8.gnuplot b/figs/plots.fig/density-8.gnuplot
new file mode 100644
index 0000000..f107ff6
--- /dev/null
+++ b/figs/plots.fig/density-8.gnuplot
@@ -0,0 +1,44 @@
+## can also set the following options
+set title "$E=8\\ \\mathrm{V}\\cdot\\mathrm{nm}^{-1}$"
+set ylabel "$|\\psi|^2$" norotate
+set xlabel "$t$ (fs)"
+#
+## 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
+
+# default output canvas size: 12.5cm x 8.75cm
+set term lua tikz size 12.5,8.75 standalone
+
+set key spacing 1.5
+
+# 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
+
+# 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 xrange [:12]
+
+plot \
+     "current_density-8.dat" using 1:2 every ::1 with lines linestyle 1 title "$x=x_0$",\
+     "current_density-8.dat" using 5:6 every ::1 with lines linestyle 2 title "$x=10x_0$"
diff --git a/figs/plots.fig/integrated_current-4.gnuplot b/figs/plots.fig/integrated_current-4.gnuplot
new file mode 100644
index 0000000..7e9dd6b
--- /dev/null
+++ b/figs/plots.fig/integrated_current-4.gnuplot
@@ -0,0 +1,43 @@
+## can also set the following options
+set title "$E=4\\ \\mathrm{V}\\cdot\\mathrm{nm}^{-1}$"
+set ylabel "$\\displaystyle\\frac J{k^2}$" norotate
+set xlabel "$t$ (fs)"
+#
+## 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
+
+# default output canvas size: 12.5cm x 8.75cm
+set term lua tikz size 12.5,8.75 standalone
+
+set key  spacing 1.5
+
+# 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
+
+# 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 xrange [:12]
+
+plot "current_density-4.dat" using 1:4 every ::1 with lines linestyle 1 title "$x=x_0$",\
+     "current_density-4.dat" using 5:8 every ::1 with lines linestyle 2 title "$x=10x_0$"
diff --git a/figs/plots.fig/integrated_current-8.gnuplot b/figs/plots.fig/integrated_current-8.gnuplot
new file mode 100644
index 0000000..3d5d501
--- /dev/null
+++ b/figs/plots.fig/integrated_current-8.gnuplot
@@ -0,0 +1,43 @@
+## can also set the following options
+set title "$E=4\\ \\mathrm{V}\\cdot\\mathrm{nm}^{-1}$"
+set ylabel "$\\displaystyle\\frac J{k^2}$" norotate
+set xlabel "$t$ (fs)"
+#
+## 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
+
+# default output canvas size: 12.5cm x 8.75cm
+set term lua tikz size 12.5,8.75 standalone
+
+set key  spacing 1.5
+
+# 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
+
+# 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 xrange [:12]
+
+plot "current_density-8.dat" using 1:4 every ::1 with lines linestyle 1 title "$x=x_0$",\
+     "current_density-8.dat" using 5:8 every ::1 with lines linestyle 2 title "$x=10x_0$"