diff options
Diffstat (limited to 'src/simpleq-hardcore.jl')
| -rw-r--r-- | src/simpleq-hardcore.jl | 573 |
1 files changed, 573 insertions, 0 deletions
diff --git a/src/simpleq-hardcore.jl b/src/simpleq-hardcore.jl new file mode 100644 index 0000000..ca64f78 --- /dev/null +++ b/src/simpleq-hardcore.jl @@ -0,0 +1,573 @@ +## 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. + +# compute energy as a function of rho +function simpleq_hardcore_energy_rho(rhos,taus,P,N,J,maxiter,tolerance) + ## spawn workers + # number of workers + nw=nworkers() + # split jobs among workers + work=Array{Array{Int64,1},1}(undef,nw) + # init empty arrays + for p in 1:nw + work[p]=zeros(0) + end + for j in 1:length(rhos) + append!(work[j%nw+1],j) + end + + # initialize vectors + (weights,weights_gL,r,T)=simpleq_hardcore_init(taus,P,N,J) + + # initial guess + u0s=Array{Array{Float64}}(undef,length(rhos)) + e0s=Array{Float64}(undef,length(rhos)) + for j in 1:length(rhos) + u0s[j]=Array{Float64}(undef,N*J) + for i in 1:N*J + u0s[j][i]=1/(1+r[i]^2)^2 + end + e0s[j]=pi*rhos[j] + end + + + # save result from each task + us=Array{Array{Float64}}(undef,length(rhos)) + es=Array{Float64}(undef,length(rhos)) + err=Array{Float64}(undef,length(rhos)) + + count=0 + # for each worker + @sync for p in 1:nw + # for each task + @async for j in work[p] + count=count+1 + if count>=nw + progress(count,length(rhos),10000) + end + # run the task + (us[j],es[j],err[j])=remotecall_fetch(simpleq_hardcore_hatu,workers()[p],u0s[j],e0s[j],rhos[j],r,taus,T,weights,weights_gL,P,N,J,maxiter,tolerance) + end + end + + for j in 1:length(rhos) + @printf("% .15e % .15e % .15e\n",rhos[j],es[j],err[j]) + end +end + +# compute u(x) +function simpleq_hardcore_ux(rho,taus,P,N,J,maxiter,tolerance) + # initialize vectors + (weights,weights_gL,r,T)=simpleq_hardcore_init(taus,P,N,J) + + # initial guess + u0=Array{Float64}(undef,N*J) + for i in 1:N*J + u0[i]=1/(1+r[i]^2)^2 + end + e0=pi*rho + + (u,e,err)=simpleq_hardcore_hatu(u0,e0,rho,r,taus,T,weights,weights_gL,P,N,J,maxiter,tolerance) + + for i in 1:N*J + @printf("% .15e % .15e\n",r[i],u[i]) + end +end + +# compute condensate fraction as a function of rho +function simpleq_hardcore_condensate_fraction_rho(rhos,taus,P,N,J,maxiter,tolerance) + ## spawn workers + # number of workers + nw=nworkers() + # split jobs among workers + work=Array{Array{Int64,1},1}(undef,nw) + # init empty arrays + for p in 1:nw + work[p]=zeros(0) + end + for j in 1:length(rhos) + append!(work[j%nw+1],j) + end + + # initialize vectors + (weights,weights_gL,r,T)=simpleq_hardcore_init(taus,P,N,J) + + # initial guess + u0s=Array{Array{Float64}}(undef,length(rhos)) + e0s=Array{Float64}(undef,length(rhos)) + for j in 1:length(rhos) + u0s[j]=Array{Float64}(undef,N*J) + for i in 1:N*J + u0s[j][i]=1/(1+r[i]^2)^2 + end + e0s[j]=pi*rhos[j] + end + + + # save result from each task + us=Array{Array{Float64}}(undef,length(rhos)) + es=Array{Float64}(undef,length(rhos)) + err=Array{Float64}(undef,length(rhos)) + + count=0 + # for each worker + @sync for p in 1:nw + # for each task + @async for j in work[p] + count=count+1 + if count>=nw + progress(count,length(rhos),10000) + end + # run the task + (us[j],es[j],err[j])=remotecall_fetch(simpleq_hardcore_hatu,workers()[p],u0s[j],e0s[j],rhos[j],r,taus,T,weights,weights_gL,P,N,J,maxiter,tolerance) + end + end + + for j in 1:length(rhos) + du=-inv(simpleq_hardcore_DXi(us[j],es[j],rhos[j],r,taus,T,weights,weights_gL,P,N,J))*simpleq_hardcore_dXidmu(us[j],es[j],rhos[j],r,taus,T,weights,weights_gL,P,N,J) + @printf("% .15e % .15e % .15e\n",rhos[j],du[N*J+1],err[j]) + end +end + + +# initialize computation +@everywhere function simpleq_hardcore_init(taus,P,N,J) + # Gauss-Legendre weights + weights=gausslegendre(N) + weights_gL=gausslaguerre(N) + + # r + r=Array{Float64}(undef,J*N) + for zeta in 0:J-1 + for j in 1:N + xj=weights[1][j] + # set kj + r[zeta*N+j]=(2+(taus[zeta+2]-taus[zeta+1])*sin(pi*xj/2)-(taus[zeta+2]+taus[zeta+1]))/(2-(taus[zeta+2]-taus[zeta+1])*sin(pi*xj/2)+taus[zeta+2]+taus[zeta+1]) + end + end + + # Chebyshev polynomials + T=chebyshev_polynomials(P) + + return (weights,weights_gL,r,T) +end + +# compute u using chebyshev expansions +@everywhere function simpleq_hardcore_hatu(u0,e0,rho,r,taus,T,weights,weights_gL,P,N,J,maxiter,tolerance) + # init + vec=Array{Float64}(undef,J*N+1) + for i in 1:J*N + vec[i]=u0[i] + end + vec[J*N+1]=e0 + + # quantify relative error + error=Inf + + # run Newton algorithm + for i in 1:maxiter-1 + u=vec[1:J*N] + e=vec[J*N+1] + new=vec-inv(simpleq_hardcore_DXi(u,e,rho,r,taus,T,weights,weights_gL,P,N,J))*simpleq_hardcore_Xi(u,e,rho,r,taus,T,weights,weights_gL,P,N,J) + + error=norm(new-vec)/norm(new) + if(error<tolerance) + vec=new + break + end + + vec=new + end + + return(vec[1:J*N],vec[J*N+1],error) +end + +# Xi +@everywhere function simpleq_hardcore_Xi(u,e,rho,r,taus,T,weights,weights_gL,P,N,J) + out=Array{Float64}(undef,J*N+1) + FU=chebyshev(u,taus,weights,P,N,J,4) + + #D's + d0=D0(e,rho,r,weights,N,J) + d1=D1(e,rho,r,taus,T,weights,weights_gL,P,N,J,4) + d2=D2(e,rho,r,taus,T,weights,weights_gL,P,N,J,4) + + # u + for i in 1:J*N + out[i]=d0[i]+dot(FU,d1[i])+dot(FU,d2[i]*FU)-u[i] + end + # e + out[J*N+1]=-e+2*pi*rho* + ((1+2*sqrt(abs(e)))-gamma0(e,rho,weights)-dot(FU,gamma1(e,rho,taus,T,weights,weights_gL,P,J,4))-dot(FU,gamma2(e,rho,taus,T,weights,weights_gL,P,J,4)*FU))/ + (1-8/3*pi*rho+rho^2*(gammabar0(weights)+dot(FU,gammabar1(taus,T,weights,P,J,4))+dot(FU,gammabar2(taus,T,weights,P,J,4)*FU))) + + return out +end +# DXi +@everywhere function simpleq_hardcore_DXi(u,e,rho,r,taus,T,weights,weights_gL,P,N,J) + out=Array{Float64,2}(undef,J*N+1,J*N+1) + FU=chebyshev(u,taus,weights,P,N,J,4) + + #D's + d0=D0(e,rho,r,weights,N,J) + d1=D1(e,rho,r,taus,T,weights,weights_gL,P,N,J,4) + d2=D2(e,rho,r,taus,T,weights,weights_gL,P,N,J,4) + dsed0=dseD0(e,rho,r,weights,N,J) + dsed1=dseD1(e,rho,r,taus,T,weights,weights_gL,P,N,J,4) + dsed2=dseD2(e,rho,r,taus,T,weights,weights_gL,P,N,J,4) + + # denominator of e + denom=1-8/3*pi*rho+rho^2*(gammabar0(weights)+dot(FU,gammabar1(taus,T,weights,P,J,4))+dot(FU,gammabar2(taus,T,weights,P,J,4)*FU)) + + for zetapp in 0:J-1 + for n in 1:N + one=zeros(Int64,J*N) + one[zetapp*N+n]=1 + Fone=chebyshev(one,taus,weights,P,N,J,4) + + for i in 1:J*N + # du/du + out[i,zetapp*N+n]=dot(Fone,d1[i])+2*dot(FU,d2[i]*Fone)-(zetapp*N+n==i ? 1 : 0) + # du/de + out[i,J*N+1]=(dsed0[i]+dot(FU,dsed1[i])+dot(FU,dsed2[i]*FU))/(2*sqrt(abs(e)))*(e>=0 ? 1 : -1) + end + # de/du + out[J*N+1,zetapp*N+n]=2*pi*rho* + (-dot(Fone,gamma1(e,rho,taus,T,weights,weights_gL,P,J,4))-2*dot(FU,gamma2(e,rho,taus,T,weights,weights_gL,P,J,4)*Fone))/denom- + 2*pi*rho* + ((1+2*sqrt(abs(e)))-gamma0(e,rho,weights)-dot(FU,gamma1(e,rho,taus,T,weights,weights_gL,P,J,4))-dot(FU,gamma2(e,rho,taus,T,weights,weights_gL,P,J,4)*FU))* + rho^2*(dot(Fone,gammabar1(taus,T,weights,P,J,4))+2*dot(FU,gammabar2(taus,T,weights,P,J,4)*Fone))/denom^2 + end + end + #de/de + out[J*N+1,J*N+1]=-1+2*pi*rho* + (2-dsedgamma0(e,rho,weights)-dot(FU,dsedgamma1(e,rho,taus,T,weights,weights_gL,P,J,4))-dot(FU,dsedgamma2(e,rho,taus,T,weights,weights_gL,P,J,4)*FU))/denom/ + (2*sqrt(abs(e)))*(e>=0 ? 1 : -1) + + return out +end + +# dXi/dmu +@everywhere function simpleq_hardcore_dXidmu(u,e,rho,r,taus,T,weights,weights_gL,P,N,J) + out=Array{Float64}(undef,J*N+1) + FU=chebyshev(u,taus,weights,P,N,J,4) + + #D's + dsmud0=dsmuD0(e,rho,r,weights,N,J) + dsmud1=dsmuD1(e,rho,r,taus,T,weights,weights_gL,P,N,J,4) + dsmud2=dsmuD2(e,rho,r,taus,T,weights,weights_gL,P,N,J,4) + + # u + for i in 1:J*N + out[i]=(dsmud0[i]+dot(FU,dsmud1[i])+dot(FU,dsmud2[i]*FU))/(4*sqrt(abs(e))) + end + # e + out[J*N+1]=2*pi*rho*(2/3+1/(2*sqrt(abs(e)))- + (dsmudgamma0(e,rho,weights)+dot(FU,dsmudgamma1(e,rho,taus,T,weights,weights_gL,P,J,4))+dot(FU,dsmudgamma2(e,rho,taus,T,weights,weights_gL,P,J,4)*FU))/(4*sqrt(abs(e))) + )/(1-8/3*pi*rho+rho^2*(gammabar0(weights)+dot(FU,gammabar1(taus,T,weights,P,J,4))+dot(FU,gammabar2(taus,T,weights,P,J,4)*FU))) + + return out +end + +# B's +@everywhere function B0(r) + return pi/12*(r-1)^2*(r+5) +end +@everywhere function B1(r,zeta,n,taus,T,weights,nu) + return (taus[zeta+1]>=(2-r)/r || taus[zeta+2]<=-r/(r+2) ? 0 : + 8*pi/(r+1)*integrate_legendre(tau-> + (1-(r-(1-tau)/(1+tau))^2)/(1+tau)^(3-nu)*T[n+1]((2*tau-(taus[zeta+1]+taus[zeta+2]))/(taus[zeta+2]-taus[zeta+1])),max(taus[zeta+1],-r/(r+2)) + ,min(taus[zeta+2],(2-r)/r),weights) + ) +end +@everywhere function B2(r,zeta,n,zetap,m,taus,T,weights,nu) + return 32*pi/(r+1)*integrate_legendre(tau-> + 1/(1+tau)^(3-nu)*T[n+1]((2*tau-(taus[zeta+1]+taus[zeta+2]))/(taus[zeta+2]-taus[zeta+1]))* + (taus[zetap+1]>=alphap(abs(r-(1-tau)/(1+tau))-2*tau/(1+tau),tau) || taus[zetap+2]<=alpham(1+r,tau) ? 0 : + integrate_legendre(s-> + 1/(1+s)^(3-nu)*T[m+1]((2*s-(taus[zetap+1]+taus[zetap+2]))/(taus[zetap+2]-taus[zetap+1])),max(taus[zetap+1] + ,alpham(1+r,tau)),min(taus[zetap+2],alphap(abs(r-(1-tau)/(1+tau))-2*tau/(1+tau),tau)),weights) + ) + ,taus[zeta+1],taus[zeta+2],weights) +end + +# D's +@everywhere function D0(e,rho,r,weights,N,J) + out=Array{Float64}(undef,J*N) + for i in 1:J*N + out[i]=exp(-2*sqrt(abs(e))*r[i])/(r[i]+1)+ + rho*sqrt(abs(e))/(r[i]+1)*integrate_legendre(s-> + (s+1)*B0(s)*(exp(-2*sqrt(abs(e))*(r[i]-s))-exp(-2*sqrt(abs(e))*(r[i]+s)))/2 + ,0,min(1,r[i]),weights)+ + (r[i]>=1 ? 0 : + rho*sqrt(abs(e))/(2*(r[i]+1))*(1-exp(-4*sqrt(abs(e))*r[i]))*integrate_legendre(s-> + (s+r[i]+1)*B0(s+r[i])*exp(-2*sqrt(abs(e))*s) + ,0,1-r[i],weights) + ) + end + return out +end +@everywhere function D1(e,rho,r,taus,T,weights,weights_gL,P,N,J,nu) + out=Array{Array{Float64}}(undef,J*N) + for i in 1:J*N + out[i]=Array{Float64}(undef,(P+1)*J) + for zeta in 0:J-1 + for n in 0:P + out[i][zeta*(P+1)+n+1]=rho*sqrt(abs(e))/(r[i]+1)*integrate_legendre(s-> + (s+1)*B1(s,zeta,n,taus,T,weights,nu)*(exp(-2*sqrt(abs(e))*(r[i]-s))-exp(-2*sqrt(abs(e))*(r[i]+s)))/2 + ,0,r[i],weights)+ + rho*sqrt(abs(e))/(2*(r[i]+1))*(1-exp(-4*sqrt(abs(e))*r[i]))*integrate_laguerre(s-> + (s+r[i]+1)*B1(s+r[i],zeta,n,taus,T,weights,nu) + ,2*sqrt(abs(e)),weights_gL) + end + end + end + + return out +end +@everywhere function D2(e,rho,r,taus,T,weights,weights_gL,P,N,J,nu) + out=Array{Array{Float64,2}}(undef,J*N) + for i in 1:J*N + out[i]=Array{Float64,2}(undef,(P+1)*J,(P+1)*J) + 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]=rho*sqrt(abs(e))/(r[i]+1)*integrate_legendre(s-> + (s+1)*B2(s,zeta,n,zetap,m,taus,T,weights,nu)*(exp(-2*sqrt(abs(e))*(r[i]-s))-exp(-2*sqrt(abs(e))*(r[i]+s)))/2 + ,0,r[i],weights)+ + rho*sqrt(abs(e))/(2*(r[i]+1))*(1-exp(-4*sqrt(abs(e))*r[i]))*integrate_laguerre(s-> + (s+r[i]+1)*B2(s+r[i],zeta,n,zetap,m,taus,T,weights,nu) + ,2*sqrt(abs(e)),weights_gL) + end + end + end + end + end + return out +end + +# dD/d sqrt(abs(e))'s +@everywhere function dseD0(e,rho,r,weights,N,J) + out=Array{Float64}(undef,J*N) + for i in 1:J*N + out[i]=-2*r[i]*exp(-2*sqrt(abs(e))*r[i])/(r[i]+1)+ + rho/(r[i]+1)*integrate_legendre(s-> + (s+1)*B0(s)*((1-2*sqrt(abs(e))*r[i])*(exp(-2*sqrt(abs(e))*(r[i]-s))-exp(-2*sqrt(abs(e))*(r[i]+s)))/2+2*sqrt(abs(e))*s*(exp(-2*sqrt(abs(e))*(r[i]-s))+exp(-2*sqrt(abs(e))*(r[i]+s)))/2) + ,0,min(1,r[i]),weights)+ + (r[i]>=1 ? 0 : + rho/(2*(r[i]+1))*integrate_legendre(s->(s+r[i]+1)*B0(s+r[i])*((1-2*sqrt(abs(e))*s)*(1-exp(-4*sqrt(abs(e))*r[i]))*exp(-2*sqrt(abs(e))*s)+4*sqrt(abs(e))*r[i]*exp(-2*sqrt(abs(e))*(2*r[i]+s))),0,1-r[i],weights) + ) + end + return out +end +@everywhere function dseD1(e,rho,r,taus,T,weights,weights_gL,P,N,J,nu) + out=Array{Array{Float64}}(undef,J*N) + for i in 1:J*N + out[i]=Array{Float64}(undef,(P+1)*J) + for zeta in 0:J-1 + for n in 0:P + out[i][zeta*(P+1)+n+1]=rho/(r[i]+1)*integrate_legendre(s-> + (s+1)*B1(s,zeta,n,taus,T,weights,nu)*((1-2*sqrt(abs(e))*r[i])*(exp(-2*sqrt(abs(e))*(r[i]-s))-exp(-2*sqrt(abs(e))*(r[i]+s)))/2+2*sqrt(abs(e))*s*(exp(-2*sqrt(abs(e))*(r[i]-s))+exp(-2*sqrt(abs(e))*(r[i]+s)))/2) + ,0,r[i],weights)+ + rho/(2*(r[i]+1))*integrate_laguerre(s-> + (s+r[i]+1)*B1(s+r[i],zeta,n,taus,T,weights,nu)*((1-2*sqrt(abs(e))*s)*(1-exp(-4*sqrt(abs(e))*r[i]))+4*sqrt(abs(e))*r[i]*exp(-4*sqrt(abs(e))*r[i])) + ,2*sqrt(abs(e)),weights_gL) + end + end + end + return out +end +@everywhere function dseD2(e,rho,r,taus,T,weights,weights_gL,P,N,J,nu) + out=Array{Array{Float64,2}}(undef,J*N) + for i in 1:J*N + out[i]=Array{Float64,2}(undef,(P+1)*J,(P+1)*J) + 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]=rho/(r[i]+1)*integrate_legendre(s-> + (s+1)*B2(s,zeta,n,zetap,m,taus,T,weights,nu)*((1-2*sqrt(abs(e))*r[i])*(exp(-2*sqrt(abs(e))*(r[i]-s))-exp(-2*sqrt(abs(e))*(r[i]+s)))/2+2*sqrt(abs(e))*s*(exp(-2*sqrt(abs(e))*(r[i]-s))+exp(-2*sqrt(abs(e))*(r[i]+s)))/2) + ,0,r[i],weights)+ + rho/(2*(r[i]+1))*integrate_laguerre(s-> + (s+r[i]+1)*B2(s+r[i],zeta,n,zetap,m,taus,T,weights,nu)*((1-2*sqrt(abs(e))*s)*(1-exp(-4*sqrt(abs(e))*r[i]))+4*sqrt(abs(e))*r[i]*exp(-4*sqrt(abs(e))*r[i])) + ,2*sqrt(abs(e)),weights_gL) + end + end + end + end + end + return out +end + +# dD/d sqrt(abs(e+mu/2))'s +@everywhere function dsmuD0(e,rho,r,weights,N,J) + out=Array{Float64}(undef,J*N) + for i in 1:J*N + out[i]=-2*r[i]*exp(-2*sqrt(abs(e))*r[i])/(r[i]+1)+ + rho/(r[i]+1)*integrate_legendre(s-> + (s+1)*B0(s)*((-1-2*sqrt(abs(e))*r[i])*(exp(-2*sqrt(abs(e))*(r[i]-s))-exp(-2*sqrt(abs(e))*(r[i]+s)))/2+2*sqrt(abs(e))*s*(exp(-2*sqrt(abs(e))*(r[i]-s))+exp(-2*sqrt(abs(e))*(r[i]+s)))/2) + ,0,min(1,r[i]),weights)+ + (r[i]>=1 ? 0 : + rho/(2*(r[i]+1))*integrate_legendre(s->(s+r[i]+1)*B0(s+r[i])*((-1-2*sqrt(abs(e))*s)*(1-exp(-4*sqrt(abs(e))*r[i]))*exp(-2*sqrt(abs(e))*s)+4*sqrt(abs(e))*r[i]*exp(-2*sqrt(abs(e))*(2*r[i]+s))),0,1-r[i],weights) + ) + end + return out +end +@everywhere function dsmuD1(e,rho,r,taus,T,weights,weights_gL,P,N,J,nu) + out=Array{Array{Float64}}(undef,J*N) + for i in 1:J*N + out[i]=Array{Float64}(undef,(P+1)*J) + for zeta in 0:J-1 + for n in 0:P + out[i][zeta*(P+1)+n+1]=rho/(r[i]+1)*integrate_legendre(s-> + (s+1)*B1(s,zeta,n,taus,T,weights,nu)*((-1-2*sqrt(abs(e))*r[i])*(exp(-2*sqrt(abs(e))*(r[i]-s))-exp(-2*sqrt(abs(e))*(r[i]+s)))/2+2*sqrt(abs(e))*s*(exp(-2*sqrt(abs(e))*(r[i]-s))+exp(-2*sqrt(abs(e))*(r[i]+s)))/2) + ,0,r[i],weights)+ + rho/(2*(r[i]+1))*integrate_laguerre(s-> + (s+r[i]+1)*B1(s+r[i],zeta,n,taus,T,weights,nu)*((-1-2*sqrt(abs(e))*s)*(1-exp(-4*sqrt(abs(e))*r[i]))+4*sqrt(abs(e))*r[i]*exp(-4*sqrt(abs(e))*r[i])) + ,2*sqrt(abs(e)),weights_gL) + end + end + end + return out +end +@everywhere function dsmuD2(e,rho,r,taus,T,weights,weights_gL,P,N,J,nu) + out=Array{Array{Float64,2}}(undef,J*N) + for i in 1:J*N + out[i]=Array{Float64,2}(undef,(P+1)*J,(P+1)*J) + 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]=rho/(r[i]+1)*integrate_legendre(s-> + (s+1)*B2(s,zeta,n,zetap,m,taus,T,weights,nu)*((-1-2*sqrt(abs(e))*r[i])*(exp(-2*sqrt(abs(e))*(r[i]-s))-exp(-2*sqrt(abs(e))*(r[i]+s)))/2+2*sqrt(abs(e))*s*(exp(-2*sqrt(abs(e))*(r[i]-s))+exp(-2*sqrt(abs(e))*(r[i]+s)))/2) + ,0,r[i],weights)+ + rho/(2*(r[i]+1))*integrate_laguerre(s-> + (s+r[i]+1)*B2(s+r[i],zeta,n,zetap,m,taus,T,weights,nu)*((-1-2*sqrt(abs(e))*s)*(1-exp(-4*sqrt(abs(e))*r[i]))+4*sqrt(abs(e))*r[i]*exp(-4*sqrt(abs(e))*r[i])) + ,2*sqrt(abs(e)),weights_gL) + end + end + end + end + end + return out +end + +# gamma's +@everywhere function gamma0(e,rho,weights) + return 2*rho*e*integrate_legendre(s->(s+1)*B0(s)*exp(-2*sqrt(abs(e))*s),0,1,weights) +end +@everywhere function gamma1(e,rho,taus,T,weights,weights_gL,P,J,nu) + out=Array{Float64}(undef,J*(P+1)) + for zeta in 0:J-1 + for n in 0:P + out[zeta*(P+1)+n+1]=2*rho*e*integrate_laguerre(s->(s+1)*B1(s,zeta,n,taus,T,weights,nu),2*sqrt(abs(e)),weights_gL) + end + end + return out +end +@everywhere function gamma2(e,rho,taus,T,weights,weights_gL,P,J,nu) + out=Array{Float64,2}(undef,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[zeta*(P+1)+n+1,zetap*(P+1)+m+1]=2*rho*e*integrate_laguerre(s->(s+1)*B2(s,zeta,n,zetap,m,taus,T,weights,nu),2*sqrt(abs(e)),weights_gL) + end + end + end + end + return out +end + +# dgamma/d sqrt(abs(e))'s +@everywhere function dsedgamma0(e,rho,weights) + return 4*rho*sqrt(abs(e))*integrate_legendre(s->(s+1)*B0(s)*(1-sqrt(abs(e))*s)*exp(-2*sqrt(abs(e))*s),0,1,weights) +end +@everywhere function dsedgamma1(e,rho,taus,T,weights,weights_gL,P,J,nu) + out=Array{Float64}(undef,J*(P+1)) + for zeta in 0:J-1 + for n in 0:P + out[zeta*(P+1)+n+1]=4*rho*e*integrate_laguerre(s->(s+1)*B1(s,zeta,n,taus,T,weights,nu)*(1-sqrt(abs(e))*s),2*sqrt(abs(e)),weights_gL) + end + end + return out +end +@everywhere function dsedgamma2(e,rho,taus,T,weights,weights_gL,P,J,nu) + out=Array{Float64,2}(undef,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[zeta*(P+1)+n+1,zetap*(P+1)+m+1]=4*rho*e*integrate_laguerre(s->(s+1)*B2(s,zeta,n,zetap,m,taus,T,weights,nu)*(1-sqrt(abs(e))*s),2*sqrt(abs(e)),weights_gL) + end + end + end + end + return out +end + +# dgamma/d sqrt(e+mu/2)'s +@everywhere function dsmudgamma0(e,rho,weights) + return -4*rho*e*integrate_legendre(s->(s+1)*s*B0(s)*exp(-2*sqrt(abs(e))*s),0,1,weights) +end +@everywhere function dsmudgamma1(e,rho,taus,T,weights,weights_gL,P,J,nu) + out=Array{Float64}(undef,J*(P+1)) + for zeta in 0:J-1 + for n in 0:P + out[zeta*(P+1)+n+1]=-4*rho*e*integrate_laguerre(s->s*(s+1)*B1(s,zeta,n,taus,T,weights,nu),2*sqrt(abs(e)),weights_gL) + end + end + return out +end +@everywhere function dsmudgamma2(e,rho,taus,T,weights,weights_gL,P,J,nu) + out=Array{Float64,2}(undef,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[zeta*(P+1)+n+1,zetap*(P+1)+m+1]=-4*rho*e*integrate_laguerre(s->s*(s+1)*B2(s,zeta,n,zetap,m,taus,T,weights,nu),2*sqrt(abs(e)),weights_gL) + end + end + end + end + return out +end + +# \bar gamma's +@everywhere function gammabar0(weights) + return 4*pi*integrate_legendre(s->s^2*B0(s-1),0,1,weights) +end +@everywhere function gammabar1(taus,T,weights,P,J,nu) + out=Array{Float64}(undef,J*(P+1)) + for zeta in 0:J-1 + for n in 0:P + out[zeta*(P+1)+n+1]=4*pi*integrate_legendre(s->s^2*B1(s-1,zeta,n,taus,T,weights,nu),0,1,weights) + end + end + return out +end +@everywhere function gammabar2(taus,T,weights,P,J,nu) + out=Array{Float64,2}(undef,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[zeta*(P+1)+n+1,zetap*(P+1)+m+1]=4*pi*integrate_legendre(s->s^2*B2(s-1,zeta,n,zetap,m,taus,T,weights,nu),0,1,weights) + end + end + end + end + return out +end |