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

# interpolation
@everywhere mutable struct Anyeq_approx
  aK::Float64
  bK::Float64
  gK::Float64
  aL1::Float64
  bL1::Float64
  aL2::Float64
  bL2::Float64
  gL2::Float64
  aL3::Float64
  bL3::Float64
  gL3::Float64
end

# compute energy for a given rho
# use minlrho, nlrho to incrementally compute the solution to medeq (to initialize the Newton algorithm)
function anyeq_energy(rho,minlrho,nlrho,taus,P,N,J,a0,v,maxiter,tolerance,approx,savefile)
  # initialize vectors
  (weights,T,k,V,V0,A,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v)
  # init Abar
  if savefile!=""
    Abar=anyeq_read_Abar(savefile,P,N,J)
  else
    Abar=anyeq_Abar_multithread(k,weights,taus,T,P,N,J,approx)
  end

  # 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,approx)

  @printf("% .15e % .15e\n",E,error)
end

# compute energy as a function of rho
function anyeq_energy_rho(rhos,taus,P,N,J,a0,v,maxiter,tolerance,approx,savefile)
  # initialize vectors
  (weights,T,k,V,V0,A,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v)

  # init Abar
  if savefile!=""
    Abar=anyeq_read_Abar(savefile,P,N,J)
  else
    Abar=anyeq_Abar_multithread(k,weights,taus,T,P,N,J,approx)
  end

  # compute initial guess from medeq
  u0s=anyeq_init_medeq(rhos,N,J,k,a0,v,maxiter,tolerance)

  # save result from each task
  es=Array{Float64,1}(undef,length(rhos))
  err=Array{Float64,1}(undef,length(rhos))

  ## spawn workers
  # number of workers
  nw=nworkers()
  # split jobs among workers
  work=Array{Array{Int64,1},1}(undef,nw)
  # init empty arrays
  for p in 1:nw
    work[p]=zeros(0)
  end
  for j in 1:length(rhos)
    append!(work[(j-1)%nw+1],j)
  end

  count=0
  # for each worker
  @sync for p in 1:nw
    # for each task
    @async for j in work[p]
      count=count+1
      if count>=nw
        progress(count,length(rhos),10000)
      end
      # run the task
      (u,es[j],err[j])=remotecall_fetch(anyeq_hatu,workers()[p],u0s[j],P,N,J,rhos[j],a0,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,v,maxiter,tolerance,approx)
    end
  end

  for j in 1:length(rhos)
    @printf("% .15e % .15e % .15e\n",rhos[j],es[j],err[j])
  end
end

# compute energy as a function of rho
# initialize with previous rho
function anyeq_energy_rho_init_prevrho(rhos,taus,P,N,J,a0,v,maxiter,tolerance,approx,savefile)
  # initialize vectors
  (weights,T,k,V,V0,A,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v)

  # init Abar
  if savefile!=""
    Abar=anyeq_read_Abar(savefile,P,N,J)
  else
    Abar=anyeq_Abar_multithread(k,weights,taus,T,P,N,J,approx)
  end

  # compute initial guess from medeq
  u0s=anyeq_init_medeq([rhos[1]],N,J,k,a0,v,maxiter,tolerance)
  u=u0s[1]

  for j in 1:length(rhos)
    progress(j,length(rhos),10000)
    # run the task
    # init Newton from previous rho
    (u,E,error)=anyeq_hatu(u,P,N,J,rhos[j],a0,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,v,maxiter,tolerance,approx)
    @printf("% .15e % .15e % .15e\n",rhos[j],E,error)
    # abort when the error gets too big
    if error>tolerance
      break
    end
  end
end
# compute energy as a function of rho
# initialize with next rho
function anyeq_energy_rho_init_nextrho(rhos,taus,P,N,J,a0,v,maxiter,tolerance,approx,savefile)
  # initialize vectors
  (weights,T,k,V,V0,A,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v)

  # init Abar
  if savefile!=""
    Abar=anyeq_read_Abar(savefile,P,N,J)
  else
    Abar=anyeq_Abar_multithread(k,weights,taus,T,P,N,J,approx)
  end

  # compute initial guess from medeq
  u0s=anyeq_init_medeq([rhos[length(rhos)]],N,J,k,a0,v,maxiter,tolerance)
  u=u0s[1]

  for j in 1:length(rhos)
    progress(j,length(rhos),10000)
    # run the task
    # init Newton from previous rho
    (u,E,error)=anyeq_hatu(u,P,N,J,rhos[length(rhos)+1-j],a0,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,v,maxiter,tolerance,approx)
    @printf("% .15e % .15e % .15e\n",rhos[length(rhos)+1-j],real(E),error)
    # abort when the error gets too big
    if error>tolerance
      break
    end
  end
end

# compute u(k)
# use minlrho, nlrho to incrementally compute the solution to medeq (to initialize the Newton algorithm)
function anyeq_uk(minlrho,nlrho,taus,P,N,J,rho,a0,v,maxiter,tolerance,approx,savefile)
  # init vectors
  (weights,T,k,V,V0,A,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v)
  # init Abar
  if savefile!=""
    Abar=anyeq_read_Abar(savefile,P,N,J)
  else
    Abar=anyeq_Abar_multithread(k,weights,taus,T,P,N,J,approx)
  end
  # 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,approx)

  for zeta in 0:J-1
    for j in 1:N
      # order k's in increasing order
      @printf("% .15e % .15e\n",k[(J-1-zeta)*N+j],u[(J-1-zeta)*N+j])
    end
  end
end

# compute u(x)
# use minlrho, nlrho to incrementally compute the solution to medeq (to initialize the Newton algorithm)
function anyeq_ux(minlrho,nlrho,taus,P,N,J,rho,a0,v,maxiter,tolerance,xmin,xmax,nx,approx,savefile)
  # init vectors
  (weights,T,k,V,V0,A,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v)
  # init Abar
  if savefile!=""
    Abar=anyeq_read_Abar(savefile,P,N,J)
  else
    Abar=anyeq_Abar_multithread(k,weights,taus,T,P,N,J,approx)
  end

  # 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,approx)

  for i in 1:nx
    x=xmin+(xmax-xmin)*i/nx
    ux=0.
    for zeta in 0:J-1
      for j in 1:N
	ux+=(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]*u[zeta*N+j]*sin(k[zeta*N+j]*x)
      end
    end
    @printf("% .15e % .15e % .15e\n",x,real(ux),imag(ux))
  end
end

# compute condensate fraction for a given rho
# use minlrho, nlrho to incrementally compute the solution to medeq (to initialize the Newton algorithm)
function anyeq_condensate_fraction(rho,minlrho,nlrho,taus,P,N,J,a0,v,maxiter,tolerance,approx,savefile)
  # initialize vectors
  (weights,T,k,V,V0,A,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v)
  # init Abar
  if savefile!=""
    Abar=anyeq_read_Abar(savefile,P,N,J)
  else
    Abar=anyeq_Abar_multithread(k,weights,taus,T,P,N,J,approx)
  end

  # 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,approx)

  # compute eta
  eta=anyeq_eta(u,P,N,J,rho,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,approx)

  @printf("% .15e % .15e\n",eta,error)
end

# condensate fraction as a function of rho
function anyeq_condensate_fraction_rho(rhos,taus,P,N,J,a0,v,maxiter,tolerance,approx,savefile)
  ## spawn workers
  # number of workers
  nw=nworkers()
  # split jobs among workers
  work=Array{Array{Int64,1},1}(undef,nw)
  # init empty arrays
  for p in 1:nw
    work[p]=zeros(0)
  end
  for j in 1:length(rhos)
    append!(work[(j-1)%nw+1],j)
  end

  # initialize vectors
  (weights,T,k,V,V0,A,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v)
  # init Abar
  if savefile!=""
    Abar=anyeq_read_Abar(savefile,P,N,J)
  else
    Abar=anyeq_Abar_multithread(k,weights,taus,T,P,N,J,approx)
  end

  # compute initial guess from medeq
  u0s=anyeq_init_medeq(rhos,N,J,k,a0,v,maxiter,tolerance)

  # compute u
  us=Array{Array{Float64,1}}(undef,length(rhos))
  errs=Array{Float64,1}(undef,length(rhos))
  count=0
  # for each worker
  @sync for p in 1:nw
    # for each task
    @async for j in work[p]
      count=count+1
      if count>=nw
	progress(count,length(rhos),10000)
      end
      # run the task
      (us[j],E,errs[j])=remotecall_fetch(anyeq_hatu,workers()[p],u0s[j],P,N,J,rhos[j],a0,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,v,maxiter,tolerance,approx)
    end
  end

  # compute eta
  etas=Array{Float64}(undef,length(rhos))
  count=0
  # for each worker
  @sync for p in 1:nw
    # for each task
    @async for j in work[p]
      count=count+1
      if count>=nw
	progress(count,length(rhos),10000)
      end
      # run the task
      etas[j]=anyeq_eta(us[j],P,N,J,rhos[j],weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,approx)
    end
  end

  for j in 1:length(rhos)
    @printf("% .15e % .15e % .15e\n",rhos[j],etas[j],errs[j])
  end
end

# compute the momentum distribution for a given rho
# use minlrho, nlrho to incrementally compute the solution to medeq (to initialize the Newton algorithm)
function anyeq_momentum_distribution(rho,minlrho,nlrho,taus,P,N,J,a0,v,maxiter,tolerance,approx,savefile)
  # initialize vectors
  (weights,T,k,V,V0,A,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v)
  # init Abar
  if savefile!=""
    Abar=anyeq_read_Abar(savefile,P,N,J)
  else
    Abar=anyeq_Abar_multithread(k,weights,taus,T,P,N,J,approx)
  end

  # 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,approx)

  # compute M
  M=anyeq_momentum(u,P,N,J,rho,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,approx)

  for zeta in 0:J-1
    for j in 1:N
      # order k's in increasing order
      @printf("% .15e % .15e\n",k[(J-1-zeta)*N+j],M[(J-1-zeta)*N+j])
    end
  end
end

# 2 point correlation function
function anyeq_2pt_correlation(minlrho,nlrho,taus,P,N,J,windowL,rho,a0,v,maxiter,tolerance,xmin,xmax,nx,approx,savefile)
  # init vectors
  (weights,T,k,V,V0,A,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v)
  # init Abar
  if savefile!=""
    Abar=anyeq_read_Abar(savefile,P,N,J)
  else
    Abar=anyeq_Abar_multithread(k,weights,taus,T,P,N,J,approx)
  end

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

  # compute u and some useful integrals
  (u,E,error)=anyeq_hatu(u0,P,N,J,rho,a0,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,v,maxiter,tolerance,approx)
  (S,E,II,JJ,X,Y,sL1,sK,G)=anyeq_SEIJGXY(rho*u,rho,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,weights,P,N,J,approx)

  for i in 1:nx
    x=xmin+(xmax-xmin)*i/nx
    C2=anyeq_2pt(x,u,P,N,J,windowL,rho,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,approx,S,E,II,JJ,X,Y,sL1,sK,G)
    @printf("% .15e % .15e\n",x,C2)
  end
end

# compute Abar
function anyeq_Abar_multithread(k,weights,taus,T,P,N,J,approx)
  if approx.bL3==0.
    return []
  end

  out=Array{Array{Float64,5}}(undef,J*N)

  ## spawn workers
  # number of workers
  nw=nworkers()
  # split jobs among workers
  work=Array{Array{Int64,1},1}(undef,nw)
  # init empty arrays
  for p in 1:nw
    work[p]=zeros(0)
  end
  for j in 1:N*J
    append!(work[(j-1)%nw+1],j)
  end

  count=0
  # for each worker
  @sync for p in 1:nw
    # for each task
    @async for j in work[p]
      count=count+1
      if count>=nw
        progress(count,N*J,10000)
      end
      # run the task
      out[j]=remotecall_fetch(barAmat,workers()[p],k[j],weights,taus,T,P,N,J,2,2,2,2,2)
    end
  end

  return out
end

# initialize computation
@everywhere function anyeq_init(taus,P,N,J,v)
  # Gauss-Legendre weights
  weights=gausslegendre(N)

  # initialize vectors V,k
  V=Array{Float64}(undef,J*N)
  k=Array{Float64}(undef,J*N)
  for zeta in 0:J-1
    for j in 1:N
      xj=weights[1][j]
      # set kj
      k[zeta*N+j]=(2+(taus[zeta+2]-taus[zeta+1])*sin(pi*xj/2)-(taus[zeta+2]+taus[zeta+1]))/(2-(taus[zeta+2]-taus[zeta+1])*sin(pi*xj/2)+taus[zeta+2]+taus[zeta+1])
      # set v
      V[zeta*N+j]=v(k[zeta*N+j])
    end
  end
  # potential at 0
  V0=v(0)

  # initialize matrix A
  T=chebyshev_polynomials(P)
  A=Amat(k,weights,taus,T,P,N,J,2,2)

  # compute Upsilon
  # Upsilonmat does not use splines, so increase precision
  weights_plus=gausslegendre(N*J)
  Upsilon=Upsilonmat(k,v,weights_plus)
  Upsilon0=Upsilon0mat(k,v,weights_plus)

  return(weights,T,k,V,V0,A,Upsilon,Upsilon0)
end

# compute initial guess from medeq
@everywhere function anyeq_init_medeq(rhos,N,J,k,a0,v,maxiter,tolerance)
  us_medeq=Array{Array{Float64,1}}(undef,length(rhos))
  u0s=Array{Array{Float64,1}}(undef,length(rhos))

  weights_medeq=gausslegendre(N*J)

  (us_medeq[1],E,err)=easyeq_hatu(easyeq_init_u(a0,J*N,weights_medeq),J*N,rhos[1],v,maxiter,tolerance,weights_medeq,Easyeq_approx(1.,1.))
  u0s[1]=easyeq_to_anyeq(us_medeq[1],weights_medeq,k,N,J)
  if err>tolerance
    print(stderr,"warning: computation of initial Ansatz failed for rho=",rhos[1],"\n")
  end

  for j in 2:length(rhos)
    (us_medeq[j],E,err)=easyeq_hatu(us_medeq[j-1],J*N,rhos[j],v,maxiter,tolerance,weights_medeq,Easyeq_approx(1.,1.))
    u0s[j]=easyeq_to_anyeq(us_medeq[j],weights_medeq,k,N,J)

    if err>tolerance
      print(stderr,"warning: computation of initial Ansatz failed for rho=",rhos[j],"\n")
    end
  end

  return u0s
end

# interpolate the solution of medeq to an input for anyeq
@everywhere function easyeq_to_anyeq(u_simple,weights,k,N,J)
  # reorder u_simple, which is evaluated at (1-x_j)/(1+x_j) with x_j\in[-1,1]
  u_s=zeros(Float64,length(u_simple))
  k_s=Array{Float64}(undef,length(u_simple))
  for j in 1:length(u_simple)
    xj=weights[1][j]
    k_s[length(u_simple)-j+1]=(1-xj)/(1+xj)
    u_s[length(u_simple)-j+1]=u_simple[j]
  end

  # initialize U
  u=zeros(Float64,J*N)
  for zeta in 0:J-1
    for j in 1:N
      u[zeta*N+j]=linear_interpolation(k[zeta*N+j],k_s,u_s)
    end
  end

  return u
end


# compute u using chebyshev expansions
@everywhere function anyeq_hatu(u0,P,N,J,rho,a0,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,v,maxiter,tolerance,approx)
  # init
  # rescale by rho (that's how u is defined)
  u=rho*u0

  # quantify relative error
  error=-1.

  # run Newton algorithm
  for i in 1:maxiter-1
    (S,E,II,JJ,X,Y,sL1,sK,G)=anyeq_SEIJGXY(u,rho,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,weights,P,N,J,approx)
    new=u-inv(anyeq_DXi(u,rho,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,weights,P,N,J,approx,S,E,II,JJ,X,Y,sL1,sK))*anyeq_Xi(u,X,Y)

    error=norm(new-u)/norm(u)
    if(error<tolerance)
      u=new
      break
    end

    u=new
  end

  energy=rho/2*V0-avg_v_chebyshev(u,Upsilon0,k,taus,weights,N,J)/2
  return(u/rho,energy,error)
end


# save Abar
function anyeq_save_Abar(taus,P,N,J,v,approx)
  # initialize vectors
  (weights,T,k,V,V0,A,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v)

  # init Abar
  Abar=anyeq_Abar_multithread(k,weights,taus,T,P,N,J,approx)

  # print params
  @printf("## P=%d N=%d J=%d\n",P,N,J)

  for i in 1:N*J
    for j1 in 1:(P+1)*J
      for j2 in 1:(P+1)*J
	for j3 in 1:(P+1)*J
	  for j4 in 1:(P+1)*J
	    for j5 in 1:(P+1)*J
	      @printf("% .15e\n",Abar[i][j1,j2,j3,j4,j5])
	    end
	  end
	end
      end
    end
  end
end

# read Abar
function anyeq_read_Abar(savefile,P,N,J)
   # open file
  ff=open(savefile,"r")
  # read all lines
  lines=readlines(ff)
  close(ff)

  # init
  Abar=Array{Array{Float64,5}}(undef,N*J)
  for i in 1:N*J
    Abar[i]=Array{Float64,5}(undef,(P+1)*J,(P+1)*J,(P+1)*J,(P+1)*J,(P+1)*J)
  end

  i=1
  j1=1
  j2=1
  j3=1
  j4=1
  j5=0
  for l in 1:length(lines)
    # drop comments
    if lines[l]!="" && lines[l][1]!='#'
      # increment counters
      if j5<(P+1)*J
	j5+=1
      else
	j5=1
	if j4<(P+1)*J
	  j4+=1
	else
	  j4=1
	  if j3<(P+1)*J
	    j3+=1
	  else
	    j3=1
	    if j2<(P+1)*J
	      j2+=1
	    else
	      j2=1
	      if j1<(P+1)*J
		j1+=1
	      else
		j1=1
		if i<N*J
		  i+=1
		else
		  print(stderr,"error: too many lines in savefile\n")
		  exit()
		end
	      end
	    end
	  end
	end
      end

      Abar[i][j1,j2,j3,j4,j5]=parse(Float64,lines[l])

    end
  end

  return Abar

end



# Xi
# takes the vector of kj's and xn's as input
@everywhere function anyeq_Xi(U,X,Y)
  return U-(Y.+1)./(2*(X.+1)).*dotPhi((Y.+1)./((X.+1).^2))
end

# DXi
# takes the vector of kj's as input
@everywhere function anyeq_DXi(U,rho,k,taus,v,v0,A,Abar,Upsilon,Upsilon0,weights,P,N,J,approx,S,E,II,JJ,X,Y,sL1,sK)
  out=Array{Float64,2}(undef,N*J,N*J)
  for zetapp in 0:J-1
    for n in 1:N
      one=zeros(Int64,J*N)
      one[zetapp*N+n]=1

      # Chebyshev expansion of U
      FU=chebyshev(U,taus,weights,P,N,J,2)

      dS=-conv_one_v_chebyshev(n,zetapp,Upsilon,k,taus,weights,N,J)/rho

      dE=-avg_one_v_chebyshev(n,zetapp,Upsilon0,k,taus,weights,N)/rho

      UU=conv_chebyshev(FU,FU,A)

      dII=zeros(Float64,N*J)
      if approx.gL2!=0.
	if approx.bL2!=0.
	  dII+=approx.gL2*approx.bL2*(conv_one_chebyshev(n,zetapp,chebyshev(U.*S,taus,weights,P,N,J,2),A,taus,weights,P,N,J,2)/rho+conv_chebyshev(FU,chebyshev(one.*S+U.*dS,taus,weights,P,N,J,2),A)/rho)
	end
	if approx.bL2!=1.
	  dII+=approx.gL2*(1-approx.bL2)*(dE/rho*UU+2*E/rho*conv_one_chebyshev(n,zetapp,FU,A,taus,weights,P,N,J,2))
	end
      end

      dJJ=zeros(Float64,J*N)
      if approx.gL3!=0.
	if approx.bL3!=0.
	  FS=chebyshev(S,taus,weights,P,N,J,2)
	  dFU=chebyshev(one,taus,weights,P,N,J,2)
	  dFS=chebyshev(dS,taus,weights,P,N,J,2)
	  dJJ+=approx.gL3*approx.bL3*(4*double_conv_S_chebyshev(FU,FU,FU,dFU,FS,Abar)+double_conv_S_chebyshev(FU,FU,FU,FU,dFS,Abar))
	end
	if approx.bL3!=1.
	  dJJ+=approx.gL3*(1-approx.bL3)*(dE*(UU/rho).^2+4*E*conv_one_chebyshev(n,zetapp,FU,A,taus,weights,P,N,J,2).*UU/rho^2)
	end
      end

      dG=zeros(Float64,N*J)
      if approx.aK!=0. && approx.gK!=0.
	if approx.bK!=0.
	  dG+=approx.aK*approx.gK*approx.bK*(conv_one_chebyshev(n,zetapp,chebyshev(2*S.*U,taus,weights,P,N,J,2),A,taus,weights,P,N,J,2)/rho+conv_chebyshev(FU,chebyshev(2*S.*one+2*dS.*U,taus,weights,P,N,J,2),A)/rho)
	end
	if approx.bK!=1.
	  dG+=approx.aK*approx.gK*(1-approx.bK)*(2*dE*UU/rho+4*E*conv_one_chebyshev(n,zetapp,FU,A,taus,weights,P,N,J,2)/rho)
	end
      end
      if approx.aL1!=0.
	if approx.bL1!=0.
	  dG-=approx.aL1*approx.bL1*(conv_one_chebyshev(n,zetapp,chebyshev(S.*(U.^2),taus,weights,P,N,J,2),A,taus,weights,P,N,J,2)/rho+conv_chebyshev(FU,chebyshev(2*S.*U.*one+dS.*(U.^2),taus,weights,P,N,J,2),A)/rho)
	end
	if approx.bL1!=1.
	  dG-=approx.aL1*(1-approx.bL1)*(E/rho*conv_one_chebyshev(n,zetapp,chebyshev((U.^2),taus,weights,P,N,J,2),taus,A,weights,P,N,J,2)+conv_chebyshev(FU,chebyshev(2*E*U.*one+dE*(U.^2),taus,weights,P,N,J,2),A)/rho)
	end
      end
      if approx.aL2!=0. && approx.gL2!=0.
	dG+=approx.aL2*(conv_one_chebyshev(n,zetapp,chebyshev(2*II.*U,taus,weights,P,N,J,2),A,taus,weights,P,N,J,2)/rho+conv_chebyshev(FU,chebyshev(2*dII.*U+2*II.*one,taus,weights,P,N,J,2),A)/rho)
      end
      if approx.aL3!=0. && approx.gL3!=0.
	dG-=approx.aL3*(conv_one_chebyshev(n,zetapp,chebyshev(JJ/2,taus,weights,P,N,J,2),A,taus,weights,P,N,J,2)/rho+conv_chebyshev(FU,chebyshev(dJJ/2,taus,weights,P,N,J,2),A)/rho)
      end

      dsK=zeros(Float64,N*J)
      if approx.gK!=0.
	if approx.bK!=0.
	  dsK+=approx.gK*approx.bK*dS
	end
	if approx.bK!=1.
	  dsK+=approx.gK*(1-approx.bK)*dE*ones(N*J)
	end
      else
      end
      dsL1=zeros(Float64,N*J)
      if approx.bL1!=0.
	dsL1+=approx.bL1*dS
      end
      if approx.bL1!=1.
	dsL1+=(1-approx.bL1)*dE*ones(Float64,N*J)
      end

      dX=(dsK-dsL1+dII)./sL1-X./sL1.*dsL1
      dY=(dS-dsL1+dG+dJJ/2)./sL1-Y./sL1.*dsL1

      out[:,zetapp*N+n]=one+((Y.+1).*dX./(X.+1)-dY)./(2*(X.+1)).*dotPhi((Y.+1)./((X.+1).^2))+(Y.+1)./(2*(X.+1).^3).*(2*(Y.+1)./(X.+1).*dX-dY).*dotdPhi((Y.+1)./(X.+1).^2)
    end
  end
  return out
end

# compute S,E,I,J,X and Y
@everywhere function anyeq_SEIJGXY(U,rho,k,taus,v,v0,A,Abar,Upsilon,Upsilon0,weights,P,N,J,approx)
  # Chebyshev expansion of U
  FU=chebyshev(U,taus,weights,P,N,J,2)

  S=v-conv_v_chebyshev(U,Upsilon,k,taus,weights,N,J)/rho
  E=v0-avg_v_chebyshev(U,Upsilon0,k,taus,weights,N,J)/rho

  UU=conv_chebyshev(FU,FU,A)

  II=zeros(Float64,N*J)
  if approx.gL2!=0.
    if approx.bL2!=0.
      II+=approx.gL2*approx.bL2*conv_chebyshev(FU,chebyshev(U.*S,taus,weights,P,N,J,2),A)/rho
    end
    if approx.bL2!=1.
      II+=approx.gL2*(1-approx,bL2)*E/rho*UU
    end
  end

  JJ=zeros(Float64,N*J)
  if approx.gL3!=0.
    if approx.bL3!=0.
      FS=chebyshev(S,taus,weights,P,N,J,2)
      JJ+=approx.gL3*approx.bL3*double_conv_S_chebyshev(FU,FU,FU,FU,FS,Abar)
    end
    if approx.bL3!=1.
      JJ+=approx.gL3*(1-approx.bL3)*E*(UU/rho).^2
    end
  end

  G=zeros(Float64,N*J)
  if approx.aK!=0. && approx.gK!=0.
    if approx.bK!=0.
      G+=approx.aK*approx.gK*approx.bK*conv_chebyshev(FU,chebyshev(2*S.*U,taus,weights,P,N,J,2),A)/rho
    end
    if approx.bK!=1.
      G+=approx.aK*approx.gK*(1-approx.bK)*2*E*UU/rho
    end
  end
  if approx.aL1!=0.
    if approx.bL1!=0.
      G-=approx.aL1*approx.bL1*conv_chebyshev(FU,chebyshev(S.*(U.^2),taus,weights,P,N,J,2),A)/rho
    end
    if approx.bL1!=1.
      G-=approx.aL1*(1-approx.bL1)*E/rho*conv_chebyshev(FU,chebyshev((U.^2),taus,weights,P,N,J,2),A)
    end
  end
  if approx.aL2!=0. && approx.gL2!=0.
    G+=approx.aL2*conv_chebyshev(FU,chebyshev(2*II.*U,taus,weights,P,N,J,2),A)/rho
  end
  if approx.aL3!=0 && approx.gL3!=0.
    G-=approx.aL3*conv_chebyshev(FU,chebyshev(JJ/2,taus,weights,P,N,J,2),A)/rho
  end

  sK=zeros(Float64,N*J)
  if approx.gK!=0.
    if approx.bK!=0.
      sK+=approx.gK*approx.bK*S
    end
    if approx.bK!=1.
      sK+=approx.gK*(1-approx.bK)*E*ones(Float64,N*J)
    end
  end

  sL1=zeros(Float64,N*J)
  if approx.bL1!=0.
    sL1+=approx.bL1*S
  end
  if approx.bL1!=1.
    sL1+=(1-approx.bL1)*E*ones(Float64,N*J)
  end

  X=(k.^2/2+rho*(sK-sL1+II))./sL1/rho
  Y=(S-sL1+G+JJ/2)./sL1

  return(S,E,II,JJ,X,Y,sL1,sK,G)
end

# condensate fraction
@everywhere function anyeq_eta(u,P,N,J,rho,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,approx)
  # compute dXi/dmu
  (S,E,II,JJ,X,Y,sL1,sK,G)=anyeq_SEIJGXY(rho*u,rho,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,weights,P,N,J,approx)
  dXidmu=(Y.+1)./(rho*sL1)./(2*(X.+1).^2).*dotPhi((Y.+1)./((X.+1).^2))+(Y.+1).^2 ./((X.+1).^4)./(rho*sL1).*dotdPhi((Y.+1)./(X.+1).^2)

  # compute eta
  du=-inv(anyeq_DXi(rho*u,rho,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,weights,P,N,J,approx,S,E,II,JJ,X,Y,sL1,sK))*dXidmu
  eta=-avg_v_chebyshev(du,Upsilon0,k,taus,weights,N,J)/2

  return eta
end

# momentum distribution
@everywhere function anyeq_momentum(u,P,N,J,rho,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,approx)
  # compute dXi/dlambda (without delta functions)
  (S,E,II,JJ,X,Y,sL1,sK,G)=anyeq_SEIJGXY(rho*u,rho,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,weights,P,N,J,approx)
  dXidlambda=-(2*pi)^3*2*u./sL1.*(dotPhi((Y.+1)./((X.+1).^2))./(2*(X.+1))+(Y.+1)./(2*(X.+1).^3).*dotdPhi((Y.+1)./(X.+1).^2))

  # approximation for delta function (without Kronecker deltas)
  delta=Array{Float64}(undef,J*N)
  for zeta in 0:J-1
    for n in 1:N
      delta[zeta*N+n]=2/pi^2/((taus[zeta+2]-taus[zeta+1])*weights[2][n]*cos(pi*weights[1][n]/2)*(1+k[zeta*N+n])^2*k[zeta*N+n]^2)
    end
  end
  
  # compute dXidu
  dXidu=inv(anyeq_DXi(rho*u,rho,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,weights,P,N,J,approx,S,E,II,JJ,X,Y,sL1,sK))

  M=Array{Float64}(undef,J*N)
  for i in 1:J*N
    # du/dlambda
    du=dXidu[:,i]*dXidlambda[i]*delta[i]

    # compute M
    M[i]=-avg_v_chebyshev(du,Upsilon0,k,taus,weights,N,J)/2
  end

  return M
end


# correlation function
@everywhere function anyeq_2pt(x,u,P,N,J,windowL,rho,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,approx,S,E,II,JJ,X,Y,sL1,sK,G)
  # initialize dV
  dV=Array{Float64}(undef,J*N)
  for i in 1:J*N
    if x>0
      dV[i]=sin(k[i]*x)/(k[i]*x)*hann(k[i],windowL)
    else
      dV[i]=hann(k[i],windowL)
    end
  end
  dV0=1.

  # compute dUpsilon
  # Upsilonmat does not use splines, so increase precision
  weights_plus=gausslegendre(N*J)
  dUpsilon=Upsilonmat(k,r->sin(r*x)/(r*x)*hann(r,windowL),weights_plus)
  dUpsilon0=Upsilon0mat(k,r->sin(r*x)/(r*x)*hann(r,windowL),weights_plus)

  du=-inv(anyeq_DXi(rho*u,rho,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,weights,P,N,J,approx,S,E,II,JJ,X,Y,sL1,sK))*anyeq_dXidv(x,rho*u,rho,k,taus,dV,dV0,A,Abar,dUpsilon,dUpsilon0,weights,P,N,J,approx,S,E,II,JJ,X,Y,sL1,sK)
  # rescale rho
  du=du/rho

  C2=rho^2*(1-integrate_f_chebyshev(s->1.,u.*dV+V.*du,k,taus,weights,N,J))

  return C2
end

# derivative of Xi with respect to v in the direction sin(kx)/kx
@everywhere function anyeq_dXidv(x,U,rho,k,taus,dv,dv0,A,Abar,dUpsilon,dUpsilon0,weights,P,N,J,approx,S,E,II,JJ,X,Y,sL1,sK)
  # Chebyshev expansion of U
  FU=chebyshev(U,taus,weights,P,N,J,2)

  dS=dv-conv_v_chebyshev(U,dUpsilon,k,taus,weights,N,J)/rho
  dE=dv0-avg_v_chebyshev(U,dUpsilon0,k,taus,weights,N,J)/rho

  UU=conv_chebyshev(FU,FU,A)

  dII=zeros(Float64,N*J)
  if approx.gL2!=0.
    if approx.bL2!=0.
      dII+=approx.gL2*approx.bL2*conv_chebyshev(FU,chebyshev(U.*dS,taus,weights,P,N,J,2),A)/rho
    end
    if approx.bL2!=1.
      dII+=approx.gL2*(1-approx,bL2)*dE/rho*UU
    end
  end

  dJJ=zeros(Float64,J*N)
  if approx.gL3!=0.
    if approx.bL3!=0.
      dFS=chebyshev(dS,taus,weights,P,N,J,2)
      dJJ+=approx.gL3*approx.bL3*double_conv_S_chebyshev(FU,FU,FU,FU,dFS,Abar)
    end
    if approx.bL3!=1.
      dJJ=approx.gL3*(1-approx.bL3)*dE*(UU/rho).^2
    end
  end

  dG=zeros(Float64,N*J)
  if approx.aK!=0. && approx.gK!=0.
    if approx.bK!=0.
      dG+=approx.aK*approx.gK*approx.bK*conv_chebyshev(FU,chebyshev(2*dS.*U,taus,weights,P,N,J,2),A)/rho
    end
    if approx.bK!=1.
      dG+=approx.aK*approx.gK*(1-approx.bK)*2*dE*UU/rho
    end
  end
  if approx.aL1!=0.
    if approx.bL1!=0.
      dG-=approx.aL1*approx.bL1*conv_chebyshev(FU,chebyshev(dS.*(U.^2),taus,weights,P,N,J,2),A)/rho
    end
    if approx.bL1!=1.
      dG-=approx.aL1*(1-approx.bL1)*dE/rho*conv_chebyshev(FU,chebyshev((U.^2),taus,weights,P,N,J,2),A)
    end
  end
  if approx.aL2!=0. && approx.gL2!=0.
    dG+=approx.aL2*conv_chebyshev(FU,chebyshev(2*dII.*U,taus,weights,P,N,J,2),A)/rho
  end
  if approx.aL3!=0. && approx.gL3!=0.
    dG-=approx.aL3*conv_chebyshev(FU,chebyshev(dJJ/2,taus,weights,P,N,J,2),A)/rho
  end

  dsK=zeros(Float64,N*J)
  if approx.gK!=0.
    if approx.bK!=0.
      dsK+=approx.gK*approx.bK*dS
    end
    if approx.bK!=1.
      dsK+=approx.gK*(1-approx.bK)*dE*ones(N*J)
    end
  end
  dsL1=zeros(Float64,N*J)
  if approx.bL1!=0.
    dsL1+=approx.bL1*dS
  end
  if approx.bL1!=1.
    dsL1+=(1-approx.bL1)*dE*ones(N*J)
  end

  dX=(dsK-dsL1+dII)./sL1-X./sL1.*dsL1
  dY=(dS-dsL1+dG+dJJ/2)./sL1-Y./sL1.*dsL1

  out=((Y.+1).*dX./(X.+1)-dY)./(2*(X.+1)).*dotPhi((Y.+1)./((X.+1).^2))+(Y.+1)./(2*(X.+1).^3).*(2*(Y.+1)./(X.+1).*dX-dY).*dotdPhi((Y.+1)./(X.+1).^2)
  return out
end