# HG changeset patch
# User Mike Miller <mtmiller@ieee.org>
# Date 1412217544 14400
# Node ID 068a3e51b7b834eeafbfa962199eb0ffee3ec1dc
# Parent  3a6fd52e1458d4f6876bd7a0b56720c1a04fc3e3
ellipj: Fix continuity of dn output when cn is near zero (bug #43344)

* lo-specfun.cc (ellipj): Fix continuity of the dn output, the current
method results in a discontinuity when cn is near zero.
* ellipj.cc: Add test for bug #43344.

diff -r 3a6fd52e1458 -r 068a3e51b7b8 libinterp/corefcn/ellipj.cc
--- a/libinterp/corefcn/ellipj.cc	Wed Oct 01 19:36:05 2014 +0100
+++ b/libinterp/corefcn/ellipj.cc	Wed Oct 01 22:39:04 2014 -0400
@@ -919,6 +919,15 @@
 %! assert (cn, C, 8*eps);
 %! assert (dn, D, 8*eps);
 
+%!test
+%! ## Test continuity of dn when cn is near zero (bug #43344)
+%! m = 0.5;
+%! u = ellipke (0.5);
+%! x = [-1e-3, -1e-12, 0, 1e-12, 1e-3];
+%! [~, ~, dn] = ellipj (u + x, m);
+%! D = 1/sqrt (2) * ones (size (x));
+%! assert (dn, D, 1e-6);
+
 %!error ellipj ()
 %!error ellipj (1)
 %!error ellipj (1,2,3,4)
diff -r 3a6fd52e1458 -r 068a3e51b7b8 liboctave/numeric/lo-specfun.cc
--- a/liboctave/numeric/lo-specfun.cc	Wed Oct 01 19:36:05 2014 +0100
+++ b/liboctave/numeric/lo-specfun.cc	Wed Oct 01 22:39:04 2014 -0400
@@ -3680,12 +3680,11 @@
       phi = ii*a[Nn]*u;
       for (n = Nn; n > 0; --n)
         {
-          t = phi;
           phi = (asin ((c[n]/a[n])* sin (phi)) + phi)/2;
         }
       sn = sin (phi);
       cn = cos (phi);
-      dn = cn/cos (t - phi);
+      dn = sqrt (1 - m*sn*sn);
     }
 }