initial conversion to templates of nmie.h

src/nmie.cc

@@ -45,8 +45,14 @@
 #include <cstdlib>
 #include <stdexcept>
 #include <vector>
+#include <boost/multiprecision/cpp_bin_float.hpp>
+//#include <boost/math/special_functions/gamma.hpp>
 
 namespace nmie {
+  //typedef float FloatType;
+  //typedef double FloatType;
+  typedef boost::multiprecision::cpp_bin_float_100 FloatType;
+  
   //**********************************************************************************//
   // This function emulates a C call to calculate the scattering coefficients         //
   // required to calculate both the near- and far-field parameters.                   //
@@ -71,16 +77,16 @@ namespace nmie {
     if (x.size() != L || m.size() != L)
         throw std::invalid_argument("Declared number of layers do not fit x and m!");
     try {
-      MultiLayerMie<double> ml_mie;
-      ml_mie.SetLayersSize(x);
-      ml_mie.SetLayersIndex(m);
+      MultiLayerMie<FloatType> ml_mie;
+      ml_mie.SetLayersSize(ConvertVector<FloatType>(x));
+      ml_mie.SetLayersIndex(ConvertComplexVector<FloatType>(m));
       ml_mie.SetPECLayer(pl);
       ml_mie.SetMaxTerms(nmax);
 
       ml_mie.calcScattCoeffs();
 
-      an = ml_mie.GetAn();
-      bn = ml_mie.GetBn();
+      an = ConvertComplexVector<double>(ml_mie.GetAn());
+      bn = ConvertComplexVector<double>(ml_mie.GetBn());
 
       return ml_mie.GetMaxTerms();
     } catch(const std::invalid_argument& ia) {
@@ -128,25 +134,24 @@ namespace nmie {
     if (Theta.size() != nTheta)
         throw std::invalid_argument("Declared number of sample for Theta is not correct!");
     try {
-      typedef double FloatType;
-      MultiLayerMie<double> ml_mie;
-      ml_mie.SetLayersSize(x);
-      ml_mie.SetLayersIndex(m);
-      ml_mie.SetAngles(Theta);
+      MultiLayerMie<FloatType> ml_mie;
+      ml_mie.SetLayersSize(ConvertVector<FloatType>(x));
+      ml_mie.SetLayersIndex(ConvertComplexVector<FloatType>(m));
+      ml_mie.SetAngles(ConvertVector<FloatType>(Theta));
       ml_mie.SetPECLayer(pl);
       ml_mie.SetMaxTerms(nmax);
 
       ml_mie.RunMieCalculation();
 
-      *Qext = ml_mie.GetQext();
-      *Qsca = ml_mie.GetQsca();
-      *Qabs = ml_mie.GetQabs();
-      *Qbk = ml_mie.GetQbk();
-      *Qpr = ml_mie.GetQpr();
-      *g = ml_mie.GetAsymmetryFactor();
-      *Albedo = ml_mie.GetAlbedo();
-      S1 = ml_mie.GetS1();
-      S2 = ml_mie.GetS2();
+      *Qext = static_cast<double>(ml_mie.GetQext());
+      *Qsca = static_cast<double>(ml_mie.GetQsca());
+      *Qabs = static_cast<double>(ml_mie.GetQabs());
+      *Qbk = static_cast<double>(ml_mie.GetQbk());
+      *Qpr = static_cast<double>(ml_mie.GetQpr());
+      *g = static_cast<double>(ml_mie.GetAsymmetryFactor());
+      *Albedo = static_cast<double>(ml_mie.GetAlbedo());
+      S1 = ConvertComplexVector<double>(ml_mie.GetS1());
+      S2 = ConvertComplexVector<double>(ml_mie.GetS2());
 
       return ml_mie.GetMaxTerms();
     } catch(const std::invalid_argument& ia) {
@@ -292,14 +297,16 @@ namespace nmie {
       if (f.size() != 3)
         throw std::invalid_argument("Field H is not 3D!");
     try {
-      MultiLayerMie<double> ml_mie;
+      MultiLayerMie<FloatType> ml_mie;
+      ml_mie.SetLayersSize(ConvertVector<FloatType>(x));
+      ml_mie.SetLayersIndex(ConvertComplexVector<FloatType>(m));
       ml_mie.SetPECLayer(pl);
-      ml_mie.SetLayersSize(x);
-      ml_mie.SetLayersIndex(m);
-      ml_mie.SetFieldCoords({Xp_vec, Yp_vec, Zp_vec});
+      ml_mie.SetFieldCoords({ConvertVector<FloatType>(Xp_vec),
+	    ConvertVector<FloatType>(Yp_vec),
+	    ConvertVector<FloatType>(Zp_vec) });
       ml_mie.RunFieldCalculation();
-      E = ml_mie.GetFieldE();
-      H = ml_mie.GetFieldH();
+      E = ConvertComplexVectorVector<double>(ml_mie.GetFieldE());
+      H = ConvertComplexVectorVector<double>(ml_mie.GetFieldH());
 
       return ml_mie.GetMaxTerms();
     } catch(const std::invalid_argument& ia) {

src/nmie.hpp

@@ -33,14 +33,8 @@
 #include <cstdlib>
 #include <iostream>
 #include <vector>
-
+#include <boost/math/constants/constants.hpp>
 namespace nmie {
-  //Used constants TODO! Change to boost PI
-  const double PI_=3.14159265358979323846;
-  // light speed [m s-1]
-  const double cc_ = 2.99792458e8;
-  // assume non-magnetic (MU=MU0=const) [N A-2]
-  const double mu_ = 4.0*PI_*1.0e-7;
   int ScattCoeffs(const unsigned int L, const int pl, std::vector<double>& x, std::vector<std::complex<double> >& m, const int nmax, std::vector<std::complex<double> >& an, std::vector<std::complex<double> >& bn);
   int nMie(const unsigned int L, const int pl, std::vector<double>& x, std::vector<std::complex<double> >& m, const unsigned int nTheta, std::vector<double>& Theta, const int nmax, double *Qext, double *Qsca, double *Qabs, double *Qbk, double *Qpr, double *g, double *Albedo, std::vector<std::complex<double> >& S1, std::vector<std::complex<double> >& S2);
   int nMie(const unsigned int L, std::vector<double>& x, std::vector<std::complex<double> >& m, const unsigned int nTheta, std::vector<double>& Theta, double *Qext, double *Qsca, double *Qabs, double *Qbk, double *Qpr, double *g, double *Albedo, std::vector<std::complex<double> >& S1, std::vector<std::complex<double> >& S2);
@@ -49,8 +43,14 @@ namespace nmie {
   int nField(const unsigned int L, const int pl, const std::vector<double>& x, const std::vector<std::complex<double> >& m, const int nmax, const unsigned int ncoord, const std::vector<double>& Xp, const std::vector<double>& Yp, const std::vector<double>& Zp, std::vector<std::vector<std::complex<double> > >& E, std::vector<std::vector<std::complex<double> > >& H);
 
   template <typename FloatType = double>
-  class MultiLayerMie {
+  class MultiLayerMie {    
    public:
+    //Used constants TODO! Change to boost PI
+    const double PI_=3.14159265358979323846;
+    // light speed [m s-1]
+    const double cc_ = 2.99792458e8;
+    // assume non-magnetic (MU=MU0=const) [N A-2]
+    const double mu_ = 4.0*PI_*1.0e-7;
     // Run calculation
     void RunMieCalculation();
     void RunFieldCalculation();

tests/c++/speed-test.cc

@@ -39,7 +39,7 @@
 //sudo aptitude install libgoogle-perftools-dev
 //#include <google/heap-profiler.h>
 #include "../../src/nmie.hpp"
-#include "../../src/nmie-impl.hpp"
+//#include "../../src/nmie-impl.hpp"
 
 timespec diff(timespec start, timespec end);
 const double PI=3.14159265358979323846;