## 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 Easyeq_approx
  bK::Float64
  bL::Float64
end

# compute energy
function easyeq_energy(
  minlrho_init::Float64,
  nlrho_init::Int64,
  order::Int64,
  rho::Float64,
  a0::Float64,
  v::Function,
  maxiter::Int64,
  tolerance::Float64,
  approx::Easyeq_approx
)
  # compute gaussian quadrature weights
  weights=gausslegendre(order)

  # compute u
  (u,E,err)= easyeq_compute_u_prevrho([rho],minlrho_init,nlrho_init,a0,order,v,maxiter,tolerance,weights,approx)

  # print energy
  @printf("% .15e % .15e\n",real(E),err)
end

# compute energy as a function of rho
function easyeq_energy_rho(
  rhos::Array{Float64,1},
  minlrho_init::Float64,
  nlrho_init::Int64,
  order::Int64,
  a0::Float64,
  v::Function,
  maxiter::Int64,
  tolerance::Float64,
  approx::Easyeq_approx
)
  # compute gaussian quadrature weights
  weights=gausslegendre(order)
  # compute u
  (us,es,errs)= easyeq_compute_u_prevrho(rhos,minlrho_init,nlrho_init,a0,order,v,maxiter,tolerance,weights,approx)

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

# compute u(k)
function easyeq_uk(
  minlrho_init::Float64,
  nlrho_init::Int64,
  order::Int64,
  rho::Float64,
  a0::Float64,
  v::Function,
  maxiter::Int64,
  tolerance::Float64,
  approx::Easyeq_approx
)
  weights=gausslegendre(order)

  # compute u
  (u,E,err)= easyeq_compute_u_prevrho_error([rho],minlrho_init,nlrho_init,a0,order,v,maxiter,tolerance,weights,approx)

  for i in 1:order
    k=(1-weights[1][i])/(1+weights[1][i])
    @printf("% .15e % .15e\n",k,real(u[i]))
  end
end

# compute u(x)
function easyeq_ux(
  minlrho_init::Float64,
  nlrho_init::Int64,
  order::Int64,
  rho::Float64,
  a0::Float64,
  v::Function,
  maxiter::Int64,
  tolerance::Float64,
  xmin::Float64,
  xmax::Float64,
  nx::Int64,
  approx::Easyeq_approx
)
  weights=gausslegendre(order)

  # compute u
  (u,E,err)= easyeq_compute_u_prevrho_error([rho],minlrho_init,nlrho_init,a0,order,v,maxiter,tolerance,weights,approx)

  for i in 1:nx
    x=xmin+(xmax-xmin)*i/nx
    @printf("% .15e % .15e\n",x,real(easyeq_u_x(x,u,weights)))
  end
end

# compute 2u(x)-rho u*u(x)
function easyeq_uux(
  minlrho_init::Float64,
  nlrho_init::Int64,
  order::Int64,
  rho::Float64,
  a0::Float64,
  v::Function,
  maxiter::Int64,
  tolerance::Float64,
  xmin::Float64,
  xmax::Float64,
  nx::Int64,
  approx::Easyeq_approx
)
  weights=gausslegendre(order)

  # compute u
  (u,E,err)= easyeq_compute_u_prevrho_error([rho],minlrho_init,nlrho_init,a0,order,v,maxiter,tolerance,weights,approx)

  for i in 1:nx
    x=xmin+(xmax-xmin)*i/nx
    @printf("% .15e % .15e\n",x,real(easyeq_u_x(x,2*u-rho*u.*u,weights)))
  end
end

# condensate fraction
function easyeq_condensate_fraction(
  minlrho_init::Float64,
  nlrho_init::Int64,
  order::Int64,
  rho::Float64,
  a0::Float64,
  v::Function,
  maxiter::Int64,
  tolerance::Float64,
  approx::Easyeq_approx
)
  # compute gaussian quadrature weights
  weights=gausslegendre(order)

  # compute u
  (u,E,err)= easyeq_compute_u_prevrho([rho],minlrho_init,nlrho_init,a0,order,v,maxiter,tolerance,weights,approx)

  # compute eta
  eta=easyeq_eta(u,order,rho,v,maxiter,tolerance,weights,approx)

  # print energy
  @printf("% .15e % .15e\n",eta,err)
end

# condensate fraction as a function of rho
function easyeq_condensate_fraction_rho(
  rhos::Array{Float64,1},
  minlrho_init::Float64,
  nlrho_init::Int64,
  order::Int64,
  a0::Float64,
  v::Function,
  maxiter::Int64,
  tolerance::Float64,
  approx::Easyeq_approx
)
  weights=gausslegendre(order)
  # compute u
  (us,es,errs)= easyeq_compute_u_prevrho(rhos,minlrho_init,nlrho_init,a0,order,v,maxiter,tolerance,weights,approx)

  for j in 1:length(rhos)
    # compute eta
    eta=easyeq_eta(us[j],order,rhos[j],v,maxiter,tolerance,weights,approx)
    @printf("% .15e % .15e % .15e\n",rhos[j],eta,errs[j])
  end
end

# 2 pt correlation function
function easyeq_2pt(
  xmin::Float64,
  xmax::Float64,
  nx::Int64,
  minlrho_init::Float64,
  nlrho_init::Int64,
  order::Int64,
  windowL::Float64,
  rho::Float64,
  a0::Float64,
  v::Function,
  maxiter::Int64,
  tolerance::Float64,
  approx::Easyeq_approx
)
  # compute gaussian quadrature weights
  weights=gausslegendre(order)

  # compute u
  (u,E,err)= easyeq_compute_u_prevrho_error([rho],minlrho_init,nlrho_init,a0,order,v,maxiter,tolerance,weights,approx)

  # compute useful terms
  (V,V0)=easyeq_init_v(weights,v)
  (Eta,Eta0)=easyeq_init_H(weights,v)
  (E,S,A,T,B,X)=easyeq_ESATBX(rho*u,V,V0,Eta,Eta0,weights,rho,approx)

  # compute C2
  for j in 1:nx
    x=xmin+(xmax-xmin)/nx*j
    C2=easyeq_C2(x,u,windowL,rho,weights,Eta,Eta0,approx,E,S,A,T,B,X)
    @printf("% .15e % .15e\n",x,C2)
  end
end

# maximum of 2 point correlation function
function easyeq_2pt_max(
  dx::Float64,
  x0::Float64, # initial guess is x0/rho^(1/3)
  minlrho_init::Float64,
  nlrho_init::Int64,
  order::Int64,
  windowL::Float64,
  rho::Float64,
  a0::Float64,
  v::Function,
  maxstep::Float64,
  maxiter::Int64,
  tolerance::Float64,
  tolerance_max::Float64,
  approx::Easyeq_approx
)
  # compute gaussian quadrature weights
  weights=gausslegendre(order)

  # compute u
  (u,E,err)= easyeq_compute_u_prevrho_error([rho],minlrho_init,nlrho_init,a0,order,v,maxiter,tolerance,weights,approx)

  # compute useful terms
  (V,V0)=easyeq_init_v(weights,v)
  (Eta,Eta0)=easyeq_init_H(weights,v)

  (x,f)=easyeq_C2_max(u,x0/rho^(1/3),dx,maxstep,maxiter,tolerance_max,windowL,rho,weights,V,V0,Eta,Eta0,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 easyeq_2pt_max_rho(
  rhos::Array{Float64,1},
  dx::Float64,
  x0::Float64, # initial guess is x0/rho^(1/3)
  minlrho_init::Float64,
  nlrho_init::Int64,
  order::Int64,
  windowL::Float64,
  a0::Float64,
  v::Function,
  maxstep::Float64,
  maxiter::Int64,
  tolerance::Float64,
  tolerance_max::Float64,
  approx::Easyeq_approx
)
  # compute gaussian quadrature weights
  weights=gausslegendre(order)

  # compute u
  (us,es,errs)= easyeq_compute_u_prevrho_error(rhos,minlrho_init,nlrho_init,a0,order,v,maxiter,tolerance,weights,approx)

  # compute useful terms
  (V,V0)=easyeq_init_v(weights,v)
  (Eta,Eta0)=easyeq_init_H(weights,v)

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

  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(easyeq_C2_max,workers()[p],us[j],x0/rhos[j]^(1/3),dx,maxstep,maxiter,tolerance_max,windowL,rhos[j],weights,V,V0,Eta,Eta0,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


# Fourier transform of 2 pt correlation function
function easyeq_2pt_fourier(
  kmin::Float64,
  kmax::Float64,
  nk::Int64,
  minlrho_init::Float64,
  nlrho_init::Int64,
  order::Int64,
  windowL::Float64,
  rho::Float64,
  a0::Float64,
  v::Function,
  maxiter::Int64,
  tolerance::Float64,
  approx::Easyeq_approx
)
  # compute gaussian quadrature weights
  weights=gausslegendre(order)

  # compute u
  (u,E,err)= easyeq_compute_u_prevrho_error([rho],minlrho_init,nlrho_init,a0,order,v,maxiter,tolerance,weights,approx)

  # compute useful terms
  (V,V0)=easyeq_init_v(weights,v)
  (Eta,Eta0)=easyeq_init_H(weights,v)
  (E,S,A,T,B,X)=easyeq_ESATBX(rho*u,V,V0,Eta,Eta0,weights,rho,approx)

  # compute C2
  for j in 1:nk
    k=kmin+(kmax-kmin)/nk*j
    C2=easyeq_C2_fourier(k,u,windowL,rho,weights,Eta,Eta0,approx,E,S,A,T,B,X)
    @printf("% .15e % .15e\n",k,C2)
  end
end

# maximum of Fourier transform of 2 point correlation function
function easyeq_2pt_fourier_max(
  dk::Float64,
  k0::Float64, # initial guess is k0*rho^(1/3)
  minlrho_init::Float64,
  nlrho_init::Int64,
  order::Int64,
  windowL::Float64,
  rho::Float64,
  a0::Float64,
  v::Function,
  maxstep::Float64,
  maxiter::Int64,
  tolerance::Float64,
  tolerance_max::Float64,
  approx::Easyeq_approx
)
  # compute gaussian quadrature weights
  weights=gausslegendre(order)

  # compute u
  (u,E,err)= easyeq_compute_u_prevrho_error([rho],minlrho_init,nlrho_init,a0,order,v,maxiter,tolerance,weights,approx)

  # compute useful terms
  (V,V0)=easyeq_init_v(weights,v)
  (Eta,Eta0)=easyeq_init_H(weights,v)

  (k,f)=easyeq_C2_fourier_max(u,k0*rho^(1/3),dk,maxstep,maxiter,tolerance_max,windowL,rho,weights,V,V0,Eta,Eta0,approx)

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

# maximum of Fourier transform of 2 point correlation function as a function of rho
function easyeq_2pt_fourier_max_rho(
  rhos::Array{Float64,1},
  dk::Float64,
  k0::Float64, # initial guess is k0*rho^(1/3)
  minlrho_init::Float64,
  nlrho_init::Int64,
  order::Int64,
  windowL::Float64,
  a0::Float64,
  v::Function,
  maxstep::Float64,
  maxiter::Int64,
  tolerance::Float64,
  tolerance_max::Float64,
  approx::Easyeq_approx
)
  # compute gaussian quadrature weights
  weights=gausslegendre(order)

  # compute u
  (us,es,errs)= easyeq_compute_u_prevrho_error(rhos,minlrho_init,nlrho_init,a0,order,v,maxiter,tolerance,weights,approx)

  # compute useful terms
  (V,V0)=easyeq_init_v(weights,v)
  (Eta,Eta0)=easyeq_init_H(weights,v)

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

  # spawn workers
  work=spawn_workers(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(easyeq_C2_fourier_max,workers()[p],us[j],k0*rhos[j]^(1/3),dk,maxstep,maxiter,tolerance_max,windowL,rhos[j],weights,V,V0,Eta,Eta0,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

# momentum distribution
function easyeq_momentum_distribution(
  kmin::Float64,
  kmax::Float64,
  minlrho_init::Float64,
  nlrho_init::Int64,
  order::Int64,
  windowL::Float64, #L=windowL/k^2
  rho::Float64,
  a0::Float64,
  v::Function,
  maxiter::Int64,
  tolerance::Float64,
  approx::Easyeq_approx
)
  # compute gaussian quadrature weights
  weights=gausslegendre(order)

  # compute u
  (u,E,err)= easyeq_compute_u_prevrho_error([rho],minlrho_init,nlrho_init,a0,order,v,maxiter,tolerance,weights,approx)

  # compute useful terms
  (V,V0)=easyeq_init_v(weights,v)
  (Eta,Eta0)=easyeq_init_H(weights,v)
  (E,S,A,T,B,X)=easyeq_ESATBX(rho*u,V,V0,Eta,Eta0,weights,rho,approx)
  # dXi/dlambda without the delta function and u
  dXidlambda=(dotPhi(B.*T./((X.+1).^2))./(2*(X.+1))+B.*T./(2*(X.+1).^3).*dotdPhi(B.*T./(X.+1).^2))./A*(16*pi^3)

  dXi=inv(easyeq_dXi(rho*u,Eta,Eta0,weights,rho,approx,E,S,A,T,B,X))

  # compute momentum distribution
  for j in 1:order
    q=(1-weights[1][j])/(1+weights[1][j])
    # drop if not in k interval
    if q<kmin || q>kmax
      continue
    end

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

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

    # compute M
    M=-integrate_legendre_sampled(y->(1-y)/y^3,Eta0.*du,0.,1.,weights)*rho/(16*pi^3)

    @printf("% .15e % .15e\n",q,M)
  end
end


# initialize u from scattering solution
@everywhere function easyeq_init_u(
  a0::Float64,
  order::Int64,
  weights::Tuple{Array{Float64,1},Array{Float64,1}}
)
  u=zeros(Float64,order)
  for j in 1:order
    # transformed k
    k=(1-weights[1][j])/(1+weights[1][j])
    u[j]=4*pi*a0/k^2
  end

  return u
end

# compute u for an array of rhos
# use scattering solution for the first one, and the previous rho for the others
@everywhere function easyeq_compute_u_rho(
  rhos::Array{Float64,1},
  a0::Float64,
  order::Int64,
  v::Function,
  maxiter::Int64,
  tolerance::Float64,
  weights::Tuple{Array{Float64,1},Array{Float64,1}},
  approx::Easyeq_approx
)
  us=Array{Array{Float64,1}}(undef,length(rhos))
  es=Array{Float64,1}(undef,length(rhos))
  errs=Array{Float64,1}(undef,length(rhos))

  (us[1],es[1],errs[1])=easyeq_hatu(easyeq_init_u(a0,order,weights),order,rhos[1],v,maxiter,tolerance,weights,approx)
  for j in 2:length(rhos)
    (us[j],es[j],errs[j])=easyeq_hatu(us[j-1],order,rhos[j],v,maxiter,tolerance,weights,approx)
  end

  return (us,es,errs)
end

# compute u for an array of rhos
# start from a smaller rho and work up to rhos[1]
@everywhere function easyeq_compute_u_prevrho(
  rhos::Array{Float64,1},
  minlrho_init::Float64,
  nlrho_init::Int64,
  a0::Float64,
  order::Int64,
  v::Function,
  maxiter::Int64,
  tolerance::Float64,
  weights::Tuple{Array{Float64,1},Array{Float64,1}},
  approx::Easyeq_approx
)

  # only work up to rhos[1] if nlrho_init>0
  if nlrho_init>0
    rhos_init=Array{Float64,1}(undef,nlrho_init)
    for j in 0:nlrho_init-1
      rhos_init[j+1]=(nlrho_init==1 ? 10^minlrho_init : 10^(minlrho_init+(log10(rhos[1])-minlrho_init)/(nlrho_init-1)*j))
    end
    append!(rhos_init,rhos)
  # start from rhos[1] if nlrho_init=0
  else
    rhos_init=rhos
  end
  (us,es,errs)=easyeq_compute_u_rho(rhos_init,a0,order,v,maxiter,tolerance,weights,approx)

  # return a single value if there was a single input
  if length(rhos)==1
    return (us[nlrho_init+1],es[nlrho_init+1],errs[nlrho_init+1])
  else
    return (us[nlrho_init+1:length(us)],es[nlrho_init+1:length(es)],errs[nlrho_init+1:length(errs)])
  end
end
# with error message if the computation failed to be accurate enough
@everywhere function easyeq_compute_u_prevrho_error(
  rhos::Array{Float64,1},
  minlrho_init::Float64,
  nlrho_init::Int64,
  a0::Float64,
  order::Int64,
  v::Function,
  maxiter::Int64,
  tolerance::Float64,
  weights::Tuple{Array{Float64,1},Array{Float64,1}},
  approx::Easyeq_approx
)
  (us,es,errs)=easyeq_compute_u_prevrho(rhos,minlrho_init,nlrho_init,a0,order,v,maxiter,tolerance,weights,approx)
  # check errs
  for j in 1:length(errs)
    if errs[j]>tolerance
      print(stderr,"warning: computation of u failed for rho=",rhos[j],"\n")
    end
  end
  return (us,es,errs)
end


# \hat u(k) computed using Newton
@everywhere function easyeq_hatu(
  u0::Array{Float64,1},
  order::Int64,
  rho::Float64,
  v::Function,
  maxiter::Int64,
  tolerance::Float64,
  weights::Tuple{Array{Float64,1},Array{Float64,1}},
  approx::Easyeq_approx
)
  # initialize V and Eta
  (V,V0)=easyeq_init_v(weights,v)
  (Eta,Eta0)=easyeq_init_H(weights,v)

  # init u
  u=rho*u0

  # iterate
  err=Inf
  for i in 1:maxiter-1
    (E,S,A,T,B,X)=easyeq_ESATBX(u,V,V0,Eta,Eta0,weights,rho,approx)
    new=u-inv(easyeq_dXi(u,Eta,Eta0,weights,rho,approx,E,S,A,T,B,X))*easyeq_Xi(u,order,S,A,T,B,X)

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

  return (u/rho,easyeq_en(u,V0,Eta0,rho,weights)*rho/2,err)
end

# \Eta
@everywhere function easyeq_H(
  x::Float64,
  t::Float64,
  weights::Tuple{Array{Float64,1},Array{Float64,1}},
  v::Function
)
  return (x>t ? 2*t/x : 2)* integrate_legendre(y->2*pi*((x+t)*y+abs(x-t)*(1-y))*v((x+t)*y+abs(x-t)*(1-y)),0.,1.,weights)
end

# initialize V
@everywhere function easyeq_init_v(
  weights::Tuple{Array{Float64,1},Array{Float64,1}},
  v::Function
)
  order=length(weights[1])
  V=Array{Float64,1}(undef,order)
  V0=v(0.)
  for i in 1:order
    k=(1-weights[1][i])/(1+weights[1][i])
    V[i]=v(k)
  end
  return(V,V0)
end

# initialize Eta
@everywhere function easyeq_init_H(
  weights::Tuple{Array{Float64,1},Array{Float64,1}},
  v::Function
)
  order=length(weights[1])
  Eta=Array{Array{Float64,1},1}(undef,order)
  Eta0=Array{Float64,1}(undef,order)
  for i in 1:order
    k=(1-weights[1][i])/(1+weights[1][i])
    Eta[i]=Array{Float64,1}(undef,order)
    for j in 1:order
      y=(weights[1][j]+1)/2
      Eta[i][j]=easyeq_H(k,(1-y)/y,weights,v)
    end
    y=(weights[1][i]+1)/2
    Eta0[i]=easyeq_H(0.,(1-y)/y,weights,v)
  end
  return(Eta,Eta0)
end

# Xi(u)
@everywhere function easyeq_Xi(
  u::Array{Float64,1},
  order::Int64,
  S::Array{Float64,1},
  A::Array{Float64,1},
  T::Array{Float64,1},
  B::Array{Float64,1},
  X::Array{Float64,1}
)
  return u.-T./(2*(X.+1)).*dotPhi(B.*T./(X.+1).^2)
end

# compute E,S,A,T,B,X
@everywhere function easyeq_ESATBX(
  u::Array{Float64,1},
  V::Array{Float64,1},
  V0::Float64,
  Eta::Array{Array{Float64,1},1},
  Eta0::Array{Float64,1},
  weights::Tuple{Array{Float64,1},Array{Float64,1}},
  rho::Float64,
  approx::Easyeq_approx
)
  order=length(weights[1])

  # init
  S=zeros(Float64,order)
  A=zeros(Float64,order)
  T=zeros(Float64,order)
  B=zeros(Float64,order)
  X=zeros(Float64,order)

  # compute E before running the loop
  E=easyeq_en(u,V0,Eta0,rho,weights)

  for i in 1:order
    # k_i
    k=(1-weights[1][i])/(1+weights[1][i])
    # S_i
    S[i]=V[i]-1/(rho*(2*pi)^3)*integrate_legendre_sampled(y->(1-y)/y^3,Eta[i].*u,0.,1.,weights)

    # A_K,i
    A[i]=0.
    if approx.bK!=0.
      A[i]+=approx.bK*S[i]
    end
    if approx.bK!=1.
      A[i]+=(1-approx.bK)*E
    end

    # T_i
    if approx.bK==1.
      T[i]=1.
    else
      T[i]=S[i]/A[i]
    end

    # B_i
    if approx.bK==approx.bL
      B[i]=1.
    else
      B[i]=(approx.bL*S[i]+(1-approx.bL*E))/(approx.bK*S[i]+(1-approx.bK*E))
    end

    # X_i
    X[i]=k^2/(2*A[i]*rho)
  end

  return (E,S,A,T,B,X)
end

# derivative of Xi
@everywhere function easyeq_dXi(
  u::Array{Float64,1},
  Eta::Array{Array{Float64,1},1},
  Eta0::Array{Float64,1},
  weights::Tuple{Array{Float64,1},Array{Float64,1}},
  rho::Float64,
  approx::Easyeq_approx,
  E::Float64,
  S::Array{Float64,1},
  A::Array{Float64,1},
  T::Array{Float64,1},
  B::Array{Float64,1},
  X::Array{Float64,1}
)
  order=length(weights[1])

  # init
  out=zeros(Float64,order,order)

  for i in 1:order
    # k_i
    k=(1-weights[1][i])/(1+weights[1][i])

    for j in 1:order
      y=(weights[1][j]+1)/2
      dS=-1/rho*(1-y)*Eta[i][j]/(2*(2*pi)^3*y^3)*weights[2][j]
      dE=-1/rho*(1-y)*Eta0[j]/(2*(2*pi)^3*y^3)*weights[2][j]
      dU=(i==j ? 1. : 0.)

      out[i,j]=easyeq_dXi_of_dSdEdU(k,dS,dE,dU,E,S[i],A[i],T[i],B[i],X[i],rho,approx)
    end
  end
  
  return out
end

# dXi given dS, dE and dU
@everywhere function easyeq_dXi_of_dSdEdU(
  k::Float64,
  dS::Float64,
  dE::Float64,
  dU::Float64,
  E::Float64,
  S::Float64,
  A::Float64,
  T::Float64,
  B::Float64,
  X::Float64,
  rho::Float64,
  approx::Easyeq_approx
)
  # dA
  dA=0.
  if approx.bK!=0.
    dA+=approx.bK*dS
  end
  if approx.bK!=1.
    dA+=(1-approx.bK)*dE
  end

  # dT,dB
  # nothing to do if bK=bL=1
  if approx.bK!=1. || approx.bK!=approx.bL
    dB=(E*dS-S*dE)/A^2
  end
  if approx.bK==1.
    dT=0.
  else
    dT=(1-approx.bK)*dB
  end
  if approx.bK==approx.bL
    dB=0.
  else
    dB=(approx.bL*(1-approx.bK)-approx.bK*(1-approx.bL))*dB
  end

  dX=-k^2/(2*A^2*rho)*dA

  return dU-(dT-T*dX/(X+1))/(2*(X+1))*Phi(B*T/(X+1)^2)-T/(2*(X+1)^3)*(B*dT+T*dB-2*B*T*dX/(X+1))*dPhi(B*T/(X+1)^2)
end

# derivative of Xi with respect to mu
@everywhere function easyeq_dXidmu(
  u::Array{Float64,1},
  order::Int64,
  rho::Float64,
  A::Array{Float64,1},
  T::Array{Float64,1},
  B::Array{Float64,1},
  X::Array{Float64,1}
)
  # init
  out=zeros(Float64,order)
  for i in 1:order
    out[i]=T[i]/(2*rho*A[i]*(X[i]+1)^2)*Phi(B[i]*T[i]/(X[i]+1)^2)+B[i]*T[i]^2/(rho*A[i]*(X[i]+1)^4)*dPhi(B[i]*T[i]/(X[i]+1)^2)
  end

  return out
end

# energy
@everywhere function easyeq_en(
  u::Array{Float64,1},
  V0::Float64,
  Eta0::Array{Float64,1},
  rho::Float64,
  weights::Tuple{Array{Float64,1},Array{Float64,1}}
)
  return V0-1/(rho*(2*pi)^3)*integrate_legendre_sampled(y->(1-y)/y^3,Eta0.*u,0.,1.,weights)
end

# condensate fraction
@everywhere function easyeq_eta(
  u::Array{Float64,1},
  order::Int64,
  rho::Float64,
  v::Function,
  maxiter::Int64,
  tolerance::Float64,
  weights::Tuple{Array{Float64,1},Array{Float64,1}},
  approx::Easyeq_approx
)
  (V,V0)=easyeq_init_v(weights,v)
  (Eta,Eta0)=easyeq_init_H(weights,v)

  (E,S,A,T,B,X)=easyeq_ESATBX(rho*u,V,V0,Eta,Eta0,weights,rho,approx)
  du=-inv(easyeq_dXi(rho*u,Eta,Eta0,weights,rho,approx,E,S,A,T,B,X))*easyeq_dXidmu(rho*u,order,rho,A,T,B,X)

  eta=-1/(2*(2*pi)^3)*integrate_legendre_sampled(y->(1-y)/y^3,Eta0.*du,0.,1.,weights)

  return eta
end

# inverse Fourier transform
@everywhere function easyeq_u_x(
  x::Float64,
  u::Array{Float64,1},
  weights::Tuple{Array{Float64,1},Array{Float64,1}}
)
  order=length(weights[1])
  out=integrate_legendre_sampled(y->(1-y)/y^3*sin(x*(1-y)/y)/x/(2*pi^2),u,0.,1.,weights)
  return out
end


# correlation function
@everywhere function easyeq_C2(
  x::Float64,
  u::Array{Float64,1},
  windowL::Float64,
  rho::Float64,
  weights::Tuple{Array{Float64,1},Array{Float64,1}},
  Eta::Array{Array{Float64,1},1},
  Eta0::Array{Float64,1},
  approx::Easyeq_approx,
  E::Float64,
  S::Array{Float64,1},
  A::Array{Float64,1},
  T::Array{Float64,1},
  B::Array{Float64,1},
  X::Array{Float64,1}
)
  g=(r,x)->(r>0. ? sin(r*x)/(r*x)*hann(r,windowL) : hann(r,windowL))
  (dEta,dEta0)=easyeq_init_H(weights,k->g(k,x))
  du=easyeq_dudv(g,x,u,rho,weights,Eta,Eta0,dEta,dEta0,approx,E,S,A,T,B,X)

  return rho^2*(1-integrate_legendre_sampled(y->(1-y)/y^3,dEta0.*u+Eta0.*du,0.,1.,weights)/(8*pi^3))
end

# derivative of u with respect to v in direction g
@everywhere function easyeq_dudv(
  g::Function,# should be of the form g(k,x) where x is a parameter
  x::Float64,
  u::Array{Float64,1},
  rho::Float64,
  weights::Tuple{Array{Float64,1},Array{Float64,1}},
  Eta::Array{Array{Float64,1},1},
  Eta0::Array{Float64,1},
  dEta::Array{Array{Float64,1},1},
  dEta0::Array{Float64,1},
  approx::Easyeq_approx,
  E::Float64,
  S::Array{Float64,1},
  A::Array{Float64,1},
  T::Array{Float64,1},
  B::Array{Float64,1},
  X::Array{Float64,1}
)
  # initialize dV and dEta
  (dV,dV0)=easyeq_init_v(weights,k->g(k,x))

  du=-inv(easyeq_dXi(rho*u,Eta,Eta0,weights,rho,approx,E,S,A,T,B,X))*easyeq_dXidv(rho*u,dV,dV0,dEta,dEta0,weights,rho,approx,E,S,A,T,B,X)
  # rescale rho
  du=du/rho

  return du
end

# derivative of Xi with respect to potential
@everywhere function easyeq_dXidv(
  u::Array{Float64,1},
  dv::Array{Float64,1},
  dv0::Float64,
  dEta::Array{Array{Float64,1},1},
  dEta0::Array{Float64,1},
  weights::Tuple{Array{Float64,1},Array{Float64,1}},
  rho::Float64,
  approx::Easyeq_approx,
  E::Float64,
  S::Array{Float64,1},
  A::Array{Float64,1},
  T::Array{Float64,1},
  B::Array{Float64,1},
  X::Array{Float64,1}
)

  order=length(weights[1])

  # init
  out=zeros(Float64,order)

  for i in 1:order
    # k_i
    k=(1-weights[1][i])/(1+weights[1][i])

    dS=dv[i]
    dE=dv0
    for j in 1:order
      y=(weights[1][j]+1)/2
      dS+=-1/rho*(1-y)*u[j]*dEta[i][j]/(2*(2*pi)^3*y^3)*weights[2][j]
      dE+=-1/rho*(1-y)*u[j]*dEta0[j]/(2*(2*pi)^3*y^3)*weights[2][j]
    end

    out[i]=easyeq_dXi_of_dSdEdU(k,dS,dE,0.,E,S[i],A[i],T[i],B[i],X[i],rho,approx)
  end

  return out
end

# maximum of 2 point correlation function
@everywhere function easyeq_C2_max(
  u::Array{Float64,1},
  x0::Float64,
  dx::Float64,
  maxstep::Float64,
  maxiter::Int64,
  tolerance::Float64,
  windowL::Float64,
  rho::Float64,
  weights::Tuple{Array{Float64,1},Array{Float64,1}},
  V::Array{Float64,1},
  V0::Float64,
  Eta::Array{Array{Float64,1},1},
  Eta0::Array{Float64,1},
  approx::Easyeq_approx
)
  # compute some useful terms
  (E,S,A,T,B,X)=easyeq_ESATBX(rho*u,V,V0,Eta,Eta0,weights,rho,approx)

  (x,f)=newton_maximum(y->easyeq_C2(y,u,windowL,rho,weights,Eta,Eta0,approx,E,S,A,T,B,X),x0,dx,maxiter,tolerance,maxstep)

  return(x,f)
end

# Fourier transform of correlation function
@everywhere function easyeq_C2_fourier(
  q::Float64,
  u::Array{Float64,1},
  windowL::Float64,
  rho::Float64,
  weights::Tuple{Array{Float64,1},Array{Float64,1}},
  Eta::Array{Array{Float64,1},1},
  Eta0::Array{Float64,1},
  approx::Easyeq_approx,
  E::Float64,
  S::Array{Float64,1},
  A::Array{Float64,1},
  T::Array{Float64,1},
  B::Array{Float64,1},
  X::Array{Float64,1}
)
  # direction in which to differentiate u
  g=(r,x)->(r>0. ? (1.)/(2*x*r*windowL)*(gaussian(x-r,(1.)/windowL)-gaussian(x+r,(1.)/windowL)) : gaussian(x,(1.)/windowL))
  (dEta,dEta0)=easyeq_init_H(weights,k->g(k,q))
  du=easyeq_dudv(g,q,u,rho,weights,Eta,Eta0,dEta,dEta0,approx,E,S,A,T,B,X)

  return rho^2*(-integrate_legendre_sampled(y->(1-y)/y^3,dEta0.*u+Eta0.*du,0.,1.,weights)/(8*pi^3))
end


# maximum of Fourier transform of 2 point correlation function
@everywhere function easyeq_C2_fourier_max(
  u::Array{Float64,1},
  k0::Float64,
  dk::Float64,
  maxstep::Float64,
  maxiter::Int64,
  tolerance::Float64,
  windowL::Float64,
  rho::Float64,
  weights::Tuple{Array{Float64,1},Array{Float64,1}},
  V::Array{Float64,1},
  V0::Float64,
  Eta::Array{Array{Float64,1},1},
  Eta0::Array{Float64,1},
  approx::Easyeq_approx
)
  # compute some useful terms
  (E,S,A,T,B,X)=easyeq_ESATBX(rho*u,V,V0,Eta,Eta0,weights,rho,approx)

  (k,f)=newton_maximum(y->easyeq_C2_fourier(y,u,windowL,rho,weights,Eta,Eta0,approx,E,S,A,T,B,X),k0,dk,maxiter,tolerance,maxstep)

  return (k,f)
end


@everywhere function barEta(
  q::Float64,
  y::Float64,
  t::Float64
)
  return (t>=abs(y-q) && t<=y+q ? pi/y/q : 0.)
end