Backward order calculation of expansion coeffs

compare.cc

@@ -252,18 +252,18 @@ int main(int argc, char *argv[]) {
     int samples = range.size();
     std::vector<double> zero(samples, 0.0);
     std::vector<double> Xp, Yp, Zp;
-    // // X line
-    // Xp.insert(Xp.end(), range.begin(), range.end());
-    // Yp.insert(Yp.end(), zero.begin(), zero.end());
-    // Zp.insert(Zp.end(), zero.begin(), zero.end());
-    // Y line
+     // X line
+    Xp.insert(Xp.end(), range.begin(), range.end());
+    Yp.insert(Yp.end(), zero.begin(), zero.end());
+    Zp.insert(Zp.end(), zero.begin(), zero.end());
+    //Y line
     Xp.insert(Xp.end(), zero.begin(), zero.end());
     Yp.insert(Yp.end(), range.begin(), range.end());
     Zp.insert(Zp.end(), zero.begin(), zero.end());
     // Z line
-    // Xp.insert(Xp.end(), zero.begin(), zero.end());
-    // Yp.insert(Yp.end(), zero.begin(), zero.end());
-    // Zp.insert(Zp.end(), range.begin(), range.end());
+    Xp.insert(Xp.end(), zero.begin(), zero.end());
+    Yp.insert(Yp.end(), zero.begin(), zero.end());
+    Zp.insert(Zp.end(), range.begin(), range.end());
     int ncoord = Xp.size();
     // Test solid sphere
     x = {size};

go.sh

@@ -99,6 +99,7 @@ time  ASAN_SYMBOLIZER_PATH=/usr/bin/llvm-symbolizer-3.4  $PROGRAM -l 5 0.4642 1.
 # python setup.py build_ext --inplace
 # cp scattnlay.so tests/python/
 # cd tests/python/
+# ./field-dielectric-sphere..py
 # ./field.py
 # ./test01.py
 # ./test01-wrapper.py

nmie.cc

@@ -1243,23 +1243,17 @@ namespace nmie {
 
     // Yang, paragraph under eq. A3
     // a^(L + 1)_n = a_n, d^(L + 1) = 1 ...
-    for (int n = 0; n < nmax_; n++) {
-      aln_[L][n] = an_[n];
-      bln_[L][n] = bn_[n];
-      cln_[L][n] = c_one;
-      dln_[L][n] = c_one;
-
-      printf("aln_[%02i, %02i] = %g,%g; bln_[%02i, %02i] = %g,%g; cln_[%02i, %02i] = %g,%g; dln_[%02i, %02i] = %g,%g\n", L, n, std::real(aln_[L][n]), std::imag(aln_[L][n]), L, n, std::real(bln_[L][n]), std::imag(bln_[L][n]), L, n, std::real(cln_[L][n]), std::imag(cln_[L][n]), L, n, real(dln_[L][n]), std::imag(dln_[L][n]));
+    for (int i = 0; i < nmax_; ++i) {
+      aln_[L][i] = an_[i];
+      bln_[L][i] = bn_[i];
+      cln_[L][i] = c_one;
+      dln_[L][i] = c_one;
     }
 
     std::vector<std::complex<double> > D1z(nmax_ + 1), D1z1(nmax_ + 1), D3z(nmax_ + 1), D3z1(nmax_ + 1);
     std::vector<std::complex<double> > Psiz(nmax_ + 1), Psiz1(nmax_ + 1), Zetaz(nmax_ + 1), Zetaz1(nmax_ + 1);
     std::complex<double> denomZeta, denomPsi, T1, T2, T3, T4;
 
-    std::vector<std::vector<std::complex<double> > > a(2);
-    a[0].resize(3);
-    a[1].resize(3);
-
     auto& m = refractive_index_;
     std::vector< std::complex<double> > m1(L);
 
@@ -1279,53 +1273,33 @@ namespace nmie {
       for (int n = 0; n < nmax_; n++) {
         int n1 = n + 1;
 
-        a[0][0] = m1[l]*D3z[n1]*Zetaz[n1];
-        a[0][1] = -m1[l]*D1z[n1]*Psiz[n1];
-        a[0][2] = aln_[l + 1][n]*m[l]*D3z1[n1]*Zetaz1[n1];
-        a[0][2] -= dln_[l + 1][n]*m[l]*D1z1[n1]*Psiz1[n1];
-
-        a[1][0] = Zetaz[n1];
-        a[1][1] = -Psiz[n1];
-        a[1][2] = aln_[l + 1][n]*Zetaz1[n1] - dln_[l + 1][n]*Psiz1[n1];
-
-        // aln
-        aln_[l][n] = (a[0][2]*a[1][1] - a[0][1]*a[1][2])/(a[0][0]*a[1][1] - a[0][1]*a[1][0]);
-        // dln
-        dln_[l][n] = (a[0][2]*a[1][0] - a[0][0]*a[1][2])/(a[0][1]*a[1][0] - a[0][0]*a[0][1]);
-
-        /*for (int i = 0; i < 2; i++) {
-          for (int j = 0; j < 3; j++) {
-            printf("a[%i, %i] = %g,%g ", i, j, real(a[i][j]), imag(a[i][j]));
-          }
-          printf("\n");
-        }
-        printf("aln_[%i, %i] = %g,%g; dln_[%i, %i] = %g,%g\n\n", l, n, real(aln_[l][n]), imag(aln_[l][n]), l, n, real(dln_[l][n]), imag(dln_[l][n]));*/
+        denomZeta = m1[l]*Zetaz[n1]*(D1z[n1] - D3z[n1]);
+        denomPsi  =  m1[l]*Psiz[n1]*(D1z[n1] - D3z[n1]);
 
-        a[0][0] = D3z[n1]*Zetaz[n1];
-        a[0][1] = -D1z[n1]*Psiz[n1];
-        a[0][2] = bln_[l + 1][n]*D3z1[n1]*Zetaz1[n1];
-        a[0][2] -= cln_[l + 1][n]*D1z1[n1]*Psiz1[n1];
+        T1 = aln_[l + 1][n]*Zetaz1[n1] - dln_[l + 1][n]*Psiz1[n1];
+        T2 = bln_[l + 1][n]*Zetaz1[n1] - cln_[l + 1][n]*Psiz1[n1];
 
-        a[1][0] = m1[l]*Zetaz[n1];
-        a[1][1] = -m1[l]*Psiz[n1];
-        a[1][2] = bln_[l + 1][n]*m[l]*Zetaz1[n1] - cln_[l + 1][n]*m[l]*Psiz1[n1];
+        T3 = D1z1[n1]*dln_[l + 1][n]*Psiz1[n1] - D3z1[n1]*aln_[l + 1][n]*Zetaz1[n1];
+        T4 = D1z1[n1]*cln_[l + 1][n]*Psiz1[n1] - D3z1[n1]*bln_[l + 1][n]*Zetaz1[n1];
 
+        // aln
+        aln_[l][n] = (D1z[n1]*m1[l]*T1 + m[l]*T3)/denomZeta;
         // bln
-        bln_[l][n] = (a[0][2]*a[1][1] - a[0][1]*a[1][2])/(a[0][0]*a[1][1] - a[0][1]*a[1][0]);
+        bln_[l][n] = (D1z[n1]*m[l]*T2 + m1[l]*T4)/denomZeta;
         // cln
-        cln_[l][n] = (a[0][2]*a[1][0] - a[0][0]*a[1][2])/(a[0][1]*a[1][0] - a[0][0]*a[0][1]);
-
-        printf("aln_[%02i, %02i] = %g,%g; bln_[%02i, %02i] = %g,%g; cln_[%02i, %02i] = %g,%g; dln_[%02i, %02i] = %g,%g\n", l, n, real(aln_[l][n]), imag(aln_[l][n]), l, n, real(bln_[l][n]), imag(bln_[l][n]), l, n, real(cln_[l][n]), imag(cln_[l][n]), l, n, real(dln_[l][n]), imag(dln_[l][n]));
+        cln_[l][n] = (D3z[n1]*m[l]*T2 + m1[l]*T4)/denomPsi;
+        // dln
+        dln_[l][n] = (D3z[n1]*m1[l]*T1 + m[l]*T3)/denomPsi;
       }  // end of all n
     }  // end of all l
 
     // Check the result and change  aln_[0][n] and aln_[0][n] for exact zero
     for (int n = 0; n < nmax_; ++n) {
-//      printf("n=%d, aln_=%g,%g,   bln_=%g,%g \n", n, real(aln_[0][n]), imag(aln_[0][n]),
-//	     real(bln_[0][n]), imag(bln_[0][n]));
-      if (std::abs(aln_[0][n]) < 1e-1) aln_[0][n] = 0.0;
+      printf("n=%d, aln_=%g,%g,   bln_=%g,%g \n", n, real(aln_[0][n]), imag(aln_[0][n]),
+	     real(bln_[0][n]), imag(bln_[0][n]));
+      if (std::abs(aln_[0][n]) < 1e-10) aln_[0][n] = 0.0;
       else throw std::invalid_argument("Unstable calculation of aln_[0][n]!");
-      if (std::abs(bln_[0][n]) < 1e-1) bln_[0][n] = 0.0;
+      if (std::abs(bln_[0][n]) < 1e-10) bln_[0][n] = 0.0;
       else throw std::invalid_argument("Unstable calculation of bln_[0][n]!");
     }
 
@@ -1517,7 +1491,7 @@ namespace nmie {
     for (int i = 0; i < 3; i++) {
       H[i] = hffact*H[i];
       Hi[i] *= hffact;
-      printf("E[%i] = %10.5er%+10.5ei; Ei[%i] = %10.5er%+10.5ei; H[%i] = %10.5er%+10.5ei; Hi[%i] = %10.5er%+10.5ei\n", i, std::real(E[i]), std::imag(E[i]), i, std::real(Ei[i]), std::imag(Ei[i]), i, std::real(H[i]), std::imag(H[i]), i, std::real(Hi[i]), std::imag(Hi[i]));
+      //      printf("E[%i] = %10.5er%+10.5ei; Ei[%i] = %10.5er%+10.5ei; H[%i] = %10.5er%+10.5ei; Hi[%i] = %10.5er%+10.5ei\n", i, std::real(E[i]), std::imag(E[i]), i, std::real(Ei[i]), std::imag(Ei[i]), i, std::real(H[i]), std::imag(H[i]), i, std::real(Hi[i]), std::imag(Hi[i]));
     }
    }  // end of MultiLayerMie::calcField(...)
 
@@ -1557,8 +1531,12 @@ namespace nmie {
     //   printf("bn_[%i] = %11.4er%+10.5ei;  bn_bulk_[%i] = %11.4er%+10.5ei\n", n, std::real(bn_[n]), std::imag(bn_[n]), n, std::real(bn1[n]), std::imag(bn1[n]));
     // }
 
+
     // Calculate expansion coefficients aln_,  bln_, cln_, and dln_
     ExpanCoeffs();
+    //ExpanCoeffsV2();
+
+
     // for (int i = 0; i < nmax_; ++i) {
     //   printf("cln_[%i] = %11.4er%+10.5ei;  dln_[%i] = %11.4er%+10.5ei\n", i, std::real(cln_[0][i]), std::imag(cln_[0][i]),
     // 	     i, std::real(dln_[0][i]), std::imag(dln_[0][i]));

tests/python/field-dielectric-sphere.py

@@ -76,8 +76,8 @@ try:
     from matplotlib import cm
     from matplotlib.colors import LogNorm
 
-    min_tick = 0.1
-    max_tick = 1.0
+    min_tick = 0.16
+    max_tick = 0.18
 
     edata = np.resize(Eh, (npts, npts))
 
@@ -88,8 +88,8 @@ try:
     scale_y = np.linspace(min(coordY), max(coordY), npts)
 
     # Define scale ticks
-    min_tick = max(0.1, min(min_tick, np.amin(edata)))
-    max_tick = max(max_tick, np.amax(edata))
+    # min_tick = max(0.1, min(min_tick, np.amin(edata)))
+    # max_tick = max(max_tick, np.amax(edata))
     #scale_ticks = np.power(10.0, np.linspace(np.log10(min_tick), np.log10(max_tick), 6))
     scale_ticks = np.linspace(min_tick,max_tick, 11)
     #scale_ticks = np.linspace(0, 2, 11)
@@ -97,7 +97,7 @@ try:
     # Interpolation can be 'nearest', 'bilinear' or 'bicubic'
     cax = ax.imshow(edata, interpolation = 'nearest', cmap = cm.afmhot,
                     origin = 'lower', vmin = min_tick, vmax = max_tick,
-                    #origin = 'lower', vmin = 0.25, vmax = 1,
+                    #origin = 'lower', vmin = 0.16, vmax = 0.18,
                     extent = (min(scale_x), max(scale_x), min(scale_y), max(scale_y))
                     #,norm = LogNorm()
                     )