removed warnings about unsigned int comparison

.gitignore

@@ -41,4 +41,4 @@ cov-int
 build*
 .ipynb_checkpoints
 tests/python/field.txt
-
+scattnlay.cpp

Makefile

@@ -44,7 +44,7 @@ builddeb:
 	export CFLAGS='-std=c++11' && dpkg-buildpackage -i -I -rfakeroot
 
 standalone: standalone.cc nmie.cc
-	export CFLAGS='-std=c++11' && c++ -DNDEBUG -O2 -std=c++11 standalone.cc nmie.cc -lm -o scattnlay
+	export CFLAGS='-std=c++11' && c++ -DNDEBUG -O2 -Wall -std=c++11 standalone.cc nmie.cc -lm -o scattnlay
 	mv scattnlay ../
 
 clean:

compare.cc

@@ -38,7 +38,7 @@
 //sudo aptitude install libgoogle-perftools-dev
 #include <google/heap-profiler.h>
 #include "nmie.h"
-#include "nmie-wrapper.h"
+#include "nmie-old.h"
 
 timespec diff(timespec start, timespec end);
 const double PI=3.14159265358979323846;
@@ -223,7 +223,7 @@ int main(int argc, char *argv[]) {
     // do {
     //   clock_gettime(CLOCK_PROCESS_CPUTIME_ID, &time1);
     //   for (int i = 0; i<repeats; ++i) {
-    // 	nmie::nMie_wrapper(L, x, m, nt, Theta, &Qextw, &Qscaw,
+    // 	nmie::nMie(L, x, m, nt, Theta, &Qextw, &Qscaw,
     // 			   &Qabsw, &Qbkw, &Qprw, &gw, &Albedow, S1w, S2w);
     //   }
     //   clock_gettime(CLOCK_PROCESS_CPUTIME_ID, &time2);
@@ -252,15 +252,15 @@ int main(int argc, char *argv[]) {
     // } while (cpptime_nsec < 1e8 && ctime_nsec < 1e8);
 
     nMie(L, x, m, nt, Theta, &Qext, &Qsca, &Qabs, &Qbk, &Qpr, &g, &Albedo, S1, S2);
-    nmie::nMie_wrapper(L, x, m, nt, Theta, &Qextw, &Qscaw, &Qabsw, &Qbkw, &Qprw, &gw, &Albedow, S1w, S2w);
+    nmie::nMie(L, x, m, nt, Theta, &Qextw, &Qscaw, &Qabsw, &Qbkw, &Qprw, &gw, &Albedow, S1w, S2w);
         printf("\n");
     
     if (has_comment) {
       printf("%6s, %+.5e, %+.5e, %+.5e, %+.5e, %+.5e, %+.5e, %+.5e  old\n", comment.c_str(), Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo);
-      printf("%6s, %+.5e, %+.5e, %+.5e, %+.5e, %+.5e, %+.5e, %+.5e  wrapper\n", comment.c_str(), Qextw, Qscaw, Qabsw, Qbkw, Qprw, gw, Albedow);
+      printf("%6s, %+.5e, %+.5e, %+.5e, %+.5e, %+.5e, %+.5e, %+.5e  \n", comment.c_str(), Qextw, Qscaw, Qabsw, Qbkw, Qprw, gw, Albedow);
     } else {
       printf("%+.5e, %+.5e, %+.5e, %+.5e, %+.5e, %+.5e, %+.5e  old\n", Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo);
-      printf("%+.5e, %+.5e, %+.5e, %+.5e, %+.5e, %+.5e, %+.5e  wrapper\n", Qextw, Qscaw, Qabsw, Qbkw, Qprw, gw, Albedow);
+      printf("%+.5e, %+.5e, %+.5e, %+.5e, %+.5e, %+.5e, %+.5e  \n", Qextw, Qscaw, Qabsw, Qbkw, Qprw, gw, Albedow);
     }
     
     if (nt > 0) {
@@ -268,21 +268,21 @@ int main(int argc, char *argv[]) {
       
       for (i = 0; i < nt; i++) {
         printf("%6.2f, %+.5e, %+.5e, %+.5e, %+.5e  old\n", Theta[i]*180.0/PI, S1[i].real(), S1[i].imag(), S2[i].real(), S2[i].imag());
-        printf("%6.2f, %+.5e, %+.5e, %+.5e, %+.5e  wrapper\n", Theta[i]*180.0/PI, S1w[i].real(), S1w[i].imag(), S2w[i].real(), S2w[i].imag());
+        printf("%6.2f, %+.5e, %+.5e, %+.5e, %+.5e  \n", Theta[i]*180.0/PI, S1w[i].real(), S1w[i].imag(), S2w[i].real(), S2w[i].imag());
       }
     }
     // Field testing
     double from_coord = -3.0, to_coord = 3000.0;
     //double size=2.0*PI*1.0/6.0;
-    double size=0.001;
+    double size=0.01;
     std::vector<double> range;
     // for (int i = 0; i < samples; ++i) {
       //range.push_back( from_coord + (to_coord-from_coord)/(static_cast<double>(samples)-1)*i );
     //range.push_back(size*0.01);
     //range.push_back(size*0.99999);
     range.push_back(size/2.0);
-    //  range.push_back(size*1.00001);
-    // range.push_back(3);
+    //range.push_back(size*1.00001);
+     range.push_back(3);
     //printf("r=%g  ", range.back());
     //}
     int samples = range.size();
@@ -302,9 +302,9 @@ int main(int argc, char *argv[]) {
     Zp.insert(Zp.end(), range.begin(), range.end());
     int ncoord = Xp.size();
     x = {size};
-    //m = {std::complex<double>(1.0000002,0.00)};
+    m = {std::complex<double>(2.000000,1.00)};
     // x = {size};
-    m = {std::complex<double>(1.33,0.0)};
+    //m = {std::complex<double>(1.33,0.0)};
     // x = {1.017, 2.016};
     // m = {std::complex<double>(1.5016,1.023), std::complex<double>(2.014,2.012)};
     L = x.size();
@@ -315,23 +315,10 @@ int main(int argc, char *argv[]) {
     for (auto& f:H) f.resize(3);
     double free_impedance = 376.73031;
     //double free_impedance = 1.0;
-    nmax =  nField( L,  pl,  x,  m, nmax,  ncoord,  Xp,  Yp,  Zp, E, H);
-    double sum_e = 0.0, sum_h = 0.0;
-    printf ("Field total sum (old)\n");
-    for (auto f:E) {
-      sum_e = 0.0;
-      for (auto c:f) sum_e+=std::abs(pow2(c));
-      printf("Field E=%g\n", std::sqrt(std::abs(sum_e)));
-    }
-    for (auto f:H) {
-      sum_h = 0.0;
-      for (auto c:f) sum_h+=std::abs(pow2(c));
-      printf("Field H=%g\n", std::sqrt(std::abs(sum_h))*free_impedance);
-    }
 
     nmie::nField( L,  pl,  x,  m, nmax,  ncoord,  Xp,  Yp,  Zp, E, H);
-    sum_e = 0.0, sum_h = 0.0;
-    printf ("Field total sum (wrapper)\n");
+    double sum_e = 0.0, sum_h = 0.0;
+    printf ("Field total sum ()\n");
     for (auto f:E) {
       sum_e = 0.0;
       for (auto c:f) sum_e+=std::abs(pow2(c));

go.sh

@@ -10,7 +10,7 @@ rm -f ../scattnlay
 
 #google profiler  ######## FAST!!!
 echo Uncomment next line to compile compare.cc
-g++ -Ofast -std=c++11 $file nmie.cc nmie-wrapper.cc -lm -lrt -o scattnlay.bin /usr/lib/libtcmalloc.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 nmie.cc nmie-old.cc -lm -lrt -o scattnlay.bin /usr/lib/libtcmalloc.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 compare.cc nmie.cc nmie-wrapper.cc -lm -lrt -o scattnlay-g.bin -ltcmalloc -fno-builtin-malloc -fno-builtin-calloc -fno-builtin-realloc -fno-builtin-free -g
 
@@ -27,8 +27,8 @@ cp scattnlay.bin ../scattnlay
 #     fi
 # done
 
-PROGRAM='../../../scattnlay'
-#time  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  -t 0.0 90.0 5 -c test01
+PROGRAM='./scattnlay.bin'
+time  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  -t 0.0 90.0 5 -c test01
 # #result
 # # test01, +1.41154e+00, +4.17695e-01, +9.93844e-01, +1.59427e-01, +1.25809e+00, +3.67376e-01, +2.95915e-01
 
@@ -93,11 +93,12 @@ PROGRAM='../../../scattnlay'
 
 # ./$file.bin
 
-rm scattnlay.so
-export CFLAGS='-std=c++11'
-python setup.py build_ext --inplace
-cp scattnlay.so tests/python/
-cd tests/python/
-./field.py
-./test01.py
-./test01-wrapper.py
+# ### Cython
+# rm scattnlay.so
+# export CFLAGS='-std=c++11'
+# python setup.py build_ext --inplace
+# cp scattnlay.so tests/python/
+# cd tests/python/
+# ./field.py
+# ./test01.py
+# ./test01-wrapper.py

nmie.cc

@@ -82,7 +82,7 @@ namespace nmie {
   // Return value:                                                                    //
   //   Number of multipolar expansion terms used for the calculations                 //
   //**********************************************************************************//
-  int nMie(const int L, const int pl, std::vector<double>& x, std::vector<std::complex<double> >& m, const 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, 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) {
 
     if (x.size() != L || m.size() != L)
         throw std::invalid_argument("Declared number of layers do not fit x and m!");
@@ -144,7 +144,7 @@ namespace nmie {
   // Return value:                                                                    //
   //   Number of multipolar expansion terms used for the calculations                 //
   //**********************************************************************************//
-  int nMie(const int L, std::vector<double>& x, std::vector<std::complex<double> >& m, const 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) {
+  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) {
     return nmie::nMie(L, -1, x, m, nTheta, Theta, -1, Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo, S1, S2);
   }
 
@@ -176,7 +176,7 @@ namespace nmie {
   // Return value:                                                                    //
   //   Number of multipolar expansion terms used for the calculations                 //
   //**********************************************************************************//
-  int nMie(const int L, const int pl, std::vector<double>& x, std::vector<std::complex<double> >& m, const 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) {
+  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, 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) {
     return nmie::nMie(L, pl, x, m, nTheta, Theta, -1, Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo, S1, S2);
   }
 
@@ -210,7 +210,7 @@ namespace nmie {
   // Return value:                                                                    //
   //   Number of multipolar expansion terms used for the calculations                 //
   //**********************************************************************************//
-  int nMie(const int L, std::vector<double>& x, std::vector<std::complex<double> >& m, const 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, 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) {
     return nmie::nMie(L, -1, x, m, nTheta, Theta, nmax, Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo, S1, S2);
   }
 
@@ -237,7 +237,7 @@ namespace nmie {
   // Return value:                                                                    //
   //   Number of multipolar expansion terms used for the calculations                 //
   //**********************************************************************************//
-  int nField(const int L, const int pl, const std::vector<double>& x, const std::vector<std::complex<double> >& m, const int nmax, const int ncoord, const std::vector<double>& Xp_vec, const std::vector<double>& Yp_vec, const std::vector<double>& Zp_vec, std::vector<std::vector<std::complex<double> > >& E, std::vector<std::vector<std::complex<double> > >& H) {
+  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_vec, const std::vector<double>& Yp_vec, const std::vector<double>& Zp_vec, std::vector<std::vector<std::complex<double> > >& E, std::vector<std::vector<std::complex<double> > >& H) {
     if (x.size() != L || m.size() != L)
       throw std::invalid_argument("Declared number of layers do not fit x and m!");
     if (Xp_vec.size() != ncoord || Yp_vec.size() != ncoord || Zp_vec.size() != ncoord
@@ -420,7 +420,7 @@ namespace nmie {
     isIntCoeffsCalc_ = false;
     isExtCoeffsCalc_ = false;
     isMieCalculated_ = false;
-    if (layer_position < 0)
+    if (layer_position < 0 && layer_position != -1)
       throw std::invalid_argument("Error! Layers are numbered from 0!");
     PEC_layer_position_ = layer_position;
   }
@@ -487,19 +487,19 @@ namespace nmie {
   // ********************************************************************** //
   // Maximum number of terms required for the calculation                   //
   // ********************************************************************** //
-  void MultiLayerMie::calcNmax(int first_layer) {
+  void MultiLayerMie::calcNmax(unsigned int first_layer) {
     int ri, riM1;
     const std::vector<double>& x = size_param_;
     const std::vector<std::complex<double> >& m = refractive_index_;
     calcNstop();  // Set initial nmax_ value
-    for (int i = first_layer; i < x.size(); i++) {
-      if (i > PEC_layer_position_)
+    for (unsigned int i = first_layer; i < x.size(); i++) {
+      if (static_cast<int>(i) > PEC_layer_position_)  // static_cast used to avoid warning
         ri = round(std::abs(x[i]*m[i]));
       else
         ri = 0;
       nmax_ = std::max(nmax_, ri);
       // first layer is pec, if pec is present
-      if ((i > first_layer) && ((i - 1) > PEC_layer_position_))
+      if ((i > first_layer) && (static_cast<int>(i - 1) > PEC_layer_position_))
         riM1 = round(std::abs(x[i - 1]* m[i]));
       else
         riM1 = 0;
@@ -993,7 +993,7 @@ namespace nmie {
       Qbktmp += (double)(n + n + 1)*(1 - 2*(n % 2))*(an_[i]- bn_[i]);
       // Calculate the scattering amplitudes (S1 and S2)    //
       // Equations (25a) - (25b)                            //
-      for (int t = 0; t < theta_.size(); t++) {
+      for (unsigned int t = 0; t < theta_.size(); t++) {
         calcPiTau(std::cos(theta_[t]), Pi, Tau);
 
         S1_[t] += calc_S1(n, an_[i], bn_[i], Pi[i], Tau[i]);
@@ -1313,12 +1313,11 @@ namespace nmie {
       std::vector<std::complex<double> > Es(3), Hs(3);
 
       // Firstly the easiest case: the field outside the particle
-      if (Rho >= GetSizeParameter()) {
-        fieldInt(Rho, Phi, Theta, Es, Hs);
-//        fieldExt(Rho, Phi, Theta, Es, Hs);
-      } else {
-        fieldInt(Rho, Phi, Theta, Es, Hs);
-      }
+      //      if (Rho >= GetSizeParameter()) {
+      //        fieldExt(Rho, Phi, Theta, Es, Hs);
+      // } else {
+      fieldInt(Rho, Phi, Theta, Es, Hs);  //Should work fine both: inside and outside the particle
+      //}
 
       { //Now, convert the fields back to cartesian coordinates
         using std::sin;

nmie.h

@@ -33,12 +33,12 @@
 #include <vector>
 
 namespace nmie {
-  int ScattCoeffs(const 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 int L, const int pl, std::vector<double>& x, std::vector<std::complex<double> >& m, const 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 int L, std::vector<double>& x, std::vector<std::complex<double> >& m, const 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);
-  int nMie(const int L, const int pl, std::vector<double>& x, std::vector<std::complex<double> >& m, const 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);
-  int nMie(const int L, std::vector<double>& x, std::vector<std::complex<double> >& m, const 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 nField(const int L, const int pl, const std::vector<double>& x, const std::vector<std::complex<double> >& m, const int nmax, const 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);
+  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);
+  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, 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, 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 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);
 
   class MultiLayerMie {
    public:
@@ -97,7 +97,7 @@ namespace nmie {
     std::vector<std::vector< std::complex<double> > > GetFieldH(){return H_;};
   private:
     void calcNstop();
-    void calcNmax(int first_layer);
+    void calcNmax(unsigned int first_layer);
 
     std::complex<double> calc_an(int n, double XL, std::complex<double> Ha, std::complex<double> mL,
                                  std::complex<double> PsiXL, std::complex<double> ZetaXL,

standalone.cc

@@ -68,7 +68,8 @@ int main(int argc, char *argv[]) {
     // for (auto arg : args) std::cout<< arg <<std::endl;
     std::string comment;
     int has_comment = 0;
-    int i, l, L = 0;
+    int i;
+    unsigned int L = 0;
     std::vector<double> x, Theta;
     std::vector<std::complex<double> > m, S1, S2;
     double Qext, Qabs, Qsca, Qbk, Qpr, g, Albedo;