Start of conversion to std::vector and std::complex

Makefile

@@ -31,7 +31,7 @@ builddeb:
 	dpkg-buildpackage -i -I -rfakeroot
 
 standalone: standalone.cc nmie.cc ucomplex.cc
-	cc standalone.cc nmie.cc ucomplex.cc -lm -o scattnlay
+	c++ -DNDEBUG -O2 -std=c++11 standalone.cc nmie.cc ucomplex.cc -lm -o scattnlay
 	mv scattnlay ../
 
 clean:

go.sh

@@ -1,13 +1,19 @@
 #!/bin/bash
 echo Compile with gcc -O2
 rm -rf *.bin
-g++ -O2 standalone.cc nmie.cc ucomplex.cc -lm -o scattnlay.bin
+
+clang++ -g -O1 -fsanitize=address  -fno-optimize-sibling-calls -fno-omit-frame-pointer -std=c++11 standalone.cc nmie.cc ucomplex.cc -lm -o scattnlay.bin
+#g++ -O2 -std=c++11 standalone.cc nmie.cc ucomplex.cc -lm -o scattnlay.bin
 cp scattnlay.bin ../scattnlay
 cd tests/shell
 # for file in `ls *.sh`; do ./$file; done
-repeats=30
-echo Run test for $repeats times
-time for i in `seq $repeats`; do ./test01.sh; done
-echo Run test with original binary for $repeats times
-cp ../../../scattnlay-0.3.0 ../../../scattnlay
-time for i in `seq $repeats`; do ./test01.sh; done
+PROGRAM='../../../scattnlay'
+ASAN_SYMBOLIZER_PATH=/usr/bin/llvm-symbolizer-3.4 $PROGRAM -l 5 0.4642 1.8000 1.7000 0.7114 0.8000 0.7000 0.7393 1.2000 0.0900 0.9168 2.8000 0.2000 1.0000 1.5000 0.4000 -c test01
+
+
+# repeats=30
+# echo Run test for $repeats times
+# time for i in `seq $repeats`; do ./test01.sh; done
+# echo Run test with original binary for $repeats times
+# cp ../../../scattnlay-0.3.0 ../../../scattnlay
+# time for i in `seq $repeats`; do ./test01.sh; done

nmie-wrapper.cc

@@ -0,0 +1,148 @@
+///
+/// @file   nmie-wrapper.cc
+/// @author Ladutenko Konstantin <kostyfisik at gmail (.) com>
+/// @date   Tue Sep  3 00:38:27 2013
+/// @copyright 2013 Ladutenko Konstantin
+///
+/// nmie-wrapper 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.
+///
+/// nmie-wrapper 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.
+///
+/// You should have received a copy of the GNU General Public License
+/// along with nmie-wrapper.  If not, see <http://www.gnu.org/licenses/>.
+///
+/// nmie-wrapper uses nmie.c from scattnlay by Ovidio Pena
+/// <ovidio@bytesfall.com> as a linked library. He has an additional condition to 
+/// his library:
+//    The only additional condition 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.                                       //
+///
+/// @brief  Wrapper class around nMie function for ease of use
+///
+#include "ucomplex.h"
+#include "nmie-wrapper.h"
+#include "nmie.h"
+#include <cstdio>
+#include <cstdlib>
+#include <stdexcept>
+#include <vector>
+namespace nmie {  
+  // ********************************************************************** //
+  // ********************************************************************** //
+  // ********************************************************************** //
+    void MultiLayerMie::SetSizeParameter(std::vector<double> size_parameter) {
+    size_parameter_.clear();
+    for (auto layer_size_parameter : size_parameter)
+      if (layer_size_parameter <= 0.0)
+        throw std::invalid_argument("Size parameter should be positive!");
+      else size_parameter_.push_back(layer_size_parameter);
+  }
+  // end of void MultiLayerMie::SetSizeParameter(...);
+  // ********************************************************************** //
+  // ********************************************************************** //
+  // ********************************************************************** //
+  void MultiLayerMie::SetIndex(std::vector< std::complex<double> > index) {
+    index_.clear();
+    for (auto value : index) index_.push_back(value);
+  }  // end of void MultiLayerMie::SetIndex(...);
+  // ********************************************************************** //
+  // ********************************************************************** //
+  // ********************************************************************** //
+  /// nMie starts layer numeration from 1 (no separation of target
+  /// and coating). So the first elment (zero-indexed) is not used
+  /// and has some unused value.
+  void MultiLayerMie::GenerateSizeParameter() {
+    // size_parameter_.clear();
+    // size_parameter_.push_back(0.0);
+    // double radius = 0.0;
+    // for (auto width : target_thickness_) {
+    //   radius += width;
+    //   size_parameter_.push_back(2*PI*radius / wavelength_);
+    // }
+    // for (auto width : coating_thickness_) {
+    //   radius += width;
+    //   size_parameter_.push_back(2*PI*radius / wavelength_);
+    // }
+    // total_radius_ = radius;
+  }  // end of void MultiLayerMie::GenerateSizeParameter();
+  // ********************************************************************** //
+  // ********************************************************************** //
+  // ********************************************************************** //
+  void MultiLayerMie::RunMie(double *Qext_out, double *Qsca_out,
+                             double *Qabs_out, double *Qbk_out) {
+    if (size_parameter_.size() != index_.size())
+      throw std::invalid_argument("Each size parameter should have only one index!");
+    int L = static_cast<int>(size_parameter_.size()) - 1;
+    double Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo;
+    int nt = 0;
+    double *Theta = (double *)malloc(nt * sizeof(double));
+    complex *S1 = (complex *)malloc(nt * sizeof(complex *));
+    complex *S2 = (complex *)malloc(nt * sizeof(complex *));
+    double *x = &(size_parameter_.front());
+    complex *m = &(index_.front());
+    int terms = 0;
+    terms = nMieFast(L, x, m, nt, Theta,
+                 &Qext, &Qsca, &Qabs, &Qbk, &Qpr, &g, &Albedo,
+                 S1,S2);
+    free(Theta);
+    free(S1);
+    free(S2);
+    if (terms == 0) {
+      *Qext_out = Qfaild_;
+      *Qsca_out = Qfaild_;
+      *Qabs_out = Qfaild_;
+      *Qbk_out = Qfaild_;
+      throw std::invalid_argument("Failed to evaluate Q!");
+    }
+    *Qext_out = Qext;
+    *Qsca_out = Qsca;
+    *Qabs_out = Qabs;
+    *Qbk_out = Qbk;
+  }  // end of void MultiLayerMie::RunMie();
+  // ********************************************************************** //
+  // ********************************************************************** //
+  // ********************************************************************** //
+  double MultiLayerMie::GetTotalRadius() {
+    if (total_radius_ == 0) GenerateSizeParameter();
+    return total_radius_;      
+  }  // end of double MultiLayerMie::GetTotalRadius();
+  // ********************************************************************** //
+  // ********************************************************************** //
+  // ********************************************************************** //
+  std::vector< std::vector<double> >
+  MultiLayerMie::GetSpectra(double from_WL, double to_WL, int samples) {
+    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) {
+      double Qext, Qsca, Qabs, Qbk;
+      wavelength_ = WL;
+      try {
+        RunMie(&Qext, &Qsca, &Qabs, &Qbk);
+      } catch( const std::invalid_argument& ia ) {
+        fails++;
+        continue;
+      }
+      //printf("%3.1f ",WL);
+      spectra.push_back({wavelength_, Qext, Qsca, Qabs, Qbk});
+    }  // end of for each WL in spectra
+    printf("fails %li\n",fails);
+    wavelength_ = wavelength_backup;
+    return spectra;
+  }
+  // ********************************************************************** //
+  // ********************************************************************** //
+  // ********************************************************************** //
+///MultiLayerMie::
+}  // end of namespace nmie

nmie-wrapper.h

@@ -0,0 +1,77 @@
+#ifndef SRC_NMIE_NMIE_WRAPPER_H_
+#define SRC_NMIE_NMIE_WRAPPER_H_
+///
+/// @file   nmie-wrapper.h
+/// @author Ladutenko Konstantin <kostyfisik at gmail (.) com>
+/// @date   Tue Sep  3 00:40:47 2013
+/// @copyright 2013 Ladutenko Konstantin
+///
+/// nmie-wrapper 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.
+///
+/// nmie-wrapper 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.
+///
+/// You should have received a copy of the GNU General Public License
+/// along with nmie-wrapper.  If not, see <http://www.gnu.org/licenses/>.
+///
+/// nmie-wrapper uses nmie.c from scattnlay by Ovidio Pena
+/// <ovidio@bytesfall.com> as a linked library. He has an additional condition to 
+/// his library:
+//    The only additional condition 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.                                       //
+///
+/// @brief  Wrapper class around nMie function for ease of use
+/// 
+///
+#include <complex>
+#include <cstdlib>
+#include <iostream>
+#include <vector>
+
+#ifndef NDEBUG
+#   define ASSERT(condition, message) \
+    do { \
+        if (! (condition)) { \
+            std::cerr << "Assertion `" #condition "` failed in " << __FILE__ \
+                      << " line " << __LINE__ << ": " << message << std::endl; \
+            std::exit(EXIT_FAILURE); \
+        } \
+    } while (false)
+#else
+#   define ASSERT(condition, message) do { } while (false)
+#endif
+
+namespace nmie {
+  class MultiLayerMie {
+   public:
+    void SetWavelength(double wavelength) {wavelength_ = wavelength;};
+    void SetSizeParameter(std::vector<double> size_parameter);
+    void SetIndex(std::vector< std::complex<double> > index);
+    void RunMie(double *Qext, double *Qsca, double *Qabs, double *Qbk);
+    std::vector< std::vector<double> >  GetSpectra(double from_WL, double to_WL,
+                                                   int samples);
+    /// Disabled functions or unused untill the end of C++ conversion
+    double GetTotalRadius();
+  private:
+    /// Size parameter for each layer.
+    std::vector<double> size_parameter_;
+    /// Complex index value for each layer.
+    std::vector< std::complex<double> > index_;
+    const double PI=3.14159265358979323846;
+    /// Disabled functions or unused untill the end of C++ conversion
+    void GenerateSizeParameter();
+    double wavelength_ = 1.0;
+    double total_radius_ = 0.0;
+
+  };  // end of class MultiLayerMie
+}  // end of namespace nmie
+#endif  // SRC_NMIE_NMIE_WRAPPER_H_

nmie.cc

@@ -503,6 +503,232 @@ int ScattCoeff(int L, int pl, double x[], complex m[], int n_max, complex **an,
 
   return n_max;
 }
+//**********************************************************************************//
+//**********************************************************************************//
+//**********************************************************************************//
+int ScattCoeff_std(double x[], complex m[], complex **an, complex **bn, 
+		   int L, int pl, std::vector<double> x_std,
+		   std::vector<std::complex<double> > m_std, int n_max,
+		   std::vector< std::complex<double> > an_std, 
+		   std::vector< std::complex<double> > bn_std){
+  //************************************************************************//
+  // Calculate the index of the first layer. It can be either 0 (default)   //
+  // or the index of the outermost PEC layer. In the latter case all layers //
+  // below the PEC are discarded.                                           //
+  //************************************************************************//
+  int fl = firstLayer(L, pl);
+
+  if (n_max <= 0) {
+    n_max = Nmax(L, fl, pl, x, m);
+  }
+  
+  complex z1, z2, cn;
+  complex Num, Denom;
+  complex G1, G2;
+  complex Temp;
+  double Tmp;
+
+  int n, l, t;
+
+  //**************************************************************************//
+  // Note that since Fri, Nov 14, 2014 all arrays start from 0 (zero), which  //
+  // means that index = layer number - 1 or index = n - 1. The only exception //
+  // are the arrays for representing D1, D3 and Q because they need a value   //
+  // for the index 0 (zero), hence it is important to consider this shift     //
+  // between different arrays. The change was done to optimize memory usage.  //
+  //**************************************************************************//
+
+  // Allocate memory to the arrays
+  complex **D1_mlxl = (complex **) malloc(L*sizeof(complex *));
+  complex **D1_mlxlM1 = (complex **) malloc(L*sizeof(complex *));
+
+  complex **D3_mlxl = (complex **) malloc(L*sizeof(complex *));
+  complex **D3_mlxlM1 = (complex **) malloc(L*sizeof(complex *));
+
+  complex **Q = (complex **) malloc(L*sizeof(complex *));
+
+  complex **Ha = (complex **) malloc(L*sizeof(complex *));
+  complex **Hb = (complex **) malloc(L*sizeof(complex *));
+
+  for (l = 0; l < L; l++) {
+    D1_mlxl[l] = (complex *) malloc((n_max + 1)*sizeof(complex));
+    D1_mlxlM1[l] = (complex *) malloc((n_max + 1)*sizeof(complex));
+
+    D3_mlxl[l] = (complex *) malloc((n_max + 1)*sizeof(complex));
+    D3_mlxlM1[l] = (complex *) malloc((n_max + 1)*sizeof(complex));
+
+    Q[l] = (complex *) malloc((n_max + 1)*sizeof(complex));
+
+    Ha[l] = (complex *) malloc(n_max*sizeof(complex));
+    Hb[l] = (complex *) malloc(n_max*sizeof(complex));
+  }
+
+  (*an) = (complex *) malloc(n_max*sizeof(complex));
+  (*bn) = (complex *) malloc(n_max*sizeof(complex));
+
+  complex *D1XL = (complex *) malloc((n_max + 1)*sizeof(complex));
+  complex *D3XL = (complex *) malloc((n_max + 1)*sizeof(complex));
+
+  complex *PsiXL = (complex *) malloc((n_max + 1)*sizeof(complex));
+  complex *ZetaXL = (complex *) malloc((n_max + 1)*sizeof(complex));
+
+  //*************************************************//
+  // Calculate D1 and D3 for z1 in the first layer   //
+  //*************************************************//
+  if (fl == pl) {  // PEC layer
+    for (n = 0; n <= n_max; n++) {
+      D1_mlxl[fl][n] = Complex(0, -1);
+      D3_mlxl[fl][n] = C_I;
+    }
+  } else { // Regular layer
+    z1 = RCmul(x[fl], m[fl]);
+
+    // Calculate D1 and D3
+    calcD1D3(z1, n_max, &(D1_mlxl[fl]), &(D3_mlxl[fl]));
+  }
+
+  //******************************************************************//
+  // Calculate Ha and Hb in the first layer - equations (7a) and (8a) //
+  //******************************************************************//
+  for (n = 0; n < n_max; n++) {
+    Ha[fl][n] = D1_mlxl[fl][n + 1];
+    Hb[fl][n] = D1_mlxl[fl][n + 1];
+  }
+
+  //*****************************************************//
+  // Iteration from the second layer to the last one (L) //
+  //*****************************************************//
+  for (l = fl + 1; l < L; l++) {
+    //************************************************************//
+    //Calculate D1 and D3 for z1 and z2 in the layers fl+1..L     //
+    //************************************************************//
+    z1 = RCmul(x[l], m[l]);
+    z2 = RCmul(x[l - 1], m[l]);
+
+    //Calculate D1 and D3 for z1
+    calcD1D3(z1, n_max, &(D1_mlxl[l]), &(D3_mlxl[l]));
+
+    //Calculate D1 and D3 for z2
+    calcD1D3(z2, n_max, &(D1_mlxlM1[l]), &(D3_mlxlM1[l]));
+
+    //*********************************************//
+    //Calculate Q, Ha and Hb in the layers fl+1..L //
+    //*********************************************//
+
+    // Upward recurrence for Q - equations (19a) and (19b)
+    Num = RCmul(exp(-2*(z1.i - z2.i)), Complex(cos(-2*z2.r) - exp(-2*z2.i), sin(-2*z2.r)));
+    Denom = Complex(cos(-2*z1.r) - exp(-2*z1.i), sin(-2*z1.r));
+    Q[l][0] = Cdiv(Num, Denom);
+
+    for (n = 1; n <= n_max; n++) {
+      cn = Complex(n, 0);
+      Num = Cmul(Cadd(Cmul(z1, D1_mlxl[l][n]), cn), Csub(cn, Cmul(z1, D3_mlxl[l][n - 1])));
+      Denom = Cmul(Cadd(Cmul(z2, D1_mlxlM1[l][n]), cn), Csub(cn, Cmul(z2, D3_mlxlM1[l][n - 1])));
+
+      Q[l][n] = Cdiv(Cmul(RCmul((x[l - 1]*x[l - 1])/(x[l]*x[l]), Q[l][n - 1]), Num), Denom);
+    }
+
+    // Upward recurrence for Ha and Hb - equations (7b), (8b) and (12) - (15)
+    for (n = 1; n <= n_max; n++) {
+      //Ha
+      if ((l - 1) == pl) { // The layer below the current one is a PEC layer
+        G1 = RCmul(-1.0, D1_mlxlM1[l][n]);
+        G2 = RCmul(-1.0, D3_mlxlM1[l][n]);
+      } else {
+        G1 = Csub(Cmul(m[l], Ha[l - 1][n - 1]), Cmul(m[l - 1], D1_mlxlM1[l][n]));
+        G2 = Csub(Cmul(m[l], Ha[l - 1][n - 1]), Cmul(m[l - 1], D3_mlxlM1[l][n]));
+      }
+
+      Temp = Cmul(Q[l][n], G1);
+
+      Num = Csub(Cmul(G2, D1_mlxl[l][n]), Cmul(Temp, D3_mlxl[l][n]));
+      Denom = Csub(G2, Temp);
+
+      Ha[l][n - 1] = Cdiv(Num, Denom);
+
+      //Hb
+      if ((l - 1) == pl) { // The layer below the current one is a PEC layer
+        G1 = Hb[l - 1][n - 1];
+        G2 = Hb[l - 1][n - 1];
+      } else {
+        G1 = Csub(Cmul(m[l - 1], Hb[l - 1][n - 1]), Cmul(m[l], D1_mlxlM1[l][n]));
+        G2 = Csub(Cmul(m[l - 1], Hb[l - 1][n - 1]), Cmul(m[l], D3_mlxlM1[l][n]));
+      }
+
+      Temp = Cmul(Q[l][n], G1);
+
+      Num = Csub(Cmul(G2, D1_mlxl[l][n]), Cmul(Temp, D3_mlxl[l][n]));
+      Denom = Csub(G2, Temp);
+
+      Hb[l][n - 1] = Cdiv(Num, Denom);
+    }
+  }
+
+  //**************************************//
+  //Calculate D1, D3, Psi and Zeta for XL //
+  //**************************************//
+  z1 = Complex(x[L - 1], 0);
+
+  // Calculate D1XL and D3XL
+  calcD1D3(z1, n_max, &D1XL, &D3XL);
+
+  // Calculate PsiXL and ZetaXL
+  calcPsiZeta(z1, n_max, D1XL, D3XL, &PsiXL, &ZetaXL);
+
+  //*********************************************************************//
+  // Finally, we calculate the scattering coefficients (an and bn) and   //
+  // the angular functions (Pi and Tau). Note that for these arrays the  //
+  // first layer is 0 (zero), in future versions all arrays will follow  //
+  // this convention to save memory. (13 Nov, 2014)                      //
+  //*********************************************************************//
+  for (n = 0; n < n_max; n++) {
+    //********************************************************************//
+    //Expressions for calculating an and bn coefficients are not valid if //
+    //there is only one PEC layer (ie, for a simple PEC sphere).          //
+    //********************************************************************//
+    if (pl < (L - 1)) {
+      (*an)[n] = calc_an(n + 1, x[L - 1], Ha[L - 1][n], m[L - 1], PsiXL[n + 1], ZetaXL[n + 1], PsiXL[n], ZetaXL[n]);
+      (*bn)[n] = calc_bn(n + 1, x[L - 1], Hb[L - 1][n], m[L - 1], PsiXL[n + 1], ZetaXL[n + 1], PsiXL[n], ZetaXL[n]);
+    } else {
+      (*an)[n] = calc_an(n + 1, x[L - 1], C_ZERO, C_ONE, PsiXL[n + 1], ZetaXL[n + 1], PsiXL[n], ZetaXL[n]);
+      (*bn)[n] = Cdiv(PsiXL[n + 1], ZetaXL[n + 1]);
+    }
+  }
+
+  // Free the memory used for the arrays
+  for (l = 0; l < L; l++) {
+    free(D1_mlxl[l]);
+    free(D1_mlxlM1[l]);
+
+    free(D3_mlxl[l]);
+    free(D3_mlxlM1[l]);
+
+    free(Q[l]);
+
+    free(Ha[l]);
+    free(Hb[l]);
+  }
+
+  free(D1_mlxl);
+  free(D1_mlxlM1);
+
+  free(D3_mlxl);
+  free(D3_mlxlM1);
+
+  free(Q);
+
+  free(Ha);
+  free(Hb);
+
+  free(D1XL);
+  free(D3XL);
+
+  free(PsiXL);
+  free(ZetaXL);
+
+  return n_max;
+}
+
 
 //**********************************************************************************//
 // This function is just a wrapper to call the function 'nMieScatt' with fewer      //
@@ -516,6 +742,13 @@ int nMie(int L, double x[], complex m[], int nTheta, double Theta[], double *Qex
   return nMieScatt(L, -1, x, m, nTheta, Theta, -1, Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo, S1, S2);
 }
 
+
+int nMie_std(double x[], complex m[], double Theta[], complex S1[], complex S2[],
+int L, std::vector<double> &x_std, std::vector<std::complex<double> > &m_std, int nTheta, std::vector<double> &Theta_std, double *Qext, double *Qsca, double *Qabs, double *Qbk, double *Qpr, double *g, double *Albedo, std::vector< std::complex<double> > &S1_std, std::vector< std::complex<double> >  &S2_std) {
+  return nMieScatt_std(x, m, Theta, S1, S2, L, -1, x_std, m_std, nTheta, Theta_std, -1, Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo, S1_std, S2_std);
+}
+
+
 //**********************************************************************************//
 // This function is just a wrapper to call the function 'nMieScatt' with fewer      //
 // parameters, it is useful if you want to include a PEC layer but not a limit      //
@@ -698,6 +931,105 @@ int nMieScatt(int L, int pl, double x[], complex m[], int nTheta, double Theta[]
   return n_max;
 }
 
+
+
+int nMieScatt_std(double x[], complex m[], double Theta[], complex S1[], complex S2[],
+		  int L, int pl,
+		  std::vector<double> &x_std, std::vector<std::complex<double> > &m_std,
+		  int nTheta, std::vector<double> &Theta_std,
+		  int n_max, double *Qext, double *Qsca, double *Qabs, double *Qbk,
+		  double *Qpr, double *g, double *Albedo,
+		  std::vector< std::complex<double> > &S1_std,
+		  std::vector< std::complex<double> >  &S2_std)  {
+  int i, n, t;
+  double **Pi, **Tau;
+  std::vector< std::vector<double> > Pi_std, Tau_std;
+  complex *an, *bn;
+  std::vector< std::complex<double> > an_std, bn_std;
+  complex Qbktmp;
+  std::complex<double> Qbktmp_std;
+  {
+    int tmp_n_max = ScattCoeff(L, pl, x, m, n_max, &an, &bn);
+    n_max = ScattCoeff(L, pl, x, m, n_max, &an, &bn);
+  }
+  Pi = (double **) malloc(n_max*sizeof(double *));
+  Tau = (double **) malloc(n_max*sizeof(double *));
+  for (n = 0; n < n_max; n++) {
+    Pi[n] = (double *) malloc(nTheta*sizeof(double));
+    Tau[n] = (double *) malloc(nTheta*sizeof(double));
+  }
+
+  calcPiTau(n_max, nTheta, Theta, &Pi, &Tau);
+
+  double x2 = x[L - 1]*x[L - 1];
+
+  // Initialize the scattering parameters
+  *Qext = 0;
+  *Qsca = 0;
+  *Qabs = 0;
+  *Qbk = 0;
+  Qbktmp = C_ZERO;
+  *Qpr = 0;
+  *g = 0;
+  *Albedo = 0;
+
+  // Initialize the scattering amplitudes
+  for (t = 0; t < nTheta; t++) {
+    S1[t] = C_ZERO;
+    S2[t] = C_ZERO;
+  }
+
+  // By using downward recurrence we avoid loss of precision due to float rounding errors
+  // See: https://docs.oracle.com/cd/E19957-01/806-3568/ncg_goldberg.html
+  //      http://en.wikipedia.org/wiki/Loss_of_significance
+  for (i = n_max - 2; i >= 0; i--) {
+    n = i + 1;
+    // Equation (27)
+    *Qext = *Qext + (double)(n + n + 1)*(an[i].r + bn[i].r);
+    // Equation (28)
+    *Qsca = *Qsca + (double)(n + n + 1)*(an[i].r*an[i].r + an[i].i*an[i].i + bn[i].r*bn[i].r + bn[i].i*bn[i].i);
+    // Equation (29)
+    *Qpr = *Qpr + ((n*(n + 2)/(n + 1))*((Cadd(Cmul(an[i], Conjg(an[n])), Cmul(bn[i], Conjg(bn[n])))).r) + ((double)(n + n + 1)/(n*(n + 1)))*(Cmul(an[i], Conjg(bn[i])).r));
+
+    // Equation (33)
+    Qbktmp = Cadd(Qbktmp, RCmul((double)((n + n + 1)*(1 - 2*(n % 2))), Csub(an[i], bn[i])));
+
+    //****************************************************//
+    // Calculate the scattering amplitudes (S1 and S2)    //
+    // Equations (25a) - (25b)                            //
+    //****************************************************//
+    for (t = 0; t < nTheta; t++) {
+      S1[t] = Cadd(S1[t], calc_S1_n(n, an[i], bn[i], Pi[i][t], Tau[i][t]));
+      S2[t] = Cadd(S2[t], calc_S2_n(n, an[i], bn[i], Pi[i][t], Tau[i][t]));
+    }
+  }
+
+  *Qext = 2*(*Qext)/x2;                                 // Equation (27)
+  *Qsca = 2*(*Qsca)/x2;                                 // Equation (28)
+  *Qpr = *Qext - 4*(*Qpr)/x2;                           // Equation (29)
+
+  *Qabs = *Qext - *Qsca;                                // Equation (30)
+  *Albedo = *Qsca / *Qext;                              // Equation (31)
+  *g = (*Qext - *Qpr) / *Qsca;                          // Equation (32)
+
+  *Qbk = (Qbktmp.r*Qbktmp.r + Qbktmp.i*Qbktmp.i)/x2;    // Equation (33)
+
+  // Free the memory used for the arrays
+  for (n = 0; n < n_max; n++) {
+    free(Pi[n]);
+    free(Tau[n]);
+  }
+
+  free(Pi);
+  free(Tau);
+
+  free(an);
+  free(bn);
+
+  return n_max;
+}
+
+
 //**********************************************************************************//
 // This function calculates complex electric and magnetic field in the surroundings //
 // and inside (TODO) the particle.                                                  //

nmie.h

@@ -24,7 +24,10 @@
 //    along with this program.  If not, see <http://www.gnu.org/licenses/>.         //
 //**********************************************************************************//
 
-#define VERSION "0.3.0"
+#define VERSION "0.3.1"
+#include <complex>
+#include <string>
+#include <vector>
 
 int nMie(int L, double x[], complex m[], int nTheta, double Theta[], double *Qext, double *Qsca, double *Qabs, double *Qbk, double *Qpr, double *g, double *Albedo, complex S1[], complex S2[]);
 
@@ -36,3 +39,21 @@ int nMieScatt(int L, int pl, double x[], complex m[], int nTheta, double Theta[]
 
 int nMieField(int L, int pl, double x[], complex m[], int n_max, int nCoords, double Xp[], double Yp[], double Zp[], complex E[], complex H[]);
 
+
+int nMie_std(double x[], complex m[], double Theta[], complex S1[], complex S2[],
+int L, std::vector<double> &x_std, std::vector<std::complex<double> > &m_std, int nTheta, std::vector<double> &Theta_std, double *Qext, double *Qsca, double *Qabs, double *Qbk, double *Qpr, double *g, double *Albedo, std::vector< std::complex<double> > &S1_std, std::vector< std::complex<double> >  &S2_std);
+
+
+int nMieScatt_std(double x[], complex m[], double Theta[], complex S1[], complex S2[],
+		  int L, int pl,
+		  std::vector<double> &x_std, std::vector<std::complex<double> > &m_std,
+		  int nTheta, std::vector<double> &Theta_std,
+		  int n_max, double *Qext, double *Qsca, double *Qabs, double *Qbk,
+		  double *Qpr, double *g, double *Albedo,
+		  std::vector< std::complex<double> > &S1_std,
+		  std::vector< std::complex<double> >  &S2_std);
+
+
+
+
+

push-to-github.sh

@@ -1,4 +1,4 @@
 #!/bin/bash
 # after adding origin with:
 #   git remote add origin git@github.com:kostyfisik/scattnlay.git
-git push -u origin master
+git push --tags -u origin master

standalone.cc

@@ -24,14 +24,21 @@
 //    along with this program.  If not, see <http://www.gnu.org/licenses/>.         //
 //**********************************************************************************//
 
-#include <vector>
+#include <algorithm>
+#include <complex>
+#include <functional>
+#include <iostream>
+#include <stdexcept>
 #include <string>
+#include <vector>
 #include <stdlib.h>
 #include <stdio.h>
 #include <time.h>
 #include <string.h>
 #include "ucomplex.h"
 #include "nmie.h"
+#include "nmie-wrapper.h"
+
 //#define MAXLAYERS 1100
 const int MAXLAYERS = 1100;
 //#define MAXTHETA 800
@@ -58,87 +65,134 @@ const double PI=3.14159265358979323846;
 //        'comment, Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo'                           //
 //***********************************************************************************//
 int main(int argc, char *argv[]) {
-  char comment[200];
-  std::string comment_std;
-  int has_comment = 0;
-  int i, j, L;
-  double *x, *Theta;
-  std::vector<double> x_std, Theta_std;
-  complex *m, *S1, *S2;
-  double Qext, Qabs, Qsca, Qbk, Qpr, g, Albedo;
-  double ti = 0.0, tf = 90.0;
-  int nt = 0;
-  
-  if (argc < 5) {
-    printf("Insufficient parameters.\nUsage: %s -l Layers x1 m1.r m1.i [x2 m2.r m2.i ...] [-t ti tf nt] [-c comment]\n", argv[0]);
-    return -1;
-  }
-
-  strcpy(comment, "");
-  for (i = 1; i < argc; i++) {
-    if (strcmp(argv[i], "-l") == 0) {
-      i++;
-      L = atoi(argv[i]);
-      x = (double *) malloc(L*sizeof(double));
-      m = (complex *) malloc(L*sizeof(complex));
-      if (argc < 3*(L + 1)) {
-        printf("Insufficient parameters.\nUsage: %s -l Layers x1 m1.r m1.i [x2 m2.r m2.i ...] [-t ti tf nt] [-c comment]\n", argv[0]);
-        return -1;
-      } else {
-        for (j = 0; j < L; j++) {
-          i++;
-          x[j] = atof(argv[i]);
-          i++;
-          m[j].r = atof(argv[i]);
-          i++;
-          m[j].i = atof(argv[i]);
-        }
+  try {
+    char comment[200];
+    std::string comment_std;
+    int has_comment = 0;
+    int i, j, L;
+    double *x, *Theta;
+    std::vector<double> x_std, Theta_std;
+    complex *m, *S1, *S2;
+    std::vector<std::complex<double> > m_std, S1_std, S2_std;
+    double Qext, Qabs, Qsca, Qbk, Qpr, g, Albedo;
+    double ti = 0.0, tf = 90.0;
+    int nt = 0;
+    
+    if (argc < 5) {
+      printf("Insufficient parameters.\n");
+      printf("Usage: %s -l Layers x1 m1.r m1.i [x2 m2.r m2.i ...] [-t ti tf nt] [-c comment]\n", argv[0]);
+      return -1;
+    }
+    
+    strcpy(comment, "");
+    for (i = 1; i < argc; i++) {
+      if (strcmp(argv[i], "-l") == 0) {
+	i++;
+	L = atoi(argv[i]);
+	x = (double *) malloc(L*sizeof(double));
+	x_std.resize(L);
+	m = (complex *) malloc(L*sizeof(complex));
+	m_std.resize(L);
+	if (argc < 3*(L + 1)) {
+	  printf("Insufficient parameters.\nUsage: %s -l Layers x1 m1.r m1.i [x2 m2.r m2.i ...] [-t ti tf nt] [-c comment]\n", argv[0]);
+	  return -1;
+	} else {
+	  for (j = 0; j < L; j++) {
+	    i++;
+	    x[j] = atof(argv[i]);
+	    x_std[j] = x[j];
+	    i++;
+	    m[j].r = atof(argv[i]);
+	    i++;
+	    m[j].i = atof(argv[i]);
+	    m_std[j] = std::complex<double>(m[j].r, m[j].i);
+	  }
+	}
+      } else if (strcmp(argv[i], "-t") == 0) {
+	i++;
+	ti = atof(argv[i]);
+	i++;
+	tf = atof(argv[i]);
+	i++;
+	nt = atoi(argv[i]);
+	Theta = (double *) malloc(nt*sizeof(double));
+	Theta_std.resize(nt);
+	S1 = (complex *) malloc(nt*sizeof(complex));
+	S2 = (complex *) malloc(nt*sizeof(complex));
+	S1_std.resize(nt);
+	S2_std.resize(nt);
+      } else if (strcmp(argv[i], "-c") == 0) {
+	i++;
+	strcpy(comment, argv[i]);
+	comment_std = std::string(comment);
+	has_comment = 1;
+      } else { i++; }
+    }
+    
+    if (nt < 0) {
+      printf("Error reading Theta.\n");
+      return -1;
+    } else if (nt == 1) {
+      Theta[0] = ti*PI/180.0;
+      Theta_std[0] = Theta[0];
+    } else {
+      for (i = 0; i < nt; i++) {
+	Theta[i] = (ti + (double)i*(tf - ti)/(nt - 1))*PI/180.0;
+	Theta_std[i] = Theta[i];
       }
-    } else if (strcmp(argv[i], "-t") == 0) {
-      i++;
-      ti = atof(argv[i]);
-      i++;
-      tf = atof(argv[i]);
-      i++;
-      nt = atoi(argv[i]);
-      Theta = (double *) malloc(nt*sizeof(double));
-      S1 = (complex *) malloc(nt*sizeof(complex));
-      S2 = (complex *) malloc(nt*sizeof(complex));
-    } else if (strcmp(argv[i], "-c") == 0) {
-      i++;
-      strcpy(comment, argv[i]);
-      has_comment = 1;
-    } else { i++; }
-  }
-
-  if (nt < 0) {
-    printf("Error reading Theta.\n");
-    return -1;
-  } else if (nt == 1) {
-    Theta[0] = ti*PI/180.0;
-  } else {
-    for (i = 0; i < nt; i++) {
-      Theta[i] = (ti + (double)i*(tf - ti)/(nt - 1))*PI/180.0;
     }
-  }
-
-  //  printf("Debug run!\n");
-  nMie(L, x, m, nt, Theta, &Qext, &Qsca, &Qabs, &Qbk, &Qpr, &g, &Albedo, S1, S2);
-
-  // if (has_comment) {
-  //   printf("%6s, %+.5e, %+.5e, %+.5e, %+.5e, %+.5e, %+.5e, %+.5e\n", comment, Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo);
-  // } else {
-  //   printf("%+.5e, %+.5e, %+.5e, %+.5e, %+.5e, %+.5e, %+.5e\n", Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo);
-  // }
-
-  // if (nt > 0) {
-  //   printf(" Theta,         S1.r,         S1.i,         S2.r,         S2.i\n");
+    
+    //  printf("Debug run!\n");
+    nMie(L, x, m, nt, Theta, &Qext, &Qsca, &Qabs, &Qbk, &Qpr, &g, &Albedo, S1, S2);
+    std::vector<double> old_result({Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo});
+    // if (has_comment) {
+    //   printf("%6s, %+.5e, %+.5e, %+.5e, %+.5e, %+.5e, %+.5e, %+.5e\n", comment, Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo);
+    // } else {
+    //   printf("%+.5e, %+.5e, %+.5e, %+.5e, %+.5e, %+.5e, %+.5e\n", Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo);
+    // }
+    
+    // if (nt > 0) {
+    //   printf(" Theta,         S1.r,         S1.i,         S2.r,         S2.i\n");
+      
+    //   for (i = 0; i < nt; i++) {
+    //     printf("%6.2f, %+.5e, %+.5e, %+.5e, %+.5e\n", Theta[i]*180.0/PI, S1[i].r, S1[i].i, S2[i].r, S2[i].i);
+    //   }
+    // }
+    /////////////////////////////////////////////////////////////////////////
+    // After conversion
+    nMie_std( x, m, Theta, S1, S2,
+	      L, x_std, m_std, nt, Theta_std, &Qext, &Qsca, &Qabs,
+	      &Qbk, &Qpr, &g, &Albedo, S1_std, S2_std);
+    std::vector<double> new_result({Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo});
+    // if (has_comment) {
+    //   printf("%6s, %+.5e, %+.5e, %+.5e, %+.5e, %+.5e, %+.5e, %+.5e **After\n", comment, Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo);
+    // } else {
+    //   printf("%+.5e, %+.5e, %+.5e, %+.5e, %+.5e, %+.5e, %+.5e **After\n", Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo);
+    // }
+    
+    // if (nt > 0) {
+    //   printf(" Theta,         S1.r,         S1.i,         S2.r,         S2.i **After\n");
+      
+    //   for (i = 0; i < nt; i++) {
+    //     printf("%6.2f, %+.5e, %+.5e, %+.5e, %+.5e **After\n", Theta[i]*180.0/PI, S1[i].r, S1[i].i, S2[i].r, S2[i].i);
+    //   }
+    // }
+    for (int i =0; i < old_result.size(); ++i) {
+      double diff = new_result[i] - old_result[i];
+      // printf("%g ", diff);
+      if (std::abs(diff) > 1e-16) printf(" ********* WARNING!!! diff = %g ********* \n", diff);
+    }
+    // std::vector<double> diff_result(old_result.size(), 0.0);
+    // std::transform(new_result.begin(), new_result.end(), old_result.begin(),
+    // 		   std::back_inserter(diff_result), std::minus<double>());
+    // std::cout << "diff: " <<  diff_result << std::endl;
 
-  //   for (i = 0; i < nt; i++) {
-  //     printf("%6.2f, %+.5e, %+.5e, %+.5e, %+.5e\n", Theta[i]*180.0/PI, S1[i].r, S1[i].i, S2[i].r, S2[i].i);
-  //   }
-  // }
-  return 0;
+  } 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;
 }
 
 

tests/dev/test-complex-lib.cc

@@ -11,7 +11,7 @@ void output(std::string operation, complex c1, std::complex<double> c2) {
 	   c1.r, c1.i,
 	   c2.real(), c2.imag());
     // Epsilon for double is 1e-16, check it 1.0 + 1e-16 = 1.0
-    if (diff_r > 1e-16||diff_i > 1e-16) 
+    if (std::abs(diff_r) > 1e-16||std::abs(diff_i) > 1e-16) 
 	//if (diff_r > 1e-15||diff_i > 1e-15) 
 	printf("\n********\n\t\tWARNING!! Non-zero diff!!!\n********\n");
 }
@@ -27,7 +27,7 @@ void output(std::string operation, complex c1, complex c2) {
 	   c1.r, c1.i,
 	   c2.r, c2.i);
     // Epsilon for double is 1e-16, check it 1.0 + 1e-16 = 1.0
-    if (diff_r > 1e-16||diff_i > 1e-16) 
+    if (std::abs(diff_r) > 1e-16||std::abs(diff_i) > 1e-16) 
 	//if (diff_r > 1e-15||diff_i > 1e-15) 
 	printf("\n********\n\t\tWARNING!! Non-zero diff!!!\n********\n");
 }
@@ -43,7 +43,7 @@ void output(std::string operation, std::complex<double> c1, std::complex<double>
 	   c1.real(), c1.imag(),
 	   c2.real(), c2.imag());
     // Epsilon for double is 1e-16, check it 1.0 + 1e-16 = 1.0
-    if (diff_r > 1e-16||diff_i > 1e-16) 
+    if (std::abs(diff_r) > 1e-16||std::abs(diff_i) > 1e-16) 
 	//if (diff_r > 1e-15||diff_i > 1e-15) 
 	printf("\n********\n\t\tWARNING!! Non-zero diff!!!\n********\n");
 }
@@ -56,7 +56,7 @@ void output_double(std::string operation, double c1, double c2) {
 	   operation.c_str(),
 	   diff, c1, c2);
     // Epsilon for double is 1e-16, check it 1.0 + 1e-16 = 1.0
-    if (diff > 1e-16)
+    if (std::abs(diff) > 1e-16)
 	printf("\n********\n\t\tWARNING!! Non-zero diff!!!\n********\n");
 }
 /********************************************************************/