Mercurial > octave-nkf

changeset 19234:068a3e51b7b8

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
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.
author Mike Miller <mtmiller@ieee.org>
date Wed, 01 Oct 2014 22:39:04 -0400
parents 3a6fd52e1458
children 566ea8a8683b
files libinterp/corefcn/ellipj.cc liboctave/numeric/lo-specfun.cc
diffstat 2 files changed, 10 insertions(+), 2 deletions(-) [+]
line wrap: on
line diff
--- 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)
--- 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);
     }
 }