1 files changed, 119 insertions, 0 deletions
diff --git a/src/simpleq-Kv.jl b/src/simpleq-Kv.jl
new file mode 100644
index 0000000..8789656
--- /dev/null
+++ b/src/simpleq-Kv.jl
@@ -0,0 +1,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