diff options
| -rw-r--r-- | Jauslin_SMM121_2019.tex | 148 | ||||
| -rw-r--r-- | Makefile | 50 | ||||
| -rw-r--r-- | README | 42 | ||||
| -rw-r--r-- | figs/anim_laser/Makefile | 12 | ||||
| -rw-r--r-- | figs/anim_laser/animate.py | 103 | ||||
| -rw-r--r-- | figs/anim_laser/laser.jl | 391 | ||||
| -rw-r--r-- | figs/libs/Makefile | 28 | ||||
| l--------- | figs/potential.fig/Makefile | 1 | ||||
| -rw-r--r-- | figs/potential.fig/potential.tikz.tex | 17 | ||||
| -rw-r--r-- | figs/potential.fig/potential_square.tikz.tex | 17 | ||||
| -rw-r--r-- | figs/potential.fig/potential_square_photonic.tikz.tex | 22 | ||||
| -rw-r--r-- | figs/potential.fig/potential_square_thermal.tikz.tex | 17 | ||||
| -rw-r--r-- | libs/ian-presentation.cls | 187 |
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 @@ -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 |