add tests from Hong Du paper (#34)

CMakeLists.txt

@@ -98,7 +98,7 @@ endif ()
 
 #include_directories(src)
 add_subdirectory(src)
-
+add_subdirectory(examples)
 #
 # Copy all python scripts to the build directory.
 #
@@ -111,3 +111,6 @@ foreach (_script ${Python_SCRIPTS})
             COPYONLY
     )
 endforeach ()
+
+enable_testing()
+add_subdirectory(tests)

examples/CMakeLists.txt

@@ -0,0 +1 @@
+add_executable(bug1 bug1.cc)

examples/bug1.cc

@@ -0,0 +1,64 @@
+//**********************************************************************************//
+//    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/>.         //
+//**********************************************************************************//
+//   This program evaluates absorption of a triple layered nanoparticle
+#include <complex>
+#include <cstdio>
+#include <string>
+#include "../src/nmie-applied.hpp"
+#include "../src/nmie-applied-impl.hpp"
+
+int main(int argc, char *argv[]) {
+  try {
+    nmie::MultiLayerMieApplied<double> multi_layer_mie;
+    const std::complex<double> index_Si(4, 0.01);
+    double delta = 1e-5;
+    double core_r = 100; //nm Si
+    double host = 2.;
+    multi_layer_mie.AddTargetLayer(core_r*host, index_Si/host);
+    multi_layer_mie.SetWavelength(400-delta);
+    multi_layer_mie.RunMieCalculation();
+    double Qsca = multi_layer_mie.GetQsca();
+    printf("at WL = 400-(%g) the result is (Qsca - ref)=%15.14g\n", delta, Qsca-2.382076221);
+
+    multi_layer_mie.SetWavelength(400);
+    multi_layer_mie.RunMieCalculation();
+    Qsca = multi_layer_mie.GetQsca();
+    printf("at WL = 400 the result is (Qsca - ref)=%15.14g\n", Qsca-2.382076221);
+
+//    multi_layer_mie.SetWavelength(400+delta);
+//    multi_layer_mie.RunMieCalculation();
+//    Qsca = multi_layer_mie.GetQsca();
+//    printf("Qsca = %15.14g\n", Qsca-2.382076221);
+  } catch( const std::invalid_argument& ia ) {
+    // Will catch if  multi_layer_mie fails or other errors.
+    std::cerr << "Invalid argument: " << ia.what() << std::endl;
+    return -1;
+  }
+    return 0;
+}
+
+

src/CMakeLists.txt

@@ -5,7 +5,7 @@ set(_scattnlay_python_sources
         ${CMAKE_CURRENT_LIST_DIR}/nmie-pybind11.hpp
         ${CMAKE_CURRENT_LIST_DIR}/nmie-pybind11.cc
         ${CMAKE_CURRENT_LIST_DIR}/nmie-precision.hpp
-        ${CMAKE_CURRENT_LIST_DIR}/nmie-impl.cc
+        ${CMAKE_CURRENT_LIST_DIR}/nmie-impl.hpp
         ${CMAKE_CURRENT_LIST_DIR}/pb11_wrapper.cc)
 
 # Define python extension

src/nmie-applied.hpp

@@ -37,7 +37,7 @@
 #include <iostream>
 #include <vector>
 #include "nmie.hpp"
-#include "nmie-impl.cc"
+#include "nmie-impl.hpp"
 
 
 namespace nmie {
@@ -60,11 +60,11 @@ namespace nmie {
       do {
 	if (!output) break;
 	++iformat;
-	printf("%23.13e",var);	     
+	printf("%23.13e",var);
 	if (iformat%4 == 0) printf("\n");
       } while (false);
     }
-    // Set parameters in applied units 
+    // Set parameters in applied units
     void SetWavelength(FloatType wavelength) {
       this->MultiLayerMie<FloatType>::MarkUncalculated();
       wavelength_ = wavelength;};
@@ -93,7 +93,7 @@ namespace nmie {
         this->MultiLayerMie<FloatType>::SetModeNmaxAndType(mode_n, mode_type);};
     void SetAnglesForPattern(FloatType from_angle, FloatType to_angle, int samples);
     std::vector<FloatType> GetAngles();
-    
+
     void ClearTarget();
     void ClearCoating();
     void ClearLayers();
@@ -124,7 +124,7 @@ namespace nmie {
     // Size parameter units
     std::vector<FloatType> GetLayerWidthSP();
     // Same as to get target and coating index
-    std::vector< std::complex<FloatType> > GetLayerIndex();  
+    std::vector< std::complex<FloatType> > GetLayerIndex();
     std::vector< std::array<FloatType,3> > GetFieldPointsSP();
     // Do we need normalize field to size parameter?
     /* std::vector<std::vector<std::complex<FloatType> > >  GetFieldESP(); */
@@ -163,7 +163,7 @@ namespace nmie {
     void sphericalBessel(std::complex<FloatType> z, std::vector<std::complex<FloatType> >& bj,
 			             std::vector<std::complex<FloatType> >& by, std::vector<std::complex<FloatType> >& bd);
     std::complex<FloatType> calcD1confra(int N, const std::complex<FloatType> z);
-    
+
     FloatType wavelength_ = 1.0;
     FloatType total_radius_ = 0.0;
     /// Width and index for each layer of the structure

src/nmie-impl.cc → src/nmie-impl.hpp

@@ -350,7 +350,8 @@ namespace nmie {
   // Calculate an - equation (5)                                            //
   // ********************************************************************** //
   template <typename FloatType>
-  std::complex<FloatType> MultiLayerMie<FloatType>::calc_an(int n, FloatType XL, std::complex<FloatType> Ha, std::complex<FloatType> mL,
+  std::complex<FloatType> MultiLayerMie<FloatType>::
+      calc_an(int n, FloatType XL, std::complex<FloatType> Ha, std::complex<FloatType> mL,
                                               std::complex<FloatType> PsiXL, std::complex<FloatType> ZetaXL,
                                               std::complex<FloatType> PsiXLM1, std::complex<FloatType> ZetaXLM1) {
 

src/nmie.cc

@@ -51,7 +51,7 @@
 
 #include "nmie.hpp"
 #include "nmie-precision.hpp"
-#include "nmie-impl.cc"
+#include "nmie-impl.hpp"
 
 namespace nmie {
 

tests/CMakeLists.txt

@@ -0,0 +1,28 @@
+cmake_minimum_required(VERSION 3.15)
+project(scattnlay_tests C CXX)
+# -- Dependency (Google Test)
+find_package(GTest REQUIRED)
+include_directories(${GTEST_INCLUDE_DIRS})
+set(LIBS ${LIBS} ${GTEST_LIBRARIES})
+set(LIBS ${LIBS} pthread)
+
+# -- Output tests in directory
+add_executable("test_bulk_sphere" test_bulk_sphere.cc)
+target_link_libraries("test_bulk_sphere" ${LIBS})
+add_test(NAME "test_bulk_sphere"
+        COMMAND "test_bulk_sphere")
+
+add_executable("test_bulk_sphere_multi_precision" test_bulk_sphere_multi_precision.cc)
+target_compile_options("test_bulk_sphere_multi_precision"
+        PRIVATE -DMULTI_PRECISION=100)
+target_link_libraries("test_bulk_sphere_multi_precision" ${LIBS})
+add_test(NAME "test_bulk_sphere_multi_precision"
+        COMMAND "test_bulk_sphere_multi_precision")
+
+
+add_executable("test_Riccati_Bessel_logarithmic_derivative"
+        test_Riccati_Bessel_logarithmic_derivative.cc)
+target_link_libraries("test_Riccati_Bessel_logarithmic_derivative" ${LIBS})
+add_test(NAME "test_Riccati_Bessel_logarithmic_derivative"
+        COMMAND "test_Riccati_Bessel_logarithmic_derivative")
+

tests/test_Riccati_Bessel_logarithmic_derivative.cc

@@ -0,0 +1,10 @@
+#include "gtest/gtest.h"
+#include "../src/nmie-impl.hpp"
+TEST(RiccatiBesselTest, HandlesInput) {
+  nmie::MultiLayerMie<double> nmie;
+  EXPECT_EQ(1, 1)<<"Should be equal";
+}
+int main(int argc, char **argv) {
+  testing::InitGoogleTest(&argc, argv);
+  return RUN_ALL_TESTS();
+}

tests/test_bulk_sphere.cc

@@ -0,0 +1,43 @@
+#include "gtest/gtest.h"
+#include "../src/nmie-impl.hpp"
+TEST(BulkSphere, HandlesInput) {
+  nmie::MultiLayerMie<double> nmie;
+  // A list of tests for a bulk sphere from
+  // Hong Du, "Mie-scattering calculation," Appl. Opt. 43, 1951-1956 (2004)
+  // table 1: sphere size and refractive index
+  // followed by resulting extinction and scattering efficiencies
+  std::vector< std::tuple< double, std::complex<double>, double, double, char > >
+      parameters_and_results
+      {
+          // x, {Re(m), Im(m)}, Qext, Qsca, test_name
+          {0.099, {0.75,0}, 7.417859e-06, 7.417859e-06, 'a'},
+          {0.101, {0.75,0}, 8.033538e-06, 8.033538e-06, 'b'},
+          {10,    {0.75,0},     2.232265, 2.232265, 'c'},
+          {1000,  {0.75,0},     1.997908, 1.997908, 'd'},
+//          {100,   {1.33,-1e-5}, 2.101321, 2.096594, 'e'},
+//          {10000, {1.33,-1e-5}, 2.004089, 1.723857, 'f'},
+//          {0.055, {1.5, -1},    0.101491, 1.131687e-05, 'g'},
+//          {0.056, {1.5, -1},   0.1033467, 1.216311e-05, 'h'},
+//          {100,   {1.5, -1},    2.097502, 1.283697, 'i'},
+//          {10000, {1.5, -1},    2.004368, 1.236574, 'j'},
+//          {1,     {10,  -10},   2.532993, 2.049405, 'k'},
+//          {100,   {10,  -10,},  2.071124, 1.836785, 'l'},
+//          {10000, {10,  -10},   2.005914, 1.795393, 'm'},
+      };
+  for (const auto &data : parameters_and_results) {
+    nmie.SetLayersSize({std::get<0>(data)});
+    nmie.SetLayersIndex({std::get<1>(data)});
+    nmie.RunMieCalculation();
+    double Qext = nmie.GetQext();
+    double Qsca = nmie.GetQsca();
+    EXPECT_FLOAT_EQ(std::get<2>(data), Qext)
+              << "Extinction of the bulk sphere, test case:" << std::get<4>(data);
+    EXPECT_FLOAT_EQ(std::get<3>(data), Qsca)
+              << "Scattering of the bulk sphere, test case:" << std::get<4>(data);
+  }
+
+}
+int main(int argc, char **argv) {
+  testing::InitGoogleTest(&argc, argv);
+  return RUN_ALL_TESTS();
+}

tests/test_bulk_sphere_multi_precision.cc

@@ -0,0 +1,44 @@
+#include "gtest/gtest.h"
+#include "../src/nmie-impl.hpp"
+#include "../src/nmie-precision.hpp"
+TEST(BulkSphereMultiPrecision, HandlesInput) {
+  nmie::MultiLayerMie<nmie::FloatType> nmie;
+  // A list of tests for a bulk sphere from
+  // Hong Du, "Mie-scattering calculation," Appl. Opt. 43, 1951-1956 (2004)
+  // table 1: sphere size and refractive index
+  // followed by resulting extinction and scattering efficiencies
+  std::vector< std::tuple< double, std::complex<double>, double, double, char > >
+      parameters_and_results
+      {
+          // x, {Re(m), Im(m)}, Qext, Qsca, test_name
+          {0.099, {0.75,0}, 7.417859e-06, 7.417859e-06, 'a'},
+          {0.101, {0.75,0}, 8.033538e-06, 8.033538e-06, 'b'},
+          {10,    {0.75,0},     2.232265, 2.232265, 'c'},
+          {1000,  {0.75,0},     1.997908, 1.997908, 'd'},
+//          {100,   {1.33,-1e-5}, 2.101321, 2.096594, 'e'},
+//          {10000, {1.33,-1e-5}, 2.004089, 1.723857, 'f'},
+//          {0.055, {1.5, -1},    0.101491, 1.131687e-05, 'g'},
+//          {0.056, {1.5, -1},   0.1033467, 1.216311e-05, 'h'},
+//          {100,   {1.5, -1},    2.097502, 1.283697, 'i'},
+//          {10000, {1.5, -1},    2.004368, 1.236574, 'j'},
+//          {1,     {10,  -10},   2.532993, 2.049405, 'k'},
+//          {100,   {10,  -10,},  2.071124, 1.836785, 'l'},
+//          {10000, {10,  -10},   2.005914, 1.795393, 'm'},
+      };
+  for (const auto &data : parameters_and_results) {
+    nmie.SetLayersSize({std::get<0>(data)});
+    nmie.SetLayersIndex({std::get<1>(data)});
+    nmie.RunMieCalculation();
+    double Qext = static_cast<double>(nmie.GetQext());
+    double Qsca = static_cast<double>(nmie.GetQsca());
+    EXPECT_FLOAT_EQ(std::get<2>(data), Qext)
+              << "Extinction of the bulk sphere, test case:" << std::get<4>(data);
+    EXPECT_FLOAT_EQ(std::get<3>(data), Qsca)
+              << "Scattering of the bulk sphere, test case:" << std::get<4>(data);
+  }
+
+}
+int main(int argc, char **argv) {
+  testing::InitGoogleTest(&argc, argv);
+  return RUN_ALL_TESTS();
+}