Prelimany version of spherical Bessel functions.

CMakeLists.txt

@@ -1,7 +1,7 @@
 # --------------------------------------------------------------------------- #
 # Author:       Joey Dumont                   <joey.dumont@gmail.com>         #
 # Author:       Denis Gagnon                  <gagnon88@gmail.com>            #
-# Date:         2015-02-26                                                    #
+# Date:         2015-06-01                                                    #
 # Description:  CMake compilation instructions.                               #
 # ----------------------------------------------------------------------------#
 
@@ -31,7 +31,7 @@ set (CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -O3 -Wall -march=native -std=c++11")
 # --                  Compilation Instructions                   -- #
 # ----------------------------------------------------------------- #
 # -- Included files
-include_directories (${CMAKE_CURRENT_SOURCE_DIR}/include INC_LIST)
+include_directories (${CMAKE_CURRENT_SOURCE_DIR}/include)
 
 # -- Source files
 aux_source_directory (${CMAKE_CURRENT_SOURCE_DIR}/src SRC_LIST)

TODO

@@ -1,4 +1,7 @@
-July 2014
- [] Add support for real orders by carefully coding the reflection formulae.
- [] Add support for real arguments via a static cast
+June 2015
+ [✓] Add support for real orders by carefully coding the reflection formulae.
+ [ ] Add support for real arguments via a static cast
 	from double to complex<double>.
+ [ ] Add support for spherical Bessel functions.
+ [ ] Test for more values of order and argument.
+ [ ] Add support for derivatives of spherical Bessel functions.

include/complex_bessel.h

@@ -1,5 +1,5 @@
-#ifndef CBESSEL_
-#define CBESSEL_
+#ifndef COMPLEX_BESSEL
+#define COMPLEX_BESSEL
 
 /*! \file complex_bessel.h
  *
@@ -23,5 +23,6 @@
 #include "complex_bessel_bits/fortranLinkage.h"
 #include "complex_bessel_bits/utilities.h"
 #include "complex_bessel_bits/besselFunctions.h"
+#include "complex_bessel_bits/sph_besselFunctions.h"
 
-#endif
+#endif // COMPLEX_BESSEL

include/complex_bessel_bits/besselFunctions.h

@@ -24,8 +24,11 @@ namespace sp_bessel {
 
  /*! @name Evaluation of Bessel functions.
  * We implement Amos' Fortran subroutines in C++.
- * \todo Provide error detection and signaling.  */
+ * \todo Provide error detection and signaling.  
+ */
+
 ///@{
+
 /*! Using a function as a template parameter, we define a function that
  * computes the derivative of the Bessel functions \f$J,\,Y,\,H^{(1,2)},\,I,\,K\f$
  * using the recurrence relations \cite ABR65 (Sects. 9.1.27/9.6.26). */
@@ -34,7 +37,7 @@ inline std::complex<double> diffBessel(double order, std::complex<double> z, int
 {
     // For J, Y, H1 and H2, phase = -1. 
     // For I, e^(order*pi*i)K, phase = 1. 
-    // First term of the serie. 
+    // First term of the series. 
     double p = 1.0;
     std::complex<double> s = T(order-n, z);
 
@@ -225,17 +228,17 @@ inline std::complex<double> besselK(double order, std::complex<double> z)
   zbesk_wrap(zr,zi,nu,kode,N,&cyr,&cyi,&nz,&ierr); // Call Fortran subroutine.
 
   // In passing from C++ to FORTRAN, the exact zero becomes the numerical zero (10^(-14)). 
-    // The limiting form of K_nu(z) for high order, Gamma(nu)/2*(z/2)^(-nu),
-    // leads to product of the form zero*infinity for the imaginary part, which destroys numerical precision. We hence
-    // manually set the imaginary part of the answer to zero is the imaginary part of the input
-    // is zero.
-    if (zi == 0.0 && zr >= 0.0) cyi = 0.0;
-    std::complex<double> answer(cyr,cyi); 
+  // The limiting form of K_nu(z) for high order, Gamma(nu)/2*(z/2)^(-nu),
+  // leads to product of the form zero*infinity for the imaginary part, which destroys numerical precision. We hence
+  // manually set the imaginary part of the answer to zero is the imaginary part of the input
+  // is zero.
+  if (zi == 0.0 && zr >= 0.0) cyi = 0.0;
+  std::complex<double> answer(cyr,cyi); 
 
-    // In case of error, we print the error code.
-    if (ierr!=0) std::cout << "besselK: Error code " << ierr << "." << std::endl;
+  // In case of error, we print the error code.
+  if (ierr!=0) std::cout << "besselK: Error code " << ierr << "." << std::endl;
 
-    return answer;
+  return answer;
 }
 
 inline std::complex<double> expBesselK(double order, std::complex<double> z)
@@ -369,6 +372,6 @@ inline std::complex<double> biryp(std::complex<double> z)
 
 ///@}
 
-}
+} // namespace sp_bessel
 
 #endif // BESSELFUNCTIONS_H

include/complex_bessel_bits/fortranLinkage.h

@@ -62,6 +62,6 @@ extern void zbiry_wrap(double,double,int,int,double*,double*,int*,int*);
 }
 ///@}
 
-}
+} // namespace sp_bessel
 
 #endif // FORTRANLINKAGE_H

include/complex_bessel_bits/sph_besselFunctions.h

@@ -0,0 +1,60 @@
+/*! \file sph_besselFunctions.h
+ *
+ * \author Joey Dumont <joey.dumont@gmail.com>
+ * \author Denis Gagnon <gagnon88@gmail.com>
+ *
+ * \brief Functions computing the spherical Bessel functions.
+ *
+ * We compute the spherical Bessel functions by evaluating the
+ * Bessel functions with half-integer order with D.E. Amos' FORTRAN
+ * implementation. 
+ *
+ * \copyright LGPL
+ */
+
+#ifndef SPH_BESSELFUNCTIONS_H
+#define SPH_BESSELFUNCTIONS_H
+
+#include <iostream>
+#include <complex>
+
+#include "besselFunctions.h"
+#include "utilities.h"
+
+namespace sp_bessel {
+
+/// Near \f$z\rightarrow0\f$, the floating point division necessary would
+/// destroy precision. We thus use an ascending series to compute sph_besselJ.
+/// See \cite ABR65 Sec. 10.1.2.
+inline std::complex<double> sph_besselJ_asc_series(double order, std::complex<double> z)
+{
+  return 0;
+}
+
+/// We compute the spherical Bessel function of the first kind.
+inline std::complex<double> sph_besselJ(double order, std::complex<double> z)
+{
+  return std::sqrt(constants::pi/(2.0*z))*besselJ(order+0.5, z);
+}
+
+/// We compute the spherical Bessel function of the second kind.
+inline std::complex<double> sph_besselY(double order, std::complex<double> z)
+{
+  return std::sqrt(constants::pi/(2.0*z))*besselY(order+0.5,z);
+}
+
+/// We compute the spherical Hankel function of the first kind.
+inline std::complex<double> sph_hankelH1(double order, std::complex<double> z)
+{
+  return std::sqrt(constants::pi/(2.0*z))*hankelH1(order+0.5,z);
+}
+
+/// We compute the spherical Hankel function of the second kind.
+inline std::complex<double> sph_hankelH2(double order, std::complex<double> z)
+{
+  return std::sqrt(constants::pi/(2.0*z))*hankelH2(order+0.5,z);
+}
+
+} // namespace sp_bessel
+
+#endif // SPH_BESSELFUNCTIONS_H

include/complex_bessel_bits/utilities.h

@@ -51,6 +51,6 @@ inline double sin_pi(double nu)
   return std::sin(constants::pi*nu);
 }
 
-}
+} // namespace sp_bessel
 
 #endif // UTILITIES_H

tests/simpleTest.cpp

@@ -1,11 +1,13 @@
 #include <complex_bessel.h>
+//#include <complex_bessel_bits/sph_besselFunctions.h>
 #include <iostream>
 
 using namespace sp_bessel;
 
 int main(int argc, char* argv[])
 {
 	std::complex<double> z(1.0,0.0);
+	std::complex<double> nearOrigin(1.0e-10,0);
 	std::cout << besselJ(2,z) << std::endl;
 	std::cout << besselJp(2,z) << std::endl;
 	std::cout << besselJp(2,z,2) << std::endl;
@@ -16,5 +18,9 @@ int main(int argc, char* argv[])
 	std::cout << besselKp(2,z,2) << std::endl;
 	std::cout << hankelH1(2,z) << std::endl;
 	std::cout << hankelH1p(2,z) << std::endl;
+	std::cout << sph_besselJ(-constants::pi,nearOrigin) << std::endl;
+	std::cout << sph_besselY(-constants::pi,nearOrigin) << std::endl;
+	std::cout << sph_hankelH1(-constants::pi,nearOrigin) << std::endl;
+	std::cout << sph_hankelH2(-constants::pi,nearOrigin) << std::endl;
 	return 0;
 }