1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
|
# \eta
function eta(x,t,weights,d,v)
if d==2
return integrate_chebyshev(y->4*((x+t)*y+abs(x-t)*(1-y))*v((x+t)*y+abs(x-t)*(1-y))/sqrt(((x+t)*y+abs(x-t)*(2-y))*((x+t)*(1+y)+abs(x-t)*(1-y))),0,1,length(weights))
elseif d==3
return (x>t ? 2*t/x : 2)* integrate_legendre(y->2*pi*((x+t)*y+abs(x-t)*(1-y))*v((x+t)*y+abs(x-t)*(1-y)),0,1,weights)
end
end
# initialize V and Eta
function init_veta(weights,d,v)
order=length(weights[1])
V=Array{Complex{Float64}}(undef,order)
Eta=Array{Array{Complex{Float64}}}(undef,order)
Eta0=Array{Complex{Float64}}(undef,order)
V0=v(0)
for i in 1:order
k=(1-weights[1][i])/(1+weights[1][i])
V[i]=v(k)
Eta[i]=Array{Complex{Float64}}(undef,order)
for j in 1:order
y=(weights[1][j]+1)/2
Eta[i][j]=eta(k,(1-y)/y,weights,d,v)
end
y=(weights[1][i]+1)/2
Eta0[i]=eta(0,(1-y)/y,weights,d,v)
end
return(V,V0,Eta,Eta0)
end
# inverse Fourier transform
function u_x(x,u,weights,d)
order=length(weights[1])
if d==2
out=integrate_legendre_sampled(y->(1-y)/y^3*besselj(0,x*(1-y)/y)/(2*pi),u,0,1,weights)
elseif d==3
out=integrate_legendre_sampled(y->(1-y)/y^3*sin(x*(1-y)/y)/x/(2*pi^2),u,0,1,weights)
end
return out
end
|
