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

# exponential potential in 3 dimensions
@everywhere function v_exp(k,a)
  return 8*pi/(1+k^2)^2*a
end
@everywhere function a0_exp(a)
  if a>0.
    return log(a)+2*MathConstants.eulergamma+2*besselk(0,2*sqrt(a))/besseli(0,2*sqrt(a))
  elseif a<0.
    return log(-a)+2*MathConstants.eulergamma-pi*bessely(0,2*sqrt(-a))/besselj(0,2*sqrt(-a))
  else
    return 0.
  end
end

# exp(-x)-a*exp(-b*x) in 3 dimensions
@everywhere function v_expcry(k,a,b)
  return 8*pi*((1+k^2)^(-2)-a*b*(b^2+k^2)^(-2))
end
@everywhere function a0_expcry(a,b)
  return 1.21751642717932720441274114683413710125487579284827462 #ish
end

# x^2*exp(-|x|) in 3 dimensions
@everywhere function v_npt(k)
  return 96*pi*(1-k^2)/(1+k^2)^4
end
@everywhere function a0_npt()
  return 1. #ish
end

# 1/(1+x^4/4) potential in 3 dimensions
@everywhere function v_alg(k)
  if(k==0)
    return 4*pi^2
  else
    return 4*pi^2*exp(-k)*sin(k)/k
  end
end
a0_alg=1. #ish

# (1+a x^4)/(1+x^2)^4 potential in 3 dimensions
@everywhere function v_algwell(k)
  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)
  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)
  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)
  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)
  return b #ish
end