instability example

compare.cc

@@ -228,7 +228,7 @@ int main(int argc, char *argv[]) {
       // if (ctime_nsec < best_c) best_c = ctime_nsec;
       
       // clock_gettime(CLOCK_PROCESS_CPUTIME_ID, &time1);
-       nmie::nMie_wrapper(L, x, m, nt, Theta, &Qextw, &Qscaw, &Qabsw, &Qbkw, &Qprw, &gw, &Albedow, S1w, S2w);
+    //nmie::nMie_wrapper(L, x, m, nt, Theta, &Qextw, &Qscaw, &Qabsw, &Qbkw, &Qprw, &gw, &Albedow, S1w, S2w);
       // clock_gettime(CLOCK_PROCESS_CPUTIME_ID, &time2);
       // //cpptime_nsec = std::min(cpptime_nsec, diff(time1,time2).tv_nsec);
       // cpptime_nsec = diff(time1,time2).tv_nsec;
@@ -244,13 +244,14 @@ int main(int argc, char *argv[]) {
     //printf("--best C/C++ time ratio: %Lg\n", static_cast<long double>(best_c)/static_cast<long double>(best_cpp));
 
     // clock_gettime(CLOCK_PROCESS_CPUTIME_ID, &time1);
-       nMie(L, x, m, nt, Theta, &Qext, &Qsca, &Qabs, &Qbk, &Qpr, &g, &Albedo, S1, S2);
+       // nMie(L, x, m, nt, Theta, &Qext, &Qsca, &Qabs, &Qbk, &Qpr, &g, &Albedo, S1, S2);
     // clock_gettime(CLOCK_PROCESS_CPUTIME_ID, &time2);
     //  ctime_nsec = diff(time1,time2).tv_nsec;
     // printf("-- C time consumed %ld sec : %ld nsec\n",diff(time1,time2).tv_sec, ctime_nsec);
 
     // clock_gettime(CLOCK_PROCESS_CPUTIME_ID, &time1);
        nmie::nMie_wrapper(L, x, m, nt, Theta, &Qextw, &Qscaw, &Qabsw, &Qbkw, &Qprw, &gw, &Albedow, S1w, S2w);
+        printf("\n");
     // clock_gettime(CLOCK_PROCESS_CPUTIME_ID, &time2);
     // cpptime_nsec = diff(time1,time2).tv_nsec;
     // printf("-- C++ time consumed %ld sec : %ld nsec\n",diff(time1,time2).tv_sec,cpptime_nsec);
@@ -261,7 +262,7 @@ int main(int argc, char *argv[]) {
       printf("%6s, %+.5e, %+.5e, %+.5e, %+.5e, %+.5e, %+.5e, %+.5e\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);
     } else {
-      printf("%+.5e, %+.5e, %+.5e, %+.5e, %+.5e, %+.5e, %+.5e\n", Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo);
+      //printf("%+.5e, %+.5e, %+.5e, %+.5e, %+.5e, %+.5e, %+.5e\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);
     }
     

go.sh

@@ -25,15 +25,15 @@ cd tests/shell
 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  -t 0.0 90.0 5 -c test01
 
-echo BUG -- All 3 designs shoul give almost the same answer
-echo
-echo $PROGRAM -l 1 4.71238898038469 2 0.0001
+echo BUG -- All 3 designs should give almost the same answer
+#echo
+#echo $PROGRAM -l 1 4.71238898038469 2 0.0001
 $PROGRAM -l 1 4.71238898038469 2 0.0001
-echo
-echo $PROGRAM -l 2 4.71238898038469 2 0.0001 9.42477796076937 1 0
+#echo
+#echo $PROGRAM -l 2 4.71238898038469 2 0.0001 9.42477796076937 1 0
 $PROGRAM -l 2 4.71238898038469 2 0.0001 9.42477796076937 1 0
-echo
-echo $PROGRAM -l 2 4.71238898038469 2 0.0001 9.42477796076938 1 0
+#echo
+#echo $PROGRAM -l 2 4.71238898038469 2 0.0001 9.42477796076938 1 0
 $PROGRAM -l 2 4.71238898038469 2 0.0001 9.42477796076938 1 0
 echo
 #  #apt-get install oprofile

nmie-wrapper.cc

@@ -520,7 +520,10 @@ namespace nmie {
     const std::complex<double> zinv = std::complex<double>(1.0, 0.0)/z;
     for (int n = nmax_; n > 0; n--) {
       D1[n - 1] = double(n)*zinv - 1.0/(D1[n] + double(n)*zinv);
+      printf(" D:");prn((D1[n-1]).real()); printf("\t diff:");
+      prn((D1[n] + double(n)*zinv).real());
     }
+    printf("\n\n"); iformat=0;
     // Upward recurrence for PsiZeta and D3 - equations (18a) - (18d)
     PsiZeta_[0] = 0.5*(1.0 - std::complex<double>(cos(2.0*z.real()), sin(2.0*z.real()))*exp(-2.0*z.imag()));
     D3[0] = std::complex<double>(0.0, 1.0);
@@ -641,6 +644,10 @@ namespace nmie {
       // Calculate D1 and D3
       calcD1D3(z1, D1_mlxl, D3_mlxl);
     }
+    // do { \
+    //   ++iformat;\
+    //   if (iformat%5 == 0) printf("%24.16e",z1.real());	\
+    // } while (false);
     //******************************************************************//
     // Calculate Ha and Hb in the first layer - equations (7a) and (8a) //
     //******************************************************************//
@@ -663,6 +670,11 @@ namespace nmie {
       calcD1D3(z1, D1_mlxl, D3_mlxl);
       //Calculate D1 and D3 for z2
       calcD1D3(z2, D1_mlxlM1, D3_mlxlM1);
+      // prn(z1.real());
+      // for ( auto i : D1_mlxl) { prn(i.real());
+      //   // prn(i.imag());
+      // 	} printf("\n");
+
       //*********************************************//
       //Calculate Q, Ha and Hb in the layers fl+1..L //
       //*********************************************//
@@ -710,6 +722,7 @@ namespace nmie {
     //**************************************//
     // Calculate D1XL and D3XL
     calcD1D3(x[L - 1],  D1XL, D3XL);
+    //printf("%5.20f\n",Ha[L-1][0].real());
     // Calculate PsiXL and ZetaXL
     calcPsiZeta(x[L - 1], D1XL, D3XL, PsiXL, ZetaXL);
     //*********************************************************************//

nmie-wrapper.h

@@ -51,7 +51,9 @@
 #   define ASSERT(condition, message) do { } while (false)
 #endif
 
+
 namespace nmie {
+
   int nMie_wrapper(int L, const std::vector<double>& x, const std::vector<std::complex<double> >& m,
          int nTheta, const std::vector<double>& Theta,
          double *Qext, double *Qsca, double *Qabs, double *Qbk, double *Qpr, double *g, double *Albedo,
@@ -61,6 +63,14 @@ namespace nmie {
     // Will throw for any error!
     // SP stands for size parameter units.
    public:
+    long iformat = 0;
+    void prn(double var) {
+      do {
+	++iformat;
+	printf("%23.13e",var);	     
+	if (iformat%4 == 0) printf("\n");
+      } while (false);
+    }
     // Set parameters in applied units 
     void SetWavelength(double wavelength) {wavelength_ = wavelength;};
     // It is possible to set only a multilayer target to run calculaitons.