buggy core-shell design

compare.cc

@@ -262,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

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

nmie-wrapper.cc

@@ -618,19 +618,20 @@ c    MM       + 1  and - 1, alternately
 			       std::vector<std::complex<double> >& D1,
 			       std::vector<std::complex<double> >& D3) {
     // Downward recurrence for D1 - equations (16a) and (16b)
-    //D1[nmax_] = std::complex<double>(0.0, 0.0);
-    D1[nmax_] = calcD1confra(nmax_, z);
+    D1[nmax_] = std::complex<double>(0.0, 0.0);
+    //D1[nmax_] = calcD1confra(nmax_, z);
     const std::complex<double> zinv = std::complex<double>(1.0, 0.0)/z;
     
     // printf(" D:");prn((D1[nmax_]).real()); printf("\t diff:");
     // prn((D1[nmax_] + double(nmax_)*zinv).real());
     for (int n = nmax_; n > 0; n--) {
-      //D1[n - 1] = double(n)*zinv - 1.0/(D1[n] + double(n)*zinv);
-      D1[n-1] = calcD1confra(n-1, z);
+      D1[n - 1] = double(n)*zinv - 1.0/(D1[n] + double(n)*zinv);
+      //D1[n-1] = calcD1confra(n-1, z);
       // printf(" D:");prn((D1[n-1]).real()); printf("\t diff:");
       // prn((D1[n] + double(n)*zinv).real());
     }
     // printf("\n\n"); iformat=0;
+    //if (D1[0].real() > 100.0 ) throw std::invalid_argument("Unstable D1!\n");
     // 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);

nmie-wrapper.h

@@ -64,7 +64,7 @@ namespace nmie {
     // SP stands for size parameter units.
    public:
     long iformat = 0;
-    bool output = false;
+    bool output = true;
     void prn(double var) {
       do {
 	if (!output) break;