nmie-applied compiles with double

Makefile

@@ -2,7 +2,7 @@ PYTHON=`which python`
 CYTHON=`which cython`
 DESTDIR=/
 PROJECT=python-scattnlay
-VERSION=2.0
+VERSION=2.1
 BUILDIR=$(CURDIR)/debian/$(PROJECT)
 SRCDIR=$(CURDIR)/src
 MULTIPREC=100

src/nmie-applied.cc

@@ -1,5 +1,5 @@
 ///
-/// @file   nmie.cc
+/// @file   nmie-applied.cc
 /// @author Ladutenko Konstantin <kostyfisik at gmail (.) com>
 /// @date   Tue Sep  3 00:38:27 2013
 /// @copyright 2013,2014,2015 Ladutenko Konstantin
@@ -30,6 +30,8 @@
 /// @brief  Wrapper class around nMie function for ease of use
 ///
 #include "nmie-applied.hpp"
+#include "nmie-applied-impl.hpp"
+#include "nmie-precision.hpp"
 #include <array>
 #include <algorithm>
 #include <cstdio>
@@ -114,350 +116,5 @@ namespace nmie {
   // ********************************************************************** //
   // ********************************************************************** //
   // ********************************************************************** //
-  void MultiLayerMieApplied::GetFailed() {
-    double faild_x = 9.42477796076938;
-    //double faild_x = 9.42477796076937;
-    std::complex<double> z(faild_x, 0.0);
-    std::vector<int> nmax_local_array = {20, 100, 500, 2500};
-    for (auto nmax_local : nmax_local_array) {
-      std::vector<std::complex<double> > D1_failed(nmax_local + 1);
-      // Downward recurrence for D1 - equations (16a) and (16b)
-      D1_failed[nmax_local] = std::complex<double>(0.0, 0.0);
-      const std::complex<double> zinv = std::complex<double>(1.0, 0.0)/z;
-      for (int n = nmax_local; n > 0; n--) {
-        D1_failed[n - 1] = double(n)*zinv - 1.0/(D1_failed[n] + double(n)*zinv);
-      }
-      printf("Faild D1[0] from reccurence (z = %16.14f, nmax = %d): %g\n",
-             faild_x, nmax_local, D1_failed[0].real());
-    }
-    printf("Faild D1[0] from continued fraction (z = %16.14f): %g\n", faild_x,
-           calcD1confra(0,z).real());
-    //D1[nmax_] = calcD1confra(nmax_, z);
-  
-    
-  }
-  // ********************************************************************** //
-  // ********************************************************************** //
-  // ********************************************************************** //
-  void MultiLayerMieApplied::AddTargetLayer(double width, std::complex<double> layer_index) {
-    MarkUncalculated();
-    if (width <= 0)
-      throw std::invalid_argument("Layer width should be positive!");
-    target_width_.push_back(width);
-    target_index_.push_back(layer_index);
-  }  // end of void  MultiLayerMieApplied::AddTargetLayer(...)  
-  // ********************************************************************** //
-  // ********************************************************************** //
-  // ********************************************************************** //
-  void MultiLayerMieApplied::SetTargetPEC(double radius) {
-    MarkUncalculated();
-    if (target_width_.size() != 0 || target_index_.size() != 0)
-      throw std::invalid_argument("Error! Define PEC target radius before any other layers!");
-    // Add layer of any index...
-    AddTargetLayer(radius, std::complex<double>(0.0, 0.0));
-    // ... and mark it as PEC
-    SetPECLayer(0);
-  }
-  // ********************************************************************** //
-  // ********************************************************************** //
-  // ********************************************************************** //
-  void MultiLayerMieApplied::SetCoatingIndex(std::vector<std::complex<double> > index) {
-    MarkUncalculated();
-    coating_index_.clear();
-    for (auto value : index) coating_index_.push_back(value);
-  }  // end of void MultiLayerMieApplied::SetCoatingIndex(std::vector<complex> index);  
-  // ********************************************************************** //
-  // ********************************************************************** //
-  // ********************************************************************** //
-  void MultiLayerMieApplied::SetCoatingWidth(std::vector<double> width) {
-    MarkUncalculated();
-    coating_width_.clear();
-    for (auto w : width)
-      if (w <= 0)
-        throw std::invalid_argument("Coating width should be positive!");
-      else coating_width_.push_back(w);
-  }
-  // end of void MultiLayerMieApplied::SetCoatingWidth(...);
-  // ********************************************************************** //
-  // ********************************************************************** //
-  // ********************************************************************** //
-  void MultiLayerMieApplied::SetWidthSP(const std::vector<double>& size_parameter) {
-    MarkUncalculated();
-    size_param_.clear();
-    double prev_size_parameter = 0.0;
-    for (auto layer_size_parameter : size_parameter) {
-      if (layer_size_parameter <= 0.0)
-        throw std::invalid_argument("Size parameter should be positive!");
-      if (prev_size_parameter > layer_size_parameter) 
-        throw std::invalid_argument
-          ("Size parameter for next layer should be larger than the previous one!");
-      prev_size_parameter = layer_size_parameter;
-      size_param_.push_back(layer_size_parameter);
-    }
-  }
-  // end of void MultiLayerMieApplied::SetWidthSP(...);
-  // ********************************************************************** //
-  // ********************************************************************** //
-  // ********************************************************************** //
-  void MultiLayerMieApplied::SetIndexSP(const std::vector< std::complex<double> >& index) {
-    MarkUncalculated();
-    //refractive_index_.clear();
-    refractive_index_ = index;
-    // for (auto value : index) refractive_index_.push_back(value);
-  }  // end of void MultiLayerMieApplied::SetIndexSP(...);  
-  // ********************************************************************** //
-  // ********************************************************************** //
-  // ********************************************************************** //
-  void MultiLayerMieApplied::SetFieldPointsSP(const std::vector< std::vector<double> >& coords_sp) {
-    if (coords_sp.size() != 3)
-      throw std::invalid_argument("Error! Wrong dimension of field monitor points!");
-    if (coords_sp[0].size() != coords_sp[1].size() || coords_sp[0].size() != coords_sp[2].size())
-      throw std::invalid_argument("Error! Missing coordinates for field monitor points!");
-    coords_sp_ = coords_sp;
-    // for (int i = 0; i < coords_sp_[0].size(); ++i) {
-    //   printf("%g, %g, %g\n", coords_sp_[0][i], coords_sp_[1][i], coords_sp_[2][i]);
-    // }
-  }  // end of void MultiLayerMieApplied::SetFieldPointsSP(...)
-  // ********************************************************************** //
-  // ********************************************************************** //
-  // ********************************************************************** //
-  void MultiLayerMieApplied::GenerateSizeParameter() {
-    MarkUncalculated();
-    size_param_.clear();
-    double radius = 0.0;
-    for (auto width : target_width_) {
-      radius += width;
-      size_param_.push_back(2*PI_*radius/wavelength_);
-    }
-    for (auto width : coating_width_) {
-      radius += width;
-      size_param_.push_back(2*PI_*radius/wavelength_);
-    }
-    total_radius_ = radius;
-  }  // end of void MultiLayerMieApplied::GenerateSizeParameter();
-  // ********************************************************************** //
-  // ********************************************************************** //
-  // ********************************************************************** //
-  void MultiLayerMieApplied::GenerateIndex() {
-    MarkUncalculated();
-    refractive_index_.clear();
-    for (auto index : target_index_) refractive_index_.push_back(index);
-    for (auto index : coating_index_) refractive_index_.push_back(index);
-  }  // end of void MultiLayerMieApplied::GenerateIndex();
-  // ********************************************************************** //
-  // ********************************************************************** //
-  // ********************************************************************** //
-  double MultiLayerMieApplied::GetTotalRadius() {
-    if (!isMieCalculated())  GenerateSizeParameter();
-    return total_radius_;      
-  }  // end of double MultiLayerMieApplied::GetTotalRadius();
-  // ********************************************************************** //
-  // ********************************************************************** //
-  // ********************************************************************** //
-  std::vector< std::vector<double> >
-  MultiLayerMieApplied::GetSpectra(double from_WL, double to_WL, int samples) {
-    if (!isMieCalculated())
-      throw std::invalid_argument("You should run calculations before result request!");
-    std::vector< std::vector<double> > spectra;
-    double step_WL = (to_WL - from_WL)/static_cast<double>(samples);
-    double wavelength_backup = wavelength_;
-    long fails = 0;
-    for (double WL = from_WL; WL < to_WL; WL += step_WL) {
-      wavelength_ = WL;
-      try {
-        RunMieCalculation();
-      } catch(const std::invalid_argument& ia) {
-        fails++;
-        continue;
-      }
-      //printf("%3.1f ",WL);
-      spectra.push_back(std::vector<double>({wavelength_, GetQext(),
-	      GetQsca(), GetQabs(), GetQbk()}));
-    }  // end of for each WL in spectra
-    printf("Spectrum has %li fails\n",fails);
-    wavelength_ = wavelength_backup;
-    return spectra;
-  }
-  // ********************************************************************** //
-  // ********************************************************************** //
-  // ********************************************************************** //
-  void MultiLayerMieApplied::ClearTarget() {
-    MarkUncalculated();
-    target_width_.clear();
-    target_index_.clear();
-  }
-  // ********************************************************************** //
-  // ********************************************************************** //
-  // ********************************************************************** //
-  void MultiLayerMieApplied::ClearCoating() {
-    MarkUncalculated();
-    coating_width_.clear();
-    coating_index_.clear();
-  }
-  // ********************************************************************** //
-  // ********************************************************************** //
-  // ********************************************************************** //
-  void MultiLayerMieApplied::ClearLayers() {
-    MarkUncalculated();
-    ClearTarget();
-    ClearCoating();
-  }
-  // ********************************************************************** //
-  // ********************************************************************** //
-  // ********************************************************************** //
-  void MultiLayerMieApplied::ClearAllDesign() {
-    MarkUncalculated();
-    ClearLayers();
-    size_param_.clear();
-    refractive_index_.clear();
-  }
-  // ********************************************************************** //
-  // ********************************************************************** //
-  // ********************************************************************** //
-  //                         Computational core
-  // ********************************************************************** //
-  // ********************************************************************** //
-  // ********************************************************************** //
-  //**********************************************************************************//
-  // Function CONFRA ported from MIEV0.f (Wiscombe,1979)
-  // Ref. to NCAR Technical Notes, Wiscombe, 1979
-  /*
-c         Compute Bessel function ratio A-sub-N from its
-c         continued fraction using Lentz method
-
-c         ZINV = Reciprocal of argument of A
-
-
-c    I N T E R N A L    V A R I A B L E S
-c    ------------------------------------
-
-c    CAK      Term in continued fraction expansion of A (Eq. R25)
-c     a_k
-
-c    CAPT     Factor used in Lentz iteration for A (Eq. R27)
-c     T_k
-
-c    CNUMER   Numerator   in capT  (Eq. R28A)
-c     N_k
-c    CDENOM   Denominator in capT  (Eq. R28B)
-c     D_k
-
-c    CDTD     Product of two successive denominators of capT factors
-c                 (Eq. R34C)
-c     xi_1
-
-c    CNTN     Product of two successive numerators of capT factors
-c                 (Eq. R34B)
-c     xi_2
-
-c    EPS1     Ill-conditioning criterion
-c    EPS2     Convergence criterion
-
-c    KK       Subscript k of cAk  (Eq. R25B)
-c     k
-
-c    KOUNT    Iteration counter (used to prevent infinite looping)
-
-c    MAXIT    Max. allowed no. of iterations
-
-c    MM + 1  and - 1, alternately
-*/
-  std::complex<double> MultiLayerMieApplied::calcD1confra(const int N, const std::complex<double> z) {
-  // NTMR -> nmax_ - 1  \\TODO nmax_ ?
-    //int N = nmax_ - 1;
-    int KK, KOUNT, MAXIT = 10000, MM;
-    //    double EPS1=1.0e-2;
-    double EPS2=1.0e-8;
-    std::complex<double> CAK, CAPT, CDENOM, CDTD, CNTN, CNUMER;
-    std::complex<double> one = std::complex<double>(1.0,0.0);
-    std::complex<double> ZINV = one/z;
-// c                                 ** Eq. R25a
-    std::complex<double> CONFRA = static_cast<std::complex<double> >(N + 1)*ZINV;   //debug ZINV
-    MM = - 1; 
-    KK = 2*N +3; //debug 3
-// c                                 ** Eq. R25b, k=2
-    CAK    = static_cast<std::complex<double> >(MM*KK)*ZINV; //debug -3 ZINV
-    CDENOM = CAK;
-    CNUMER = CDENOM + one/CONFRA; //-3zinv+z
-    KOUNT  = 1;
-    //10 CONTINUE
-    do {      ++KOUNT;
-      if (KOUNT > MAXIT) {
-        printf("re(%g):im(%g)\t\n", CONFRA.real(), CONFRA.imag());
-        throw std::invalid_argument("ConFra--Iteration failed to converge!\n");
-      }
-      MM *= - 1;      KK += 2;  //debug  mm=1 kk=5
-      CAK = static_cast<std::complex<double> >(MM*KK)*ZINV; //    ** Eq. R25b //debug 5zinv
-     //  //c ** Eq. R32    Ill-conditioned case -- stride two terms instead of one
-     //  if (std::abs(CNUMER/CAK) >= EPS1 ||  std::abs(CDENOM/CAK) >= EPS1) {
-     //         //c                       ** Eq. R34
-     //         CNTN   = CAK*CNUMER + 1.0;
-     //         CDTD   = CAK*CDENOM + 1.0;
-     //         CONFRA = (CNTN/CDTD)*CONFRA; // ** Eq. R33
-     //         MM  *= - 1;        KK  += 2;
-     //         CAK = static_cast<std::complex<double> >(MM*KK)*ZINV; // ** Eq. R25b
-     //         //c                        ** Eq. R35
-     //         CNUMER = CAK + CNUMER/CNTN;
-     //         CDENOM = CAK + CDENOM/CDTD;
-     //         ++KOUNT;
-     //         //GO TO  10
-     //         continue;
-     // } else { //c                           *** Well-conditioned case
-      {
-        CAPT   = CNUMER/CDENOM; // ** Eq. R27 //debug (-3zinv + z)/(-3zinv)
-        // printf("re(%g):im(%g)**\t", CAPT.real(), CAPT.imag());
-       CONFRA = CAPT*CONFRA; // ** Eq. R26
-       //if (N == 0) {output=true;printf(" re:");prn(CONFRA.real());printf(" im:"); prn(CONFRA.imag());output=false;};
-       //c                                  ** Check for convergence; Eq. R31
-       if (std::abs(CAPT.real() - 1.0) >= EPS2 ||  std::abs(CAPT.imag()) >= EPS2) {
-//c                                        ** Eq. R30
-         CNUMER = CAK + one/CNUMER;
-         CDENOM = CAK + one/CDENOM;
-         continue;
-         //GO TO  10
-       }  // end of if < eps2
-      }
-      break;
-    } while(1);    
-    //if (N == 0)  printf(" return confra for z=(%g,%g)\n", ZINV.real(), ZINV.imag());
-    return CONFRA;
-  }
-  // ********************************************************************** //
-  // ********************************************************************** //
-  // ********************************************************************** //
-  void MultiLayerMieApplied::ConvertToSP() {
-    MarkUncalculated();
-    if (target_width_.size() + coating_width_.size() == 0)
-      return;  // Nothing to convert, we suppose that SP was set directly
-    GenerateSizeParameter();
-    GenerateIndex();
-    if (size_param_.size() != refractive_index_.size())
-      throw std::invalid_argument("Ivalid conversion of width to size parameter units!/n");
-  }
-  // ********************************************************************** //
-  // ********************************************************************** //
-  // ********************************************************************** //
-  void MultiLayerMieApplied::RunMieCalculation() {
-    ConvertToSP();
-    MultiLayerMie::RunMieCalculation(); 
-  }
-  // ********************************************************************** //
-  // ********************************************************************** //
-  // ********************************************************************** //
-  void MultiLayerMieApplied::GetExpanCoeffs( std::vector< std::vector<std::complex<double> > >& aln, std::vector< std::vector<std::complex<double> > >& bln, std::vector< std::vector<std::complex<double> > >& cln, std::vector< std::vector<std::complex<double> > >& dln) {
-    ConvertToSP();  // Case of call before running full Mie calculation.
-    // Calculate scattering coefficients an_ and bn_
-    calcScattCoeffs();
-    // Calculate expansion coefficients aln_,  bln_, cln_, and dln_
-    calcExpanCoeffs();
-    aln = aln_;
-    bln = bln_;
-    cln = cln_;
-    dln = dln_;
-    
-  }  // end of void MultiLayerMieApplied::GetExpanCoeffs( ...)
-  // ********************************************************************** //
-  // ********************************************************************** //
-  // ********************************************************************** //
 
 }  // end of namespace nmie

src/nmie-applied.hpp

@@ -33,6 +33,7 @@
 #include <iostream>
 #include <vector>
 #include "nmie.hpp"
+#include "nmie-impl.hpp"
 
 
 namespace nmie {
@@ -51,7 +52,7 @@ namespace nmie {
     void GetFailed();
     long iformat = 0;
     bool output = true;
-    void prn(double var) {
+    void prn(FloatType var) {
       do {
 	if (!output) break;
 	++iformat;
@@ -60,28 +61,28 @@ namespace nmie {
       } while (false);
     }
     // Set parameters in applied units 
-    void SetWavelength(double wavelength) {wavelength_ = wavelength;};
+    void SetWavelength(FloatType wavelength) {wavelength_ = wavelength;};
     // It is possible to set only a multilayer target to run calculaitons.
     // For many runs it can be convenient to separate target and coating layers.
     // Per layer
-    void AddTargetLayer(double layer_width, std::complex<double> layer_index);
-    void AddCoatingLayer(double layer_width, std::complex<double> layer_index);
+    void AddTargetLayer(FloatType layer_width, std::complex<FloatType> layer_index);
+    void AddCoatingLayer(FloatType layer_width, std::complex<FloatType> layer_index);
     // For all layers
-    void SetTargetWidth(std::vector<double> width);
-    void SetTargetIndex(std::vector< std::complex<double> > index);
-    void SetTargetPEC(double radius);
-    void SetCoatingWidth(std::vector<double> width);
-    void SetCoatingIndex(std::vector< std::complex<double> > index);
-    void SetFieldPoints(std::vector< std::array<double,3> > coords);
+    void SetTargetWidth(std::vector<FloatType> width);
+    void SetTargetIndex(std::vector< std::complex<FloatType> > index);
+    void SetTargetPEC(FloatType radius);
+    void SetCoatingWidth(std::vector<FloatType> width);
+    void SetCoatingIndex(std::vector< std::complex<FloatType> > index);
+    void SetFieldPoints(std::vector< std::array<FloatType,3> > coords);
 
     //Set parameters in size parameter units
-    void SetWidthSP(const std::vector<double>& width);
-    void SetIndexSP(const std::vector< std::complex<double> >& index);
-    void SetFieldPointsSP(const std::vector< std::vector<double> >& coords_sp);
+    void SetWidthSP(const std::vector<FloatType>& width);
+    void SetIndexSP(const std::vector< std::complex<FloatType> >& index);
+    void SetFieldPointsSP(const std::vector< std::vector<FloatType> >& coords_sp);
 
     // Set common parameters
-    void SetAnglesForPattern(double from_angle, double to_angle, int samples);
-    std::vector<double> GetAngles();
+    void SetAnglesForPattern(FloatType from_angle, FloatType to_angle, int samples);
+    std::vector<FloatType> GetAngles();
     
     void ClearTarget();
     void ClearCoating();
@@ -89,43 +90,43 @@ namespace nmie {
     void ClearAllDesign(); //Layers + SP + index_
 
     // Applied units requests
-    double GetTotalRadius();
-    double GetTargetRadius();
-    double GetCoatingWidth();
-    std::vector<double>                  GetTargetLayersWidth();
-    std::vector< std::complex<double> >  GetTargetLayersIndex();
-    std::vector<double>                  GetCoatingLayersWidth();
-    std::vector< std::complex<double> >  GetCoatingLayersIndex();
-    std::vector< std::vector<double> >   GetFieldPoints();
-    std::vector< std::vector<double> > GetSpectra(double from_WL, double to_WL, int samples);  // ext, sca, abs, bk
-    double GetRCSext();
-    double GetRCSsca();
-    double GetRCSabs();
-    double GetRCSbk();
-    std::vector<double> GetPatternEk();
-    std::vector<double> GetPatternHk();
-    std::vector<double> GetPatternUnpolarized();
+    FloatType GetTotalRadius();
+    FloatType GetTargetRadius();
+    FloatType GetCoatingWidth();
+    std::vector<FloatType>                  GetTargetLayersWidth();
+    std::vector< std::complex<FloatType> >  GetTargetLayersIndex();
+    std::vector<FloatType>                  GetCoatingLayersWidth();
+    std::vector< std::complex<FloatType> >  GetCoatingLayersIndex();
+    std::vector< std::vector<FloatType> >   GetFieldPoints();
+    std::vector< std::vector<FloatType> > GetSpectra(FloatType from_WL, FloatType to_WL, int samples);  // ext, sca, abs, bk
+    FloatType GetRCSext();
+    FloatType GetRCSsca();
+    FloatType GetRCSabs();
+    FloatType GetRCSbk();
+    std::vector<FloatType> GetPatternEk();
+    std::vector<FloatType> GetPatternHk();
+    std::vector<FloatType> GetPatternUnpolarized();
 
     // Size parameter units
-    std::vector<double> GetLayerWidthSP();
+    std::vector<FloatType> GetLayerWidthSP();
     // Same as to get target and coating index
-    std::vector< std::complex<double> > GetLayerIndex();  
-    std::vector< std::array<double,3> > GetFieldPointsSP();
+    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<double> > >  GetFieldESP(); */
-    /* std::vector<std::vector<std::complex<double> > >  GetFieldHSP(); */
-    std::vector< std::array<double,5> > GetSpectraSP(double from_SP, double to_SP, int samples);  // WL,ext, sca, abs, bk
+    /* std::vector<std::vector<std::complex<FloatType> > >  GetFieldESP(); */
+    /* std::vector<std::vector<std::complex<FloatType> > >  GetFieldHSP(); */
+    std::vector< std::array<FloatType,5> > GetSpectraSP(FloatType from_SP, FloatType to_SP, int samples);  // WL,ext, sca, abs, bk
 
 
-    std::vector<double> GetPatternEkSP();
-    std::vector<double> GetPatternHkSP();
-    std::vector<double> GetPatternUnpolarizedSP();
+    std::vector<FloatType> GetPatternEkSP();
+    std::vector<FloatType> GetPatternHkSP();
+    std::vector<FloatType> GetPatternUnpolarizedSP();
 
     void GetExpanCoeffs
-      (std::vector< std::vector<std::complex<double> > >& aln,
-       std::vector< std::vector<std::complex<double> > >& bln,
-       std::vector< std::vector<std::complex<double> > >& cln,
-       std::vector< std::vector<std::complex<double> > >& dln);
+      (std::vector< std::vector<std::complex<FloatType> > >& aln,
+       std::vector< std::vector<std::complex<FloatType> > >& bln,
+       std::vector< std::vector<std::complex<FloatType> > >& cln,
+       std::vector< std::vector<std::complex<FloatType> > >& dln);
 
 
     // Output results (data file + python script to plot it with matplotlib)
@@ -142,20 +143,20 @@ namespace nmie {
     void GenerateIndex();
     void InitMieCalculations();
 
-    void sbesjh(std::complex<double> z, std::vector<std::complex<double> >& jn,
-	            std::vector<std::complex<double> >& jnp, std::vector<std::complex<double> >& h1n,
-	            std::vector<std::complex<double> >& h1np);
-    void sphericalBessel(std::complex<double> z, std::vector<std::complex<double> >& bj,
-			             std::vector<std::complex<double> >& by, std::vector<std::complex<double> >& bd);
-    std::complex<double> calcD1confra(int N, const std::complex<double> z);
+    void sbesjh(std::complex<FloatType> z, std::vector<std::complex<FloatType> >& jn,
+	            std::vector<std::complex<FloatType> >& jnp, std::vector<std::complex<FloatType> >& h1n,
+	            std::vector<std::complex<FloatType> >& h1np);
+    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);
     
-    double wavelength_ = 1.0;
-    double total_radius_ = 0.0;
+    FloatType wavelength_ = 1.0;
+    FloatType total_radius_ = 0.0;
     /// Width and index for each layer of the structure
-    std::vector<double> target_width_, coating_width_;
-    std::vector< std::complex<double> > target_index_, coating_index_;
+    std::vector<FloatType> target_width_, coating_width_;
+    std::vector< std::complex<FloatType> > target_index_, coating_index_;
 
-    std::vector< std::vector<double> > coords_sp_;
+    std::vector< std::vector<FloatType> > coords_sp_;
 
 
 

tests/c++/go-speed-test.sh

@@ -14,8 +14,8 @@ rm -f $PROGRAM
 #g++ -Ofast -std=c++11 $file ../../src/nmie.cc  -lm -lrt -o $PROGRAM /usr/lib/libtcmalloc_minimal.so.4 -fno-builtin-malloc -fno-builtin-calloc -fno-builtin-realloc -fno-builtin-free -march=native -mtune=native -msse4.2
 
 file=speed-test-applied.cc
-g++ -Ofast -std=c++11 $file ../../src/nmie.cc ../../src/nmie-applied.cc -DMULTI_PRECISION=200 -lm -lrt -o $PROGRAM /usr/lib/libtcmalloc_minimal.so.4 -fno-builtin-malloc -fno-builtin-calloc -fno-builtin-realloc -fno-builtin-free -march=native -mtune=native -msse4.2
-# g++ -Ofast -std=c++11 $file ../../src/nmie.cc ../../src/nmie-applied.cc -lm -lrt -o $PROGRAM /usr/lib/libtcmalloc_minimal.so.4 -fno-builtin-malloc -fno-builtin-calloc -fno-builtin-realloc -fno-builtin-free -march=native -mtune=native -msse4.2
+# g++ -Ofast -std=c++11 $file ../../src/nmie.cc ../../src/nmie-applied.cc -DMULTI_PRECISION=200 -lm -lrt -o $PROGRAM /usr/lib/libtcmalloc_minimal.so.4 -fno-builtin-malloc -fno-builtin-calloc -fno-builtin-realloc -fno-builtin-free -march=native -mtune=native -msse4.2
+g++ -Ofast -std=c++11 $file ../../src/nmie.cc ../../src/nmie-applied.cc -lm -lrt -o $PROGRAM /usr/lib/libtcmalloc_minimal.so.4 -fno-builtin-malloc -fno-builtin-calloc -fno-builtin-realloc -fno-builtin-free -march=native -mtune=native -msse4.2
 
 echo Should be:
 echo test01, +1.41154e+00, +4.17695e-01, +9.93844e-01, +1.59427e-01, +1.25809e+00, +3.67376e-01, +2.95915e-01