## 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