diff options
Diffstat (limited to 'src/easyeq.jl')
| -rw-r--r-- | src/easyeq.jl | 1040 |
1 files changed, 854 insertions, 186 deletions
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 |