diff options
Diffstat (limited to 'src')
| -rw-r--r-- | src/anyeq.jl | 1485 | ||||
| -rw-r--r-- | src/chebyshev.jl | 323 | ||||
| -rw-r--r-- | src/easyeq.jl | 1040 | ||||
| -rw-r--r-- | src/integration.jl | 50 | ||||
| -rw-r--r-- | src/interpolation.jl | 36 | ||||
| -rw-r--r-- | src/main.jl | 168 | ||||
| -rw-r--r-- | src/multithread.jl | 16 | ||||
| -rw-r--r-- | src/optimization.jl | 94 | ||||
| -rw-r--r-- | src/potentials.jl | 58 | ||||
| -rw-r--r-- | src/print.jl | 2 | ||||
| -rw-r--r-- | src/simpleq-Kv.jl | 68 | ||||
| -rw-r--r-- | src/simpleq-hardcore.jl | 494 | ||||
| -rw-r--r-- | src/simpleq-iteration.jl | 49 | ||||
| -rw-r--r-- | src/tools.jl | 22 |
14 files changed, 3202 insertions, 703 deletions
diff --git a/src/anyeq.jl b/src/anyeq.jl index 8459cb0..fee77d8 100644 --- a/src/anyeq.jl +++ b/src/anyeq.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. @@ -28,24 +28,26 @@ end # compute energy for a given rho -# use minlrho, nlrho to incrementally compute the solution to medeq (to initialize the Newton algorithm) -function anyeq_energy(rho,minlrho,nlrho,taus,P,N,J,a0,v,maxiter,tolerance,approx,savefile) +function anyeq_energy( + rho::Float64, + minlrho_init::Float64, + nlrho_init::Int64, + taus::Array{Float64,1}, + P::Int64, + N::Int64, + J::Int64, + a0::Float64, + v::Function, + maxiter::Int64, + tolerance::Float64, + approx::Anyeq_approx, + savefile::String +) # initialize vectors - (weights,T,k,V,V0,A,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v) - # init Abar - if savefile!="" - Abar=anyeq_read_Abar(savefile,P,N,J) - else - Abar=anyeq_Abar_multithread(k,weights,taus,T,P,N,J,approx) - end + (weights,T,k,V,V0,A,Abar,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v,approx,savefile) # compute initial guess from medeq - rhos=Array{Float64}(undef,nlrho) - for j in 0:nlrho-1 - rhos[j+1]=(nlrho==1 ? rho : 10^(minlrho+(log10(rho)-minlrho)/(nlrho-1)*j)) - end - u0s=anyeq_init_medeq(rhos,N,J,k,a0,v,maxiter,tolerance) - u0=u0s[nlrho] + u0=anyeq_init_medeq([rho],minlrho_init,nlrho_init,N,J,k,a0,v,maxiter,tolerance) (u,E,error)=anyeq_hatu(u0,P,N,J,rho,a0,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,v,maxiter,tolerance,approx) @@ -53,72 +55,60 @@ function anyeq_energy(rho,minlrho,nlrho,taus,P,N,J,a0,v,maxiter,tolerance,approx end # compute energy as a function of rho -function anyeq_energy_rho(rhos,taus,P,N,J,a0,v,maxiter,tolerance,approx,savefile) +function anyeq_energy_rho( + rhos::Array{Float64,1}, + minlrho_init::Float64, + nlrho_init::Int64, + taus::Array{Float64,1}, + P::Int64, + N::Int64, + J::Int64, + a0::Float64, + v::Function, + maxiter::Int64, + tolerance::Float64, + approx::Anyeq_approx, + savefile::String +) # initialize vectors - (weights,T,k,V,V0,A,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v) - - # init Abar - if savefile!="" - Abar=anyeq_read_Abar(savefile,P,N,J) - else - Abar=anyeq_Abar_multithread(k,weights,taus,T,P,N,J,approx) - end + (weights,T,k,V,V0,A,Abar,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v,approx,savefile) # compute initial guess from medeq - u0s=anyeq_init_medeq(rhos,N,J,k,a0,v,maxiter,tolerance) - - # save result from each task - es=Array{Float64,1}(undef,length(rhos)) - err=Array{Float64,1}(undef,length(rhos)) + u0s=anyeq_init_medeq(rhos,minlrho_init,nlrho_init,N,J,k,a0,v,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-1)%nw+1],j) - end + # spawn workers + work=spawn_workers(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 - (u,es[j],err[j])=remotecall_fetch(anyeq_hatu,workers()[p],u0s[j],P,N,J,rhos[j],a0,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,v,maxiter,tolerance,approx) - end - end + # compute u + (us,es,errs)=anyeq_hatu_rho_multithread(u0s,P,N,J,rhos,a0,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,v,maxiter,tolerance,approx,work) for j in 1:length(rhos) - @printf("% .15e % .15e % .15e\n",rhos[j],es[j],err[j]) + @printf("% .15e % .15e % .15e\n",rhos[j],es[j],errs[j]) end end # compute energy as a function of rho # initialize with previous rho -function anyeq_energy_rho_init_prevrho(rhos,taus,P,N,J,a0,v,maxiter,tolerance,approx,savefile) +function anyeq_energy_rho_init_prevrho( + rhos::Array{Float64,1}, + minlrho_init::Float64, + nlrho_init::Int64, + taus::Array{Float64,1}, + P::Int64, + N::Int64, + J::Int64, + a0::Float64, + v::Function, + maxiter::Int64, + tolerance::Float64, + approx::Anyeq_approx, + savefile::String +) # initialize vectors - (weights,T,k,V,V0,A,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v) - - # init Abar - if savefile!="" - Abar=anyeq_read_Abar(savefile,P,N,J) - else - Abar=anyeq_Abar_multithread(k,weights,taus,T,P,N,J,approx) - end + (weights,T,k,V,V0,A,Abar,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v,approx,savefile) # compute initial guess from medeq - u0s=anyeq_init_medeq([rhos[1]],N,J,k,a0,v,maxiter,tolerance) - u=u0s[1] + u=anyeq_init_medeq([rhos[1]],minlrho_init,nlrho_init,N,J,k,a0,v,maxiter,tolerance) for j in 1:length(rhos) progress(j,length(rhos),10000) @@ -134,20 +124,26 @@ function anyeq_energy_rho_init_prevrho(rhos,taus,P,N,J,a0,v,maxiter,tolerance,ap end # compute energy as a function of rho # initialize with next rho -function anyeq_energy_rho_init_nextrho(rhos,taus,P,N,J,a0,v,maxiter,tolerance,approx,savefile) +function anyeq_energy_rho_init_nextrho( + rhos::Array{Float64,1}, + minlrho_init::Float64, + nlrho_init::Int64, + taus::Array{Float64,1}, + P::Int64, + N::Int64, + J::Int64, + a0::Float64, + v::Function, + maxiter::Int64, + tolerance::Float64, + approx::Anyeq_approx, + savefile::String +) # initialize vectors - (weights,T,k,V,V0,A,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v) - - # init Abar - if savefile!="" - Abar=anyeq_read_Abar(savefile,P,N,J) - else - Abar=anyeq_Abar_multithread(k,weights,taus,T,P,N,J,approx) - end + (weights,T,k,V,V0,A,Abar,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v,approx,savefile) # compute initial guess from medeq - u0s=anyeq_init_medeq([rhos[length(rhos)]],N,J,k,a0,v,maxiter,tolerance) - u=u0s[1] + u=anyeq_init_medeq([rhos[length(rhos)]],minlrho_init,nlrho_init,N,J,k,a0,v,maxiter,tolerance) for j in 1:length(rhos) progress(j,length(rhos),10000) @@ -163,23 +159,26 @@ function anyeq_energy_rho_init_nextrho(rhos,taus,P,N,J,a0,v,maxiter,tolerance,ap end # compute u(k) -# use minlrho, nlrho to incrementally compute the solution to medeq (to initialize the Newton algorithm) -function anyeq_uk(minlrho,nlrho,taus,P,N,J,rho,a0,v,maxiter,tolerance,approx,savefile) +function anyeq_uk( + minlrho_init::Float64, + nlrho_init::Int64, + taus::Array{Float64,1}, + P::Int64, + N::Int64, + J::Int64, + rho::Float64, + a0::Float64, + v::Function, + maxiter::Int64, + tolerance::Float64, + approx::Anyeq_approx, + savefile::String +) # init vectors - (weights,T,k,V,V0,A,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v) - # init Abar - if savefile!="" - Abar=anyeq_read_Abar(savefile,P,N,J) - else - Abar=anyeq_Abar_multithread(k,weights,taus,T,P,N,J,approx) - end + (weights,T,k,V,V0,A,Abar,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v,approx,savefile) + # compute initial guess from medeq - rhos=Array{Float64}(undef,nlrho) - for j in 0:nlrho-1 - rhos[j+1]=(nlrho==1 ? rho : 10^(minlrho+(log10(rho)-minlrho)/(nlrho-1)*j)) - end - u0s=anyeq_init_medeq(rhos,N,J,k,a0,v,maxiter,tolerance) - u0=u0s[nlrho] + u0=anyeq_init_medeq([rho],minlrho_init,nlrho_init,N,J,k,a0,v,maxiter,tolerance) (u,E,error)=anyeq_hatu(u0,P,N,J,rho,a0,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,v,maxiter,tolerance,approx) @@ -192,58 +191,60 @@ function anyeq_uk(minlrho,nlrho,taus,P,N,J,rho,a0,v,maxiter,tolerance,approx,sav end # compute u(x) -# use minlrho, nlrho to incrementally compute the solution to medeq (to initialize the Newton algorithm) -function anyeq_ux(minlrho,nlrho,taus,P,N,J,rho,a0,v,maxiter,tolerance,xmin,xmax,nx,approx,savefile) +function anyeq_ux( + minlrho_init::Float64, + nlrho_init::Int64, + taus::Array{Float64,1}, + P::Int64, + N::Int64, + J::Int64, + rho::Float64, + a0::Float64, + v::Function, + maxiter::Int64, + tolerance::Float64, + xmin::Float64, + xmax::Float64, + nx::Int64, + approx::Anyeq_approx, + savefile::String +) # init vectors - (weights,T,k,V,V0,A,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v) - # init Abar - if savefile!="" - Abar=anyeq_read_Abar(savefile,P,N,J) - else - Abar=anyeq_Abar_multithread(k,weights,taus,T,P,N,J,approx) - end + (weights,T,k,V,V0,A,Abar,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v,approx,savefile) # compute initial guess from medeq - rhos=Array{Float64}(undef,nlrho) - for j in 0:nlrho-1 - rhos[j+1]=(nlrho==1 ? rho : 10^(minlrho+(log10(rho)-minlrho)/(nlrho-1)*j)) - end - u0s=anyeq_init_medeq(rhos,N,J,k,a0,v,maxiter,tolerance) - u0=u0s[nlrho] + u0=anyeq_init_medeq([rho],minlrho_init,nlrho_init,N,J,k,a0,v,maxiter,tolerance) (u,E,error)=anyeq_hatu(u0,P,N,J,rho,a0,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,v,maxiter,tolerance,approx) for i in 1:nx x=xmin+(xmax-xmin)*i/nx - ux=0. - for zeta in 0:J-1 - for j in 1:N - ux+=(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 - @printf("% .15e % .15e % .15e\n",x,real(ux),imag(ux)) + ux=anyeq_u_x(x,u,k,N,J,taus,weights) + @printf("% .15e % .15e % .15e\n",x,ux,) end end # compute condensate fraction for a given rho -# use minlrho, nlrho to incrementally compute the solution to medeq (to initialize the Newton algorithm) -function anyeq_condensate_fraction(rho,minlrho,nlrho,taus,P,N,J,a0,v,maxiter,tolerance,approx,savefile) +function anyeq_condensate_fraction( + rho::Float64, + minlrho_init::Float64, + nlrho_init::Int64, + taus::Array{Float64,1}, + P::Int64, + N::Int64, + J::Int64, + a0::Float64, + v::Function, + maxiter::Int64, + tolerance::Float64, + approx::Anyeq_approx, + savefile::String +) # initialize vectors - (weights,T,k,V,V0,A,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v) - # init Abar - if savefile!="" - Abar=anyeq_read_Abar(savefile,P,N,J) - else - Abar=anyeq_Abar_multithread(k,weights,taus,T,P,N,J,approx) - end + (weights,T,k,V,V0,A,Abar,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v,approx,savefile) # compute initial guess from medeq - rhos=Array{Float64}(undef,nlrho) - for j in 0:nlrho-1 - rhos[j+1]=(nlrho==1 ? rho : 10^(minlrho+(log10(rho)-minlrho)/(nlrho-1)*j)) - end - u0s=anyeq_init_medeq(rhos,N,J,k,a0,v,maxiter,tolerance) - u0=u0s[nlrho] + u0=anyeq_init_medeq([rho],minlrho_init,nlrho_init,N,J,k,a0,v,maxiter,tolerance) (u,E,error)=anyeq_hatu(u0,P,N,J,rho,a0,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,v,maxiter,tolerance,approx) @@ -254,161 +255,571 @@ function anyeq_condensate_fraction(rho,minlrho,nlrho,taus,P,N,J,a0,v,maxiter,tol end # condensate fraction as a function of rho -function anyeq_condensate_fraction_rho(rhos,taus,P,N,J,a0,v,maxiter,tolerance,approx,savefile) - ## 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) +function anyeq_condensate_fraction_rho( + rhos::Array{Float64,1}, + minlrho_init::Float64, + nlrho_init::Int64, + taus::Array{Float64,1}, + P::Int64, + N::Int64, + J::Int64, + a0::Float64, + v::Function, + maxiter::Int64, + tolerance::Float64, + approx::Anyeq_approx, + savefile::String +) + # spawn workers + work=spawn_workers(length(rhos)) + + # initialize vectors + (weights,T,k,V,V0,A,Abar,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v,approx,savefile) + + # compute initial guess from medeq + u0s=anyeq_init_medeq(rhos,minlrho_init,nlrho_init,N,J,k,a0,v,maxiter,tolerance) + + # compute u + (us,es,errs)=anyeq_hatu_rho_multithread(u0s,P,N,J,rhos,a0,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,v,maxiter,tolerance,approx,work) + + # compute eta + etas=Array{Float64,1}(undef,length(rhos)) + count=0 + # for each worker + @sync for p in 1:length(work) + # for each task + @async for j in work[p] + count=count+1 + if count>=length(work) + progress(count,length(rhos),10000) + end + # run the task + etas[j]=remotecall_fetch(anyeq_eta,workers()[p],us[j],P,N,J,rhos[j],weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,approx) + end end + for j in 1:length(rhos) - append!(work[(j-1)%nw+1],j) + @printf("% .15e % .15e % .15e\n",rhos[j],etas[j],errs[j]) end +end +# compute the momentum distribution for a given rho +function anyeq_momentum_distribution( + kmin::Float64, + kmax::Float64, + rho::Float64, + minlrho_init::Float64, + nlrho_init::Int64, + taus::Array{Float64,1}, + P::Int64, + N::Int64, + J::Int64, + windowL::Float64, + a0::Float64, + v::Function, + maxiter::Int64, + tolerance::Float64, + approx::Anyeq_approx, + savefile::String +) # initialize vectors - (weights,T,k,V,V0,A,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v) - # init Abar - if savefile!="" - Abar=anyeq_read_Abar(savefile,P,N,J) + (weights,T,k,V,V0,A,Abar,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v,approx,savefile) + + # compute initial guess from medeq + u0=anyeq_init_medeq([rho],minlrho_init,nlrho_init,N,J,k,a0,v,maxiter,tolerance) + + (u,E,error)=anyeq_hatu(u0,P,N,J,rho,a0,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,v,maxiter,tolerance,approx) + + # compute M + if windowL==Inf + # use discrete approximation of delta function + M=anyeq_momentum_discrete_delta(kmin,kmax,u,P,N,J,rho,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,approx) else - Abar=anyeq_Abar_multithread(k,weights,taus,T,P,N,J,approx) + M=anyeq_momentum_window(kmin,kmax,u,P,N,J,windowL,rho,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,approx) end + # order k's in increasing order + for zeta in 0:J-1 + for j in 1:N + q=k[(J-1-zeta)*N+j] + # drop if not in k interval + if q<kmin || q>kmax + continue + end + @printf("% .15e % .15e\n",q,M[(J-1-zeta)*N+j]) + end + end +end + +# 2 point correlation function +function anyeq_2pt_correlation( + minlrho_init::Float64, + nlrho_init::Int64, + taus::Array{Float64,1}, + P::Int64, + N::Int64, + J::Int64, + windowL::Float64, + rho::Float64, + a0::Float64, + v::Function, + maxiter::Int64, + tolerance::Float64, + xmin::Float64, + xmax::Float64, + nx::Int64, + approx::Anyeq_approx, + savefile::String +) + # init vectors + (weights,T,k,V,V0,A,Abar,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v,approx,savefile) + # compute initial guess from medeq - u0s=anyeq_init_medeq(rhos,N,J,k,a0,v,maxiter,tolerance) + u0=anyeq_init_medeq([rho],minlrho_init,nlrho_init,N,J,k,a0,v,maxiter,tolerance) + + # compute u and some useful integrals + (u,E,error)=anyeq_hatu(u0,P,N,J,rho,a0,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,v,maxiter,tolerance,approx) + (S,E,II,JJ,X,Y,sL1,sK,G)=anyeq_SEIJGXY(rho*u,rho,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,weights,P,N,J,approx) + + for i in 1:nx + x=xmin+(xmax-xmin)*i/nx + C2=anyeq_2pt(x,u,P,N,J,windowL,rho,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,approx,S,E,II,JJ,X,Y,sL1,sK) + @printf("% .15e % .15e\n",x,C2) + end +end + +# maximum of 2 point correlation function +function anyeq_2pt_correlation_max( + rho::Float64, + minlrho_init::Float64, + nlrho_init::Int64, + dx::Float64, + x0::Float64, # initial guess is x0/rho^(1/3) + maxstep::Float64, + taus::Array{Float64,1}, + P::Int64, + N::Int64, + J::Int64, + windowL::Float64, + a0::Float64, + v::Function, + maxiter::Int64, + tolerance::Float64, + tolerance_max::Float64, + approx::Anyeq_approx, + savefile::String +) + # init vectors + (weights,T,k,V,V0,A,Abar,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v,approx,savefile) + + # compute initial guess from medeq + u0=anyeq_init_medeq([rho],minlrho_init,nlrho_init,N,J,k,a0,v,maxiter,tolerance) + + # initial guess + x0=1/rho^(1/3) + + (u,E,error)=anyeq_hatu(u0,P,N,J,rho,a0,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,v,maxiter,tolerance,approx) + (x,f)=anyeq_2pt_max(u,P,N,J,x0/rho^(1/3),dx,maxstep,maxiter,tolerance_max,windowL,rho,weights,T,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,approx) + + if(x==Inf) + @printf(stderr,"max search failed for rho=%e\n",rho) + else + @printf("% .15e % .15e\n",x,f) + end +end + +# maximum of 2 point correlation function as a function of rho +function anyeq_2pt_correlation_max_rho( + rhos::Array{Float64,1}, + minlrho_init::Float64, + nlrho_init::Int64, + dx::Float64, + x0::Float64, # initial guess is x0/rho^(1/3) + maxstep::Float64, + taus::Array{Float64,1}, + P::Int64, + N::Int64, + J::Int64, + windowL::Float64, + a0::Float64, + v::Function, + maxiter::Int64, + tolerance::Float64, + tolerance_max::Float64, + approx::Anyeq_approx, + savefile::String +) + # init vectors + (weights,T,k,V,V0,A,Abar,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v,approx,savefile) + + # compute initial guess from medeq + u0s=anyeq_init_medeq(rhos,minlrho_init,nlrho_init,N,J,k,a0,v,maxiter,tolerance) + + # save result from each task + xs=Array{Float64,1}(undef,length(rhos)) + fs=Array{Float64,1}(undef,length(rhos)) + + # spawn workers + work=spawn_workers(length(rhos)) # compute u - us=Array{Array{Float64,1}}(undef,length(rhos)) - errs=Array{Float64,1}(undef,length(rhos)) + (us,es,errs)=anyeq_hatu_rho_multithread(u0s,P,N,J,rhos,a0,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,v,maxiter,tolerance,approx,work) + 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 - progress(count,length(rhos),10000) + if count>=length(work) + progress(count,length(rhos),10000) end # run the task - (us[j],E,errs[j])=remotecall_fetch(anyeq_hatu,workers()[p],u0s[j],P,N,J,rhos[j],a0,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,v,maxiter,tolerance,approx) + (xs[j],fs[j])=remotecall_fetch(anyeq_2pt_max,workers()[p],us[j],P,N,J,x0/rhos[j]^(1/3),dx,maxstep,maxiter,tolerance_max,windowL,rhos[j],weights,T,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,approx) end end - # compute eta - etas=Array{Float64}(undef,length(rhos)) + for j in 1:length(rhos) + if(xs[j]==Inf) + @printf(stderr,"max search failed for rho=%e\n",rhos[j]) + else + @printf("% .15e % .15e % .15e\n",rhos[j],xs[j],fs[j]) + end + end +end + +# Correlation function in Fourier space +function anyeq_2pt_correlation_fourier( + minlrho_init::Float64, + nlrho_init::Int64, + taus::Array{Float64,1}, + P::Int64, + N::Int64, + J::Int64, + windowL::Float64, + rho::Float64, + a0::Float64, + v::Function, + maxiter::Int64, + tolerance::Float64, + kmin::Float64, + kmax::Float64, + nk::Int64, + approx::Anyeq_approx, + savefile::String +) + # init vectors + (weights,T,k,V,V0,A,Abar,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v,approx,savefile) + + # compute initial guess from medeq + u0=anyeq_init_medeq([rho],minlrho_init,nlrho_init,N,J,k,a0,v,maxiter,tolerance) + + # compute u and some useful integrals + (u,E,error)=anyeq_hatu(u0,P,N,J,rho,a0,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,v,maxiter,tolerance,approx) + (S,E,II,JJ,X,Y,sL1,sK,G)=anyeq_SEIJGXY(rho*u,rho,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,weights,P,N,J,approx) + + # save result from each task + C2s=Array{Float64,1}(undef,nk) + + # spawn workers + work=spawn_workers(nk) + 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 - progress(count,length(rhos),10000) + if count>=length(work) + progress(count,nk,10000) end # run the task - etas[j]=anyeq_eta(us[j],P,N,J,rhos[j],weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,approx) + C2s[j]=remotecall_fetch(anyeq_2pt_fourier,workers()[p],kmin+(kmax-kmin)*j/nk,u,P,N,J,windowL,rho,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,approx,S,E,II,JJ,X,Y,sL1,sK) end end - for j in 1:length(rhos) - @printf("% .15e % .15e % .15e\n",rhos[j],etas[j],errs[j]) + for j in 1:nk + @printf("% .15e % .15e\n",kmin+(kmax-kmin)*j/nk,C2s[j]) end end -# compute the momentum distribution for a given rho -# use minlrho, nlrho to incrementally compute the solution to medeq (to initialize the Newton algorithm) -function anyeq_momentum_distribution(rho,minlrho,nlrho,taus,P,N,J,a0,v,maxiter,tolerance,approx,savefile) - # initialize vectors - (weights,T,k,V,V0,A,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v) - # init Abar - if savefile!="" - Abar=anyeq_read_Abar(savefile,P,N,J) +# maximum of Fourier transform of 2 point correlation function +function anyeq_2pt_correlation_fourier_max( + rho::Float64, + minlrho_init::Float64, + nlrho_init::Int64, + dk::Float64, + k0::Float64, # initial guess is k0*rho^(1/3) + maxstep::Float64, + taus::Array{Float64,1}, + P::Int64, + N::Int64, + J::Int64, + windowL::Float64, + a0::Float64, + v::Function, + maxiter::Int64, + tolerance::Float64, + tolerance_max::Float64, + approx::Anyeq_approx, + savefile::String +) + # init vectors + (weights,T,k,V,V0,A,Abar,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v,approx,savefile) + + # compute initial guess from medeq + u0=anyeq_init_medeq([rho],minlrho_init,nlrho_init,N,J,k,a0,v,maxiter,tolerance) + + # compute u + (u,E,error)=anyeq_hatu(u0,P,N,J,rho,a0,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,v,maxiter,tolerance,approx) + (ko,f)=anyeq_2pt_fourier_max(u,P,N,J,k0*rho^(1/3),dk,maxstep,maxiter,tolerance_max,windowL,rho,weights,T,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,approx) + + if(ko==Inf) + @printf(stderr,"max search failed for rho=%e\n",rho) else - Abar=anyeq_Abar_multithread(k,weights,taus,T,P,N,J,approx) + @printf("% .15e % .15e\n",ko,f) end +end + +# maximum of fourier transform of 2 point correlation function as a function of rho +function anyeq_2pt_correlation_fourier_max_rho( + rhos::Array{Float64,1}, + minlrho_init::Float64, + nlrho_init::Int64, + dk::Float64, + k0::Float64, # initial guess is k0*rho^(1/3) + maxstep::Float64, + taus::Array{Float64,1}, + P::Int64, + N::Int64, + J::Int64, + windowL::Float64, + a0::Float64, + v::Function, + maxiter::Int64, + tolerance::Float64, + tolerance_max::Float64, + approx::Anyeq_approx, + savefile::String +) + # init vectors + (weights,T,k,V,V0,A,Abar,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v,approx,savefile) # compute initial guess from medeq - rhos=Array{Float64}(undef,nlrho) - for j in 0:nlrho-1 - rhos[j+1]=(nlrho==1 ? rho : 10^(minlrho+(log10(rho)-minlrho)/(nlrho-1)*j)) - end - u0s=anyeq_init_medeq(rhos,N,J,k,a0,v,maxiter,tolerance) - u0=u0s[nlrho] + u0s=anyeq_init_medeq(rhos,minlrho_init,nlrho_init,N,J,k,a0,v,maxiter,tolerance) - (u,E,error)=anyeq_hatu(u0,P,N,J,rho,a0,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,v,maxiter,tolerance,approx) + # spawn workers + work=spawn_workers(length(rhos)) - # compute M - M=anyeq_momentum(u,P,N,J,rho,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,approx) + # compute u + (us,es,errs)=anyeq_hatu_rho_multithread(u0s,P,N,J,rhos,a0,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,v,maxiter,tolerance,approx,work) - for zeta in 0:J-1 - for j in 1:N - # order k's in increasing order - @printf("% .15e % .15e\n",k[(J-1-zeta)*N+j],M[(J-1-zeta)*N+j]) + # save result from each task + ks=Array{Float64,1}(undef,length(rhos)) + fs=Array{Float64,1}(undef,length(rhos)) + + count=0 + # for each worker + @sync for p in 1:length(work) + # for each task + @async for j in work[p] + count=count+1 + if count>=length(work) + progress(count,length(rhos),10000) + end + # run the task + (ks[j],fs[j])=remotecall_fetch(anyeq_2pt_fourier_max,workers()[p],us[j],P,N,J,k0*rhos[j]^(1/3),dk,maxstep,maxiter,tolerance_max,windowL,rhos[j],weights,T,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,approx) + end + end + + for j in 1:length(rhos) + if(ks[j]==Inf) + @printf(stderr,"max search failed for rho=%e\n",rhos[j]) + else + @printf("% .15e % .15e % .15e\n",rhos[j],ks[j],fs[j]) end end end -# 2 point correlation function -function anyeq_2pt_correlation(minlrho,nlrho,taus,P,N,J,windowL,rho,a0,v,maxiter,tolerance,xmin,xmax,nx,approx,savefile) +# uncondensed 2 point correlation function +function anyeq_uncondensed_2pt_correlation( + minlrho_init::Float64, + nlrho_init::Int64, + taus::Array{Float64,1}, + P::Int64, + N::Int64, + J::Int64, + windowL::Float64, + rho::Float64, + a0::Float64, + v::Function, + maxiter::Int64, + tolerance::Float64, + xmin::Float64, + xmax::Float64, + nx::Int64, + approx::Anyeq_approx, + savefile::String +) # init vectors - (weights,T,k,V,V0,A,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v) - # init Abar - if savefile!="" - Abar=anyeq_read_Abar(savefile,P,N,J) - else - Abar=anyeq_Abar_multithread(k,weights,taus,T,P,N,J,approx) + (weights,T,k,V,V0,A,Abar,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v,approx,savefile) + + # compute initial guess from medeq + u0=anyeq_init_medeq([rho],minlrho_init,nlrho_init,N,J,k,a0,v,maxiter,tolerance) + + # compute u and some useful integrals + (u,E,error)=anyeq_hatu(u0,P,N,J,rho,a0,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,v,maxiter,tolerance,approx) + (S,E,II,JJ,X,Y,sL1,sK,G)=anyeq_SEIJGXY(rho*u,rho,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,weights,P,N,J,approx) + + # compute u(xi) + uxis=Array{Float64,1}(undef,nx) + for i in 1:nx + uxis[i]=anyeq_u_x(xmin+(xmax-xmin)*i/nx,u,k,N,J,taus,weights) + end + + + # spawn workers + work=spawn_workers(nx) + gamma2s=Array{Float64,1}(undef,nx) + count=0 + # for each worker + @sync for p in 1:length(work) + # for each task + @async for j in work[p] + count=count+1 + if count>=length(work) + progress(count,nx,10000) + end + # run the task + gamma2s[j]=remotecall_fetch(anyeq_uncondensed_2pt,workers()[p],xmin+(xmax-xmin)*j/nx,uxis[j],u,P,N,J,windowL,rho,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,approx,S,E,II,JJ,X,Y,sL1,sK) + end end + for i in 1:nx + @printf("% .15e % .15e\n",xmin+(xmax-xmin)*i/nx,gamma2s[i]) + end +end + +# compute the compressibility +function anyeq_compressibility_rho( + rhos::Array{Float64,1}, + minlrho_init::Float64, + nlrho_init::Int64, + taus::Array{Float64,1}, + P::Int64, + N::Int64, + J::Int64, + a0::Float64, + v::Function, + maxiter::Int64, + tolerance::Float64, + approx::Anyeq_approx, + savefile::String +) + # initialize vectors + (weights,T,k,V,V0,A,Abar,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v,approx,savefile) + # compute initial guess from medeq - rhos=Array{Float64}(undef,nlrho) - for j in 0:nlrho-1 - rhos[j+1]=(nlrho==1 ? rho : 10^(minlrho+(log10(rho)-minlrho)/(nlrho-1)*j)) + u0s=anyeq_init_medeq(rhos,minlrho_init,nlrho_init,N,J,k,a0,v,maxiter,tolerance) + + # save result from each task + es=Array{Float64,1}(undef,length(rhos)) + err=Array{Float64,1}(undef,length(rhos)) + + # spawn workers + work=spawn_workers(length(rhos)) + + # compute u + (us,es,errs)=anyeq_hatu_rho_multithread(u0s,P,N,J,rhos,a0,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,v,maxiter,tolerance,approx,work) + + for j in 1:length(rhos)-2 + # compute \rho^2\partial^2(\rho e)=\partial^2_{\log\rho}(\rho e)-\partial_{\log\rho}(\rho e) as a function of rho + # \partial^2_{\log\rho}(\rho e) + p2=(rhos[j+2]*es[j+2]-2*rhos[j+1]*es[j+1]+rhos[j]*es[j])/(log(rhos[j+2])-log(rhos[j+1]))/(log(rhos[j+1])-log(rhos[j])) + # \partial_{\log\rho}(\rho e) (take average of front and back) + p11=(rhos[j+2]*es[j+2]-rhos[j+1]*es[j+1])/(log(rhos[j+2])-log(rhos[j+1])) + p12=(rhos[j+1]*es[j+1]-rhos[j]*es[j])/(log(rhos[j+1])-log(rhos[j])) + # compressibility is 1/this + @printf("% .15e % .15e\n",rhos[j+1],1/(p2-(p11+p12)/2)) end - u0s=anyeq_init_medeq(rhos,N,J,k,a0,v,maxiter,tolerance) - u0=u0s[nlrho] +end + +# Fourier transform 2 point correlation function test: compute by transforming anyeq_2pt +function anyeq_2pt_correlation_fourier_test( + minlrho_init::Float64, + nlrho_init::Int64, + taus::Array{Float64,1}, + P::Int64, + N::Int64, + J::Int64, + windowL::Float64, + rho::Float64, + a0::Float64, + v::Function, + maxiter::Int64, + tolerance::Float64, + xmax::Float64, + kmin::Float64, + kmax::Float64, + nk::Int64, + approx::Anyeq_approx, + savefile::String +) + # init vectors + (weights,T,k,V,V0,A,Abar,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v,approx,savefile) + + # compute initial guess from medeq + u0=anyeq_init_medeq([rho],minlrho_init,nlrho_init,N,J,k,a0,v,maxiter,tolerance) # compute u and some useful integrals (u,E,error)=anyeq_hatu(u0,P,N,J,rho,a0,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,v,maxiter,tolerance,approx) (S,E,II,JJ,X,Y,sL1,sK,G)=anyeq_SEIJGXY(rho*u,rho,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,weights,P,N,J,approx) - for i in 1:nx - x=xmin+(xmax-xmin)*i/nx - C2=anyeq_2pt(x,u,P,N,J,windowL,rho,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,approx,S,E,II,JJ,X,Y,sL1,sK,G) - @printf("% .15e % .15e\n",x,C2) + # compute C2 in real space + C2=Array{Float64,1}(undef,N*J) + for i in 1:N*J + if k[i]<xmax + C2[i]=anyeq_2pt(k[i],u,P,N,J,windowL,rho,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,approx,S,E,II,JJ,X,Y,sL1,sK)-rho^2 + else + C2[i]=0. + end + end + # invert Fourier transform + for i in 1:nk + k0=kmin+(kmax-kmin)*i/nk + hatC2=inverse_fourier_chebyshev(C2,k0,k,taus,weights,N,J) + @printf("% .15e % .15e\n",k0,hatC2) end end # compute Abar -function anyeq_Abar_multithread(k,weights,taus,T,P,N,J,approx) +function anyeq_Abar_multithread( + k::Array{Float64,1}, + weights::Tuple{Array{Float64,1},Array{Float64,1}}, + taus::Array{Float64,1}, + T::Array{Polynomial,1}, + P::Int64, + N::Int64, + J::Int64, + approx::Anyeq_approx +) if approx.bL3==0. - return [] + # empty array + return Array{Array{Float64,5}}(undef,0) end out=Array{Array{Float64,5}}(undef,J*N) - ## 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:N*J - append!(work[(j-1)%nw+1],j) - end + # spawn workers + work=spawn_workers(N*J) 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,N*J,10000) end # run the task @@ -419,14 +830,23 @@ function anyeq_Abar_multithread(k,weights,taus,T,P,N,J,approx) return out end + # initialize computation -@everywhere function anyeq_init(taus,P,N,J,v) +@everywhere function anyeq_init( + taus::Array{Float64,1}, + P::Int64, + N::Int64, + J::Int64, + v::Function, + approx::Anyeq_approx, + savefile::String +) # Gauss-Legendre weights weights=gausslegendre(N) # initialize vectors V,k - V=Array{Float64}(undef,J*N) - k=Array{Float64}(undef,J*N) + V=Array{Float64,1}(undef,J*N) + k=Array{Float64,1}(undef,J*N) for zeta in 0:J-1 for j in 1:N xj=weights[1][j] @@ -437,7 +857,7 @@ end end end # potential at 0 - V0=v(0) + V0=v(0.) # initialize matrix A T=chebyshev_polynomials(P) @@ -449,39 +869,62 @@ end Upsilon=Upsilonmat(k,v,weights_plus) Upsilon0=Upsilon0mat(k,v,weights_plus) - return(weights,T,k,V,V0,A,Upsilon,Upsilon0) + # init Abar + if savefile!="" + Abar=anyeq_read_Abar(savefile,P,N,J) + else + Abar=anyeq_Abar_multithread(k,weights,taus,T,P,N,J,approx) + end + + return(weights,T,k,V,V0,A,Abar,Upsilon,Upsilon0) end # compute initial guess from medeq -@everywhere function anyeq_init_medeq(rhos,N,J,k,a0,v,maxiter,tolerance) - us_medeq=Array{Array{Float64,1}}(undef,length(rhos)) - u0s=Array{Array{Float64,1}}(undef,length(rhos)) - +@everywhere function anyeq_init_medeq( + rhos::Array{Float64,1}, + minlrho_init::Float64, + nlrho_init::Int64, + N::Int64, + J::Int64, + k::Array{Float64,1}, + a0::Float64, + v::Function, + maxiter::Int64, + tolerance::Float64 +) weights_medeq=gausslegendre(N*J) + (u0s_medeq,es,errs)=easyeq_compute_u_prevrho(rhos,minlrho_init,nlrho_init,a0,N*J,v,maxiter,tolerance,weights_medeq,Easyeq_approx(1.,1.)) - (us_medeq[1],E,err)=easyeq_hatu(easyeq_init_u(a0,J*N,weights_medeq),J*N,rhos[1],v,maxiter,tolerance,weights_medeq,Easyeq_approx(1.,1.)) - u0s[1]=easyeq_to_anyeq(us_medeq[1],weights_medeq,k,N,J) - if err>tolerance - print(stderr,"warning: computation of initial Ansatz failed for rho=",rhos[1],"\n") - end - - for j in 2:length(rhos) - (us_medeq[j],E,err)=easyeq_hatu(us_medeq[j-1],J*N,rhos[j],v,maxiter,tolerance,weights_medeq,Easyeq_approx(1.,1.)) - u0s[j]=easyeq_to_anyeq(us_medeq[j],weights_medeq,k,N,J) - - if err>tolerance + # check errs + for j in 1:length(errs) + if errs[j]>tolerance print(stderr,"warning: computation of initial Ansatz failed for rho=",rhos[j],"\n") end end - return u0s + # return a single vector if there is a single rho + if length(rhos)>1 + u0s=Array{Array{Float64,1}}(undef,length(rhos)) + for j in 1:length(u0s_medeq) + u0s[j]=easyeq_to_anyeq(u0s_medeq[j],weights_medeq,k,N,J) + end + return u0s + else + return easyeq_to_anyeq(u0s_medeq,weights_medeq,k,N,J) + end end # interpolate the solution of medeq to an input for anyeq -@everywhere function easyeq_to_anyeq(u_simple,weights,k,N,J) +@everywhere function easyeq_to_anyeq( + u_simple::Array{Float64,1}, + weights::Tuple{Array{Float64,1},Array{Float64,1}}, + k::Array{Float64,1}, + N::Int64, + J::Int64 +) # reorder u_simple, which is evaluated at (1-x_j)/(1+x_j) with x_j\in[-1,1] u_s=zeros(Float64,length(u_simple)) - k_s=Array{Float64}(undef,length(u_simple)) + k_s=Array{Float64,1}(undef,length(u_simple)) for j in 1:length(u_simple) xj=weights[1][j] k_s[length(u_simple)-j+1]=(1-xj)/(1+xj) @@ -501,7 +944,27 @@ end # compute u using chebyshev expansions -@everywhere function anyeq_hatu(u0,P,N,J,rho,a0,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,v,maxiter,tolerance,approx) +@everywhere function anyeq_hatu( + u0::Array{Float64,1}, + P::Int64, + N::Int64, + J::Int64, + rho::Float64, + a0::Float64, + weights::Tuple{Array{Float64,1},Array{Float64,1}}, + k::Array{Float64,1}, + taus::Array{Float64,1}, + V::Array{Float64,1}, + V0::Float64, + A::Array{Array{Float64,2},1}, + Abar::Array{Array{Float64,5},1}, + Upsilon::Array{Array{Float64,1},1}, + Upsilon0::Array{Float64,1}, + v::Function, + maxiter::Int64, + tolerance::Float64, + approx::Anyeq_approx +) # init # rescale by rho (that's how u is defined) u=rho*u0 @@ -512,7 +975,7 @@ end # run Newton algorithm for i in 1:maxiter-1 (S,E,II,JJ,X,Y,sL1,sK,G)=anyeq_SEIJGXY(u,rho,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,weights,P,N,J,approx) - new=u-inv(anyeq_DXi(u,rho,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,weights,P,N,J,approx,S,E,II,JJ,X,Y,sL1,sK))*anyeq_Xi(u,X,Y) + new=u-inv(anyeq_DXi(u,rho,k,taus,A,Abar,Upsilon,Upsilon0,weights,P,N,J,approx,S,E,II,JJ,X,Y,sL1,sK))*anyeq_Xi(u,X,Y) error=norm(new-u)/norm(u) if(error<tolerance) @@ -528,8 +991,60 @@ end end +# compute u for various rho +function anyeq_hatu_rho_multithread( + u0s::Array{Array{Float64,1},1}, + P::Int64, + N::Int64, + J::Int64, + rhos::Array{Float64,1}, + a0::Float64, + weights::Tuple{Array{Float64,1},Array{Float64,1}}, + k::Array{Float64,1}, + taus::Array{Float64,1}, + V::Array{Float64,1}, + V0::Float64, + A::Array{Array{Float64,2},1}, + Abar::Array{Array{Float64,5},1}, + Upsilon::Array{Array{Float64,1},1}, + Upsilon0::Array{Float64,1}, + v::Function, + maxiter::Int64, + tolerance::Float64, + approx::Anyeq_approx, + work::Array{Array{Int64,1},1} +) + # compute u + us=Array{Array{Float64,1}}(undef,length(rhos)) + es=Array{Float64,1}(undef,length(rhos)) + errs=Array{Float64,1}(undef,length(rhos)) + count=0 + # for each worker + @sync for p in 1:length(work) + # for each task + @async for j in work[p] + count=count+1 + if count>=length(work) + progress(count,length(rhos),10000) + end + # run the task + (us[j],es[j],errs[j])=remotecall_fetch(anyeq_hatu,workers()[p],u0s[j],P,N,J,rhos[j],a0,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,v,maxiter,tolerance,approx) + end + end + + return (us,es,errs) +end + + # save Abar -function anyeq_save_Abar(taus,P,N,J,v,approx) +function anyeq_save_Abar( + taus::Array{Float64,1}, + P::Int64, + N::Int64, + J::Int64, + v::Function, + approx::Anyeq_approx +) # initialize vectors (weights,T,k,V,V0,A,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v) @@ -555,7 +1070,12 @@ function anyeq_save_Abar(taus,P,N,J,v,approx) end # read Abar -function anyeq_read_Abar(savefile,P,N,J) +function anyeq_read_Abar( + savefile::String, + P::Int64, + N::Int64, + J::Int64 +) # open file ff=open(savefile,"r") # read all lines @@ -623,13 +1143,39 @@ end # Xi # takes the vector of kj's and xn's as input -@everywhere function anyeq_Xi(U,X,Y) +@everywhere function anyeq_Xi( + U::Array{Float64,1}, + X::Array{Float64,1}, + Y::Array{Float64,1} +) return U-(Y.+1)./(2*(X.+1)).*dotPhi((Y.+1)./((X.+1).^2)) end # DXi # takes the vector of kj's as input -@everywhere function anyeq_DXi(U,rho,k,taus,v,v0,A,Abar,Upsilon,Upsilon0,weights,P,N,J,approx,S,E,II,JJ,X,Y,sL1,sK) +@everywhere function anyeq_DXi( + U::Array{Float64,1}, + rho::Float64, + k::Array{Float64,1}, + taus::Array{Float64,1}, + A::Array{Array{Float64,2},1}, + Abar::Array{Array{Float64,5},1}, + Upsilon::Array{Array{Float64,1},1}, + Upsilon0::Array{Float64,1}, + weights::Tuple{Array{Float64,1},Array{Float64,1}}, + P::Int64, + N::Int64, + J::Int64, + approx::Anyeq_approx, + S::Array{Float64,1}, + E::Float64, + II::Array{Float64,1}, + JJ::Array{Float64,1}, + X::Array{Float64,1}, + Y::Array{Float64,1}, + sL1::Array{Float64,1}, + sK::Array{Float64,1} +) out=Array{Float64,2}(undef,N*J,N*J) for zetapp in 0:J-1 for n in 1:N @@ -720,7 +1266,23 @@ end end # compute S,E,I,J,X and Y -@everywhere function anyeq_SEIJGXY(U,rho,k,taus,v,v0,A,Abar,Upsilon,Upsilon0,weights,P,N,J,approx) +@everywhere function anyeq_SEIJGXY( + U::Array{Float64,1}, + rho::Float64, + k::Array{Float64,1}, + taus::Array{Float64,1}, + v::Array{Float64,1}, + v0::Float64, + A::Array{Array{Float64,2},1}, + Abar::Array{Array{Float64,5},1}, + Upsilon::Array{Array{Float64,1},1}, + Upsilon0::Array{Float64,1}, + weights::Tuple{Array{Float64,1},Array{Float64,1}}, + P::Int64, + N::Int64, + J::Int64, + approx::Anyeq_approx +) # Chebyshev expansion of U FU=chebyshev(U,taus,weights,P,N,J,2) @@ -798,79 +1360,207 @@ end return(S,E,II,JJ,X,Y,sL1,sK,G) end + +# u(x) +@everywhere function anyeq_u_x( + x::Float64, + u::Array{Float64,1}, + k::Array{Float64,1}, + N::Int64, + J::Int64, + taus::Array{Float64,1}, + weights::Tuple{Array{Float64,1},Array{Float64,1}} +) + ux=0. + for zeta in 0:J-1 + for j in 1:N + ux+=(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 real(ux) +end + # condensate fraction -@everywhere function anyeq_eta(u,P,N,J,rho,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,approx) +@everywhere function anyeq_eta( + u::Array{Float64,1}, + P::Int64, + N::Int64, + J::Int64, + rho::Float64, + weights::Tuple{Array{Float64,1},Array{Float64,1}}, + k::Array{Float64,1}, + taus::Array{Float64,1}, + V::Array{Float64,1}, + V0::Float64, + A::Array{Array{Float64,2},1}, + Abar::Array{Array{Float64,5},1}, + Upsilon::Array{Array{Float64,1},1}, + Upsilon0::Array{Float64,1}, + approx::Anyeq_approx +) # compute dXi/dmu (S,E,II,JJ,X,Y,sL1,sK,G)=anyeq_SEIJGXY(rho*u,rho,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,weights,P,N,J,approx) dXidmu=(Y.+1)./(rho*sL1)./(2*(X.+1).^2).*dotPhi((Y.+1)./((X.+1).^2))+(Y.+1).^2 ./((X.+1).^4)./(rho*sL1).*dotdPhi((Y.+1)./(X.+1).^2) # compute eta - du=-inv(anyeq_DXi(rho*u,rho,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,weights,P,N,J,approx,S,E,II,JJ,X,Y,sL1,sK))*dXidmu + du=-inv(anyeq_DXi(rho*u,rho,k,taus,A,Abar,Upsilon,Upsilon0,weights,P,N,J,approx,S,E,II,JJ,X,Y,sL1,sK))*dXidmu eta=-avg_v_chebyshev(du,Upsilon0,k,taus,weights,N,J)/2 return eta end -# momentum distribution -@everywhere function anyeq_momentum(u,P,N,J,rho,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,approx) - # compute dXi/dlambda (without delta functions) - (S,E,II,JJ,X,Y,sL1,sK,G)=anyeq_SEIJGXY(rho*u,rho,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,weights,P,N,J,approx) - dXidlambda=-(2*pi)^3*2*u./sL1.*(dotPhi((Y.+1)./((X.+1).^2))./(2*(X.+1))+(Y.+1)./(2*(X.+1).^3).*dotdPhi((Y.+1)./(X.+1).^2)) - - # approximation for delta function (without Kronecker deltas) - delta=Array{Float64}(undef,J*N) - for zeta in 0:J-1 - for n in 1:N - delta[zeta*N+n]=2/pi^2/((taus[zeta+2]-taus[zeta+1])*weights[2][n]*cos(pi*weights[1][n]/2)*(1+k[zeta*N+n])^2*k[zeta*N+n]^2) - end - end - - # compute dXidu - dXidu=inv(anyeq_DXi(rho*u,rho,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,weights,P,N,J,approx,S,E,II,JJ,X,Y,sL1,sK)) - M=Array{Float64}(undef,J*N) - for i in 1:J*N - # du/dlambda - du=dXidu[:,i]*dXidlambda[i]*delta[i] +# correlation function +@everywhere function anyeq_2pt( + x::Float64, + u::Array{Float64,1}, + P::Int64, + N::Int64, + J::Int64, + windowL::Float64, + rho::Float64, + weights::Tuple{Array{Float64,1},Array{Float64,1}}, + k::Array{Float64,1}, + taus::Array{Float64,1}, + V::Array{Float64,1}, + V0::Float64, + A::Array{Array{Float64,2},1}, + Abar::Array{Array{Float64,5},1}, + Upsilon::Array{Array{Float64,1},1}, + Upsilon0::Array{Float64,1}, + approx::Anyeq_approx, + S::Array{Float64,1}, + E::Float64, + II::Array{Float64,1}, + JJ::Array{Float64,1}, + X::Array{Float64,1}, + Y::Array{Float64,1}, + sL1::Array{Float64,1}, + sK::Array{Float64,1} +) + g=(r,x)->sinc(r*x)*hann(r,windowL) + du=anyeq_dudv(g, x, u, P, N, J, rho, weights, k, taus, V, V0, A, Abar, Upsilon, Upsilon0, approx, S, E, II, JJ, X, Y, sL1, sK) + + C2=rho^2*(1-integrate_f_chebyshev(s->g(s,x),u,k,taus,weights,N,J)-integrate_f_chebyshev(s->1.,V.*du,k,taus,weights,N,J)) +end - # compute M - M[i]=-avg_v_chebyshev(du,Upsilon0,k,taus,weights,N,J)/2 +# uncondensed correlation function +@everywhere function anyeq_uncondensed_2pt( + xi::Float64, + uxi::Float64, + u::Array{Float64,1}, + P::Int64, + N::Int64, + J::Int64, + windowL::Float64, + rho::Float64, + weights::Tuple{Array{Float64,1},Array{Float64,1}}, + k::Array{Float64,1}, + taus::Array{Float64,1}, + V::Array{Float64,1}, + V0::Float64, + A::Array{Array{Float64,2},1}, + Abar::Array{Array{Float64,5},1}, + Upsilon::Array{Array{Float64,1},1}, + Upsilon0::Array{Float64,1}, + approx::Anyeq_approx, + S::Array{Float64,1}, + E::Float64, + II::Array{Float64,1}, + JJ::Array{Float64,1}, + X::Array{Float64,1}, + Y::Array{Float64,1}, + sL1::Array{Float64,1}, + sK::Array{Float64,1} +) + # compute dXi/dmu + g=Array{Float64,1}(undef,length(k)) + for i in 1:length(k) + g[i]=-2*rho*uxi*sinc(k[i]*xi)*hann(k[i],windowL) end + dXidmu=-g./sL1./(2*(X.+1)).*dotPhi((Y.+1)./((X.+1).^2))-(Y.+1).*g./sL1./(2(X.+1).^3).*dotdPhi((Y.+1)./(X.+1).^2) - return M -end + # compute gamma2 + du=-inv(anyeq_DXi(rho*u,rho,k,taus,A,Abar,Upsilon,Upsilon0,weights,P,N,J,approx,S,E,II,JJ,X,Y,sL1,sK))*dXidmu + gamma2=-avg_v_chebyshev(du,Upsilon0,k,taus,weights,N,J)/2 + return gamma2 +end -# correlation function -@everywhere function anyeq_2pt(x,u,P,N,J,windowL,rho,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,approx,S,E,II,JJ,X,Y,sL1,sK,G) +# compute the directional derivative of u with respect to v in direction g +@everywhere function anyeq_dudv( + g::Function,# should be of the form g(k,x) where x is a parameter + x::Float64, + u::Array{Float64,1}, + P::Int64, + N::Int64, + J::Int64, + rho::Float64, + weights::Tuple{Array{Float64,1},Array{Float64,1}}, + k::Array{Float64,1}, + taus::Array{Float64,1}, + V::Array{Float64,1}, + V0::Float64, + A::Array{Array{Float64,2},1}, + Abar::Array{Array{Float64,5},1}, + Upsilon::Array{Array{Float64,1},1}, + Upsilon0::Array{Float64,1}, + approx::Anyeq_approx, + S::Array{Float64,1}, + E::Float64, + II::Array{Float64,1}, + JJ::Array{Float64,1}, + X::Array{Float64,1}, + Y::Array{Float64,1}, + sL1::Array{Float64,1}, + sK::Array{Float64,1} +) # initialize dV - dV=Array{Float64}(undef,J*N) + dV=Array{Float64,1}(undef,J*N) for i in 1:J*N - if x>0 - dV[i]=sin(k[i]*x)/(k[i]*x)*hann(k[i],windowL) - else - dV[i]=hann(k[i],windowL) - end + dV[i]=g(k[i],x) end - dV0=1. + dV0=g(0.,x) # compute dUpsilon # Upsilonmat does not use splines, so increase precision weights_plus=gausslegendre(N*J) - dUpsilon=Upsilonmat(k,r->sin(r*x)/(r*x)*hann(r,windowL),weights_plus) - dUpsilon0=Upsilon0mat(k,r->sin(r*x)/(r*x)*hann(r,windowL),weights_plus) + dUpsilon=Upsilonmat(k,s->g(s,x),weights_plus) + dUpsilon0=Upsilon0mat(k,s->g(s,x),weights_plus) - du=-inv(anyeq_DXi(rho*u,rho,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,weights,P,N,J,approx,S,E,II,JJ,X,Y,sL1,sK))*anyeq_dXidv(x,rho*u,rho,k,taus,dV,dV0,A,Abar,dUpsilon,dUpsilon0,weights,P,N,J,approx,S,E,II,JJ,X,Y,sL1,sK) + du=-inv(anyeq_DXi(rho*u,rho,k,taus,A,Abar,Upsilon,Upsilon0,weights,P,N,J,approx,S,E,II,JJ,X,Y,sL1,sK))*anyeq_dXidv(rho*u,rho,k,taus,dV,dV0,A,Abar,dUpsilon,dUpsilon0,weights,P,N,J,approx,S,E,II,JJ,X,Y,sL1,sK) # rescale rho du=du/rho - C2=rho^2*(1-integrate_f_chebyshev(s->1.,u.*dV+V.*du,k,taus,weights,N,J)) - - return C2 + return du end -# derivative of Xi with respect to v in the direction sin(kx)/kx -@everywhere function anyeq_dXidv(x,U,rho,k,taus,dv,dv0,A,Abar,dUpsilon,dUpsilon0,weights,P,N,J,approx,S,E,II,JJ,X,Y,sL1,sK) +# derivative of Xi with respect to v in the specified by dUpsilon and dUpsilon0 +@everywhere function anyeq_dXidv( + U::Array{Float64,1}, + rho::Float64, + k::Array{Float64,1}, + taus::Array{Float64,1}, + dv::Array{Float64,1}, + dv0::Float64, + A::Array{Array{Float64,2},1}, + Abar::Array{Array{Float64,5},1}, + dUpsilon::Array{Array{Float64,1},1}, + dUpsilon0::Array{Float64,1}, + weights::Tuple{Array{Float64,1},Array{Float64,1}}, + P::Int64, + N::Int64, + J::Int64, + approx::Anyeq_approx, + S::Array{Float64,1}, + E::Float64, + II::Array{Float64,1}, + JJ::Array{Float64,1}, + X::Array{Float64,1}, + Y::Array{Float64,1}, + sL1::Array{Float64,1}, + sK::Array{Float64,1} +) # Chebyshev expansion of U FU=chebyshev(U,taus,weights,P,N,J,2) @@ -896,7 +1586,7 @@ end dJJ+=approx.gL3*approx.bL3*double_conv_S_chebyshev(FU,FU,FU,FU,dFS,Abar) end if approx.bL3!=1. - dJJ=approx.gL3*(1-approx.bL3)*dE*(UU/rho).^2 + dJJ+=approx.gL3*(1-approx.bL3)*dE*(UU/rho).^2 end end @@ -947,3 +1637,220 @@ end out=((Y.+1).*dX./(X.+1)-dY)./(2*(X.+1)).*dotPhi((Y.+1)./((X.+1).^2))+(Y.+1)./(2*(X.+1).^3).*(2*(Y.+1)./(X.+1).*dX-dY).*dotdPhi((Y.+1)./(X.+1).^2) return out end + +# maximum of 2 point correlation function +@everywhere function anyeq_2pt_max( + u::Array{Float64,1}, + P::Int64, + N::Int64, + J::Int64, + x0::Float64, + dx::Float64, + maxstep::Float64, + maxiter::Int64, + tolerance::Float64, + windowL::Float64, + rho::Float64, + weights::Tuple{Array{Float64,1},Array{Float64,1}}, + T::Array{Polynomial,1}, + k::Array{Float64,1}, + taus::Array{Float64,1}, + V::Array{Float64,1}, + V0::Float64, + A::Array{Array{Float64,2},1}, + Abar::Array{Array{Float64,5},1}, + Upsilon::Array{Array{Float64,1},1}, + Upsilon0::Array{Float64,1}, + approx::Anyeq_approx +) + # compute some useful integrals + (S,E,II,JJ,X,Y,sL1,sK,G)=anyeq_SEIJGXY(rho*u,rho,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,weights,P,N,J,approx) + + (x,f)=newton_maximum(y->anyeq_2pt(y,u,P,N,J,windowL,rho,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,approx,S,E,II,JJ,X,Y,sL1,sK),x0,dx,maxiter,tolerance,maxstep) + + return(x,f) +end + + +# Fourier transform of 2pt correlation function at q +@everywhere function anyeq_2pt_fourier( + q::Float64, + u::Array{Float64,1}, + P::Int64, + N::Int64, + J::Int64, + windowL::Float64, + rho::Float64, + weights::Tuple{Array{Float64,1},Array{Float64,1}}, + k::Array{Float64,1}, + taus::Array{Float64,1}, + V::Array{Float64,1}, + V0::Float64, + A::Array{Array{Float64,2},1}, + Abar::Array{Array{Float64,5},1}, + Upsilon::Array{Array{Float64,1},1}, + Upsilon0::Array{Float64,1}, + approx::Anyeq_approx, + S::Array{Float64,1}, + E::Float64, + II::Array{Float64,1}, + JJ::Array{Float64,1}, + X::Array{Float64,1}, + Y::Array{Float64,1}, + sL1::Array{Float64,1}, + sK::Array{Float64,1} +) + # direction in which to differentiate u + weights_plus=gausslegendre(N*J) + #g=(r,x)->(r>0. ? (1.)/(2*x*r)*integrate_legendre(s->s*gaussian(s,(1.)/windowL),abs(x-r),x+r,weights_plus) : gaussian(x,(1.)/windowL)) + g=(r,x)->(r>0. ? (1.)/(2*x*r*windowL)*(gaussian(x-r,(1.)/windowL)-gaussian(x+r,(1.)/windowL)) : gaussian(x,(1.)/windowL)) + + du=anyeq_dudv(g,q,u,P,N,J,rho,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,approx,S,E,II,JJ,X,Y,sL1,sK) + + C2=rho^2*(-integrate_f_chebyshev(s->g(s,q),u,k,taus,weights,N,J)-integrate_f_chebyshev(s->1.,V.*du,k,taus,weights,N,J)) + + return C2 +end + +# maximum of Fourier transform of 2 point correlation function +@everywhere function anyeq_2pt_fourier_max( + u::Array{Float64,1}, + P::Int64, + N::Int64, + J::Int64, + k0::Float64, + dk::Float64, + maxstep::Float64, + maxiter::Int64, + tolerance::Float64, + windowL::Float64, + rho::Float64, + weights::Tuple{Array{Float64,1},Array{Float64,1}}, + T::Array{Polynomial,1}, + k::Array{Float64,1}, + taus::Array{Float64,1}, + V::Array{Float64,1}, + V0::Float64, + A::Array{Array{Float64,2},1}, + Abar::Array{Array{Float64,5},1}, + Upsilon::Array{Array{Float64,1},1}, + Upsilon0::Array{Float64,1}, + approx::Anyeq_approx +) + # compute some useful integrals + (S,E,II,JJ,X,Y,sL1,sK,G)=anyeq_SEIJGXY(rho*u,rho,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,weights,P,N,J,approx) + + (ko,f)=newton_maximum(y->anyeq_2pt_fourier(y,u,P,N,J,windowL,rho,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,approx,S,E,II,JJ,X,Y,sL1,sK),k0,dk,maxiter,tolerance,maxstep) + + return(ko,f) +end + +# momentum distribution, computed using a Gaussian window +@everywhere function anyeq_momentum_window( + kmin::Float64, + kmax::Float64, + u::Array{Float64,1}, + P::Int64, + N::Int64, + J::Int64, + windowL::Float64, # L is windowL/k^2 + rho::Float64, + weights::Tuple{Array{Float64,1},Array{Float64,1}}, + k::Array{Float64,1}, + taus::Array{Float64,1}, + V::Array{Float64,1}, + V0::Float64, + A::Array{Array{Float64,2},1}, + Abar::Array{Array{Float64,5},1}, + Upsilon::Array{Array{Float64,1},1}, + Upsilon0::Array{Float64,1}, + approx::Anyeq_approx +) + + # compute dXi/dlambda without the delta function of u(q) + (S,E,II,JJ,X,Y,sL1,sK,G)=anyeq_SEIJGXY(rho*u,rho,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,weights,P,N,J,approx) + dXidlambda=(16*pi^3)*(dotPhi((Y.+1)./((X.+1).^2))./(2*(X.+1))+(Y.+1)./(2*(X.+1).^3).*dotdPhi((Y.+1)./(X.+1).^2))./sL1 + + # compute dXidu + dXidu=inv(anyeq_DXi(rho*u,rho,k,taus,A,Abar,Upsilon,Upsilon0,weights,P,N,J,approx,S,E,II,JJ,X,Y,sL1,sK)) + + M=Array{Float64,1}(undef,N*J) + for j in 1:J*N + # the Chebyshev polynomial expansion is often not good enough to compute M away from k[i] + q=k[j] + + # drop the computation if not in k interval + if q<kmin || q>kmax + continue + end + + # delta function + delta=Array{Float64,1}(undef,J*N) + L=windowL/q^2 + for i in 1:J*N + delta[i]=(1.)/(2*k[i]*q*L)*(gaussian(k[i]-q,(1.)/L)-gaussian(k[i]+q,(1.)/L)) + end + + # du/dlambda + du=-dXidu*(dXidlambda.*delta*u[j]) + # rescale u + du=du/rho + + # compute M + M[j]=-avg_v_chebyshev(du,Upsilon0,k,taus,weights,N,J)*rho/2 + end + + return M +end + +# momentum distribution, computed using a discrete approximation of the delta function +@everywhere function anyeq_momentum_discrete_delta( + kmin::Float64, + kmax::Float64, + u::Array{Float64,1}, + P::Int64, + N::Int64, + J::Int64, + rho::Float64, + weights::Tuple{Array{Float64,1},Array{Float64,1}}, + k::Array{Float64,1}, + taus::Array{Float64,1}, + V::Array{Float64,1}, + V0::Float64, + A::Array{Array{Float64,2},1}, + Abar::Array{Array{Float64,5},1}, + Upsilon::Array{Array{Float64,1},1}, + Upsilon0::Array{Float64,1}, + approx::Anyeq_approx +) + # compute dXi/dlambda (without delta functions) + (S,E,II,JJ,X,Y,sL1,sK,G)=anyeq_SEIJGXY(rho*u,rho,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,weights,P,N,J,approx) + dXidlambda=(2*pi)^3*2*u./sL1.*(dotPhi((Y.+1)./((X.+1).^2))./(2*(X.+1))+(Y.+1)./(2*(X.+1).^3).*dotdPhi((Y.+1)./(X.+1).^2)) + + # approximation for delta function (without Kronecker deltas) + delta=Array{Float64,1}(undef,J*N) + for zeta in 0:J-1 + for n in 1:N + delta[zeta*N+n]=8/pi/((taus[zeta+2]-taus[zeta+1])*weights[2][n]*cos(pi*weights[1][n]/2)*(1+k[zeta*N+n])^2) + end + end + + # compute dXidu + dXidu=inv(anyeq_DXi(rho*u,rho,k,taus,A,Abar,Upsilon,Upsilon0,weights,P,N,J,approx,S,E,II,JJ,X,Y,sL1,sK)) + + M=Array{Float64,1}(undef,J*N) + for i in 1:J*N + # drop the computation if not in k interval + if k[i]<kmin || k[i]>kmax + continue + end + + # du/dlambda + du=-dXidu[:,i]*dXidlambda[i]*delta[i] + + # compute M + M[i]=-avg_v_chebyshev(du,Upsilon0,k,taus,weights,N,J)/2 + end + + return M +end diff --git a/src/chebyshev.jl b/src/chebyshev.jl index 28c8f1f..af4be40 100644 --- a/src/chebyshev.jl +++ b/src/chebyshev.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,7 +13,15 @@ ## limitations under the License. # Chebyshev expansion -@everywhere function chebyshev(a,taus,weights,P,N,J,nu) +@everywhere function chebyshev( + a::Array{Float64,1}, + taus::Array{Float64,1}, + weights::Tuple{Array{Float64,1},Array{Float64,1}}, + P::Int64, + N::Int64, + J::Int64, + nu::Int64 +) out=zeros(Float64,J*(P+1)) for zeta in 0:J-1 for n in 0:P @@ -27,7 +35,15 @@ end # evaluate function from Chebyshev expansion -@everywhere function chebyshev_eval(Fa,x,taus,chebyshev,P,J,nu) +@everywhere function chebyshev_eval( + Fa::Array{Float64,1}, + x::Float64, + taus::Array{Float64,1}, + chebyshev::Array{Polynomial,1}, + P::Int64, + J::Int64, + nu::Int64 +) # change variable tau=(1-x)/(1+x) @@ -48,7 +64,11 @@ end # convolution # input the Chebyshev expansion of a and b, as well as the A matrix -@everywhere function conv_chebyshev(Fa,Fb,A) +@everywhere function conv_chebyshev( + Fa::Array{Float64,1}, + Fb::Array{Float64,1}, + A::Array{Array{Float64,2},1} +) out=zeros(Float64,length(A)) for i in 1:length(A) out[i]=dot(Fa,A[i]*Fb) @@ -57,12 +77,27 @@ end end # <ab> -@everywhere function avg_chebyshev(Fa,Fb,A0) +@everywhere function avg_chebyshev( + Fa::Array{Float64,1}, + Fb::Array{Float64,1}, + A0::Float64 +) return dot(Fa,A0*Fb) end # 1_n * a -@everywhere function conv_one_chebyshev(n,zetapp,Fa,A,taus,weights,P,N,J,nu1) +@everywhere function conv_one_chebyshev( + n::Int64, + zetapp::Int64, + Fa::Array{Float64,1}, + A::Array{Array{Float64,2},1}, + taus::Array{Float64,1}, + weights::Tuple{Array{Float64,1},Array{Float64,1}}, + P::Int64, + N::Int64, + J::Int64, + nu1::Int64 +) out=zeros(Float64,N*J) for m in 1:N*J for l in 0:P @@ -74,7 +109,15 @@ end return out end # a * v -@everywhere function conv_v_chebyshev(a,Upsilon,k,taus,weights,N,J) +@everywhere function conv_v_chebyshev( + a::Array{Float64,1}, + Upsilon::Array{Array{Float64,1},1}, + k::Array{Float64,1}, + taus::Array{Float64,1}, + weights::Tuple{Array{Float64,1},Array{Float64,1}}, + N::Int64, + J::Int64 +) out=zeros(Float64,J*N) for i in 1:J*N for zetap in 0:J-1 @@ -86,7 +129,16 @@ end end return out end -@everywhere function conv_one_v_chebyshev(n,zetap,Upsilon,k,taus,weights,N,J) +@everywhere function conv_one_v_chebyshev( + n::Int64, + zetap::Int64, + Upsilon::Array{Array{Float64,1},1}, + k::Array{Float64,1}, + taus::Array{Float64,1}, + weights::Tuple{Array{Float64,1},Array{Float64,1}}, + N::Int64, + J::Int64 +) out=zeros(Float64,J*N) xj=weights[1][n] for i in 1:J*N @@ -96,7 +148,14 @@ end end # <av> -@everywhere function avg_v_chebyshev(a,Upsilon0,k,taus,weights,N,J) +@everywhere function avg_v_chebyshev(a, + Upsilon0::Array{Float64,1}, + k::Array{Float64,1}, + taus::Array{Float64,1}, + weights::Tuple{Array{Float64,1},Array{Float64,1}}, + N::Int64, + J::Int64 +) out=0. for zetap in 0:J-1 for j in 1:N @@ -107,13 +166,28 @@ end return out end # <1_nv> -@everywhere function avg_one_v_chebyshev(n,zetap,Upsilon0,k,taus,weights,N) +@everywhere function avg_one_v_chebyshev( + n::Int64, + zetap::Int64, + Upsilon0::Array{Float64,1}, + k::Array{Float64,1}, + taus::Array{Float64,1}, + weights::Tuple{Array{Float64,1},Array{Float64,1}}, + N::Int64 +) 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) +@everywhere function double_conv_S_chebyshev( + FU1::Array{Float64,1}, + FU2::Array{Float64,1}, + FU3::Array{Float64,1}, + FU4::Array{Float64,1}, + FU5::Array{Float64,1}, + Abar::Array{Float64,5} +) out=zeros(Float64,length(Abar)) for i in 1:length(Abar) for j1 in 1:length(FU1) @@ -133,7 +207,17 @@ end # compute A -@everywhere function Amat(k,weights,taus,T,P,N,J,nua,nub) +@everywhere function Amat( + k::Array{Float64,1}, + weights::Tuple{Array{Float64,1},Array{Float64,1}}, + taus::Array{Float64,1}, + T::Array{Polynomial,1}, + P::Int64, + N::Int64, + J::Int64, + nua::Int64, + nub::Int64 +) 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)) @@ -152,7 +236,11 @@ end end # compute Upsilon -@everywhere function Upsilonmat(k,v,weights) +@everywhere function Upsilonmat( + k::Array{Float64,1}, + v::Function, + weights::Tuple{Array{Float64,1},Array{Float64,1}} +) out=Array{Array{Float64,1},1}(undef,length(k)) for i in 1:length(k) out[i]=Array{Float64,1}(undef,length(k)) @@ -162,7 +250,11 @@ end end return out end -@everywhere function Upsilon0mat(k,v,weights) +@everywhere function Upsilon0mat( + k::Array{Float64,1}, + v::Function, + weights::Tuple{Array{Float64,1},Array{Float64,1}} +) out=Array{Float64,1}(undef,length(k)) for j in 1:length(k) out[j]=2*k[j]*v(k[j]) @@ -171,17 +263,36 @@ end end # alpha_- -@everywhere function alpham(k,t) +@everywhere function alpham( + k::Float64, + t::Float64 +) return (1-k-(1-t)/(1+t))/(1+k+(1-t)/(1+t)) end # alpha_+ -@everywhere function alphap(k,t) +@everywhere function alphap( + k::Float64, + t::Float64 +) 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) +@everywhere function barAmat( + k::Float64, + weights::Tuple{Array{Float64,1},Array{Float64,1}}, + taus::Array{Float64,1}, + T::Array{Polynomial,1}, + P::Int64, + N::Int64, + J::Int64, + nu1::Int64, + nu2::Int64, + nu3::Int64, + nu4::Int64, + nu5::Int64 +) 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 @@ -211,27 +322,107 @@ 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) +@everywhere function barAmat_int1(tau, + k::Float64, + taus::Array{Float64,1}, + T::Array{Polynomial,1}, + weights::Tuple{Array{Float64,1},Array{Float64,1}}, + nu1::Int64, + nu2::Int64, + nu3::Int64, + nu4::Int64, + nu5::Int64, + zeta1::Int64, + zeta2::Int64, + zeta3::Int64, + zeta4::Int64, + zeta5::Int64, + n1::Int64, + n2::Int64, + n3::Int64, + n4::Int64, + n5::Int64 +) 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) +@everywhere function barAmat_int2(tau, + sigma::Float64, + k::Float64, + taus::Array{Float64,1}, + T::Array{Polynomial,1}, + weights::Tuple{Array{Float64,1},Array{Float64,1}}, + nu2::Int64, + nu3::Int64, + nu4::Int64, + nu5::Int64, + zeta2::Int64, + zeta3::Int64, + zeta4::Int64, + zeta5::Int64, + n2::Int64, + n3::Int64, + n4::Int64, + n5::Int64 +) 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) +@everywhere function barAmat_int3(tau, + sigma::Float64, + taup::Float64, + k::Float64, + taus::Array{Float64,1}, + T::Array{Polynomial,1}, + weights::Tuple{Array{Float64,1},Array{Float64,1}}, + nu3::Int64, + nu4::Int64, + nu5::Int64, + zeta3::Int64, + zeta4::Int64, + zeta5::Int64, + n3::Int64, + n4::Int64, + n5::Int64 +) 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) +@everywhere function barAmat_int4(tau, + sigma::Float64, + taup::Float64, + sigmap::Float64, + k::Float64, + taus::Array{Float64,1}, + T::Array{Polynomial,1}, + weights::Tuple{Array{Float64,1},Array{Float64,1}}, + nu4::Int64, + nu5::Int64, + zeta4::Int64, + zeta5::Int64, + n4::Int64, + n5::Int64 +) 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) +@everywhere function barAmat_int5(tau, + sigma::Float64, + taup::Float64, + sigmap::Float64, + theta::Float64, + k::Float64, + taus::Array{Float64,1}, + T::Array{Polynomial,1}, + weights::Tuple{Array{Float64,1},Array{Float64,1}}, + nu5::Int64, + zeta5::Int64, + n5::Int64 +) 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])) @@ -240,13 +431,22 @@ end 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) +@everywhere function barAmat_R( + s::Float64, + t::Float64, + sp::Float64, + tp::Float64, + theta::Float64, + k::Float64 +) + 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) +@everywhere function chebyshev_polynomials( + P::Int64 +) + T=Array{Polynomial,1}(undef,P+1) T[1]=Polynomial([1]) T[2]=Polynomial([0,1]) for n in 1:P-1 @@ -258,7 +458,15 @@ end end # compute \int f*u dk/(2*pi)^3 -@everywhere function integrate_f_chebyshev(f,u,k,taus,weights,N,J) +@everywhere function integrate_f_chebyshev( + f::Function, + u::Array{Float64,1}, + k::Array{Float64,1}, + taus::Array{Float64,1}, + weights::Tuple{Array{Float64,1},Array{Float64,1}}, + N::Int64, + J::Int64 +) out=0. for zeta in 0:J-1 for i in 1:N @@ -268,7 +476,15 @@ end return out end -@everywhere function inverse_fourier_chebyshev(u,x,k,taus,weights,N,J) +@everywhere function inverse_fourier_chebyshev( + u::Array{Float64,1}, + x::Float64, + k::Array{Float64,1}, + taus::Array{Float64,1}, + weights::Tuple{Array{Float64,1},Array{Float64,1}}, + N::Int64, + J::Int64 +) out=0. for zeta in 0:J-1 for j in 1:N @@ -277,3 +493,54 @@ end end return out end + +# compute B (for the computation of the fourier transform of the two-point correlation) +@everywhere function Bmat( + q::Float64, + k::Array{Float64,1}, + weights::Tuple{Array{Float64,1},Array{Float64,1}}, + taus::Array{Float64,1}, + T::Array{Polynomial,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]=zeros(Float64,J*(P+1)) + for zeta in 0:J-1 + for n in 0:P + out[i][zeta*(P+1)+n+1]=1/(8*pi^3*k[i]*q)*(betam(k[i],q)>taus[zeta+2] || betap(k[i],q)<taus[zeta+1] ? 0. : integrate_legendre(sigma->(1-sigma)/(1+sigma)^(3-nu)*T[n+1]((2*sigma-(taus[zeta+1]+taus[zeta+2]))/(taus[zeta+2]-taus[zeta+1])),max(taus[zeta+1],betam(k[i],q)),min(taus[zeta+2],betap(k[i],q)),weights)) + end + end + end + + return out +end +# beta_- +@everywhere function betam( + k::Float64, + q::Float64 +) + return (1-k-q)/(1+k+q) +end +# beta_+ +@everywhere function betap( + k::Float64, + q::Float64 +) + return (1-abs(k-q))/(1+abs(k-q)) +end + +# mathfrak S (for the computation of the fourier transform of the two-point correlation) +@everywhere function chebyshev_frakS( + Ff::Array{Float64,1}, + B::Array{Array{Float64,1},1} +) + out=zeros(Float64,length(B)) + for i in 1:length(B) + out[i]=dot(Ff,B[i]) + end + return out +end diff --git a/src/easyeq.jl b/src/easyeq.jl index 0bde3ab..2dbbd1a 100644 --- a/src/easyeq.jl +++ b/src/easyeq.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. @@ -19,47 +19,65 @@ end # compute energy -function easyeq_energy(minlrho_init,nlrho_init,order,rho,a0,v,maxiter,tolerance,approx) +function easyeq_energy( + minlrho_init::Float64, + nlrho_init::Int64, + order::Int64, + rho::Float64, + a0::Float64, + v::Function, + maxiter::Int64, + tolerance::Float64, + approx::Easyeq_approx +) # compute gaussian quadrature weights weights=gausslegendre(order) - # compute initial guess from previous rho - (u,E,err)=easyeq_hatu(easyeq_init_u(a0,order,weights),order,(10.)^minlrho_init,v,maxiter,tolerance,weights,approx) - for j in 2:nlrho_init - rho_tmp=10^(minlrho_init+(log10(rho)-minlrho_init)*(j-1)/(nlrho_init-1)) - (u,E,err)=easyeq_hatu(u,order,rho_tmp,v,maxiter,tolerance,weights,approx) - end + # compute u + (u,E,err)= easyeq_compute_u_prevrho([rho],minlrho_init,nlrho_init,a0,order,v,maxiter,tolerance,weights,approx) # print energy @printf("% .15e % .15e\n",real(E),err) end # compute energy as a function of rho -function easyeq_energy_rho(rhos,order,a0,v,maxiter,tolerance,approx) +function easyeq_energy_rho( + rhos::Array{Float64,1}, + minlrho_init::Float64, + nlrho_init::Int64, + order::Int64, + a0::Float64, + v::Function, + maxiter::Int64, + tolerance::Float64, + approx::Easyeq_approx +) # compute gaussian quadrature weights weights=gausslegendre(order) - # init u - u=easyeq_init_u(a0,order,weights) + # compute u + (us,es,errs)= easyeq_compute_u_prevrho(rhos,minlrho_init,nlrho_init,a0,order,v,maxiter,tolerance,weights,approx) for j in 1:length(rhos) - # compute u (init newton with previously computed u) - (u,E,err)=easyeq_hatu(u,order,rhos[j],v,maxiter,tolerance,weights,approx) - - @printf("% .15e % .15e % .15e\n",rhos[j],real(E),err) - + @printf("% .15e % .15e % .15e\n",rhos[j],real(es[j]),errs[j]) end end # compute u(k) -function easyeq_uk(minlrho_init,nlrho_init,order,rho,a0,v,maxiter,tolerance,approx) +function easyeq_uk( + minlrho_init::Float64, + nlrho_init::Int64, + order::Int64, + rho::Float64, + a0::Float64, + v::Function, + maxiter::Int64, + tolerance::Float64, + approx::Easyeq_approx +) weights=gausslegendre(order) - # compute initial guess from previous rho - (u,E,err)=easyeq_hatu(easyeq_init_u(a0,order,weights),order,(10.)^minlrho_init,v,maxiter,tolerance,weights,approx) - for j in 2:nlrho_init - rho_tmp=10^(minlrho_init+(log10(rho)-minlrho_init)*(j-1)/(nlrho_init-1)) - (u,E,err)=easyeq_hatu(u,order,rho_tmp,v,maxiter,tolerance,weights,approx) - end + # compute u + (u,E,err)= easyeq_compute_u_prevrho_error([rho],minlrho_init,nlrho_init,a0,order,v,maxiter,tolerance,weights,approx) for i in 1:order k=(1-weights[1][i])/(1+weights[1][i]) @@ -68,15 +86,24 @@ function easyeq_uk(minlrho_init,nlrho_init,order,rho,a0,v,maxiter,tolerance,appr end # compute u(x) -function easyeq_ux(minlrho_init,nlrho_init,order,rho,a0,v,maxiter,tolerance,xmin,xmax,nx,approx) +function easyeq_ux( + minlrho_init::Float64, + nlrho_init::Int64, + order::Int64, + rho::Float64, + a0::Float64, + v::Function, + maxiter::Int64, + tolerance::Float64, + xmin::Float64, + xmax::Float64, + nx::Int64, + approx::Easyeq_approx +) weights=gausslegendre(order) - # compute initial guess from previous rho - (u,E,err)=easyeq_hatu(easyeq_init_u(a0,order,weights),order,(10.)^minlrho_init,v,maxiter,tolerance,weights,approx) - for j in 2:nlrho_init - rho_tmp=10^(minlrho_init+(log10(rho)-minlrho_init)*(j-1)/(nlrho_init-1)) - (u,E,err)=easyeq_hatu(u,order,rho_tmp,v,maxiter,tolerance,weights,approx) - end + # compute u + (u,E,err)= easyeq_compute_u_prevrho_error([rho],minlrho_init,nlrho_init,a0,order,v,maxiter,tolerance,weights,approx) for i in 1:nx x=xmin+(xmax-xmin)*i/nx @@ -85,15 +112,24 @@ function easyeq_ux(minlrho_init,nlrho_init,order,rho,a0,v,maxiter,tolerance,xmin end # compute 2u(x)-rho u*u(x) -function easyeq_uux(minlrho_init,nlrho_init,order,rho,a0,v,maxiter,tolerance,xmin,xmax,nx,approx) +function easyeq_uux( + minlrho_init::Float64, + nlrho_init::Int64, + order::Int64, + rho::Float64, + a0::Float64, + v::Function, + maxiter::Int64, + tolerance::Float64, + xmin::Float64, + xmax::Float64, + nx::Int64, + approx::Easyeq_approx +) weights=gausslegendre(order) - # compute initial guess from previous rho - (u,E,err)=easyeq_hatu(easyeq_init_u(a0,order,weights),order,(10.)^minlrho_init,v,maxiter,tolerance,weights,approx) - for j in 2:nlrho_init - rho_tmp=10^(minlrho_init+(log10(rho)-minlrho_init)*(j-1)/(nlrho_init-1)) - (u,E,err)=easyeq_hatu(u,order,rho_tmp,v,maxiter,tolerance,weights,approx) - end + # compute u + (u,E,err)= easyeq_compute_u_prevrho_error([rho],minlrho_init,nlrho_init,a0,order,v,maxiter,tolerance,weights,approx) for i in 1:nx x=xmin+(xmax-xmin)*i/nx @@ -102,16 +138,22 @@ function easyeq_uux(minlrho_init,nlrho_init,order,rho,a0,v,maxiter,tolerance,xmi end # condensate fraction -function easyeq_condensate_fraction(minlrho_init,nlrho_init,order,rho,a0,v,maxiter,tolerance,approx) +function easyeq_condensate_fraction( + minlrho_init::Float64, + nlrho_init::Int64, + order::Int64, + rho::Float64, + a0::Float64, + v::Function, + maxiter::Int64, + tolerance::Float64, + approx::Easyeq_approx +) # compute gaussian quadrature weights weights=gausslegendre(order) - # compute initial guess from previous rho - (u,E,err)=easyeq_hatu(easyeq_init_u(a0,order,weights),order,(10.)^minlrho_init,v,maxiter,tolerance,weights,approx) - for j in 2:nlrho_init - rho_tmp=10^(minlrho_init+(log10(rho)-minlrho_init)*(j-1)/(nlrho_init-1)) - (u,E,err)=easyeq_hatu(u,order,rho_tmp,v,maxiter,tolerance,weights,approx) - end + # compute u + (u,E,err)= easyeq_compute_u_prevrho([rho],minlrho_init,nlrho_init,a0,order,v,maxiter,tolerance,weights,approx) # compute eta eta=easyeq_eta(u,order,rho,v,maxiter,tolerance,weights,approx) @@ -121,25 +163,350 @@ function easyeq_condensate_fraction(minlrho_init,nlrho_init,order,rho,a0,v,maxit end # condensate fraction as a function of rho -function easyeq_condensate_fraction_rho(rhos,order,a0,v,maxiter,tolerance,approx) +function easyeq_condensate_fraction_rho( + rhos::Array{Float64,1}, + minlrho_init::Float64, + nlrho_init::Int64, + order::Int64, + a0::Float64, + v::Function, + maxiter::Int64, + tolerance::Float64, + approx::Easyeq_approx +) weights=gausslegendre(order) - # init u - u=easyeq_init_u(a0,order,weights) + # compute u + (us,es,errs)= easyeq_compute_u_prevrho(rhos,minlrho_init,nlrho_init,a0,order,v,maxiter,tolerance,weights,approx) for j in 1:length(rhos) - # compute u (init newton with previously computed u) - (u,E,err)=easyeq_hatu(u,order,rhos[j],v,maxiter,tolerance,weights,approx) - # compute eta - eta=easyeq_eta(u,order,rhos[j],v,maxiter,tolerance,weights,approx) + eta=easyeq_eta(us[j],order,rhos[j],v,maxiter,tolerance,weights,approx) + @printf("% .15e % .15e % .15e\n",rhos[j],eta,errs[j]) + end +end + +# 2 pt correlation function +function easyeq_2pt( + xmin::Float64, + xmax::Float64, + nx::Int64, + minlrho_init::Float64, + nlrho_init::Int64, + order::Int64, + windowL::Float64, + rho::Float64, + a0::Float64, + v::Function, + maxiter::Int64, + tolerance::Float64, + approx::Easyeq_approx +) + # compute gaussian quadrature weights + weights=gausslegendre(order) + + # compute u + (u,E,err)= easyeq_compute_u_prevrho_error([rho],minlrho_init,nlrho_init,a0,order,v,maxiter,tolerance,weights,approx) - @printf("% .15e % .15e % .15e\n",rhos[j],eta,err) + # compute useful terms + (V,V0)=easyeq_init_v(weights,v) + (Eta,Eta0)=easyeq_init_H(weights,v) + (E,S,A,T,B,X)=easyeq_ESATBX(rho*u,V,V0,Eta,Eta0,weights,rho,approx) + + # compute C2 + for j in 1:nx + x=xmin+(xmax-xmin)/nx*j + C2=easyeq_C2(x,u,windowL,rho,weights,Eta,Eta0,approx,E,S,A,T,B,X) + @printf("% .15e % .15e\n",x,C2) end end +# maximum of 2 point correlation function +function easyeq_2pt_max( + dx::Float64, + x0::Float64, # initial guess is x0/rho^(1/3) + minlrho_init::Float64, + nlrho_init::Int64, + order::Int64, + windowL::Float64, + rho::Float64, + a0::Float64, + v::Function, + maxstep::Float64, + maxiter::Int64, + tolerance::Float64, + tolerance_max::Float64, + approx::Easyeq_approx +) + # compute gaussian quadrature weights + weights=gausslegendre(order) + + # compute u + (u,E,err)= easyeq_compute_u_prevrho_error([rho],minlrho_init,nlrho_init,a0,order,v,maxiter,tolerance,weights,approx) + + # compute useful terms + (V,V0)=easyeq_init_v(weights,v) + (Eta,Eta0)=easyeq_init_H(weights,v) -# initialize u -@everywhere function easyeq_init_u(a0,order,weights) + (x,f)=easyeq_C2_max(u,x0/rho^(1/3),dx,maxstep,maxiter,tolerance_max,windowL,rho,weights,V,V0,Eta,Eta0,approx) + + if(x==Inf) + @printf(stderr,"max search failed for rho=%e\n",rho) + else + @printf("% .15e % .15e\n",x,f) + end +end + +# maximum of 2 point correlation function as a function of rho +function easyeq_2pt_max_rho( + rhos::Array{Float64,1}, + dx::Float64, + x0::Float64, # initial guess is x0/rho^(1/3) + minlrho_init::Float64, + nlrho_init::Int64, + order::Int64, + windowL::Float64, + a0::Float64, + v::Function, + maxstep::Float64, + maxiter::Int64, + tolerance::Float64, + tolerance_max::Float64, + approx::Easyeq_approx +) + # compute gaussian quadrature weights + weights=gausslegendre(order) + + # compute u + (us,es,errs)= easyeq_compute_u_prevrho_error(rhos,minlrho_init,nlrho_init,a0,order,v,maxiter,tolerance,weights,approx) + + # compute useful terms + (V,V0)=easyeq_init_v(weights,v) + (Eta,Eta0)=easyeq_init_H(weights,v) + + # save result from each task + xs=Array{Float64,1}(undef,length(rhos)) + fs=Array{Float64,1}(undef,length(rhos)) + + # spawn workers + work=spawn_workers(length(rhos)) + + count=0 + # for each worker + @sync for p in 1:length(work) + # for each task + @async for j in work[p] + count=count+1 + if count>=length(work) + progress(count,length(rhos),10000) + end + # run the task + (xs[j],fs[j])=remotecall_fetch(easyeq_C2_max,workers()[p],us[j],x0/rhos[j]^(1/3),dx,maxstep,maxiter,tolerance_max,windowL,rhos[j],weights,V,V0,Eta,Eta0,approx) + end + end + + for j in 1:length(rhos) + if(xs[j]==Inf) + @printf(stderr,"max search failed for rho=%e\n",rhos[j]) + else + @printf("% .15e % .15e % .15e\n",rhos[j],xs[j],fs[j]) + end + end +end + + +# Fourier transform of 2 pt correlation function +function easyeq_2pt_fourier( + kmin::Float64, + kmax::Float64, + nk::Int64, + minlrho_init::Float64, + nlrho_init::Int64, + order::Int64, + windowL::Float64, + rho::Float64, + a0::Float64, + v::Function, + maxiter::Int64, + tolerance::Float64, + approx::Easyeq_approx +) + # compute gaussian quadrature weights + weights=gausslegendre(order) + + # compute u + (u,E,err)= easyeq_compute_u_prevrho_error([rho],minlrho_init,nlrho_init,a0,order,v,maxiter,tolerance,weights,approx) + + # compute useful terms + (V,V0)=easyeq_init_v(weights,v) + (Eta,Eta0)=easyeq_init_H(weights,v) + (E,S,A,T,B,X)=easyeq_ESATBX(rho*u,V,V0,Eta,Eta0,weights,rho,approx) + + # compute C2 + for j in 1:nk + k=kmin+(kmax-kmin)/nk*j + C2=easyeq_C2_fourier(k,u,windowL,rho,weights,Eta,Eta0,approx,E,S,A,T,B,X) + @printf("% .15e % .15e\n",k,C2) + end +end + +# maximum of Fourier transform of 2 point correlation function +function easyeq_2pt_fourier_max( + dk::Float64, + k0::Float64, # initial guess is k0*rho^(1/3) + minlrho_init::Float64, + nlrho_init::Int64, + order::Int64, + windowL::Float64, + rho::Float64, + a0::Float64, + v::Function, + maxstep::Float64, + maxiter::Int64, + tolerance::Float64, + tolerance_max::Float64, + approx::Easyeq_approx +) + # compute gaussian quadrature weights + weights=gausslegendre(order) + + # compute u + (u,E,err)= easyeq_compute_u_prevrho_error([rho],minlrho_init,nlrho_init,a0,order,v,maxiter,tolerance,weights,approx) + + # compute useful terms + (V,V0)=easyeq_init_v(weights,v) + (Eta,Eta0)=easyeq_init_H(weights,v) + + (k,f)=easyeq_C2_fourier_max(u,k0*rho^(1/3),dk,maxstep,maxiter,tolerance_max,windowL,rho,weights,V,V0,Eta,Eta0,approx) + + if(k==Inf) + @printf(stderr,"max search failed for rho=%e\n",rho) + else + @printf("% .15e % .15e\n",k,f) + end +end + +# maximum of Fourier transform of 2 point correlation function as a function of rho +function easyeq_2pt_fourier_max_rho( + rhos::Array{Float64,1}, + dk::Float64, + k0::Float64, # initial guess is k0*rho^(1/3) + minlrho_init::Float64, + nlrho_init::Int64, + order::Int64, + windowL::Float64, + a0::Float64, + v::Function, + maxstep::Float64, + maxiter::Int64, + tolerance::Float64, + tolerance_max::Float64, + approx::Easyeq_approx +) + # compute gaussian quadrature weights + weights=gausslegendre(order) + + # compute u + (us,es,errs)= easyeq_compute_u_prevrho_error(rhos,minlrho_init,nlrho_init,a0,order,v,maxiter,tolerance,weights,approx) + + # compute useful terms + (V,V0)=easyeq_init_v(weights,v) + (Eta,Eta0)=easyeq_init_H(weights,v) + + # save result from each task + ks=Array{Float64,1}(undef,length(rhos)) + fs=Array{Float64,1}(undef,length(rhos)) + + # spawn workers + work=spawn_workers(length(rhos)) + + count=0 + # for each worker + @sync for p in 1:length(work) + # for each task + @async for j in work[p] + count=count+1 + if count>=length(work) + progress(count,length(rhos),10000) + end + # run the task + (ks[j],fs[j])=remotecall_fetch(easyeq_C2_fourier_max,workers()[p],us[j],k0*rhos[j]^(1/3),dk,maxstep,maxiter,tolerance_max,windowL,rhos[j],weights,V,V0,Eta,Eta0,approx) + end + end + + for j in 1:length(rhos) + if(ks[j]==Inf) + @printf(stderr,"max search failed for rho=%e\n",rhos[j]) + else + @printf("% .15e % .15e % .15e\n",rhos[j],ks[j],fs[j]) + end + end +end + +# momentum distribution +function easyeq_momentum_distribution( + kmin::Float64, + kmax::Float64, + minlrho_init::Float64, + nlrho_init::Int64, + order::Int64, + windowL::Float64, #L=windowL/k^2 + rho::Float64, + a0::Float64, + v::Function, + maxiter::Int64, + tolerance::Float64, + approx::Easyeq_approx +) + # compute gaussian quadrature weights + weights=gausslegendre(order) + + # compute u + (u,E,err)= easyeq_compute_u_prevrho_error([rho],minlrho_init,nlrho_init,a0,order,v,maxiter,tolerance,weights,approx) + + # compute useful terms + (V,V0)=easyeq_init_v(weights,v) + (Eta,Eta0)=easyeq_init_H(weights,v) + (E,S,A,T,B,X)=easyeq_ESATBX(rho*u,V,V0,Eta,Eta0,weights,rho,approx) + # dXi/dlambda without the delta function and u + dXidlambda=(dotPhi(B.*T./((X.+1).^2))./(2*(X.+1))+B.*T./(2*(X.+1).^3).*dotdPhi(B.*T./(X.+1).^2))./A*(16*pi^3) + + dXi=inv(easyeq_dXi(rho*u,Eta,Eta0,weights,rho,approx,E,S,A,T,B,X)) + + # compute momentum distribution + for j in 1:order + q=(1-weights[1][j])/(1+weights[1][j]) + # drop if not in k interval + if q<kmin || q>kmax + continue + end + + # delta(k_i,q) + delta=Array{Float64,1}(undef,order) + L=windowL/q^2 + for i in 1:order + k=(1-weights[1][i])/(1+weights[1][i]) + delta[i]=(1.)/(2*k*q*L)*(gaussian(k-q,(1.)/L)-gaussian(k+q,(1.)/L)) + end + + # du/dlambda + du=-dXi*(dXidlambda.*delta*u[j]) + # rescale u + du=du/rho + + # compute M + M=-integrate_legendre_sampled(y->(1-y)/y^3,Eta0.*du,0.,1.,weights)*rho/(16*pi^3) + + @printf("% .15e % .15e\n",q,M) + end +end + + +# initialize u from scattering solution +@everywhere function easyeq_init_u( + a0::Float64, + order::Int64, + weights::Tuple{Array{Float64,1},Array{Float64,1}} +) u=zeros(Float64,order) for j in 1:order # transformed k @@ -150,8 +517,100 @@ end return u end +# compute u for an array of rhos +# use scattering solution for the first one, and the previous rho for the others +@everywhere function easyeq_compute_u_rho( + rhos::Array{Float64,1}, + a0::Float64, + order::Int64, + v::Function, + maxiter::Int64, + tolerance::Float64, + weights::Tuple{Array{Float64,1},Array{Float64,1}}, + approx::Easyeq_approx +) + us=Array{Array{Float64,1}}(undef,length(rhos)) + es=Array{Float64,1}(undef,length(rhos)) + errs=Array{Float64,1}(undef,length(rhos)) + + (us[1],es[1],errs[1])=easyeq_hatu(easyeq_init_u(a0,order,weights),order,rhos[1],v,maxiter,tolerance,weights,approx) + for j in 2:length(rhos) + (us[j],es[j],errs[j])=easyeq_hatu(us[j-1],order,rhos[j],v,maxiter,tolerance,weights,approx) + end + + return (us,es,errs) +end + +# compute u for an array of rhos +# start from a smaller rho and work up to rhos[1] +@everywhere function easyeq_compute_u_prevrho( + rhos::Array{Float64,1}, + minlrho_init::Float64, + nlrho_init::Int64, + a0::Float64, + order::Int64, + v::Function, + maxiter::Int64, + tolerance::Float64, + weights::Tuple{Array{Float64,1},Array{Float64,1}}, + approx::Easyeq_approx +) + + # only work up to rhos[1] if nlrho_init>0 + if nlrho_init>0 + rhos_init=Array{Float64,1}(undef,nlrho_init) + for j in 0:nlrho_init-1 + rhos_init[j+1]=(nlrho_init==1 ? 10^minlrho_init : 10^(minlrho_init+(log10(rhos[1])-minlrho_init)/(nlrho_init-1)*j)) + end + append!(rhos_init,rhos) + # start from rhos[1] if nlrho_init=0 + else + rhos_init=rhos + end + (us,es,errs)=easyeq_compute_u_rho(rhos_init,a0,order,v,maxiter,tolerance,weights,approx) + + # return a single value if there was a single input + if length(rhos)==1 + return (us[nlrho_init+1],es[nlrho_init+1],errs[nlrho_init+1]) + else + return (us[nlrho_init+1:length(us)],es[nlrho_init+1:length(es)],errs[nlrho_init+1:length(errs)]) + end +end +# with error message if the computation failed to be accurate enough +@everywhere function easyeq_compute_u_prevrho_error( + rhos::Array{Float64,1}, + minlrho_init::Float64, + nlrho_init::Int64, + a0::Float64, + order::Int64, + v::Function, + maxiter::Int64, + tolerance::Float64, + weights::Tuple{Array{Float64,1},Array{Float64,1}}, + approx::Easyeq_approx +) + (us,es,errs)=easyeq_compute_u_prevrho(rhos,minlrho_init,nlrho_init,a0,order,v,maxiter,tolerance,weights,approx) + # check errs + for j in 1:length(errs) + if errs[j]>tolerance + print(stderr,"warning: computation of u failed for rho=",rhos[j],"\n") + end + end + return (us,es,errs) +end + + # \hat u(k) computed using Newton -@everywhere function easyeq_hatu(u0,order,rho,v,maxiter,tolerance,weights,approx) +@everywhere function easyeq_hatu( + u0::Array{Float64,1}, + order::Int64, + rho::Float64, + v::Function, + maxiter::Int64, + tolerance::Float64, + weights::Tuple{Array{Float64,1},Array{Float64,1}}, + approx::Easyeq_approx +) # initialize V and Eta (V,V0)=easyeq_init_v(weights,v) (Eta,Eta0)=easyeq_init_H(weights,v) @@ -162,7 +621,8 @@ end # iterate err=Inf for i in 1:maxiter-1 - new=u-inv(easyeq_dXi(u,V,V0,Eta,Eta0,weights,rho,approx))*easyeq_Xi(u,V,V0,Eta,Eta0,weights,rho,approx) + (E,S,A,T,B,X)=easyeq_ESATBX(u,V,V0,Eta,Eta0,weights,rho,approx) + new=u-inv(easyeq_dXi(u,Eta,Eta0,weights,rho,approx,E,S,A,T,B,X))*easyeq_Xi(u,order,S,A,T,B,X) err=norm(new-u)/norm(u) if(err<tolerance) @@ -176,15 +636,23 @@ end end # \Eta -@everywhere function easyeq_H(x,t,weights,v) - 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) +@everywhere function easyeq_H( + x::Float64, + t::Float64, + weights::Tuple{Array{Float64,1},Array{Float64,1}}, + v::Function +) + 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 # initialize V -@everywhere function easyeq_init_v(weights,v) +@everywhere function easyeq_init_v( + weights::Tuple{Array{Float64,1},Array{Float64,1}}, + v::Function +) order=length(weights[1]) - V=Array{Float64}(undef,order) - V0=v(0) + V=Array{Float64,1}(undef,order) + V0=v(0.) for i in 1:order k=(1-weights[1][i])/(1+weights[1][i]) V[i]=v(k) @@ -193,29 +661,58 @@ end end # initialize Eta -@everywhere function easyeq_init_H(weights,v) +@everywhere function easyeq_init_H( + weights::Tuple{Array{Float64,1},Array{Float64,1}}, + v::Function +) order=length(weights[1]) - Eta=Array{Array{Float64}}(undef,order) - Eta0=Array{Float64}(undef,order) + Eta=Array{Array{Float64,1},1}(undef,order) + Eta0=Array{Float64,1}(undef,order) for i in 1:order k=(1-weights[1][i])/(1+weights[1][i]) - Eta[i]=Array{Float64}(undef,order) + Eta[i]=Array{Float64,1}(undef,order) for j in 1:order y=(weights[1][j]+1)/2 Eta[i][j]=easyeq_H(k,(1-y)/y,weights,v) end y=(weights[1][i]+1)/2 - Eta0[i]=easyeq_H(0,(1-y)/y,weights,v) + Eta0[i]=easyeq_H(0.,(1-y)/y,weights,v) end return(Eta,Eta0) end # Xi(u) -@everywhere function easyeq_Xi(u,V,V0,Eta,Eta0,weights,rho,approx) +@everywhere function easyeq_Xi( + u::Array{Float64,1}, + order::Int64, + S::Array{Float64,1}, + A::Array{Float64,1}, + T::Array{Float64,1}, + B::Array{Float64,1}, + X::Array{Float64,1} +) + return u.-T./(2*(X.+1)).*dotPhi(B.*T./(X.+1).^2) +end + +# compute E,S,A,T,B,X +@everywhere function easyeq_ESATBX( + u::Array{Float64,1}, + V::Array{Float64,1}, + V0::Float64, + Eta::Array{Array{Float64,1},1}, + Eta0::Array{Float64,1}, + weights::Tuple{Array{Float64,1},Array{Float64,1}}, + rho::Float64, + approx::Easyeq_approx +) order=length(weights[1]) # init - out=zeros(Float64,order) + S=zeros(Float64,order) + A=zeros(Float64,order) + T=zeros(Float64,order) + B=zeros(Float64,order) + X=zeros(Float64,order) # compute E before running the loop E=easyeq_en(u,V0,Eta0,rho,weights) @@ -224,188 +721,359 @@ end # k_i k=(1-weights[1][i])/(1+weights[1][i]) # S_i - S=V[i]-1/(rho*(2*pi)^3)*integrate_legendre_sampled(y->(1-y)/y^3,Eta[i].*u,0,1,weights) + S[i]=V[i]-1/(rho*(2*pi)^3)*integrate_legendre_sampled(y->(1-y)/y^3,Eta[i].*u,0.,1.,weights) # A_K,i - A=0. + A[i]=0. if approx.bK!=0. - A+=approx.bK*S + A[i]+=approx.bK*S[i] end if approx.bK!=1. - A+=(1-approx.bK)*E + A[i]+=(1-approx.bK)*E end - # T + # T_i if approx.bK==1. - T=1. + T[i]=1. else - T=S/A + T[i]=S[i]/A[i] end - # B + # B_i if approx.bK==approx.bL - B=1. + B[i]=1. else - B=(approx.bL*S+(1-approx.bL*E))/(approx.bK*S+(1-approx.bK*E)) + B[i]=(approx.bL*S[i]+(1-approx.bL*E))/(approx.bK*S[i]+(1-approx.bK*E)) end # X_i - X=k^2/(2*A*rho) - - # U_i - out[i]=u[i]-T/(2*(X+1))*Phi(B*T/(X+1)^2) + X[i]=k^2/(2*A[i]*rho) end - return out + return (E,S,A,T,B,X) end # derivative of Xi -@everywhere function easyeq_dXi(u,V,V0,Eta,Eta0,weights,rho,approx) +@everywhere function easyeq_dXi( + u::Array{Float64,1}, + Eta::Array{Array{Float64,1},1}, + Eta0::Array{Float64,1}, + weights::Tuple{Array{Float64,1},Array{Float64,1}}, + rho::Float64, + approx::Easyeq_approx, + E::Float64, + S::Array{Float64,1}, + A::Array{Float64,1}, + T::Array{Float64,1}, + B::Array{Float64,1}, + X::Array{Float64,1} +) order=length(weights[1]) # init out=zeros(Float64,order,order) - # compute E before the loop - E=easyeq_en(u,V0,Eta0,rho,weights) - for i in 1:order # k_i k=(1-weights[1][i])/(1+weights[1][i]) - # S_i - S=V[i]-1/(rho*(2*pi)^3)*integrate_legendre_sampled(y->(1-y)/y^3,Eta[i].*u,0,1,weights) - - # A_K,i - A=0. - if approx.bK!=0. - A+=approx.bK*S - end - if approx.bK!=1. - A+=(1-approx.bK)*E - end - - # T - if approx.bK==1. - T=1. - else - T=S/A - end - - # B - if approx.bK==approx.bL - B=1. - else - B=(approx.bL*S+(1-approx.bL*E))/(approx.bK*S+(1-approx.bK*E)) - end - - # X_i - X=k^2/(2*A*rho) for j in 1:order y=(weights[1][j]+1)/2 dS=-1/rho*(1-y)*Eta[i][j]/(2*(2*pi)^3*y^3)*weights[2][j] dE=-1/rho*(1-y)*Eta0[j]/(2*(2*pi)^3*y^3)*weights[2][j] + dU=(i==j ? 1. : 0.) - # dA - dA=0. - if approx.bK!=0. - dA+=approx.bK*dS - end - if approx.bK!=1. - dA+=(1-approx.bK)*dE - end + out[i,j]=easyeq_dXi_of_dSdEdU(k,dS,dE,dU,E,S[i],A[i],T[i],B[i],X[i],rho,approx) + end + end + + return out +end - # dT - if approx.bK==1. - dT=0. - else - dT=(1-approx.bK)*(E*dS-S*dE)/A^2 - end +# dXi given dS, dE and dU +@everywhere function easyeq_dXi_of_dSdEdU( + k::Float64, + dS::Float64, + dE::Float64, + dU::Float64, + E::Float64, + S::Float64, + A::Float64, + T::Float64, + B::Float64, + X::Float64, + rho::Float64, + approx::Easyeq_approx +) + # dA + dA=0. + if approx.bK!=0. + dA+=approx.bK*dS + end + if approx.bK!=1. + dA+=(1-approx.bK)*dE + end - # dB - if approx.bK==approx.bL - dB=0. - else - dB=(approx.bL*(1-approx.bK)-approx.bK*(1-approx.bL))*(E*dS-S*dE)/(approx.bK*S+(1-approx.bK*E))^2 - end + # dT,dB + # nothing to do if bK=bL=1 + if approx.bK!=1. || approx.bK!=approx.bL + dB=(E*dS-S*dE)/A^2 + end + if approx.bK==1. + dT=0. + else + dT=(1-approx.bK)*dB + end + if approx.bK==approx.bL + dB=0. + else + dB=(approx.bL*(1-approx.bK)-approx.bK*(1-approx.bL))*dB + end - dX=-k^2/(2*A^2*rho)*dA + dX=-k^2/(2*A^2*rho)*dA - out[i,j]=(i==j ? 1 : 0)-(dT-T*dX/(X+1))/(2*(X+1))*Phi(B*T/(X+1)^2)-T/(2*(X+1)^3)*(B*dT+T*dB-2*B*T*dX/(X+1))*dPhi(B*T/(X+1)^2) - end + return dU-(dT-T*dX/(X+1))/(2*(X+1))*Phi(B*T/(X+1)^2)-T/(2*(X+1)^3)*(B*dT+T*dB-2*B*T*dX/(X+1))*dPhi(B*T/(X+1)^2) +end + +# derivative of Xi with respect to mu +@everywhere function easyeq_dXidmu( + u::Array{Float64,1}, + order::Int64, + rho::Float64, + A::Array{Float64,1}, + T::Array{Float64,1}, + B::Array{Float64,1}, + X::Array{Float64,1} +) + # init + out=zeros(Float64,order) + for i in 1:order + out[i]=T[i]/(2*rho*A[i]*(X[i]+1)^2)*Phi(B[i]*T[i]/(X[i]+1)^2)+B[i]*T[i]^2/(rho*A[i]*(X[i]+1)^4)*dPhi(B[i]*T[i]/(X[i]+1)^2) end return out end -# derivative of Xi with respect to mu -@everywhere function easyeq_dXidmu(u,V,V0,Eta,Eta0,weights,rho,approx) +# energy +@everywhere function easyeq_en( + u::Array{Float64,1}, + V0::Float64, + Eta0::Array{Float64,1}, + rho::Float64, + weights::Tuple{Array{Float64,1},Array{Float64,1}} +) + return V0-1/(rho*(2*pi)^3)*integrate_legendre_sampled(y->(1-y)/y^3,Eta0.*u,0.,1.,weights) +end + +# condensate fraction +@everywhere function easyeq_eta( + u::Array{Float64,1}, + order::Int64, + rho::Float64, + v::Function, + maxiter::Int64, + tolerance::Float64, + weights::Tuple{Array{Float64,1},Array{Float64,1}}, + approx::Easyeq_approx +) + (V,V0)=easyeq_init_v(weights,v) + (Eta,Eta0)=easyeq_init_H(weights,v) + + (E,S,A,T,B,X)=easyeq_ESATBX(rho*u,V,V0,Eta,Eta0,weights,rho,approx) + du=-inv(easyeq_dXi(rho*u,Eta,Eta0,weights,rho,approx,E,S,A,T,B,X))*easyeq_dXidmu(rho*u,order,rho,A,T,B,X) + + eta=-1/(2*(2*pi)^3)*integrate_legendre_sampled(y->(1-y)/y^3,Eta0.*du,0.,1.,weights) + + return eta +end + +# inverse Fourier transform +@everywhere function easyeq_u_x( + x::Float64, + u::Array{Float64,1}, + weights::Tuple{Array{Float64,1},Array{Float64,1}} +) + order=length(weights[1]) + out=integrate_legendre_sampled(y->(1-y)/y^3*sin(x*(1-y)/y)/x/(2*pi^2),u,0.,1.,weights) + return out +end + + +# correlation function +@everywhere function easyeq_C2( + x::Float64, + u::Array{Float64,1}, + windowL::Float64, + rho::Float64, + weights::Tuple{Array{Float64,1},Array{Float64,1}}, + Eta::Array{Array{Float64,1},1}, + Eta0::Array{Float64,1}, + approx::Easyeq_approx, + E::Float64, + S::Array{Float64,1}, + A::Array{Float64,1}, + T::Array{Float64,1}, + B::Array{Float64,1}, + X::Array{Float64,1} +) + g=(r,x)->(r>0. ? sin(r*x)/(r*x)*hann(r,windowL) : hann(r,windowL)) + (dEta,dEta0)=easyeq_init_H(weights,k->g(k,x)) + du=easyeq_dudv(g,x,u,rho,weights,Eta,Eta0,dEta,dEta0,approx,E,S,A,T,B,X) + + return rho^2*(1-integrate_legendre_sampled(y->(1-y)/y^3,dEta0.*u+Eta0.*du,0.,1.,weights)/(8*pi^3)) +end + +# derivative of u with respect to v in direction g +@everywhere function easyeq_dudv( + g::Function,# should be of the form g(k,x) where x is a parameter + x::Float64, + u::Array{Float64,1}, + rho::Float64, + weights::Tuple{Array{Float64,1},Array{Float64,1}}, + Eta::Array{Array{Float64,1},1}, + Eta0::Array{Float64,1}, + dEta::Array{Array{Float64,1},1}, + dEta0::Array{Float64,1}, + approx::Easyeq_approx, + E::Float64, + S::Array{Float64,1}, + A::Array{Float64,1}, + T::Array{Float64,1}, + B::Array{Float64,1}, + X::Array{Float64,1} +) + # initialize dV and dEta + (dV,dV0)=easyeq_init_v(weights,k->g(k,x)) + + du=-inv(easyeq_dXi(rho*u,Eta,Eta0,weights,rho,approx,E,S,A,T,B,X))*easyeq_dXidv(rho*u,dV,dV0,dEta,dEta0,weights,rho,approx,E,S,A,T,B,X) + # rescale rho + du=du/rho + + return du +end + +# derivative of Xi with respect to potential +@everywhere function easyeq_dXidv( + u::Array{Float64,1}, + dv::Array{Float64,1}, + dv0::Float64, + dEta::Array{Array{Float64,1},1}, + dEta0::Array{Float64,1}, + weights::Tuple{Array{Float64,1},Array{Float64,1}}, + rho::Float64, + approx::Easyeq_approx, + E::Float64, + S::Array{Float64,1}, + A::Array{Float64,1}, + T::Array{Float64,1}, + B::Array{Float64,1}, + X::Array{Float64,1} +) + order=length(weights[1]) # init out=zeros(Float64,order) - # compute E before running the loop - E=easyeq_en(u,V0,Eta0,rho,weights) - for i in 1:order # k_i k=(1-weights[1][i])/(1+weights[1][i]) - # S_i - S=V[i]-1/(rho*(2*pi)^3)*integrate_legendre_sampled(y->(1-y)/y^3,Eta[i].*u,0,1,weights) - # A_K,i - A=0. - if approx.bK!=0. - A+=approx.bK*S - end - if approx.bK!=1. - A+=(1-approx.bK)*E - end - - # T - if approx.bK==1. - T=1. - else - T=S/A - end - - # B - if approx.bK==approx.bL - B=1. - else - B=(approx.bL*S+(1-approx.bL*E))/(approx.bK*S+(1-approx.bK*E)) + dS=dv[i] + dE=dv0 + for j in 1:order + y=(weights[1][j]+1)/2 + dS+=-1/rho*(1-y)*u[j]*dEta[i][j]/(2*(2*pi)^3*y^3)*weights[2][j] + dE+=-1/rho*(1-y)*u[j]*dEta0[j]/(2*(2*pi)^3*y^3)*weights[2][j] end - # X_i - X=k^2/(2*A*rho) - - out[i]=T/(2*rho*A*(X+1)^2)*Phi(B*T/(X+1)^2)+B*T^2/(rho*A*(X+1)^4)*dPhi(B*T/(X+1)^2) + out[i]=easyeq_dXi_of_dSdEdU(k,dS,dE,0.,E,S[i],A[i],T[i],B[i],X[i],rho,approx) end return out end -# energy -@everywhere function easyeq_en(u,V0,Eta0,rho,weights) - return V0-1/(rho*(2*pi)^3)*integrate_legendre_sampled(y->(1-y)/y^3,Eta0.*u,0,1,weights) +# maximum of 2 point correlation function +@everywhere function easyeq_C2_max( + u::Array{Float64,1}, + x0::Float64, + dx::Float64, + maxstep::Float64, + maxiter::Int64, + tolerance::Float64, + windowL::Float64, + rho::Float64, + weights::Tuple{Array{Float64,1},Array{Float64,1}}, + V::Array{Float64,1}, + V0::Float64, + Eta::Array{Array{Float64,1},1}, + Eta0::Array{Float64,1}, + approx::Easyeq_approx +) + # compute some useful terms + (E,S,A,T,B,X)=easyeq_ESATBX(rho*u,V,V0,Eta,Eta0,weights,rho,approx) + + (x,f)=newton_maximum(y->easyeq_C2(y,u,windowL,rho,weights,Eta,Eta0,approx,E,S,A,T,B,X),x0,dx,maxiter,tolerance,maxstep) + + return(x,f) end -# condensate fraction -@everywhere function easyeq_eta(u,order,rho,v,maxiter,tolerance,weights,approx) - (V,V0)=easyeq_init_v(weights,v) - (Eta,Eta0)=easyeq_init_H(weights,v) - - du=-inv(easyeq_dXi(rho*u,V,V0,Eta,Eta0,weights,rho,approx))*easyeq_dXidmu(rho*u,V,V0,Eta,Eta0,weights,rho,approx) +# Fourier transform of correlation function +@everywhere function easyeq_C2_fourier( + q::Float64, + u::Array{Float64,1}, + windowL::Float64, + rho::Float64, + weights::Tuple{Array{Float64,1},Array{Float64,1}}, + Eta::Array{Array{Float64,1},1}, + Eta0::Array{Float64,1}, + approx::Easyeq_approx, + E::Float64, + S::Array{Float64,1}, + A::Array{Float64,1}, + T::Array{Float64,1}, + B::Array{Float64,1}, + X::Array{Float64,1} +) + # direction in which to differentiate u + g=(r,x)->(r>0. ? (1.)/(2*x*r*windowL)*(gaussian(x-r,(1.)/windowL)-gaussian(x+r,(1.)/windowL)) : gaussian(x,(1.)/windowL)) + (dEta,dEta0)=easyeq_init_H(weights,k->g(k,q)) + du=easyeq_dudv(g,q,u,rho,weights,Eta,Eta0,dEta,dEta0,approx,E,S,A,T,B,X) + + return rho^2*(-integrate_legendre_sampled(y->(1-y)/y^3,dEta0.*u+Eta0.*du,0.,1.,weights)/(8*pi^3)) +end - eta=-1/(2*(2*pi)^3)*integrate_legendre_sampled(y->(1-y)/y^3,Eta0.*du,0,1,weights) - return eta +# maximum of Fourier transform of 2 point correlation function +@everywhere function easyeq_C2_fourier_max( + u::Array{Float64,1}, + k0::Float64, + dk::Float64, + maxstep::Float64, + maxiter::Int64, + tolerance::Float64, + windowL::Float64, + rho::Float64, + weights::Tuple{Array{Float64,1},Array{Float64,1}}, + V::Array{Float64,1}, + V0::Float64, + Eta::Array{Array{Float64,1},1}, + Eta0::Array{Float64,1}, + approx::Easyeq_approx +) + # compute some useful terms + (E,S,A,T,B,X)=easyeq_ESATBX(rho*u,V,V0,Eta,Eta0,weights,rho,approx) + + (k,f)=newton_maximum(y->easyeq_C2_fourier(y,u,windowL,rho,weights,Eta,Eta0,approx,E,S,A,T,B,X),k0,dk,maxiter,tolerance,maxstep) + + return (k,f) end -# inverse Fourier transform -@everywhere function easyeq_u_x(x,u,weights) - order=length(weights[1]) - out=integrate_legendre_sampled(y->(1-y)/y^3*sin(x*(1-y)/y)/x/(2*pi^2),u,0,1,weights) - return out + +@everywhere function barEta( + q::Float64, + y::Float64, + t::Float64 +) + return (t>=abs(y-q) && t<=y+q ? pi/y/q : 0.) end diff --git a/src/integration.jl b/src/integration.jl index 9be4641..223d6cc 100644 --- a/src/integration.jl +++ b/src/integration.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,7 +13,12 @@ ## limitations under the License. # approximate \int_a^b f using Gauss-Legendre quadratures -@everywhere function integrate_legendre(f,a,b,weights) +@everywhere function integrate_legendre( + f::Function, + a::Float64, + b::Float64, + weights::Tuple{Array{Float64,1},Array{Float64,1}} +) out=0 for i in 1:length(weights[1]) out+=(b-a)/2*weights[2][i]*f((b-a)/2*weights[1][i]+(b+a)/2) @@ -21,7 +26,13 @@ return out end # \int f*g where g is sampled at the Legendre nodes -@everywhere function integrate_legendre_sampled(f,g,a,b,weights) +@everywhere function integrate_legendre_sampled( + f::Function, + g::Array{Float64,1}, + a::Float64, + b::Float64, + weights::Tuple{Array{Float64,1},Array{Float64,1}} +) out=0 for i in 1:length(weights[1]) out+=(b-a)/2*weights[2][i]*f((b-a)/2*weights[1][i]+(b+a)/2)*g[i] @@ -31,7 +42,12 @@ end # approximate \int_a^b f/sqrt((b-x)(x-a)) using Gauss-Chebyshev quadratures -@everywhere function integrate_chebyshev(f,a,b,N) +@everywhere function integrate_chebyshev( + f::Function, + a::Float64, + b::Float64, + N::Int64 +) out=0 for i in 1:N out=out+pi/N*f((b-a)/2*cos((2*i-1)/(2*N)*pi)+(b+a)/2) @@ -40,7 +56,11 @@ end end # approximate \int_0^\infty dr f(r)*exp(-a*r) using Gauss-Chebyshev quadratures -@everywhere function integrate_laguerre(f,a,weights_gL) +@everywhere function integrate_laguerre( + f::Function, + a::Float64, + weights_gL::Tuple{Array{Float64,1},Array{Float64,1}} +) out=0. for i in 1:length(weights_gL[1]) out+=1/a*f(weights_gL[1][i]/a)*weights_gL[2][i] @@ -49,10 +69,28 @@ end end # Hann window -@everywhere function hann(x,L) +@everywhere function hann( + x::Float64, + L::Float64 +) if abs(x)<L/2 return cos(pi*x/L)^2 else return 0. end end +# Fourier transform (in 3d) +@everywhere function hann_fourier( + k::Float64, + L::Float64 +) + return L^2*4*pi^3/k*(((k*L)^3-4*k*L*pi^2)*cos(k*L/2)-2*(3*(k*L)^2-4*pi^2)*sin(k*L/2))/((k*L)^3-4*k*L*pi^2)^2 +end + +# normalized Gaussian (in 3d) +@everywhere function gaussian( + k::Float64, + L::Float64 +) + return exp(-k^2/(2*L))/sqrt(8*pi^3*L^3) +end diff --git a/src/interpolation.jl b/src/interpolation.jl index fa3bcdb..066bc20 100644 --- a/src/interpolation.jl +++ b/src/interpolation.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. @@ -14,7 +14,11 @@ # linear interpolation: given vectors x,y, compute a linear interpolation for y(x0) # assume x is ordered -@everywhere function linear_interpolation(x0,x,y) +@everywhere function linear_interpolation( + x0::Float64, + x::Array{Float64,1}, + y::Array{Float64,1} +) # if x0 is beyond all x's, then return the corresponding boundary value. if x0>x[length(x)] return y[length(y)] @@ -28,7 +32,12 @@ # interpolate return y[i]+(y[i+1]-y[i])*(x0-x[i])/(x[i+1]-x[i]) end -@everywhere function bracket(x0,x,a,b) +@everywhere function bracket( + x0::Float64, + x::Array{Float64,1}, + a::Int64, + b::Int64 +) i=floor(Int64,(a+b)/2) if x0<x[i] return bracket(x0,x,a,i) @@ -41,15 +50,18 @@ end # polynomial interpolation of a family of points -@everywhere function poly_interpolation(x,y) +@everywhere function poly_interpolation( + x::Array{Float64,1}, + y::Array{Float64,1} +) # init for recursion - rec=Array{Polynomial{Float64}}(undef,length(x)) + rec=Array{Polynomial{Float64},1}(undef,length(x)) for i in 1:length(x) rec[i]=Polynomial([1.]) end # compute \prod (x-x_i) - poly_interpolation_rec(rec,x,1,length(x)) + poly_interpolation_rec(rec,x,1.,length(x)) # sum terms together out=0. @@ -60,7 +72,12 @@ end return out end # recursive helper function -@everywhere function poly_interpolation_rec(out,x,a,b) +@everywhere function poly_interpolation_rec( + out::Array{Float64,1}, + x::Array{Float64,1}, + a::Float64, + b::Float64 +) if a==b return end @@ -91,7 +108,10 @@ end return end ## the following does the same, but has complexity N^2, the function above has N*log(N) -#@everywhere function poly_interpolation(x,y) +#@everywhere function poly_interpolation( +# x::Array{Float64,1}, +# y::Array{Float64,1} +#) # out=Polynomial([0.]) # for i in 1:length(x) # prod=Polynomial([1.]) diff --git a/src/main.jl b/src/main.jl index 382fe6b..28fc2be 100644 --- a/src/main.jl +++ b/src/main.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. @@ -22,6 +22,8 @@ include("chebyshev.jl") include("integration.jl") include("interpolation.jl") include("tools.jl") +include("multithread.jl") +include("optimization.jl") include("potentials.jl") include("print.jl") include("easyeq.jl") @@ -36,9 +38,6 @@ function main() rho=1e-6 e=1e-4 - # incrementally initialize Newton algorithm - nlrho_init=1 - # potential v=k->v_exp(k,1.) a0=a0_exp(1.) @@ -49,19 +48,29 @@ function main() v_param_e=1. # plot range when plotting in rho - minlrho=-6 - maxlrho=2 + # linear + minrho=1e-6 + maxrho=1e2 + nrho=0 + # logarithmic + minlrho=-6. + maxlrho=2. nlrho=100 - rhos=Array{Float64}(undef,0) + # list + rhos=Array{Float64,1}(undef,0) # plot range when plotting in e - minle=-6 - maxle=2 + minle=-6. + maxle=2. nle=100 - es=Array{Float64}(undef,0) + es=Array{Float64,1}(undef,0) # plot range when plotting in x - xmin=0 - xmax=100 + xmin=0. + xmax=100. nx=100 + # plot range when plotting in k + kmin=0. + kmax=10. + nk=100 # cutoffs tolerance=1e-11 @@ -77,8 +86,8 @@ function main() J=10 # starting rho from which to incrementally initialize Newton algorithm - # default must be set after reading rho, if not set explicitly - minlrho_init=nothing + minlrho_init=-6. + nlrho_init=0 # Hann window for Fourier transforms windowL=1e3 @@ -98,6 +107,18 @@ function main() anyeq_compleq_approx=Anyeq_approx(1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.) anyeq_approx=anyeq_bigeq_approx + # numerical approximations of derivatives + dx=1e-7 + dk=1e-7 + + # initial guess for 2pt_max + x0=1. + k0=1. + # maximum step in 2pt_max + maxstep=Inf + # tolerance for max search + tolerance_max=Inf + # read cli arguments (params,potential,method,savefile,command)=read_args(ARGS) @@ -109,8 +130,8 @@ function main() print(stderr,"error: could not read parameter '",param,"'.\n") exit(-1) end - lhs=terms[1] - rhs=terms[2] + lhs=string(terms[1]) + rhs=string(terms[2]) if lhs=="rho" rho=parse(Float64,rhs) elseif lhs=="minlrho_init" @@ -133,6 +154,12 @@ function main() maxlrho=parse(Float64,rhs) elseif lhs=="nlrho" nlrho=parse(Int64,rhs) + elseif lhs=="minrho" + minrho=parse(Float64,rhs) + elseif lhs=="maxrho" + maxrho=parse(Float64,rhs) + elseif lhs=="nrho" + nrho=parse(Int64,rhs) elseif lhs=="es" es=parse_list(rhs) elseif lhs=="minle" @@ -147,6 +174,12 @@ function main() xmax=parse(Float64,rhs) elseif lhs=="nx" nx=parse(Int64,rhs) + elseif lhs=="kmin" + kmin=parse(Float64,rhs) + elseif lhs=="kmax" + kmax=parse(Float64,rhs) + elseif lhs=="nk" + nk=parse(Int64,rhs) elseif lhs=="P" P=parse(Int64,rhs) elseif lhs=="N" @@ -155,6 +188,18 @@ function main() J=parse(Int64,rhs) elseif lhs=="window_L" windowL=parse(Float64,rhs) + elseif lhs=="dx" + dx=parse(Float64,rhs) + elseif lhs=="x0" + x0=parse(Float64,rhs) + elseif lhs=="dk" + dk=parse(Float64,rhs) + elseif lhs=="k0" + k0=parse(Float64,rhs) + elseif lhs=="maxstep" + maxstep=parse(Float64,rhs) + elseif lhs=="tolerance_max" + tolerance_max=parse(Float64,rhs) elseif lhs=="aK" anyeq_approx.aK=parse(Float64,rhs) elseif lhs=="bK" @@ -242,37 +287,48 @@ function main() ## set parameters # rhos if length(rhos)==0 - rhos=Array{Float64}(undef,nlrho) - for j in 0:nlrho-1 - rhos[j+1]=(nlrho==1 ? 10^minlrho : 10^(minlrho+(maxlrho-minlrho)/(nlrho-1)*j)) + # linear only if nrho is specified + if nrho>0 + rhos=Array{Float64,1}(undef,nrho) + for j in 0:nrho-1 + rhos[j+1]=(nrho==1 ? minrho : minrho+(maxrho-minrho)/(nrho-1)*j) + end + else + rhos=Array{Float64,1}(undef,nlrho) + for j in 0:nlrho-1 + rhos[j+1]=(nlrho==1 ? 10^minlrho : 10^(minlrho+(maxlrho-minlrho)/(nlrho-1)*j)) + end end end + # es if length(es)==0 - es=Array{Float64}(undef,nle) + es=Array{Float64,1}(undef,nle) for j in 0:nle-1 es[j+1]=(nle==1 ? 10^minle : 10^(minle+(maxle-minle)/(nle-1)*j)) end end - # default minlrho_init - if (minlrho_init==nothing) - minlrho_init=log10(rho) - end - # splines - taus=Array{Float64}(undef,J+1) + taus=Array{Float64,1}(undef,J+1) for j in 0:J taus[j+1]=-1+2*j/J end + # tolerance_max + if tolerance_max==Inf + tolerance_max=tolerance + end + ## run command - if method=="easyeq" + if command=="scattering_length" + @printf("% .15e\n",a0) + elseif method=="easyeq" if command=="energy" easyeq_energy(minlrho_init,nlrho_init,order,rho,a0,v,maxiter,tolerance,easyeq_approx) # e(rho) elseif command=="energy_rho" - easyeq_energy_rho(rhos,order,a0,v,maxiter,tolerance,easyeq_approx) + easyeq_energy_rho(rhos,minlrho_init,nlrho_init,order,a0,v,maxiter,tolerance,easyeq_approx) # u(k) elseif command=="uk" easyeq_uk(minlrho_init,nlrho_init,order,rho,a0,v,maxiter,tolerance,easyeq_approx) @@ -285,7 +341,28 @@ function main() elseif command=="condensate_fraction" easyeq_condensate_fraction(minlrho_init,nlrho_init,order,rho,a0,v,maxiter,tolerance,easyeq_approx) elseif command=="condensate_fraction_rho" - easyeq_condensate_fraction_rho(rhos,order,a0,v,maxiter,tolerance,easyeq_approx) + easyeq_condensate_fraction_rho(rhos,minlrho_init,nlrho_init,order,a0,v,maxiter,tolerance,easyeq_approx) + # 2pt correlation + elseif command=="2pt" + easyeq_2pt(xmin,xmax,nx,minlrho_init,nlrho_init,order,windowL,rho,a0,v,maxiter,tolerance,easyeq_approx) + # max of 2pt correlation + elseif command=="2pt_max" + easyeq_2pt_max(dx,x0,minlrho_init,nlrho_init,order,windowL,rho,a0,v,maxstep,maxiter,tolerance,tolerance_max,easyeq_approx) + # max of 2pt correlation as a function of rho + elseif command=="2pt_max_rho" + easyeq_2pt_max_rho(rhos,dx,x0,minlrho_init,nlrho_init,order,windowL,a0,v,maxstep,maxiter,tolerance,tolerance_max,easyeq_approx) + # fourier transform of 2pt correlation + elseif command=="2pt_fourier" + easyeq_2pt_fourier(kmin,kmax,nk,minlrho_init,nlrho_init,order,windowL,rho,a0,v,maxiter,tolerance,easyeq_approx) + # max of fourier transform of 2pt correlation + elseif command=="2pt_fourier_max" + easyeq_2pt_fourier_max(dk,k0,minlrho_init,nlrho_init,order,windowL,rho,a0,v,maxstep,maxiter,tolerance,tolerance_max,easyeq_approx) + # max of 2pt correlation as a function of rho + elseif command=="2pt_fourier_max_rho" + easyeq_2pt_fourier_max_rho(rhos,dk,k0,minlrho_init,nlrho_init,order,windowL,a0,v,maxstep,maxiter,tolerance,tolerance_max,easyeq_approx) + # momentum distribution + elseif command=="momentum_distribution" + easyeq_momentum_distribution(kmin,kmax,minlrho_init,nlrho_init,order,windowL,rho,a0,v,maxiter,tolerance,easyeq_approx) else print(stderr,"unrecognized command '",command,"'.\n") exit(-1) @@ -316,7 +393,7 @@ function main() anyeq_energy(rho,minlrho_init,nlrho_init,taus,P,N,J,a0,v,maxiter,tolerance,anyeq_approx,savefile) # e(rho) elseif command=="energy_rho" - anyeq_energy_rho(rhos,taus,P,N,J,a0,v,maxiter,tolerance,anyeq_approx,savefile) + anyeq_energy_rho(rhos,minlrho_init,nlrho_init,taus,P,N,J,a0,v,maxiter,tolerance,anyeq_approx,savefile) elseif command=="energy_rho_init_prevrho" anyeq_energy_rho_init_prevrho(rhos,taus,P,N,J,a0,v,maxiter,tolerance,anyeq_approx,savefile) elseif command=="energy_rho_init_nextrho" @@ -331,12 +408,29 @@ function main() elseif command=="condensate_fraction" anyeq_condensate_fraction(rho,minlrho_init,nlrho_init,taus,P,N,J,a0,v,maxiter,tolerance,anyeq_approx,savefile) elseif command=="condensate_fraction_rho" - anyeq_condensate_fraction_rho(rhos,taus,P,N,J,a0,v,maxiter,tolerance,anyeq_approx,savefile) + anyeq_condensate_fraction_rho(rhos,minlrho_init,nlrho_init,taus,P,N,J,a0,v,maxiter,tolerance,anyeq_approx,savefile) # momentum distribution elseif command=="momentum_distribution" - anyeq_momentum_distribution(rho,minlrho_init,nlrho_init,taus,P,N,J,a0,v,maxiter,tolerance,anyeq_approx,savefile) + anyeq_momentum_distribution(kmin,kmax,rho,minlrho_init,nlrho_init,taus,P,N,J,windowL,a0,v,maxiter,tolerance,anyeq_approx,savefile) elseif command=="2pt" anyeq_2pt_correlation(minlrho_init,nlrho_init,taus,P,N,J,windowL,rho,a0,v,maxiter,tolerance,xmin,xmax,nx,anyeq_approx,savefile) + elseif command=="2pt_max" + anyeq_2pt_correlation_max(rho,minlrho_init,nlrho_init,dx,x0,maxstep,taus,P,N,J,windowL,a0,v,maxiter,tolerance,tolerance_max,anyeq_approx,savefile) + elseif command=="2pt_max_rho" + anyeq_2pt_correlation_max_rho(rhos,minlrho_init,nlrho_init,dx,x0,maxstep,taus,P,N,J,windowL,a0,v,maxiter,tolerance,tolerance_max,anyeq_approx,savefile) + elseif command=="2pt_fourier" + anyeq_2pt_correlation_fourier(minlrho_init,nlrho_init,taus,P,N,J,windowL,rho,a0,v,maxiter,tolerance,kmin,kmax,nk,anyeq_approx,savefile) + elseif command=="2pt_fourier_test" + anyeq_2pt_correlation_fourier_test(minlrho_init,nlrho_init,taus,P,N,J,windowL,rho,a0,v,maxiter,tolerance,xmax,kmin,kmax,nk,anyeq_approx,savefile) + elseif command=="2pt_fourier_max" + anyeq_2pt_correlation_fourier_max(rho,minlrho_init,nlrho_init,dk,k0,maxstep,taus,P,N,J,windowL,a0,v,maxiter,tolerance,tolerance_max,anyeq_approx,savefile) + elseif command=="2pt_fourier_max_rho" + anyeq_2pt_correlation_fourier_max_rho(rhos,minlrho_init,nlrho_init,dk,k0,maxstep,taus,P,N,J,windowL,a0,v,maxiter,tolerance,tolerance_max,anyeq_approx,savefile) + elseif command=="uncondensed_2pt" + anyeq_uncondensed_2pt_correlation(minlrho_init,nlrho_init,taus,P,N,J,windowL,rho,a0,v,maxiter,tolerance,xmin,xmax,nx,anyeq_approx,savefile) + # compressibility + elseif command=="compressibility_rho" + anyeq_compressibility_rho(rhos,minlrho_init,nlrho_init,taus,P,N,J,a0,v,maxiter,tolerance,anyeq_approx,savefile) else print(stderr,"unrecognized command: '",command,"'\n") exit(-1) @@ -360,9 +454,11 @@ function main() end # parse a comma separated list as an array of Float64 -function parse_list(str) +function parse_list( + str::String +) elems=split(str,",") - out=Array{Float64}(undef,length(elems)) + out=Array{Float64,1}(undef,length(elems)) for i in 1:length(elems) out[i]=parse(Float64,elems[i]) end @@ -370,7 +466,9 @@ function parse_list(str) end # read cli arguments -function read_args(ARGS) +function read_args( + ARGS +) # flag flag="" diff --git a/src/multithread.jl b/src/multithread.jl new file mode 100644 index 0000000..f61cdd7 --- /dev/null +++ b/src/multithread.jl @@ -0,0 +1,16 @@ +# split up 1...n among workers +function spawn_workers(n::Int64) + # 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]=Int64[] + end + for i in 1:n + append!(work[(i-1)%nw+1],[i]) + end + + return work +end diff --git a/src/optimization.jl b/src/optimization.jl new file mode 100644 index 0000000..bacead7 --- /dev/null +++ b/src/optimization.jl @@ -0,0 +1,94 @@ +# gradient descent: find local minimum of function of one variable from initial guess +# numerically estimate the derivative +@everywhere function gradient_descent( + f::Function, + x0::Float64, + delta::Float64, # shift is delta*df + dx::Float64, # finite difference for numerical derivative evaluation + maxiter::Int64 # interrupt and fail after maxiter steps +) + counter=0 + + # init + x=x0 + + while counter<maxiter + # value at x and around + val=f(x) + valm=f(x-dx) + valp=f(x+dx) + # quit if minimum + if(val<valm && val<valp) + return(x,val) + end + + # derivative + df=(valp-val)/dx + # step + x=x-delta*df + counter+=1 + end + + # fail + return(Inf,Inf) +end + +# Newton algorithm to compute extrema +# numerically estimate the derivatives +@everywhere function newton_maximum( + f::Function, + x0::Float64, + dx::Float64, # finite difference for numerical derivative evaluation + maxiter::Int64, + tolerance::Float64, + maxstep::Float64 # maximal size of step +) + counter=0 + + # init + x=x0 + + while counter<maxiter + # value at x and around + val=f(x) + valm=f(x-dx) + valp=f(x+dx) + + # derivative + dfp=(valp-val)/dx + dfm=(val-valm)/dx + # second derivative + ddf=(dfp-dfm)/dx + + #@printf(stderr,"% .15e % .15e % .15e % .15e\n",x,val,dfp,ddf) + + if abs(dfp/ddf)<tolerance + # check it is a local maximum + if ddf<0 + return(x,val) + else + return(Inf,Inf) + end + end + + # step + step=dfp/abs(ddf) + + # step too large + if abs(step)>maxstep + step=maxstep*sign(step) + end + + x=x+step + + # fail if off to infinity + if x==Inf || x==-Inf + return(x,val) + end + counter+=1 + end + + # fail + return(Inf,Inf) +end + diff --git a/src/potentials.jl b/src/potentials.jl index 46cafc0..b480db0 100644 --- a/src/potentials.jl +++ b/src/potentials.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,10 +13,15 @@ ## limitations under the License. # exponential potential in 3 dimensions -@everywhere function v_exp(k,a) +@everywhere function v_exp( + k::Float64, + a::Float64 +) return 8*pi/(1+k^2)^2*a end -@everywhere function a0_exp(a) +@everywhere function a0_exp( + a::Float64 +) if a>0. return log(a)+2*MathConstants.eulergamma+2*besselk(0,2*sqrt(a))/besseli(0,2*sqrt(a)) elseif a<0. @@ -27,15 +32,24 @@ end end # exp(-x)-a*exp(-b*x) in 3 dimensions -@everywhere function v_expcry(k,a,b) +@everywhere function v_expcry( + k::Float64, + a::Float64, + b::Float64 +) return 8*pi*((1+k^2)^(-2)-a*b*(b^2+k^2)^(-2)) end -@everywhere function a0_expcry(a,b) +@everywhere function a0_expcry( + a::Float64, + b::Float64 +) return 1.21751642717932720441274114683413710125487579284827462 #ish end # x^2*exp(-|x|) in 3 dimensions -@everywhere function v_npt(k) +@everywhere function v_npt( + k::Float64 +) return 96*pi*(1-k^2)/(1+k^2)^4 end @everywhere function a0_npt() @@ -43,7 +57,9 @@ end end # 1/(1+x^4/4) potential in 3 dimensions -@everywhere function v_alg(k) +@everywhere function v_alg( + k::Float64 +) if(k==0) return 4*pi^2 else @@ -53,32 +69,50 @@ end a0_alg=1. #ish # (1+a x^4)/(1+x^2)^4 potential in 3 dimensions -@everywhere function v_algwell(k) +@everywhere function v_algwell( + k::Float64 +) a=4 return pi^2/24*exp(-k)*(a*(k^2-9*k+15)+k^2+3*k+3) end a0_algwell=1. #ish # potential corresponding to the exact solution c/(1+b^2x^2)^2 -@everywhere function v_exact(k,b,c,e) +@everywhere function v_exact( + k::Float64, + b::Float64, + c::Float64, + e::Float64 +) if k!=0 return 48*pi^2*((18+3*sqrt(c)-(4-3*e/b^2)*c-(1-2*e/b^2)*c^1.5)/(4*(3+sqrt(c))^2*sqrt(c))*exp(-sqrt(1-sqrt(c))*k/b)+(-18+3*sqrt(c)+(4-3*e/b^2)*c-(1-2*e/b^2)*c^1.5)/(4*(3-sqrt(c))^2*sqrt(c))*exp(-sqrt(1+sqrt(c))*k/b)+(1-k/b)/2*exp(-k/b)-c*e/b^2*(3*(9-c)*k/b+8*c)/(8*(9-c)^2)*exp(-2*k/b))/k else return 48*pi^2*(-sqrt(1-sqrt(c))/b*(18+3*sqrt(c)-(4-3*e/b^2)*c-(1-2*e/b^2)*c^1.5)/(4*(3+sqrt(c))^2*sqrt(c))-sqrt(1+sqrt(c))/b*(-18+3*sqrt(c)+(4-3*e/b^2)*c-(1-2*e/b^2)*c^1.5)/(4*(3-sqrt(c))^2*sqrt(c))-1/b-c*e/b^2*(27-19*c)/(8*(9-c)^2)) end end -@everywhere function a0_exact(b,c,e) +@everywhere function a0_exact( + b::Float64, + c::Float64, + e::Float64 +) return 1. #ish end # tent potential (convolution of soft sphere with itself): a*pi/12*(2*|x|/b-2)^2*(2*|x|/b+4) for |x|<b -@everywhere function v_tent(k,a,b) +@everywhere function v_tent( + k::Float64, + a::Float64, + b::Float64 +) if k!=0 return (b/2)^3*a*(4*pi*(sin(k*b/2)-k*b/2*cos(k*b/2))/(k*b/2)^3)^2 else return (b/2)^3*a*(4*pi/3)^2 end end -@everywhere function a0_tent(a,b) +@everywhere function a0_tent( + a::Float64, + b::Float64 +) return b #ish end diff --git a/src/print.jl b/src/print.jl index bef1c4d..3587728 100644 --- a/src/print.jl +++ b/src/print.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. diff --git a/src/simpleq-Kv.jl b/src/simpleq-Kv.jl index 8789656..5a6579c 100644 --- a/src/simpleq-Kv.jl +++ b/src/simpleq-Kv.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,12 +13,27 @@ ## limitations under the License. # Compute Kv=(-\Delta+v+4e(1-\rho u*))^{-1}v -function anyeq_Kv(minlrho,nlrho,taus,P,N,J,rho,a0,v,maxiter,tolerance,xmin,xmax,nx) +function anyeq_Kv( + minlrho::Float64, + nlrho::Int64, + taus::Array{Float64,1}, + P::Int64, + N::Int64, + J::Int64, + rho::Float64, + a0::Float64, + v::Function, + maxiter::Int64, + tolerance::Float64, + xmin::Float64, + xmax::Float64, + nx::Int64 +) # init vectors (weights,T,k,V,V0,A,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v) # compute initial guess from medeq - rhos=Array{Float64}(undef,nlrho) + rhos=Array{Float64,1}(undef,nlrho) for j in 0:nlrho-1 rhos[j+1]=(nlrho==1 ? rho : 10^(minlrho+(log10(rho)-minlrho)/(nlrho-1)*j)) end @@ -44,7 +59,17 @@ function anyeq_Kv(minlrho,nlrho,taus,P,N,J,rho,a0,v,maxiter,tolerance,xmin,xmax, end # Compute the condensate fraction for simpleq using Kv -function simpleq_Kv_condensate_fraction(rhos,taus,P,N,J,a0,v,maxiter,tolerance) +function simpleq_Kv_condensate_fraction( + rhos::Array{Float64,1}, + taus::Array{Float64,1}, + P::Int64, + N::Int64, + J::Int64, + a0::Float64, + v::Function, + maxiter::Int64, + tolerance::Float64 +) # init vectors (weights,T,k,V,V0,A,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v) @@ -67,19 +92,35 @@ function simpleq_Kv_condensate_fraction(rhos,taus,P,N,J,a0,v,maxiter,tolerance) end # Compute the two-point correlation function for simpleq using Kv -function simpleq_Kv_2pt(minlrho,nlrho,taus,P,N,J,rho,a0,v,maxiter,tolerance,xmin,xmax,nx) +function simpleq_Kv_2pt( + minlrho::Float64, + nlrho::Int64, + taus::Array{Float64,1}, + P::Int64, + N::Int64, + J::Int64, + rho::Float64, + a0::Float64, + v::Function, + maxiter::Int64, + tolerance::Float64, + xmin::Float64, + xmax::Float64, + nx::Int64 +) # init vectors (weights,T,k,V,V0,A,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v) # compute initial guess from medeq - rhos=Array{Float64}(undef,nlrho) + rhos=Array{Float64,1}(undef,nlrho) for j in 0:nlrho-1 rhos[j+1]=(nlrho==1 ? rho : 10^(minlrho+(log10(rho)-minlrho)/(nlrho-1)*j)) end u0s=anyeq_init_medeq(rhos,N,J,k,a0,v,maxiter,tolerance) u0=u0s[nlrho] - (u,E,error)=anyeq_hatu(u0,P,N,J,rho,a0,weights,k,taus,V,V0,A,nothing,Upsilon,Upsilon0,v,maxiter,tolerance,Anyeq_approx(0.,0.,1.,0.,0.,0.,0.,0.,0.,0.,0.)) + Abar=Array{Float64,5}(undef,0,0,0,0,0) + (u,E,error)=anyeq_hatu(u0,P,N,J,rho,a0,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,v,maxiter,tolerance,Anyeq_approx(0.,0.,1.,0.,0.,0.,0.,0.,0.,0.,0.)) # Kv in Fourier space Kvk=simpleq_Kv_Kvk(u,V,E,rho,Upsilon,k,taus,weights,N,J) @@ -103,7 +144,18 @@ function simpleq_Kv_2pt(minlrho,nlrho,taus,P,N,J,rho,a0,v,maxiter,tolerance,xmin end # Kv -function simpleq_Kv_Kvk(u,V,E,rho,Upsilon,k,taus,weights,N,J) +function simpleq_Kv_Kvk( + u::Array{Float64,1}, + V::Array{Float64,1}, + E::Float64, + rho::Float64, + Upsilon::Array{Array{Float64,1},1}, + k::Array{Float64,1}, + taus::Array{Float64,1}, + weights::Tuple{Array{Float64,1},Array{Float64,1}}, + N::Int64, + J::Int64 +) # (-Delta+v+4e(1-\rho u*)) in Fourier space M=Array{Float64,2}(undef,N*J,N*J) for zetapp in 0:J-1 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 diff --git a/src/simpleq-iteration.jl b/src/simpleq-iteration.jl index 98977b8..4c4de07 100644 --- a/src/simpleq-iteration.jl +++ b/src/simpleq-iteration.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,7 +13,12 @@ ## limitations under the License. # compute rho(e) using the iteration -function simpleq_iteration_rho_e(es,order,v,maxiter) +function simpleq_iteration_rho_e( + es::Array{Float64,1}, + order::Int64, + v::Function, + maxiter::Int64 +) for j in 1:length(es) (u,rho)=simpleq_iteration_hatun(es[j],order,v,maxiter) @printf("% .15e % .15e\n",es[j],real(rho[maxiter+1])) @@ -21,7 +26,15 @@ function simpleq_iteration_rho_e(es,order,v,maxiter) end # compute u(x) using the iteration and print at every step -function simpleq_iteration_ux(order,e,v,maxiter,xmin,xmax,nx) +function simpleq_iteration_ux( + order::Int64, + e::Float64, + v::Function, + maxiter::Int64, + xmin::Float64, + xmax::Float64, + nx::Int64 +) (u,rho)=simpleq_iteration_hatun(e,order,v,maxiter) weights=gausslegendre(order) @@ -37,7 +50,12 @@ end # \hat u_n -function simpleq_iteration_hatun(e,order,v,maxiter) +function simpleq_iteration_hatun( + e::Float64, + order::Int64, + v::Function, + maxiter::Int64 +) # gauss legendre weights weights=gausslegendre(order) @@ -46,7 +64,7 @@ function simpleq_iteration_hatun(e,order,v,maxiter) (Eta,Eta0)=easyeq_init_H(weights,v) # init u and rho - u=Array{Array{Float64}}(undef,maxiter+1) + u=Array{Array{Float64,1},1}(undef,maxiter+1) u[1]=zeros(Float64,order) rho=zeros(Float64,maxiter+1) @@ -60,7 +78,11 @@ function simpleq_iteration_hatun(e,order,v,maxiter) end # A -function simpleq_iteration_A(e,weights,Eta) +function simpleq_iteration_A( + e::Float64, + weights::Tuple{Array{Float64,1},Array{Float64,1}}, + Eta::Array{Float64,1} +) N=length(weights[1]) out=zeros(Float64,N,N) for i in 1:N @@ -75,7 +97,12 @@ function simpleq_iteration_A(e,weights,Eta) end # b -function simpleq_iteration_b(u,e,rho,V) +function simpleq_iteration_b( + u::Array{Float64,1}, + e::Float64, + rho::Float64, + V::Array{Float64,1} +) out=zeros(Float64,length(V)) for i in 1:length(V) out[i]=V[i]+2*e*rho*u[i]^2 @@ -84,7 +111,13 @@ function simpleq_iteration_b(u,e,rho,V) end # rho_n -function simpleq_iteration_rhon(u,e,weights,V0,Eta0) +function simpleq_iteration_rhon( + u::Array{Float64,1}, + e::Float64, + weights::Tuple{Array{Float64,1},Array{Float64,1}}, + V0::Float64, + Eta0::Float64 +) S=V0 for i in 1:length(weights[1]) y=(weights[1][i]+1)/2 diff --git a/src/tools.jl b/src/tools.jl index 0d3dc7f..5d1678d 100644 --- a/src/tools.jl +++ b/src/tools.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,7 +13,9 @@ ## limitations under the License. # \Phi(x)=2*(1-sqrt(1-x))/x -@everywhere function Phi(x) +@everywhere function Phi( + x::Float64 +) if abs(x)>1e-5 return 2*(1-sqrt(abs(1-x)))/x else @@ -21,7 +23,9 @@ end end # \partial\Phi -@everywhere function dPhi(x) +@everywhere function dPhi( + x::Float64 +) #if abs(x-1)<1e-5 # @printf(stderr,"warning: dPhi is singular at 1, and evaluating it at (% .8e+i% .8e)\n",real(x),imag(x)) #end @@ -33,17 +37,25 @@ end end # apply Phi to every element of a vector -@everywhere function dotPhi(v) +@everywhere function dotPhi( + v::Array{Float64,1} +) out=zeros(Float64,length(v)) for i in 1:length(v) out[i]=Phi(v[i]) end return out end -@everywhere function dotdPhi(v) +@everywhere function dotdPhi( + v::Array{Float64,1} +) out=zeros(Float64,length(v)) for i in 1:length(v) out[i]=dPhi(v[i]) end return out end + +@everywhere function sinc(x::Float64) + return (x == 0 ? 1 : sin(x)/x) +end |