module Math.EllipticIntegrals.Elliptic
  where
import Math.EllipticIntegrals.Carlson
import Data.Complex
import Math.EllipticIntegrals.Internal

-- | Elliptic integral of the first kind.
ellipticF' :: 
     Double -- ^ bound on the relative error passed to `carlsonRF'`
  -> Cplx   -- ^ amplitude
  -> Cplx   -- ^ parameter
  -> Cplx
ellipticF' :: Double -> Cplx -> Cplx -> Cplx
ellipticF' Double
err Cplx
phi Cplx
m
  | Cplx
phi forall a. Eq a => a -> a -> Bool
== Cplx
0 =
    Double -> Cplx
toCplx Double
0
  | Cplx
m forall a. Eq a => a -> a -> Bool
== Cplx
1 Bool -> Bool -> Bool
&& forall a. Num a => a -> a
abs(forall a. Complex a -> a
realPart Cplx
phi) forall a. Eq a => a -> a -> Bool
== forall a. Floating a => a
piforall a. Fractional a => a -> a -> a
/Double
2 =
    Double -> Cplx
toCplx (Double
0forall a. Fractional a => a -> a -> a
/Double
0)
  | Cplx
m forall a. Eq a => a -> a -> Bool
== Cplx
1 Bool -> Bool -> Bool
&& forall a. Num a => a -> a
abs(forall a. Complex a -> a
realPart Cplx
phi) forall a. Ord a => a -> a -> Bool
< forall a. Floating a => a
piforall a. Fractional a => a -> a -> a
/Double
2 =
    forall a. Floating a => a -> a
atanh(forall a. Floating a => a -> a
sin Cplx
phi)
  | forall a. Num a => a -> a
abs(forall a. Complex a -> a
realPart Cplx
phi) forall a. Ord a => a -> a -> Bool
<= forall a. Floating a => a
piforall a. Fractional a => a -> a -> a
/Double
2 =
    if Cplx
m forall a. Eq a => a -> a -> Bool
== Cplx
0
      then
        Cplx
phi
      else
        let sine :: Cplx
sine = forall a. Floating a => a -> a
sin Cplx
phi in
        let sine2 :: Cplx
sine2 = Cplx
sineforall a. Num a => a -> a -> a
*Cplx
sine in
        let (Cplx
cosine2, Cplx
oneminusmsine2) = (Cplx
1 forall a. Num a => a -> a -> a
- Cplx
sine2, Cplx
1 forall a. Num a => a -> a -> a
- Cplx
mforall a. Num a => a -> a -> a
*Cplx
sine2) in
        Cplx
sine forall a. Num a => a -> a -> a
* Double -> Cplx -> Cplx -> Cplx -> Cplx
carlsonRF' Double
err Cplx
cosine2 Cplx
oneminusmsine2 Cplx
1
  | Bool
otherwise =
    let (Cplx
phi', Int
k) = Cplx -> (Cplx, Int)
getPhiK Cplx
phi in
    Cplx
2 forall a. Num a => a -> a -> a
* forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
k forall a. Num a => a -> a -> a
* Double -> Cplx -> Cplx -> Cplx
ellipticF' Double
err (forall a. Floating a => a
piforall a. Fractional a => a -> a -> a
/Cplx
2) Cplx
m forall a. Num a => a -> a -> a
+ Double -> Cplx -> Cplx -> Cplx
ellipticF' Double
err Cplx
phi' Cplx
m

-- | Elliptic integral of the first kind.
ellipticF :: 
     Cplx -- ^ amplitude
  -> Cplx -- ^ parameter
  -> Cplx
ellipticF :: Cplx -> Cplx -> Cplx
ellipticF = Double -> Cplx -> Cplx -> Cplx
ellipticF' Double
1e-15

-- | Elliptic integral of the second kind.
ellipticE' :: 
     Double -- ^ bound on the relative error passed to `carlsonRF'` and `carlsonRD'` 
  -> Cplx   -- ^ amplitude
  -> Cplx   -- ^ parameter
  -> Cplx
ellipticE' :: Double -> Cplx -> Cplx -> Cplx
ellipticE' Double
err Cplx
phi Cplx
m
  | Cplx
phi forall a. Eq a => a -> a -> Bool
== Cplx
0 =
    Double -> Cplx
toCplx Double
0
  | forall a. Num a => a -> a
abs(forall a. Complex a -> a
realPart Cplx
phi) forall a. Ord a => a -> a -> Bool
<= forall a. Floating a => a
piforall a. Fractional a => a -> a -> a
/Double
2 =
    case Cplx
m of
      Cplx
0 -> Cplx
phi
      Cplx
1 -> forall a. Floating a => a -> a
sin Cplx
phi
      Cplx
_ ->
        let sine :: Cplx
sine = forall a. Floating a => a -> a
sin Cplx
phi in
        let sine2 :: Cplx
sine2 = Cplx
sineforall a. Num a => a -> a -> a
*Cplx
sine in
        let (Cplx
cosine2, Cplx
oneminusmsine2) = (Cplx
1 forall a. Num a => a -> a -> a
- Cplx
sine2, Cplx
1 forall a. Num a => a -> a -> a
- Cplx
mforall a. Num a => a -> a -> a
*Cplx
sine2) in
        Cplx
sine forall a. Num a => a -> a -> a
* (Double -> Cplx -> Cplx -> Cplx -> Cplx
carlsonRF' Double
err Cplx
cosine2 Cplx
oneminusmsine2 Cplx
1 forall a. Num a => a -> a -> a
-
          Cplx
m forall a. Num a => a -> a -> a
* Cplx
sine2 forall a. Fractional a => a -> a -> a
/ Cplx
3 forall a. Num a => a -> a -> a
* Double -> Cplx -> Cplx -> Cplx -> Cplx
carlsonRD' Double
err Cplx
cosine2 Cplx
oneminusmsine2 Cplx
1)
  | Bool
otherwise =
    let (Cplx
phi', Int
k) = Cplx -> (Cplx, Int)
getPhiK Cplx
phi in
    Cplx
2 forall a. Num a => a -> a -> a
* forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
k forall a. Num a => a -> a -> a
* Double -> Cplx -> Cplx -> Cplx
ellipticE' Double
err (forall a. Floating a => a
piforall a. Fractional a => a -> a -> a
/Cplx
2) Cplx
m forall a. Num a => a -> a -> a
+ Double -> Cplx -> Cplx -> Cplx
ellipticE' Double
err Cplx
phi' Cplx
m

-- | Elliptic integral of the second kind.
ellipticE :: 
     Cplx -- ^ amplitude
  -> Cplx -- ^ parameter
  -> Cplx
ellipticE :: Cplx -> Cplx -> Cplx
ellipticE = Double -> Cplx -> Cplx -> Cplx
ellipticE' Double
1e-15

-- | Elliptic integral of the third kind.
ellipticPI' :: 
     Double -- ^ bound on the relative error passed to `carlsonRF'` and `carlsonRJ'` 
  -> Cplx   -- ^ amplitude
  -> Cplx   -- ^ characteristic
  -> Cplx   -- ^ parameter
  -> Cplx
ellipticPI' :: Double -> Cplx -> Cplx -> Cplx -> Cplx
ellipticPI' Double
err Cplx
phi Cplx
n Cplx
m
  | Cplx
phi forall a. Eq a => a -> a -> Bool
== Cplx
0 =
    Double -> Cplx
toCplx Double
0
  | Cplx
phi forall a. Eq a => a -> a -> Bool
== forall a. Floating a => a
piforall a. Fractional a => a -> a -> a
/Cplx
2 Bool -> Bool -> Bool
&& Cplx
n forall a. Eq a => a -> a -> Bool
== Cplx
1 =
    Cplx
0forall a. Fractional a => a -> a -> a
/Cplx
0
  | Cplx
phi forall a. Eq a => a -> a -> Bool
== forall a. Floating a => a
piforall a. Fractional a => a -> a -> a
/Cplx
2 Bool -> Bool -> Bool
&& Cplx
m forall a. Eq a => a -> a -> Bool
== Cplx
0 =
    forall a. Floating a => a
piforall a. Fractional a => a -> a -> a
/Cplx
2forall a. Fractional a => a -> a -> a
/forall a. Floating a => a -> a
sqrt(Cplx
1forall a. Num a => a -> a -> a
-Cplx
n)
  | Cplx
phi forall a. Eq a => a -> a -> Bool
== forall a. Floating a => a
piforall a. Fractional a => a -> a -> a
/Cplx
2 Bool -> Bool -> Bool
&& Cplx
m forall a. Eq a => a -> a -> Bool
== Cplx
n =
    Double -> Cplx -> Cplx -> Cplx
ellipticE' Double
err (forall a. Floating a => a
piforall a. Fractional a => a -> a -> a
/Cplx
2) Cplx
m forall a. Fractional a => a -> a -> a
/ (Cplx
1forall a. Num a => a -> a -> a
-Cplx
m)
  | Cplx
phi forall a. Eq a => a -> a -> Bool
== forall a. Floating a => a
piforall a. Fractional a => a -> a -> a
/Cplx
2 Bool -> Bool -> Bool
&& Cplx
n forall a. Eq a => a -> a -> Bool
== Cplx
0 =
    Double -> Cplx -> Cplx -> Cplx
ellipticF' Double
err (forall a. Floating a => a
piforall a. Fractional a => a -> a -> a
/Cplx
2) Cplx
m
  | forall a. Num a => a -> a
abs(forall a. Complex a -> a
realPart Cplx
phi) forall a. Ord a => a -> a -> Bool
<= forall a. Floating a => a
piforall a. Fractional a => a -> a -> a
/Double
2 =
    let sine :: Cplx
sine = forall a. Floating a => a -> a
sin Cplx
phi in
    let sine2 :: Cplx
sine2 = Cplx
sineforall a. Num a => a -> a -> a
*Cplx
sine in
    let (Cplx
cosine2, Cplx
oneminusmsine2) = (Cplx
1 forall a. Num a => a -> a -> a
- Cplx
sine2, Cplx
1 forall a. Num a => a -> a -> a
- Cplx
mforall a. Num a => a -> a -> a
*Cplx
sine2) in
    Cplx
sine forall a. Num a => a -> a -> a
* (Double -> Cplx -> Cplx -> Cplx -> Cplx
carlsonRF' Double
err Cplx
cosine2 Cplx
oneminusmsine2 Cplx
1 forall a. Num a => a -> a -> a
+
      Cplx
n forall a. Num a => a -> a -> a
* Cplx
sine2 forall a. Fractional a => a -> a -> a
/ Cplx
3 forall a. Num a => a -> a -> a
* Double -> Cplx -> Cplx -> Cplx -> Cplx -> Cplx
carlsonRJ' Double
err Cplx
cosine2 Cplx
oneminusmsine2 Cplx
1 (Cplx
1forall a. Num a => a -> a -> a
-Cplx
nforall a. Num a => a -> a -> a
*Cplx
sine2))
  | Bool
otherwise =
    let (Cplx
phi', Int
k) = Cplx -> (Cplx, Int)
getPhiK Cplx
phi in
    Cplx
2 forall a. Num a => a -> a -> a
* forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
k forall a. Num a => a -> a -> a
* Double -> Cplx -> Cplx -> Cplx -> Cplx
ellipticPI' Double
err (forall a. Floating a => a
piforall a. Fractional a => a -> a -> a
/Cplx
2) Cplx
n Cplx
m forall a. Num a => a -> a -> a
+ Double -> Cplx -> Cplx -> Cplx -> Cplx
ellipticPI' Double
err Cplx
phi' Cplx
n Cplx
m

-- | Elliptic integral of the third kind.
ellipticPI ::
     Cplx -- ^ amplitude
  -> Cplx -- ^ characteristic
  -> Cplx -- ^ parameter
  -> Cplx
ellipticPI :: Cplx -> Cplx -> Cplx -> Cplx
ellipticPI = Double -> Cplx -> Cplx -> Cplx -> Cplx
ellipticPI' Double
1e-15

-- | Jacobi zeta function.
jacobiZeta' ::
     Double -- ^ bound on the relative error passed to `ellipticF'` and `ellipticE'` 
  -> Cplx   -- ^ amplitude
  -> Cplx   -- ^ parameter
  -> Cplx
jacobiZeta' :: Double -> Cplx -> Cplx -> Cplx
jacobiZeta' Double
err Cplx
phi Cplx
m =
  if Cplx
m forall a. Eq a => a -> a -> Bool
== Cplx
1
    then
      if forall a. Num a => a -> a
abs(forall a. Complex a -> a
realPart Cplx
phi) forall a. Ord a => a -> a -> Bool
<= forall a. Floating a => a
piforall a. Fractional a => a -> a -> a
/Double
2
        then forall a. Floating a => a -> a
sin Cplx
phi
        else let (Cplx
phi',Int
_) = Cplx -> (Cplx, Int)
getPhiK Cplx
phi in forall a. Floating a => a -> a
sin Cplx
phi'
    else
      Double -> Cplx -> Cplx -> Cplx
ellipticE' Double
err Cplx
phi Cplx
m forall a. Num a => a -> a -> a
-
        Double -> Cplx -> Cplx -> Cplx
ellipticE' Double
err (forall a. Floating a => a
piforall a. Fractional a => a -> a -> a
/Cplx
2) Cplx
m forall a. Fractional a => a -> a -> a
/ Double -> Cplx -> Cplx -> Cplx
ellipticF' Double
err (forall a. Floating a => a
piforall a. Fractional a => a -> a -> a
/Cplx
2) Cplx
m forall a. Num a => a -> a -> a
*
        Double -> Cplx -> Cplx -> Cplx
ellipticF' Double
err Cplx
phi Cplx
m

-- | Jacobi zeta function.
jacobiZeta ::
     Cplx -- ^ amplitude
  -> Cplx -- ^ parameter
  -> Cplx
jacobiZeta :: Cplx -> Cplx -> Cplx
jacobiZeta = Double -> Cplx -> Cplx -> Cplx
jacobiZeta' Double
1e-15