Update to v0.4v0.4

  feature: compute the 2-point correlation function in easyeq.

  feature: compute the Fourier transform of the 2-point correlation function
           in anyeq and easyeq.

  feature: compute the local maximum of the 2-point correlation function and
           its Fourier transform.

  feature: compute the compressibility for anyeq.

  feature: allow for linear spacing of rho's.

  feature: print the scattering length.

  change: ux and uk now return real numbers.

  fix: error in the computation of the momentum distribution: wrong
       definition of delta functions.

  fix: various minor bugs.

  optimization: assign explicit types to variables.

14 files changed, 3195 insertions, 696 deletions

diff --git a/src/anyeq.jl b/src/anyeq.jl
index 8459cb0..fee77d8 100644
--- a/src/anyeq.jl
+++ b/src/anyeq.jl
@@ -1,4 +1,4 @@
-## Copyright 2021 Ian Jauslin
+## 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.
@@ -28,24 +28,26 @@
 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)
+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,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
+  (weights,T,k,V,V0,A,Abar,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v,approx,savefile)
 
   # 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]
+  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)
 
@@ -53,72 +55,60 @@ function anyeq_energy(rho,minlrho,nlrho,taus,P,N,J,a0,v,maxiter,tolerance,approx
 end
 
 # compute energy as a function of rho
-function anyeq_energy_rho(rhos,taus,P,N,J,a0,v,maxiter,tolerance,approx,savefile)
+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,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
+  (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,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))
+  u0s=anyeq_init_medeq(rhos,minlrho_init,nlrho_init,N,J,k,a0,v,maxiter,tolerance)
 
-  ## 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
+  # spawn workers
+  work=spawn_workers(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
-      (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
+  # 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],err[j])
+    @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,taus,P,N,J,a0,v,maxiter,tolerance,approx,savefile)
+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,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
+  (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[1]],N,J,k,a0,v,maxiter,tolerance)
-  u=u0s[1]
+  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)
@@ -134,20 +124,26 @@ function anyeq_energy_rho_init_prevrho(rhos,taus,P,N,J,a0,v,maxiter,tolerance,ap
 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)
+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,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
+  (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[length(rhos)]],N,J,k,a0,v,maxiter,tolerance)
-  u=u0s[1]
+  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)
@@ -163,23 +159,26 @@ function anyeq_energy_rho_init_nextrho(rhos,taus,P,N,J,a0,v,maxiter,tolerance,ap
 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)
+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,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
+  (weights,T,k,V,V0,A,Abar,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v,approx,savefile)
+
   # 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]
+  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)
 
@@ -192,58 +191,60 @@ function anyeq_uk(minlrho,nlrho,taus,P,N,J,rho,a0,v,maxiter,tolerance,approx,sav
 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)
+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,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
+  (weights,T,k,V,V0,A,Abar,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v,approx,savefile)
 
   # 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]
+  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=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))
+    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
-# 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)
+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,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
+  (weights,T,k,V,V0,A,Abar,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v,approx,savefile)
 
   # 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]
+  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)
 
@@ -254,161 +255,571 @@ function anyeq_condensate_fraction(rho,minlrho,nlrho,taus,P,N,J,a0,v,maxiter,tol
 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)
+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)
-    append!(work[(j-1)%nw+1],j)
+    @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,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v)
-  # init Abar
-  if savefile!=""
-    Abar=anyeq_read_Abar(savefile,P,N,J)
+  (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
-    Abar=anyeq_Abar_multithread(k,weights,taus,T,P,N,J,approx)
+    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
-  u0s=anyeq_init_medeq(rhos,N,J,k,a0,v,maxiter,tolerance)
+  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=Array{Array{Float64,1}}(undef,length(rhos))
-  errs=Array{Float64,1}(undef,length(rhos))
+  (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:nw
+  @sync for p in 1:length(work)
     # for each task
     @async for j in work[p]
       count=count+1
-      if count>=nw
-	progress(count,length(rhos),10000)
+      if count>=length(work)
+        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)
+      (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
 
-  # compute eta
-  etas=Array{Float64}(undef,length(rhos))
+  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:nw
+  @sync for p in 1:length(work)
     # for each task
     @async for j in work[p]
       count=count+1
-      if count>=nw
-	progress(count,length(rhos),10000)
+      if count>=length(work)
+        progress(count,nk,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)
+      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:length(rhos)
-    @printf("% .15e % .15e % .15e\n",rhos[j],etas[j],errs[j])
+  for j in 1:nk
+    @printf("% .15e % .15e\n",kmin+(kmax-kmin)*j/nk,C2s[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)
+# 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
-    Abar=anyeq_Abar_multithread(k,weights,taus,T,P,N,J,approx)
+    @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
-  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]
+  u0s=anyeq_init_medeq(rhos,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)
+  # spawn workers
+  work=spawn_workers(length(rhos))
 
-  # compute M
-  M=anyeq_momentum(u,P,N,J,rho,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,approx)
+  # 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 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])
+  # 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
 
-# 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)
+# 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,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)
+  (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
-  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))
+  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
-  u0s=anyeq_init_medeq(rhos,N,J,k,a0,v,maxiter,tolerance)
-  u0=u0s[nlrho]
+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)
 
-  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)
+  # 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,weights,taus,T,P,N,J,approx)
+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.
-    return []
+    # empty array
+    return Array{Array{Float64,5}}(undef,0)
   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
+  # spawn workers
+  work=spawn_workers(N*J)
 
   count=0
   # for each worker
-  @sync for p in 1:nw
+  @sync for p in 1:length(work)
     # for each task
     @async for j in work[p]
       count=count+1
-      if count>=nw
+      if count>=length(work)
         progress(count,N*J,10000)
       end
       # run the task
@@ -419,14 +830,23 @@ function anyeq_Abar_multithread(k,weights,taus,T,P,N,J,approx)
   return out
 end
 
+
 # initialize computation
-@everywhere function anyeq_init(taus,P,N,J,v)
+@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}(undef,J*N)
-  k=Array{Float64}(undef,J*N)
+  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]
@@ -437,7 +857,7 @@ end
     end
   end
   # potential at 0
-  V0=v(0)
+  V0=v(0.)
 
   # initialize matrix A
   T=chebyshev_polynomials(P)
@@ -449,39 +869,62 @@ end
   Upsilon=Upsilonmat(k,v,weights_plus)
   Upsilon0=Upsilon0mat(k,v,weights_plus)
 
-  return(weights,T,k,V,V0,A,Upsilon,Upsilon0)
+  # 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,N,J,k,a0,v,maxiter,tolerance)
-  us_medeq=Array{Array{Float64,1}}(undef,length(rhos))
-  u0s=Array{Array{Float64,1}}(undef,length(rhos))
-
+@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.))
 
-  (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
+  # 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 u0s
+  # 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,weights,k,N,J)
+@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}(undef,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)
@@ -501,7 +944,27 @@ 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)
+@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
@@ -512,7 +975,7 @@ end
   # 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)
+    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)
@@ -528,8 +991,60 @@ end
 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,P,N,J,v,approx)
+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)
 
@@ -555,7 +1070,12 @@ function anyeq_save_Abar(taus,P,N,J,v,approx)
 end
 
 # read Abar
-function anyeq_read_Abar(savefile,P,N,J)
+function anyeq_read_Abar(
+  savefile::String,
+  P::Int64,
+  N::Int64,
+  J::Int64
+)
    # open file
   ff=open(savefile,"r")
   # read all lines
@@ -623,13 +1143,39 @@ end
 
 # Xi
 # takes the vector of kj's and xn's as input
-@everywhere function anyeq_Xi(U,X,Y)
+@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,rho,k,taus,v,v0,A,Abar,Upsilon,Upsilon0,weights,P,N,J,approx,S,E,II,JJ,X,Y,sL1,sK)
+@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
@@ -720,7 +1266,23 @@ end
 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)
+@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)
 
@@ -798,79 +1360,207 @@ end
   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,P,N,J,rho,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,approx)
+@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,V,V0,A,Abar,Upsilon,Upsilon0,weights,P,N,J,approx,S,E,II,JJ,X,Y,sL1,sK))*dXidmu
+  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
 
-# 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))
+# 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)
 
-  M=Array{Float64}(undef,J*N)
-  for i in 1:J*N
-    # du/dlambda
-    du=dXidu[:,i]*dXidlambda[i]*delta[i]
+  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
 
-    # compute M
-    M[i]=-avg_v_chebyshev(du,Upsilon0,k,taus,weights,N,J)/2
+# 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)
 
-  return M
-end
+  # 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
 
-# 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)
+# 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}(undef,J*N)
+  dV=Array{Float64,1}(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
+    dV[i]=g(k[i],x)
   end
-  dV0=1.
+  dV0=g(0.,x)
 
   # 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)
+  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,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)
+  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
 
-  C2=rho^2*(1-integrate_f_chebyshev(s->1.,u.*dV+V.*du,k,taus,weights,N,J))
-
-  return C2
+  return du
 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)
+# 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)
 
@@ -896,7 +1586,7 @@ end
       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
+      dJJ+=approx.gL3*(1-approx.bL3)*dE*(UU/rho).^2
     end
   end
 
@@ -947,3 +1637,220 @@ end
   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
diff --git a/src/chebyshev.jl b/src/chebyshev.jl
index 28c8f1f..af4be40 100644
--- a/src/chebyshev.jl
+++ b/src/chebyshev.jl
@@ -1,4 +1,4 @@
-## Copyright 2021 Ian Jauslin
+## 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.
@@ -13,7 +13,15 @@
 ## limitations under the License.
 
 # Chebyshev expansion
-@everywhere function chebyshev(a,taus,weights,P,N,J,nu)
+@everywhere function chebyshev(
+  a::Array{Float64,1},
+  taus::Array{Float64,1},
+  weights::Tuple{Array{Float64,1},Array{Float64,1}},
+  P::Int64,
+  N::Int64,
+  J::Int64,
+  nu::Int64
+)
   out=zeros(Float64,J*(P+1))
   for zeta in 0:J-1
     for n in 0:P
@@ -27,7 +35,15 @@
 end
 
 # evaluate function from Chebyshev expansion
-@everywhere function chebyshev_eval(Fa,x,taus,chebyshev,P,J,nu)
+@everywhere function chebyshev_eval(
+  Fa::Array{Float64,1},
+  x::Float64,
+  taus::Array{Float64,1},
+  chebyshev::Array{Polynomial,1},
+  P::Int64,
+  J::Int64,
+  nu::Int64
+)
   # change variable
   tau=(1-x)/(1+x)
 
@@ -48,7 +64,11 @@ end
 
 # convolution 
 # input the Chebyshev expansion of a and b, as well as the A matrix
-@everywhere function conv_chebyshev(Fa,Fb,A)
+@everywhere function conv_chebyshev(
+  Fa::Array{Float64,1},
+  Fb::Array{Float64,1},
+  A::Array{Array{Float64,2},1}
+)
   out=zeros(Float64,length(A))
   for i in 1:length(A)
     out[i]=dot(Fa,A[i]*Fb)
@@ -57,12 +77,27 @@ end
 end
 
 # <ab>
-@everywhere function avg_chebyshev(Fa,Fb,A0)
+@everywhere function avg_chebyshev(
+  Fa::Array{Float64,1},
+  Fb::Array{Float64,1},
+  A0::Float64
+)
   return dot(Fa,A0*Fb)
 end
 
 # 1_n * a
-@everywhere function conv_one_chebyshev(n,zetapp,Fa,A,taus,weights,P,N,J,nu1)
+@everywhere function conv_one_chebyshev(
+  n::Int64,
+  zetapp::Int64,
+  Fa::Array{Float64,1},
+  A::Array{Array{Float64,2},1},
+  taus::Array{Float64,1},
+  weights::Tuple{Array{Float64,1},Array{Float64,1}},
+  P::Int64,
+  N::Int64,
+  J::Int64,
+  nu1::Int64
+)
   out=zeros(Float64,N*J)
   for m in 1:N*J
     for l in 0:P
@@ -74,7 +109,15 @@ end
   return out
 end
 # a * v
-@everywhere function conv_v_chebyshev(a,Upsilon,k,taus,weights,N,J)
+@everywhere function conv_v_chebyshev(
+  a::Array{Float64,1},
+  Upsilon::Array{Array{Float64,1},1},
+  k::Array{Float64,1},
+  taus::Array{Float64,1},
+  weights::Tuple{Array{Float64,1},Array{Float64,1}},
+  N::Int64,
+  J::Int64
+)
   out=zeros(Float64,J*N)
   for i in 1:J*N
     for zetap in 0:J-1
@@ -86,7 +129,16 @@ end
   end
   return out
 end
-@everywhere function conv_one_v_chebyshev(n,zetap,Upsilon,k,taus,weights,N,J)
+@everywhere function conv_one_v_chebyshev(
+  n::Int64,
+  zetap::Int64,
+  Upsilon::Array{Array{Float64,1},1},
+  k::Array{Float64,1},
+  taus::Array{Float64,1},
+  weights::Tuple{Array{Float64,1},Array{Float64,1}},
+  N::Int64,
+  J::Int64
+)
   out=zeros(Float64,J*N)
   xj=weights[1][n]
   for i in 1:J*N
@@ -96,7 +148,14 @@ end
 end
 
 # <av>
-@everywhere function avg_v_chebyshev(a,Upsilon0,k,taus,weights,N,J)
+@everywhere function avg_v_chebyshev(a,
+  Upsilon0::Array{Float64,1},
+  k::Array{Float64,1},
+  taus::Array{Float64,1},
+  weights::Tuple{Array{Float64,1},Array{Float64,1}},
+  N::Int64,
+  J::Int64
+)
   out=0.
   for zetap in 0:J-1
     for j in 1:N
@@ -107,13 +166,28 @@ end
   return out
 end
 # <1_nv>
-@everywhere function avg_one_v_chebyshev(n,zetap,Upsilon0,k,taus,weights,N)
+@everywhere function avg_one_v_chebyshev(
+  n::Int64,
+  zetap::Int64,
+  Upsilon0::Array{Float64,1},
+  k::Array{Float64,1},
+  taus::Array{Float64,1},
+  weights::Tuple{Array{Float64,1},Array{Float64,1}},
+  N::Int64
+)
   xj=weights[1][n]
   return (taus[zetap+2]-taus[zetap+1])/(32*pi)*weights[2][n]*cos(pi*xj/2)*(1+k[zetap*N+n])^2*k[zetap*N+n]*Upsilon0[zetap*N+n]
 end
 
 # compute \int dq dxi u1(k-xi)u2(q)u3(xi)u4(k-q)u5(xi-q)
-@everywhere function double_conv_S_chebyshev(FU1,FU2,FU3,FU4,FU5,Abar)
+@everywhere function double_conv_S_chebyshev(
+  FU1::Array{Float64,1},
+  FU2::Array{Float64,1},
+  FU3::Array{Float64,1},
+  FU4::Array{Float64,1},
+  FU5::Array{Float64,1},
+  Abar::Array{Float64,5}
+)
   out=zeros(Float64,length(Abar))
   for i in 1:length(Abar)
     for j1 in 1:length(FU1)
@@ -133,7 +207,17 @@ end
 
 
 # compute A
-@everywhere function Amat(k,weights,taus,T,P,N,J,nua,nub)
+@everywhere function Amat(
+  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,
+  nua::Int64,
+  nub::Int64
+)
   out=Array{Array{Float64,2},1}(undef,J*N)
   for i in 1:J*N
     out[i]=zeros(Float64,J*(P+1),J*(P+1))
@@ -152,7 +236,11 @@ end
 end
 
 # compute Upsilon
-@everywhere function Upsilonmat(k,v,weights)
+@everywhere function Upsilonmat(
+  k::Array{Float64,1},
+  v::Function,
+  weights::Tuple{Array{Float64,1},Array{Float64,1}}
+)
   out=Array{Array{Float64,1},1}(undef,length(k))
   for i in 1:length(k)
     out[i]=Array{Float64,1}(undef,length(k))
@@ -162,7 +250,11 @@ end
   end
   return out
 end
-@everywhere function Upsilon0mat(k,v,weights)
+@everywhere function Upsilon0mat(
+  k::Array{Float64,1},
+  v::Function,
+  weights::Tuple{Array{Float64,1},Array{Float64,1}}
+)
   out=Array{Float64,1}(undef,length(k))
   for j in 1:length(k)
     out[j]=2*k[j]*v(k[j])
@@ -171,17 +263,36 @@ end
 end
 
 # alpha_-
-@everywhere function alpham(k,t)
+@everywhere function alpham(
+  k::Float64,
+  t::Float64
+)
   return (1-k-(1-t)/(1+t))/(1+k+(1-t)/(1+t))
 end
 # alpha_+
-@everywhere function alphap(k,t)
+@everywhere function alphap(
+  k::Float64,
+  t::Float64
+)
   return (1-abs(k-(1-t)/(1+t)))/(1+abs(k-(1-t)/(1+t)))
 end
 
 
 # compute \bar A
-@everywhere function barAmat(k,weights,taus,T,P,N,J,nu1,nu2,nu3,nu4,nu5)
+@everywhere function barAmat(
+  k::Float64,
+  weights::Tuple{Array{Float64,1},Array{Float64,1}},
+  taus::Array{Float64,1},
+  T::Array{Polynomial,1},
+  P::Int64,
+  N::Int64,
+  J::Int64,
+  nu1::Int64,
+  nu2::Int64,
+  nu3::Int64,
+  nu4::Int64,
+  nu5::Int64
+)
   out=zeros(Float64,J*(P+1),J*(P+1),J*(P+1),J*(P+1),J*(P+1))
   for zeta1 in 0:J-1
     for n1 in 0:P
@@ -211,27 +322,107 @@ end
 
   return out
 end
-@everywhere function barAmat_int1(tau,k,taus,T,weights,nu1,nu2,nu3,nu4,nu5,zeta1,zeta2,zeta3,zeta4,zeta5,n1,n2,n3,n4,n5)
+@everywhere function barAmat_int1(tau,
+  k::Float64,
+  taus::Array{Float64,1},
+  T::Array{Polynomial,1},
+  weights::Tuple{Array{Float64,1},Array{Float64,1}},
+  nu1::Int64,
+  nu2::Int64,
+  nu3::Int64,
+  nu4::Int64,
+  nu5::Int64,
+  zeta1::Int64,
+  zeta2::Int64,
+  zeta3::Int64,
+  zeta4::Int64,
+  zeta5::Int64,
+  n1::Int64,
+  n2::Int64,
+  n3::Int64,
+  n4::Int64,
+  n5::Int64
+)
   if(alpham(k,tau)<taus[zeta2+2] && alphap(k,tau)>taus[zeta2+1])
     return 2*(1-tau)/(1+tau)^(3-nu1)*T[n1+1]((2*tau-(taus[zeta1+1]+taus[zeta1+2]))/(taus[zeta1+2]-taus[zeta1+1]))*integrate_legendre(sigma->barAmat_int2(tau,sigma,k,taus,T,weights,nu2,nu3,nu4,nu5,zeta2,zeta3,zeta4,zeta5,n2,n3,n4,n5),max(taus[zeta2+1],alpham(k,tau)),min(taus[zeta2+2],alphap(k,tau)),weights)
   else
     return 0.
   end
 end
-@everywhere function barAmat_int2(tau,sigma,k,taus,T,weights,nu2,nu3,nu4,nu5,zeta2,zeta3,zeta4,zeta5,n2,n3,n4,n5)
+@everywhere function barAmat_int2(tau,
+  sigma::Float64,
+  k::Float64,
+  taus::Array{Float64,1},
+  T::Array{Polynomial,1},
+  weights::Tuple{Array{Float64,1},Array{Float64,1}},
+  nu2::Int64,
+  nu3::Int64,
+  nu4::Int64,
+  nu5::Int64,
+  zeta2::Int64,
+  zeta3::Int64,
+  zeta4::Int64,
+  zeta5::Int64,
+  n2::Int64,
+  n3::Int64,
+  n4::Int64,
+  n5::Int64
+)
   return 2*(1-sigma)/(1+sigma)^(3-nu2)*T[n2+1]((2*sigma-(taus[zeta2+1]+taus[zeta2+2]))/(taus[zeta2+2]-taus[zeta2+1]))*integrate_legendre(taup->barAmat_int3(tau,sigma,taup,k,taus,T,weights,nu3,nu4,nu5,zeta3,zeta4,zeta5,n3,n4,n5),taus[zeta3+1],taus[zeta3+2],weights)
 end
-@everywhere function barAmat_int3(tau,sigma,taup,k,taus,T,weights,nu3,nu4,nu5,zeta3,zeta4,zeta5,n3,n4,n5)
+@everywhere function barAmat_int3(tau,
+  sigma::Float64,
+  taup::Float64,
+  k::Float64,
+  taus::Array{Float64,1},
+  T::Array{Polynomial,1},
+  weights::Tuple{Array{Float64,1},Array{Float64,1}},
+  nu3::Int64,
+  nu4::Int64,
+  nu5::Int64,
+  zeta3::Int64,
+  zeta4::Int64,
+  zeta5::Int64,
+  n3::Int64,
+  n4::Int64,
+  n5::Int64
+)
   if(alpham(k,taup)<taus[zeta4+2] && alphap(k,taup)>taus[zeta4+1])
     return 2*(1-taup)/(1+taup)^(3-nu3)*T[n3+1]((2*taup-(taus[zeta3+1]+taus[zeta3+2]))/(taus[zeta3+2]-taus[zeta3+1]))*integrate_legendre(sigmap->barAmat_int4(tau,sigma,taup,sigmap,k,taus,T,weights,nu4,nu5,zeta4,zeta5,n4,n5),max(taus[zeta4+1],alpham(k,taup)),min(taus[zeta4+2],alphap(k,taup)),weights)
   else
     return 0.
   end
 end
-@everywhere function barAmat_int4(tau,sigma,taup,sigmap,k,taus,T,weights,nu4,nu5,zeta4,zeta5,n4,n5)
+@everywhere function barAmat_int4(tau,
+  sigma::Float64,
+  taup::Float64,
+  sigmap::Float64,
+  k::Float64,
+  taus::Array{Float64,1},
+  T::Array{Polynomial,1},
+  weights::Tuple{Array{Float64,1},Array{Float64,1}},
+  nu4::Int64,
+  nu5::Int64,
+  zeta4::Int64,
+  zeta5::Int64,
+  n4::Int64,
+  n5::Int64
+)
   return 2*(1-sigmap)/(1+sigmap)^(3-nu4)*T[n4+1]((2*sigma-(taus[zeta4+1]+taus[zeta4+2]))/(taus[zeta4+2]-taus[zeta4+1]))*integrate_legendre(theta->barAmat_int5(tau,sigma,taup,sigmap,theta,k,taus,T,weights,nu5,zeta5,n5),0.,2*pi,weights)
 end
-@everywhere function barAmat_int5(tau,sigma,taup,sigmap,theta,k,taus,T,weights,nu5,zeta5,n5)
+@everywhere function barAmat_int5(tau,
+  sigma::Float64,
+  taup::Float64,
+  sigmap::Float64,
+  theta::Float64,
+  k::Float64,
+  taus::Array{Float64,1},
+  T::Array{Polynomial,1},
+  weights::Tuple{Array{Float64,1},Array{Float64,1}},
+  nu5::Int64,
+  zeta5::Int64,
+  n5::Int64
+)
   R=barAmat_R((1-sigma)/(1+sigma),(1-tau)/(1+tau),(1-sigmap)/(1+sigmap),(1-taup)/(1+taup),theta,k)
   if((1-R)/(1+R)<taus[zeta5+2] && (1-R)/(1+R)>taus[zeta5+1])
     return (2/(2+R))^nu5*T[n5+1]((2*(1-R)/(1+R)-(taus[zeta5+1]+taus[zeta5+2]))/(taus[zeta5+2]-taus[zeta5+1]))
@@ -240,13 +431,22 @@ end
   end
 end
 # R(s,t,s',t,theta,k)
-@everywhere function barAmat_R(s,t,sp,tp,theta,k)
-  return sqrt(k^2*(s^2+t^2+sp^2+tp^2)-k^4-(s^2-t^2)*(sp^2-tp^2)-sqrt((4*k^2*s^2-(k^2+s^2-t^2)^2)*(4*k^2*sp^2-(k^2+sp^2-tp^2)^2))*cos(theta))/(sqrt(2)*k)
+@everywhere function barAmat_R(
+  s::Float64,
+  t::Float64,
+  sp::Float64,
+  tp::Float64,
+  theta::Float64,
+  k::Float64
+)
+  return sqrt(k^2*(s^2+t^2+sp^2+tp^2)-k^4-(s^2-t^2)*(sp^2-tp^2)-sqrt((4*k^2*s^2-(k^2+s^2-t^2)^2)*(4*k^2*sp^2-(k^2+sp^2-tp^2)^2))*cos(theta))/(sqrt(2.)*k)
 end
 
 # compute Chebyshev polynomials
-@everywhere function chebyshev_polynomials(P)
-  T=Array{Polynomial}(undef,P+1)
+@everywhere function chebyshev_polynomials(
+  P::Int64
+)
+  T=Array{Polynomial,1}(undef,P+1)
   T[1]=Polynomial([1])
   T[2]=Polynomial([0,1])
   for n in 1:P-1
@@ -258,7 +458,15 @@ end
 end
 
 # compute \int f*u dk/(2*pi)^3
-@everywhere function integrate_f_chebyshev(f,u,k,taus,weights,N,J)
+@everywhere function integrate_f_chebyshev(
+  f::Function,
+  u::Array{Float64,1},
+  k::Array{Float64,1},
+  taus::Array{Float64,1},
+  weights::Tuple{Array{Float64,1},Array{Float64,1}},
+  N::Int64,
+  J::Int64
+)
   out=0.
   for zeta in 0:J-1
     for i in 1:N
@@ -268,7 +476,15 @@ end
   return out
 end
 
-@everywhere function inverse_fourier_chebyshev(u,x,k,taus,weights,N,J)
+@everywhere function inverse_fourier_chebyshev(
+  u::Array{Float64,1},
+  x::Float64,
+  k::Array{Float64,1},
+  taus::Array{Float64,1},
+  weights::Tuple{Array{Float64,1},Array{Float64,1}},
+  N::Int64,
+  J::Int64
+)
   out=0.
   for zeta in 0:J-1
     for j in 1:N
@@ -277,3 +493,54 @@ end
   end
   return out
 end
+
+# compute B (for the computation of the fourier transform of the two-point correlation)
+@everywhere function Bmat(
+  q::Float64,
+  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,
+  nu::Int64
+)
+  out=Array{Array{Float64,1},1}(undef,J*N)
+  for i in 1:J*N
+    out[i]=zeros(Float64,J*(P+1))
+    for zeta in 0:J-1
+      for n in 0:P
+	out[i][zeta*(P+1)+n+1]=1/(8*pi^3*k[i]*q)*(betam(k[i],q)>taus[zeta+2] || betap(k[i],q)<taus[zeta+1] ? 0. : integrate_legendre(sigma->(1-sigma)/(1+sigma)^(3-nu)*T[n+1]((2*sigma-(taus[zeta+1]+taus[zeta+2]))/(taus[zeta+2]-taus[zeta+1])),max(taus[zeta+1],betam(k[i],q)),min(taus[zeta+2],betap(k[i],q)),weights))
+      end
+    end
+  end
+
+  return out
+end
+# beta_-
+@everywhere function betam(
+  k::Float64,
+  q::Float64
+)
+  return (1-k-q)/(1+k+q)
+end
+# beta_+
+@everywhere function betap(
+  k::Float64,
+  q::Float64
+)
+  return (1-abs(k-q))/(1+abs(k-q))
+end
+
+# mathfrak S (for the computation of the fourier transform of the two-point correlation)
+@everywhere function chebyshev_frakS(
+  Ff::Array{Float64,1},
+  B::Array{Array{Float64,1},1}
+)
+  out=zeros(Float64,length(B))
+  for i in 1:length(B)
+    out[i]=dot(Ff,B[i])
+  end
+  return out
+end
diff --git a/src/easyeq.jl b/src/easyeq.jl
index 0bde3ab..2dbbd1a 100644
--- a/src/easyeq.jl
+++ b/src/easyeq.jl
@@ -1,4 +1,4 @@
-## Copyright 2021 Ian Jauslin
+## 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.
@@ -19,47 +19,65 @@
 end
 
 # compute energy
-function easyeq_energy(minlrho_init,nlrho_init,order,rho,a0,v,maxiter,tolerance,approx)
+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 initial guess from previous rho
-  (u,E,err)=easyeq_hatu(easyeq_init_u(a0,order,weights),order,(10.)^minlrho_init,v,maxiter,tolerance,weights,approx)
-  for j in 2:nlrho_init
-    rho_tmp=10^(minlrho_init+(log10(rho)-minlrho_init)*(j-1)/(nlrho_init-1))
-    (u,E,err)=easyeq_hatu(u,order,rho_tmp,v,maxiter,tolerance,weights,approx)
-  end
+  # 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,order,a0,v,maxiter,tolerance,approx)
+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)
-  # init u
-  u=easyeq_init_u(a0,order,weights)
+  # 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 u (init newton with previously computed u)
-    (u,E,err)=easyeq_hatu(u,order,rhos[j],v,maxiter,tolerance,weights,approx)
-
-    @printf("% .15e % .15e % .15e\n",rhos[j],real(E),err)
-
+    @printf("% .15e % .15e % .15e\n",rhos[j],real(es[j]),errs[j])
   end
 end
 
 # compute u(k)
-function easyeq_uk(minlrho_init,nlrho_init,order,rho,a0,v,maxiter,tolerance,approx)
+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 initial guess from previous rho
-  (u,E,err)=easyeq_hatu(easyeq_init_u(a0,order,weights),order,(10.)^minlrho_init,v,maxiter,tolerance,weights,approx)
-  for j in 2:nlrho_init
-    rho_tmp=10^(minlrho_init+(log10(rho)-minlrho_init)*(j-1)/(nlrho_init-1))
-    (u,E,err)=easyeq_hatu(u,order,rho_tmp,v,maxiter,tolerance,weights,approx)
-  end
+  # 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])
@@ -68,15 +86,24 @@ function easyeq_uk(minlrho_init,nlrho_init,order,rho,a0,v,maxiter,tolerance,appr
 end
 
 # compute u(x)
-function easyeq_ux(minlrho_init,nlrho_init,order,rho,a0,v,maxiter,tolerance,xmin,xmax,nx,approx)
+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 initial guess from previous rho
-  (u,E,err)=easyeq_hatu(easyeq_init_u(a0,order,weights),order,(10.)^minlrho_init,v,maxiter,tolerance,weights,approx)
-  for j in 2:nlrho_init
-    rho_tmp=10^(minlrho_init+(log10(rho)-minlrho_init)*(j-1)/(nlrho_init-1))
-    (u,E,err)=easyeq_hatu(u,order,rho_tmp,v,maxiter,tolerance,weights,approx)
-  end
+  # 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
@@ -85,15 +112,24 @@ function easyeq_ux(minlrho_init,nlrho_init,order,rho,a0,v,maxiter,tolerance,xmin
 end
 
 # compute 2u(x)-rho u*u(x)
-function easyeq_uux(minlrho_init,nlrho_init,order,rho,a0,v,maxiter,tolerance,xmin,xmax,nx,approx)
+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 initial guess from previous rho
-  (u,E,err)=easyeq_hatu(easyeq_init_u(a0,order,weights),order,(10.)^minlrho_init,v,maxiter,tolerance,weights,approx)
-  for j in 2:nlrho_init
-    rho_tmp=10^(minlrho_init+(log10(rho)-minlrho_init)*(j-1)/(nlrho_init-1))
-    (u,E,err)=easyeq_hatu(u,order,rho_tmp,v,maxiter,tolerance,weights,approx)
-  end
+  # 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
@@ -102,16 +138,22 @@ function easyeq_uux(minlrho_init,nlrho_init,order,rho,a0,v,maxiter,tolerance,xmi
 end
 
 # condensate fraction
-function easyeq_condensate_fraction(minlrho_init,nlrho_init,order,rho,a0,v,maxiter,tolerance,approx)
+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 initial guess from previous rho
-  (u,E,err)=easyeq_hatu(easyeq_init_u(a0,order,weights),order,(10.)^minlrho_init,v,maxiter,tolerance,weights,approx)
-  for j in 2:nlrho_init
-    rho_tmp=10^(minlrho_init+(log10(rho)-minlrho_init)*(j-1)/(nlrho_init-1))
-    (u,E,err)=easyeq_hatu(u,order,rho_tmp,v,maxiter,tolerance,weights,approx)
-  end
+  # 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)
@@ -121,25 +163,350 @@ function easyeq_condensate_fraction(minlrho_init,nlrho_init,order,rho,a0,v,maxit
 end
 
 # condensate fraction as a function of rho
-function easyeq_condensate_fraction_rho(rhos,order,a0,v,maxiter,tolerance,approx)
+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)
-  # init u
-  u=easyeq_init_u(a0,order,weights)
+  # 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 u (init newton with previously computed u)
-    (u,E,err)=easyeq_hatu(u,order,rhos[j],v,maxiter,tolerance,weights,approx)
-
     # compute eta
-    eta=easyeq_eta(u,order,rhos[j],v,maxiter,tolerance,weights,approx)
+    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)
 
-    @printf("% .15e % .15e % .15e\n",rhos[j],eta,err)
+  # 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)
 
-# initialize u
-@everywhere function easyeq_init_u(a0,order,weights)
+  # 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
@@ -150,8 +517,100 @@ 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,order,rho,v,maxiter,tolerance,weights,approx)
+@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)
@@ -162,7 +621,8 @@ end
   # iterate
   err=Inf
   for i in 1:maxiter-1
-    new=u-inv(easyeq_dXi(u,V,V0,Eta,Eta0,weights,rho,approx))*easyeq_Xi(u,V,V0,Eta,Eta0,weights,rho,approx)
+    (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)
@@ -176,15 +636,23 @@ end
 end
 
 # \Eta
-@everywhere function easyeq_H(x,t,weights,v)
-  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)
+@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,v)
+@everywhere function easyeq_init_v(
+  weights::Tuple{Array{Float64,1},Array{Float64,1}},
+  v::Function
+)
   order=length(weights[1])
-  V=Array{Float64}(undef,order)
-  V0=v(0)
+  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)
@@ -193,29 +661,58 @@ end
 end
 
 # initialize Eta
-@everywhere function easyeq_init_H(weights,v)
+@everywhere function easyeq_init_H(
+  weights::Tuple{Array{Float64,1},Array{Float64,1}},
+  v::Function
+)
   order=length(weights[1])
-  Eta=Array{Array{Float64}}(undef,order)
-  Eta0=Array{Float64}(undef,order)
+  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}(undef,order)
+    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)
+    Eta0[i]=easyeq_H(0.,(1-y)/y,weights,v)
   end
   return(Eta,Eta0)
 end
 
 # Xi(u)
-@everywhere function easyeq_Xi(u,V,V0,Eta,Eta0,weights,rho,approx)
+@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
-  out=zeros(Float64,order)
+  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)
@@ -224,188 +721,359 @@ end
     # k_i
     k=(1-weights[1][i])/(1+weights[1][i])
     # S_i
-    S=V[i]-1/(rho*(2*pi)^3)*integrate_legendre_sampled(y->(1-y)/y^3,Eta[i].*u,0,1,weights)
+    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=0.
+    A[i]=0.
     if approx.bK!=0.
-      A+=approx.bK*S
+      A[i]+=approx.bK*S[i]
     end
     if approx.bK!=1.
-      A+=(1-approx.bK)*E
+      A[i]+=(1-approx.bK)*E
     end
 
-    # T
+    # T_i
     if approx.bK==1.
-      T=1.
+      T[i]=1.
     else
-      T=S/A
+      T[i]=S[i]/A[i]
     end
 
-    # B
+    # B_i
     if approx.bK==approx.bL
-      B=1.
+      B[i]=1.
     else
-      B=(approx.bL*S+(1-approx.bL*E))/(approx.bK*S+(1-approx.bK*E))
+      B[i]=(approx.bL*S[i]+(1-approx.bL*E))/(approx.bK*S[i]+(1-approx.bK*E))
     end
 
     # X_i
-    X=k^2/(2*A*rho)
-
-    # U_i
-    out[i]=u[i]-T/(2*(X+1))*Phi(B*T/(X+1)^2)
+    X[i]=k^2/(2*A[i]*rho)
   end
 
-  return out
+  return (E,S,A,T,B,X)
 end
 
 # derivative of Xi
-@everywhere function easyeq_dXi(u,V,V0,Eta,Eta0,weights,rho,approx)
+@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)
 
-  # compute E before 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=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=0.
-    if approx.bK!=0.
-      A+=approx.bK*S
-    end
-    if approx.bK!=1.
-      A+=(1-approx.bK)*E
-    end
-
-    # T
-    if approx.bK==1.
-      T=1.
-    else
-      T=S/A
-    end
-
-    # B
-    if approx.bK==approx.bL
-      B=1.
-    else
-      B=(approx.bL*S+(1-approx.bL*E))/(approx.bK*S+(1-approx.bK*E))
-    end
-
-    # X_i
-    X=k^2/(2*A*rho)
 
     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.)
 
-      # dA
-      dA=0.
-      if approx.bK!=0.
-	dA+=approx.bK*dS
-      end
-      if approx.bK!=1.
-	dA+=(1-approx.bK)*dE
-      end
+      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
 
-      # dT
-      if approx.bK==1.
-	dT=0.
-      else
-	dT=(1-approx.bK)*(E*dS-S*dE)/A^2
-      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
 
-      # dB
-      if approx.bK==approx.bL
-	dB=0.
-      else
-	dB=(approx.bL*(1-approx.bK)-approx.bK*(1-approx.bL))*(E*dS-S*dE)/(approx.bK*S+(1-approx.bK*E))^2
-      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
+  dX=-k^2/(2*A^2*rho)*dA
 
-      out[i,j]=(i==j ? 1 : 0)-(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
+  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
 
-# derivative of Xi with respect to mu
-@everywhere function easyeq_dXidmu(u,V,V0,Eta,Eta0,weights,rho,approx)
+# 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)
 
-  # 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=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=0.
-    if approx.bK!=0.
-      A+=approx.bK*S
-    end
-    if approx.bK!=1.
-      A+=(1-approx.bK)*E
+    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
 
-    # T
-    if approx.bK==1.
-      T=1.
-    else
-      T=S/A
-    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
 
-    # B
-    if approx.bK==approx.bL
-      B=1.
-    else
-      B=(approx.bL*S+(1-approx.bL*E))/(approx.bK*S+(1-approx.bK*E))
-    end
+  return out
+end
 
-    # X_i
-    X=k^2/(2*A*rho)
+# 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)
 
-    out[i]=T/(2*rho*A*(X+1)^2)*Phi(B*T/(X+1)^2)+B*T^2/(rho*A*(X+1)^4)*dPhi(B*T/(X+1)^2)
-  end
+  (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 out
+  return(x,f)
 end
 
-# energy
-@everywhere function easyeq_en(u,V0,Eta0,rho,weights)
-  return V0-1/(rho*(2*pi)^3)*integrate_legendre_sampled(y->(1-y)/y^3,Eta0.*u,0,1,weights)
+# 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
 
-# condensate fraction
-@everywhere function easyeq_eta(u,order,rho,v,maxiter,tolerance,weights,approx)
-  (V,V0)=easyeq_init_v(weights,v)
-  (Eta,Eta0)=easyeq_init_H(weights,v)
 
-  du=-inv(easyeq_dXi(rho*u,V,V0,Eta,Eta0,weights,rho,approx))*easyeq_dXidmu(rho*u,V,V0,Eta,Eta0,weights,rho,approx)
+# 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)
 
-  eta=-1/(2*(2*pi)^3)*integrate_legendre_sampled(y->(1-y)/y^3,Eta0.*du,0,1,weights)
+  (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 eta
+  return (k,f)
 end
 
-# inverse Fourier transform
-@everywhere function easyeq_u_x(x,u,weights)
-  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
+
+@everywhere function barEta(
+  q::Float64,
+  y::Float64,
+  t::Float64
+)
+  return (t>=abs(y-q) && t<=y+q ? pi/y/q : 0.)
 end
diff --git a/src/integration.jl b/src/integration.jl
index 9be4641..223d6cc 100644
--- a/src/integration.jl
+++ b/src/integration.jl
@@ -1,4 +1,4 @@
-## Copyright 2021 Ian Jauslin
+## 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.
@@ -13,7 +13,12 @@
 ## limitations under the License.
 
 # approximate \int_a^b f using Gauss-Legendre quadratures
-@everywhere function integrate_legendre(f,a,b,weights)
+@everywhere function integrate_legendre(
+  f::Function,
+  a::Float64,
+  b::Float64,
+  weights::Tuple{Array{Float64,1},Array{Float64,1}}
+)
   out=0
   for i in 1:length(weights[1])
     out+=(b-a)/2*weights[2][i]*f((b-a)/2*weights[1][i]+(b+a)/2)
@@ -21,7 +26,13 @@
   return out
 end
 # \int f*g where g is sampled at the Legendre nodes
-@everywhere function integrate_legendre_sampled(f,g,a,b,weights)
+@everywhere function integrate_legendre_sampled(
+  f::Function,
+  g::Array{Float64,1},
+  a::Float64,
+  b::Float64,
+  weights::Tuple{Array{Float64,1},Array{Float64,1}}
+)
   out=0
   for i in 1:length(weights[1])
     out+=(b-a)/2*weights[2][i]*f((b-a)/2*weights[1][i]+(b+a)/2)*g[i]
@@ -31,7 +42,12 @@ end
 
 
 # approximate \int_a^b f/sqrt((b-x)(x-a)) using Gauss-Chebyshev quadratures
-@everywhere function integrate_chebyshev(f,a,b,N)
+@everywhere function integrate_chebyshev(
+  f::Function,
+  a::Float64,
+  b::Float64,
+  N::Int64
+)
   out=0
   for i in 1:N
     out=out+pi/N*f((b-a)/2*cos((2*i-1)/(2*N)*pi)+(b+a)/2)
@@ -40,7 +56,11 @@ end
 end
 
 # approximate \int_0^\infty dr f(r)*exp(-a*r) using Gauss-Chebyshev quadratures
-@everywhere function integrate_laguerre(f,a,weights_gL)
+@everywhere function integrate_laguerre(
+  f::Function,
+  a::Float64,
+  weights_gL::Tuple{Array{Float64,1},Array{Float64,1}}
+)
   out=0.
   for i in 1:length(weights_gL[1])
     out+=1/a*f(weights_gL[1][i]/a)*weights_gL[2][i]
@@ -49,10 +69,28 @@ end
 end
 
 # Hann window
-@everywhere function hann(x,L)
+@everywhere function hann(
+  x::Float64,
+  L::Float64
+)
   if abs(x)<L/2
     return cos(pi*x/L)^2
   else
     return 0.
   end
 end
+# Fourier transform (in 3d)
+@everywhere function hann_fourier(
+  k::Float64,
+  L::Float64
+)
+  return L^2*4*pi^3/k*(((k*L)^3-4*k*L*pi^2)*cos(k*L/2)-2*(3*(k*L)^2-4*pi^2)*sin(k*L/2))/((k*L)^3-4*k*L*pi^2)^2
+end
+
+# normalized Gaussian (in 3d)
+@everywhere function gaussian(
+  k::Float64,
+  L::Float64
+)
+  return exp(-k^2/(2*L))/sqrt(8*pi^3*L^3)
+end
diff --git a/src/interpolation.jl b/src/interpolation.jl
index fa3bcdb..066bc20 100644
--- a/src/interpolation.jl
+++ b/src/interpolation.jl
@@ -1,4 +1,4 @@
-## Copyright 2021 Ian Jauslin
+## 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.
@@ -14,7 +14,11 @@
 
 # linear interpolation: given vectors x,y, compute a linear interpolation for y(x0)
 # assume x is ordered
-@everywhere function linear_interpolation(x0,x,y)
+@everywhere function linear_interpolation(
+  x0::Float64,
+  x::Array{Float64,1},
+  y::Array{Float64,1}
+)
   # if x0 is beyond all x's, then return the corresponding boundary value.
   if x0>x[length(x)]
     return y[length(y)]
@@ -28,7 +32,12 @@
   # interpolate
   return y[i]+(y[i+1]-y[i])*(x0-x[i])/(x[i+1]-x[i])
 end
-@everywhere function bracket(x0,x,a,b)
+@everywhere function bracket(
+  x0::Float64,
+  x::Array{Float64,1},
+  a::Int64,
+  b::Int64
+)
   i=floor(Int64,(a+b)/2)
   if x0<x[i]
     return bracket(x0,x,a,i)
@@ -41,15 +50,18 @@ end
 
 
 # polynomial interpolation of a family of points
-@everywhere function poly_interpolation(x,y)
+@everywhere function poly_interpolation(
+  x::Array{Float64,1},
+  y::Array{Float64,1}
+)
   # init for recursion
-  rec=Array{Polynomial{Float64}}(undef,length(x))
+  rec=Array{Polynomial{Float64},1}(undef,length(x))
   for i in 1:length(x)
     rec[i]=Polynomial([1.])
   end
 
   # compute \prod (x-x_i)
-  poly_interpolation_rec(rec,x,1,length(x))
+  poly_interpolation_rec(rec,x,1.,length(x))
 
   # sum terms together
   out=0.
@@ -60,7 +72,12 @@ end
   return out
 end
 # recursive helper function
-@everywhere function poly_interpolation_rec(out,x,a,b)
+@everywhere function poly_interpolation_rec(
+  out::Array{Float64,1},
+  x::Array{Float64,1},
+  a::Float64,
+  b::Float64
+)
   if a==b
     return
   end
@@ -91,7 +108,10 @@ end
   return
 end
 ## the following does the same, but has complexity N^2, the function above has N*log(N)
-#@everywhere function poly_interpolation(x,y)
+#@everywhere function poly_interpolation(
+#  x::Array{Float64,1},
+#  y::Array{Float64,1}
+#)
 #  out=Polynomial([0.])
 #  for i in 1:length(x)
 #    prod=Polynomial([1.])
diff --git a/src/main.jl b/src/main.jl
index 382fe6b..28fc2be 100644
--- a/src/main.jl
+++ b/src/main.jl
@@ -1,4 +1,4 @@
-## Copyright 2021 Ian Jauslin
+## 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.
@@ -22,6 +22,8 @@ include("chebyshev.jl")
 include("integration.jl")
 include("interpolation.jl")
 include("tools.jl")
+include("multithread.jl")
+include("optimization.jl")
 include("potentials.jl")
 include("print.jl")
 include("easyeq.jl")
@@ -36,9 +38,6 @@ function main()
   rho=1e-6
   e=1e-4
 
-  # incrementally initialize Newton algorithm
-  nlrho_init=1
-
   # potential
   v=k->v_exp(k,1.)
   a0=a0_exp(1.)
@@ -49,19 +48,29 @@ function main()
   v_param_e=1.
 
   # plot range when plotting in rho
-  minlrho=-6
-  maxlrho=2
+  # linear
+  minrho=1e-6
+  maxrho=1e2
+  nrho=0
+  # logarithmic
+  minlrho=-6.
+  maxlrho=2.
   nlrho=100
-  rhos=Array{Float64}(undef,0)
+  # list
+  rhos=Array{Float64,1}(undef,0)
   # plot range when plotting in e
-  minle=-6
-  maxle=2
+  minle=-6.
+  maxle=2.
   nle=100
-  es=Array{Float64}(undef,0)
+  es=Array{Float64,1}(undef,0)
   # plot range when plotting in x
-  xmin=0
-  xmax=100
+  xmin=0.
+  xmax=100.
   nx=100
+  # plot range when plotting in k
+  kmin=0.
+  kmax=10.
+  nk=100
 
   # cutoffs
   tolerance=1e-11
@@ -77,8 +86,8 @@ function main()
   J=10
 
   # starting rho from which to incrementally initialize Newton algorithm
-  # default must be set after reading rho, if not set explicitly
-  minlrho_init=nothing
+  minlrho_init=-6.
+  nlrho_init=0
 
   # Hann window for Fourier transforms
   windowL=1e3
@@ -98,6 +107,18 @@ function main()
   anyeq_compleq_approx=Anyeq_approx(1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.)
   anyeq_approx=anyeq_bigeq_approx
 
+  # numerical approximations of derivatives
+  dx=1e-7
+  dk=1e-7
+
+  # initial guess for 2pt_max
+  x0=1.
+  k0=1.
+  # maximum step in 2pt_max
+  maxstep=Inf
+  # tolerance for max search
+  tolerance_max=Inf
+
   # read cli arguments
   (params,potential,method,savefile,command)=read_args(ARGS)
 
@@ -109,8 +130,8 @@ function main()
 	print(stderr,"error: could not read parameter '",param,"'.\n")
 	exit(-1)
       end
-      lhs=terms[1]
-      rhs=terms[2]
+      lhs=string(terms[1])
+      rhs=string(terms[2])
       if lhs=="rho"
 	rho=parse(Float64,rhs)
       elseif lhs=="minlrho_init"
@@ -133,6 +154,12 @@ function main()
 	maxlrho=parse(Float64,rhs)
       elseif lhs=="nlrho"
 	nlrho=parse(Int64,rhs)
+      elseif lhs=="minrho"
+	minrho=parse(Float64,rhs)
+      elseif lhs=="maxrho"
+	maxrho=parse(Float64,rhs)
+      elseif lhs=="nrho"
+	nrho=parse(Int64,rhs)
       elseif lhs=="es"
 	es=parse_list(rhs)
       elseif lhs=="minle"
@@ -147,6 +174,12 @@ function main()
 	xmax=parse(Float64,rhs)
       elseif lhs=="nx"
 	nx=parse(Int64,rhs)
+      elseif lhs=="kmin"
+	kmin=parse(Float64,rhs)
+      elseif lhs=="kmax"
+	kmax=parse(Float64,rhs)
+      elseif lhs=="nk"
+	nk=parse(Int64,rhs)
       elseif lhs=="P"
 	P=parse(Int64,rhs)
       elseif lhs=="N"
@@ -155,6 +188,18 @@ function main()
 	J=parse(Int64,rhs)
       elseif lhs=="window_L"
 	windowL=parse(Float64,rhs)
+      elseif lhs=="dx"
+	dx=parse(Float64,rhs)
+      elseif lhs=="x0"
+	x0=parse(Float64,rhs)
+      elseif lhs=="dk"
+	dk=parse(Float64,rhs)
+      elseif lhs=="k0"
+	k0=parse(Float64,rhs)
+      elseif lhs=="maxstep"
+	maxstep=parse(Float64,rhs)
+      elseif lhs=="tolerance_max"
+	tolerance_max=parse(Float64,rhs)
       elseif lhs=="aK"
 	anyeq_approx.aK=parse(Float64,rhs)
       elseif lhs=="bK"
@@ -242,37 +287,48 @@ function main()
   ## set parameters
   # rhos
   if length(rhos)==0
-    rhos=Array{Float64}(undef,nlrho)
-    for j in 0:nlrho-1
-      rhos[j+1]=(nlrho==1 ? 10^minlrho : 10^(minlrho+(maxlrho-minlrho)/(nlrho-1)*j))
+    # linear only if nrho is specified
+    if nrho>0
+      rhos=Array{Float64,1}(undef,nrho)
+      for j in 0:nrho-1
+	rhos[j+1]=(nrho==1 ? minrho : minrho+(maxrho-minrho)/(nrho-1)*j)
+      end
+    else
+      rhos=Array{Float64,1}(undef,nlrho)
+      for j in 0:nlrho-1
+	rhos[j+1]=(nlrho==1 ? 10^minlrho : 10^(minlrho+(maxlrho-minlrho)/(nlrho-1)*j))
+      end
     end
   end
+
   # es
   if length(es)==0
-    es=Array{Float64}(undef,nle)
+    es=Array{Float64,1}(undef,nle)
     for j in 0:nle-1
       es[j+1]=(nle==1 ? 10^minle : 10^(minle+(maxle-minle)/(nle-1)*j))
     end
   end
 
-  # default minlrho_init
-  if (minlrho_init==nothing)
-    minlrho_init=log10(rho)
-  end
-
   # splines
-  taus=Array{Float64}(undef,J+1)
+  taus=Array{Float64,1}(undef,J+1)
   for j in 0:J
     taus[j+1]=-1+2*j/J
   end
 
+  # tolerance_max
+  if tolerance_max==Inf
+    tolerance_max=tolerance
+  end
+
   ## run command
-  if method=="easyeq"
+  if command=="scattering_length"
+    @printf("% .15e\n",a0)
+  elseif method=="easyeq"
     if command=="energy"
       easyeq_energy(minlrho_init,nlrho_init,order,rho,a0,v,maxiter,tolerance,easyeq_approx)
     # e(rho)
     elseif command=="energy_rho"
-      easyeq_energy_rho(rhos,order,a0,v,maxiter,tolerance,easyeq_approx)
+      easyeq_energy_rho(rhos,minlrho_init,nlrho_init,order,a0,v,maxiter,tolerance,easyeq_approx)
     # u(k)
     elseif command=="uk"
       easyeq_uk(minlrho_init,nlrho_init,order,rho,a0,v,maxiter,tolerance,easyeq_approx)
@@ -285,7 +341,28 @@ function main()
     elseif command=="condensate_fraction"
       easyeq_condensate_fraction(minlrho_init,nlrho_init,order,rho,a0,v,maxiter,tolerance,easyeq_approx)
     elseif command=="condensate_fraction_rho"
-      easyeq_condensate_fraction_rho(rhos,order,a0,v,maxiter,tolerance,easyeq_approx)
+      easyeq_condensate_fraction_rho(rhos,minlrho_init,nlrho_init,order,a0,v,maxiter,tolerance,easyeq_approx)
+    # 2pt correlation
+    elseif command=="2pt"
+      easyeq_2pt(xmin,xmax,nx,minlrho_init,nlrho_init,order,windowL,rho,a0,v,maxiter,tolerance,easyeq_approx)
+    # max of 2pt correlation
+    elseif command=="2pt_max"
+      easyeq_2pt_max(dx,x0,minlrho_init,nlrho_init,order,windowL,rho,a0,v,maxstep,maxiter,tolerance,tolerance_max,easyeq_approx)
+    # max of 2pt correlation as a function of rho
+    elseif command=="2pt_max_rho"
+      easyeq_2pt_max_rho(rhos,dx,x0,minlrho_init,nlrho_init,order,windowL,a0,v,maxstep,maxiter,tolerance,tolerance_max,easyeq_approx)
+    # fourier transform of 2pt correlation
+    elseif command=="2pt_fourier"
+      easyeq_2pt_fourier(kmin,kmax,nk,minlrho_init,nlrho_init,order,windowL,rho,a0,v,maxiter,tolerance,easyeq_approx)
+    # max of fourier transform of 2pt correlation
+    elseif command=="2pt_fourier_max"
+      easyeq_2pt_fourier_max(dk,k0,minlrho_init,nlrho_init,order,windowL,rho,a0,v,maxstep,maxiter,tolerance,tolerance_max,easyeq_approx)
+    # max of 2pt correlation as a function of rho
+    elseif command=="2pt_fourier_max_rho"
+      easyeq_2pt_fourier_max_rho(rhos,dk,k0,minlrho_init,nlrho_init,order,windowL,a0,v,maxstep,maxiter,tolerance,tolerance_max,easyeq_approx)
+    # momentum distribution
+    elseif command=="momentum_distribution"
+      easyeq_momentum_distribution(kmin,kmax,minlrho_init,nlrho_init,order,windowL,rho,a0,v,maxiter,tolerance,easyeq_approx)
     else
       print(stderr,"unrecognized command '",command,"'.\n")
       exit(-1)
@@ -316,7 +393,7 @@ function main()
       anyeq_energy(rho,minlrho_init,nlrho_init,taus,P,N,J,a0,v,maxiter,tolerance,anyeq_approx,savefile)
     # e(rho)
     elseif command=="energy_rho"
-      anyeq_energy_rho(rhos,taus,P,N,J,a0,v,maxiter,tolerance,anyeq_approx,savefile)
+      anyeq_energy_rho(rhos,minlrho_init,nlrho_init,taus,P,N,J,a0,v,maxiter,tolerance,anyeq_approx,savefile)
     elseif command=="energy_rho_init_prevrho"
       anyeq_energy_rho_init_prevrho(rhos,taus,P,N,J,a0,v,maxiter,tolerance,anyeq_approx,savefile)
     elseif command=="energy_rho_init_nextrho"
@@ -331,12 +408,29 @@ function main()
     elseif command=="condensate_fraction"
       anyeq_condensate_fraction(rho,minlrho_init,nlrho_init,taus,P,N,J,a0,v,maxiter,tolerance,anyeq_approx,savefile)
     elseif command=="condensate_fraction_rho"
-      anyeq_condensate_fraction_rho(rhos,taus,P,N,J,a0,v,maxiter,tolerance,anyeq_approx,savefile)
+      anyeq_condensate_fraction_rho(rhos,minlrho_init,nlrho_init,taus,P,N,J,a0,v,maxiter,tolerance,anyeq_approx,savefile)
     # momentum distribution
     elseif command=="momentum_distribution"
-      anyeq_momentum_distribution(rho,minlrho_init,nlrho_init,taus,P,N,J,a0,v,maxiter,tolerance,anyeq_approx,savefile)
+      anyeq_momentum_distribution(kmin,kmax,rho,minlrho_init,nlrho_init,taus,P,N,J,windowL,a0,v,maxiter,tolerance,anyeq_approx,savefile)
     elseif command=="2pt"
       anyeq_2pt_correlation(minlrho_init,nlrho_init,taus,P,N,J,windowL,rho,a0,v,maxiter,tolerance,xmin,xmax,nx,anyeq_approx,savefile)
+    elseif command=="2pt_max"
+      anyeq_2pt_correlation_max(rho,minlrho_init,nlrho_init,dx,x0,maxstep,taus,P,N,J,windowL,a0,v,maxiter,tolerance,tolerance_max,anyeq_approx,savefile)
+    elseif command=="2pt_max_rho"
+      anyeq_2pt_correlation_max_rho(rhos,minlrho_init,nlrho_init,dx,x0,maxstep,taus,P,N,J,windowL,a0,v,maxiter,tolerance,tolerance_max,anyeq_approx,savefile)
+    elseif command=="2pt_fourier"
+      anyeq_2pt_correlation_fourier(minlrho_init,nlrho_init,taus,P,N,J,windowL,rho,a0,v,maxiter,tolerance,kmin,kmax,nk,anyeq_approx,savefile)
+    elseif command=="2pt_fourier_test"
+      anyeq_2pt_correlation_fourier_test(minlrho_init,nlrho_init,taus,P,N,J,windowL,rho,a0,v,maxiter,tolerance,xmax,kmin,kmax,nk,anyeq_approx,savefile)
+    elseif command=="2pt_fourier_max"
+      anyeq_2pt_correlation_fourier_max(rho,minlrho_init,nlrho_init,dk,k0,maxstep,taus,P,N,J,windowL,a0,v,maxiter,tolerance,tolerance_max,anyeq_approx,savefile)
+    elseif command=="2pt_fourier_max_rho"
+      anyeq_2pt_correlation_fourier_max_rho(rhos,minlrho_init,nlrho_init,dk,k0,maxstep,taus,P,N,J,windowL,a0,v,maxiter,tolerance,tolerance_max,anyeq_approx,savefile)
+    elseif command=="uncondensed_2pt"
+      anyeq_uncondensed_2pt_correlation(minlrho_init,nlrho_init,taus,P,N,J,windowL,rho,a0,v,maxiter,tolerance,xmin,xmax,nx,anyeq_approx,savefile)
+    # compressibility
+    elseif command=="compressibility_rho"
+      anyeq_compressibility_rho(rhos,minlrho_init,nlrho_init,taus,P,N,J,a0,v,maxiter,tolerance,anyeq_approx,savefile)
     else
       print(stderr,"unrecognized command: '",command,"'\n")
       exit(-1)
@@ -360,9 +454,11 @@ function main()
 end
 
 # parse a comma separated list as an array of Float64
-function parse_list(str)
+function parse_list(
+  str::String
+)
   elems=split(str,",")
-  out=Array{Float64}(undef,length(elems))
+  out=Array{Float64,1}(undef,length(elems))
   for i in 1:length(elems)
     out[i]=parse(Float64,elems[i])
   end
@@ -370,7 +466,9 @@ function parse_list(str)
 end
 
 # read cli arguments
-function read_args(ARGS)
+function read_args(
+  ARGS
+)
   # flag
   flag=""
 
diff --git a/src/multithread.jl b/src/multithread.jl
new file mode 100644
index 0000000..f61cdd7
--- /dev/null
+++ b/src/multithread.jl
@@ -0,0 +1,16 @@
+# split up 1...n among workers
+function spawn_workers(n::Int64)
+  # 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]=Int64[]
+  end
+  for i in 1:n
+    append!(work[(i-1)%nw+1],[i])
+  end
+
+  return work
+end
diff --git a/src/optimization.jl b/src/optimization.jl
new file mode 100644
index 0000000..bacead7
--- /dev/null
+++ b/src/optimization.jl
@@ -0,0 +1,94 @@
+# gradient descent: find local minimum of function of one variable from initial guess
+# numerically estimate the derivative
+@everywhere function gradient_descent(
+  f::Function,
+  x0::Float64,
+  delta::Float64, # shift is delta*df
+  dx::Float64, # finite difference for numerical derivative evaluation
+  maxiter::Int64 # interrupt and fail after maxiter steps
+)
+  counter=0
+
+  # init
+  x=x0
+
+  while counter<maxiter
+    # value at x and around
+    val=f(x)
+    valm=f(x-dx)
+    valp=f(x+dx)
+    # quit if minimum
+    if(val<valm && val<valp)
+      return(x,val)
+    end
+
+    # derivative
+    df=(valp-val)/dx
+    # step
+    x=x-delta*df
+    counter+=1
+  end
+
+  # fail
+  return(Inf,Inf)
+end
+
+# Newton algorithm to compute extrema
+# numerically estimate the derivatives
+@everywhere function newton_maximum(
+  f::Function,
+  x0::Float64,
+  dx::Float64, # finite difference for numerical derivative evaluation
+  maxiter::Int64,
+  tolerance::Float64,
+  maxstep::Float64 # maximal size of step
+)
+  counter=0
+
+  # init
+  x=x0
+
+  while counter<maxiter
+    # value at x and around
+    val=f(x)
+    valm=f(x-dx)
+    valp=f(x+dx)
+
+    # derivative
+    dfp=(valp-val)/dx
+    dfm=(val-valm)/dx
+    # second derivative
+    ddf=(dfp-dfm)/dx
+
+    #@printf(stderr,"% .15e % .15e % .15e % .15e\n",x,val,dfp,ddf)
+
+    if abs(dfp/ddf)<tolerance
+      # check it is a local maximum
+      if ddf<0
+	return(x,val)
+      else
+	return(Inf,Inf)
+      end
+    end
+
+    # step
+    step=dfp/abs(ddf)
+
+    # step too large
+    if abs(step)>maxstep
+      step=maxstep*sign(step)
+    end
+
+    x=x+step
+
+    # fail if off to infinity
+    if x==Inf || x==-Inf
+      return(x,val)
+    end
+    counter+=1
+  end
+
+  # fail
+  return(Inf,Inf)
+end
+
diff --git a/src/potentials.jl b/src/potentials.jl
index 46cafc0..b480db0 100644
--- a/src/potentials.jl
+++ b/src/potentials.jl
@@ -1,4 +1,4 @@
-## Copyright 2021 Ian Jauslin
+## 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.
@@ -13,10 +13,15 @@
 ## limitations under the License.
 
 # exponential potential in 3 dimensions
-@everywhere function v_exp(k,a)
+@everywhere function v_exp(
+  k::Float64,
+  a::Float64
+)
   return 8*pi/(1+k^2)^2*a
 end
-@everywhere function a0_exp(a)
+@everywhere function a0_exp(
+  a::Float64
+)
   if a>0.
     return log(a)+2*MathConstants.eulergamma+2*besselk(0,2*sqrt(a))/besseli(0,2*sqrt(a))
   elseif a<0.
@@ -27,15 +32,24 @@ end
 end
 
 # exp(-x)-a*exp(-b*x) in 3 dimensions
-@everywhere function v_expcry(k,a,b)
+@everywhere function v_expcry(
+  k::Float64,
+  a::Float64,
+  b::Float64
+)
   return 8*pi*((1+k^2)^(-2)-a*b*(b^2+k^2)^(-2))
 end
-@everywhere function a0_expcry(a,b)
+@everywhere function a0_expcry(
+  a::Float64,
+  b::Float64
+)
   return 1.21751642717932720441274114683413710125487579284827462 #ish
 end
 
 # x^2*exp(-|x|) in 3 dimensions
-@everywhere function v_npt(k)
+@everywhere function v_npt(
+  k::Float64
+)
   return 96*pi*(1-k^2)/(1+k^2)^4
 end
 @everywhere function a0_npt()
@@ -43,7 +57,9 @@ end
 end
 
 # 1/(1+x^4/4) potential in 3 dimensions
-@everywhere function v_alg(k)
+@everywhere function v_alg(
+  k::Float64
+)
   if(k==0)
     return 4*pi^2
   else
@@ -53,32 +69,50 @@ end
 a0_alg=1. #ish
 
 # (1+a x^4)/(1+x^2)^4 potential in 3 dimensions
-@everywhere function v_algwell(k)
+@everywhere function v_algwell(
+  k::Float64
+)
   a=4
   return pi^2/24*exp(-k)*(a*(k^2-9*k+15)+k^2+3*k+3)
 end
 a0_algwell=1. #ish
 
 # potential corresponding to the exact solution c/(1+b^2x^2)^2
-@everywhere function v_exact(k,b,c,e)
+@everywhere function v_exact(
+  k::Float64,
+  b::Float64,
+  c::Float64,
+  e::Float64
+)
   if k!=0
     return 48*pi^2*((18+3*sqrt(c)-(4-3*e/b^2)*c-(1-2*e/b^2)*c^1.5)/(4*(3+sqrt(c))^2*sqrt(c))*exp(-sqrt(1-sqrt(c))*k/b)+(-18+3*sqrt(c)+(4-3*e/b^2)*c-(1-2*e/b^2)*c^1.5)/(4*(3-sqrt(c))^2*sqrt(c))*exp(-sqrt(1+sqrt(c))*k/b)+(1-k/b)/2*exp(-k/b)-c*e/b^2*(3*(9-c)*k/b+8*c)/(8*(9-c)^2)*exp(-2*k/b))/k
   else
     return 48*pi^2*(-sqrt(1-sqrt(c))/b*(18+3*sqrt(c)-(4-3*e/b^2)*c-(1-2*e/b^2)*c^1.5)/(4*(3+sqrt(c))^2*sqrt(c))-sqrt(1+sqrt(c))/b*(-18+3*sqrt(c)+(4-3*e/b^2)*c-(1-2*e/b^2)*c^1.5)/(4*(3-sqrt(c))^2*sqrt(c))-1/b-c*e/b^2*(27-19*c)/(8*(9-c)^2))
   end
 end
-@everywhere function a0_exact(b,c,e)
+@everywhere function a0_exact(
+  b::Float64,
+  c::Float64,
+  e::Float64
+)
   return 1. #ish
 end
 
 # tent potential (convolution of soft sphere with itself): a*pi/12*(2*|x|/b-2)^2*(2*|x|/b+4) for |x|<b
-@everywhere function v_tent(k,a,b)
+@everywhere function v_tent(
+  k::Float64,
+  a::Float64,
+  b::Float64
+)
   if k!=0
     return (b/2)^3*a*(4*pi*(sin(k*b/2)-k*b/2*cos(k*b/2))/(k*b/2)^3)^2
   else
     return (b/2)^3*a*(4*pi/3)^2
   end
 end
-@everywhere function a0_tent(a,b)
+@everywhere function a0_tent(
+  a::Float64,
+  b::Float64
+)
   return b #ish
 end
diff --git a/src/print.jl b/src/print.jl
index bef1c4d..3587728 100644
--- a/src/print.jl
+++ b/src/print.jl
@@ -1,4 +1,4 @@
-## Copyright 2021 Ian Jauslin
+## 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.
diff --git a/src/simpleq-Kv.jl b/src/simpleq-Kv.jl
index 8789656..5a6579c 100644
--- a/src/simpleq-Kv.jl
+++ b/src/simpleq-Kv.jl
@@ -1,4 +1,4 @@
-## Copyright 2021 Ian Jauslin
+## 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.
@@ -13,12 +13,27 @@
 ## limitations under the License.
 
 # Compute Kv=(-\Delta+v+4e(1-\rho u*))^{-1}v
-function anyeq_Kv(minlrho,nlrho,taus,P,N,J,rho,a0,v,maxiter,tolerance,xmin,xmax,nx)
+function anyeq_Kv(
+  minlrho::Float64,
+  nlrho::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
+)
   # init vectors
   (weights,T,k,V,V0,A,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v)
 
   # compute initial guess from medeq
-  rhos=Array{Float64}(undef,nlrho)
+  rhos=Array{Float64,1}(undef,nlrho)
   for j in 0:nlrho-1
     rhos[j+1]=(nlrho==1 ? rho : 10^(minlrho+(log10(rho)-minlrho)/(nlrho-1)*j))
   end
@@ -44,7 +59,17 @@ function anyeq_Kv(minlrho,nlrho,taus,P,N,J,rho,a0,v,maxiter,tolerance,xmin,xmax,
 end
 
 # Compute the condensate fraction for simpleq using Kv
-function simpleq_Kv_condensate_fraction(rhos,taus,P,N,J,a0,v,maxiter,tolerance)
+function simpleq_Kv_condensate_fraction(
+  rhos::Array{Float64,1},
+  taus::Array{Float64,1},
+  P::Int64,
+  N::Int64,
+  J::Int64,
+  a0::Float64,
+  v::Function,
+  maxiter::Int64,
+  tolerance::Float64
+)
   # init vectors
   (weights,T,k,V,V0,A,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v)
 
@@ -67,19 +92,35 @@ function simpleq_Kv_condensate_fraction(rhos,taus,P,N,J,a0,v,maxiter,tolerance)
 end
 
 # Compute the two-point correlation function for simpleq using Kv
-function simpleq_Kv_2pt(minlrho,nlrho,taus,P,N,J,rho,a0,v,maxiter,tolerance,xmin,xmax,nx)
+function simpleq_Kv_2pt(
+  minlrho::Float64,
+  nlrho::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
+)
   # init vectors
   (weights,T,k,V,V0,A,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v)
 
   # compute initial guess from medeq
-  rhos=Array{Float64}(undef,nlrho)
+  rhos=Array{Float64,1}(undef,nlrho)
   for j in 0:nlrho-1
     rhos[j+1]=(nlrho==1 ? rho : 10^(minlrho+(log10(rho)-minlrho)/(nlrho-1)*j))
   end
   u0s=anyeq_init_medeq(rhos,N,J,k,a0,v,maxiter,tolerance)
   u0=u0s[nlrho]
 
-  (u,E,error)=anyeq_hatu(u0,P,N,J,rho,a0,weights,k,taus,V,V0,A,nothing,Upsilon,Upsilon0,v,maxiter,tolerance,Anyeq_approx(0.,0.,1.,0.,0.,0.,0.,0.,0.,0.,0.))
+  Abar=Array{Float64,5}(undef,0,0,0,0,0)
+  (u,E,error)=anyeq_hatu(u0,P,N,J,rho,a0,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,v,maxiter,tolerance,Anyeq_approx(0.,0.,1.,0.,0.,0.,0.,0.,0.,0.,0.))
 
   # Kv in Fourier space
   Kvk=simpleq_Kv_Kvk(u,V,E,rho,Upsilon,k,taus,weights,N,J)
@@ -103,7 +144,18 @@ function simpleq_Kv_2pt(minlrho,nlrho,taus,P,N,J,rho,a0,v,maxiter,tolerance,xmin
 end
 
 # Kv
-function simpleq_Kv_Kvk(u,V,E,rho,Upsilon,k,taus,weights,N,J)
+function simpleq_Kv_Kvk(
+  u::Array{Float64,1},
+  V::Array{Float64,1},
+  E::Float64,
+  rho::Float64,
+  Upsilon::Array{Array{Float64,1},1},
+  k::Array{Float64,1},
+  taus::Array{Float64,1},
+  weights::Tuple{Array{Float64,1},Array{Float64,1}},
+  N::Int64,
+  J::Int64
+)
   # (-Delta+v+4e(1-\rho u*)) in Fourier space
   M=Array{Float64,2}(undef,N*J,N*J)
   for zetapp in 0:J-1
diff --git a/src/simpleq-hardcore.jl b/src/simpleq-hardcore.jl
index ca64f78..398ac06 100644
--- a/src/simpleq-hardcore.jl
+++ b/src/simpleq-hardcore.jl
@@ -1,4 +1,4 @@
-## Copyright 2021 Ian Jauslin
+## 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.
@@ -13,28 +13,26 @@
 ## limitations under the License.
 
 # compute energy as a function of rho
-function simpleq_hardcore_energy_rho(rhos,taus,P,N,J,maxiter,tolerance)
-  ## 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%nw+1],j)
-  end
+function simpleq_hardcore_energy_rho(
+  rhos::Array{Float64,1},
+  taus::Array{Float64,1},
+  P::Int64,
+  N::Int64,
+  J::Int64,
+  maxiter::Int64,
+  tolerance::Float64
+)
+  # spawn workers
+  work=spawn_workers(length(rhos))
 
   # initialize vectors
   (weights,weights_gL,r,T)=simpleq_hardcore_init(taus,P,N,J)
 
   # initial guess
-  u0s=Array{Array{Float64}}(undef,length(rhos))
-  e0s=Array{Float64}(undef,length(rhos))
+  u0s=Array{Array{Float64,1}}(undef,length(rhos))
+  e0s=Array{Float64,1}(undef,length(rhos))
   for j in 1:length(rhos)
-    u0s[j]=Array{Float64}(undef,N*J)
+    u0s[j]=Array{Float64,1}(undef,N*J)
     for i in 1:N*J
       u0s[j][i]=1/(1+r[i]^2)^2
     end
@@ -43,17 +41,17 @@ function simpleq_hardcore_energy_rho(rhos,taus,P,N,J,maxiter,tolerance)
 
 
   # save result from each task
-  us=Array{Array{Float64}}(undef,length(rhos))
-  es=Array{Float64}(undef,length(rhos))
-  err=Array{Float64}(undef,length(rhos))
+  us=Array{Array{Float64,1}}(undef,length(rhos))
+  es=Array{Float64,1}(undef,length(rhos))
+  err=Array{Float64,1}(undef,length(rhos))
 
   count=0
   # for each worker
-  @sync for p in 1:nw
+  @sync for p in 1:length(work)
     # for each task
     @async for j in work[p]
       count=count+1
-      if count>=nw
+      if count>=length(work)
         progress(count,length(rhos),10000)
       end
       # run the task
@@ -67,12 +65,20 @@ function simpleq_hardcore_energy_rho(rhos,taus,P,N,J,maxiter,tolerance)
 end
 
 # compute u(x)
-function simpleq_hardcore_ux(rho,taus,P,N,J,maxiter,tolerance)
+function simpleq_hardcore_ux(
+  rho::Float64,
+  taus::Array{Float64,1},
+  P::Int64,
+  N::Int64,
+  J::Int64,
+  maxiter::Int64,
+  tolerance::Float64
+)
   # initialize vectors
   (weights,weights_gL,r,T)=simpleq_hardcore_init(taus,P,N,J)
 
   # initial guess
-  u0=Array{Float64}(undef,N*J)
+  u0=Array{Float64,1}(undef,N*J)
   for i in 1:N*J
     u0[i]=1/(1+r[i]^2)^2
   end
@@ -86,28 +92,26 @@ function simpleq_hardcore_ux(rho,taus,P,N,J,maxiter,tolerance)
 end
 
 # compute condensate fraction as a function of rho
-function simpleq_hardcore_condensate_fraction_rho(rhos,taus,P,N,J,maxiter,tolerance)
-  ## 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%nw+1],j)
-  end
+function simpleq_hardcore_condensate_fraction_rho(
+  rhos::Array{Float64,1},
+  taus::Array{Float64,1},
+  P::Int64,
+  N::Int64,
+  J::Int64,
+  maxiter::Int64,
+  tolerance::Float64
+)
+  # spawn workers
+  work=spawn_workers(length(rhos))
 
   # initialize vectors
   (weights,weights_gL,r,T)=simpleq_hardcore_init(taus,P,N,J)
 
   # initial guess
-  u0s=Array{Array{Float64}}(undef,length(rhos))
-  e0s=Array{Float64}(undef,length(rhos))
+  u0s=Array{Array{Float64,1}}(undef,length(rhos))
+  e0s=Array{Float64,1}(undef,length(rhos))
   for j in 1:length(rhos)
-    u0s[j]=Array{Float64}(undef,N*J)
+    u0s[j]=Array{Float64,1}(undef,N*J)
     for i in 1:N*J
       u0s[j][i]=1/(1+r[i]^2)^2
     end
@@ -116,17 +120,17 @@ function simpleq_hardcore_condensate_fraction_rho(rhos,taus,P,N,J,maxiter,tolera
 
 
   # save result from each task
-  us=Array{Array{Float64}}(undef,length(rhos))
-  es=Array{Float64}(undef,length(rhos))
-  err=Array{Float64}(undef,length(rhos))
+  us=Array{Array{Float64,1}}(undef,length(rhos))
+  es=Array{Float64,1}(undef,length(rhos))
+  err=Array{Float64,1}(undef,length(rhos))
 
   count=0
   # for each worker
-  @sync for p in 1:nw
+  @sync for p in 1:length(work)
     # for each task
     @async for j in work[p]
       count=count+1
-      if count>=nw
+      if count>=length(work)
         progress(count,length(rhos),10000)
       end
       # run the task
@@ -142,13 +146,18 @@ end
 
 
 # initialize computation
-@everywhere function simpleq_hardcore_init(taus,P,N,J)
+@everywhere function simpleq_hardcore_init(
+  taus::Array{Float64,1},
+  P::Int64,
+  N::Int64,
+  J::Int64
+)
   # Gauss-Legendre weights
   weights=gausslegendre(N)
   weights_gL=gausslaguerre(N)
 
   # r
-  r=Array{Float64}(undef,J*N)
+  r=Array{Float64,1}(undef,J*N)
   for zeta in 0:J-1
     for j in 1:N
       xj=weights[1][j]
@@ -164,9 +173,23 @@ end
 end
 
 # compute u using chebyshev expansions
-@everywhere function simpleq_hardcore_hatu(u0,e0,rho,r,taus,T,weights,weights_gL,P,N,J,maxiter,tolerance)
+@everywhere function simpleq_hardcore_hatu(
+  u0::Array{Float64,1},
+  e0::Float64,
+  rho::Float64,
+  r::Array{Float64,1},
+  taus::Array{Float64,1},
+  T::Array{Polynomial,1},
+  weights::Tuple{Array{Float64,1},Array{Float64,1}},
+  weights_gL::Tuple{Array{Float64,1},Array{Float64,1}},
+  P::Int64,
+  N::Int64,
+  J::Int64,
+  maxiter::Int64,
+  tolerance::Float64
+)
   # init
-  vec=Array{Float64}(undef,J*N+1)
+  vec=Array{Float64,1}(undef,J*N+1)
   for i in 1:J*N
     vec[i]=u0[i]
   end
@@ -194,8 +217,20 @@ end
 end
 
 # Xi
-@everywhere function simpleq_hardcore_Xi(u,e,rho,r,taus,T,weights,weights_gL,P,N,J)
-  out=Array{Float64}(undef,J*N+1)
+@everywhere function simpleq_hardcore_Xi(
+  u::Array{Float64,1},
+  e::Float64,
+  rho::Float64,
+  r::Array{Float64,1},
+  taus::Array{Float64,1},
+  T::Array{Polynomial,1},
+  weights::Tuple{Array{Float64,1},Array{Float64,1}},
+  weights_gL::Tuple{Array{Float64,1},Array{Float64,1}},
+  P::Int64,
+  N::Int64,
+  J::Int64
+)
+  out=Array{Float64,1}(undef,J*N+1)
   FU=chebyshev(u,taus,weights,P,N,J,4)
 
   #D's
@@ -215,7 +250,19 @@ end
   return out
 end
 # DXi
-@everywhere function simpleq_hardcore_DXi(u,e,rho,r,taus,T,weights,weights_gL,P,N,J)
+@everywhere function simpleq_hardcore_DXi(
+  u::Array{Float64,1},
+  e::Float64,
+  rho::Float64,
+  r::Array{Float64,1},
+  taus::Array{Float64,1},
+  T::Array{Polynomial,1},
+  weights::Tuple{Array{Float64,1},Array{Float64,1}},
+  weights_gL::Tuple{Array{Float64,1},Array{Float64,1}},
+  P::Int64,
+  N::Int64,
+  J::Int64
+)
   out=Array{Float64,2}(undef,J*N+1,J*N+1)
   FU=chebyshev(u,taus,weights,P,N,J,4)
 
@@ -232,15 +279,15 @@ end
 
   for zetapp in 0:J-1
     for n in 1:N
-      one=zeros(Int64,J*N)
-      one[zetapp*N+n]=1
+      one=zeros(Float64,J*N)
+      one[zetapp*N+n]=1.
       Fone=chebyshev(one,taus,weights,P,N,J,4)
 
       for i in 1:J*N
 	# du/du
 	out[i,zetapp*N+n]=dot(Fone,d1[i])+2*dot(FU,d2[i]*Fone)-(zetapp*N+n==i ? 1 : 0)
 	# du/de
-	out[i,J*N+1]=(dsed0[i]+dot(FU,dsed1[i])+dot(FU,dsed2[i]*FU))/(2*sqrt(abs(e)))*(e>=0 ? 1 : -1)
+	out[i,J*N+1]=(dsed0[i]+dot(FU,dsed1[i])+dot(FU,dsed2[i]*FU))/(2*sqrt(abs(e)))*(e>=0. ? 1. : -1.)
       end
       # de/du
       out[J*N+1,zetapp*N+n]=2*pi*rho*
@@ -253,14 +300,26 @@ end
   #de/de
   out[J*N+1,J*N+1]=-1+2*pi*rho*
     (2-dsedgamma0(e,rho,weights)-dot(FU,dsedgamma1(e,rho,taus,T,weights,weights_gL,P,J,4))-dot(FU,dsedgamma2(e,rho,taus,T,weights,weights_gL,P,J,4)*FU))/denom/
-    (2*sqrt(abs(e)))*(e>=0 ? 1 : -1)
+    (2*sqrt(abs(e)))*(e>=0. ? 1. : -1.)
 
   return out
 end
 
 # dXi/dmu
-@everywhere function simpleq_hardcore_dXidmu(u,e,rho,r,taus,T,weights,weights_gL,P,N,J)
-  out=Array{Float64}(undef,J*N+1)
+@everywhere function simpleq_hardcore_dXidmu(
+  u::Array{Float64,1},
+  e::Float64,
+  rho::Float64,
+  r::Array{Float64,1},
+  taus::Array{Float64,1},
+  T::Array{Polynomial,1},
+  weights::Tuple{Array{Float64,1},Array{Float64,1}},
+  weights_gL::Tuple{Array{Float64,1},Array{Float64,1}},
+  P::Int64,
+  N::Int64,
+  J::Int64
+)
+  out=Array{Float64,1}(undef,J*N+1)
   FU=chebyshev(u,taus,weights,P,N,J,4)
 
   #D's
@@ -281,17 +340,37 @@ end
 end
 
 # B's
-@everywhere function B0(r)
+@everywhere function B0(
+  r::Float64
+)
   return pi/12*(r-1)^2*(r+5)
 end
-@everywhere function B1(r,zeta,n,taus,T,weights,nu)
+@everywhere function B1(
+  r::Float64,
+  zeta::Int64,
+  n::Int64,
+  taus::Array{Float64,1},
+  T::Array{Polynomial,1},
+  weights::Tuple{Array{Float64,1},Array{Float64,1}},
+  nu::Int64
+)
   return (taus[zeta+1]>=(2-r)/r || taus[zeta+2]<=-r/(r+2) ? 0 :
     8*pi/(r+1)*integrate_legendre(tau->
       (1-(r-(1-tau)/(1+tau))^2)/(1+tau)^(3-nu)*T[n+1]((2*tau-(taus[zeta+1]+taus[zeta+2]))/(taus[zeta+2]-taus[zeta+1])),max(taus[zeta+1],-r/(r+2))
     ,min(taus[zeta+2],(2-r)/r),weights)
   )
 end
-@everywhere function B2(r,zeta,n,zetap,m,taus,T,weights,nu)
+@everywhere function B2(
+  r::Float64,
+  zeta::Int64,
+  n::Int64,
+  zetap::Int64,
+  m::Int64,
+  taus::Array{Float64,1},
+  T::Array{Polynomial,1},
+  weights::Tuple{Array{Float64,1},Array{Float64,1}},
+  nu::Int64
+)
   return 32*pi/(r+1)*integrate_legendre(tau->
     1/(1+tau)^(3-nu)*T[n+1]((2*tau-(taus[zeta+1]+taus[zeta+2]))/(taus[zeta+2]-taus[zeta+1]))*
     (taus[zetap+1]>=alphap(abs(r-(1-tau)/(1+tau))-2*tau/(1+tau),tau) || taus[zetap+2]<=alpham(1+r,tau) ? 0 :
@@ -303,30 +382,49 @@ end
 end
 
 # D's
-@everywhere function D0(e,rho,r,weights,N,J)
-  out=Array{Float64}(undef,J*N)
+@everywhere function D0(
+  e::Float64,
+  rho::Float64,
+  r::Array{Float64,1},
+  weights::Tuple{Array{Float64,1},Array{Float64,1}},
+  N::Int64,
+  J::Int64
+)
+  out=Array{Float64,1}(undef,J*N)
   for i in 1:J*N
     out[i]=exp(-2*sqrt(abs(e))*r[i])/(r[i]+1)+
       rho*sqrt(abs(e))/(r[i]+1)*integrate_legendre(s->
 	(s+1)*B0(s)*(exp(-2*sqrt(abs(e))*(r[i]-s))-exp(-2*sqrt(abs(e))*(r[i]+s)))/2
-      ,0,min(1,r[i]),weights)+
+      ,0.,min(1.,r[i]),weights)+
       (r[i]>=1 ? 0 :
 	rho*sqrt(abs(e))/(2*(r[i]+1))*(1-exp(-4*sqrt(abs(e))*r[i]))*integrate_legendre(s->
 	  (s+r[i]+1)*B0(s+r[i])*exp(-2*sqrt(abs(e))*s)
-	,0,1-r[i],weights)
+	,0.,1. -r[i],weights)
       )
   end
   return out
 end
-@everywhere function D1(e,rho,r,taus,T,weights,weights_gL,P,N,J,nu)
-  out=Array{Array{Float64}}(undef,J*N)
+@everywhere function D1(
+  e::Float64,
+  rho::Float64,
+  r::Array{Float64,1},
+  taus::Array{Float64,1},
+  T::Array{Polynomial,1},
+  weights::Tuple{Array{Float64,1},Array{Float64,1}},
+  weights_gL::Tuple{Array{Float64,1},Array{Float64,1}},
+  P::Int64,
+  N::Int64,
+  J::Int64,
+  nu::Int64
+)
+  out=Array{Array{Float64,1},1}(undef,J*N)
   for i in 1:J*N
-    out[i]=Array{Float64}(undef,(P+1)*J)
+    out[i]=Array{Float64,1}(undef,(P+1)*J)
     for zeta in 0:J-1
       for n in 0:P
 	out[i][zeta*(P+1)+n+1]=rho*sqrt(abs(e))/(r[i]+1)*integrate_legendre(s->
 	    (s+1)*B1(s,zeta,n,taus,T,weights,nu)*(exp(-2*sqrt(abs(e))*(r[i]-s))-exp(-2*sqrt(abs(e))*(r[i]+s)))/2
-	  ,0,r[i],weights)+
+	  ,0.,r[i],weights)+
 	  rho*sqrt(abs(e))/(2*(r[i]+1))*(1-exp(-4*sqrt(abs(e))*r[i]))*integrate_laguerre(s->
 	    (s+r[i]+1)*B1(s+r[i],zeta,n,taus,T,weights,nu)
 	  ,2*sqrt(abs(e)),weights_gL)
@@ -336,7 +434,19 @@ end
 
   return out
 end
-@everywhere function D2(e,rho,r,taus,T,weights,weights_gL,P,N,J,nu)
+@everywhere function D2(
+  e::Float64,
+  rho::Float64,
+  r::Array{Float64,1},
+  taus::Array{Float64,1},
+  T::Array{Polynomial,1},
+  weights::Tuple{Array{Float64,1},Array{Float64,1}},
+  weights_gL::Tuple{Array{Float64,1},Array{Float64,1}},
+  P::Int64,
+  N::Int64,
+  J::Int64,
+  nu::Int64
+)
   out=Array{Array{Float64,2}}(undef,J*N)
   for i in 1:J*N
     out[i]=Array{Float64,2}(undef,(P+1)*J,(P+1)*J)
@@ -346,7 +456,7 @@ end
 	  for m in 0:P
 	    out[i][zeta*(P+1)+n+1,zetap*(P+1)+m+1]=rho*sqrt(abs(e))/(r[i]+1)*integrate_legendre(s->
 		(s+1)*B2(s,zeta,n,zetap,m,taus,T,weights,nu)*(exp(-2*sqrt(abs(e))*(r[i]-s))-exp(-2*sqrt(abs(e))*(r[i]+s)))/2
-	      ,0,r[i],weights)+
+	      ,0.,r[i],weights)+
 	      rho*sqrt(abs(e))/(2*(r[i]+1))*(1-exp(-4*sqrt(abs(e))*r[i]))*integrate_laguerre(s->
 		(s+r[i]+1)*B2(s+r[i],zeta,n,zetap,m,taus,T,weights,nu)
 	      ,2*sqrt(abs(e)),weights_gL)
@@ -359,28 +469,47 @@ end
 end
 
 # dD/d sqrt(abs(e))'s
-@everywhere function dseD0(e,rho,r,weights,N,J)
-  out=Array{Float64}(undef,J*N)
+@everywhere function dseD0(
+  e::Float64,
+  rho::Float64,
+  r::Array{Float64,1},
+  weights::Tuple{Array{Float64,1},Array{Float64,1}},
+  N::Int64,
+  J::Int64
+)
+  out=Array{Float64,1}(undef,J*N)
   for i in 1:J*N
     out[i]=-2*r[i]*exp(-2*sqrt(abs(e))*r[i])/(r[i]+1)+
       rho/(r[i]+1)*integrate_legendre(s->
 	(s+1)*B0(s)*((1-2*sqrt(abs(e))*r[i])*(exp(-2*sqrt(abs(e))*(r[i]-s))-exp(-2*sqrt(abs(e))*(r[i]+s)))/2+2*sqrt(abs(e))*s*(exp(-2*sqrt(abs(e))*(r[i]-s))+exp(-2*sqrt(abs(e))*(r[i]+s)))/2)
-      ,0,min(1,r[i]),weights)+
+      ,0.,min(1.,r[i]),weights)+
       (r[i]>=1 ? 0 :
-	rho/(2*(r[i]+1))*integrate_legendre(s->(s+r[i]+1)*B0(s+r[i])*((1-2*sqrt(abs(e))*s)*(1-exp(-4*sqrt(abs(e))*r[i]))*exp(-2*sqrt(abs(e))*s)+4*sqrt(abs(e))*r[i]*exp(-2*sqrt(abs(e))*(2*r[i]+s))),0,1-r[i],weights)
+	rho/(2*(r[i]+1))*integrate_legendre(s->(s+r[i]+1)*B0(s+r[i])*((1-2*sqrt(abs(e))*s)*(1-exp(-4*sqrt(abs(e))*r[i]))*exp(-2*sqrt(abs(e))*s)+4*sqrt(abs(e))*r[i]*exp(-2*sqrt(abs(e))*(2*r[i]+s))),0.,1. -r[i],weights)
       )
   end
   return out
 end
-@everywhere function dseD1(e,rho,r,taus,T,weights,weights_gL,P,N,J,nu)
-  out=Array{Array{Float64}}(undef,J*N)
+@everywhere function dseD1(
+  e::Float64,
+  rho::Float64,
+  r::Array{Float64,1},
+  taus::Array{Float64,1},
+  T::Array{Polynomial,1},
+  weights::Tuple{Array{Float64,1},Array{Float64,1}},
+  weights_gL::Tuple{Array{Float64,1},Array{Float64,1}},
+  P::Int64,
+  N::Int64,
+  J::Int64,
+  nu::Int64
+)
+  out=Array{Array{Float64,1},1}(undef,J*N)
   for i in 1:J*N
-    out[i]=Array{Float64}(undef,(P+1)*J)
+    out[i]=Array{Float64,1}(undef,(P+1)*J)
     for zeta in 0:J-1
       for n in 0:P
 	out[i][zeta*(P+1)+n+1]=rho/(r[i]+1)*integrate_legendre(s->
 	  (s+1)*B1(s,zeta,n,taus,T,weights,nu)*((1-2*sqrt(abs(e))*r[i])*(exp(-2*sqrt(abs(e))*(r[i]-s))-exp(-2*sqrt(abs(e))*(r[i]+s)))/2+2*sqrt(abs(e))*s*(exp(-2*sqrt(abs(e))*(r[i]-s))+exp(-2*sqrt(abs(e))*(r[i]+s)))/2)
-	,0,r[i],weights)+
+	,0.,r[i],weights)+
 	rho/(2*(r[i]+1))*integrate_laguerre(s->
 	  (s+r[i]+1)*B1(s+r[i],zeta,n,taus,T,weights,nu)*((1-2*sqrt(abs(e))*s)*(1-exp(-4*sqrt(abs(e))*r[i]))+4*sqrt(abs(e))*r[i]*exp(-4*sqrt(abs(e))*r[i]))
 	,2*sqrt(abs(e)),weights_gL)
@@ -389,7 +518,19 @@ end
   end
   return out
 end
-@everywhere function dseD2(e,rho,r,taus,T,weights,weights_gL,P,N,J,nu)
+@everywhere function dseD2(
+  e::Float64,
+  rho::Float64,
+  r::Array{Float64,1},
+  taus::Array{Float64,1},
+  T::Array{Polynomial,1},
+  weights::Tuple{Array{Float64,1},Array{Float64,1}},
+  weights_gL::Tuple{Array{Float64,1},Array{Float64,1}},
+  P::Int64,
+  N::Int64,
+  J::Int64,
+  nu::Int64
+)
   out=Array{Array{Float64,2}}(undef,J*N)
   for i in 1:J*N
     out[i]=Array{Float64,2}(undef,(P+1)*J,(P+1)*J)
@@ -399,7 +540,7 @@ end
 	  for m in 0:P
 	    out[i][zeta*(P+1)+n+1,zetap*(P+1)+m+1]=rho/(r[i]+1)*integrate_legendre(s->
 	      (s+1)*B2(s,zeta,n,zetap,m,taus,T,weights,nu)*((1-2*sqrt(abs(e))*r[i])*(exp(-2*sqrt(abs(e))*(r[i]-s))-exp(-2*sqrt(abs(e))*(r[i]+s)))/2+2*sqrt(abs(e))*s*(exp(-2*sqrt(abs(e))*(r[i]-s))+exp(-2*sqrt(abs(e))*(r[i]+s)))/2)
-	    ,0,r[i],weights)+
+	    ,0.,r[i],weights)+
 	    rho/(2*(r[i]+1))*integrate_laguerre(s->
 	      (s+r[i]+1)*B2(s+r[i],zeta,n,zetap,m,taus,T,weights,nu)*((1-2*sqrt(abs(e))*s)*(1-exp(-4*sqrt(abs(e))*r[i]))+4*sqrt(abs(e))*r[i]*exp(-4*sqrt(abs(e))*r[i]))
 	    ,2*sqrt(abs(e)),weights_gL)
@@ -412,28 +553,47 @@ end
 end
 
 # dD/d sqrt(abs(e+mu/2))'s
-@everywhere function dsmuD0(e,rho,r,weights,N,J)
-  out=Array{Float64}(undef,J*N)
+@everywhere function dsmuD0(
+  e::Float64,
+  rho::Float64,
+  r::Array{Float64,1},
+  weights::Tuple{Array{Float64,1},Array{Float64,1}},
+  N::Int64,
+  J::Int64
+)
+  out=Array{Float64,1}(undef,J*N)
   for i in 1:J*N
     out[i]=-2*r[i]*exp(-2*sqrt(abs(e))*r[i])/(r[i]+1)+
       rho/(r[i]+1)*integrate_legendre(s->
 	(s+1)*B0(s)*((-1-2*sqrt(abs(e))*r[i])*(exp(-2*sqrt(abs(e))*(r[i]-s))-exp(-2*sqrt(abs(e))*(r[i]+s)))/2+2*sqrt(abs(e))*s*(exp(-2*sqrt(abs(e))*(r[i]-s))+exp(-2*sqrt(abs(e))*(r[i]+s)))/2)
-      ,0,min(1,r[i]),weights)+
+      ,0.,min(1.,r[i]),weights)+
       (r[i]>=1 ? 0 :
-	rho/(2*(r[i]+1))*integrate_legendre(s->(s+r[i]+1)*B0(s+r[i])*((-1-2*sqrt(abs(e))*s)*(1-exp(-4*sqrt(abs(e))*r[i]))*exp(-2*sqrt(abs(e))*s)+4*sqrt(abs(e))*r[i]*exp(-2*sqrt(abs(e))*(2*r[i]+s))),0,1-r[i],weights)
+	rho/(2*(r[i]+1))*integrate_legendre(s->(s+r[i]+1)*B0(s+r[i])*((-1-2*sqrt(abs(e))*s)*(1-exp(-4*sqrt(abs(e))*r[i]))*exp(-2*sqrt(abs(e))*s)+4*sqrt(abs(e))*r[i]*exp(-2*sqrt(abs(e))*(2*r[i]+s))),0.,1. -r[i],weights)
       )
   end
   return out
 end
-@everywhere function dsmuD1(e,rho,r,taus,T,weights,weights_gL,P,N,J,nu)
-  out=Array{Array{Float64}}(undef,J*N)
+@everywhere function dsmuD1(
+  e::Float64,
+  rho::Float64,
+  r::Array{Float64,1},
+  taus::Array{Float64,1},
+  T::Array{Polynomial,1},
+  weights::Tuple{Array{Float64,1},Array{Float64,1}},
+  weights_gL::Tuple{Array{Float64,1},Array{Float64,1}},
+  P::Int64,
+  N::Int64,
+  J::Int64,
+  nu::Int64
+)
+  out=Array{Array{Float64,1},1}(undef,J*N)
   for i in 1:J*N
-    out[i]=Array{Float64}(undef,(P+1)*J)
+    out[i]=Array{Float64,1}(undef,(P+1)*J)
     for zeta in 0:J-1
       for n in 0:P
 	out[i][zeta*(P+1)+n+1]=rho/(r[i]+1)*integrate_legendre(s->
 	  (s+1)*B1(s,zeta,n,taus,T,weights,nu)*((-1-2*sqrt(abs(e))*r[i])*(exp(-2*sqrt(abs(e))*(r[i]-s))-exp(-2*sqrt(abs(e))*(r[i]+s)))/2+2*sqrt(abs(e))*s*(exp(-2*sqrt(abs(e))*(r[i]-s))+exp(-2*sqrt(abs(e))*(r[i]+s)))/2)
-	,0,r[i],weights)+
+	,0.,r[i],weights)+
 	rho/(2*(r[i]+1))*integrate_laguerre(s->
 	  (s+r[i]+1)*B1(s+r[i],zeta,n,taus,T,weights,nu)*((-1-2*sqrt(abs(e))*s)*(1-exp(-4*sqrt(abs(e))*r[i]))+4*sqrt(abs(e))*r[i]*exp(-4*sqrt(abs(e))*r[i]))
 	,2*sqrt(abs(e)),weights_gL)
@@ -442,7 +602,19 @@ end
   end
   return out
 end
-@everywhere function dsmuD2(e,rho,r,taus,T,weights,weights_gL,P,N,J,nu)
+@everywhere function dsmuD2(
+  e::Float64,
+  rho::Float64,
+  r::Array{Float64,1},
+  taus::Array{Float64,1},
+  T::Array{Polynomial,1},
+  weights::Tuple{Array{Float64,1},Array{Float64,1}},
+  weights_gL::Tuple{Array{Float64,1},Array{Float64,1}},
+  P::Int64,
+  N::Int64,
+  J::Int64,
+  nu::Int64
+)
   out=Array{Array{Float64,2}}(undef,J*N)
   for i in 1:J*N
     out[i]=Array{Float64,2}(undef,(P+1)*J,(P+1)*J)
@@ -452,7 +624,7 @@ end
 	  for m in 0:P
 	    out[i][zeta*(P+1)+n+1,zetap*(P+1)+m+1]=rho/(r[i]+1)*integrate_legendre(s->
 	      (s+1)*B2(s,zeta,n,zetap,m,taus,T,weights,nu)*((-1-2*sqrt(abs(e))*r[i])*(exp(-2*sqrt(abs(e))*(r[i]-s))-exp(-2*sqrt(abs(e))*(r[i]+s)))/2+2*sqrt(abs(e))*s*(exp(-2*sqrt(abs(e))*(r[i]-s))+exp(-2*sqrt(abs(e))*(r[i]+s)))/2)
-	    ,0,r[i],weights)+
+	    ,0.,r[i],weights)+
 	    rho/(2*(r[i]+1))*integrate_laguerre(s->
 	      (s+r[i]+1)*B2(s+r[i],zeta,n,zetap,m,taus,T,weights,nu)*((-1-2*sqrt(abs(e))*s)*(1-exp(-4*sqrt(abs(e))*r[i]))+4*sqrt(abs(e))*r[i]*exp(-4*sqrt(abs(e))*r[i]))
 	    ,2*sqrt(abs(e)),weights_gL)
@@ -465,11 +637,25 @@ end
 end
 
 # gamma's
-@everywhere function gamma0(e,rho,weights)
-  return 2*rho*e*integrate_legendre(s->(s+1)*B0(s)*exp(-2*sqrt(abs(e))*s),0,1,weights)
+@everywhere function gamma0(
+  e::Float64,
+  rho::Float64,
+  weights::Tuple{Array{Float64,1},Array{Float64,1}}
+)
+  return 2*rho*e*integrate_legendre(s->(s+1)*B0(s)*exp(-2*sqrt(abs(e))*s),0.,1.,weights)
 end
-@everywhere function gamma1(e,rho,taus,T,weights,weights_gL,P,J,nu)
-  out=Array{Float64}(undef,J*(P+1))
+@everywhere function gamma1(
+  e::Float64,
+  rho::Float64,
+  taus::Array{Float64,1},
+  T::Array{Polynomial,1},
+  weights::Tuple{Array{Float64,1},Array{Float64,1}},
+  weights_gL::Tuple{Array{Float64,1},Array{Float64,1}},
+  P::Int64,
+  J::Int64,
+  nu::Int64
+)
+  out=Array{Float64,1}(undef,J*(P+1))
   for zeta in 0:J-1
     for n in 0:P
       out[zeta*(P+1)+n+1]=2*rho*e*integrate_laguerre(s->(s+1)*B1(s,zeta,n,taus,T,weights,nu),2*sqrt(abs(e)),weights_gL)
@@ -477,7 +663,17 @@ end
   end
   return out
 end
-@everywhere function gamma2(e,rho,taus,T,weights,weights_gL,P,J,nu)
+@everywhere function gamma2(
+  e::Float64,
+  rho::Float64,
+  taus::Array{Float64,1},
+  T::Array{Polynomial,1},
+  weights::Tuple{Array{Float64,1},Array{Float64,1}},
+  weights_gL::Tuple{Array{Float64,1},Array{Float64,1}},
+  P::Int64,
+  J::Int64,
+  nu::Int64
+)
   out=Array{Float64,2}(undef,J*(P+1),J*(P+1))
   for zeta in 0:J-1
     for n in 0:P
@@ -492,11 +688,25 @@ end
 end
 
 # dgamma/d sqrt(abs(e))'s
-@everywhere function dsedgamma0(e,rho,weights)
-  return 4*rho*sqrt(abs(e))*integrate_legendre(s->(s+1)*B0(s)*(1-sqrt(abs(e))*s)*exp(-2*sqrt(abs(e))*s),0,1,weights)
+@everywhere function dsedgamma0(
+  e::Float64,
+  rho::Float64,
+  weights::Tuple{Array{Float64,1},Array{Float64,1}}
+)
+  return 4*rho*sqrt(abs(e))*integrate_legendre(s->(s+1)*B0(s)*(1-sqrt(abs(e))*s)*exp(-2*sqrt(abs(e))*s),0.,1.,weights)
 end
-@everywhere function dsedgamma1(e,rho,taus,T,weights,weights_gL,P,J,nu)
-  out=Array{Float64}(undef,J*(P+1))
+@everywhere function dsedgamma1(
+  e::Float64,
+  rho::Float64,
+  taus::Array{Float64,1},
+  T::Array{Polynomial,1},
+  weights::Tuple{Array{Float64,1},Array{Float64,1}},
+  weights_gL::Tuple{Array{Float64,1},Array{Float64,1}},
+  P::Int64,
+  J::Int64,
+  nu::Int64
+)
+  out=Array{Float64,1}(undef,J*(P+1))
   for zeta in 0:J-1
     for n in 0:P
       out[zeta*(P+1)+n+1]=4*rho*e*integrate_laguerre(s->(s+1)*B1(s,zeta,n,taus,T,weights,nu)*(1-sqrt(abs(e))*s),2*sqrt(abs(e)),weights_gL)
@@ -504,7 +714,17 @@ end
   end
   return out
 end
-@everywhere function dsedgamma2(e,rho,taus,T,weights,weights_gL,P,J,nu)
+@everywhere function dsedgamma2(
+  e::Float64,
+  rho::Float64,
+  taus::Array{Float64,1},
+  T::Array{Polynomial,1},
+  weights::Tuple{Array{Float64,1},Array{Float64,1}},
+  weights_gL::Tuple{Array{Float64,1},Array{Float64,1}},
+  P::Int64,
+  J::Int64,
+  nu::Int64
+)
   out=Array{Float64,2}(undef,J*(P+1),J*(P+1))
   for zeta in 0:J-1
     for n in 0:P
@@ -519,11 +739,25 @@ end
 end
 
 # dgamma/d sqrt(e+mu/2)'s
-@everywhere function dsmudgamma0(e,rho,weights)
-  return -4*rho*e*integrate_legendre(s->(s+1)*s*B0(s)*exp(-2*sqrt(abs(e))*s),0,1,weights)
+@everywhere function dsmudgamma0(
+  e::Float64,
+  rho::Float64,
+  weights::Tuple{Array{Float64,1},Array{Float64,1}}
+)
+  return -4*rho*e*integrate_legendre(s->(s+1)*s*B0(s)*exp(-2*sqrt(abs(e))*s),0.,1.,weights)
 end
-@everywhere function dsmudgamma1(e,rho,taus,T,weights,weights_gL,P,J,nu)
-  out=Array{Float64}(undef,J*(P+1))
+@everywhere function dsmudgamma1(
+  e::Float64,
+  rho::Float64,
+  taus::Array{Float64,1},
+  T::Array{Polynomial,1},
+  weights::Tuple{Array{Float64,1},Array{Float64,1}},
+  weights_gL::Tuple{Array{Float64,1},Array{Float64,1}},
+  P::Int64,
+  J::Int64,
+  nu::Int64
+)
+  out=Array{Float64,1}(undef,J*(P+1))
   for zeta in 0:J-1
     for n in 0:P
       out[zeta*(P+1)+n+1]=-4*rho*e*integrate_laguerre(s->s*(s+1)*B1(s,zeta,n,taus,T,weights,nu),2*sqrt(abs(e)),weights_gL)
@@ -531,7 +765,17 @@ end
   end
   return out
 end
-@everywhere function dsmudgamma2(e,rho,taus,T,weights,weights_gL,P,J,nu)
+@everywhere function dsmudgamma2(
+  e::Float64,
+  rho::Float64,
+  taus::Array{Float64,1},
+  T::Array{Polynomial,1},
+  weights::Tuple{Array{Float64,1},Array{Float64,1}},
+  weights_gL::Tuple{Array{Float64,1},Array{Float64,1}},
+  P::Int64,
+  J::Int64,
+  nu::Int64
+)
   out=Array{Float64,2}(undef,J*(P+1),J*(P+1))
   for zeta in 0:J-1
     for n in 0:P
@@ -546,25 +790,41 @@ end
 end
 
 # \bar gamma's
-@everywhere function gammabar0(weights)
-  return 4*pi*integrate_legendre(s->s^2*B0(s-1),0,1,weights)
+@everywhere function gammabar0(
+  weights::Tuple{Array{Float64,1},Array{Float64,1}}
+)
+  return 4*pi*integrate_legendre(s->s^2*B0(s-1),0.,1.,weights)
 end
-@everywhere function gammabar1(taus,T,weights,P,J,nu)
-  out=Array{Float64}(undef,J*(P+1))
+@everywhere function gammabar1(
+  taus::Array{Float64,1},
+  T::Array{Polynomial,1},
+  weights::Tuple{Array{Float64,1},Array{Float64,1}},
+  P::Int64,
+  J::Int64,
+  nu::Int64
+)
+  out=Array{Float64,1}(undef,J*(P+1))
   for zeta in 0:J-1
     for n in 0:P
-      out[zeta*(P+1)+n+1]=4*pi*integrate_legendre(s->s^2*B1(s-1,zeta,n,taus,T,weights,nu),0,1,weights)
+      out[zeta*(P+1)+n+1]=4*pi*integrate_legendre(s->s^2*B1(s-1,zeta,n,taus,T,weights,nu),0.,1.,weights)
     end
   end
   return out
 end
-@everywhere function gammabar2(taus,T,weights,P,J,nu)
+@everywhere function gammabar2(
+  taus::Array{Float64,1},
+  T::Array{Polynomial,1},
+  weights::Tuple{Array{Float64,1},Array{Float64,1}},
+  P::Int64,
+  J::Int64,
+  nu::Int64
+)
   out=Array{Float64,2}(undef,J*(P+1),J*(P+1))
   for zeta in 0:J-1
     for n in 0:P
       for zetap in 0:J-1
 	for m in 0:P
-	  out[zeta*(P+1)+n+1,zetap*(P+1)+m+1]=4*pi*integrate_legendre(s->s^2*B2(s-1,zeta,n,zetap,m,taus,T,weights,nu),0,1,weights)
+	  out[zeta*(P+1)+n+1,zetap*(P+1)+m+1]=4*pi*integrate_legendre(s->s^2*B2(s-1,zeta,n,zetap,m,taus,T,weights,nu),0.,1.,weights)
 	end
       end
     end
diff --git a/src/simpleq-iteration.jl b/src/simpleq-iteration.jl
index 98977b8..4c4de07 100644
--- a/src/simpleq-iteration.jl
+++ b/src/simpleq-iteration.jl
@@ -1,4 +1,4 @@
-## Copyright 2021 Ian Jauslin
+## 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.
@@ -13,7 +13,12 @@
 ## limitations under the License.
 
 # compute rho(e) using the iteration
-function simpleq_iteration_rho_e(es,order,v,maxiter)
+function simpleq_iteration_rho_e(
+  es::Array{Float64,1},
+  order::Int64,
+  v::Function,
+  maxiter::Int64
+)
   for j in 1:length(es)
     (u,rho)=simpleq_iteration_hatun(es[j],order,v,maxiter)
     @printf("% .15e % .15e\n",es[j],real(rho[maxiter+1]))
@@ -21,7 +26,15 @@ function simpleq_iteration_rho_e(es,order,v,maxiter)
 end
 
 # compute u(x) using the iteration and print at every step
-function simpleq_iteration_ux(order,e,v,maxiter,xmin,xmax,nx)
+function simpleq_iteration_ux(
+  order::Int64,
+  e::Float64,
+  v::Function,
+  maxiter::Int64,
+  xmin::Float64,
+  xmax::Float64,
+  nx::Int64
+)
   (u,rho)=simpleq_iteration_hatun(e,order,v,maxiter)
 
   weights=gausslegendre(order)
@@ -37,7 +50,12 @@ end
 
 
 # \hat u_n
-function simpleq_iteration_hatun(e,order,v,maxiter)
+function simpleq_iteration_hatun(
+  e::Float64,
+  order::Int64,
+  v::Function,
+  maxiter::Int64
+)
   # gauss legendre weights
   weights=gausslegendre(order)
 
@@ -46,7 +64,7 @@ function simpleq_iteration_hatun(e,order,v,maxiter)
   (Eta,Eta0)=easyeq_init_H(weights,v)
 
   # init u and rho
-  u=Array{Array{Float64}}(undef,maxiter+1)
+  u=Array{Array{Float64,1},1}(undef,maxiter+1)
   u[1]=zeros(Float64,order)
   rho=zeros(Float64,maxiter+1)
 
@@ -60,7 +78,11 @@ function simpleq_iteration_hatun(e,order,v,maxiter)
 end
 
 # A
-function simpleq_iteration_A(e,weights,Eta)
+function simpleq_iteration_A(
+  e::Float64,
+  weights::Tuple{Array{Float64,1},Array{Float64,1}},
+  Eta::Array{Float64,1}
+)
   N=length(weights[1])
   out=zeros(Float64,N,N)
   for i in 1:N
@@ -75,7 +97,12 @@ function simpleq_iteration_A(e,weights,Eta)
 end
 
 # b
-function simpleq_iteration_b(u,e,rho,V)
+function simpleq_iteration_b(
+  u::Array{Float64,1},
+  e::Float64,
+  rho::Float64,
+  V::Array{Float64,1}
+)
   out=zeros(Float64,length(V))
   for i in 1:length(V)
     out[i]=V[i]+2*e*rho*u[i]^2
@@ -84,7 +111,13 @@ function simpleq_iteration_b(u,e,rho,V)
 end
 
 # rho_n
-function simpleq_iteration_rhon(u,e,weights,V0,Eta0)
+function simpleq_iteration_rhon(
+  u::Array{Float64,1},
+  e::Float64,
+  weights::Tuple{Array{Float64,1},Array{Float64,1}},
+  V0::Float64,
+  Eta0::Float64
+)
   S=V0
   for i in 1:length(weights[1])
     y=(weights[1][i]+1)/2
diff --git a/src/tools.jl b/src/tools.jl
index 0d3dc7f..5d1678d 100644
--- a/src/tools.jl
+++ b/src/tools.jl
@@ -1,4 +1,4 @@
-## Copyright 2021 Ian Jauslin
+## 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.
@@ -13,7 +13,9 @@
 ## limitations under the License.
 
 # \Phi(x)=2*(1-sqrt(1-x))/x
-@everywhere function Phi(x)
+@everywhere function Phi(
+  x::Float64
+)
   if abs(x)>1e-5
     return 2*(1-sqrt(abs(1-x)))/x
   else
@@ -21,7 +23,9 @@
   end
 end
 # \partial\Phi
-@everywhere function dPhi(x)
+@everywhere function dPhi(
+  x::Float64
+)
   #if abs(x-1)<1e-5
   #  @printf(stderr,"warning: dPhi is singular at 1, and evaluating it at (% .8e+i% .8e)\n",real(x),imag(x))
   #end
@@ -33,17 +37,25 @@ end
 end
 
 # apply Phi to every element of a vector
-@everywhere function dotPhi(v)
+@everywhere function dotPhi(
+  v::Array{Float64,1}
+)
   out=zeros(Float64,length(v))
   for i in 1:length(v)
     out[i]=Phi(v[i])
   end
   return out
 end
-@everywhere function dotdPhi(v)
+@everywhere function dotdPhi(
+  v::Array{Float64,1}
+)
   out=zeros(Float64,length(v))
   for i in 1:length(v)
     out[i]=dPhi(v[i])
   end
   return out
 end
+
+@everywhere function sinc(x::Float64)
+  return (x == 0 ? 1 : sin(x)/x)
+end