diff options
Diffstat (limited to 'src/chebyshev.jl')
| -rw-r--r-- | src/chebyshev.jl | 279 |
1 files changed, 279 insertions, 0 deletions
diff --git a/src/chebyshev.jl b/src/chebyshev.jl new file mode 100644 index 0000000..28c8f1f --- /dev/null +++ b/src/chebyshev.jl @@ -0,0 +1,279 @@ +## Copyright 2021 Ian Jauslin +## +## Licensed under the Apache License, Version 2.0 (the "License"); +## you may not use this file except in compliance with the License. +## You may obtain a copy of the License at +## +## http://www.apache.org/licenses/LICENSE-2.0 +## +## Unless required by applicable law or agreed to in writing, software +## distributed under the License is distributed on an "AS IS" BASIS, +## WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +## See the License for the specific language governing permissions and +## limitations under the License. + +# Chebyshev expansion +@everywhere function chebyshev(a,taus,weights,P,N,J,nu) + out=zeros(Float64,J*(P+1)) + for zeta in 0:J-1 + for n in 0:P + for j in 1:N + out[zeta*(P+1)+n+1]+=(2-(n==0 ? 1 : 0))/2*weights[2][j]*cos(n*pi*(1+weights[1][j])/2)*a[zeta*N+j]/(1-(taus[zeta+2]-taus[zeta+1])/2*sin(pi*weights[1][j]/2)+(taus[zeta+2]+taus[zeta+1])/2)^nu + end + end + end + + return out +end + +# evaluate function from Chebyshev expansion +@everywhere function chebyshev_eval(Fa,x,taus,chebyshev,P,J,nu) + # change variable + tau=(1-x)/(1+x) + + out=0. + for zeta in 0:J-1 + # check that the spline is right + if tau<taus[zeta+2] && tau>=taus[zeta+1] + for n in 0:P + out+=Fa[zeta*(P+1)+n+1]*chebyshev[n+1]((2*tau-(taus[zeta+1]+taus[zeta+2]))/(taus[zeta+2]-taus[zeta+1])) + end + break + end + end + + return (1+tau)^nu*out +end + + +# convolution +# input the Chebyshev expansion of a and b, as well as the A matrix +@everywhere function conv_chebyshev(Fa,Fb,A) + out=zeros(Float64,length(A)) + for i in 1:length(A) + out[i]=dot(Fa,A[i]*Fb) + end + return out +end + +# <ab> +@everywhere function avg_chebyshev(Fa,Fb,A0) + return dot(Fa,A0*Fb) +end + +# 1_n * a +@everywhere function conv_one_chebyshev(n,zetapp,Fa,A,taus,weights,P,N,J,nu1) + out=zeros(Float64,N*J) + for m in 1:N*J + for l in 0:P + for p in 1:J*(P+1) + out[m]+=(2-(l==0 ? 1 : 0))/2*weights[2][n]*cos(l*pi*(1+weights[1][n])/2)/(1-(taus[zetapp+2]-taus[zetapp+1])/2*sin(pi*weights[1][n]/2)+(taus[zetapp+2]+taus[zetapp+1])/2)^nu1*A[m][zetapp*(P+1)+l+1,p]*Fa[p] + end + end + end + return out +end +# a * v +@everywhere function conv_v_chebyshev(a,Upsilon,k,taus,weights,N,J) + out=zeros(Float64,J*N) + for i in 1:J*N + for zetap in 0:J-1 + for j in 1:N + xj=weights[1][j] + out[i]+=(taus[zetap+2]-taus[zetap+1])/(32*pi)*weights[2][j]*cos(pi*xj/2)*(1+k[zetap*N+j])^2*k[zetap*N+j]*a[zetap*N+j]*Upsilon[zetap*N+j][i] + end + end + end + return out +end +@everywhere function conv_one_v_chebyshev(n,zetap,Upsilon,k,taus,weights,N,J) + out=zeros(Float64,J*N) + xj=weights[1][n] + for i in 1:J*N + out[i]=(taus[zetap+2]-taus[zetap+1])/(32*pi)*weights[2][n]*cos(pi*xj/2)*(1+k[zetap*N+n])^2*k[zetap*N+n]*Upsilon[zetap*N+n][i] + end + return out +end + +# <av> +@everywhere function avg_v_chebyshev(a,Upsilon0,k,taus,weights,N,J) + out=0. + for zetap in 0:J-1 + for j in 1:N + xj=weights[1][j] + out+=(taus[zetap+2]-taus[zetap+1])/(32*pi)*weights[2][j]*cos(pi*xj/2)*(1+k[zetap*N+j])^2*k[zetap*N+j]*a[zetap*N+j]*Upsilon0[zetap*N+j] + end + end + return out +end +# <1_nv> +@everywhere function avg_one_v_chebyshev(n,zetap,Upsilon0,k,taus,weights,N) + xj=weights[1][n] + return (taus[zetap+2]-taus[zetap+1])/(32*pi)*weights[2][n]*cos(pi*xj/2)*(1+k[zetap*N+n])^2*k[zetap*N+n]*Upsilon0[zetap*N+n] +end + +# compute \int dq dxi u1(k-xi)u2(q)u3(xi)u4(k-q)u5(xi-q) +@everywhere function double_conv_S_chebyshev(FU1,FU2,FU3,FU4,FU5,Abar) + out=zeros(Float64,length(Abar)) + for i in 1:length(Abar) + for j1 in 1:length(FU1) + for j2 in 1:length(FU2) + for j3 in 1:length(FU3) + for j4 in 1:length(FU4) + for j5 in 1:length(FU5) + out[i]+=Abar[i][j1,j2,j3,j4,j5]*FU1[j1]*FU2[j2]*FU3[j3]*FU4[j4]*FU5[j5] + end + end + end + end + end + end + return out +end + + +# compute A +@everywhere function Amat(k,weights,taus,T,P,N,J,nua,nub) + out=Array{Array{Float64,2},1}(undef,J*N) + for i in 1:J*N + out[i]=zeros(Float64,J*(P+1),J*(P+1)) + for zeta in 0:J-1 + for n in 0:P + for zetap in 0:J-1 + for m in 0:P + out[i][zeta*(P+1)+n+1,zetap*(P+1)+m+1]=1/(pi^2*k[i])*integrate_legendre(tau->(1-tau)/(1+tau)^(3-nua)*T[n+1]((2*tau-(taus[zeta+1]+taus[zeta+2]))/(taus[zeta+2]-taus[zeta+1]))*(alpham(k[i],tau)>taus[zetap+2] || alphap(k[i],tau)<taus[zetap+1] ? 0. : integrate_legendre(sigma->(1-sigma)/(1+sigma)^(3-nub)*T[m+1]((2*sigma-(taus[zetap+1]+taus[zetap+2]))/(taus[zetap+2]-taus[zetap+1])),max(taus[zetap+1],alpham(k[i],tau)),min(taus[zetap+2],alphap(k[i],tau)),weights)),taus[zeta+1],taus[zeta+2],weights) + end + end + end + end + end + + return out +end + +# compute Upsilon +@everywhere function Upsilonmat(k,v,weights) + out=Array{Array{Float64,1},1}(undef,length(k)) + for i in 1:length(k) + out[i]=Array{Float64,1}(undef,length(k)) + for j in 1:length(k) + out[i][j]=integrate_legendre(s->s*v(s)/k[j],abs(k[j]-k[i]),k[j]+k[i],weights) + end + end + return out +end +@everywhere function Upsilon0mat(k,v,weights) + out=Array{Float64,1}(undef,length(k)) + for j in 1:length(k) + out[j]=2*k[j]*v(k[j]) + end + return out +end + +# alpha_- +@everywhere function alpham(k,t) + return (1-k-(1-t)/(1+t))/(1+k+(1-t)/(1+t)) +end +# alpha_+ +@everywhere function alphap(k,t) + return (1-abs(k-(1-t)/(1+t)))/(1+abs(k-(1-t)/(1+t))) +end + + +# compute \bar A +@everywhere function barAmat(k,weights,taus,T,P,N,J,nu1,nu2,nu3,nu4,nu5) + out=zeros(Float64,J*(P+1),J*(P+1),J*(P+1),J*(P+1),J*(P+1)) + for zeta1 in 0:J-1 + for n1 in 0:P + for zeta2 in 0:J-1 + for n2 in 0:P + for zeta3 in 0:J-1 + for n3 in 0:P + for zeta4 in 0:J-1 + for n4 in 0:P + for zeta5 in 0:J-1 + for n5 in 0:P + out[zeta1*(P+1)+n1+1, + zeta2*(P+1)+n2+1, + zeta3*(P+1)+n3+1, + zeta4*(P+1)+n4+1, + zeta5*(P+1)+n5+1]=1/((2*pi)^5*k^2)*integrate_legendre(tau->barAmat_int1(tau,k,taus,T,weights,nu1,nu2,nu3,nu4,nu5,zeta1,zeta2,zeta3,zeta4,zeta5,n1,n2,n3,n4,n5),taus[zeta1+1],taus[zeta1+2],weights) + end + end + end + end + end + end + end + end + end + end + + return out +end +@everywhere function barAmat_int1(tau,k,taus,T,weights,nu1,nu2,nu3,nu4,nu5,zeta1,zeta2,zeta3,zeta4,zeta5,n1,n2,n3,n4,n5) + if(alpham(k,tau)<taus[zeta2+2] && alphap(k,tau)>taus[zeta2+1]) + return 2*(1-tau)/(1+tau)^(3-nu1)*T[n1+1]((2*tau-(taus[zeta1+1]+taus[zeta1+2]))/(taus[zeta1+2]-taus[zeta1+1]))*integrate_legendre(sigma->barAmat_int2(tau,sigma,k,taus,T,weights,nu2,nu3,nu4,nu5,zeta2,zeta3,zeta4,zeta5,n2,n3,n4,n5),max(taus[zeta2+1],alpham(k,tau)),min(taus[zeta2+2],alphap(k,tau)),weights) + else + return 0. + end +end +@everywhere function barAmat_int2(tau,sigma,k,taus,T,weights,nu2,nu3,nu4,nu5,zeta2,zeta3,zeta4,zeta5,n2,n3,n4,n5) + return 2*(1-sigma)/(1+sigma)^(3-nu2)*T[n2+1]((2*sigma-(taus[zeta2+1]+taus[zeta2+2]))/(taus[zeta2+2]-taus[zeta2+1]))*integrate_legendre(taup->barAmat_int3(tau,sigma,taup,k,taus,T,weights,nu3,nu4,nu5,zeta3,zeta4,zeta5,n3,n4,n5),taus[zeta3+1],taus[zeta3+2],weights) +end +@everywhere function barAmat_int3(tau,sigma,taup,k,taus,T,weights,nu3,nu4,nu5,zeta3,zeta4,zeta5,n3,n4,n5) + if(alpham(k,taup)<taus[zeta4+2] && alphap(k,taup)>taus[zeta4+1]) + return 2*(1-taup)/(1+taup)^(3-nu3)*T[n3+1]((2*taup-(taus[zeta3+1]+taus[zeta3+2]))/(taus[zeta3+2]-taus[zeta3+1]))*integrate_legendre(sigmap->barAmat_int4(tau,sigma,taup,sigmap,k,taus,T,weights,nu4,nu5,zeta4,zeta5,n4,n5),max(taus[zeta4+1],alpham(k,taup)),min(taus[zeta4+2],alphap(k,taup)),weights) + else + return 0. + end +end +@everywhere function barAmat_int4(tau,sigma,taup,sigmap,k,taus,T,weights,nu4,nu5,zeta4,zeta5,n4,n5) + return 2*(1-sigmap)/(1+sigmap)^(3-nu4)*T[n4+1]((2*sigma-(taus[zeta4+1]+taus[zeta4+2]))/(taus[zeta4+2]-taus[zeta4+1]))*integrate_legendre(theta->barAmat_int5(tau,sigma,taup,sigmap,theta,k,taus,T,weights,nu5,zeta5,n5),0.,2*pi,weights) +end +@everywhere function barAmat_int5(tau,sigma,taup,sigmap,theta,k,taus,T,weights,nu5,zeta5,n5) + R=barAmat_R((1-sigma)/(1+sigma),(1-tau)/(1+tau),(1-sigmap)/(1+sigmap),(1-taup)/(1+taup),theta,k) + if((1-R)/(1+R)<taus[zeta5+2] && (1-R)/(1+R)>taus[zeta5+1]) + return (2/(2+R))^nu5*T[n5+1]((2*(1-R)/(1+R)-(taus[zeta5+1]+taus[zeta5+2]))/(taus[zeta5+2]-taus[zeta5+1])) + else + return 0. + end +end +# R(s,t,s',t,theta,k) +@everywhere function barAmat_R(s,t,sp,tp,theta,k) + return sqrt(k^2*(s^2+t^2+sp^2+tp^2)-k^4-(s^2-t^2)*(sp^2-tp^2)-sqrt((4*k^2*s^2-(k^2+s^2-t^2)^2)*(4*k^2*sp^2-(k^2+sp^2-tp^2)^2))*cos(theta))/(sqrt(2)*k) +end + +# compute Chebyshev polynomials +@everywhere function chebyshev_polynomials(P) + T=Array{Polynomial}(undef,P+1) + T[1]=Polynomial([1]) + T[2]=Polynomial([0,1]) + for n in 1:P-1 + # T_n + T[n+2]=2*T[2]*T[n+1]-T[n] + end + + return T +end + +# compute \int f*u dk/(2*pi)^3 +@everywhere function integrate_f_chebyshev(f,u,k,taus,weights,N,J) + out=0. + for zeta in 0:J-1 + for i in 1:N + out+=(taus[zeta+2]-taus[zeta+1])/(16*pi)*weights[2][i]*cos(pi*weights[1][i]/2)*(1+k[zeta*N+i])^2*k[zeta*N+i]^2*u[zeta*N+i]*f(k[zeta*N+i]) + end + end + return out +end + +@everywhere function inverse_fourier_chebyshev(u,x,k,taus,weights,N,J) + out=0. + for zeta in 0:J-1 + for j in 1:N + out+=(taus[zeta+2]-taus[zeta+1])/(16*pi*x)*weights[2][j]*cos(pi*weights[1][j]/2)*(1+k[zeta*N+j])^2*k[zeta*N+j]*u[zeta*N+j]*sin(k[zeta*N+j]*x) + end + end + return out +end |