From c1b477a1b2b796617c4e345a7296a8d429d7a067 Mon Sep 17 00:00:00 2001 From: Ian Jauslin Date: Sun, 26 Feb 2023 18:36:05 -0500 Subject: Update to v0.4 feature: compute the 2-point correlation function in easyeq. feature: compute the Fourier transform of the 2-point correlation function in anyeq and easyeq. feature: compute the local maximum of the 2-point correlation function and its Fourier transform. feature: compute the compressibility for anyeq. feature: allow for linear spacing of rho's. feature: print the scattering length. change: ux and uk now return real numbers. fix: error in the computation of the momentum distribution: wrong definition of delta functions. fix: various minor bugs. optimization: assign explicit types to variables. --- src/simpleq-hardcore.jl | 494 ++++++++++++++++++++++++++++++++++++------------ 1 file changed, 377 insertions(+), 117 deletions(-) (limited to 'src/simpleq-hardcore.jl') diff --git a/src/simpleq-hardcore.jl b/src/simpleq-hardcore.jl index ca64f78..398ac06 100644 --- a/src/simpleq-hardcore.jl +++ b/src/simpleq-hardcore.jl @@ -1,4 +1,4 @@ -## Copyright 2021 Ian Jauslin +## Copyright 2021-2023 Ian Jauslin ## ## Licensed under the Apache License, Version 2.0 (the "License"); ## you may not use this file except in compliance with the License. @@ -13,28 +13,26 @@ ## 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 +function simpleq_hardcore_energy_rho( + rhos::Array{Float64,1}, + taus::Array{Float64,1}, + P::Int64, + N::Int64, + J::Int64, + maxiter::Int64, + tolerance::Float64 +) + # spawn workers + work=spawn_workers(length(rhos)) # 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)) + u0s=Array{Array{Float64,1}}(undef,length(rhos)) + e0s=Array{Float64,1}(undef,length(rhos)) for j in 1:length(rhos) - u0s[j]=Array{Float64}(undef,N*J) + u0s[j]=Array{Float64,1}(undef,N*J) for i in 1:N*J u0s[j][i]=1/(1+r[i]^2)^2 end @@ -43,17 +41,17 @@ function simpleq_hardcore_energy_rho(rhos,taus,P,N,J,maxiter,tolerance) # save result from each task - us=Array{Array{Float64}}(undef,length(rhos)) - es=Array{Float64}(undef,length(rhos)) - err=Array{Float64}(undef,length(rhos)) + us=Array{Array{Float64,1}}(undef,length(rhos)) + es=Array{Float64,1}(undef,length(rhos)) + err=Array{Float64,1}(undef,length(rhos)) count=0 # for each worker - @sync for p in 1:nw + @sync for p in 1:length(work) # for each task @async for j in work[p] count=count+1 - if count>=nw + if count>=length(work) progress(count,length(rhos),10000) end # run the task @@ -67,12 +65,20 @@ function simpleq_hardcore_energy_rho(rhos,taus,P,N,J,maxiter,tolerance) end # compute u(x) -function simpleq_hardcore_ux(rho,taus,P,N,J,maxiter,tolerance) +function simpleq_hardcore_ux( + rho::Float64, + taus::Array{Float64,1}, + P::Int64, + N::Int64, + J::Int64, + maxiter::Int64, + tolerance::Float64 +) # initialize vectors (weights,weights_gL,r,T)=simpleq_hardcore_init(taus,P,N,J) # initial guess - u0=Array{Float64}(undef,N*J) + u0=Array{Float64,1}(undef,N*J) for i in 1:N*J u0[i]=1/(1+r[i]^2)^2 end @@ -86,28 +92,26 @@ function simpleq_hardcore_ux(rho,taus,P,N,J,maxiter,tolerance) 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 +function simpleq_hardcore_condensate_fraction_rho( + rhos::Array{Float64,1}, + taus::Array{Float64,1}, + P::Int64, + N::Int64, + J::Int64, + maxiter::Int64, + tolerance::Float64 +) + # spawn workers + work=spawn_workers(length(rhos)) # 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)) + u0s=Array{Array{Float64,1}}(undef,length(rhos)) + e0s=Array{Float64,1}(undef,length(rhos)) for j in 1:length(rhos) - u0s[j]=Array{Float64}(undef,N*J) + u0s[j]=Array{Float64,1}(undef,N*J) for i in 1:N*J u0s[j][i]=1/(1+r[i]^2)^2 end @@ -116,17 +120,17 @@ function simpleq_hardcore_condensate_fraction_rho(rhos,taus,P,N,J,maxiter,tolera # save result from each task - us=Array{Array{Float64}}(undef,length(rhos)) - es=Array{Float64}(undef,length(rhos)) - err=Array{Float64}(undef,length(rhos)) + us=Array{Array{Float64,1}}(undef,length(rhos)) + es=Array{Float64,1}(undef,length(rhos)) + err=Array{Float64,1}(undef,length(rhos)) count=0 # for each worker - @sync for p in 1:nw + @sync for p in 1:length(work) # for each task @async for j in work[p] count=count+1 - if count>=nw + if count>=length(work) progress(count,length(rhos),10000) end # run the task @@ -142,13 +146,18 @@ end # initialize computation -@everywhere function simpleq_hardcore_init(taus,P,N,J) +@everywhere function simpleq_hardcore_init( + taus::Array{Float64,1}, + P::Int64, + N::Int64, + J::Int64 +) # Gauss-Legendre weights weights=gausslegendre(N) weights_gL=gausslaguerre(N) # r - r=Array{Float64}(undef,J*N) + r=Array{Float64,1}(undef,J*N) for zeta in 0:J-1 for j in 1:N xj=weights[1][j] @@ -164,9 +173,23 @@ end 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) +@everywhere function simpleq_hardcore_hatu( + u0::Array{Float64,1}, + e0::Float64, + rho::Float64, + r::Array{Float64,1}, + taus::Array{Float64,1}, + T::Array{Polynomial,1}, + weights::Tuple{Array{Float64,1},Array{Float64,1}}, + weights_gL::Tuple{Array{Float64,1},Array{Float64,1}}, + P::Int64, + N::Int64, + J::Int64, + maxiter::Int64, + tolerance::Float64 +) # init - vec=Array{Float64}(undef,J*N+1) + vec=Array{Float64,1}(undef,J*N+1) for i in 1:J*N vec[i]=u0[i] end @@ -194,8 +217,20 @@ end 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) +@everywhere function simpleq_hardcore_Xi( + u::Array{Float64,1}, + e::Float64, + rho::Float64, + r::Array{Float64,1}, + taus::Array{Float64,1}, + T::Array{Polynomial,1}, + weights::Tuple{Array{Float64,1},Array{Float64,1}}, + weights_gL::Tuple{Array{Float64,1},Array{Float64,1}}, + P::Int64, + N::Int64, + J::Int64 +) + out=Array{Float64,1}(undef,J*N+1) FU=chebyshev(u,taus,weights,P,N,J,4) #D's @@ -215,7 +250,19 @@ end return out end # DXi -@everywhere function simpleq_hardcore_DXi(u,e,rho,r,taus,T,weights,weights_gL,P,N,J) +@everywhere function simpleq_hardcore_DXi( + u::Array{Float64,1}, + e::Float64, + rho::Float64, + r::Array{Float64,1}, + taus::Array{Float64,1}, + T::Array{Polynomial,1}, + weights::Tuple{Array{Float64,1},Array{Float64,1}}, + weights_gL::Tuple{Array{Float64,1},Array{Float64,1}}, + P::Int64, + N::Int64, + J::Int64 +) out=Array{Float64,2}(undef,J*N+1,J*N+1) FU=chebyshev(u,taus,weights,P,N,J,4) @@ -232,15 +279,15 @@ end for zetapp in 0:J-1 for n in 1:N - one=zeros(Int64,J*N) - one[zetapp*N+n]=1 + one=zeros(Float64,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) + 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* @@ -253,14 +300,26 @@ 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) + (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) +@everywhere function simpleq_hardcore_dXidmu( + u::Array{Float64,1}, + e::Float64, + rho::Float64, + r::Array{Float64,1}, + taus::Array{Float64,1}, + T::Array{Polynomial,1}, + weights::Tuple{Array{Float64,1},Array{Float64,1}}, + weights_gL::Tuple{Array{Float64,1},Array{Float64,1}}, + P::Int64, + N::Int64, + J::Int64 +) + out=Array{Float64,1}(undef,J*N+1) FU=chebyshev(u,taus,weights,P,N,J,4) #D's @@ -281,17 +340,37 @@ end end # B's -@everywhere function B0(r) +@everywhere function B0( + r::Float64 +) return pi/12*(r-1)^2*(r+5) end -@everywhere function B1(r,zeta,n,taus,T,weights,nu) +@everywhere function B1( + r::Float64, + zeta::Int64, + n::Int64, + taus::Array{Float64,1}, + T::Array{Polynomial,1}, + weights::Tuple{Array{Float64,1},Array{Float64,1}}, + nu::Int64 +) 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) +@everywhere function B2( + r::Float64, + zeta::Int64, + n::Int64, + zetap::Int64, + m::Int64, + taus::Array{Float64,1}, + T::Array{Polynomial,1}, + weights::Tuple{Array{Float64,1},Array{Float64,1}}, + nu::Int64 +) 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 : @@ -303,30 +382,49 @@ end end # D's -@everywhere function D0(e,rho,r,weights,N,J) - out=Array{Float64}(undef,J*N) +@everywhere function D0( + e::Float64, + rho::Float64, + r::Array{Float64,1}, + weights::Tuple{Array{Float64,1},Array{Float64,1}}, + N::Int64, + J::Int64 +) + out=Array{Float64,1}(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)+ + ,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) + ,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) +@everywhere function D1( + e::Float64, + rho::Float64, + r::Array{Float64,1}, + taus::Array{Float64,1}, + T::Array{Polynomial,1}, + weights::Tuple{Array{Float64,1},Array{Float64,1}}, + weights_gL::Tuple{Array{Float64,1},Array{Float64,1}}, + P::Int64, + N::Int64, + J::Int64, + nu::Int64 +) + out=Array{Array{Float64,1},1}(undef,J*N) for i in 1:J*N - out[i]=Array{Float64}(undef,(P+1)*J) + out[i]=Array{Float64,1}(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)+ + ,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) @@ -336,7 +434,19 @@ end return out end -@everywhere function D2(e,rho,r,taus,T,weights,weights_gL,P,N,J,nu) +@everywhere function D2( + e::Float64, + rho::Float64, + r::Array{Float64,1}, + taus::Array{Float64,1}, + T::Array{Polynomial,1}, + weights::Tuple{Array{Float64,1},Array{Float64,1}}, + weights_gL::Tuple{Array{Float64,1},Array{Float64,1}}, + P::Int64, + N::Int64, + J::Int64, + nu::Int64 +) 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) @@ -346,7 +456,7 @@ end 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)+ + ,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) @@ -359,28 +469,47 @@ end end # dD/d sqrt(abs(e))'s -@everywhere function dseD0(e,rho,r,weights,N,J) - out=Array{Float64}(undef,J*N) +@everywhere function dseD0( + e::Float64, + rho::Float64, + r::Array{Float64,1}, + weights::Tuple{Array{Float64,1},Array{Float64,1}}, + N::Int64, + J::Int64 +) + out=Array{Float64,1}(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)+ + ,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) + 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) +@everywhere function dseD1( + e::Float64, + rho::Float64, + r::Array{Float64,1}, + taus::Array{Float64,1}, + T::Array{Polynomial,1}, + weights::Tuple{Array{Float64,1},Array{Float64,1}}, + weights_gL::Tuple{Array{Float64,1},Array{Float64,1}}, + P::Int64, + N::Int64, + J::Int64, + nu::Int64 +) + out=Array{Array{Float64,1},1}(undef,J*N) for i in 1:J*N - out[i]=Array{Float64}(undef,(P+1)*J) + out[i]=Array{Float64,1}(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)+ + ,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) @@ -389,7 +518,19 @@ end end return out end -@everywhere function dseD2(e,rho,r,taus,T,weights,weights_gL,P,N,J,nu) +@everywhere function dseD2( + e::Float64, + rho::Float64, + r::Array{Float64,1}, + taus::Array{Float64,1}, + T::Array{Polynomial,1}, + weights::Tuple{Array{Float64,1},Array{Float64,1}}, + weights_gL::Tuple{Array{Float64,1},Array{Float64,1}}, + P::Int64, + N::Int64, + J::Int64, + nu::Int64 +) 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) @@ -399,7 +540,7 @@ end 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)+ + ,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) @@ -412,28 +553,47 @@ end end # dD/d sqrt(abs(e+mu/2))'s -@everywhere function dsmuD0(e,rho,r,weights,N,J) - out=Array{Float64}(undef,J*N) +@everywhere function dsmuD0( + e::Float64, + rho::Float64, + r::Array{Float64,1}, + weights::Tuple{Array{Float64,1},Array{Float64,1}}, + N::Int64, + J::Int64 +) + out=Array{Float64,1}(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)+ + ,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) + 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) +@everywhere function dsmuD1( + e::Float64, + rho::Float64, + r::Array{Float64,1}, + taus::Array{Float64,1}, + T::Array{Polynomial,1}, + weights::Tuple{Array{Float64,1},Array{Float64,1}}, + weights_gL::Tuple{Array{Float64,1},Array{Float64,1}}, + P::Int64, + N::Int64, + J::Int64, + nu::Int64 +) + out=Array{Array{Float64,1},1}(undef,J*N) for i in 1:J*N - out[i]=Array{Float64}(undef,(P+1)*J) + out[i]=Array{Float64,1}(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)+ + ,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) @@ -442,7 +602,19 @@ end end return out end -@everywhere function dsmuD2(e,rho,r,taus,T,weights,weights_gL,P,N,J,nu) +@everywhere function dsmuD2( + e::Float64, + rho::Float64, + r::Array{Float64,1}, + taus::Array{Float64,1}, + T::Array{Polynomial,1}, + weights::Tuple{Array{Float64,1},Array{Float64,1}}, + weights_gL::Tuple{Array{Float64,1},Array{Float64,1}}, + P::Int64, + N::Int64, + J::Int64, + nu::Int64 +) 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) @@ -452,7 +624,7 @@ end 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)+ + ,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) @@ -465,11 +637,25 @@ end 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)) +@everywhere function gamma0( + e::Float64, + rho::Float64, + weights::Tuple{Array{Float64,1},Array{Float64,1}} +) + 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::Float64, + rho::Float64, + taus::Array{Float64,1}, + T::Array{Polynomial,1}, + weights::Tuple{Array{Float64,1},Array{Float64,1}}, + weights_gL::Tuple{Array{Float64,1},Array{Float64,1}}, + P::Int64, + J::Int64, + nu::Int64 +) + out=Array{Float64,1}(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) @@ -477,7 +663,17 @@ end end return out end -@everywhere function gamma2(e,rho,taus,T,weights,weights_gL,P,J,nu) +@everywhere function gamma2( + e::Float64, + rho::Float64, + taus::Array{Float64,1}, + T::Array{Polynomial,1}, + weights::Tuple{Array{Float64,1},Array{Float64,1}}, + weights_gL::Tuple{Array{Float64,1},Array{Float64,1}}, + P::Int64, + J::Int64, + nu::Int64 +) out=Array{Float64,2}(undef,J*(P+1),J*(P+1)) for zeta in 0:J-1 for n in 0:P @@ -492,11 +688,25 @@ end 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)) +@everywhere function dsedgamma0( + e::Float64, + rho::Float64, + weights::Tuple{Array{Float64,1},Array{Float64,1}} +) + 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::Float64, + rho::Float64, + taus::Array{Float64,1}, + T::Array{Polynomial,1}, + weights::Tuple{Array{Float64,1},Array{Float64,1}}, + weights_gL::Tuple{Array{Float64,1},Array{Float64,1}}, + P::Int64, + J::Int64, + nu::Int64 +) + out=Array{Float64,1}(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) @@ -504,7 +714,17 @@ end end return out end -@everywhere function dsedgamma2(e,rho,taus,T,weights,weights_gL,P,J,nu) +@everywhere function dsedgamma2( + e::Float64, + rho::Float64, + taus::Array{Float64,1}, + T::Array{Polynomial,1}, + weights::Tuple{Array{Float64,1},Array{Float64,1}}, + weights_gL::Tuple{Array{Float64,1},Array{Float64,1}}, + P::Int64, + J::Int64, + nu::Int64 +) out=Array{Float64,2}(undef,J*(P+1),J*(P+1)) for zeta in 0:J-1 for n in 0:P @@ -519,11 +739,25 @@ end 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)) +@everywhere function dsmudgamma0( + e::Float64, + rho::Float64, + weights::Tuple{Array{Float64,1},Array{Float64,1}} +) + 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::Float64, + rho::Float64, + taus::Array{Float64,1}, + T::Array{Polynomial,1}, + weights::Tuple{Array{Float64,1},Array{Float64,1}}, + weights_gL::Tuple{Array{Float64,1},Array{Float64,1}}, + P::Int64, + J::Int64, + nu::Int64 +) + out=Array{Float64,1}(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) @@ -531,7 +765,17 @@ end end return out end -@everywhere function dsmudgamma2(e,rho,taus,T,weights,weights_gL,P,J,nu) +@everywhere function dsmudgamma2( + e::Float64, + rho::Float64, + taus::Array{Float64,1}, + T::Array{Polynomial,1}, + weights::Tuple{Array{Float64,1},Array{Float64,1}}, + weights_gL::Tuple{Array{Float64,1},Array{Float64,1}}, + P::Int64, + J::Int64, + nu::Int64 +) out=Array{Float64,2}(undef,J*(P+1),J*(P+1)) for zeta in 0:J-1 for n in 0:P @@ -546,25 +790,41 @@ end 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)) +@everywhere function gammabar0( + weights::Tuple{Array{Float64,1},Array{Float64,1}} +) + return 4*pi*integrate_legendre(s->s^2*B0(s-1),0.,1.,weights) +end +@everywhere function gammabar1( + taus::Array{Float64,1}, + T::Array{Polynomial,1}, + weights::Tuple{Array{Float64,1},Array{Float64,1}}, + P::Int64, + J::Int64, + nu::Int64 +) + out=Array{Float64,1}(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) + 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) +@everywhere function gammabar2( + taus::Array{Float64,1}, + T::Array{Polynomial,1}, + weights::Tuple{Array{Float64,1},Array{Float64,1}}, + P::Int64, + J::Int64, + nu::Int64 +) 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) + 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 -- cgit v1.2.3-70-g09d2