## 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 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
function anyeq_energy(
  rho::Float64,
  minlrho_init::Float64,
  nlrho_init::Int64,
  taus::Array{Float64,1},
  P::Int64,
  N::Int64,
  J::Int64,
  a0::Float64,
  v::Function,
  maxiter::Int64,
  tolerance::Float64,
  approx::Anyeq_approx,
  savefile::String
)
  # initialize vectors
  (weights,T,k,V,V0,A,Abar,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v,approx,savefile)

  # compute initial guess from medeq
  u0=anyeq_init_medeq([rho],minlrho_init,nlrho_init,N,J,k,a0,v,maxiter,tolerance)

  (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::Array{Float64,1},
  minlrho_init::Float64,
  nlrho_init::Int64,
  taus::Array{Float64,1},
  P::Int64,
  N::Int64,
  J::Int64,
  a0::Float64,
  v::Function,
  maxiter::Int64,
  tolerance::Float64,
  approx::Anyeq_approx,
  savefile::String
)
  # initialize vectors
  (weights,T,k,V,V0,A,Abar,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v,approx,savefile)

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

  # spawn workers
  work=spawn_workers(length(rhos))

  # compute u
  (us,es,errs)=anyeq_hatu_rho_multithread(u0s,P,N,J,rhos,a0,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,v,maxiter,tolerance,approx,work)

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

# compute energy as a function of rho
# initialize with previous rho
function anyeq_energy_rho_init_prevrho(
  rhos::Array{Float64,1},
  minlrho_init::Float64,
  nlrho_init::Int64,
  taus::Array{Float64,1},
  P::Int64,
  N::Int64,
  J::Int64,
  a0::Float64,
  v::Function,
  maxiter::Int64,
  tolerance::Float64,
  approx::Anyeq_approx,
  savefile::String
)
  # initialize vectors
  (weights,T,k,V,V0,A,Abar,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v,approx,savefile)

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

  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::Array{Float64,1},
  minlrho_init::Float64,
  nlrho_init::Int64,
  taus::Array{Float64,1},
  P::Int64,
  N::Int64,
  J::Int64,
  a0::Float64,
  v::Function,
  maxiter::Int64,
  tolerance::Float64,
  approx::Anyeq_approx,
  savefile::String
)
  # initialize vectors
  (weights,T,k,V,V0,A,Abar,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v,approx,savefile)

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

  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)
function anyeq_uk(
  minlrho_init::Float64,
  nlrho_init::Int64,
  taus::Array{Float64,1},
  P::Int64,
  N::Int64,
  J::Int64,
  rho::Float64,
  a0::Float64,
  v::Function,
  maxiter::Int64,
  tolerance::Float64,
  approx::Anyeq_approx,
  savefile::String
)
  # init vectors
  (weights,T,k,V,V0,A,Abar,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v,approx,savefile)

  # compute initial guess from medeq
  u0=anyeq_init_medeq([rho],minlrho_init,nlrho_init,N,J,k,a0,v,maxiter,tolerance)

  (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)
function anyeq_ux(
  minlrho_init::Float64,
  nlrho_init::Int64,
  taus::Array{Float64,1},
  P::Int64,
  N::Int64,
  J::Int64,
  rho::Float64,
  a0::Float64,
  v::Function,
  maxiter::Int64,
  tolerance::Float64,
  xmin::Float64,
  xmax::Float64,
  nx::Int64,
  approx::Anyeq_approx,
  savefile::String
)
  # init vectors
  (weights,T,k,V,V0,A,Abar,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v,approx,savefile)

  # compute initial guess from medeq
  u0=anyeq_init_medeq([rho],minlrho_init,nlrho_init,N,J,k,a0,v,maxiter,tolerance)

  (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=anyeq_u_x(x,u,k,N,J,taus,weights)
    @printf("% .15e % .15e % .15e\n",x,ux,)
  end
end

# compute condensate fraction for a given rho
function anyeq_condensate_fraction(
  rho::Float64,
  minlrho_init::Float64,
  nlrho_init::Int64,
  taus::Array{Float64,1},
  P::Int64,
  N::Int64,
  J::Int64,
  a0::Float64,
  v::Function,
  maxiter::Int64,
  tolerance::Float64,
  approx::Anyeq_approx,
  savefile::String
)
  # initialize vectors
  (weights,T,k,V,V0,A,Abar,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v,approx,savefile)

  # compute initial guess from medeq
  u0=anyeq_init_medeq([rho],minlrho_init,nlrho_init,N,J,k,a0,v,maxiter,tolerance)

  (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::Array{Float64,1},
  minlrho_init::Float64,
  nlrho_init::Int64,
  taus::Array{Float64,1},
  P::Int64,
  N::Int64,
  J::Int64,
  a0::Float64,
  v::Function,
  maxiter::Int64,
  tolerance::Float64,
  approx::Anyeq_approx,
  savefile::String
)
  # spawn workers
  work=spawn_workers(length(rhos))

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

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

  # compute u
  (us,es,errs)=anyeq_hatu_rho_multithread(u0s,P,N,J,rhos,a0,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,v,maxiter,tolerance,approx,work)

  # compute eta
  etas=Array{Float64,1}(undef,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
      etas[j]=remotecall_fetch(anyeq_eta,workers()[p],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
function anyeq_momentum_distribution(
  kmin::Float64,
  kmax::Float64,
  rho::Float64,
  minlrho_init::Float64,
  nlrho_init::Int64,
  taus::Array{Float64,1},
  P::Int64,
  N::Int64,
  J::Int64,
  windowL::Float64,
  a0::Float64,
  v::Function,
  maxiter::Int64,
  tolerance::Float64,
  approx::Anyeq_approx,
  savefile::String
)
  # initialize vectors
  (weights,T,k,V,V0,A,Abar,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v,approx,savefile)

  # compute initial guess from medeq
  u0=anyeq_init_medeq([rho],minlrho_init,nlrho_init,N,J,k,a0,v,maxiter,tolerance)

  (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
  if windowL==Inf
    # use discrete approximation of delta function
    M=anyeq_momentum_discrete_delta(kmin,kmax,u,P,N,J,rho,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,approx)
  else
    M=anyeq_momentum_window(kmin,kmax,u,P,N,J,windowL,rho,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,approx)
  end

  # order k's in increasing order
  for zeta in 0:J-1
    for j in 1:N
      q=k[(J-1-zeta)*N+j]
      # drop if not in k interval
      if q<kmin || q>kmax
	continue
      end
      @printf("% .15e % .15e\n",q,M[(J-1-zeta)*N+j])
    end
  end
end

# 2 point correlation function
function anyeq_2pt_correlation(
  minlrho_init::Float64,
  nlrho_init::Int64,
  taus::Array{Float64,1},
  P::Int64,
  N::Int64,
  J::Int64,
  windowL::Float64,
  rho::Float64,
  a0::Float64,
  v::Function,
  maxiter::Int64,
  tolerance::Float64,
  xmin::Float64,
  xmax::Float64,
  nx::Int64,
  approx::Anyeq_approx,
  savefile::String
)
  # init vectors
  (weights,T,k,V,V0,A,Abar,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v,approx,savefile)

  # compute initial guess from medeq
  u0=anyeq_init_medeq([rho],minlrho_init,nlrho_init,N,J,k,a0,v,maxiter,tolerance)

  # 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)
    @printf("% .15e % .15e\n",x,C2)
  end
end

# maximum of 2 point correlation function
function anyeq_2pt_correlation_max(
  rho::Float64,
  minlrho_init::Float64,
  nlrho_init::Int64,
  dx::Float64,
  x0::Float64, # initial guess is x0/rho^(1/3)
  maxstep::Float64,
  taus::Array{Float64,1},
  P::Int64,
  N::Int64,
  J::Int64,
  windowL::Float64,
  a0::Float64,
  v::Function,
  maxiter::Int64,
  tolerance::Float64,
  tolerance_max::Float64,
  approx::Anyeq_approx,
  savefile::String
)
  # init vectors
  (weights,T,k,V,V0,A,Abar,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v,approx,savefile)

  # compute initial guess from medeq
  u0=anyeq_init_medeq([rho],minlrho_init,nlrho_init,N,J,k,a0,v,maxiter,tolerance)

  # initial guess
  x0=1/rho^(1/3)

  (u,E,error)=anyeq_hatu(u0,P,N,J,rho,a0,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,v,maxiter,tolerance,approx)
  (x,f)=anyeq_2pt_max(u,P,N,J,x0/rho^(1/3),dx,maxstep,maxiter,tolerance_max,windowL,rho,weights,T,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,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 anyeq_2pt_correlation_max_rho(
  rhos::Array{Float64,1},
  minlrho_init::Float64,
  nlrho_init::Int64,
  dx::Float64,
  x0::Float64, # initial guess is x0/rho^(1/3)
  maxstep::Float64,
  taus::Array{Float64,1},
  P::Int64,
  N::Int64,
  J::Int64,
  windowL::Float64,
  a0::Float64,
  v::Function,
  maxiter::Int64,
  tolerance::Float64,
  tolerance_max::Float64,
  approx::Anyeq_approx,
  savefile::String
)
  # init vectors
  (weights,T,k,V,V0,A,Abar,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v,approx,savefile)

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

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

  # compute u
  (us,es,errs)=anyeq_hatu_rho_multithread(u0s,P,N,J,rhos,a0,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,v,maxiter,tolerance,approx,work)

  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(anyeq_2pt_max,workers()[p],us[j],P,N,J,x0/rhos[j]^(1/3),dx,maxstep,maxiter,tolerance_max,windowL,rhos[j],weights,T,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,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

# Correlation function in Fourier space
function anyeq_2pt_correlation_fourier(
  minlrho_init::Float64,
  nlrho_init::Int64,
  taus::Array{Float64,1},
  P::Int64,
  N::Int64,
  J::Int64,
  windowL::Float64,
  rho::Float64,
  a0::Float64,
  v::Function,
  maxiter::Int64,
  tolerance::Float64,
  kmin::Float64,
  kmax::Float64,
  nk::Int64,
  approx::Anyeq_approx,
  savefile::String
)
  # init vectors
  (weights,T,k,V,V0,A,Abar,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v,approx,savefile)

  # compute initial guess from medeq
  u0=anyeq_init_medeq([rho],minlrho_init,nlrho_init,N,J,k,a0,v,maxiter,tolerance)

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

  # save result from each task
  C2s=Array{Float64,1}(undef,nk)

  # spawn workers
  work=spawn_workers(nk)

  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,nk,10000)
      end
      # run the task
      C2s[j]=remotecall_fetch(anyeq_2pt_fourier,workers()[p],kmin+(kmax-kmin)*j/nk,u,P,N,J,windowL,rho,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,approx,S,E,II,JJ,X,Y,sL1,sK)
    end
  end

  for j in 1:nk
    @printf("% .15e % .15e\n",kmin+(kmax-kmin)*j/nk,C2s[j])
  end
end

# maximum of Fourier transform of 2 point correlation function
function anyeq_2pt_correlation_fourier_max(
  rho::Float64,
  minlrho_init::Float64,
  nlrho_init::Int64,
  dk::Float64,
  k0::Float64, # initial guess is k0*rho^(1/3)
  maxstep::Float64,
  taus::Array{Float64,1},
  P::Int64,
  N::Int64,
  J::Int64,
  windowL::Float64,
  a0::Float64,
  v::Function,
  maxiter::Int64,
  tolerance::Float64,
  tolerance_max::Float64,
  approx::Anyeq_approx,
  savefile::String
)
  # init vectors
  (weights,T,k,V,V0,A,Abar,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v,approx,savefile)

  # compute initial guess from medeq
  u0=anyeq_init_medeq([rho],minlrho_init,nlrho_init,N,J,k,a0,v,maxiter,tolerance)

  # compute u
  (u,E,error)=anyeq_hatu(u0,P,N,J,rho,a0,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,v,maxiter,tolerance,approx)
  (ko,f)=anyeq_2pt_fourier_max(u,P,N,J,k0*rho^(1/3),dk,maxstep,maxiter,tolerance_max,windowL,rho,weights,T,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,approx)

  if(ko==Inf)
    @printf(stderr,"max search failed for rho=%e\n",rho)
  else
    @printf("% .15e % .15e\n",ko,f)
  end
end

# maximum of fourier transform of 2 point correlation function as a function of rho
function anyeq_2pt_correlation_fourier_max_rho(
  rhos::Array{Float64,1},
  minlrho_init::Float64,
  nlrho_init::Int64,
  dk::Float64,
  k0::Float64, # initial guess is k0*rho^(1/3)
  maxstep::Float64,
  taus::Array{Float64,1},
  P::Int64,
  N::Int64,
  J::Int64,
  windowL::Float64,
  a0::Float64,
  v::Function,
  maxiter::Int64,
  tolerance::Float64,
  tolerance_max::Float64,
  approx::Anyeq_approx,
  savefile::String
)
  # init vectors
  (weights,T,k,V,V0,A,Abar,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v,approx,savefile)

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

  # spawn workers
  work=spawn_workers(length(rhos))

  # compute u
  (us,es,errs)=anyeq_hatu_rho_multithread(u0s,P,N,J,rhos,a0,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,v,maxiter,tolerance,approx,work)

  # save result from each task
  ks=Array{Float64,1}(undef,length(rhos))
  fs=Array{Float64,1}(undef,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(anyeq_2pt_fourier_max,workers()[p],us[j],P,N,J,k0*rhos[j]^(1/3),dk,maxstep,maxiter,tolerance_max,windowL,rhos[j],weights,T,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,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

# uncondensed 2 point correlation function
function anyeq_uncondensed_2pt_correlation(
  minlrho_init::Float64,
  nlrho_init::Int64,
  taus::Array{Float64,1},
  P::Int64,
  N::Int64,
  J::Int64,
  windowL::Float64,
  rho::Float64,
  a0::Float64,
  v::Function,
  maxiter::Int64,
  tolerance::Float64,
  xmin::Float64,
  xmax::Float64,
  nx::Int64,
  approx::Anyeq_approx,
  savefile::String
)
  # init vectors
  (weights,T,k,V,V0,A,Abar,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v,approx,savefile)

  # compute initial guess from medeq
  u0=anyeq_init_medeq([rho],minlrho_init,nlrho_init,N,J,k,a0,v,maxiter,tolerance)

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

  # compute u(xi)
  uxis=Array{Float64,1}(undef,nx)
  for i in 1:nx
    uxis[i]=anyeq_u_x(xmin+(xmax-xmin)*i/nx,u,k,N,J,taus,weights)
  end
    

  # spawn workers
  work=spawn_workers(nx)
  gamma2s=Array{Float64,1}(undef,nx)
  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,nx,10000)
      end
      # run the task
      gamma2s[j]=remotecall_fetch(anyeq_uncondensed_2pt,workers()[p],xmin+(xmax-xmin)*j/nx,uxis[j],u,P,N,J,windowL,rho,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,approx,S,E,II,JJ,X,Y,sL1,sK)
    end
  end

  for i in 1:nx
    @printf("% .15e % .15e\n",xmin+(xmax-xmin)*i/nx,gamma2s[i])
  end
end

# compute the compressibility
function anyeq_compressibility_rho(
  rhos::Array{Float64,1},
  minlrho_init::Float64,
  nlrho_init::Int64,
  taus::Array{Float64,1},
  P::Int64,
  N::Int64,
  J::Int64,
  a0::Float64,
  v::Function,
  maxiter::Int64,
  tolerance::Float64,
  approx::Anyeq_approx,
  savefile::String
)
  # initialize vectors
  (weights,T,k,V,V0,A,Abar,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v,approx,savefile)

  # compute initial guess from medeq
  u0s=anyeq_init_medeq(rhos,minlrho_init,nlrho_init,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
  work=spawn_workers(length(rhos))

  # compute u
  (us,es,errs)=anyeq_hatu_rho_multithread(u0s,P,N,J,rhos,a0,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,v,maxiter,tolerance,approx,work)

  for j in 1:length(rhos)-2
    # compute \rho^2\partial^2(\rho e)=\partial^2_{\log\rho}(\rho e)-\partial_{\log\rho}(\rho e) as a function of rho
    # \partial^2_{\log\rho}(\rho e)
    p2=(rhos[j+2]*es[j+2]-2*rhos[j+1]*es[j+1]+rhos[j]*es[j])/(log(rhos[j+2])-log(rhos[j+1]))/(log(rhos[j+1])-log(rhos[j]))
    # \partial_{\log\rho}(\rho e) (take average of front and back)
    p11=(rhos[j+2]*es[j+2]-rhos[j+1]*es[j+1])/(log(rhos[j+2])-log(rhos[j+1]))
    p12=(rhos[j+1]*es[j+1]-rhos[j]*es[j])/(log(rhos[j+1])-log(rhos[j]))
    # compressibility is 1/this
    @printf("% .15e % .15e\n",rhos[j+1],1/(p2-(p11+p12)/2))
  end
end

# Fourier transform 2 point correlation function test: compute by transforming anyeq_2pt
function anyeq_2pt_correlation_fourier_test(
  minlrho_init::Float64,
  nlrho_init::Int64,
  taus::Array{Float64,1},
  P::Int64,
  N::Int64,
  J::Int64,
  windowL::Float64,
  rho::Float64,
  a0::Float64,
  v::Function,
  maxiter::Int64,
  tolerance::Float64,
  xmax::Float64,
  kmin::Float64,
  kmax::Float64,
  nk::Int64,
  approx::Anyeq_approx,
  savefile::String
)
  # init vectors
  (weights,T,k,V,V0,A,Abar,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v,approx,savefile)

  # compute initial guess from medeq
  u0=anyeq_init_medeq([rho],minlrho_init,nlrho_init,N,J,k,a0,v,maxiter,tolerance)

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

  # compute C2 in real space
  C2=Array{Float64,1}(undef,N*J)
  for i in 1:N*J
    if k[i]<xmax
      C2[i]=anyeq_2pt(k[i],u,P,N,J,windowL,rho,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,approx,S,E,II,JJ,X,Y,sL1,sK)-rho^2
    else
      C2[i]=0.
    end
  end
  # invert Fourier transform
  for i in 1:nk
    k0=kmin+(kmax-kmin)*i/nk
    hatC2=inverse_fourier_chebyshev(C2,k0,k,taus,weights,N,J)
    @printf("% .15e % .15e\n",k0,hatC2)
  end
end

# compute Abar
function anyeq_Abar_multithread(
  k::Array{Float64,1},
  weights::Tuple{Array{Float64,1},Array{Float64,1}},
  taus::Array{Float64,1},
  T::Array{Polynomial,1},
  P::Int64,
  N::Int64,
  J::Int64,
  approx::Anyeq_approx
)
  if approx.bL3==0.
    # empty array
    return Array{Array{Float64,5}}(undef,0)
  end

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

  # spawn workers
  work=spawn_workers(N*J)

  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,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::Array{Float64,1},
  P::Int64,
  N::Int64,
  J::Int64,
  v::Function,
  approx::Anyeq_approx,
  savefile::String
)
  # Gauss-Legendre weights
  weights=gausslegendre(N)

  # initialize vectors V,k
  V=Array{Float64,1}(undef,J*N)
  k=Array{Float64,1}(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)

  # 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

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

# compute initial guess from medeq
@everywhere function anyeq_init_medeq(
  rhos::Array{Float64,1},
  minlrho_init::Float64,
  nlrho_init::Int64,
  N::Int64,
  J::Int64,
  k::Array{Float64,1},
  a0::Float64,
  v::Function,
  maxiter::Int64,
  tolerance::Float64
)
  weights_medeq=gausslegendre(N*J)
  (u0s_medeq,es,errs)=easyeq_compute_u_prevrho(rhos,minlrho_init,nlrho_init,a0,N*J,v,maxiter,tolerance,weights_medeq,Easyeq_approx(1.,1.))

  # check errs
  for j in 1:length(errs)
    if errs[j]>tolerance
      print(stderr,"warning: computation of initial Ansatz failed for rho=",rhos[j],"\n")
    end
  end

  # return a single vector if there is a single rho
  if length(rhos)>1
    u0s=Array{Array{Float64,1}}(undef,length(rhos))
    for j in 1:length(u0s_medeq)
      u0s[j]=easyeq_to_anyeq(u0s_medeq[j],weights_medeq,k,N,J)
    end
    return u0s
  else
    return easyeq_to_anyeq(u0s_medeq,weights_medeq,k,N,J)
  end
end

# interpolate the solution of medeq to an input for anyeq
@everywhere function easyeq_to_anyeq(
  u_simple::Array{Float64,1},
  weights::Tuple{Array{Float64,1},Array{Float64,1}},
  k::Array{Float64,1},
  N::Int64,
  J::Int64
)
  # 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,1}(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::Array{Float64,1},
  P::Int64,
  N::Int64,
  J::Int64,
  rho::Float64,
  a0::Float64,
  weights::Tuple{Array{Float64,1},Array{Float64,1}},
  k::Array{Float64,1},
  taus::Array{Float64,1},
  V::Array{Float64,1},
  V0::Float64,
  A::Array{Array{Float64,2},1},
  Abar::Array{Array{Float64,5},1},
  Upsilon::Array{Array{Float64,1},1},
  Upsilon0::Array{Float64,1},
  v::Function,
  maxiter::Int64,
  tolerance::Float64,
  approx::Anyeq_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,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


# compute u for various rho
function anyeq_hatu_rho_multithread(
  u0s::Array{Array{Float64,1},1},
  P::Int64,
  N::Int64,
  J::Int64,
  rhos::Array{Float64,1},
  a0::Float64,
  weights::Tuple{Array{Float64,1},Array{Float64,1}},
  k::Array{Float64,1},
  taus::Array{Float64,1},
  V::Array{Float64,1},
  V0::Float64,
  A::Array{Array{Float64,2},1},
  Abar::Array{Array{Float64,5},1},
  Upsilon::Array{Array{Float64,1},1},
  Upsilon0::Array{Float64,1},
  v::Function,
  maxiter::Int64,
  tolerance::Float64,
  approx::Anyeq_approx,
  work::Array{Array{Int64,1},1}
)
  # compute u
  us=Array{Array{Float64,1}}(undef,length(rhos))
  es=Array{Float64,1}(undef,length(rhos))
  errs=Array{Float64,1}(undef,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
      (us[j],es[j],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

  return (us,es,errs)
end


# save Abar
function anyeq_save_Abar(
  taus::Array{Float64,1},
  P::Int64,
  N::Int64,
  J::Int64,
  v::Function,
  approx::Anyeq_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::String,
  P::Int64,
  N::Int64,
  J::Int64
)
   # 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::Array{Float64,1},
  X::Array{Float64,1},
  Y::Array{Float64,1}
)
  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::Array{Float64,1},
  rho::Float64,
  k::Array{Float64,1},
  taus::Array{Float64,1},
  A::Array{Array{Float64,2},1},
  Abar::Array{Array{Float64,5},1},
  Upsilon::Array{Array{Float64,1},1},
  Upsilon0::Array{Float64,1},
  weights::Tuple{Array{Float64,1},Array{Float64,1}},
  P::Int64,
  N::Int64,
  J::Int64,
  approx::Anyeq_approx,
  S::Array{Float64,1},
  E::Float64,
  II::Array{Float64,1},
  JJ::Array{Float64,1},
  X::Array{Float64,1},
  Y::Array{Float64,1},
  sL1::Array{Float64,1},
  sK::Array{Float64,1}
)
  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::Array{Float64,1},
  rho::Float64,
  k::Array{Float64,1},
  taus::Array{Float64,1},
  v::Array{Float64,1},
  v0::Float64,
  A::Array{Array{Float64,2},1},
  Abar::Array{Array{Float64,5},1},
  Upsilon::Array{Array{Float64,1},1},
  Upsilon0::Array{Float64,1},
  weights::Tuple{Array{Float64,1},Array{Float64,1}},
  P::Int64,
  N::Int64,
  J::Int64,
  approx::Anyeq_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


# u(x)
@everywhere function anyeq_u_x(
  x::Float64,
  u::Array{Float64,1},
  k::Array{Float64,1},
  N::Int64,
  J::Int64,
  taus::Array{Float64,1},
  weights::Tuple{Array{Float64,1},Array{Float64,1}}
)
  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
  return real(ux)
end

# condensate fraction
@everywhere function anyeq_eta(
  u::Array{Float64,1},
  P::Int64,
  N::Int64,
  J::Int64,
  rho::Float64,
  weights::Tuple{Array{Float64,1},Array{Float64,1}},
  k::Array{Float64,1},
  taus::Array{Float64,1},
  V::Array{Float64,1},
  V0::Float64,
  A::Array{Array{Float64,2},1},
  Abar::Array{Array{Float64,5},1},
  Upsilon::Array{Array{Float64,1},1},
  Upsilon0::Array{Float64,1},
  approx::Anyeq_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,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


# correlation function
@everywhere function anyeq_2pt(
  x::Float64,
  u::Array{Float64,1},
  P::Int64,
  N::Int64,
  J::Int64,
  windowL::Float64,
  rho::Float64,
  weights::Tuple{Array{Float64,1},Array{Float64,1}},
  k::Array{Float64,1},
  taus::Array{Float64,1},
  V::Array{Float64,1},
  V0::Float64,
  A::Array{Array{Float64,2},1},
  Abar::Array{Array{Float64,5},1},
  Upsilon::Array{Array{Float64,1},1},
  Upsilon0::Array{Float64,1},
  approx::Anyeq_approx,
  S::Array{Float64,1},
  E::Float64,
  II::Array{Float64,1},
  JJ::Array{Float64,1},
  X::Array{Float64,1},
  Y::Array{Float64,1},
  sL1::Array{Float64,1},
  sK::Array{Float64,1}
)
  g=(r,x)->sinc(r*x)*hann(r,windowL)
  du=anyeq_dudv(g, x, u, P, N, J, rho, weights, k, taus, V, V0, A, Abar, Upsilon, Upsilon0, approx, S, E, II, JJ, X, Y, sL1, sK)

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

# uncondensed correlation function
@everywhere function anyeq_uncondensed_2pt(
  xi::Float64,
  uxi::Float64,
  u::Array{Float64,1},
  P::Int64,
  N::Int64,
  J::Int64,
  windowL::Float64,
  rho::Float64,
  weights::Tuple{Array{Float64,1},Array{Float64,1}},
  k::Array{Float64,1},
  taus::Array{Float64,1},
  V::Array{Float64,1},
  V0::Float64,
  A::Array{Array{Float64,2},1},
  Abar::Array{Array{Float64,5},1},
  Upsilon::Array{Array{Float64,1},1},
  Upsilon0::Array{Float64,1},
  approx::Anyeq_approx,
  S::Array{Float64,1},
  E::Float64,
  II::Array{Float64,1},
  JJ::Array{Float64,1},
  X::Array{Float64,1},
  Y::Array{Float64,1},
  sL1::Array{Float64,1},
  sK::Array{Float64,1}
)
  # compute dXi/dmu
  g=Array{Float64,1}(undef,length(k))
  for i in 1:length(k)
    g[i]=-2*rho*uxi*sinc(k[i]*xi)*hann(k[i],windowL)
  end
  dXidmu=-g./sL1./(2*(X.+1)).*dotPhi((Y.+1)./((X.+1).^2))-(Y.+1).*g./sL1./(2(X.+1).^3).*dotdPhi((Y.+1)./(X.+1).^2)

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

  return gamma2
end

# compute the directional derivative of u with respect to v in direction g
@everywhere function anyeq_dudv(
  g::Function,# should be of the form g(k,x) where x is a parameter
  x::Float64,
  u::Array{Float64,1},
  P::Int64,
  N::Int64,
  J::Int64,
  rho::Float64,
  weights::Tuple{Array{Float64,1},Array{Float64,1}},
  k::Array{Float64,1},
  taus::Array{Float64,1},
  V::Array{Float64,1},
  V0::Float64,
  A::Array{Array{Float64,2},1},
  Abar::Array{Array{Float64,5},1},
  Upsilon::Array{Array{Float64,1},1},
  Upsilon0::Array{Float64,1},
  approx::Anyeq_approx,
  S::Array{Float64,1},
  E::Float64,
  II::Array{Float64,1},
  JJ::Array{Float64,1},
  X::Array{Float64,1},
  Y::Array{Float64,1},
  sL1::Array{Float64,1},
  sK::Array{Float64,1}
)
  # initialize dV
  dV=Array{Float64,1}(undef,J*N)
  for i in 1:J*N
    dV[i]=g(k[i],x)
  end
  dV0=g(0.,x)

  # compute dUpsilon
  # Upsilonmat does not use splines, so increase precision
  weights_plus=gausslegendre(N*J)
  dUpsilon=Upsilonmat(k,s->g(s,x),weights_plus)
  dUpsilon0=Upsilon0mat(k,s->g(s,x),weights_plus)

  du=-inv(anyeq_DXi(rho*u,rho,k,taus,A,Abar,Upsilon,Upsilon0,weights,P,N,J,approx,S,E,II,JJ,X,Y,sL1,sK))*anyeq_dXidv(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

  return du
end

# derivative of Xi with respect to v in the specified by dUpsilon and dUpsilon0
@everywhere function anyeq_dXidv(
  U::Array{Float64,1},
  rho::Float64,
  k::Array{Float64,1},
  taus::Array{Float64,1},
  dv::Array{Float64,1},
  dv0::Float64,
  A::Array{Array{Float64,2},1},
  Abar::Array{Array{Float64,5},1},
  dUpsilon::Array{Array{Float64,1},1},
  dUpsilon0::Array{Float64,1},
  weights::Tuple{Array{Float64,1},Array{Float64,1}},
  P::Int64,
  N::Int64,
  J::Int64,
  approx::Anyeq_approx,
  S::Array{Float64,1},
  E::Float64,
  II::Array{Float64,1},
  JJ::Array{Float64,1},
  X::Array{Float64,1},
  Y::Array{Float64,1},
  sL1::Array{Float64,1},
  sK::Array{Float64,1}
)
  # 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

# maximum of 2 point correlation function
@everywhere function anyeq_2pt_max(
  u::Array{Float64,1},
  P::Int64,
  N::Int64,
  J::Int64,
  x0::Float64,
  dx::Float64,
  maxstep::Float64,
  maxiter::Int64,
  tolerance::Float64,
  windowL::Float64,
  rho::Float64,
  weights::Tuple{Array{Float64,1},Array{Float64,1}},
  T::Array{Polynomial,1},
  k::Array{Float64,1},
  taus::Array{Float64,1},
  V::Array{Float64,1},
  V0::Float64,
  A::Array{Array{Float64,2},1},
  Abar::Array{Array{Float64,5},1},
  Upsilon::Array{Array{Float64,1},1},
  Upsilon0::Array{Float64,1},
  approx::Anyeq_approx
)
  # compute some useful integrals
  (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)

  (x,f)=newton_maximum(y->anyeq_2pt(y,u,P,N,J,windowL,rho,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,approx,S,E,II,JJ,X,Y,sL1,sK),x0,dx,maxiter,tolerance,maxstep)

  return(x,f)
end


# Fourier transform of 2pt correlation function at q
@everywhere function anyeq_2pt_fourier(
  q::Float64,
  u::Array{Float64,1},
  P::Int64,
  N::Int64,
  J::Int64,
  windowL::Float64,
  rho::Float64,
  weights::Tuple{Array{Float64,1},Array{Float64,1}},
  k::Array{Float64,1},
  taus::Array{Float64,1},
  V::Array{Float64,1},
  V0::Float64,
  A::Array{Array{Float64,2},1},
  Abar::Array{Array{Float64,5},1},
  Upsilon::Array{Array{Float64,1},1},
  Upsilon0::Array{Float64,1},
  approx::Anyeq_approx,
  S::Array{Float64,1},
  E::Float64,
  II::Array{Float64,1},
  JJ::Array{Float64,1},
  X::Array{Float64,1},
  Y::Array{Float64,1},
  sL1::Array{Float64,1},
  sK::Array{Float64,1}
)
  # direction in which to differentiate u
  weights_plus=gausslegendre(N*J)
  #g=(r,x)->(r>0. ? (1.)/(2*x*r)*integrate_legendre(s->s*gaussian(s,(1.)/windowL),abs(x-r),x+r,weights_plus) : gaussian(x,(1.)/windowL))
  g=(r,x)->(r>0. ? (1.)/(2*x*r*windowL)*(gaussian(x-r,(1.)/windowL)-gaussian(x+r,(1.)/windowL)) : gaussian(x,(1.)/windowL))

  du=anyeq_dudv(g,q,u,P,N,J,rho,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,approx,S,E,II,JJ,X,Y,sL1,sK)

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

  return C2
end

# maximum of Fourier transform of 2 point correlation function
@everywhere function anyeq_2pt_fourier_max(
  u::Array{Float64,1},
  P::Int64,
  N::Int64,
  J::Int64,
  k0::Float64,
  dk::Float64,
  maxstep::Float64,
  maxiter::Int64,
  tolerance::Float64,
  windowL::Float64,
  rho::Float64,
  weights::Tuple{Array{Float64,1},Array{Float64,1}},
  T::Array{Polynomial,1},
  k::Array{Float64,1},
  taus::Array{Float64,1},
  V::Array{Float64,1},
  V0::Float64,
  A::Array{Array{Float64,2},1},
  Abar::Array{Array{Float64,5},1},
  Upsilon::Array{Array{Float64,1},1},
  Upsilon0::Array{Float64,1},
  approx::Anyeq_approx
)
  # compute some useful integrals
  (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)

  (ko,f)=newton_maximum(y->anyeq_2pt_fourier(y,u,P,N,J,windowL,rho,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,approx,S,E,II,JJ,X,Y,sL1,sK),k0,dk,maxiter,tolerance,maxstep)

  return(ko,f)
end

# momentum distribution, computed using a Gaussian window
@everywhere function anyeq_momentum_window(
  kmin::Float64,
  kmax::Float64,
  u::Array{Float64,1},
  P::Int64,
  N::Int64,
  J::Int64,
  windowL::Float64, # L is windowL/k^2
  rho::Float64,
  weights::Tuple{Array{Float64,1},Array{Float64,1}},
  k::Array{Float64,1},
  taus::Array{Float64,1},
  V::Array{Float64,1},
  V0::Float64,
  A::Array{Array{Float64,2},1},
  Abar::Array{Array{Float64,5},1},
  Upsilon::Array{Array{Float64,1},1},
  Upsilon0::Array{Float64,1},
  approx::Anyeq_approx
)

  # compute dXi/dlambda without the delta function of u(q)
  (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=(16*pi^3)*(dotPhi((Y.+1)./((X.+1).^2))./(2*(X.+1))+(Y.+1)./(2*(X.+1).^3).*dotdPhi((Y.+1)./(X.+1).^2))./sL1
  
  # compute dXidu
  dXidu=inv(anyeq_DXi(rho*u,rho,k,taus,A,Abar,Upsilon,Upsilon0,weights,P,N,J,approx,S,E,II,JJ,X,Y,sL1,sK))

  M=Array{Float64,1}(undef,N*J)
  for j in 1:J*N
    # the Chebyshev polynomial expansion is often not good enough to compute M away from k[i]
    q=k[j]

    # drop the computation if not in k interval
    if q<kmin || q>kmax
      continue
    end

    # delta function
    delta=Array{Float64,1}(undef,J*N)
    L=windowL/q^2
    for i in 1:J*N
      delta[i]=(1.)/(2*k[i]*q*L)*(gaussian(k[i]-q,(1.)/L)-gaussian(k[i]+q,(1.)/L))
    end

    # du/dlambda
    du=-dXidu*(dXidlambda.*delta*u[j])
    # rescale u
    du=du/rho

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

  return M
end

# momentum distribution, computed using a discrete approximation of the delta function
@everywhere function anyeq_momentum_discrete_delta(
  kmin::Float64,
  kmax::Float64,
  u::Array{Float64,1},
  P::Int64,
  N::Int64,
  J::Int64,
  rho::Float64,
  weights::Tuple{Array{Float64,1},Array{Float64,1}},
  k::Array{Float64,1},
  taus::Array{Float64,1},
  V::Array{Float64,1},
  V0::Float64,
  A::Array{Array{Float64,2},1},
  Abar::Array{Array{Float64,5},1},
  Upsilon::Array{Array{Float64,1},1},
  Upsilon0::Array{Float64,1},
  approx::Anyeq_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,1}(undef,J*N)
  for zeta in 0:J-1
    for n in 1:N
      delta[zeta*N+n]=8/pi/((taus[zeta+2]-taus[zeta+1])*weights[2][n]*cos(pi*weights[1][n]/2)*(1+k[zeta*N+n])^2)
    end
  end
  
  # compute dXidu
  dXidu=inv(anyeq_DXi(rho*u,rho,k,taus,A,Abar,Upsilon,Upsilon0,weights,P,N,J,approx,S,E,II,JJ,X,Y,sL1,sK))

  M=Array{Float64,1}(undef,J*N)
  for i in 1:J*N
    # drop the computation if not in k interval
    if k[i]<kmin || k[i]>kmax
      continue
    end

    # 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