|
|
@@ -0,0 +1,61 @@
|
|
|
+//**********************************************************************************//
|
|
|
+// 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 <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];
|
|
|
+ }
|
|
|
+ }
|
|
|
+ // 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];
|
|
|
+ }
|
|
|
+ }
|
|
|
+} // end of namespace nmie
|
