#!/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))