diff --git a/CHANGELOG.md b/CHANGELOG.md
new file mode 100644
--- /dev/null
+++ b/CHANGELOG.md
@@ -0,0 +1,5 @@
+# Changelog for `jacobi-elliptic`
+
+## 0.1.0.0 - 2023-02-20
+
+First release.
diff --git a/LICENSE b/LICENSE
new file mode 100644
--- /dev/null
+++ b/LICENSE
@@ -0,0 +1,30 @@
+Copyright Stéphane Laurent (c) 2023
+
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are met:
+
+    * Redistributions of source code must retain the above copyright
+      notice, this list of conditions and the following disclaimer.
+
+    * Redistributions in binary form must reproduce the above
+      copyright notice, this list of conditions and the following
+      disclaimer in the documentation and/or other materials provided
+      with the distribution.
+
+    * Neither the name of Stéphane Laurent nor the names of other
+      contributors may be used to endorse or promote products derived
+      from this software without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
+OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
diff --git a/README.md b/README.md
new file mode 100644
--- /dev/null
+++ b/README.md
@@ -0,0 +1,8 @@
+# jacobi-elliptic
+
+<!-- badges: start -->
+[![Stack-lts](https://github.com/stla/jacobi-elliptic/actions/workflows/Stack-lts.yml/badge.svg)](https://github.com/stla/jacobi-elliptic/actions/workflows/Stack-lts.yml)
+[![Stack-nightly](https://github.com/stla/jacobi-elliptic/actions/workflows/Stack-nightly.yml/badge.svg)](https://github.com/stla/jacobi-elliptic/actions/workflows/Stack-nightly.yml)
+<!-- badges: end -->
+
+Evaluation of the Neville theta functions and the Jacobi elliptic functions.
diff --git a/Setup.hs b/Setup.hs
new file mode 100644
--- /dev/null
+++ b/Setup.hs
@@ -0,0 +1,2 @@
+import Distribution.Simple
+main = defaultMain
diff --git a/jacobi-elliptic.cabal b/jacobi-elliptic.cabal
new file mode 100644
--- /dev/null
+++ b/jacobi-elliptic.cabal
@@ -0,0 +1,48 @@
+name:                jacobi-elliptic
+version:             0.1.0.0
+synopsis:            Neville Theta Functions and Jacobi Elliptic Functions
+description:         Evaluation of the Neville theta functions and the Jacobi elliptic functions.
+homepage:            https://github.com/stla/jacobi-elliptic#readme
+license:             BSD3
+license-file:        LICENSE
+author:              Stéphane Laurent
+maintainer:          laurent_step@outlook.fr
+copyright:           2023 Stéphane Laurent
+category:            Math, Numeric
+build-type:          Simple
+extra-source-files:  README.md
+                     CHANGELOG.md
+cabal-version:       >=1.10
+
+library
+  hs-source-dirs:      src
+  exposed-modules:     Math.NevilleTheta
+                     , Math.JacobiElliptic
+  build-depends:       base >= 4.7 && < 5
+                     , jacobi-theta >= 0.1.0.0
+                     , elliptic-integrals >= 0.1.0.0
+  default-language:    Haskell2010
+  ghc-options:         -Wall
+                       -Wcompat
+                       -Widentities
+                       -Wincomplete-record-updates
+                       -Wincomplete-uni-patterns
+                       -Wmissing-export-lists
+                       -Wmissing-home-modules
+                       -Wpartial-fields
+                       -Wredundant-constraints
+
+test-suite unit-tests
+  type:                 exitcode-stdio-1.0
+  main-is:              Main.hs
+  hs-source-dirs:       tests/
+  other-modules:        Approx
+  Build-Depends:        base >= 4.7 && < 5
+                      , tasty
+                      , tasty-hunit
+                      , jacobi-elliptic
+  Default-Language:     Haskell2010
+
+source-repository head
+  type:     git
+  location: https://github.com/stla/jacobi-elliptic
diff --git a/src/Math/JacobiElliptic.hs b/src/Math/JacobiElliptic.hs
new file mode 100644
--- /dev/null
+++ b/src/Math/JacobiElliptic.hs
@@ -0,0 +1,62 @@
+module Math.JacobiElliptic
+    ( jellip,
+      jellip'
+    ) where
+import Data.Complex ( Complex )
+import Math.NevilleTheta
+    ( theta_c,
+      theta_d,
+      theta_n,
+      theta_s,
+      theta_c',
+      theta_d',
+      theta_n',
+      theta_s' )
+
+type Cplx = Complex Double
+
+-- | Jacobi elliptic function in terms of the nome.
+jellip :: 
+     Char -- ^ a letter among 'c', 'd', 'n', 's' identifying the Neville function at the numerator
+  -> Char -- ^ a letter among 'c', 'd', 'n', 's' identifying the Neville function at the denominator
+  -> Cplx -- ^ z, the variable
+  -> Cplx -- ^ q, the nome
+  -> Cplx
+jellip p q z nome = 
+  theta_num z nome / theta_den z nome
+  where
+    theta_num = case p of
+      'c' -> theta_c
+      'd' -> theta_d
+      'n' -> theta_n
+      's' -> theta_s
+      _   -> error "Invalid numerator identifier."
+    theta_den = case q of
+      'c' -> theta_c
+      'd' -> theta_d
+      'n' -> theta_n
+      's' -> theta_s
+      _   -> error "Invalid denominator identifier."
+
+-- | Jacobi elliptic function in terms of the squared modulus.
+jellip' :: 
+     Char -- ^ a letter among 'c', 'd', 'n', 's' identifying the Neville function at the numerator
+  -> Char -- ^ a letter among 'c', 'd', 'n', 's' identifying the Neville function at the denominator
+  -> Cplx -- ^ z, the variable
+  -> Cplx -- ^ m, the squared modulus
+  -> Cplx
+jellip' p q z m = 
+  theta_num z m / theta_den z m
+  where
+    theta_num = case p of
+      'c' -> theta_c'
+      'd' -> theta_d'
+      'n' -> theta_n'
+      's' -> theta_s'
+      _   -> error "Invalid numerator identifier."
+    theta_den = case q of
+      'c' -> theta_c'
+      'd' -> theta_d'
+      'n' -> theta_n'
+      's' -> theta_s'
+      _   -> error "Invalid denominator identifier."
diff --git a/src/Math/NevilleTheta.hs b/src/Math/NevilleTheta.hs
new file mode 100644
--- /dev/null
+++ b/src/Math/NevilleTheta.hs
@@ -0,0 +1,98 @@
+module Math.NevilleTheta
+    ( theta_c, 
+      theta_d,
+      theta_n,
+      theta_s,
+      theta_c', 
+      theta_d',
+      theta_n',
+      theta_s'
+    ) where
+import Data.Complex           ( Complex(..) )
+import Math.EllipticIntegrals ( ellipticF )
+import Math.JacobiTheta
+    ( jtheta1, jtheta1Dash, jtheta2, jtheta3, jtheta4 )
+
+type Cplx = Complex Double
+
+i_ :: Cplx
+i_ = 0.0 :+ 1.0
+
+tauFromM :: Cplx -> Cplx
+tauFromM m = i_ * ellipticF (pi/2) (1 - m) / ellipticF (pi/2) m
+
+nomeFromM :: Cplx -> Cplx
+nomeFromM m = exp (i_ * pi * tauFromM m)
+
+-- | Neville theta-c function in terms of the nome.
+theta_c :: 
+     Cplx -- ^ z
+  -> Cplx -- ^ q, the nome
+  -> Cplx
+theta_c z q = 
+  jtheta2 z' q / jtheta2 0 q
+  where
+    j3 = jtheta3 0 q
+    z' = z / (j3 * j3)
+
+-- | Neville theta-d function in terms of the nome.
+theta_d :: 
+     Cplx -- ^ z
+  -> Cplx -- ^ q, the nome
+  -> Cplx
+theta_d z q = 
+  jtheta3 z' q / jtheta3 0 q
+  where
+    j3 = jtheta3 0 q
+    z' = z / (j3 * j3)
+
+-- | Neville theta-n function in terms of the nome.
+theta_n :: 
+     Cplx -- ^ z
+  -> Cplx -- ^ q, the nome
+  -> Cplx
+theta_n z q = 
+  jtheta4 z' q / jtheta4 0 q
+  where
+    j3 = jtheta3 0 q
+    z' = z / (j3 * j3)
+
+-- | Neville theta-d function in terms of the nome.
+theta_s :: 
+     Cplx -- ^ z
+  -> Cplx -- ^ q, the nome
+  -> Cplx
+theta_s z q = 
+  j3sq * jtheta1 z' q / jtheta1Dash 0 q
+  where
+    j3 = jtheta3 0 q
+    j3sq = j3 * j3
+    z' = z / j3sq
+
+-- | Neville theta-c function in terms of the squared modulus.
+theta_c' :: 
+     Cplx -- ^ z
+  -> Cplx -- ^ m, the squared modulus
+  -> Cplx
+theta_c' z m = theta_c z (nomeFromM m)
+
+-- | Neville theta-d function in terms of the squared modulus.
+theta_d' :: 
+     Cplx -- ^ z
+  -> Cplx -- ^ m, the squared modulus
+  -> Cplx
+theta_d' z m = theta_d z (nomeFromM m)
+
+-- | Neville theta-n function in terms of the squared modulus.
+theta_n' :: 
+     Cplx -- ^ z
+  -> Cplx -- ^ m, the squared modulus
+  -> Cplx
+theta_n' z m = theta_n z (nomeFromM m)
+
+-- | Neville theta-s function in terms of the squared modulus.
+theta_s' :: 
+     Cplx -- ^ z
+  -> Cplx -- ^ m, the squared modulus
+  -> Cplx
+theta_s' z m = theta_s z (nomeFromM m)
diff --git a/tests/Approx.hs b/tests/Approx.hs
new file mode 100644
--- /dev/null
+++ b/tests/Approx.hs
@@ -0,0 +1,15 @@
+module Approx (assertApproxEqual) where
+import           Data.Complex     ( imagPart, realPart, Complex(..) )
+import           Test.Tasty.HUnit ( Assertion, assertEqual )
+
+-- round x to n digits
+approx0 :: Int -> Double -> Double
+approx0 n x = fromInteger (round $ x * (10^n)) / (10.0^^n)
+
+-- round z to n digits
+approx :: Int -> Complex Double -> Complex Double
+approx n z = approx0 n (realPart z) :+ approx0 n (imagPart z)
+
+assertApproxEqual :: String -> Int -> Complex Double -> Complex Double -> Assertion
+assertApproxEqual prefix n z1 z2 = 
+  assertEqual prefix (approx n z1) (approx n z2)
diff --git a/tests/Main.hs b/tests/Main.hs
new file mode 100644
--- /dev/null
+++ b/tests/Main.hs
@@ -0,0 +1,140 @@
+module Main where
+import           Approx               ( assertApproxEqual )
+import           Data.Complex         ( Complex(..) )
+import           Math.NevilleTheta
+                                      ( theta_c,
+                                        theta_d,
+                                        theta_n,
+                                        theta_s,
+                                        theta_c',
+                                        theta_d',
+                                        theta_n',
+                                        theta_s' )
+import           Math.JacobiElliptic  ( jellip' )
+import           Test.Tasty           ( defaultMain, testGroup )
+import           Test.Tasty.HUnit     ( testCase )
+
+i_ :: Complex Double
+i_ = 0.0 :+ 1.0
+
+z :: Complex Double
+z = 1.0 :+ 1.0
+
+q :: Complex Double 
+q = exp (-pi)
+
+q' :: Complex Double 
+q' = exp (-pi/10)
+
+q'' :: Complex Double 
+q'' = exp (i_ * pi * tau)
+  where
+    tau = 2.0 :+ 2.0
+
+u :: Complex Double
+u = 0.3 :+ 0.7
+
+m :: Complex Double
+m = 0.4 :+ 0.0
+
+
+main :: IO ()
+main = defaultMain $
+  testGroup "Tests"
+  [ 
+    testCase "theta_c value 1" $ do
+      let expected = 0.902705416117337 :+ (-0.718974020880116)
+          obtained = theta_c z q
+      assertApproxEqual "" 10 expected obtained,
+
+    testCase "theta_c value 2" $ do
+      let expected = 0.997974260633626 :+ (-0.063618983904188)
+          obtained = theta_c z q'
+      assertApproxEqual "" 10 expected obtained,
+
+    testCase "theta_c value 3" $ do
+      let expected = 0.838567437919619 :+ (-0.974584266572289)
+          obtained = theta_c z q''
+      assertApproxEqual "" 10 expected obtained,
+
+    testCase "theta_d value 1" $ do
+      let expected = 0.892748081976972 :+ (-0.207593861225047)
+          obtained = theta_d z q
+      assertApproxEqual "" 10 expected obtained,
+
+    testCase "theta_d value 2" $ do
+      let expected = 0.997974260633412 :+ (-0.063618983903874)
+          obtained = theta_d z q'
+      assertApproxEqual "" 10 expected obtained,
+
+    testCase "theta_d value 3" $ do
+      let expected = 0.990723180697351 :+ (-0.012164484951676)
+          obtained = theta_d z q''
+      assertApproxEqual "" 10 expected obtained,
+
+    testCase "theta_n value 1" $ do
+      let expected = 1.12730988168993 :+ 0.2469274015421
+          obtained = theta_n z q
+      assertApproxEqual "" 10 expected obtained,
+
+    testCase "theta_n value 2" $ do
+      let expected = 0.894953772623932 :+ 0.933853399701569
+          obtained = theta_n z q'
+      assertApproxEqual "" 10 expected obtained,
+
+    testCase "theta_n value 3" $ do
+      let expected = 1.00934637387594 :+ 0.01225569246714
+          obtained = theta_n z q''
+      assertApproxEqual "" 10 expected obtained,
+
+    testCase "theta_s value 1" $ do
+      let expected = 1.22039326540444 :+ 0.75990704701835
+          obtained = theta_s z q
+      assertApproxEqual "" 10 expected obtained,
+
+    testCase "theta_s value 2" $ do
+      let expected = 0.7162841953585 :+ 1.25543148570321
+          obtained = theta_s z q'
+      assertApproxEqual "" 10 expected obtained,
+
+    testCase "theta_s value 3" $ do
+      let expected = 1.29457805665579 :+ 0.64084576896851
+          obtained = theta_s z q''
+      assertApproxEqual "" 10 expected obtained,
+
+    testCase "a value of theta_c prime" $ do
+      let expected = -0.65900466676738154967
+          obtained = theta_c' 2.5 0.3
+      assertApproxEqual "" 15 expected obtained,
+
+    testCase "a value of theta_d prime" $ do
+      let expected = 0.95182196661267561994
+          obtained = theta_d' 2.5 0.3
+      assertApproxEqual "" 15 expected obtained,
+
+    testCase "a value of theta_n prime" $ do
+      let expected = 1.0526693354651613637
+          obtained = theta_n' 2.5 0.3
+      assertApproxEqual "" 14 expected obtained,
+
+    testCase "a value of theta_s prime" $ do
+      let expected = 0.82086879524530400536
+          obtained = theta_s' 2.5 0.3
+      assertApproxEqual "" 14 expected obtained,
+
+    testCase "jellip relation 1" $ do
+      let z1 = jellip' 'c' 'n' u m 
+          z2 = jellip' 'n' 'c' (i_ * u) (1 - m) 
+      assertApproxEqual "" 14 z1 z2, 
+
+    testCase "jellip relation 2" $ do
+      let z1 = jellip' 's' 'n' u m 
+          z2 = -i_ * jellip' 's' 'c' (i_ * u) (1 - m) 
+      assertApproxEqual "" 14 z1 z2, 
+
+    testCase "jellip relation 3" $ do
+      let z1 = jellip' 'd' 'n' u m 
+          z2 = jellip' 'd' 'c' (i_ * u) (1 - m) 
+      assertApproxEqual "" 14 z1 z2
+
+  ]