[Scipy-svn] r6941 - in trunk/scipy/special: . specfun tests

[scipy-svn at scipy.org]

Author: ptvirtan
Date: 2010-11-23 15:23:35 -0600 (Tue, 23 Nov 2010)
New Revision: 6941

Modified:
   trunk/scipy/special/basic.py
   trunk/scipy/special/specfun/specfun.f
   trunk/scipy/special/tests/test_basic.py
Log:
BUG: special: initialize variables in specfun.f:JDZO and fix off-by-one error in jnjnp_zeros (#1330)

Initializing the variables avoids garbage output from JDZO.

Modified: trunk/scipy/special/basic.py
===================================================================
--- trunk/scipy/special/basic.py	2010-11-23 21:23:10 UTC (rev 6940)
+++ trunk/scipy/special/basic.py	2010-11-23 21:23:35 UTC (rev 6941)
@@ -74,7 +74,7 @@
         raise ValueError("Number must be integer <= 1200.")
     nt = int(nt)
     n,m,t,zo = specfun.jdzo(nt)
-    return zo[:nt],n[:nt],m[:nt],t[:nt]
+    return zo[1:nt+1],n[:nt],m[:nt],t[:nt]
 
 def jnyn_zeros(n,nt):
     """Compute nt zeros of the Bessel functions Jn(x), Jn'(x), Yn(x), and

Modified: trunk/scipy/special/specfun/specfun.f
===================================================================
--- trunk/scipy/special/specfun/specfun.f	2010-11-23 21:23:10 UTC (rev 6940)
+++ trunk/scipy/special/specfun/specfun.f	2010-11-23 21:23:35 UTC (rev 6941)
@@ -367,6 +367,8 @@
         INTEGER P(1400), P1(70)
         DIMENSION N(1400),M(1400),ZO(0:1400),N1(70),M1(70),
      &            ZOC(0:70),BJ(101),DJ(101),FJ(101)
+        X = 0
+        ZOC(0) = 0
         IF (NT.LT.600) THEN
            XM=-1.0+2.248485*NT**0.5-.0159382*NT+3.208775E-4
      &        *NT**1.5

Modified: trunk/scipy/special/tests/test_basic.py
===================================================================
--- trunk/scipy/special/tests/test_basic.py	2010-11-23 21:23:10 UTC (rev 6940)
+++ trunk/scipy/special/tests/test_basic.py	2010-11-23 21:23:35 UTC (rev 6941)
@@ -8,7 +8,6 @@
 #8   Means test runs forever
 
 ###  test_besselpoly
-###  test_jnjnp_zeros
 ###  test_mathieu_a
 ###  test_mathieu_even_coef
 ###  test_mathieu_odd_coef
@@ -27,7 +26,7 @@
 
 from numpy.testing import assert_equal, assert_almost_equal, assert_array_equal, \
         assert_array_almost_equal, assert_approx_equal, assert_, \
-        rand, dec, TestCase, run_module_suite
+        rand, dec, TestCase, run_module_suite, assert_allclose
 from scipy import special
 import scipy.special._cephes as cephes
 
@@ -1347,11 +1346,18 @@
                                         3101.86438139042]), rtol=1e-8)
 
     def test_jnjnp_zeros(self):
-        pass
-        #jnjp = jnjnp(3)
-        #assert_array_almost_equal(jnjp,(array([
-        #I don't think specfun jdzo is working properly the outputs do not seem to correlate
-        #to the inputs
+        jn = special.jn
+        def jnp(n, x):
+            return (jn(n-1,x) - jn(n+1,x))/2
+        for nt in range(1, 30):
+            z, n, m, t = special.jnjnp_zeros(nt)
+            for zz, nn, tt in zip(z, n, t):
+                if tt == 0:
+                    assert_allclose(jn(nn, zz), 0, atol=1e-6)
+                elif tt == 1:
+                    assert_allclose(jnp(nn, zz), 0, atol=1e-6)
+                else:
+                    raise AssertionError("Invalid t return for nt=%d" % nt)
 
     def test_jnp_zeros(self):
         jnp = special.jnp_zeros(1,5)