From d0c870377cc4e0388a343a6297f7ec5fa54adf16 Mon Sep 17 00:00:00 2001
From: Ian Jauslin <jauslin@ias.edu>
Date: Mon, 13 May 2019 22:20:51 -0400
Subject: As presented at MSU on 2019-05-02

---
 figs/contour.fig/Makefile         | 30 ++++++++++++++
 figs/contour.fig/contour.py       | 82 +++++++++++++++++++++++++++++++++++++++
 figs/contour.fig/contour.tikz.tex | 34 ++++++++++++++++
 3 files changed, 146 insertions(+)
 create mode 100644 figs/contour.fig/Makefile
 create mode 100644 figs/contour.fig/contour.py
 create mode 100644 figs/contour.fig/contour.tikz.tex

(limited to 'figs/contour.fig')

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}
-- 
cgit v1.2.3-70-g09d2