1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
|
## Copyright 2021 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.
# 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)
# 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)
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,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)
# switch to real space
for i in 1:nx
x=xmin+(xmax-xmin)*i/nx
Kv=0.
for zeta in 0:J-1
for j in 1:N
Kv+=(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]*Kvk[zeta*N+j]*sin(k[zeta*N+j]*x)
end
end
@printf("% .15e % .15e\n",x,Kv)
end
end
# Compute the condensate fraction for simpleq using Kv
function simpleq_Kv_condensate_fraction(rhos,taus,P,N,J,a0,v,maxiter,tolerance)
# init vectors
(weights,T,k,V,V0,A,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v)
# compute initial guess from medeq
u0s=anyeq_init_medeq(rhos,N,J,k,a0,v,0.,maxiter,tolerance)
for j in 1:length(rhos)
(u,E,error)=anyeq_hatu(u0s[j],0.,P,N,J,rhos[j],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)
# denominator: (1-rho*\int (Kv)(2u-rho u*u))
denom=1-rhos[j]*integrate_f_chebyshev(s->1.,Kvk.*(2*u-rhos[j]*u.^2),k,taus,weights,N,J)
eta=rhos[j]*integrate_f_chebyshev(s->1.,Kvk.*u,k,taus,weights,N,J)/denom
@printf("% .15e % .15e\n",rhos[j],eta)
end
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)
# 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)
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.))
# Kv in Fourier space
Kvk=simpleq_Kv_Kvk(u,V,E,rho,Upsilon,k,taus,weights,N,J)
# denominator: (1-rho*\int (Kv)(2u-rho u*u))
denom=1-rho*integrate_f_chebyshev(s->1.,Kvk.*(2*u-rho*u.^2),k,taus,weights,N,J)
# Kv, u, u*Kv, u*u*Kv in real space
for i in 1:nx
x=xmin+(xmax-xmin)*i/nx
ux=inverse_fourier_chebyshev(u,x,k,taus,weights,N,J)
Kv=inverse_fourier_chebyshev(Kvk,x,k,taus,weights,N,J)
uKv=inverse_fourier_chebyshev(u.*Kvk,x,k,taus,weights,N,J)
uuKv=inverse_fourier_chebyshev(u.*u.*Kvk,x,k,taus,weights,N,J)
# correlation in real space
Cx=rho^2*((1-ux)-((1-ux)*Kv-2*rho*uKv+rho^2*uuKv)/denom)
@printf("% .15e % .15e\n",x,Cx)
end
end
# Kv
function simpleq_Kv_Kvk(u,V,E,rho,Upsilon,k,taus,weights,N,J)
# (-Delta+v+4e(1-\rho u*)) in Fourier space
M=Array{Float64,2}(undef,N*J,N*J)
for zetapp in 0:J-1
for n in 1:N
M[:,zetapp*N+n]=conv_one_v_chebyshev(n,zetapp,Upsilon,k,taus,weights,N,J)
end
end
for i in 1:J*N
M[i,i]+=k[i]^2+4*E*(1-rho*u[i])
end
return inv(M)*V
end
|
