diff --git a/build.rs b/build.rs
index e52ab87..8b5daac 100644
--- a/build.rs
+++ b/build.rs
@@ -104,7 +104,7 @@ const XSF_TYPES: &[(&str, &str)] = &[
     ("digamma", "d->d"),
     ("digamma", "D->D"),
     // ellip.h
-    //  TODO: ellipj
+    ("ellipj", "dd->dddd"),
     ("ellipk", "d->d"),
     ("ellipkm1", "d->d"),
     ("ellipkinc", "dd->d"),
diff --git a/src/ellip.rs b/src/ellip.rs
index 1e4b0fd..52671bb 100644
--- a/src/ellip.rs
+++ b/src/ellip.rs
@@ -1,3 +1,4 @@
+use crate::bindings;
 use crate::bindings::xsf_impl;
 
 xsf_impl!(ellipk, (m: f64), "Complete elliptic integral of the first kind");
@@ -6,3 +7,29 @@ xsf_impl!(ellipkinc, (phi: f64, m: f64), "Incomplete elliptic integral of the fi
 
 xsf_impl!(ellipe, (m: f64), "Complete elliptic integral of the second kind");
 xsf_impl!(ellipeinc, (phi: f64, m: f64), "Incomplete elliptic integral of the second kind");
+
+/// Jacobian elliptic functions
+///
+/// # Arguments
+/// - `u` - Real argument
+/// - `m` - Parameter between 0 and 1
+///
+/// # Returns
+/// - *sn(u | m)* - sine amplitude
+/// - *cn(u | m)* - cosine amplitude
+/// - *dn(u | m)* - delta amplitude
+/// - phase *φ* s.t. *sn(u | m) = sin(φ)* and *cn(u | m) = cos(φ)*
+///
+/// # See Also
+/// - [`ellipk`] - Complete elliptic integral of the first kind
+/// - [`ellipkinc`] - Incomplete elliptic integral of the first kind
+pub fn ellipj(u: f64, m: f64) -> (f64, f64, f64, f64) {
+    let mut sn = f64::NAN;
+    let mut cn = f64::NAN;
+    let mut dn = f64::NAN;
+    let mut phi = f64::NAN;
+    unsafe {
+        bindings::ellipj(u, m, &mut sn, &mut cn, &mut dn, &mut phi);
+    }
+    (sn, cn, dn, phi)
+}
diff --git a/src/lib.rs b/src/lib.rs
index 5c80660..37c6bba 100644
--- a/src/lib.rs
+++ b/src/lib.rs
@@ -37,7 +37,7 @@ pub use digamma::digamma;
 
 // ellip.h
 mod ellip;
-pub use ellip::{ellipe, ellipeinc, ellipk, ellipkinc, ellipkm1};
+pub use ellip::{ellipe, ellipeinc, ellipj, ellipk, ellipkinc, ellipkm1};
 
 // erf.h
 mod erf;
diff --git a/tests/test_functions.rs b/tests/test_functions.rs
index 233eb69..8fd0b4c 100644
--- a/tests/test_functions.rs
+++ b/tests/test_functions.rs
@@ -615,6 +615,17 @@ macro_rules! xsref_test {
             );
         }
     };
+    (@single $f:ident, "dd->dddd") => {
+        paste::paste! {
+            _test!(
+                [<test_ $f _ d>],
+                $f,
+                "d_d-d_d_d_d",
+                |x: &[f64]| xsf::$f(x[0], x[1]),
+                (f64, f64, f64, f64)
+            );
+        }
+    };
     (@single $f:ident, "dd->dd") => {
         paste::paste! {
             _test!(
@@ -817,6 +828,7 @@ xsref_test!(ellipkm1, "d->d");
 xsref_test!(ellipkinc, "dd->d");
 xsref_test!(ellipe, "d->d");
 xsref_test!(ellipeinc, "dd->d");
+xsref_test!(ellipj, "dd->dddd");
 
 // erf.h
 xsref_test!(erf, "d->d", "D->D");