new bessel caclulation

bessel.cc

@@ -24,38 +24,355 @@
 //    You should have received a copy of the GNU General Public License             //
 //    along with this program.  If not, see <http://www.gnu.org/licenses/>.         //
 //**********************************************************************************//
+#include "bessel.h"
+#include <algorithm>
 #include <complex>
 #include <cmath>
+#include <limits>
 #include <stdexcept>
 #include <vector>
 
 namespace nmie {
-  // Implementation of Bessel functions from Bohren and Huffman book, pp. 86-87, eq 4.11
-  // Calculate all orders of function from 0 to nmax (included) for argument rho
-  std::vector< std::complex<double> > bessel_j(int nmax, std::complex<double> rho) {
-    if (nmax < 0) throw std::invalid_argument("Bessel order should be >= 0 (nmie::bessel_j)\n");
-    std::vector< std::complex<double> > j(nmax+1);
-    j[0] = std::sin(rho)/rho;
-    if (nmax == 0) return j;
-    j[1] = std::sin(rho)/(rho*rho) - std::cos(rho)/rho;
-    if (nmax == 1) return j;
-    for (int i = 2; i < n+1; ++i) {
-      int n = i - 1;
-      j[n+1] = static_cast<double>(2*n+1)/rho*j[n] - j[n-1];
+  namespace bessel {
+    // Implementation of spherical Bessel functions from
+    // L.-W. Cai / Computer Physics Communications 182 (2011) 663–668
+    std::vector< std::complex<double> > bessel_j(int nmax, std::complex<double> rho) {
+      if (nmax < 0) throw std::invalid_argument("Bessel order should be >= 0 (nmie::bessel_j)\n");
+      std::complex<double> c_zero(0.0,0.0);
+      std::vector< std::complex<double> > j(nmax+1, c_zero);
+      // Select normalization amplitude
+      int norm_n = 0;
+      std::complex<double> norm_value = std::sin(rho)/rho,
+	tmp = std::sin(rho)/(rho*rho) - std::cos(rho)/rho;
+      if (std::abs(tmp) > std::abs(norm_value)) {
+	norm_value = tmp;
+	norm_n = 1;
+      }
+      // calculate number of terms for backward recursion Cai eq(11)
+      double theta = std::arg(rho), r = std::abs(rho);
+      double M = (1.83 + 4.1* std::pow(std::sin(theta), 0.36)) *
+	std::pow(r, (0.91 - 0.43*std::pow(std::sin(theta),0.33)))
+	+ 9.0*(1.0-std::sqrt(std::sin(theta)));
+      int terms = std::max(nmax+2, static_cast<int>(std::ceil(M))+2);
+      printf("terms = %d\n", terms);
+      // Seed for eq(2)
+      std::complex<double> b_np1(0.0, 0.0), b_nm1;
+      double eps = std::numeric_limits<double>::epsilon();
+      std::complex<double> b_n(eps, 0.0);
+      //recurence
+      for (int n = terms*2; n > 0; --n) {
+	b_nm1 = (2.0*static_cast<double>(n)+1.0)/rho*b_n - b_np1;
+	if (n-1 < nmax+1) j[n-1] = b_nm1;
+      }
+      // normalize
+      norm_value /= j[norm_n];
+      for (auto& value : j) value *=norm_value;
+      
+      return j;
     }
-  }
-  // Implementation of Bessel functions from Bohren and Huffman book, pp. 86-87, eq 4.11
-  // Calculate all orders of function from 0 to nmax (included) for argument rho
-  std::vector< std::complex<double> > bessel_y(int nmax, std::complex<double> rho) {
-    if (nmax < 0) throw std::invalid_argument("Bessel order should be >= 0 (nmie::bessel_j)\n");
-    std::vector< std::complex<double> > j(nmax+1);
-    j[0] = -std::cos(rho)/rho;
-    if (nmax == 0) return j;
-    j[1] = -std::cos(rho)/(rho*rho) - std::sin(rho)/rho;
-    if (nmax == 1) return j;
-    for (int i = 2; i < n+1; ++i) {
-      int n = i - 1;
-      j[n+1] = static_cast<double>(2*n+1)/rho*j[n] - j[n-1];
+    // Implementation of Bessel functions from Bohren and Huffman book, pp. 86-87, eq 4.11
+    // Calculate all orders of function from 0 to nmax (included) for argument rho
+    std::vector< std::complex<double> > bh_bessel_j(int nmax, std::complex<double> rho) {
+      if (nmax < 0) throw std::invalid_argument("Bessel order should be >= 0 (nmie::bessel_j)\n");
+      std::vector< std::complex<double> > j(nmax+1);
+      j[0] = std::sin(rho)/rho;
+      if (nmax == 0) return j;
+      j[1] = std::sin(rho)/(rho*rho) - std::cos(rho)/rho;
+      if (nmax == 1) return j;
+      for (int i = 2; i < nmax+1; ++i) {
+	int n = i - 1;
+	j[n+1] = static_cast<double>(2*n+1)/rho*j[n] - j[n-1];
+      }
+      return j;
     }
-  }
+
+
+// !*****************************************************************************80
+//  
+//  C++ port of fortran code
+//  
+// !! CSPHJY: spherical Bessel functions jn(z) and yn(z) for complex argument.
+// !
+// !  Discussion:
+// !
+// !    This procedure computes spherical Bessel functions jn(z) and yn(z)
+// !    and their derivatives for a complex argument.
+// !
+// !  Licensing:
+// !
+// !    This routine is copyrighted by Shanjie Zhang and Jianming Jin.  However, 
+// !    they give permission to incorporate this routine into a user program 
+// !    provided that the copyright is acknowledged.
+// !
+// !  Modified:
+// !
+// !    01 August 2012
+// !
+// !  Author:
+// !
+// !    Shanjie Zhang, Jianming Jin
+// !
+// !  Reference:
+// !
+// !    Shanjie Zhang, Jianming Jin,
+// !    Computation of Special Functions,
+// !    Wiley, 1996,
+// !    ISBN: 0-471-11963-6,
+// !    LC: QA351.C45.
+// ! 
+// !  Parameters:
+// !
+// !    Input, integer ( kind = 4 ) N, the order of jn(z) and yn(z).
+// !
+// !    Input, complex ( kind = 8 ) Z, the argument.
+// !
+// !    Output, integer ( kind = 4 ) NM, the highest order computed.
+// !
+// !    Output, complex ( kind = 8 ) CSJ(0:N0, CDJ(0:N), CSY(0:N), CDY(0:N),
+// !    the values of jn(z), jn'(z), yn(z), yn'(z).
+// !
+    void csphjy (int n, std::complex<double>z, int& nm,
+		 std::vector< std::complex<double> >& csj,
+		 std::vector< std::complex<double> >& cdj,
+		 std::vector< std::complex<double> >& csy,
+		 std::vector< std::complex<double> >& cdy ) {
+      double a0;
+      csj.resize(n+1);
+      cdj.resize(n+1);
+      csy.resize(n+1);
+      cdy.resize(n+1);
+      std::complex<double> cf, cf0, cf1, cs, csa, csb;
+      int k, m;
+      a0 = std::abs(z);
+      nm = n;      
+      if (a0 < 1.0e-60) {
+	for (int k = 0; k < n+1; ++k) {
+	  csj[k] = 0.0;
+	  cdj[k] = 0.0;
+	  csy[k] = -1.0e+300;
+	  cdy[k] = 1.0e+300;
+	}
+	csj[0] = std::complex<double>( 1.0, 0.0);
+	cdj[1] =  std::complex<double>( 0.333333333333333, 0.0);
+	return;
+      }
+      csj[0] = std::sin ( z ) / z;
+      csj[1] = ( csj[0] - std::cos ( z ) ) / z;
+      
+      if ( 2 <= n ) {
+	csa = csj[0];
+	csb = csj[1];
+	m = msta1 ( a0, 200 );
+	if ( m < n ) nm = m;
+	else m = msta2 ( a0, n, 15 );
+	cf0 = 0.0;
+	cf1 = 1.0e-100;
+	for (int k = m; k>=0; --k) {
+	  cf = ( 2.0 * k + 3.0 ) * cf1 / z - cf0;
+	  if ( k <= nm ) csj[k] = cf;
+	  cf0 = cf1;
+	  cf1 = cf;
+	}
+	if ( std::abs ( csa ) <= std::abs ( csb ) ) cs = csb / cf0;
+	else  cs = csa / cf;
+	for (int k = 0; k <= nm; ++k) {
+	  csj[k] = cs * csj[k];
+	}
+      }
+      cdj[0] = ( std::cos ( z ) - std::sin ( z ) / z ) / z;
+      for (int k = 1; k <=nm; ++k) {
+	cdj[k] = csj[k-1] - ( k + 1.0 ) * csj[k] / z;
+      }
+      csy[0] = - std::cos ( z ) / z;
+      csy[1] = ( csy[0] - std::sin ( z ) ) / z;
+      cdy[0] = ( std::sin ( z ) + std::cos ( z ) / z ) / z;
+      cdy[1] = ( 2.0 * cdy[0] - std::cos ( z ) )  / z;
+      for (int k = 2; k<=nm; ++k) {
+	if ( std::abs ( csj[k-2] ) < std::abs ( csj[k-1] ) ) {
+	  csy[k] = ( csj[k] * csy[k-1] - 1.0 / ( z * z ) ) / csj[k-1];
+	} else {
+	  csy[k] = ( csj[k] * csy[k-2] - ( 2.0 * k - 1.0 ) / std::pow(z,3) )
+	    / csj[k-2];
+	}
+      }
+      for (int k = 2; k<=nm; ++k) {
+	cdy[k] = csy[k-1] - ( k + 1.0 ) * csy[k] / z;
+      }
+      
+      return;
+    }
+      // function msta2 ( x, n, mp )
+      
+      // !*****************************************************************************80
+      // !
+      // !! MSTA2 determines a backward recurrence starting point for Jn(x).
+      // !
+      // !  Discussion:
+      // !
+      // !    This procedure determines the starting point for a backward
+      // !    recurrence such that all Jn(x) has MP significant digits.
+      // !
+      // !  Licensing:
+      // !
+      // !    This routine is copyrighted by Shanjie Zhang and Jianming Jin.  However, 
+      // !    they give permission to incorporate this routine into a user program 
+      // !    provided that the copyright is acknowledged.
+      // !
+      // !  Modified:
+      // !
+      // !    08 July 2012
+      // !
+      // !  Author:
+      // !
+      // !    Shanjie Zhang, Jianming Jin
+      // !
+      // !  Reference:
+      // !
+      // !    Shanjie Zhang, Jianming Jin,
+      // !    Computation of Special Functions,
+      // !    Wiley, 1996,
+      // !    ISBN: 0-471-11963-6,
+      // !    LC: QA351.C45.
+      // !
+      // !  Parameters:
+      // !
+      // !    Input, real ( kind = 8 ) X, the argument of Jn(x).
+      // !
+      // !    Input, integer ( kind = 4 ) N, the order of Jn(x).
+      // !
+      // !    Input, integer ( kind = 4 ) MP, the number of significant digits.
+      // !
+      // !    Output, integer ( kind = 4 ) MSTA2, the starting point.
+      // !
+    int msta2 ( double x, int n, int mp ) {
+      double a0, ejn, f, f0, f1, hmp;
+      int it, n0, n1, nn;
+      double obj;
+      a0 = std::abs ( x );
+      hmp = 0.5 * mp;
+      ejn = envj ( n, a0 );
+      if ( ejn <= hmp ) {
+	obj = mp;
+	n0 = static_cast<int> ( 1.1 * a0 );
+      } else {
+	obj = hmp + ejn;
+	n0 = n;
+      }
+      f0 = envj ( n0, a0 ) - obj;
+      n1 = n0 + 5;
+      f1 = envj ( n1, a0 ) - obj;
+      for (int it = 1; it < 21; ++it) {
+	nn = n1 - ( n1 - n0 ) / ( 1.0 - f0 / f1 );
+	f = envj ( nn, a0 ) - obj;
+	if ( std::abs ( nn - n1 ) < 1 ) break;
+	n0 = n1;
+	f0 = f1;
+	n1 = nn;
+	f1 = f;
+      }
+      return  nn + 10;
+    }
+
+
+
+// !*****************************************************************************80
+// !
+// !! MSTA1 determines a backward recurrence starting point for Jn(x).
+// !
+// !  Discussion:
+// !
+// !    This procedure determines the starting point for backward  
+// !    recurrence such that the magnitude of    
+// !    Jn(x) at that point is about 10^(-MP).
+// !
+// !  Licensing:
+// !
+// !    This routine is copyrighted by Shanjie Zhang and Jianming Jin.  However, 
+// !    they give permission to incorporate this routine into a user program 
+// !    provided that the copyright is acknowledged.
+// !
+// !  Modified:
+// !
+// !    08 July 2012
+// !
+// !  Author:
+// !
+// !    Shanjie Zhang, Jianming Jin
+// !
+// !  Reference:
+// !
+// !    Shanjie Zhang, Jianming Jin,
+// !    Computation of Special Functions,
+// !    Wiley, 1996,
+// !    ISBN: 0-471-11963-6,
+// !    LC: QA351.C45.
+// !
+// !  Parameters:
+// !
+// !    Input, real ( kind = 8 ) X, the argument.
+// !
+// !    Input, integer ( kind = 4 ) MP, the negative logarithm of the 
+// !    desired magnitude.
+// !
+// !    Output, integer ( kind = 4 ) MSTA1, the starting point.
+// !
+    int  msta1 ( double x, int mp ) {
+      double a0, f, f0, f1;
+      int it, n0, n1, nn;
+      a0 = std::abs ( x );
+      n0 = static_cast<int> (1.1 * a0 ) + 1;
+      f0 = envj ( n0, a0 ) - mp;
+      n1 = n0 + 5;
+      f1 = envj ( n1, a0 ) - mp;
+      for (int it = 1; it <= 20; ++it) { 
+	nn = n1 - ( n1 - n0 ) / ( 1.0 - f0 / f1 );
+	f = envj ( nn, a0 ) - mp;
+	if ( abs ( nn - n1 ) < 1 ) break;
+	n0 = n1;
+	f0 = f1;
+	n1 = nn;
+	f1 = f;
+      }
+      return nn;
+    }
+      // function envj ( n, x )
+
+      // !*****************************************************************************80
+      // !
+      // !! ENVJ is a utility function used by MSTA1 and MSTA2.
+      // !
+      // !  Licensing:
+      // !
+      // !    This routine is copyrighted by Shanjie Zhang and Jianming Jin.  However, 
+      // !    they give permission to incorporate this routine into a user program 
+      // !    provided that the copyright is acknowledged.
+      // !
+      // !  Modified:
+      // !
+      // !    14 March 2012
+      // !
+      // !  Author:
+      // !
+      // !    Shanjie Zhang, Jianming Jin
+      // !
+      // !  Reference:
+      // !
+      // !    Shanjie Zhang, Jianming Jin,
+      // !    Computation of Special Functions,
+      // !    Wiley, 1996,
+      // !    ISBN: 0-471-11963-6,
+      // !    LC: QA351.C45.
+      // !
+      // !  Parameters:
+      // !
+      // !    Input, integer ( kind = 4 ) N, ?
+      // !
+      // !    Input, real ( kind = 8 ) X, ?
+      // !
+      // !    Output, real ( kind = 8 ) ENVJ, ?
+      // !
+    double envj (int n, double x ) {
+      return  0.5 * std::log10(6.28 * n) - n * std::log10(1.36 * x / static_cast<double>(n) );
+    }
+  }  // end of namespace bessel
 }  // end of namespace nmie

bessel.h

@@ -27,6 +27,16 @@
 #include <complex>
 #include <vector>
 namespace nmie {
-  std::vector< std::complex<double> > bessel_j(int n, std::complex<double> rho);
-  std::vector< std::complex<double> > bessel_y(int n, std::complex<double> rho);
+  namespace bessel {
+    std::vector< std::complex<double> > bessel_j(int n, std::complex<double> rho);
+    std::vector< std::complex<double> > bh_bessel_j(int n, std::complex<double> rho);
+    void csphjy (int n, std::complex<double>z, int& nm,
+		 std::vector< std::complex<double> >& csj,
+		 std::vector< std::complex<double> >& cdj,
+		 std::vector< std::complex<double> >& csy,
+		 std::vector< std::complex<double> >& cdy );
+    double envj (int n, double x );
+    int msta2 ( double x, int n, int mp );
+    int  msta1 ( double x, int mp );
+  }  // end of namespace bessel
 }  // end of namespace nmie

tests/bessel/test_bessel.cc

@@ -0,0 +1,90 @@
+//**********************************************************************************//
+//    Copyright (C) 2009-2015  Ovidio Pena <ovidio@bytesfall.com>                   //
+//    Copyright (C) 2013-2015  Konstantin Ladutenko <kostyfisik@gmail.com>          //
+//                                                                                  //
+//    This file is part of scattnlay                                                //
+//                                                                                  //
+//    This program is free software: you can redistribute it and/or modify          //
+//    it under the terms of the GNU General Public License as published by          //
+//    the Free Software Foundation, either version 3 of the License, or             //
+//    (at your option) any later version.                                           //
+//                                                                                  //
+//    This program is distributed in the hope that it will be useful,               //
+//    but WITHOUT ANY WARRANTY; without even the implied warranty of                //
+//    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the                 //
+//    GNU General Public License for more details.                                  //
+//                                                                                  //
+//    The only additional remark is that we expect that all publications            //
+//    describing work using this software, or all commercial products               //
+//    using it, cite the following reference:                                       //
+//    [1] O. Pena and U. Pal, "Scattering of electromagnetic radiation by           //
+//        a multilayered sphere," Computer Physics Communications,                  //
+//        vol. 180, Nov. 2009, pp. 2348-2354.                                       //
+//                                                                                  //
+//    You should have received a copy of the GNU General Public License             //
+//    along with this program.  If not, see <http://www.gnu.org/licenses/>.         //
+//**********************************************************************************//
+#include <complex>
+#include <cmath>
+#include <cstdio>
+#include <stdexcept>
+#include <vector>
+
+#include "../../bessel.h"
+
+int main(int argc, char *argv[]) {
+  
+  // auto result = nmie::bessel::bh_bessel_j(20, 1);
+  auto result = nmie::bessel::bessel_j(20, 1);
+  printf("Calculate and compare against Wolfram Alpha\n");
+  printf("j(0,1) = %16.18f\n", real(result[0]));
+  printf("WA j() = 0.841470984807896506652502321630\n");
+  printf("j(1,1) = %16.18f\n", real(result[1]));
+  printf("WA j() = 0.301168678939756789251565714187322395890252640\n");
+  printf("j(20,1) = %16.18f\n", real(result[20]));
+  printf("WA j() = 7.537795722236872993957557741584960855614358030 × 10^-26\n");
+  
+  result = nmie::bessel::bessel_j(4, std::complex<double>(1,2));
+  printf("Calculate and compare against Wolfram Alpha\n");
+  printf("j(0,1+i2) = re(%16.18f)\n         im(%16.18f)\n", 
+	 real(result[0]),
+	 imag(result[0]));
+  printf("WA j() = re(1.4169961192118759)\n         im(-0.874391197002)\n");
+
+  printf("j(1,1+i2) = re(%16.18f)\n         im(%16.18f)\n", 
+	 real(result[1]),
+	 imag(result[1]));
+  printf("WA j() = re(0.74785726329830368)\n         im(0.68179207555304)\n");
+
+  printf("j(4,1+i2) = re(%16.18f)\n         im(%16.18f)\n", 
+	 real(result[4]),
+	 imag(result[4]));
+  printf("WA j() = re(-0.01352410550046)\n         im(-0.027169663050653)\n");
+
+
+  int n = 10;
+  std::complex<double> z(1.0, 2.0);
+  int nm;
+  std::vector< std::complex<double> > csj, cdj, csy, cdy;
+  nmie::bessel::csphjy (n, z, nm, csj, cdj,  csy, cdy );
+
+  result =  csj;
+  printf("===========Calculate and compare against Wolfram Alpha\n");
+  printf("j(0,1+i2) = re(%16.18f)\n         im(%16.18f)\n", 
+	 real(result[0]),
+	 imag(result[0]));
+  printf("WA j() = re(1.4169961192118759)\n         im(-0.874391197002)\n");
+
+  printf("j(1,1+i2) = re(%16.18f)\n         im(%16.18f)\n", 
+	 real(result[1]),
+	 imag(result[1]));
+  printf("WA j() = re(0.74785726329830368)\n         im(0.68179207555304)\n");
+
+  printf("j(4,1+i2) = re(%16.18f)\n         im(%16.18f)\n", 
+	 real(result[4]),
+	 imag(result[4]));
+  printf("WA j() = re(-0.01352410550046)\n         im(-0.027169663050653)\n");
+
+
+}
+