replace of the reccurence with CONFRA call. Result is the same

compare.cc

@@ -244,7 +244,7 @@ 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);

go.sh

@@ -23,7 +23,9 @@ 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

nmie-wrapper.cc

@@ -545,9 +545,9 @@ c    MAXIT    Max. allowed no. of iterations
 
 c    MM       + 1  and - 1, alternately
 */
-  std::complex<double> MultiLayerMie::calcD1confra(const std::complex<double> z) {
+  std::complex<double> MultiLayerMie::calcD1confra(const int N, const std::complex<double> z) {
   // NTMR -> nmax_ -1  \\TODO nmax_ ?
-    int N = nmax_ - 1;
+    //int N = nmax_ - 1;
     int KK, KOUNT, MAXIT = 10000, MM;
     double EPS1=1.0e-2, EPS2=1.0e-8;
     std::complex<double> CAK, CAPT, CDENOM, CDTD, CNTN, CNUMER;
@@ -585,10 +585,10 @@ c    MM       + 1  and - 1, alternately
      // } else { //c                           *** Well-conditioned case
       {
 	CAPT   = CNUMER / CDENOM; // ** Eq. R27
-       printf("re(%g):im(%g)**\t", CAPT.real(), CAPT.imag());
+	// printf("re(%g):im(%g)**\t", CAPT.real(), CAPT.imag());
        CONFRA = CAPT * CONFRA; // ** Eq. R26
-       //       printf("re(%g):im(%g)\t", CONFRA.real(), CONFRA.imag());
-//c                                  ** Check for convergence; Eq. R31
+       //       if (N == 0) {output=true;printf(" re:");prn(CONFRA.real());printf(" im:"); prn(CONFRA.imag());output=false;};
+       //c                                  ** Check for convergence; Eq. R31
        if ( std::abs(CAPT.real() - 1.0) >= EPS2 ||  std::abs(CAPT.imag()) >= EPS2 ) {
 //c                                        ** Eq. R30
 	 CNUMER = CAK + one/CNUMER;
@@ -599,7 +599,7 @@ c    MM       + 1  and - 1, alternately
       }
       break;
     } while(1);    
-    printf(" return confra\n");
+    //printf(" return confra\n");
     return CONFRA;
   }
   //**********************************************************************************//
@@ -619,17 +619,18 @@ c    MM       + 1  and - 1, alternately
 			       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(z);
+    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());
+    // 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);
-      printf(" D:");prn((D1[n-1]).real()); printf("\t diff:");
-      prn((D1[n] + double(n)*zinv).real());
+      //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;
+    // 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);

nmie-wrapper.h

@@ -64,8 +64,10 @@ namespace nmie {
     // SP stands for size parameter units.
    public:
     long iformat = 0;
+    bool output = false;
     void prn(double var) {
       do {
+	if (!output) break;
 	++iformat;
 	printf("%23.13e",var);	     
 	if (iformat%4 == 0) printf("\n");
@@ -186,7 +188,7 @@ namespace nmie {
 		     std::vector<std::complex<double> > D3,
 		     std::vector<std::complex<double> >& Psi,
 		     std::vector<std::complex<double> >& Zeta);
-    std::complex<double> calcD1confra(const std::complex<double> z);
+    std::complex<double> calcD1confra(int N, const std::complex<double> z);
     void calcD1D3(std::complex<double> z,
 		  std::vector<std::complex<double> >& D1,
 		  std::vector<std::complex<double> >& D3);