As presented at tthe 121st Statistical Mechanics Meeting on 2019-05-14HEADv1.0master

13 files changed, 1035 insertions, 0 deletions

diff --git a/Jauslin_SMM121_2019.tex b/Jauslin_SMM121_2019.tex
new file mode 100644
index 0000000..974976d
--- /dev/null
+++ b/Jauslin_SMM121_2019.tex
@@ -0,0 +1,148 @@
+\documentclass{ian-presentation}
+
+\usepackage[hidelinks]{hyperref}
+\usepackage{graphicx}
+\usepackage{array}
+
+\begin{document}
+\pagestyle{empty}
+\hbox{}\vfil
+\bf\Large
+\hfil Exact solution of the time dependent Schrödinger\par
+\hfil equation for photoemission from a metal surface
+\vfil
+\large
+\hfil Ian Jauslin
+\normalsize
+\vfil
+\hfil\rm joint with {\bf Ovidiu Costin}, {\bf Rodica Costin}, and {\bf Joel L. Lebowitz}\par
+\vfil
+\hfill{\tt \href{http://ian.jauslin.org}{http://ian.jauslin.org}}
+\eject
+
+\setcounter{page}1
+\pagestyle{plain}
+
+\title{Electron emission}
+$$
+  V(x)=U\Theta(x)
+  ,\quad
+  E_{\mathrm F}=k_{\mathrm F}^2<U
+$$
+\hfil\includegraphics[height=5cm]{potential_square.pdf}
+\vfill
+\eject
+
+\title{Thermal emission}
+$$
+  V(x)=U\Theta(x)
+  ,\quad
+  k^2>U
+$$
+\hfil\includegraphics[height=5cm]{potential_square_thermal.pdf}
+\vfill
+\eject
+
+\title{Field emission}
+$$
+  V(x)=\Theta(x)(U-Ex)
+$$
+\hfil\includegraphics[height=5cm]{potential.pdf}
+\vfill
+\eject
+
+\title{Photoemission}
+$$
+  V_t(x)=\Theta(x)(U-E_tx)
+  ,\quad
+  E_t=2\epsilon\omega\cos(\omega t)
+$$
+\hfil\includegraphics[height=5cm]{potential_square_photonic.pdf}
+\vfill
+\eject
+
+\title{Photoemission}
+\begin{itemize}
+  \item Photoelectric effect: early observations: \href{https://doi.org/10.1002/andp.18872670827}{[Hertz, 1887]}, \href{https://doi.org/10.1002/andp.18882690206}{[Hallwachs, 1888]}, \href{https://doi.org/10.1002/andp.19003070611}{[Lenard, 1900]}.
+  \item When a metal is irradiated with ultra-violet light, electrons are ionized, with kinetic energies at integer multiples of $\hbar\omega$.
+  \item \href{https://doi.org/10.1002/andp.19053220607}{[Einstein, 1905]}: interpretation: quanta of light ({\it photons}) of energy $\hbar\omega$ are absorbed by the electrons, whose kinetic energy is raised by $n\hbar\omega$, and can escape the metal.
+\end{itemize}
+\vfill
+\eject
+
+\title{Photoemission}
+\vskip-10pt
+\begin{itemize}
+  \item Time dependent potential:
+  $$
+    V_t(x)=\Theta(x)(U-2\epsilon\omega\cos(\omega t)x)
+  $$
+\vskip-10pt
+  \item Schr\"odinger equation
+  $$
+    i\partial_t\psi(x,t)=-\Delta\psi(x,t)+V_t(x)\psi(x,t)
+  $$
+\vskip-10pt
+  \item Magnetic gauge:
+  $$
+    \Psi(x,t)
+    :=\psi(x,t)e^{-ix\Theta(x)A(t)}
+    ,\quad
+    A(t):=\int_0^t ds\ 2\epsilon\omega\cos(\omega s)
+  =
+  2\epsilon\sin(\omega t)
+  $$
+  satisfies
+  $$
+    i\partial_t\Psi(x,t)=\left((-i\nabla+\Theta(x)A(t))^2+\Theta(x)U\right)\Psi(x,t)
+  $$
+\end{itemize}
+\vfill
+\eject
+
+\title{Periodic solution}
+\vskip-10pt
+\begin{itemize}
+  \item \href{https://doi.org/10.1103/PhysRevA.72.023412}{[Faisal, Kami\'nski, Saczuk, 2005]}
+  $$
+    \Psi_{\mathrm{FKS}}(x,t)=\left\{\begin{array}{ll}
+      e^{ikx}\exp\left(-ik^2t\right)+\Psi_R(x,t)&\mathrm{\ if\ }x<0\\
+      \Psi_T(x,t)&\mathrm{\ if\ }x>0
+    \end{array}\right.
+  $$
+  $$
+    \Psi_R(x,t)=\sum_{M\in\mathbb Z}R_Me^{-iq_Mx}\exp\left(-iq_M^2t\right)
+    ,\quad
+    q_M=\sqrt{k^2+M\omega}
+  $$
+  $$
+    \Psi_T(x,t)=\sum_{M\in\mathbb Z}T_Me^{ip_Mx}\exp\left(-iUt-i\int_0^td\tau\ (p_M-A(\tau))^2\right)
+  $$
+  $$
+    p_M=\sqrt{k^2-U+M\omega-2\epsilon^2}
+  $$
+  \item $\Psi(x,t)$, $(-i\nabla+\Theta(x)A(t))\Psi(x,t)$ are continuous.
+\end{itemize}
+\vfill
+\eject
+
+\title{Initial value problem}
+\begin{itemize}
+  \item Initial condition:
+  $$
+    \Psi(x,0)=
+    \left\{ \begin{array}{l@{\ }l}
+      e^{ikx}+R_0e^{-ikx} & x<0\\
+      T_0 e^{-\sqrt{U-k^2}x} & x>0
+    \end{array}\right.  
+  $$
+  $R_0$ and $T_0$ ensure that $\Psi$ and $\partial\Psi$ are continuous.
+  \item In progress: $\Psi(x,t)$ behaves asymptotically like $\Psi_{\mathrm{FKS}}$:
+  $$
+    \psi(x,t)
+    =\psi_{\mathrm{FKS}}(x,t)+\left(\frac{t}{\tau_{\mathrm{FKS}}(x)}\right)^{-\frac32}+O(t^{-\frac52})
+    .
+  $$
+\end{itemize}
+
+\end{document}
diff --git a/Makefile b/Makefile
new file mode 100644
index 0000000..4ed7b60
--- /dev/null
+++ b/Makefile
@@ -0,0 +1,50 @@
+PROJECTNAME=$(basename $(wildcard *.tex))
+LIBS=$(notdir $(wildcard libs/*))
+FIGS=$(notdir $(wildcard figs/*.fig))
+
+PDFS=$(addsuffix .pdf, $(PROJECTNAME))
+SYNCTEXS=$(addsuffix .synctex.gz, $(PROJECTNAME))
+
+all: $(PROJECTNAME)
+
+$(PROJECTNAME): $(LIBS) $(FIGS) images
+	pdflatex -file-line-error $@.tex
+	pdflatex -synctex=1 $@.tex
+
+$(SYNCTEXS): $(LIBS) $(FIGS) images
+	pdflatex -synctex=1 $(patsubst %.synctex.gz, %.tex, $@)
+
+
+$(LIBS):
+	ln -fs libs/$@ ./
+
+
+$(FIGS):
+	make -C figs/$@
+	for pdf in $$(find figs/$@/ -name '*.pdf'); do ln -fs "$$pdf" ./ ; done
+
+images:
+	for jpg in figs/*.jpg; do ln -fs "$$jpg" ./ ; done
+	for png in figs/*.png; do ln -fs "$$png" ./ ; done
+
+
+clean-aux: clean-figs-aux
+	rm -f $(addsuffix .aux, $(PROJECTNAME))
+	rm -f $(addsuffix .log, $(PROJECTNAME))
+	rm -f $(addsuffix .out, $(PROJECTNAME))
+
+clean-libs:
+	rm -f $(LIBS)
+
+clean-figs:
+	$(foreach fig,$(addprefix figs/, $(FIGS)), make -C $(fig) clean; )
+	rm -f $(notdir $(wildcard figs/*.fig/*.pdf))
+
+clean-figs-aux:
+	$(foreach fig,$(addprefix figs/, $(FIGS)), make -C $(fig) clean-aux; )
+
+
+clean-tex:
+	rm -f $(PDFS) $(SYNCTEXS)
+
+clean: clean-aux clean-tex clean-libs clean-figs
diff --git a/README b/README
new file mode 100644
index 0000000..1be07e7
--- /dev/null
+++ b/README
@@ -0,0 +1,42 @@
+This directory contains the source files to typeset the presentation, and
+generate the figures. This can be accomplished by running
+  make
+
+An animation was prepared for this talk, which shows the time-evolution of the
+density and current. The source files are located in the 'figs/anim_laser'
+directory.  To run the animation, run
+  cd figs/anim_laser
+  make run
+This can take a while the first time it is run.
+
+This document uses a custom class file, located in the 'libs' directory, which
+defines a number of commands.
+
+
+* Dependencies:
+
+  pdflatex
+  TeXlive packages:
+    amsfonts
+    array
+    graphics
+    hyperref
+    latex
+    pgf
+    standalone
+  GNU make
+  julia
+  python3
+
+
+* Files:
+
+  Jauslin_SMM121_2019.tex:
+    main LaTeX file
+
+  libs:
+    custom LaTeX class file
+
+  figs:
+    source code for the figures
+
diff --git a/figs/anim_laser/Makefile b/figs/anim_laser/Makefile
new file mode 100644
index 0000000..9b66e37
--- /dev/null
+++ b/figs/anim_laser/Makefile
@@ -0,0 +1,12 @@
+PROJECTNAME=animate
+
+all: animate.dat
+
+run: animate.dat
+	python3 animate.py animate.dat &
+
+animate.dat:
+	julia laser.jl > animate.dat
+
+clean:
+	rm -f animate.dat
diff --git a/figs/anim_laser/animate.py b/figs/anim_laser/animate.py
new file mode 100644
index 0000000..1e6995a
--- /dev/null
+++ b/figs/anim_laser/animate.py
@@ -0,0 +1,103 @@
+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:
+    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')
+cn_rho, = axr.plot([],[],color='#FF0000')
+
+pl.subplot(212)
+axJ=fig.gca()
+asym_J, = axJ.plot([],[],linewidth=3.5,color='#00FF00')
+cn_J, = axJ.plot([],[],color='#FF0000')
+
+# 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]>maxyr:
+            maxyr=frame[i][3]
+        if frame[i][4]>maxyJ:
+            maxyJ=frame[i][4]
+        if frame[i][5]>maxyJ:
+            maxyJ=frame[i][5]
+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]<minyr:
+            minyr=frame[i][3]
+        if frame[i][4]<minyJ:
+            minyJ=frame[i][4]
+        if frame[i][5]<minyJ:
+            minyJ=frame[i][5]
+
+# 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")
+    axJ.hlines(0,xmin,xmax,linestyles="dotted")
+def update(frame):
+    axr.set_title("t=% .3f fs" % (frame[0][0]))
+    xdata=[]
+    ydata=[]
+    asym_data=[]
+    cn_data=[]
+    for i in range(len(frame)):
+        xdata.append(frame[i][1])
+        asym_data.append(frame[i][2])
+        cn_data.append(frame[i][3])
+    asym_rho.set_data(xdata,asym_data)
+    cn_rho.set_data(xdata,cn_data)
+
+    xdata=[]
+    ydata=[]
+    asym_data=[]
+    cn_data=[]
+    for i in range(len(frame)):
+        xdata.append(frame[i][1])
+        asym_data.append(frame[i][4])
+        cn_data.append(frame[i][5])
+    asym_J.set_data(xdata,asym_data)
+    cn_J.set_data(xdata,cn_data)
+anim = animation.FuncAnimation(fig, update, frames=frames, blit=False, interval=100, repeat=True, init_func=init_plot)
+#anim.save('laser_schrodinger.m4v', fps=10, extra_args=['-vcodec', 'libx264'])
+pl.show()
diff --git a/figs/anim_laser/laser.jl b/figs/anim_laser/laser.jl
new file mode 100644
index 0000000..78cabff
--- /dev/null
+++ b/figs/anim_laser/laser.jl
@@ -0,0 +1,391 @@
+using FFTW
+using Printf
+using SpecialFunctions
+using LinearAlgebra
+
+# numerical values
+hbar=6.58e-16 # eV.s
+m=9.11e-31 # kg
+V_=6.6 # eV
+I=1e17 # W/m^2
+EF=3.3 # eV
+omega_=1.5498*4 # eV
+
+# dimensionless quantities
+V=1
+k=sqrt(EF/V_)
+epsilon=sqrt(hbar^3*I/(8*m*V_*omega_^2))
+omega=omega_/V_
+
+#epsilon=8*epsilon
+#omega=4*omega
+
+nu=60
+ell=20
+
+# for Crank-Nicolson
+xmin_cn=-100
+xmax_cn=100
+nx_cn=4000
+tmax=10
+# mu=dt/dx^2
+mu=0.1
+nt=floor(Int,tmax/(mu*((xmax_cn-xmin_cn)/nx_cn)^2))
+# only save up to memt time steps
+memt=10000
+
+
+# asymptotic value of \varphi_t(x)
+function phi_asym(x,T,gkz,k,epsilon,omega,V,nu,ell)
+  # compute residues
+  res=residue(x,T,gkz,k,epsilon,omega,V,nu)
+  times=compute_times(omega,nu,ell)
+  N=(2*nu+1)*(2*ell+1)
+  fun=zeros(Complex{Float64},N)
+  for m in 1:N
+    for n in -nu:nu
+      fun[m]+=-1im*res[n+nu+1]*exp(-1im*(k*k+n*omega)*times[m])
+    end
+  end
+  return fun
+end
+# asymptotic value of \partial\varphi_t(x)
+function dphi_asym(x,T,gkz,k,epsilon,omega,V,nu,ell)
+  # compute residues
+  res=dresidue(x,T,gkz,k,epsilon,omega,V,nu)
+  times=compute_times(omega,nu,ell)
+  N=(2*nu+1)*(2*ell+1)
+  fun=zeros(Complex{Float64},N)
+  for m in 1:N
+    for n in -nu:nu
+      fun[m]+=-1im*res[n+nu+1]*exp(-1im*(k*k+n*omega)*times[m])
+    end
+  end
+  return fun
+end
+
+# return times
+function compute_times(omega,nu,ell)
+  N=(2*nu+1)*(2*ell+1)
+  times=zeros(N)
+  for j in 1:N
+    times[j]=pi/(omega*(nu+0.5))*(j-1)
+  end
+  return times
+end
+
+# \kappa_m^{(\sigma)}
+function K(m,sigma,V,k,epsilon,omega)
+  return pow(V+2*epsilon*epsilon-k*k-sigma-m*omega,0.5,pi/2)
+end
+# mathfrak q_k
+function Q(k,V)
+  return sqrt(V-k*k)
+end
+
+# T_0
+function T0(k,V)
+  return 2im*k/(1im*k-sqrt(V-k*k))
+end
+
+# compute g
+function g(kappa,epsilon,omega,N)
+  B=zeros(Complex{Float64},N)
+  # prepare vector of B's
+  for i in 1:N
+    t=2*pi/(omega*N)*(i-1)
+    B[i]=exp(1im*epsilon*epsilon/omega*sin(2*omega*t)-kappa*4*epsilon/omega*cos(omega*t))
+  end
+  # take the fft
+  fft!(B)
+  # correct by a factor 1/N
+  for i in 1:N
+    B[i]=B[i]/N
+  end
+  return(B)
+end
+
+# index of m-n
+function mmn(n,nu)
+  return (n>=0 ? n+1 : n+4*nu+2)
+end
+
+# compute residues at -ik^2-in\omega
+function residue(x,T,gkz,k,epsilon,omega,V,nu)
+  # compute residue
+  out=zeros(Complex{Float64},2*nu+1)
+  if x>=0
+    for n in -nu:nu
+      for m in -nu:nu
+	out[n+nu+1]+=1im*gkz[m+nu+1][mmn(n-m,nu)]*exp(-K(m,0,V,k,epsilon,omega)*x)*T[m+nu+1]
+      end
+    end
+  else
+    for n in -nu:nu
+      out[n+nu+1]=(n==0 ? 1im*(exp(1im*k*x)-exp(-1im*k*x)) : 0)
+      for m in -nu:nu
+	out[n+nu+1]+=1im*T[m+nu+1]*gkz[m+nu+1][mmn(n-m,nu)]*exp(-1im*x*pow(k*k+n*omega,0.5,-pi/2))
+      end
+    end
+  end
+
+  return out
+end
+# compute \partial residues at -ik^2-in\omega
+function dresidue(x,T,gkz,k,epsilon,omega,V,nu)
+  # compute residue
+  out=zeros(Complex{Float64},2*nu+1)
+  if x>=0
+    for n in -nu:nu
+      for m in -nu:nu
+	out[n+nu+1]+=-K(m,0,V,k,epsilon,omega)*1im*gkz[m+nu+1][mmn(n-m,nu)]*exp(-K(m,0,V,k,epsilon,omega)*x)*T[m+nu+1]
+      end
+    end
+  else
+    for n in -nu:nu
+      out[n+nu+1]=(n==0 ? -k*(exp(1im*k*x)+exp(-1im*k*x)) : 0)
+      for m in -nu:nu
+	out[n+nu+1]+=pow(k*k+n*omega,0.5,-pi/2)*T[m+nu+1]*gkz[m+nu+1][mmn(n-m,nu)]*exp(-1im*x*pow(k*k+n*omega,0.5,-pi/2))
+      end
+    end
+  end
+
+  return out
+end
+
+
+# matching condition for residue
+function residue_matching(k,epsilon,omega,V,nu)
+  # preliminary: g^{(\kappa_m^{(0)})}
+  # g^{(\kappa_m^{(0)})}_n=gkz[m+nu+1][n+1] if n>=0 and gkz[m+nu+1][n+4*nu+2] if n<0
+  gkz=Array{Array{Complex{Float64}},1}(undef,2*nu+1)
+  for m in -nu:nu
+    gkz[m+nu+1]=g(K(m,0,V,k,epsilon,omega),epsilon,omega,4*nu+1)
+  end
+
+  # solve matching condition
+  G=zeros(Complex{Float64},2*nu+1,2*nu+1)
+  v=zeros(Complex{Float64},2*nu+1)
+  for n in -nu:nu
+    for m in -nu:nu
+      G[n+nu+1,m+nu+1]=((K(m,0,V,k,epsilon,omega)-1im*pow(k*k+n*omega,0.5,-pi/2))*gkz[m+nu+1][mmn(n-m,nu)]-epsilon*(gkz[m+nu+1][mmn(n-m+1,nu)]-gkz[m+nu+1][mmn(n-m-1,nu)]))
+    end
+  end
+  v[nu+1]=-2im*k
+  
+  return (G\v,gkz)
+end
+
+
+# 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
+
+# current
+function J(phi,dphi,x,epsilon,omega,t)
+  if x>=0
+    return(2*imag(conj(phi)*(dphi+2im*epsilon*sin(omega*t)*phi)))
+  else
+    return(2*imag(conj(phi)*dphi))
+  end
+end
+
+# compute \psi_t(x) using Crank-Nicolson
+function psi_cn(V,omega,epsilon,k,xs,ts,memt)
+  nx=length(xs)
+  nt=length(ts)
+
+  # length of output vector
+  len=min(nt,memt)
+  # how often to write to output
+  freq=ceil(Int,nt/memt)
+
+  # initialize
+  psis=zeros(Complex{Float64},nx,len)
+  times=zeros(Float64,len)
+
+  # previous time step
+  psi_prev=zeros(Complex{Float64},nx)
+
+  # init
+  for i in 1:nx
+    times[1]=0
+    psis[i,1]=(xs[i]<0 ? exp(1im*k*xs[i])+(T0(k,V)-1)*exp(-1im*k*xs[i]) : T0(k,V)*exp(-Q(k,V)*xs[i]))
+    psi_prev[i]=psis[i,1]
+  end
+  # the boundary condition is that psi at the boundary is constant (nothing needs to be done for this)
+
+  # matrix structure of the Crank-Nicolson algorithm
+  # diagonals of M
+  M0=zeros(Complex{Float64},nx-2)
+  Mp=zeros(Complex{Float64},nx-3)
+  Mm=zeros(Complex{Float64},nx-3)
+  v=zeros(Complex{Float64},nx-2)
+
+  for j in 1:nt-1
+    # print progress
+    progress(j,nt-1,1000)
+
+    # t_j
+    t=ts[j]
+    # t_{j+1}
+    tp=ts[j+1]
+    dt=tp-t
+    for mu in 2:nx-1
+      # x_mu
+      x=xs[mu]
+      dx2=(xs[mu+1]-xs[mu])*(xs[mu]-xs[mu-1])
+      # tridiagonal matrix
+      M0[mu-1]=1im/dt-1/dx2-(x<0 ? 0 : (V-2*epsilon*omega*cos(omega*tp)*x)/2)
+      if mu<nx-1
+	Mp[mu-1]=1/(2*dx2)
+      end
+      if mu>2
+	Mm[mu-2]=1/(2*dx2)
+      end
+      # right side
+      v[mu-1]=1im*psi_prev[mu]/dt-(psi_prev[mu+1]+psi_prev[mu-1]-2*psi_prev[mu])/(2*dx2)+(x<0 ? 0 : (V-2*epsilon*omega*cos(omega*t)*x)*psi_prev[mu]/2)
+      # correct for boundary conditions
+      # assumes the boundary condition is constant in time!
+      if mu==2
+	v[mu-1]-=psi_prev[1]/(2*dx2)
+      end
+      if mu==nx-1
+	v[mu-1]-=psi_prev[nx]/(2*dx2)
+      end
+    end
+
+    M=Tridiagonal(Mm,M0,Mp)
+
+    # copy to psi_prev
+    psi_prev[2:nx-1]=M\v
+
+    # write to output
+    if j%freq==0
+      psis[:,Int(j/freq)+1]=psi_prev
+      times[Int(j/freq)+1]=tp
+    end
+  end
+
+  return (times,psis)
+end
+
+# returns the times in b that are closest to those in a
+# assumes the vectors are ordered
+function nearest_times(a,b)
+  out=zeros(Int,length(a))
+  pointer=1
+  for i in 1:length(a)
+    if pointer==length(b)
+      out[i]=length(b)
+    end
+    for j in pointer:length(b)-1
+      if b[j+1]>a[i]
+	out[i]=(abs(b[j+1]-a[i])<abs(b[j]-a[i]) ? j+1 : j)
+	pointer=out[i]
+	break
+      end
+      # all remaining points are beyond the range
+      if j==length(b)-1
+	out[i]=length(b)
+	pointer=out[i]
+      end
+    end
+  end
+  return out
+end
+
+# space derivative using Crank Nicolson
+function deriv_cn(phi,i,xs)
+  if i==1
+    return (phi[i+1]-phi[i])/(xs[i+1]-xs[i])
+  end
+  if i==length(phi)
+    return (phi[i]-phi[i-1])/(xs[i]-xs[i-1])
+  end
+  return (phi[i+1]-phi[i])/(xs[i+1]-xs[i])/2 + (phi[i]-phi[i-1])/(xs[i]-xs[i-1])/2
+end
+
+################################################
+
+
+# print progress
+function progress(j,tot,freq)
+  if (j-1)%ceil(Int,tot/freq)==0
+    if j>1
+      @printf(stderr,"\r")
+    end
+    @printf(stderr,"%d/%d",j,tot)
+  end
+  if j==tot
+    @printf(stderr,"%d/%d\n",j,tot)
+  end
+end
+
+# print animation data using Crank Nicolson and compare to Laplace
+function print_anim_cn_cmp()
+  @printf(stderr,"epsilon=% .8e omega=% .8e k=% .8e V=% .8e z=% .8e T=% .8e gamma=% .8e\n",epsilon,omega,k,V,4*sqrt(nu)*epsilon/sqrt(omega),2*pi/(omega_/hbar)*1e15,sqrt((V-k^2)/(4*epsilon^2)))
+  xmin=-10
+  xmax=10
+  nx=200
+
+  # array of times for the asymptote
+  times=compute_times(omega,nu,ell)
+
+  # positions at which to compute phi
+  xs=Array(0:nx_cn)*(xmax_cn-xmin_cn)/nx_cn.+xmin_cn
+  # times at which to compute phi
+  ts=Array(0:nt-1)*tmax/nt
+
+  # compute matching condition
+  (T,gkz)=residue_matching(k,epsilon,omega,V,nu)
+
+  # compute wave function
+  (t_psi,psi)=psi_cn(V,omega,epsilon,k,xs,ts,memt)
+
+  # nearest times using Crank Nicolson
+  indices_t=nearest_times(times,t_psi)
+  indices_x=nearest_times(Array(0:nx-1)*(xmax-xmin)/nx.+xmin,xs)
+
+  ps_cn=Array{Array{Complex{Float64},1}}(undef,nx)
+  ps_asym=Array{Array{Complex{Float64},1}}(undef,nx)
+  dps_cn=Array{Array{Complex{Float64},1}}(undef,nx)
+  dps_asym=Array{Array{Complex{Float64},1}}(undef,nx)
+  for i in 1:nx
+    progress(i,nx,1000)
+    x=xmin+(xmax-xmin)*(i-1)/nx
+    ps_cn[i]=Array{Complex{Float64}}(undef,length(times))
+    dps_cn[i]=Array{Complex{Float64}}(undef,length(times))
+    for j in 1:length(times)
+      # phase change in magnetic gauge
+      phase=(x<0 ? 1 : exp(-2im*epsilon*sin(omega*times[j])*x))
+      ps_cn[i][j]=phase*psi[indices_x[i],indices_t[j]]
+      dps_cn[i][j]=phase*deriv_cn(psi[:,indices_t[j]],indices_x[i],xs)-(xs[indices_x[i]]<0 ? 0 : 2im*epsilon*sin(omega*times[j])*psi[indices_x[i],indices_t[j]])
+    end
+    ps_asym[i]=phi_asym(x,T,gkz,k,epsilon,omega,V,nu,ell)
+    dps_asym[i]=dphi_asym(x,T,gkz,k,epsilon,omega,V,nu,ell)
+  end
+
+  # print values at each time
+  for j in 1:length(times)
+    if times[j]>tmax
+      break
+    end
+    for i in 1:nx
+      x=(xmin+(xmax-xmin)*i/nx)
+      # remove point at 0
+      if i!=Int(nx/2)
+	@printf("% .8e % .8e % .8e % .8e % .8e % 8e\n",times[j]*hbar/V_*1e15,x*sqrt(hbar^2*1.6e-19/(2*m*V_))*1e9,abs(ps_asym[i][j])^2,abs(ps_cn[i][j])^2,J(ps_asym[i][j],dps_asym[i][j],x,epsilon,omega,times[j]),J(ps_cn[i][j],dps_cn[i][j],x,epsilon,omega,times[j]))
+      end
+    end
+    print('\n')
+  end
+end
+
+
+print_anim_cn_cmp()
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..ac1ba36
--- /dev/null
+++ b/figs/potential.fig/potential.tikz.tex
@@ -0,0 +1,17 @@
+\documentclass{standalone}
+
+\usepackage{tikz}
+
+\begin{document}
+\begin{tikzpicture}
+
+\draw(-3,0)--(3,0);
+\draw(0,-1.5)--(0,3);
+
+\draw[color=gray,line width=1pt, densely dotted](-2.5,1.25)--++(2.5,0);
+
+\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..bd13467
--- /dev/null
+++ b/figs/potential.fig/potential_square.tikz.tex
@@ -0,0 +1,17 @@
+\documentclass{standalone}
+
+\usepackage{tikz}
+
+\begin{document}
+\begin{tikzpicture}
+
+\draw(-3,0)--(3,0);
+\draw(0,-1.5)--(0,3);
+
+\draw[color=gray,line width=1pt, densely dotted](-2.5,1.25)--++(2.5,0);
+
+\draw[line width=1.5pt](-2.5,0)--(0,0)--(0,2.5)--(2.5,2.5);
+
+\end{tikzpicture}
+\end{document}
+
diff --git a/figs/potential.fig/potential_square_photonic.tikz.tex b/figs/potential.fig/potential_square_photonic.tikz.tex
new file mode 100644
index 0000000..8a35fb4
--- /dev/null
+++ b/figs/potential.fig/potential_square_photonic.tikz.tex
@@ -0,0 +1,22 @@
+\documentclass{standalone}
+
+\usepackage{tikz}
+
+\begin{document}
+\begin{tikzpicture}
+
+\draw(-3,0)--(3,0);
+\draw(0,-1.5)--(0,3);
+
+\draw[color=gray,line width=1pt, densely dotted](-2.5,1.25)--++(2.5,0);
+\draw[color=gray,line width=1pt, densely dotted](-1,2.00)--++(1,0);
+\draw[color=gray,line width=1pt, densely dotted](-1,2.75)--++(3.5,0);
+
+\draw[->,line width=0.75pt,color=red](-0.5,1.25)--++(0,0.75);
+\draw[->,line width=0.75pt,color=red](-0.5,2.00)--++(0,0.75);
+
+\draw[line width=1.5pt](-2.5,0)--(0,0)--(0,2.5)--(2.5,2.5);
+
+\end{tikzpicture}
+\end{document}
+
diff --git a/figs/potential.fig/potential_square_thermal.tikz.tex b/figs/potential.fig/potential_square_thermal.tikz.tex
new file mode 100644
index 0000000..3e8fd22
--- /dev/null
+++ b/figs/potential.fig/potential_square_thermal.tikz.tex
@@ -0,0 +1,17 @@
+\documentclass{standalone}
+
+\usepackage{tikz}
+
+\begin{document}
+\begin{tikzpicture}
+
+\draw(-3,0)--(3,0);
+\draw(0,-1.5)--(0,3);
+
+\draw[color=gray,line width=1pt, densely dotted](-2.5,2.75)--++(5,0);
+
+\draw[line width=1.5pt](-2.5,0)--(0,0)--(0,2.5)--(2.5,2.5);
+
+\end{tikzpicture}
+\end{document}
+
diff --git a/libs/ian-presentation.cls b/libs/ian-presentation.cls
new file mode 100644
index 0000000..91bd487
--- /dev/null
+++ b/libs/ian-presentation.cls
@@ -0,0 +1,187 @@
+%%
+%% Ian's presentation class
+%%
+
+%% TeX format
+\NeedsTeXFormat{LaTeX2e}[1995/12/01]
+
+%% class name
+\ProvidesClass{ian-presentation}[2017/09/29]
+
+\def\ian@defaultoptions{
+  \pagestyle{plain}
+  \RequirePackage{color}
+  \RequirePackage{amssymb}
+}
+
+%% paper dimensions
+\setlength\paperheight{240pt}
+\setlength\paperwidth{427pt}
+
+%% fonts
+\input{size11.clo}
+\DeclareOldFontCommand{\rm}{\normalfont\rmfamily}{\mathrm}
+\DeclareOldFontCommand{\sf}{\normalfont\sffamily}{\mathsf}
+\DeclareOldFontCommand{\tt}{\normalfont\ttfamily}{\mathtt}
+\DeclareOldFontCommand{\bf}{\normalfont\bfseries}{\mathbf}
+\DeclareOldFontCommand{\it}{\normalfont\itshape}{\mathit}
+
+%% text dimensions
+\textheight=208pt
+\textwidth=370pt
+\hoffset=-1in
+\voffset=-1in
+\oddsidemargin=24pt
+\evensidemargin=24pt
+\topmargin=8pt
+\headheight=0pt
+\headsep=0pt
+\marginparsep=0pt
+\marginparwidth=0pt
+\footskip=16pt
+
+
+%% remove default skips
+\parindent=0pt
+\parskip=0pt
+\baselineskip=0pt
+
+%% something is wrong with \thepage, redefine it
+\gdef\thepage{\the\c@page}
+
+%% correct vertical alignment at the end of a document
+\AtEndDocument{
+  % save total slide count
+  \immediate\write\@auxout{\noexpand\gdef\noexpand\slidecount{\thepage}}
+  \vfill
+  \eject
+}
+
+
+%% footer
+\def\ps@plain{
+  \def\@oddhead{}
+  \def\@evenhead{\@oddhead}
+  \def\@oddfoot{\tiny\hfill\thepage/\safe\slidecount\hfill}
+  \def\@evenfoot{\@oddfoot}
+}
+\def\ps@empty{
+  \def\@oddhead{}
+  \def\@evenhead{\@oddhead}
+  \def\@oddfoot{}
+  \def\@evenfoot{\@oddfoot}
+}
+
+
+%% title of slide
+\def\title#1{
+  \hfil{\bf\large #1}\par
+  \hfil\vrule width0.75\textwidth height0.3pt\par
+  \vskip5pt
+}
+
+
+%% hyperlinks
+% hyperlinkcounter
+\newcounter{lncount}
+% hyperref anchor
+\def\hrefanchor{%
+\stepcounter{lncount}%
+\hypertarget{ln.\thelncount}{}%
+}
+
+%% define a command and write it to aux file
+\def\outdef#1#2{%
+  % define command%
+  \expandafter\xdef\csname #1\endcsname{#2}%
+  % hyperlink number%
+  \expandafter\xdef\csname #1@hl\endcsname{\thelncount}%
+  % write command to aux%
+  \immediate\write\@auxout{\noexpand\expandafter\noexpand\gdef\noexpand\csname #1\endcsname{\csname #1\endcsname}}%
+  \immediate\write\@auxout{\noexpand\expandafter\noexpand\gdef\noexpand\csname #1@hl\endcsname{\thelncount}}%
+}
+
+%% can call commands even when they are not defined
+\def\safe#1{%
+  \ifdefined#1%
+    #1%
+  \else%
+    {\color{red}\bf?}%
+  \fi%
+}
+
+
+%% itemize
+\newlength\itemizeskip
+% left margin for items
+\setlength\itemizeskip{20pt}
+\newlength\itemizeseparator
+% space between the item symbol and the text
+\setlength\itemizeseparator{5pt}
+% penalty preceding an itemize
+\def\itemizepenalty{0}
+% counter counting the itemize level
+\newcounter{itemizecount}
+
+% item symbol
+\def\itemizept#1{
+  \ifnum#1=1
+    \textbullet
+  \else
+    $\scriptstyle\blacktriangleright$
+  \fi
+}
+
+\newlength\current@itemizeskip
+\setlength\current@itemizeskip{0pt}
+\def\itemize{
+  \par\penalty\itemizepenalty\medskip\penalty\itemizepenalty
+  \addtocounter{itemizecount}{1}
+  \addtolength\current@itemizeskip{\itemizeskip}
+  \leftskip\current@itemizeskip
+}
+\def\enditemize{
+  \addtocounter{itemizecount}{-1}
+  \addtolength\current@itemizeskip{-\itemizeskip}
+  \par\leftskip\current@itemizeskip
+  \medskip
+}
+\newlength\itempt@total
+\def\item{
+  \settowidth\itempt@total{\itemizept\theitemizecount}
+  \addtolength\itempt@total{\itemizeseparator}
+  \par
+  \medskip
+  \hskip-\itempt@total\itemizept\theitemizecount\hskip\itemizeseparator
+}
+
+%% enumerate
+\newcounter{enumerate@count}
+\def\enumerate{
+  \setcounter{enumerate@count}0
+  \let\olditem\item
+  \let\olditemizept\itemizept
+  \def\item{
+    % counter
+    \stepcounter{enumerate@count}
+    % set header
+    \def\itemizept{\theenumerate@count.}
+    % hyperref anchor
+    \hrefanchor
+    % define tag (for \label)
+    \xdef\tag{\theenumerate@count}
+    \olditem
+  }
+  \itemize
+}
+\def\endenumerate{
+  \enditemize
+  \let\item\olditem
+  \let\itemizept\olditemizept
+}
+
+
+%% end
+\ian@defaultoptions
+
+\endinput