171 lines
7.3 KiB
Julia
171 lines
7.3 KiB
Julia
# 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
|
|
|
|
# asymptotic airy functions
|
|
# specify a branch cut for the fractional power
|
|
function airyai_asym(x,cut)
|
|
if(abs(real(pow(x,3/2,cut)))<airy_threshold)
|
|
return(exp(2/3*pow(x,3/2,cut))*airyai(x))
|
|
else
|
|
ret=0
|
|
for n in 0:airy_order
|
|
ret+=gamma(n+5/6)*gamma(n+1/6)*(-3/4)^n/(4*pi^(3/2)*factorial(n)*pow(x,3*n/2+1/4,cut))
|
|
end
|
|
return ret
|
|
end
|
|
end
|
|
function airyaiprime_asym(x,cut)
|
|
if(abs(real(pow(x,3/2,cut)))<airy_threshold)
|
|
return(exp(2/3*pow(x,3/2,cut))*airyaiprime(x))
|
|
else
|
|
ret=0
|
|
for n in 0:airy_order
|
|
ret+=gamma(n+5/6)*gamma(n+1/6)*(-3/4)^n/(4*pi^(3/2)*factorial(n))*(-1/pow(x,3*n/2-1/4,cut)-(3/2*n+1/4)/pow(x,3*n/2+5/4,cut))
|
|
end
|
|
return ret
|
|
end
|
|
end
|
|
|
|
# solutions of (-\Delta+U-ip)phi=0
|
|
# assume that p has an infinitesimal real part (and adjust the branch cuts appropriately)
|
|
function phi(p,x,E,U)
|
|
return(airyai_asym(2^(1/3)*exp(-1im*pi/3)*(E^(1/3)*x-E^(-2/3)*(U-1im*p)),pi))
|
|
end
|
|
function dphi(p,x,E,U)
|
|
return(2^(1/3)*exp(-1im*pi/3)*E^(1/3)*airyaiprime_asym(2^(1/3)*exp(-1im*pi/3)*(E^(1/3)*x-E^(-2/3)*(U-1im*p)),pi))
|
|
end
|
|
function eta(p,x,E,U)
|
|
return(exp(-1im*pi/3)*airyai_asym(-2^(1/3)*(E^(1/3)*x-E^(-2/3)*(U-1im*p)),pi/2))
|
|
end
|
|
function deta(p,x,E,U)
|
|
return(-2^(1/3)*exp(-1im*pi/3)*E^(1/3)*airyaiprime_asym(-2^(1/3)*(E^(1/3)*x-E^(-2/3)*(U-1im*p)),pi/2))
|
|
end
|
|
|
|
# Laplace transform of psi
|
|
# assume that p has an infinitesimal real part (and adjust the branch cuts appropriately)
|
|
# for example, (1im*p-U)^(3/2) becomes pow(1im*p-U,3/2,-pi/2) because when 1im*p is real negative, its square root should be imaginary positive
|
|
function f(p,x,k0,E,U)
|
|
T=2im*k0/(1im*k0-sqrt(2*U-k0*k0))
|
|
R=T-1
|
|
|
|
if x>=0
|
|
C2=-2im*T/(pow(-2im*p,1/2,pi/2)*phi(p,0,E,U)-dphi(p,0,E,U))*((sqrt(2*U-k0*k0)+pow(-2im*p,1/2,pi/2))/(-2im*p+k0*k0)-2im*(2*E)^(-1/3)*pi*quadgk(y -> (pow(-2im*p,1/2,pi/2)*eta(p,0,E,U)-deta(p,0,E,U))*phi(p,y,E,U)*exp(-sqrt(2*U-k0*k0)*y)*exp(sqrt(2)*2im/3*(pow(E^(1/3)*y+E^(-2/3)*(1im*p-U),3/2,-pi/2)-E^(-1)*pow(1im*p-U,3/2,-pi/2))),0,Inf)[1])
|
|
FT=4*(2*E)^(-1/3)*pi*(quadgk(y -> phi(p,x,E,U)*eta(p,y,E,U)*exp(-sqrt(2*U-k0*k0)*y)*exp(sqrt(2)*2im/3*(pow(E^(1/3)*x+E^(-2/3)*(1im*p-U),3/2,-pi/2)-pow(E^(1/3)*y+E^(-2/3)*(1im*p-U),3/2,-pi/2))),0,x)[1]+quadgk(y -> eta(p,x,E,U)*phi(p,y,E,U)*exp(-sqrt(2*U-k0*k0)*y)*exp(sqrt(2)*2im/3*(pow(E^(1/3)*y+E^(-2/3)*(1im*p-U),3/2,-pi/2)-pow(E^(1/3)*x+E^(-2/3)*(1im*p-U),3/2,-pi/2))),x,Inf)[1])
|
|
main=C2*phi(p,x,E,U)*exp(sqrt(2)*2im/3*(pow(E^(1/3)*x+E^(-2/3)*(1im*p-U),3/2,-pi/2)-E^(-1)*pow(1im*p-U,3/2,-pi/2)))+T*FT
|
|
|
|
# subtract the contribution of the pole, which will be added back in after the integration
|
|
pole=psi_pole(x,k0,E,U)/(p+1im*k0*k0/2)
|
|
return(main-pole)
|
|
else
|
|
C1=-2im*T*((sqrt(2*U-k0*k0)*phi(p,0,E,U)+dphi(p,0,E,U))/(-2im*p+k0*k0)/(pow(-2im*p,1/2,pi/2)*phi(p,0,E,U)-dphi(p,0,E,U))+quadgk(y -> phi(p,y,E,U)/(pow(-2im*p,1/2,pi/2)*phi(p,0,E,U)-dphi(p,0,E,U))*exp(-sqrt(2*U-k0*k0)*y)*exp(sqrt(2)*2im/3*(pow(E^(1/3)*y+E^(-2/3)*(1im*p-U),3/2,-pi/2)-E^(-1)*pow(1im*p-U,3/2,-pi/2))),0,Inf)[1])
|
|
FI=-2im*exp(1im*k0*x)/(-2im*p+k0*k0)
|
|
FR=-2im*exp(-1im*k0*x)/(-2im*p+k0*k0)
|
|
main=C1*exp(pow(-2im*p,1/2,pi/2)*x)+FI+R*FR
|
|
|
|
# subtract the contribution of the pole, which will be added back in after the integration
|
|
pole=psi_pole(x,k0,E,U)/(p+1im*k0*k0/2)
|
|
return(main-pole)
|
|
end
|
|
end
|
|
# its derivative
|
|
function df(p,x,k0,E,U)
|
|
T=2im*k0/(1im*k0-sqrt(2*U-k0*k0))
|
|
R=T-1
|
|
|
|
if x>=0
|
|
C2=-2im*T/(pow(-2im*p,1/2,pi/2)*phi(p,0,E,U)-dphi(p,0,E,U))*((sqrt(2*U-k0*k0)+pow(-2im*p,1/2,pi/2))/(-2im*p+k0*k0)-2im*(2*E)^(-1/3)*pi*quadgk(y -> (pow(-2im*p,1/2,pi/2)*eta(p,0,E,U)-deta(p,0,E,U))*phi(p,y,E,U)*exp(-sqrt(2*U-k0*k0)*y)*exp(sqrt(2)*2im/3*(pow(E^(1/3)*y+E^(-2/3)*(1im*p-U),3/2,-pi/2)-E^(-1)*pow(1im*p-U,3/2,-pi/2))),0,Inf)[1])
|
|
dFT=4*(2*E)^(-1/3)*pi*(quadgk(y -> dphi(p,x,E,U)*eta(p,y,E,U)*exp(-sqrt(2*U-k0*k0)*y)*exp(sqrt(2)*2im/3*(pow(E^(1/3)*x+E^(-2/3)*(1im*p-U),3/2,-pi/2)-pow(E^(1/3)*y+E^(-2/3)*(1im*p-U),3/2,-pi/2))),0,x)[1]+quadgk(y -> deta(p,x,E,U)*phi(p,y,E,U)*exp(-sqrt(2*U-k0*k0)*y)*exp(sqrt(2)*2im/3*(pow(E^(1/3)*y+E^(-2/3)*(1im*p-U),3/2,-pi/2)-pow(E^(1/3)*x+E^(-2/3)*(1im*p-U),3/2,-pi/2))),x,Inf)[1])
|
|
main=C2*dphi(p,x,E,U)*exp(sqrt(2)*2im/3*(pow(E^(1/3)*x+E^(-2/3)*(1im*p-U),3/2,-pi/2)-E^(-1)*pow(1im*p-U,3/2,-pi/2)))+T*dFT
|
|
|
|
# subtract the contribution of the pole, which will be added back in after the integration
|
|
pole=dpsi_pole(x,k0,E,U)/(p+1im*k0*k0/2)
|
|
return(main-pole)
|
|
else
|
|
C1=-2im*T*((sqrt(2*U-k0*k0)*phi(p,0,E,U)+dphi(p,0,E,U))/(-2im*p+k0*k0)/(pow(-2im*p,1/2,pi/2)*phi(p,0,E,U)-dphi(p,0,E,U))+quadgk(y -> phi(p,y,E,U)/(pow(-2im*p,1/2,pi/2)*phi(p,0,E,U)-dphi(p,0,E,U))*exp(-sqrt(2*U-k0*k0)*y)*exp(sqrt(2)*2im/3*(pow(E^(1/3)*y+E^(-2/3)*(1im*p-U),3/2,-pi/2)-E^(-1)*pow(1im*p-U,3/2,-pi/2))),0,Inf)[1])
|
|
dFI=2*k0*exp(1im*k0*x)/(-2im*p+k0*k0)
|
|
dFR=-2*k0*exp(-1im*k0*x)/(-2im*p+k0*k0)
|
|
main=C1*pow(-2im*p,1/2,pi/2)*exp(pow(-2im*p,1/2,pi/2)*x)+dFI+R*dFR
|
|
|
|
# subtract the contribution of the pole, which will be added back in after the integration
|
|
pole=dpsi_pole(x,k0,E,U)/(p+1im*k0*k0/2)
|
|
return(main-pole)
|
|
end
|
|
end
|
|
|
|
# psi (returns t,psi(x,t))
|
|
function psi(x,k0,E,U,p_npoints,p_cutoff)
|
|
fft=fourier_fft(f,x,k0,E,U,p_npoints,p_cutoff)
|
|
# add the contribution of the pole
|
|
for i in 1:p_npoints
|
|
fft[2][i]=fft[2][i]+psi_pole(x,k0,E,U)*exp(-1im*k0*k0/2*fft[1][i])
|
|
end
|
|
return(fft)
|
|
end
|
|
# its derivative
|
|
function dpsi(x,k0,E,U,p_npoints,p_cutoff)
|
|
fft=fourier_fft(df,x,k0,E,U,p_npoints,p_cutoff)
|
|
# add the contribution of the pole
|
|
for i in 1:p_npoints
|
|
fft[2][i]=fft[2][i]+dpsi_pole(x,k0,E,U)*exp(-1im*k0*k0/2*fft[1][i])
|
|
end
|
|
return(fft)
|
|
end
|
|
|
|
# compute Fourier transform by sampling and fft
|
|
function fourier_fft(A,x,k0,E,U,p_npoints,p_cutoff)
|
|
fun=zeros(Complex{Float64},p_npoints)
|
|
times=zeros(p_npoints)
|
|
|
|
# prepare fft
|
|
for i in 1:p_npoints
|
|
fun[i]=p_cutoff/pi*A(1im*(-p_cutoff+2*p_cutoff*(i-1)/p_npoints),x,k0,E,U)
|
|
times[i]=(i-1)*pi/p_cutoff
|
|
end
|
|
|
|
ifft!(fun)
|
|
|
|
# correct the phase
|
|
for i in 2:2:p_npoints
|
|
fun[i]=-fun[i]
|
|
end
|
|
return([times,fun])
|
|
end
|
|
|
|
# asymptotic value of psi
|
|
function psi_pole(x,k0,E,U)
|
|
if x>=0
|
|
return(1im*phi(-1im*k0*k0/2,x,E,U)*2*k0/(1im*k0*phi(-1im*k0*k0/2,0,E,U)+dphi(-1im*k0*k0/2,0,E,U))*exp(sqrt(2)*2im/3*(pow(E^(1/3)*x+E^(-2/3)*(k0*k0/2-U),3/2,-pi/2)-E^(-1)*pow(k0*k0/2-U,3/2,-pi/2))))
|
|
else
|
|
return((1im*k0*phi(-1im*k0*k0/2,0,E,U)-dphi(-1im*k0*k0/2,0,E,U))/(1im*k0*phi(-1im*k0*k0/2,0,E,U)+dphi(-1im*k0*k0/2,0,E,U))*exp(-1im*k0*x)+exp(1im*k0*x))
|
|
end
|
|
end
|
|
function dpsi_pole(x,k0,E,U)
|
|
if x>=0
|
|
return(1im*dphi(-1im*k0*k0/2,x,E,U)*2*k0/(1im*k0*phi(-1im*k0*k0/2,0,E,U)+dphi(-1im*k0*k0/2,0,E,U))*exp(sqrt(2)*2im/3*(pow(E^(1/3)*x+E^(-2/3)*(k0*k0/2-U),3/2,-pi/2)-E^(-1)*pow(k0*k0/2-U,3/2,-pi/2))))
|
|
else
|
|
return(-1im*k0*(1im*k0*phi(-1im*k0*k0/2,0,E,U)-dphi(-1im*k0*k0/2,0,E,U))/(1im*k0*phi(-1im*k0*k0/2,0,E,U)+dphi(-1im*k0*k0/2,0,E,U))*exp(-1im*k0*x)+1im*k0*exp(1im*k0*x))
|
|
end
|
|
end
|
|
|
|
# current
|
|
function J(ps,dps)
|
|
return(2*imag(conj(ps)*dps))
|
|
end
|
|
|
|
# complete computation of the current
|
|
function current(x,k0,E,U,p_npoints,p_cutoff)
|
|
ps=psi(x,k0,E,U,p_npoints,p_cutoff)
|
|
dps=dpsi(x,k0,E,U,p_npoints,p_cutoff)
|
|
Js=zeros(Complex{Float64},p_npoints)
|
|
for i in 1:p_npoints
|
|
Js[i]=J(ps[2][i],dps[2][i])
|
|
end
|
|
return(Js)
|
|
end
|