Ian Jauslin

Ian Jauslin

Home Publications Talks Software Teaching Contact
summaryrefslogtreecommitdiff
path: root/src/simpleq-hardcore.jl
diff options
context:
space:
mode:
Diffstat (limited to 'src/simpleq-hardcore.jl')
-rw-r--r--src/simpleq-hardcore.jl494
1 files changed, 377 insertions, 117 deletions
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
generated by cgit v1.2.3-54-g00ecf (git 2.45.1) at 2024-05-18 17:53:54 +0000