an, bn and an_bulk, bn_bulk looks to be the same

compare.cc

@@ -213,26 +213,26 @@ int main(int argc, char *argv[]) {
     //   repeats *= 10;
     // } while (cpptime_nsec < 1e8 && ctime_nsec < 1e8);
 
-    // nMie(L, x, m, nt, Theta, &Qext, &Qsca, &Qabs, &Qbk, &Qpr, &g, &Albedo, S1, S2);
-    // nmie::nMie(L, x, m, nt, Theta, &Qextw, &Qscaw, &Qabsw, &Qbkw, &Qprw, &gw, &Albedow, S1w, S2w);
-    //     printf("\n");
+    nMie(L, x, m, nt, Theta, &Qext, &Qsca, &Qabs, &Qbk, &Qpr, &g, &Albedo, S1, S2);
+    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  \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  \n", Qextw, Qscaw, Qabsw, Qbkw, Qprw, gw, Albedow);
-    // }
+    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  \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  \n", Qextw, Qscaw, Qabsw, Qbkw, Qprw, gw, Albedow);
+    }
     
-    // if (nt > 0) {
-    //   printf(" Theta,         S1.r,         S1.i,         S2.r,         S2.i\n");
+    if (nt > 0) {
+      printf(" Theta,         S1.r,         S1.i,         S2.r,         S2.i\n");
       
-    //   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  \n", Theta[i]*180.0/PI, S1w[i].real(), S1w[i].imag(), S2w[i].real(), S2w[i].imag());
-    //   }
-    // }
+      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  \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;

nmie.cc

@@ -629,32 +629,36 @@ namespace nmie {
                                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);
-    const std::complex<double> zinv = std::complex<double>(1.0, 0.0)/z;
-
-    for (int n = nmax_; n > 0; n--) {
-      D1[n - 1] = static_cast<double>(n)*zinv - 1.0/(D1[n] + static_cast<double>(n)*zinv);
-    }
-
-    if (std::abs(D1[0]) > 100000.0)
-      throw std::invalid_argument("Unstable D1! Please, try to change input parameters!\n");
-
-    // Upward recurrence for PsiZeta and D3 - equations (18a) - (18d)
-    PsiZeta_[0] = 0.5*(1.0 - std::complex<double>(std::cos(2.0*z.real()), std::sin(2.0*z.real()))
-                       *std::exp(-2.0*z.imag()));
-    D3[0] = std::complex<double>(0.0, 1.0);
-    for (int n = 1; n <= nmax_; n++) {
-      PsiZeta_[n] = PsiZeta_[n - 1]*(static_cast<double>(n)*zinv - D1[n - 1])
-                                   *(static_cast<double>(n)*zinv - D3[n - 1]);
-      D3[n] = D1[n] + std::complex<double>(0.0, 1.0)/PsiZeta_[n];
-    }
+    // // Downward recurrence for D1 - equations (16a) and (16b)
+    // D1[nmax_] = std::complex<double>(0.0, 0.0);
+    // const std::complex<double> zinv = std::complex<double>(1.0, 0.0)/z;
+
+    // for (int n = nmax_; n > 0; n--) {
+    //   D1[n - 1] = static_cast<double>(n)*zinv - 1.0/(D1[n] + static_cast<double>(n)*zinv);
+    // }
+
+    // if (std::abs(D1[0]) > 100000.0)
+    //   throw std::invalid_argument("Unstable D1! Please, try to change input parameters!\n");
+
+    // // Upward recurrence for PsiZeta and D3 - equations (18a) - (18d)
+    // PsiZeta_[0] = 0.5*(1.0 - std::complex<double>(std::cos(2.0*z.real()), std::sin(2.0*z.real()))
+    //                    *std::exp(-2.0*z.imag()));
+    // D3[0] = std::complex<double>(0.0, 1.0);
+    // for (int n = 1; n <= nmax_; n++) {
+    //   PsiZeta_[n] = PsiZeta_[n - 1]*(static_cast<double>(n)*zinv - D1[n - 1])
+    //                                *(static_cast<double>(n)*zinv - D3[n - 1]);
+    //   D3[n] = D1[n] + std::complex<double>(0.0, 1.0)/PsiZeta_[n];
+    // }
     std::vector<std::complex<double> >  ZetaZ(nmax_+1), dZetaZ(nmax_ + 1);
     bessel::calcZeta(nmax_, z, ZetaZ, dZetaZ );
     for (int n = 0; n < nmax_+1; ++n) {
       D3[n]=dZetaZ[n]/ZetaZ[n];
     }
-
+    std::vector<std::complex<double> >  PsiZ(nmax_+1), dPsiZ(nmax_ + 1);
+    bessel::calcPsi(nmax_, z, PsiZ, dPsiZ );
+    for (int n = 0; n < nmax_+1; ++n) {
+      D1[n]=dPsiZ[n]/PsiZ[n];
+    }
   }
 
 
@@ -674,24 +678,24 @@ namespace nmie {
                                   std::vector<std::complex<double> >& Psi,
                                   std::vector<std::complex<double> >& Zeta) {
 
-    std::complex<double> c_i(0.0, 1.0);
-    std::vector<std::complex<double> > D1(nmax_ + 1), D3(nmax_ + 1);
+    // std::complex<double> c_i(0.0, 1.0);
+    // std::vector<std::complex<double> > D1(nmax_ + 1), D3(nmax_ + 1);
 
-    // First, calculate the logarithmic derivatives
-    calcD1D3(z, D1, D3);
+    // // First, calculate the logarithmic derivatives
+    // calcD1D3(z, D1, D3);
 
-    // Now, use the upward recurrence to calculate Psi and Zeta - equations (20a) - (21b)
-    Psi[0] = std::sin(z);
-    Zeta[0] = std::sin(z) - c_i*std::cos(z);
-    for (int n = 1; n <= nmax_; n++) {
-      Psi[n]  =  Psi[n - 1]*(static_cast<double>(n)/z - D1[n - 1]);
-      Zeta[n] = Zeta[n - 1]*(static_cast<double>(n)/z - D3[n - 1]);
-    }
+    // // Now, use the upward recurrence to calculate Psi and Zeta - equations (20a) - (21b)
+    // Psi[0] = std::sin(z);
+    // Zeta[0] = std::sin(z) - c_i*std::cos(z);
+    // for (int n = 1; n <= nmax_; n++) {
+    //   Psi[n]  =  Psi[n - 1]*(static_cast<double>(n)/z - D1[n - 1]);
+    //   Zeta[n] = Zeta[n - 1]*(static_cast<double>(n)/z - D3[n - 1]);
+    // }
     
     std::vector<std::complex<double> >  dZetaZ(nmax_ + 1);
     bessel::calcZeta(nmax_, z, Zeta, dZetaZ);
-
-
+    std::vector<std::complex<double> > dPsiZ(nmax_ + 1);
+    bessel::calcPsi(nmax_, z, Psi, dPsiZ );
   }
 
 
@@ -739,6 +743,12 @@ namespace nmie {
         h1np[n] = h1n[n - 1] - static_cast<double>(n + 1)*h1n[n]/z;
       }
     }
+    // std::vector< std::complex<double> > yn, ynp;
+    // int nm;
+    // csphjy (nmax_, z, nm, jn, jnp,  yn, ynp );
+    // auto c_i = std::complex<double>(0.0,1.0);
+
+
   }