seem to work ok with optimizer

go.sh

@@ -86,9 +86,9 @@ echo
 
 cd $path
 file=test-negative-epsilon.cc
-g++ -Ofast -std=c++11 $file nmie.cc -lm -lrt -o $file.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 -lm -lrt -o $file.bin /usr/lib/libtcmalloc.so.4 -fno-builtin-malloc -fno-builtin-calloc -fno-builtin-realloc -fno-builtin-free -march=native -mtune=native -msse4.2
 
 #DEBUG!
-#clang++ -g -O1 -fsanitize=address  -fno-optimize-sibling-calls -fno-omit-frame-pointer -std=c++11 $file nmie.cc -lm -lrt -o $file.bin
+clang++ -g -O1 -fsanitize=address  -fno-optimize-sibling-calls -fno-omit-frame-pointer -std=c++11 $file nmie.cc -lm -lrt -o $file.bin
 
 ./$file.bin

nmie.cc

@@ -262,8 +262,10 @@ namespace nmie {
   // ********************************************************************** //
   // ********************************************************************** //
   // ********************************************************************** //
-  void MultiLayerMie::SetMaxTermsNumber(int nmax) {
+  void MultiLayerMie::SetMaxTermsNumber(int nmax) {    
     nmax_preset_ = nmax;
+    //debug
+    printf("Setting max terms: %d\n", nmax_preset_);
   }
   // ********************************************************************** //
   // ********************************************************************** //

test-negative-epsilon.cc

@@ -32,26 +32,59 @@
 #include <stdio.h>
 template<class T> inline T pow2(const T value) {return value*value;};
 const double PI=3.14159265358979323846;  
+nmie::MultiLayerMie multi_layer_mie_;  
+double lambda_work_ = 3.75; // cm
+double a_ = 0.75*lambda_work_;  // 2.8125 cm - size of PEC core
+double min_index_ = 1e-11;
+// ********************************************************************** //
+// ********************************************************************** //
+// ********************************************************************** //
+double EvaluateScatterOnlyIndex(std::vector<double> input) {
+  double Qsca;
+  std::vector<std::complex<double>> cindex;
+  cindex.clear();
+  double k = min_index_, n=min_index_;
+  for (double epsilon : input) {
+    // sqrt(epsilon) = n + i*k
+    k = min_index_;
+    n=min_index_;
+    if (epsilon > 0.0) n=std::sqrt(epsilon);
+    else k = std::sqrt(-epsilon);
+    if (n < min_index_) n = min_index_;
+    if (k < min_index_) k = min_index_;
+    //printf("eps= %g, n=%g, k=%g\n", epsilon, n, k);
+    cindex.push_back(std::complex<double>(n, k));
+  }
+  multi_layer_mie_.SetCoatingIndex(cindex);
+  try {
+    multi_layer_mie_.RunMieCalculations();
+    Qsca = multi_layer_mie_.GetQsca();
+  } catch( const std::invalid_argument& ia ) {
+    Qsca = 0;
+    printf("#");
+    // Will catch if  multi_layer_mie_ fails or other errors.
+    //std::cerr << "Invalid argument: " << ia.what() << std::endl;
+  }  
+  double total_r = multi_layer_mie_.GetTotalRadius();
+  return Qsca*PI*pow2(total_r);
+}
+
 int main(int argc, char *argv[]) {
   try {
-    double lambda_work_ = 3.75; // cm
-    double a_ = 0.75*lambda_work_;  // 2.8125 cm - size of PEC core
-    double min_index_ = 1e-11;
-    nmie::MultiLayerMie multi_layer_mie;  
     // Only PEC target
-    multi_layer_mie.SetTargetPEC(a_);
-    multi_layer_mie.SetWavelength(lambda_work_);
-    multi_layer_mie.RunMieCalculations();
-    double PEC_Qsca = multi_layer_mie.GetQsca();
-    double PEC_r = multi_layer_mie.GetTotalRadius();
+    multi_layer_mie_.SetTargetPEC(a_);
+    multi_layer_mie_.SetWavelength(lambda_work_);
+    multi_layer_mie_.RunMieCalculations();
+    double PEC_Qsca = multi_layer_mie_.GetQsca();
+    double PEC_r = multi_layer_mie_.GetTotalRadius();
     double PEC_RCS = PEC_Qsca*PI*pow2(PEC_r);
     
     // PEC target covered with with air layer
-    multi_layer_mie.SetCoatingWidth({0.1});
-    multi_layer_mie.SetCoatingIndex({{1.0,0.0}});
-    multi_layer_mie.RunMieCalculations();
-    double Qsca1 = multi_layer_mie.GetQsca();
-    double total_r1 = multi_layer_mie.GetTotalRadius();
+    multi_layer_mie_.SetCoatingWidth({0.1});
+    multi_layer_mie_.SetCoatingIndex({{1.0,0.0}});
+    multi_layer_mie_.RunMieCalculations();
+    double Qsca1 = multi_layer_mie_.GetQsca();
+    double total_r1 = multi_layer_mie_.GetTotalRadius();
     double initial_RCS1 = Qsca1*PI*pow2(total_r1);
     printf("RCS = %g cm^2 with (r=%g) and  RCS=%g cm^2 without (r=%g)air coating.\n",
 	   initial_RCS1, total_r1, 
@@ -69,24 +102,37 @@ int main(int argc, char *argv[]) {
     k=min_index_;
     cindex.push_back(std::complex<double>(n, k));
 
-    multi_layer_mie.SetCoatingWidth({0.1,0.1});
-    multi_layer_mie.SetCoatingIndex(cindex);
-    multi_layer_mie.RunMieCalculations();
-    double Qsca = multi_layer_mie.GetQsca();
-    double total_r = multi_layer_mie.GetTotalRadius();
+    multi_layer_mie_.SetCoatingWidth({0.1,0.1});
+    multi_layer_mie_.SetCoatingIndex(cindex);
+    multi_layer_mie_.RunMieCalculations();
+    double Qsca = multi_layer_mie_.GetQsca();
+    double total_r = multi_layer_mie_.GetTotalRadius();
     double initial_RCS = Qsca*PI*pow2(total_r);
     printf("RCS=%g for bi-layer coating (total R=%g).\n", initial_RCS,total_r);
 
     n=1.0;
     k=min_index_;
     cindex.push_back(std::complex<double>(n, k));
-    multi_layer_mie.SetCoatingWidth({0.1,0.1,0.1});
-    multi_layer_mie.SetCoatingIndex(cindex);
-    multi_layer_mie.RunMieCalculations();
-    Qsca = multi_layer_mie.GetQsca();
-    total_r = multi_layer_mie.GetTotalRadius();
+    multi_layer_mie_.SetCoatingWidth({0.1,0.1,0.1});
+    multi_layer_mie_.SetCoatingIndex(cindex);
+    multi_layer_mie_.RunMieCalculations();
+    Qsca = multi_layer_mie_.GetQsca();
+    total_r = multi_layer_mie_.GetTotalRadius();
     initial_RCS = Qsca*PI*pow2(total_r);
     printf("RCS=%g for bi-layer+air coating (total R=%g).\n", initial_RCS,total_r);
+
+
+    //multi_layer_mie_.SetMaxTermsNumber(15);
+    multi_layer_mie_.SetCoatingWidth({0.1,0.1});
+    printf("With %g  coating= (26.24).\n",
+	   EvaluateScatterOnlyIndex({-0.29, 24.6}));
+    multi_layer_mie_.SetCoatingWidth({0.1,0.1,0.1});
+    printf("With %g  coating> (26.24).\n",
+	   EvaluateScatterOnlyIndex({-0.29, 24.6, 1.0}));
+    //multi_layer_mie_.SetMaxTermsNumber(-1);
+
+      // 26.24: 25||   -0.29   +24.60 
+      // 28.48: 38||   -0.29   +24.60    +1.00
     
   } catch( const std::invalid_argument& ia ) {
     // Will catch if  multi_layer_mie fails or other errors.