software/simplesolv
1
Fork 0
Code Issues Pull Requests Releases Activity
simplesolv/src/easyeq.jl
Ian Jauslin c1b477a1b2 Update to v0.4 
feature: compute the 2-point correlation function in easyeq.

  feature: compute the Fourier transform of the 2-point correlation function
           in anyeq and easyeq.

  feature: compute the local maximum of the 2-point correlation function and
           its Fourier transform.

  feature: compute the compressibility for anyeq.

  feature: allow for linear spacing of rho's.

  feature: print the scattering length.

  change: ux and uk now return real numbers.

  fix: error in the computation of the momentum distribution: wrong
       definition of delta functions.

  fix: various minor bugs.

  optimization: assign explicit types to variables.
2023-02-26 18:36:05 -05:00

1080 lines
26 KiB
Julia
Raw Blame History

## 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.
## You may obtain a copy of the License at
##
## http://www.apache.org/licenses/LICENSE-2.0
##
## Unless required by applicable law or agreed to in writing, software
## distributed under the License is distributed on an "AS IS" BASIS,
## WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
## See the License for the specific language governing permissions and
## limitations under the License.
# interpolation
@everywhere mutable struct Easyeq_approx
bK::Float64
bL::Float64
end
# compute energy
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 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::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)
# 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)
@printf("% .15e % .15e % .15e\n",rhos[j],real(es[j]),errs[j])
end
end
# compute u(k)
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 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])
@printf("% .15e % .15e\n",k,real(u[i]))
end
end
# compute u(x)
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 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
@printf("% .15e % .15e\n",x,real(easyeq_u_x(x,u,weights)))
end
end
# compute 2u(x)-rho u*u(x)
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 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
@printf("% .15e % .15e\n",x,real(easyeq_u_x(x,2*u-rho*u.*u,weights)))
end
end
# condensate fraction
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 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)
# print energy
@printf("% .15e % .15e\n",eta,err)
end
# condensate fraction as a function of rho
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)
# 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 eta
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)
# 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)
(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
k=(1-weights[1][j])/(1+weights[1][j])
u[j]=4*pi*a0/k^2
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::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)
# init u
u=rho*u0
# iterate
err=Inf
for i in 1:maxiter-1
(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)
u=new
break
end
u=new
end
return (u/rho,easyeq_en(u,V0,Eta0,rho,weights)*rho/2,err)
end
# \Eta
@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::Tuple{Array{Float64,1},Array{Float64,1}},
v::Function
)
order=length(weights[1])
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)
end
return(V,V0)
end
# initialize Eta
@everywhere function easyeq_init_H(
weights::Tuple{Array{Float64,1},Array{Float64,1}},
v::Function
)
order=length(weights[1])
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,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)
end
return(Eta,Eta0)
end
# Xi(u)
@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
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)
for i in 1:order
# k_i
k=(1-weights[1][i])/(1+weights[1][i])
# S_i
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[i]=0.
if approx.bK!=0.
A[i]+=approx.bK*S[i]
end
if approx.bK!=1.
A[i]+=(1-approx.bK)*E
end
# T_i
if approx.bK==1.
T[i]=1.
else
T[i]=S[i]/A[i]
end
# B_i
if approx.bK==approx.bL
B[i]=1.
else
B[i]=(approx.bL*S[i]+(1-approx.bL*E))/(approx.bK*S[i]+(1-approx.bK*E))
end
# X_i
X[i]=k^2/(2*A[i]*rho)
end
return (E,S,A,T,B,X)
end
# derivative of Xi
@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)
for i in 1:order
# k_i
k=(1-weights[1][i])/(1+weights[1][i])
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.)
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
# 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
# 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
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
# 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)
for i in 1:order
# k_i
k=(1-weights[1][i])/(1+weights[1][i])
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
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
# 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
# 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
# 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
@everywhere function barEta(
q::Float64,
y::Float64,
t::Float64
)
return (t>=abs(y-q) && t<=y+q ? pi/y/q : 0.)
end
Reference in New Issue View Git Blame Copy Permalink