module Math.Weierstrass
    ( halfPeriods,
      ellipticInvariants,
      weierstrassP,
      weierstrassPdash,
      weierstrassPinv,
      weierstrassSigma,
      weierstrassZeta
    ) where
import           Data.Complex           ( Complex(..) )
import           Internal               ( (%^%) )
import           Math.Eisenstein        ( eisensteinE4, 
                                          eisensteinE6, 
                                          kleinJinv, 
                                          jtheta1DashDashDash0 ) 
import           Math.JacobiTheta       ( jtheta2, 
                                          jtheta3, 
                                          jtheta1, 
                                          jtheta4,
                                          jtheta1Dash )
import           Math.Gamma             ( gamma )
import           Math.EllipticIntegrals ( carlsonRF' )



i_ :: Complex Double
i_ :: Complex Double
i_ = Double
0.0 forall a. a -> a -> Complex a
:+ Double
1.0

eisensteinG4 :: Complex Double -> Complex Double
eisensteinG4 :: Complex Double -> Complex Double
eisensteinG4 Complex Double
tau = forall a. Floating a => a
pi Complex Double -> Int -> Complex Double
%^% Int
4 forall a. Fractional a => a -> a -> a
/ Complex Double
45 forall a. Num a => a -> a -> a
* Complex Double -> Complex Double
eisensteinE4 Complex Double
tau

eisensteinG6_over_eisensteinG4 :: Complex Double -> Complex Double
eisensteinG6_over_eisensteinG4 :: Complex Double -> Complex Double
eisensteinG6_over_eisensteinG4 Complex Double
tau = 
  Complex Double
2 forall a. Num a => a -> a -> a
* forall a. Floating a => a
pi forall a. Num a => a -> a -> a
* forall a. Floating a => a
pi forall a. Fractional a => a -> a -> a
/ Complex Double
21 forall a. Num a => a -> a -> a
* Complex Double -> Complex Double
eisensteinE6 Complex Double
tau forall a. Fractional a => a -> a -> a
/ Complex Double -> Complex Double
eisensteinE4 Complex Double
tau

omega1_and_tau :: 
  Complex Double -> Complex Double -> (Complex Double, Complex Double)
omega1_and_tau :: Complex Double
-> Complex Double -> (Complex Double, Complex Double)
omega1_and_tau Complex Double
g2 Complex Double
g3 = (Complex Double
omega1, Complex Double
tau)
  where
    (Complex Double
omega1, Complex Double
tau) 
      | Complex Double
g2 forall a. Eq a => a -> a -> Bool
== Complex Double
0 = 
        (
          forall a. Gamma a => a -> a
gamma (Complex Double
1forall a. Fractional a => a -> a -> a
/Complex Double
3) Complex Double -> Int -> Complex Double
%^% Int
3 forall a. Fractional a => a -> a -> a
/ (Complex Double
4 forall a. Num a => a -> a -> a
* forall a. Floating a => a
pi forall a. Num a => a -> a -> a
* Complex Double
g3 forall a. Floating a => a -> a -> a
** (Complex Double
1forall a. Fractional a => a -> a -> a
/Complex Double
6)),
          Double
0.5 forall a. a -> a -> Complex a
:+ (forall a. Floating a => a -> a
sqrt Double
3 forall a. Fractional a => a -> a -> a
/ Double
2)
        )
      | Complex Double
g3 forall a. Eq a => a -> a -> Bool
== Complex Double
0 = 
        (
          Complex Double
i_ forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
sqrt(forall a. Floating a => a -> a
sqrt(Complex Double
3.75 forall a. Num a => a -> a -> a
* Complex Double -> Complex Double
eisensteinG4 Complex Double
tau' forall a. Fractional a => a -> a -> a
/ Complex Double
g2)),
          Complex Double
tau'
        )
      | Bool
otherwise = 
        (
          forall a. Floating a => a -> a
sqrt(Complex Double
7 forall a. Num a => a -> a -> a
* Complex Double
g2 forall a. Num a => a -> a -> a
* Complex Double -> Complex Double
eisensteinG6_over_eisensteinG4 Complex Double
tau' forall a. Fractional a => a -> a -> a
/ (Complex Double
12 forall a. Num a => a -> a -> a
* Complex Double
g3)),
          Complex Double
tau'
        )
      where 
        g2cube :: Complex Double
g2cube = Complex Double
g2 forall a. Num a => a -> a -> a
* Complex Double
g2 forall a. Num a => a -> a -> a
* Complex Double
g2
        j :: Complex Double
j = Complex Double
1728 forall a. Num a => a -> a -> a
* Complex Double
g2cube forall a. Fractional a => a -> a -> a
/ (Complex Double
g2cube forall a. Num a => a -> a -> a
- Complex Double
27 forall a. Num a => a -> a -> a
* Complex Double
g3 forall a. Num a => a -> a -> a
* Complex Double
g3)
        tau' :: Complex Double
tau' = Complex Double -> Complex Double
kleinJinv Complex Double
j

-- | Half-periods from elliptic invariants.
halfPeriods :: 
    Complex Double -- ^ g2
 -> Complex Double -- ^ g3
 -> (Complex Double, Complex Double) -- ^ omega1, omega2
halfPeriods :: Complex Double
-> Complex Double -> (Complex Double, Complex Double)
halfPeriods Complex Double
g2 Complex Double
g3 = (Complex Double
omega1, Complex Double
tau forall a. Num a => a -> a -> a
* Complex Double
omega1)
  where
    (Complex Double
omega1, Complex Double
tau) = Complex Double
-> Complex Double -> (Complex Double, Complex Double)
omega1_and_tau Complex Double
g2 Complex Double
g3

g_from_omega1_and_tau :: 
  Complex Double -> Complex Double -> (Complex Double, Complex Double)
g_from_omega1_and_tau :: Complex Double
-> Complex Double -> (Complex Double, Complex Double)
g_from_omega1_and_tau Complex Double
omega1 Complex Double
tau = (Complex Double
g2, Complex Double
g3)
  where
    q :: Complex Double
q = forall a. Floating a => a -> a
exp (Complex Double
i_ forall a. Num a => a -> a -> a
* forall a. Floating a => a
pi forall a. Num a => a -> a -> a
* Complex Double
tau)
    j2 :: Complex Double
j2 = Complex Double -> Complex Double -> Complex Double
jtheta2 Complex Double
0 Complex Double
q
    j3 :: Complex Double
j3 = Complex Double -> Complex Double -> Complex Double
jtheta3 Complex Double
0 Complex Double
q
    j2pow4 :: Complex Double
j2pow4  = Complex Double
j2 Complex Double -> Int -> Complex Double
%^% Int
4
    j2pow8 :: Complex Double
j2pow8  = Complex Double
j2pow4 forall a. Num a => a -> a -> a
* Complex Double
j2pow4
    j2pow12 :: Complex Double
j2pow12 = Complex Double
j2pow4 forall a. Num a => a -> a -> a
* Complex Double
j2pow8
    j3pow4 :: Complex Double
j3pow4  = Complex Double
j3 Complex Double -> Int -> Complex Double
%^% Int
4
    j3pow8 :: Complex Double
j3pow8  = Complex Double
j3pow4 forall a. Num a => a -> a -> a
* Complex Double
j3pow4
    j3pow12 :: Complex Double
j3pow12 = Complex Double
j3pow4 forall a. Num a => a -> a -> a
* Complex Double
j3pow8
    g2 :: Complex Double
g2 = Complex Double
4forall a. Fractional a => a -> a -> a
/Complex Double
3 forall a. Num a => a -> a -> a
* (forall a. Floating a => a
pi forall a. Fractional a => a -> a -> a
/ Complex Double
2 forall a. Fractional a => a -> a -> a
/ Complex Double
omega1) Complex Double -> Int -> Complex Double
%^% Int
4 forall a. Num a => a -> a -> a
* (Complex Double
j2pow8 forall a. Num a => a -> a -> a
- Complex Double
j2pow4 forall a. Num a => a -> a -> a
* Complex Double
j3pow4 forall a. Num a => a -> a -> a
+ Complex Double
j3pow8)
    g3 :: Complex Double
g3 = Complex Double
8forall a. Fractional a => a -> a -> a
/Complex Double
27 forall a. Num a => a -> a -> a
* (forall a. Floating a => a
pi forall a. Fractional a => a -> a -> a
/ Complex Double
2 forall a. Fractional a => a -> a -> a
/ Complex Double
omega1) Complex Double -> Int -> Complex Double
%^% Int
6 forall a. Num a => a -> a -> a
*
      (Complex Double
j2pow12 forall a. Num a => a -> a -> a
- ((Complex Double
1.5 forall a. Num a => a -> a -> a
* Complex Double
j2pow8 forall a. Num a => a -> a -> a
* Complex Double
j3pow4) forall a. Num a => a -> a -> a
+ (Complex Double
1.5 forall a. Num a => a -> a -> a
* Complex Double
j2pow4 forall a. Num a => a -> a -> a
* Complex Double
j3pow8)) forall a. Num a => a -> a -> a
+ Complex Double
j3pow12)

-- | Elliptic invariants from half-periods.
ellipticInvariants :: 
    Complex Double -- ^ omega1
 -> Complex Double -- ^ omega2
 -> (Complex Double, Complex Double) -- ^ g2, g3
ellipticInvariants :: Complex Double
-> Complex Double -> (Complex Double, Complex Double)
ellipticInvariants Complex Double
omega1 Complex Double
omega2 = 
  Complex Double
-> Complex Double -> (Complex Double, Complex Double)
g_from_omega1_and_tau Complex Double
omega1 (Complex Double
omega2 forall a. Fractional a => a -> a -> a
/ Complex Double
omega1)

weierstrassP_from_tau :: Complex Double -> Complex Double -> Complex Double
weierstrassP_from_tau :: Complex Double -> Complex Double -> Complex Double
weierstrassP_from_tau Complex Double
z Complex Double
tau = 
  (forall a. Floating a => a
pi forall a. Num a => a -> a -> a
* Complex Double
j2 forall a. Num a => a -> a -> a
* Complex Double
j3 forall a. Num a => a -> a -> a
* Complex Double
j4 forall a. Fractional a => a -> a -> a
/ Complex Double
j1) Complex Double -> Int -> Complex Double
%^% Int
2 forall a. Num a => a -> a -> a
- forall a. Floating a => a
pi forall a. Num a => a -> a -> a
* forall a. Floating a => a
pi forall a. Num a => a -> a -> a
* (Complex Double
j2 Complex Double -> Int -> Complex Double
%^% Int
4 forall a. Num a => a -> a -> a
+ Complex Double
j3 Complex Double -> Int -> Complex Double
%^% Int
4) forall a. Fractional a => a -> a -> a
/ Complex Double
3
  where
    q :: Complex Double
q = forall a. Floating a => a -> a
exp (Complex Double
i_ forall a. Num a => a -> a -> a
* forall a. Floating a => a
pi forall a. Num a => a -> a -> a
* Complex Double
tau)
    j2 :: Complex Double
j2 = Complex Double -> Complex Double -> Complex Double
jtheta2 Complex Double
0 Complex Double
q
    j3 :: Complex Double
j3 = Complex Double -> Complex Double -> Complex Double
jtheta3 Complex Double
0 Complex Double
q
    z' :: Complex Double
z' = forall a. Floating a => a
pi forall a. Num a => a -> a -> a
* Complex Double
z
    j1 :: Complex Double
j1 = Complex Double -> Complex Double -> Complex Double
jtheta1 Complex Double
z' Complex Double
q
    j4 :: Complex Double
j4 = Complex Double -> Complex Double -> Complex Double
jtheta4 Complex Double
z' Complex Double
q

weierstrassP_from_omega :: 
  Complex Double -> Complex Double -> Complex Double -> Complex Double
weierstrassP_from_omega :: Complex Double
-> Complex Double -> Complex Double -> Complex Double
weierstrassP_from_omega Complex Double
z Complex Double
omega1 Complex Double
omega2 = 
  Complex Double -> Complex Double -> Complex Double
weierstrassP_from_tau 
    (Complex Double
z forall a. Fractional a => a -> a -> a
/ Complex Double
omega1 forall a. Fractional a => a -> a -> a
/ Complex Double
2) (Complex Double
omega2 forall a. Fractional a => a -> a -> a
/ Complex Double
omega1) forall a. Fractional a => a -> a -> a
/ (Complex Double
4 forall a. Num a => a -> a -> a
* Complex Double
omega1 forall a. Num a => a -> a -> a
* Complex Double
omega1)

-- | Weierstrass p-function
weierstrassP ::
    Complex Double -- ^ z
 -> Complex Double -- ^ elliptic invariant g2
 -> Complex Double -- ^ elliptic invariant g3
 -> Complex Double
weierstrassP :: Complex Double
-> Complex Double -> Complex Double -> Complex Double
weierstrassP Complex Double
z Complex Double
g2 Complex Double
g3 = Complex Double
-> Complex Double -> Complex Double -> Complex Double
weierstrassP_from_omega Complex Double
z Complex Double
omega1 Complex Double
omega2
  where
    (Complex Double
omega1, Complex Double
omega2) = Complex Double
-> Complex Double -> (Complex Double, Complex Double)
halfPeriods Complex Double
g2 Complex Double
g3

-- | Derivative of Weierstrass p-function
weierstrassPdash ::
    Complex Double -- ^ z
 -> Complex Double -- ^ elliptic invariant g2
 -> Complex Double -- ^ elliptic invariant g3
 -> Complex Double
weierstrassPdash :: Complex Double
-> Complex Double -> Complex Double -> Complex Double
weierstrassPdash Complex Double
z Complex Double
g2 Complex Double
g3 = Complex Double
2 forall a. Fractional a => a -> a -> a
/ (Complex Double
w1 Complex Double -> Int -> Complex Double
%^% Int
3) forall a. Num a => a -> a -> a
* Complex Double
j2 forall a. Num a => a -> a -> a
* Complex Double
j3 forall a. Num a => a -> a -> a
* Complex Double
j4 forall a. Num a => a -> a -> a
* Complex Double
f
  where
    (Complex Double
omega1, Complex Double
omega2) = Complex Double
-> Complex Double -> (Complex Double, Complex Double)
halfPeriods Complex Double
g2 Complex Double
g3
    w1 :: Complex Double
w1 = Complex Double
2 forall a. Num a => a -> a -> a
* Complex Double
omega1 forall a. Fractional a => a -> a -> a
/ forall a. Floating a => a
pi
    tau :: Complex Double
tau = Complex Double
omega2 forall a. Fractional a => a -> a -> a
/ Complex Double
omega1
    q :: Complex Double
q = forall a. Floating a => a -> a
exp (Complex Double
i_ forall a. Num a => a -> a -> a
* forall a. Floating a => a
pi forall a. Num a => a -> a -> a
* Complex Double
tau)
    z' :: Complex Double
z' = Complex Double
z forall a. Fractional a => a -> a -> a
/ Complex Double
w1 
    j1 :: Complex Double
j1 = Complex Double -> Complex Double -> Complex Double
jtheta1 Complex Double
z' Complex Double
q
    j2 :: Complex Double
j2 = Complex Double -> Complex Double -> Complex Double
jtheta2 Complex Double
z' Complex Double
q
    j3 :: Complex Double
j3 = Complex Double -> Complex Double -> Complex Double
jtheta3 Complex Double
z' Complex Double
q
    j4 :: Complex Double
j4 = Complex Double -> Complex Double -> Complex Double
jtheta4 Complex Double
z' Complex Double
q
    j1dash :: Complex Double
j1dash = Complex Double -> Complex Double -> Complex Double
jtheta1Dash Complex Double
0 Complex Double
q
    j2zero :: Complex Double
j2zero = Complex Double -> Complex Double -> Complex Double
jtheta2 Complex Double
0 Complex Double
q
    j3zero :: Complex Double
j3zero = Complex Double -> Complex Double -> Complex Double
jtheta3 Complex Double
0 Complex Double
q
    j4zero :: Complex Double
j4zero = Complex Double -> Complex Double -> Complex Double
jtheta4 Complex Double
0 Complex Double
q
    f :: Complex Double
f = Complex Double
j1dash Complex Double -> Int -> Complex Double
%^% Int
3 forall a. Fractional a => a -> a -> a
/ (Complex Double
j1 Complex Double -> Int -> Complex Double
%^% Int
3 forall a. Num a => a -> a -> a
* Complex Double
j2zero forall a. Num a => a -> a -> a
* Complex Double
j3zero forall a. Num a => a -> a -> a
* Complex Double
j4zero)

-- | Inverse of Weierstrass p-function
weierstrassPinv ::
    Complex Double -- ^ w
 -> Complex Double -- ^ elliptic invariant g2
 -> Complex Double -- ^ elliptic invariant g3
 -> Complex Double
weierstrassPinv :: Complex Double
-> Complex Double -> Complex Double -> Complex Double
weierstrassPinv Complex Double
w Complex Double
g2 Complex Double
g3 = Double
-> Complex Double
-> Complex Double
-> Complex Double
-> Complex Double
carlsonRF' Double
1e-14 (Complex Double
w forall a. Num a => a -> a -> a
- Complex Double
e1) (Complex Double
w forall a. Num a => a -> a -> a
- Complex Double
e2) (Complex Double
w forall a. Num a => a -> a -> a
- Complex Double
e3)
  where
    (Complex Double
omega1, Complex Double
omega2) = Complex Double
-> Complex Double -> (Complex Double, Complex Double)
halfPeriods Complex Double
g2 Complex Double
g3
    e1 :: Complex Double
e1 = Complex Double
-> Complex Double -> Complex Double -> Complex Double
weierstrassP Complex Double
omega1 Complex Double
g2 Complex Double
g3
    e2 :: Complex Double
e2 = Complex Double
-> Complex Double -> Complex Double -> Complex Double
weierstrassP Complex Double
omega2 Complex Double
g2 Complex Double
g3
    e3 :: Complex Double
e3 = Complex Double
-> Complex Double -> Complex Double -> Complex Double
weierstrassP (-Complex Double
omega1 forall a. Num a => a -> a -> a
- Complex Double
omega2) Complex Double
g2 Complex Double
g3

-- | Weierstrass sigma function
weierstrassSigma ::
    Complex Double -- ^ z
 -> Complex Double -- ^ elliptic invariant g2
 -> Complex Double -- ^ elliptic invariant g3
 -> Complex Double
weierstrassSigma :: Complex Double
-> Complex Double -> Complex Double -> Complex Double
weierstrassSigma Complex Double
z Complex Double
g2 Complex Double
g3 = Complex Double
w1 forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
exp (Complex Double
h forall a. Num a => a -> a -> a
* Complex Double
z forall a. Num a => a -> a -> a
* Complex Double
z1 forall a. Fractional a => a -> a -> a
/ forall a. Floating a => a
pi) forall a. Num a => a -> a -> a
* Complex Double
j1 forall a. Fractional a => a -> a -> a
/ Complex Double
j1dash
  where
    (Complex Double
omega1, Complex Double
omega2) = Complex Double
-> Complex Double -> (Complex Double, Complex Double)
halfPeriods Complex Double
g2 Complex Double
g3
    tau :: Complex Double
tau = Complex Double
omega2 forall a. Fractional a => a -> a -> a
/ Complex Double
omega1
    q :: Complex Double
q = forall a. Floating a => a -> a
exp (Complex Double
i_ forall a. Num a => a -> a -> a
* forall a. Floating a => a
pi forall a. Num a => a -> a -> a
* Complex Double
tau)
    w1 :: Complex Double
w1 = -Complex Double
2 forall a. Num a => a -> a -> a
* Complex Double
omega1 forall a. Fractional a => a -> a -> a
/ forall a. Floating a => a
pi
    z1 :: Complex Double
z1 = Complex Double
z forall a. Fractional a => a -> a -> a
/ Complex Double
w1
    j1 :: Complex Double
j1 = Complex Double -> Complex Double -> Complex Double
jtheta1 Complex Double
z1 Complex Double
q
    j1dash :: Complex Double
j1dash = Complex Double -> Complex Double -> Complex Double
jtheta1Dash Complex Double
0 Complex Double
q
    h :: Complex Double
h = - forall a. Floating a => a
pi forall a. Fractional a => a -> a -> a
/ (Complex Double
6 forall a. Num a => a -> a -> a
* Complex Double
w1) forall a. Num a => a -> a -> a
* Complex Double -> Complex Double
jtheta1DashDashDash0 Complex Double
tau forall a. Fractional a => a -> a -> a
/ Complex Double
j1dash

-- | Weierstrass zeta function
weierstrassZeta ::
    Complex Double -- ^ z
 -> Complex Double -- ^ elliptic invariant g2
 -> Complex Double -- ^ elliptic invariant g3
 -> Complex Double
weierstrassZeta :: Complex Double
-> Complex Double -> Complex Double -> Complex Double
weierstrassZeta Complex Double
z Complex Double
g2 Complex Double
g3 = - Complex Double
eta1 forall a. Num a => a -> a -> a
* Complex Double
z forall a. Num a => a -> a -> a
+ Complex Double
p forall a. Num a => a -> a -> a
* Complex Double
lj1dash
  where
    (Complex Double
omega1, Complex Double
omega2) = Complex Double
-> Complex Double -> (Complex Double, Complex Double)
halfPeriods Complex Double
g2 Complex Double
g3
    tau :: Complex Double
tau = Complex Double
omega2 forall a. Fractional a => a -> a -> a
/ Complex Double
omega1
    q :: Complex Double
q = forall a. Floating a => a -> a
exp (Complex Double
i_ forall a. Num a => a -> a -> a
* forall a. Floating a => a
pi forall a. Num a => a -> a -> a
* Complex Double
tau)
    w1 :: Complex Double
w1 = - Complex Double
omega1 forall a. Fractional a => a -> a -> a
/ forall a. Floating a => a
pi
    p :: Complex Double
p = Complex Double
0.5 forall a. Fractional a => a -> a -> a
/ Complex Double
w1
    j1dash :: Complex Double
j1dash = Complex Double -> Complex Double -> Complex Double
jtheta1Dash Complex Double
0 Complex Double
q
    eta1 :: Complex Double
eta1 = Complex Double
p forall a. Num a => a -> a -> a
* Complex Double -> Complex Double
jtheta1DashDashDash0 Complex Double
tau forall a. Fractional a => a -> a -> a
/ (Complex Double
6 forall a. Num a => a -> a -> a
* Complex Double
w1 forall a. Num a => a -> a -> a
* Complex Double
j1dash)
    pz :: Complex Double
pz = Complex Double
p forall a. Num a => a -> a -> a
* Complex Double
z
    lj1dash :: Complex Double
lj1dash = Complex Double -> Complex Double -> Complex Double
jtheta1Dash Complex Double
pz Complex Double
q forall a. Fractional a => a -> a -> a
/ Complex Double -> Complex Double -> Complex Double
jtheta1 Complex Double
pz Complex Double
q