Ian Jauslin
summaryrefslogtreecommitdiff
path: root/figs
diff options
context:
space:
mode:
authorIan Jauslin <jauslin@ias.edu>2019-05-14 11:28:19 -0400
committerIan Jauslin <jauslin@ias.edu>2019-05-14 11:28:19 -0400
commitf34f94a8fd01975cf2e17a8a8a7e04665aed1cca (patch)
tree102901914bcfd8f117afa2d9127c429d7ad64573 /figs
As presented at tthe 121st Statistical Mechanics Meeting on 2019-05-14HEADv1.0master
Diffstat (limited to 'figs')
-rw-r--r--figs/anim_laser/Makefile12
-rw-r--r--figs/anim_laser/animate.py103
-rw-r--r--figs/anim_laser/laser.jl391
-rw-r--r--figs/libs/Makefile28
l---------figs/potential.fig/Makefile1
-rw-r--r--figs/potential.fig/potential.tikz.tex17
-rw-r--r--figs/potential.fig/potential_square.tikz.tex17
-rw-r--r--figs/potential.fig/potential_square_photonic.tikz.tex22
-rw-r--r--figs/potential.fig/potential_square_thermal.tikz.tex17
9 files changed, 608 insertions, 0 deletions
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}
+