{-|
Module      : Math.JacobiTheta
Description : Jacobi theta functions.
Copyright   : (c) Stéphane Laurent, 2023
License     : BSD3
Maintainer  : laurent_step@outlook.fr

Provides the four usual Jacobi theta functions, the Jacobi theta function 
with characteristics, the derivative of the first Jacobi theta function, 
as well as a function for the derivative at @0@ only of the first Jacobi 
theta function.
-}
module Math.JacobiTheta
  (
    jtheta1,
    jtheta1',
    jtheta2,
    jtheta2',
    jtheta3,
    jtheta3',
    jtheta4,
    jtheta4',
    jthetaAB,
    jthetaAB',
    jtheta1Dash0,
    jtheta1Dash 
  )
  where
import Data.Complex ( imagPart, magnitude, realPart, Complex(..) )

type Cplx = Complex Double

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

machinePrecision :: Double
machinePrecision :: Double
machinePrecision = Double
2forall a. Floating a => a -> a -> a
**(-Double
52)

areClose :: Cplx -> Cplx -> Bool
areClose :: Cplx -> Cplx -> Bool
areClose Cplx
z1 Cplx
z2 = forall a. RealFloat a => Complex a -> a
magnitude (Cplx
z1 forall a. Num a => a -> a -> a
- Cplx
z2) forall a. Ord a => a -> a -> Bool
< Double
epsilon forall a. Num a => a -> a -> a
* Double
h
  where
    epsilon :: Double
epsilon = Double
2.0 forall a. Num a => a -> a -> a
* Double
machinePrecision
    magn2 :: Double
magn2 = forall a. RealFloat a => Complex a -> a
magnitude Cplx
z2
    h :: Double
h = if Double
magn2 forall a. Ord a => a -> a -> Bool
< Double
epsilon then Double
1.0 else forall a. Ord a => a -> a -> a
max (forall a. RealFloat a => Complex a -> a
magnitude Cplx
z1) Double
magn2

modulo :: Double -> Int -> Double
modulo :: Double -> Int -> Double
modulo Double
a Int
p = 
  let p' :: Double
p' = forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
p
  in
  if Double
a forall a. Ord a => a -> a -> Bool
> Double
0 
    then Double
a forall a. Num a => a -> a -> a
- forall a b. (Integral a, Num b) => a -> b
fromIntegral(Int
p forall a. Num a => a -> a -> a
* forall a b. (RealFrac a, Integral b) => a -> b
floor(Double
aforall a. Fractional a => a -> a -> a
/Double
p'))  
    else Double
a forall a. Num a => a -> a -> a
- forall a b. (Integral a, Num b) => a -> b
fromIntegral(Int
p forall a. Num a => a -> a -> a
* forall a b. (RealFrac a, Integral b) => a -> b
ceiling(Double
aforall a. Fractional a => a -> a -> a
/Double
p'))

dologtheta4 :: Cplx -> Cplx -> Int -> Int -> Cplx
dologtheta4 :: Cplx -> Cplx -> Int -> Int -> Cplx
dologtheta4 Cplx
z Cplx
tau Int
passes Int
maxiter = 
  Cplx -> Cplx -> Int -> Int -> Cplx
dologtheta3 (Cplx
z forall a. Num a => a -> a -> a
+ Cplx
0.5) Cplx
tau (Int
passesforall a. Num a => a -> a -> a
+Int
1) Int
maxiter

dologtheta3 :: Cplx -> Cplx -> Int -> Int -> Cplx
dologtheta3 :: Cplx -> Cplx -> Int -> Int -> Cplx
dologtheta3 Cplx
z Cplx
tau Int
passes Int
maxiterloc
  | forall a. Complex a -> a
realPart Cplx
tau2 forall a. Ord a => a -> a -> Bool
> Double
0.6  = Cplx -> Cplx -> Int -> Int -> Cplx
dologtheta4 Cplx
z (Cplx
tau2 forall a. Num a => a -> a -> a
- Cplx
1) (Int
passes forall a. Num a => a -> a -> a
+ Int
1) Int
maxiterloc
  | forall a. Complex a -> a
realPart Cplx
tau2 forall a. Ord a => a -> a -> Bool
< -Double
0.6 = Cplx -> Cplx -> Int -> Int -> Cplx
dologtheta4 Cplx
z (Cplx
tau2 forall a. Num a => a -> a -> a
+ Cplx
1) (Int
passes forall a. Num a => a -> a -> a
+ Int
1) Int
maxiterloc
  | forall a. RealFloat a => Complex a -> a
magnitude Cplx
tau2 forall a. Ord a => a -> a -> Bool
< Double
0.98 Bool -> Bool -> Bool
&& forall a. Complex a -> a
imagPart Cplx
tau2 forall a. Ord a => a -> a -> Bool
< Double
0.98 = 
      Cplx
i_ forall a. Num a => a -> a -> a
* forall a. Floating a => a
pi forall a. Num a => a -> a -> a
* Cplx
tauprime forall a. Num a => a -> a -> a
* Cplx
z forall a. Num a => a -> a -> a
* Cplx
z 
      forall a. Num a => a -> a -> a
+ Cplx -> Cplx -> Int -> Int -> Cplx
dologtheta3 (Cplx
z forall a. Num a => a -> a -> a
* Cplx
tauprime) Cplx
tauprime (Int
passes forall a. Num a => a -> a -> a
+ Int
1) Int
maxiterloc
      forall a. Num a => a -> a -> a
- forall a. Floating a => a -> a
log(forall a. Floating a => a -> a
sqrt Cplx
tau2 forall a. Fractional a => a -> a -> a
/ forall a. Floating a => a -> a
sqrt Cplx
i_) 
  | Bool
otherwise = Cplx -> Cplx -> Int -> Int -> Cplx
argtheta3 Cplx
z Cplx
tau2 Int
0 Int
maxiterloc
    where
      rPtau :: Double
rPtau = forall a. Complex a -> a
realPart Cplx
tau
      rPtau2 :: Double
rPtau2 = if Double
rPtau forall a. Ord a => a -> a -> Bool
> Double
0
        then Double -> Int -> Double
modulo (Double
rPtau forall a. Num a => a -> a -> a
+ Double
1) Int
2 forall a. Num a => a -> a -> a
- Double
1
        else Double -> Int -> Double
modulo (Double
rPtau forall a. Num a => a -> a -> a
- Double
1) Int
2 forall a. Num a => a -> a -> a
+ Double
1
      tau2 :: Cplx
tau2 = Double
rPtau2 forall a. a -> a -> Complex a
:+ forall a. Complex a -> a
imagPart Cplx
tau
      tauprime :: Cplx
tauprime = -Cplx
1 forall a. Fractional a => a -> a -> a
/ Cplx
tau2

argtheta3 :: Cplx -> Cplx -> Int -> Int -> Cplx
argtheta3 :: Cplx -> Cplx -> Int -> Int -> Cplx
argtheta3 Cplx
z Cplx
tau Int
passes Int
maxiterloc
  | Int
passes forall a. Ord a => a -> a -> Bool
> Int
maxiterloc = forall a. HasCallStack => [Char] -> a
error [Char]
"Reached maximal iteration."
  | Double
iPz forall a. Ord a => a -> a -> Bool
< -Double
iPtau forall a. Fractional a => a -> a -> a
/ Double
2 = Cplx -> Cplx -> Int -> Int -> Cplx
argtheta3 (-Cplx
zuse) Cplx
tau (Int
passes forall a. Num a => a -> a -> a
+ Int
1) Int
maxiterloc
  | Double
iPz forall a. Ord a => a -> a -> Bool
>= Double
iPtau forall a. Fractional a => a -> a -> a
/ Double
2 = 
      -Cplx
2 forall a. Num a => a -> a -> a
* forall a. Floating a => a
pi forall a. Num a => a -> a -> a
* Cplx
quotient forall a. Num a => a -> a -> a
* Cplx
i_ forall a. Num a => a -> a -> a
* Cplx
zmin 
      forall a. Num a => a -> a -> a
+ Cplx -> Cplx -> Int -> Int -> Cplx
argtheta3 Cplx
zmin Cplx
tau (Int
passes forall a. Num a => a -> a -> a
+ Int
1) Int
maxiterloc
      forall a. Num a => a -> a -> a
- Cplx
i_ forall a. Num a => a -> a -> a
* forall a. Floating a => a
pi forall a. Num a => a -> a -> a
* Cplx
tau forall a. Num a => a -> a -> a
* Cplx
quotient forall a. Num a => a -> a -> a
* Cplx
quotient
  | Bool
otherwise = Cplx -> Cplx -> Cplx
calctheta3 Cplx
zuse Cplx
tau
    where
      iPz :: Double
iPz = forall a. Complex a -> a
imagPart Cplx
z
      iPtau :: Double
iPtau = forall a. Complex a -> a
imagPart Cplx
tau
      zuse :: Cplx
zuse = Double -> Int -> Double
modulo (forall a. Complex a -> a
realPart Cplx
z) Int
1 forall a. a -> a -> Complex a
:+ Double
iPz
      quotient :: Cplx
quotient = Int -> Cplx
fromInt forall a b. (a -> b) -> a -> b
$ forall a b. (RealFrac a, Integral b) => a -> b
floor(Double
iPz forall a. Fractional a => a -> a -> a
/ Double
iPtau forall a. Num a => a -> a -> a
+ Double
0.5)
      zmin :: Cplx
zmin = Cplx
zuse forall a. Num a => a -> a -> a
- Cplx
tau forall a. Num a => a -> a -> a
* Cplx
quotient
      fromInt :: Int -> Cplx
      fromInt :: Int -> Cplx
fromInt = forall a b. (Integral a, Num b) => a -> b
fromIntegral

calctheta3 :: Cplx -> Cplx -> Cplx
calctheta3 :: Cplx -> Cplx -> Cplx
calctheta3 Cplx
z Cplx
tau = 
    Int -> Cplx -> Cplx
go Int
1 Cplx
1
    where
      qw :: Int -> Cplx
      qw :: Int -> Cplx
qw Int
n = forall a. Floating a => a -> a
exp(Cplx
inpi forall a. Num a => a -> a -> a
* (Cplx
taun forall a. Num a => a -> a -> a
+ Cplx
2 forall a. Num a => a -> a -> a
* Cplx
z)) forall a. Num a => a -> a -> a
+ forall a. Floating a => a -> a
exp(Cplx
inpi forall a. Num a => a -> a -> a
* (Cplx
taun forall a. Num a => a -> a -> a
- Cplx
2 forall a. Num a => a -> a -> a
* Cplx
z))
        where
          n' :: Cplx
n' = forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n 
          inpi :: Cplx
inpi = Cplx
i_ forall a. Num a => a -> a -> a
* Cplx
n' forall a. Num a => a -> a -> a
* forall a. Floating a => a
pi 
          taun :: Cplx
taun = Cplx
n' forall a. Num a => a -> a -> a
* Cplx
tau     
      go :: Int -> Cplx -> Cplx
go Int
n Cplx
res
        | forall a. RealFloat a => a -> Bool
isNaN Double
modulus = forall a. HasCallStack => [Char] -> a
error [Char]
"NaN has occured in the summation."
        | forall a. RealFloat a => a -> Bool
isInfinite Double
modulus = forall a. HasCallStack => [Char] -> a
error [Char]
"Infinity reached in the summation."
--        | modulus == 0 = error "Zero has occured in the summation."

        | Int
n forall a. Ord a => a -> a -> Bool
>= Int
3 Bool -> Bool -> Bool
&& Cplx -> Cplx -> Bool
areClose Cplx
res Cplx
resnew = forall a. Floating a => a -> a
log Cplx
res
        | Bool
otherwise = Int -> Cplx -> Cplx
go (Int
n forall a. Num a => a -> a -> a
+ Int
1) Cplx
resnew
          where
            modulus :: Double
modulus = forall a. RealFloat a => Complex a -> a
magnitude Cplx
res
            resnew :: Cplx
resnew = Cplx
res forall a. Num a => a -> a -> a
+ Int -> Cplx
qw Int
n

-------------------------------------------------------------------------------

tauFromQ :: Cplx -> Cplx
tauFromQ :: Cplx -> Cplx
tauFromQ Cplx
q = if forall a. Complex a -> a
imagPart Cplx
q forall a. Eq a => a -> a -> Bool
== Double
0 Bool -> Bool -> Bool
&& forall a. Complex a -> a
realPart Cplx
q forall a. Ord a => a -> a -> Bool
< Double
0 
  then Double
1 forall a. a -> a -> Complex a
:+ (-forall a. Floating a => a -> a
log(-(forall a. Complex a -> a
realPart Cplx
q)) forall a. Fractional a => a -> a -> a
/ forall a. Floating a => a
pi)
  else -Cplx
i_ forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
log Cplx
q forall a. Fractional a => a -> a -> a
/ forall a. Floating a => a
pi

checkQ :: Cplx -> Cplx
checkQ :: Cplx -> Cplx
checkQ Cplx
q
  | forall a. RealFloat a => Complex a -> a
magnitude Cplx
q forall a. Ord a => a -> a -> Bool
>= Double
1 = 
    forall a. HasCallStack => [Char] -> a
error [Char]
"The modulus of the nome must be smaller than one."
  | Cplx
q forall a. Eq a => a -> a -> Bool
== Cplx
0 = 
    forall a. HasCallStack => [Char] -> a
error [Char]
"The nome cannot be zero."
  | Bool
otherwise = Cplx
q

getTauFromQ :: Cplx -> Cplx
getTauFromQ :: Cplx -> Cplx
getTauFromQ = Cplx -> Cplx
tauFromQ forall b c a. (b -> c) -> (a -> b) -> a -> c
. Cplx -> Cplx
checkQ

funM :: Cplx -> Cplx -> Cplx
funM :: Cplx -> Cplx -> Cplx
funM Cplx
z Cplx
tau = Cplx
i_ forall a. Num a => a -> a -> a
* forall a. Floating a => a
pi forall a. Num a => a -> a -> a
* (Cplx
z forall a. Num a => a -> a -> a
+ Cplx
tauforall a. Fractional a => a -> a -> a
/Cplx
4)

ljtheta1 :: Cplx -> Cplx -> Cplx
ljtheta1 :: Cplx -> Cplx -> Cplx
ljtheta1 Cplx
z Cplx
tau = Cplx -> Cplx -> Cplx
ljtheta2 (Cplx
z forall a. Num a => a -> a -> a
- Cplx
0.5) Cplx
tau

-- | First Jacobi theta function in function of the nome.

jtheta1 ::
     Complex Double -- ^ z

  -> Complex Double -- ^ q, the nome

  -> Complex Double
jtheta1 :: Cplx -> Cplx -> Cplx
jtheta1 Cplx
z Cplx
q = forall a. Floating a => a -> a
exp(Cplx -> Cplx -> Cplx
ljtheta1 (Cplx
zforall a. Fractional a => a -> a -> a
/forall a. Floating a => a
pi) Cplx
tau)
  where
    tau :: Cplx
tau = Cplx -> Cplx
getTauFromQ Cplx
q

-- | First Jacobi theta function in function of @tau@.

jtheta1' ::
     Complex Double -- ^ z

  -> Complex Double -- ^ tau

  -> Complex Double
jtheta1' :: Cplx -> Cplx -> Cplx
jtheta1' Cplx
z Cplx
tau
  | forall a. Complex a -> a
imagPart Cplx
tau forall a. Ord a => a -> a -> Bool
<= Double
0 = forall a. HasCallStack => [Char] -> a
error [Char]
"`tau` must have a nonnegative imaginary part."
  | Bool
otherwise = forall a. Floating a => a -> a
exp(Cplx -> Cplx -> Cplx
ljtheta1 (Cplx
zforall a. Fractional a => a -> a -> a
/forall a. Floating a => a
pi) Cplx
tau)

ljtheta2 :: Cplx -> Cplx -> Cplx
ljtheta2 :: Cplx -> Cplx -> Cplx
ljtheta2 Cplx
z Cplx
tau = 
  Cplx -> Cplx -> Cplx
funM Cplx
z Cplx
tau forall a. Num a => a -> a -> a
+ Cplx -> Cplx -> Int -> Int -> Cplx
dologtheta3 (Cplx
z forall a. Num a => a -> a -> a
+ Cplx
0.5 forall a. Num a => a -> a -> a
* Cplx
tau) Cplx
tau Int
0 Int
1000

-- | Second Jacobi theta function in function of the nome.

jtheta2 ::
     Complex Double -- ^ z

  -> Complex Double -- ^ q, the nome

  -> Complex Double
jtheta2 :: Cplx -> Cplx -> Cplx
jtheta2 Cplx
z Cplx
q = forall a. Floating a => a -> a
exp(Cplx -> Cplx -> Cplx
ljtheta2 (Cplx
zforall a. Fractional a => a -> a -> a
/forall a. Floating a => a
pi) Cplx
tau)
  where
    tau :: Cplx
tau = Cplx -> Cplx
getTauFromQ Cplx
q

-- | Second Jacobi theta function in function of @tau@.

jtheta2' ::
     Complex Double -- ^ z

  -> Complex Double -- ^ tau

  -> Complex Double
jtheta2' :: Cplx -> Cplx -> Cplx
jtheta2' Cplx
z Cplx
tau
  | forall a. Complex a -> a
imagPart Cplx
tau forall a. Ord a => a -> a -> Bool
<= Double
0 = forall a. HasCallStack => [Char] -> a
error [Char]
"`tau` must have a nonnegative imaginary part."
  | Bool
otherwise = forall a. Floating a => a -> a
exp(Cplx -> Cplx -> Cplx
ljtheta2 (Cplx
zforall a. Fractional a => a -> a -> a
/forall a. Floating a => a
pi) Cplx
tau)

-- | Third Jacobi theta function in function of the nome.

jtheta3 ::
     Complex Double -- ^ z

  -> Complex Double -- ^ q, the nome

  -> Complex Double
jtheta3 :: Cplx -> Cplx -> Cplx
jtheta3 Cplx
z Cplx
q = forall a. Floating a => a -> a
exp(Cplx -> Cplx -> Int -> Int -> Cplx
dologtheta3 (Cplx
zforall a. Fractional a => a -> a -> a
/forall a. Floating a => a
pi) Cplx
tau Int
0 Int
1000)
  where
    tau :: Cplx
tau = Cplx -> Cplx
getTauFromQ Cplx
q

-- | Third Jacobi theta function in function of @tau@.

jtheta3' ::
     Complex Double -- ^ z

  -> Complex Double -- ^ tau

  -> Complex Double
jtheta3' :: Cplx -> Cplx -> Cplx
jtheta3' Cplx
z Cplx
tau
  | forall a. Complex a -> a
imagPart Cplx
tau forall a. Ord a => a -> a -> Bool
<= Double
0 = forall a. HasCallStack => [Char] -> a
error [Char]
"`tau` must have a nonnegative imaginary part."
  | Bool
otherwise = forall a. Floating a => a -> a
exp(Cplx -> Cplx -> Int -> Int -> Cplx
dologtheta3 (Cplx
zforall a. Fractional a => a -> a -> a
/forall a. Floating a => a
pi) Cplx
tau Int
0 Int
1000)

-- | Fourth Jacobi theta function in function of the nome.

jtheta4 ::
     Complex Double -- ^ z

  -> Complex Double -- ^ q, the nome

  -> Complex Double
jtheta4 :: Cplx -> Cplx -> Cplx
jtheta4 Cplx
z Cplx
q = forall a. Floating a => a -> a
exp(Cplx -> Cplx -> Int -> Int -> Cplx
dologtheta4 (Cplx
zforall a. Fractional a => a -> a -> a
/forall a. Floating a => a
pi) Cplx
tau Int
0 Int
1000)
  where
    tau :: Cplx
tau = Cplx -> Cplx
getTauFromQ Cplx
q

-- | Fourth Jacobi theta function in function of @tau@.

jtheta4' ::
     Complex Double -- ^ z

  -> Complex Double -- ^ tau

  -> Complex Double
jtheta4' :: Cplx -> Cplx -> Cplx
jtheta4' Cplx
z Cplx
tau
  | forall a. Complex a -> a
imagPart Cplx
tau forall a. Ord a => a -> a -> Bool
<= Double
0 = forall a. HasCallStack => [Char] -> a
error [Char]
"`tau` must have a nonnegative imaginary part."
  | Bool
otherwise = forall a. Floating a => a -> a
exp(Cplx -> Cplx -> Int -> Int -> Cplx
dologtheta4 (Cplx
zforall a. Fractional a => a -> a -> a
/forall a. Floating a => a
pi) Cplx
tau Int
0 Int
1000)

-- | Jacobi theta function with characteristics. This is a family of functions, 

--  containing the opposite of the first Jacobi theta function (@a=b=0.5@), 

--  the second Jacobi theta function (@a=0.5, b=0@), the third Jacobi theta 

--  function (@a=b=0@) and the fourth Jacobi theta function (@a=0, b=0.5@). 

--  The examples below show the periodicity-like properties of these functions:

--  

-- >>> import Data.Complex

-- >>> a = 2 :+ 0.3

-- >>> b = 1 :+ (-0.6)

-- >>> z = 0.1 :+ 0.4

-- >>> tau = 0.2 :+ 0.3

-- >>> im = 0 :+ 1 

-- >>> q = exp(im * pi * tau)

-- >>> jab = jthetaAB a b z q

-- >>> jthetaAB a b (z + pi) q

-- (-5.285746223832433e-3) :+ 0.1674462628348814

-- 

-- >>> jab * exp(2 * im * pi * a)

-- (-5.285746223831987e-3) :+ 0.16744626283488154

-- 

-- >>> jtheta_ab a b (z + pi*tau) q

-- 0.10389127606987271 :+ 0.10155646232306936

-- 

-- >>> jab * exp(-im * (pi*tau + 2*z + 2*pi*b))

-- 0.10389127606987278 :+ 0.10155646232306961

jthetaAB ::
     Complex Double -- ^ characteristic a

  -> Complex Double -- ^ characteristic b

  -> Complex Double -- ^ z

  -> Complex Double -- ^ q, the nome

  -> Complex Double
jthetaAB :: Cplx -> Cplx -> Cplx -> Cplx -> Cplx
jthetaAB Cplx
a Cplx
b Cplx
z Cplx
q = Cplx
c forall a. Num a => a -> a -> a
* Cplx -> Cplx -> Cplx
jtheta3 (Cplx
alpha forall a. Num a => a -> a -> a
+ Cplx
beta) Cplx
q
  where
    tau :: Cplx
tau = Cplx -> Cplx
getTauFromQ Cplx
q
    alpha :: Cplx
alpha = forall a. Floating a => a
pi forall a. Num a => a -> a -> a
* Cplx
a forall a. Num a => a -> a -> a
* Cplx
tau 
    beta :: Cplx
beta  = Cplx
z forall a. Num a => a -> a -> a
+ forall a. Floating a => a
pi forall a. Num a => a -> a -> a
* Cplx
b
    c :: Cplx
c     = forall a. Floating a => a -> a
exp(Cplx
i_ forall a. Num a => a -> a -> a
* Cplx
a forall a. Num a => a -> a -> a
* (Cplx
alpha forall a. Num a => a -> a -> a
+ Cplx
2forall a. Num a => a -> a -> a
*Cplx
beta)) 
    -- c     = q**(a*a) * exp(2 * i_ * a * beta)


-- | Jacobi theta function with characteristics in function of @tau@.

jthetaAB' ::
     Complex Double -- ^ characteristic a

  -> Complex Double -- ^ characteristic b

  -> Complex Double -- ^ z

  -> Complex Double -- ^ tau

  -> Complex Double
jthetaAB' :: Cplx -> Cplx -> Cplx -> Cplx -> Cplx
jthetaAB' Cplx
a Cplx
b Cplx
z Cplx
tau = if forall a. Complex a -> a
imagPart Cplx
tau forall a. Ord a => a -> a -> Bool
<= Double
0
  then forall a. HasCallStack => [Char] -> a
error [Char]
"`tau` must have a nonnegative imaginary part."
  else Cplx
c forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
exp(Cplx -> Cplx -> Int -> Int -> Cplx
dologtheta3 (Cplx
alphaforall a. Num a => a -> a -> a
+Cplx
beta) Cplx
tau Int
0 Int
1000)
  where
    alpha :: Cplx
alpha = Cplx
a forall a. Num a => a -> a -> a
* Cplx
tau 
    beta :: Cplx
beta  = Cplx
zforall a. Fractional a => a -> a -> a
/forall a. Floating a => a
pi forall a. Num a => a -> a -> a
+ Cplx
b
    c :: Cplx
c     = forall a. Floating a => a -> a
exp(Cplx
i_ forall a. Num a => a -> a -> a
* forall a. Floating a => a
pi forall a. Num a => a -> a -> a
* Cplx
a forall a. Num a => a -> a -> a
* (Cplx
alpha forall a. Num a => a -> a -> a
+ Cplx
2forall a. Num a => a -> a -> a
*Cplx
beta)) 

-- | Derivative at 0 of the first Jacobi theta function. This is much more 

--  efficient than evaluating @jtheta1Dash@ at @0@.

jtheta1Dash0 :: 
     Complex Double -- ^ q, the nome

  -> Complex Double
jtheta1Dash0 :: Cplx -> Cplx
jtheta1Dash0 Cplx
q = 
  -Cplx
2 forall a. Num a => a -> a -> a
* Cplx
i_ forall a. Num a => a -> a -> a
* Cplx
jab forall a. Num a => a -> a -> a
* Cplx
jab forall a. Num a => a -> a -> a
* Cplx
jab
  where
    tau :: Cplx
tau = Cplx -> Cplx
getTauFromQ Cplx
q
    jab :: Cplx
jab = Cplx -> Cplx -> Cplx -> Cplx -> Cplx
jthetaAB' (Cplx
1forall a. Fractional a => a -> a -> a
/Cplx
6) Cplx
0.5 Cplx
0 (Cplx
3forall a. Num a => a -> a -> a
*Cplx
tau)

-- | Derivative of the first Jacobi theta function.

jtheta1Dash :: 
     Complex Double -- ^ z

  -> Complex Double -- ^ q, the nome

  -> Complex Double
jtheta1Dash :: Cplx -> Cplx -> Cplx
jtheta1Dash Cplx
z Cplx
q = 
  Int -> Cplx -> Cplx -> Cplx -> Cplx -> Cplx
go Int
0 (Double
0.0 forall a. a -> a -> Complex a
:+ Double
0.0) Cplx
1.0 (Cplx
1.0 forall a. Fractional a => a -> a -> a
/ Cplx
qsq) Cplx
1.0
  where 
    q' :: Cplx
q' = Cplx -> Cplx
checkQ Cplx
q
    qsq :: Cplx
qsq = Cplx
q' forall a. Num a => a -> a -> a
* Cplx
q'
    go :: Int -> Cplx -> Cplx -> Cplx -> Cplx -> Cplx
    go :: Int -> Cplx -> Cplx -> Cplx -> Cplx -> Cplx
go Int
n Cplx
out Cplx
alt Cplx
q_2n Cplx
q_n_np1 
      | Int
n forall a. Ord a => a -> a -> Bool
> Int
3000 = forall a. HasCallStack => [Char] -> a
error [Char]
"Reached 3000 iterations."
      | Cplx -> Cplx -> Bool
areClose Cplx
out Cplx
outnew = Cplx
2.0 forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
sqrt (forall a. Floating a => a -> a
sqrt Cplx
q) forall a. Num a => a -> a -> a
* Cplx
out
      | Bool
otherwise = Int -> Cplx -> Cplx -> Cplx -> Cplx -> Cplx
go (Int
n forall a. Num a => a -> a -> a
+ Int
1) Cplx
outnew (-Cplx
alt) Cplx
q_2np1 Cplx
q_np1_np2
        where
          q_2np1 :: Cplx
q_2np1 = Cplx
q_2n forall a. Num a => a -> a -> a
* Cplx
qsq
          q_np1_np2 :: Cplx
q_np1_np2 = Cplx
q_n_np1 forall a. Num a => a -> a -> a
* Cplx
q_2np1
          n' :: Cplx
n' = forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n 
          k :: Cplx
k = Cplx
2.0 forall a. Num a => a -> a -> a
* Cplx
n' forall a. Num a => a -> a -> a
+ Cplx
1.0
          outnew :: Cplx
outnew = Cplx
out forall a. Num a => a -> a -> a
+ Cplx
k forall a. Num a => a -> a -> a
* Cplx
alt forall a. Num a => a -> a -> a
* Cplx
q_np1_np2 forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
cos (Cplx
k forall a. Num a => a -> a -> a
* Cplx
z)