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

scipy-svn@scip... scipy-svn@scip...
Sun Jan 25 12:28:13 CST 2009


Author: ptvirtan
Date: 2009-01-25 12:28:00 -0600 (Sun, 25 Jan 2009)
New Revision: 5521

Modified:
   trunk/scipy/special/amos_wrappers.c
   trunk/scipy/special/tests/test_basic.py
Log:
Fixing #853: also I_v: use correct symmetry for negative orders (complex x)

Modified: trunk/scipy/special/amos_wrappers.c
===================================================================
--- trunk/scipy/special/amos_wrappers.c	2009-01-24 10:43:46 UTC (rev 5520)
+++ trunk/scipy/special/amos_wrappers.c	2009-01-25 18:28:00 UTC (rev 5521)
@@ -102,6 +102,24 @@
     return 1;
 }
 
+static int
+reflect_i(Py_complex *ik, double v)
+{
+    if (v != floor(v))
+        return 0;
+    return 1; /* I is symmetric for integer v */
+}
+
+static Py_complex
+rotate_i(Py_complex i, Py_complex k, double v)
+{
+    Py_complex w;
+    double s = sin(v * M_PI)*(2.0/M_PI);
+    w.real = i.real + s*k.real;
+    w.imag = i.imag + s*k.imag;
+    return w;
+}
+
 int cairy_wrap(Py_complex z, Py_complex *ai, Py_complex *aip, Py_complex *bi, Py_complex *bip) {
   int id = 0;
   int ierr = 0;
@@ -142,28 +160,57 @@
 Py_complex cbesi_wrap( double v, Py_complex z) {
   int n = 1;
   int kode = 1;
+  int sign = 1;
   int nz, ierr;
-  Py_complex cy;
+  Py_complex cy, cy_k;
 
   if (v < 0) {
     v = -v;
+    sign = -1;
   }
   F_FUNC(zbesi,ZBESI)(CADDR(z), &v,  &kode, &n, CADDR(cy), &nz, &ierr);
   DO_MTHERR("iv:");
+
+  if (sign == -1) {
+    if (!reflect_i(&cy, v)) {
+      F_FUNC(zbesk,ZBESK)(CADDR(z), &v,  &kode, &n, CADDR(cy_k), &nz, &ierr);
+      DO_MTHERR("iv(kv):");
+      cy = rotate_i(cy, cy_k, v);
+    }
+  }
+      
   return cy;
 }
 
 Py_complex cbesi_wrap_e( double v, Py_complex z) {
   int n = 1;
   int kode = 2;
+  int sign = 1;
   int nz, ierr;
-  Py_complex cy;
+  Py_complex cy, cy_k;
 
   if (v < 0) {
     v = -v;
+    sign = -1;
   }
   F_FUNC(zbesi,ZBESI)(CADDR(z), &v,  &kode, &n, CADDR(cy), &nz, &ierr);
   DO_MTHERR("ive:");
+
+  if (sign == -1) {
+    if (!reflect_i(&cy, v)) {
+      F_FUNC(zbesk,ZBESK)(CADDR(z), &v,  &kode, &n, CADDR(cy_k), &nz, &ierr);
+      DO_MTHERR("ive(kv):");
+      /* adjust scaling to match zbesi */
+      cy_k = rotate(cy_k, -z.imag/M_PI);
+      if (z.real > 0) {
+          cy_k.real *= exp(-2*z.real);
+          cy_k.imag *= exp(-2*z.real);
+      }
+      /* v -> -v */
+      cy = rotate_i(cy, cy_k, v);
+    }
+  }
+
   return cy;
 }
 
@@ -181,6 +228,7 @@
   }
   F_FUNC(zbesj,ZBESJ)(CADDR(z), &v,  &kode, &n, CADDR(cy_j), &nz, &ierr);
   DO_MTHERR("jv:");
+
   if (sign == -1) {
     if (!reflect_jy(&cy_j, v)) {
       F_FUNC(zbesy,ZBESY)(CADDR(z), &v,  &kode, &n, CADDR(cy_y), &nz, CADDR(cwork), &ierr);
@@ -270,6 +318,7 @@
   Py_complex cy;
 
   if (v < 0) {
+    /* K_v == K_{-v} even for non-integer v */
     v = -v;
   }
   F_FUNC(zbesk,ZBESK)(CADDR(z), &v,  &kode, &n, CADDR(cy), &nz, &ierr);
@@ -284,6 +333,7 @@
   Py_complex cy;
 
   if (v < 0) {
+    /* K_v == K_{-v} even for non-integer v */
     v = -v;
   }
   F_FUNC(zbesk,ZBESK)(CADDR(z), &v, &kode, &n, CADDR(cy), &nz, &ierr);

Modified: trunk/scipy/special/tests/test_basic.py
===================================================================
--- trunk/scipy/special/tests/test_basic.py	2009-01-24 10:43:46 UTC (rev 5520)
+++ trunk/scipy/special/tests/test_basic.py	2009-01-25 18:28:00 UTC (rev 5521)
@@ -1650,25 +1650,45 @@
         assert_tol_equal(jv(301, 1296.0682), -0.0224174325312048)
 
     def test_ticket_853(self):
+        """Negative-order Bessels"""
         # cephes
         assert_tol_equal(jv(-1,   1   ), -0.4400505857449335)
         assert_tol_equal(jv(-2,   1   ), 0.1149034849319005)
         assert_tol_equal(yv(-1,   1   ), 0.7812128213002887)
         assert_tol_equal(yv(-2,   1   ), -1.650682606816255)
+        assert_tol_equal(iv(-1,   1   ), 0.5651591039924851)
+        assert_tol_equal(iv(-2,   1   ), 0.1357476697670383)
+        assert_tol_equal(kv(-1,   1   ), 0.6019072301972347)
+        assert_tol_equal(kv(-2,   1   ), 1.624838898635178)
         assert_tol_equal(jv(-0.5, 1   ), 0.43109886801837607952)
         assert_tol_equal(yv(-0.5, 1   ), 0.6713967071418031)
+        assert_tol_equal(iv(-0.5, 1   ), 1.231200214592967)
+        assert_tol_equal(kv(-0.5, 1   ), 0.4610685044478945)
         # amos
         assert_tol_equal(jv(-1,   1+0j), -0.4400505857449335)
         assert_tol_equal(jv(-2,   1+0j), 0.1149034849319005)
         assert_tol_equal(yv(-1,   1+0j), 0.7812128213002887)
         assert_tol_equal(yv(-2,   1+0j), -1.650682606816255)
+
+        assert_tol_equal(iv(-1,   1+0j), 0.5651591039924851)
+        assert_tol_equal(iv(-2,   1+0j), 0.1357476697670383)
+        assert_tol_equal(kv(-1,   1+0j), 0.6019072301972347)
+        assert_tol_equal(kv(-2,   1+0j), 1.624838898635178)
+        
         assert_tol_equal(jv(-0.5, 1+0j), 0.43109886801837607952)
         assert_tol_equal(jv(-0.5, 1+1j), 0.2628946385649065-0.827050182040562j)
         assert_tol_equal(yv(-0.5, 1+0j), 0.6713967071418031)
         assert_tol_equal(yv(-0.5, 1+1j), 0.967901282890131+0.0602046062142816j)
+        
+        assert_tol_equal(iv(-0.5, 1+0j), 1.231200214592967)
+        assert_tol_equal(iv(-0.5, 1+1j), 0.77070737376928+0.39891821043561j)
+        assert_tol_equal(kv(-0.5, 1+0j), 0.4610685044478945)
+        assert_tol_equal(kv(-0.5, 1+1j), 0.06868578341999-0.38157825981268j)
 
-        assert_tol_equal(jve(-0.5,1+0j), jv(-0.5, 1+0j))
-        assert_tol_equal(yve(-0.5,1+0j), yv(-0.5, 1+0j))
+        assert_tol_equal(jve(-0.5,1+0.3j), jv(-0.5, 1+0.3j)*exp(-0.3))
+        assert_tol_equal(yve(-0.5,1+0.3j), yv(-0.5, 1+0.3j)*exp(-0.3))
+        assert_tol_equal(ive(-0.5,0.3+1j), iv(-0.5, 0.3+1j)*exp(-0.3))
+        assert_tol_equal(kve(-0.5,0.3+1j), kv(-0.5, 0.3+1j)*exp(0.3+1j))
 
         assert_tol_equal(hankel1(-0.5, 1+1j), jv(-0.5, 1+1j) + 1j*yv(-0.5,1+1j))
         assert_tol_equal(hankel2(-0.5, 1+1j), jv(-0.5, 1+1j) - 1j*yv(-0.5,1+1j))



More information about the Scipy-svn mailing list