## Copyright 2021 Ian Jauslin
## 
## Licensed under the Apache License, Version 2.0 (the "License");
## you may not use this file except in compliance with the License.
## You may obtain a copy of the License at
## 
##     http://www.apache.org/licenses/LICENSE-2.0
## 
## Unless required by applicable law or agreed to in writing, software
## distributed under the License is distributed on an "AS IS" BASIS,
## WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
## See the License for the specific language governing permissions and
## limitations under the License.

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

  # 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))
  for j in 1:length(rhos)
    u0s[j]=Array{Float64}(undef,N*J)
    for i in 1:N*J
      u0s[j][i]=1/(1+r[i]^2)^2
    end
    e0s[j]=pi*rhos[j]
  end


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

  count=0
  # for each worker
  @sync for p in 1:nw
    # for each task
    @async for j in work[p]
      count=count+1
      if count>=nw
        progress(count,length(rhos),10000)
      end
      # run the task
      (us[j],es[j],err[j])=remotecall_fetch(simpleq_hardcore_hatu,workers()[p],u0s[j],e0s[j],rhos[j],r,taus,T,weights,weights_gL,P,N,J,maxiter,tolerance)
    end
  end

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

# compute u(x)
function simpleq_hardcore_ux(rho,taus,P,N,J,maxiter,tolerance)
  # initialize vectors
  (weights,weights_gL,r,T)=simpleq_hardcore_init(taus,P,N,J)

  # initial guess
  u0=Array{Float64}(undef,N*J)
  for i in 1:N*J
    u0[i]=1/(1+r[i]^2)^2
  end
  e0=pi*rho

  (u,e,err)=simpleq_hardcore_hatu(u0,e0,rho,r,taus,T,weights,weights_gL,P,N,J,maxiter,tolerance)

  for i in 1:N*J
    @printf("% .15e % .15e\n",r[i],u[i])
  end
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

  # 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))
  for j in 1:length(rhos)
    u0s[j]=Array{Float64}(undef,N*J)
    for i in 1:N*J
      u0s[j][i]=1/(1+r[i]^2)^2
    end
    e0s[j]=pi*rhos[j]
  end


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

  count=0
  # for each worker
  @sync for p in 1:nw
    # for each task
    @async for j in work[p]
      count=count+1
      if count>=nw
        progress(count,length(rhos),10000)
      end
      # run the task
      (us[j],es[j],err[j])=remotecall_fetch(simpleq_hardcore_hatu,workers()[p],u0s[j],e0s[j],rhos[j],r,taus,T,weights,weights_gL,P,N,J,maxiter,tolerance)
    end
  end

  for j in 1:length(rhos)
    du=-inv(simpleq_hardcore_DXi(us[j],es[j],rhos[j],r,taus,T,weights,weights_gL,P,N,J))*simpleq_hardcore_dXidmu(us[j],es[j],rhos[j],r,taus,T,weights,weights_gL,P,N,J)
    @printf("% .15e % .15e % .15e\n",rhos[j],du[N*J+1],err[j])
  end
end


# initialize computation
@everywhere function simpleq_hardcore_init(taus,P,N,J)
  # Gauss-Legendre weights
  weights=gausslegendre(N)
  weights_gL=gausslaguerre(N)

  # r
  r=Array{Float64}(undef,J*N)
  for zeta in 0:J-1
    for j in 1:N
      xj=weights[1][j]
      # set kj
      r[zeta*N+j]=(2+(taus[zeta+2]-taus[zeta+1])*sin(pi*xj/2)-(taus[zeta+2]+taus[zeta+1]))/(2-(taus[zeta+2]-taus[zeta+1])*sin(pi*xj/2)+taus[zeta+2]+taus[zeta+1])
    end
  end

  # Chebyshev polynomials
  T=chebyshev_polynomials(P)

  return (weights,weights_gL,r,T)
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)
  # init
  vec=Array{Float64}(undef,J*N+1)
  for i in 1:J*N
    vec[i]=u0[i]
  end
  vec[J*N+1]=e0

  # quantify relative error
  error=Inf

  # run Newton algorithm
  for i in 1:maxiter-1
    u=vec[1:J*N]
    e=vec[J*N+1]
    new=vec-inv(simpleq_hardcore_DXi(u,e,rho,r,taus,T,weights,weights_gL,P,N,J))*simpleq_hardcore_Xi(u,e,rho,r,taus,T,weights,weights_gL,P,N,J)

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

    vec=new
  end

  return(vec[1:J*N],vec[J*N+1],error)
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)
  FU=chebyshev(u,taus,weights,P,N,J,4)

  #D's
  d0=D0(e,rho,r,weights,N,J)
  d1=D1(e,rho,r,taus,T,weights,weights_gL,P,N,J,4)
  d2=D2(e,rho,r,taus,T,weights,weights_gL,P,N,J,4)

  # u
  for i in 1:J*N
    out[i]=d0[i]+dot(FU,d1[i])+dot(FU,d2[i]*FU)-u[i]
  end
  # e
  out[J*N+1]=-e+2*pi*rho*
    ((1+2*sqrt(abs(e)))-gamma0(e,rho,weights)-dot(FU,gamma1(e,rho,taus,T,weights,weights_gL,P,J,4))-dot(FU,gamma2(e,rho,taus,T,weights,weights_gL,P,J,4)*FU))/
    (1-8/3*pi*rho+rho^2*(gammabar0(weights)+dot(FU,gammabar1(taus,T,weights,P,J,4))+dot(FU,gammabar2(taus,T,weights,P,J,4)*FU)))

  return out
end
# DXi
@everywhere function simpleq_hardcore_DXi(u,e,rho,r,taus,T,weights,weights_gL,P,N,J)
  out=Array{Float64,2}(undef,J*N+1,J*N+1)
  FU=chebyshev(u,taus,weights,P,N,J,4)

  #D's
  d0=D0(e,rho,r,weights,N,J)
  d1=D1(e,rho,r,taus,T,weights,weights_gL,P,N,J,4)
  d2=D2(e,rho,r,taus,T,weights,weights_gL,P,N,J,4)
  dsed0=dseD0(e,rho,r,weights,N,J)
  dsed1=dseD1(e,rho,r,taus,T,weights,weights_gL,P,N,J,4)
  dsed2=dseD2(e,rho,r,taus,T,weights,weights_gL,P,N,J,4)

  # denominator of e
  denom=1-8/3*pi*rho+rho^2*(gammabar0(weights)+dot(FU,gammabar1(taus,T,weights,P,J,4))+dot(FU,gammabar2(taus,T,weights,P,J,4)*FU))

  for zetapp in 0:J-1
    for n in 1:N
      one=zeros(Int64,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)
      end
      # de/du
      out[J*N+1,zetapp*N+n]=2*pi*rho*
	  (-dot(Fone,gamma1(e,rho,taus,T,weights,weights_gL,P,J,4))-2*dot(FU,gamma2(e,rho,taus,T,weights,weights_gL,P,J,4)*Fone))/denom-
        2*pi*rho*
	  ((1+2*sqrt(abs(e)))-gamma0(e,rho,weights)-dot(FU,gamma1(e,rho,taus,T,weights,weights_gL,P,J,4))-dot(FU,gamma2(e,rho,taus,T,weights,weights_gL,P,J,4)*FU))*
	  rho^2*(dot(Fone,gammabar1(taus,T,weights,P,J,4))+2*dot(FU,gammabar2(taus,T,weights,P,J,4)*Fone))/denom^2
    end
  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)

  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)
  FU=chebyshev(u,taus,weights,P,N,J,4)

  #D's
  dsmud0=dsmuD0(e,rho,r,weights,N,J)
  dsmud1=dsmuD1(e,rho,r,taus,T,weights,weights_gL,P,N,J,4)
  dsmud2=dsmuD2(e,rho,r,taus,T,weights,weights_gL,P,N,J,4)

  # u
  for i in 1:J*N
    out[i]=(dsmud0[i]+dot(FU,dsmud1[i])+dot(FU,dsmud2[i]*FU))/(4*sqrt(abs(e)))
  end
  # e
  out[J*N+1]=2*pi*rho*(2/3+1/(2*sqrt(abs(e)))-
    (dsmudgamma0(e,rho,weights)+dot(FU,dsmudgamma1(e,rho,taus,T,weights,weights_gL,P,J,4))+dot(FU,dsmudgamma2(e,rho,taus,T,weights,weights_gL,P,J,4)*FU))/(4*sqrt(abs(e)))
  )/(1-8/3*pi*rho+rho^2*(gammabar0(weights)+dot(FU,gammabar1(taus,T,weights,P,J,4))+dot(FU,gammabar2(taus,T,weights,P,J,4)*FU)))

  return out
end

# B's
@everywhere function B0(r)
  return pi/12*(r-1)^2*(r+5)
end
@everywhere function B1(r,zeta,n,taus,T,weights,nu)
  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)
  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 :
      integrate_legendre(s->
	1/(1+s)^(3-nu)*T[m+1]((2*s-(taus[zetap+1]+taus[zetap+2]))/(taus[zetap+2]-taus[zetap+1])),max(taus[zetap+1]
      ,alpham(1+r,tau)),min(taus[zetap+2],alphap(abs(r-(1-tau)/(1+tau))-2*tau/(1+tau),tau)),weights)
    )
  ,taus[zeta+1],taus[zeta+2],weights)
end

# D's
@everywhere function D0(e,rho,r,weights,N,J)
  out=Array{Float64}(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)+
      (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)
      )
  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)
  for i in 1:J*N
    out[i]=Array{Float64}(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)+
	  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)
      end
    end
  end

  return out
end
@everywhere function D2(e,rho,r,taus,T,weights,weights_gL,P,N,J,nu)
  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)
    for zeta in 0:J-1
      for n in 0:P
	for zetap in 0:J-1
	  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)+
	      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)
	  end
	end
      end
    end
  end
  return out
end

# dD/d sqrt(abs(e))'s
@everywhere function dseD0(e,rho,r,weights,N,J)
  out=Array{Float64}(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)+
      (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)
      )
  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)
  for i in 1:J*N
    out[i]=Array{Float64}(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)+
	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)
      end
    end
  end
  return out
end
@everywhere function dseD2(e,rho,r,taus,T,weights,weights_gL,P,N,J,nu)
  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)
    for zeta in 0:J-1
      for n in 0:P
	for zetap in 0:J-1
	  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)+
	    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)
	  end
	end
      end
    end
  end
  return out
end

# dD/d sqrt(abs(e+mu/2))'s
@everywhere function dsmuD0(e,rho,r,weights,N,J)
  out=Array{Float64}(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)+
      (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)
      )
  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)
  for i in 1:J*N
    out[i]=Array{Float64}(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)+
	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)
      end
    end
  end
  return out
end
@everywhere function dsmuD2(e,rho,r,taus,T,weights,weights_gL,P,N,J,nu)
  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)
    for zeta in 0:J-1
      for n in 0:P
	for zetap in 0:J-1
	  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)+
	    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)
	  end
	end
      end
    end
  end
  return out
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)
end
@everywhere function gamma1(e,rho,taus,T,weights,weights_gL,P,J,nu)
  out=Array{Float64}(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)
    end
  end
  return out
end
@everywhere function gamma2(e,rho,taus,T,weights,weights_gL,P,J,nu)
  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]=2*rho*e*integrate_laguerre(s->(s+1)*B2(s,zeta,n,zetap,m,taus,T,weights,nu),2*sqrt(abs(e)),weights_gL)
	end
      end
    end
  end
  return out
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)
end
@everywhere function dsedgamma1(e,rho,taus,T,weights,weights_gL,P,J,nu)
  out=Array{Float64}(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)
    end
  end
  return out
end
@everywhere function dsedgamma2(e,rho,taus,T,weights,weights_gL,P,J,nu)
  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*rho*e*integrate_laguerre(s->(s+1)*B2(s,zeta,n,zetap,m,taus,T,weights,nu)*(1-sqrt(abs(e))*s),2*sqrt(abs(e)),weights_gL)
	end
      end
    end
  end
  return out
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)
end
@everywhere function dsmudgamma1(e,rho,taus,T,weights,weights_gL,P,J,nu)
  out=Array{Float64}(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)
    end
  end
  return out
end
@everywhere function dsmudgamma2(e,rho,taus,T,weights,weights_gL,P,J,nu)
  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*rho*e*integrate_laguerre(s->s*(s+1)*B2(s,zeta,n,zetap,m,taus,T,weights,nu),2*sqrt(abs(e)),weights_gL)
	end
      end
    end
  end
  return out
end

# \bar gamma's
@everywhere function gammabar0(weights)
  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))
  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)
    end
  end
  return out
end
@everywhere function gammabar2(taus,T,weights,P,J,nu)
  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)
	end
      end
    end
  end
  return out
end