From 0425ca56937e60d194c475c29b3f145c0bac30bc Mon Sep 17 00:00:00 2001 From: Ian Jauslin Date: Fri, 12 Oct 2018 20:24:46 +0000 Subject: Initial commit --- figs/plots.fig/FN_base.jl | 170 ++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 170 insertions(+) create mode 100644 figs/plots.fig/FN_base.jl (limited to 'figs/plots.fig/FN_base.jl') diff --git a/figs/plots.fig/FN_base.jl b/figs/plots.fig/FN_base.jl new file mode 100644 index 0000000..cff6d8e --- /dev/null +++ b/figs/plots.fig/FN_base.jl @@ -0,0 +1,170 @@ +# 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)))=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 -- cgit v1.2.3-70-g09d2