calc an bn bulk sphere

compare.cc

@@ -134,22 +134,6 @@ int main(int argc, char *argv[]) {
 	mode = read_x;
 	continue;
       }
-      // if (strcmp(argv[i], "-l") == 0) {
-      //   i++;
-      //   L = atoi(argv[i]);
-      //   x.resize(L);
-      //   m.resize(L);
-      //   if (argc < 3*(L + 1)) {
-      // 	  throw std::invalid_argument(error_msg);
-      //   } else {
-      //     for (l = 0; l < L; l++) {
-      //       i++;
-      //       x[l] = atof(argv[i]);
-      //       i++;
-      //       m[l] = std::complex<double>(atof(argv[i]), atof(argv[i + 1]));
-      //       i++;
-      //     }
-      //   }
       if (mode == read_ti) {
 	ti = std::stod(arg);
 	mode = read_tf;
@@ -169,28 +153,11 @@ int main(int argc, char *argv[]) {
         S2w.resize(nt);
 	continue;
       }
-      //} else if (strcmp(argv[i], "-t") == 0) {
-        // i++;
-        // ti = atof(argv[i]);
-        // i++;
-        // tf = atof(argv[i]);
-        // i++;
-        // nt = atoi(argv[i]);
-
-        // Theta.resize(nt);
-        // S1.resize(nt);
-        // S2.resize(nt);
       if (mode ==  read_comment) {
 	comment = arg;
         has_comment = 1;
 	continue;
       }
-      // } else if (strcmp(argv[i], "-c") == 0) {
-      //   i++;
-      // 	comment = args[i];
-      //   //strcpy(comment, argv[i]);
-      //   has_comment = 1;
-      // } else { i++; }
     }
     if ( (x.size() != m.size()) || (L != x.size()) ) 
       throw std::invalid_argument(std::string("Broken structure!\n")
@@ -209,12 +176,7 @@ int main(int argc, char *argv[]) {
       Theta[i] = (ti + (double)i*(tf - ti)/(nt - 1))*PI/180.0;
       }
     }
-
-
-
-
-    timespec time1, time2;
-
+    // timespec time1, time2;
     // long cpptime_nsec, best_cpp;
     // long ctime_nsec, best_c;
     // long cpptime_sec, ctime_sec;
@@ -251,36 +213,36 @@ 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;
-    double size=0.01;
+    double size=0.001;
     std::vector<double> range;
     // for (int i = 0; i < samples; ++i) {
       //range.push_back( from_coord + (to_coord-from_coord)/(static_cast<double>(samples)-1)*i );
     //range.push_back(size*0.01);
     //range.push_back(size*0.99999);
-    range.push_back(size/2.0);
+    range.push_back(size/20.0);
     //range.push_back(size*1.00001);
     //range.push_back(3);
     //printf("r=%g  ", range.back());
@@ -288,25 +250,22 @@ 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());
+    // // 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());
+    // // // 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());
     int ncoord = Xp.size();
+    // Test solid sphere
     x = {size};
-    m = {std::complex<double>(4.000000,0.00)};
-    // x = {size};
-    //m = {std::complex<double>(1.33,0.0)};
-    // x = {1.017, 2.016};
-    // m = {std::complex<double>(1.5016,1.023), std::complex<double>(2.014,2.012)};
+    m = {std::complex<double>(2.000000,0.00)};
     L = x.size();
     int pl = 0;
     int nmax = 0;
@@ -321,14 +280,17 @@ int main(int argc, char *argv[]) {
     printf ("Field total sum ()\n");
     for (auto f:E) {
       sum_e = 0.0;
-      for (auto c:f) sum_e+=std::abs(pow2(c));
+      for (auto c:f) {
+	sum_e+=std::abs(pow2(c));
+	printf("component: %g + %g i\n", std::real(c), std::imag(c));
+      }
       printf("Field E=%g\n", std::sqrt(std::abs(sum_e)));
     }
-    for (auto f:H) {
-      sum_h = 0.0;
-      for (auto c:f) sum_h+=std::abs(pow2(c));
-      printf("Field H=%g\n", std::sqrt(std::abs(sum_h))*free_impedance);
-    }
+    // for (auto f:H) {
+    //   sum_h = 0.0;
+    //   for (auto c:f) sum_h+=std::abs(pow2(c));
+    //   printf("Field H=%g\n", std::sqrt(std::abs(sum_h))*free_impedance);
+    // }
 
   } catch( const std::invalid_argument& ia ) {
     // Will catch if  multi_layer_mie fails or other errors.

go.sh

@@ -10,12 +10,12 @@ rm -f ../scattnlay
 
 #google profiler  ######## FAST!!!
 echo Uncomment next line to compile compare.cc
-g++ -Ofast -std=c++11 $file nmie.cc nmie-old.cc -lm -lrt -o scattnlay.bin /usr/lib/libtcmalloc.so.4 -fno-builtin-malloc -fno-builtin-calloc -fno-builtin-realloc -fno-builtin-free -march=native -mtune=native -msse4.2
+# g++ -Ofast -std=c++11 $file nmie.cc nmie-old.cc -lm -lrt -o scattnlay.bin /usr/lib/libtcmalloc.so.4 -fno-builtin-malloc -fno-builtin-calloc -fno-builtin-realloc -fno-builtin-free -march=native -mtune=native -msse4.2
 
 #  g++ -Ofast -std=c++11 compare.cc nmie.cc nmie-wrapper.cc -lm -lrt -o scattnlay-g.bin -ltcmalloc -fno-builtin-malloc -fno-builtin-calloc -fno-builtin-realloc -fno-builtin-free -g
 
 #DEBUG!
-#clang++ -g -O1 -fsanitize=address  -fno-optimize-sibling-calls -fno-omit-frame-pointer -std=c++11 compare.cc nmie.cc nmie-wrapper.cc -lm -lrt -o scattnlay.bin
+clang++ -g -O1 -fsanitize=address  -fno-optimize-sibling-calls -fno-omit-frame-pointer -std=c++11 compare.cc nmie.cc nmie-old.cc -lm -lrt -o scattnlay.bin
 
 cp scattnlay.bin ../scattnlay
 # cp scattnlay-g.bin ../scattnlay-g

nmie.cc

@@ -536,6 +536,76 @@ namespace nmie {
     return Num/Denom;
   }
 
+  // ********************************************************************** //
+  // Calculate an and bn for bulk sphere size x and index m                 //
+  // equation (4.56) and (4.57) BH                                          //
+  // ********************************************************************** //
+  void MultiLayerMie::calc_an_bn_bulk(std::vector<std::complex<double> >& an,
+				      std::vector<std::complex<double> >& bn,
+				      double x, std::complex<double> m) {
+
+    printf("==========\n m = %g,%g,    x= %g\n", std::real(m), std::imag(m), x);
+    std::vector<std::complex<double> > PsiX(nmax_ + 1), ZetaX(nmax_ + 1);
+    std::vector<std::complex<double> > PsiMX(nmax_ + 1), ZetaMX(nmax_ + 1);
+    // First, calculate the Riccati-Bessel functions
+    calcPsiZeta(x, PsiX, ZetaX);
+    calcPsiZeta(m*x, PsiMX, ZetaMX);
+    std::vector<std::complex<double> > D1X(nmax_ + 1),  D3X(nmax_ + 1);
+    std::vector<std::complex<double> > D1MX(nmax_ + 1),  D3MX(nmax_ + 1);
+    // Calculate the logarithmic derivatives
+    calcD1D3(x, D1X, D3X);
+    calcD1D3(m*x, D1MX, D3MX);
+    std::vector<std::complex<double> > dPsiX(nmax_ + 1), dZetaX(nmax_ + 1);
+    std::vector<std::complex<double> > dPsiMX(nmax_ + 1);
+    for (int i = 0; i < nmax_ + 1; ++i) {
+      dPsiX[i] = D1X[i]*PsiX[i];
+      dPsiMX[i] = D1MX[i]*PsiMX[i];
+      dZetaX[i] = D3X[i]*ZetaX[i];
+    }
+
+    an.resize(nmax_);
+    bn.resize(nmax_);
+    for (int i = 0; i < nmax_; i++) {
+      int n = i+1;
+      std::complex<double> Num = m*PsiMX[n]*dPsiX[n] - PsiX[n]*dPsiMX[n];
+      std::complex<double> Denom = m*PsiMX[n]*dZetaX[n] - ZetaX[n]*dPsiMX[n];      
+      an[i] = Num/Denom;
+      Num = PsiMX[n]*dPsiX[n] - m*PsiX[n]*dPsiMX[n];
+      Denom = PsiMX[n]*dZetaX[n] - m*ZetaX[n]*dPsiMX[n];      
+      bn[i] = Num/Denom;      
+    }
+    printf("dPsiX\n");
+    for (auto a: dPsiX) printf("%10.5g + %10.5gi  ",std::real(a), std::imag(a));
+    printf("\ndPsiMX\n");
+    for (auto a: dPsiMX) printf("%10.5g + %10.5gi  ",std::real(a), std::imag(a));
+    printf("\nPsiX\n");
+    for (auto a: PsiX) printf("%10.5g + %10.5gi  ",std::real(a), std::imag(a));
+    printf("\nPsiMX\n");
+    for (auto a: PsiMX) printf("%10.5g + %10.5gi  ",std::real(a), std::imag(a));
+    printf("\nZetaX\n");
+    for (auto a: ZetaX) printf("%10.5e + %10.5ei  ",std::real(a), std::imag(a));
+    //   It should be
+    // Columns 1 through 3:
+    //   3.3333e-07 - 1.0000e+03i   2.6645e-10 - 3.0000e+06i   5.1499e-07 - 1.5000e+10i
+    // Columns 4 through 6:
+    //  -3.3479e-03 - 1.0500e+14i   1.1890e+01 - 9.4500e+17i   1.1000e+06 - 1.0395e+22i
+    // Columns 7 through 9:
+    //   1.9937e+10 - 1.3514e+26i   3.6291e+14 - 2.0270e+30i   4.2838e+18 - 3.4459e+34i
+    
+    // We have ZetaX
+    // 1.00000e-03 + -1.00000e+00i
+    // 3.33333e-07 + -1.00000e+03i  -5.26776e-08 + -3.00000e+06i  -2.63721e-04 + -1.50000e+10i
+    // -1.84605e+00 + -1.05000e+14i  -1.66144e+04 + -9.45000e+17i  -1.82759e+08 + -1.03950e+22i
+    // -2.37587e+12 + -1.35135e+26i  -3.56380e+16 + -2.02703e+30i  -6.05846e+20 + -3.44594e+34i
+
+    printf("\ndZetaX\n");
+    for (auto a: dZetaX) printf("%10.5g + %10.5gi  ",std::real(a), std::imag(a));
+
+    printf("\nsize param: %g\n", x);
+
+
+  }
+
 
   // ********************************************************************** //
   // Calculates S1 - equation (25a)                                         //
@@ -1251,37 +1321,6 @@ namespace nmie {
       ml = refractive_index_[l];
     }
 
-  // std::complex<double> z1(1.0, 0.0);
-  // sbesjh(z1, jn, jnp, h1n, h1np);
-  // printf("nmax_ = %d\n", nmax_);
-  // printf("$$$$ REAL $$$$$$ Ricatti-Bessel: Calculate and compare against Wolfram Alpha\n");
-  // printf("j(0,1) = %16.18f\n", real(jn[0]));
-  // printf("WA j() = 0.841470984807896506652502321630\n");
-  // printf("j(1,1) = %16.18f\n", real(jn[1]));
-  // printf("WA j() = 0.301168678939756789251565714187322395890252640\n");
-  // printf("j(10,1) = %.14g\n", real(jn[10]));
-  // printf("WA j() = 7.116552640047313023966191737248811458533809572 × 10^-11\n");
-  // printf("j(20,1) = %.14g\n", real(jn[16]));
-  // printf("WA j() = 7.537795722236872993957557741584960855614358030 × 10^-26\n");
-  // std::complex<double> z(1.0, 2.0);
-  // sbesjh(z, jn, jnp, h1n, h1np);
-  // auto result =  jn;
-  // printf("===========Calculate and compare against Wolfram Alpha\n");
-  // printf("j(0,1+i2) = re(%16.18f)\n         im(%16.18f)\n", 
-  // 	 real(result[0]),
-  // 	 imag(result[0]));
-  // printf("WA j() = re(1.4169961192118759)\n         im(-0.874391197002)\n");
-
-  // printf("j(1,1+i2) = re(%16.18f)\n         im(%16.18f)\n", 
-  // 	 real(result[1]),
-  // 	 imag(result[1]));
-  // printf("WA j() = re(0.74785726329830368)\n         im(0.68179207555304)\n");
-
-  // printf("j(4,1+i2) = re(%16.18f)\n         im(%16.18f)\n", 
-  // 	 real(result[4]),
-  // 	 imag(result[4]));
-  // printf("WA j() = re(-0.01352410550046)\n         im(-0.027169663050653)\n");
-
     // Calculate spherical Bessel and Hankel functions and their derivatives
     sbesjh(Rho*ml, jn, jnp, h1n, h1np);
 
@@ -1346,6 +1385,20 @@ namespace nmie {
   void MultiLayerMie::RunFieldCalculation() {
     // Calculate external scattering coefficients an_ and bn_
     ExternalScattCoeffs();
+  
+    // printf("an\n");
+    // for (auto a: an_) printf("%g + %g i\n",std::real(a), std::imag(a));
+    // printf("bn\n");
+    // for (auto a: bn_) printf("%g + %g i\n",std::real(a), std::imag(a));
+    // printf("size param: %g\n", size_param_.back());
+
+    calc_an_bn_bulk(an_, bn_, size_param_.back(), refractive_index_.back());
+    // printf("bulk an\n");
+    // for (auto a: an_) printf("%g + %g i\n",std::real(a), std::imag(a));
+    // printf("bulk bn\n");
+    // for (auto a: bn_) printf("%g + %g i\n",std::real(a), std::imag(a));
+    // printf("size param: %g\n", size_param_.back());
+
     // Calculate internal scattering coefficients aln_,  bln_, cln_, and dln_
     InternalScattCoeffs();
 

nmie.h

@@ -105,6 +105,10 @@ namespace nmie {
     std::complex<double> calc_bn(int n, double XL, std::complex<double> Hb, std::complex<double> mL,
                                  std::complex<double> PsiXL, std::complex<double> ZetaXL,
                                  std::complex<double> PsiXLM1, std::complex<double> ZetaXLM1);
+    void calc_an_bn_bulk(std::vector<std::complex<double> >& an,
+			 std::vector<std::complex<double> >& bn,
+			 double x, std::complex<double> m);
+
     std::complex<double> calc_S1(int n, std::complex<double> an, std::complex<double> bn,
                                  double Pi, double Tau);
     std::complex<double> calc_S2(int n, std::complex<double> an, std::complex<double> bn,

tests/python/field-dielectric-sphere.py

@@ -47,7 +47,7 @@ x[0, 0] = 0.001
 #x[0, 1] = 2.0*np.pi*0.06/1.064
 
 m = np.ones((1, 1), dtype = np.complex128)
-m[0, 0] = 2.0
+m[0, 0] = 4.0
 #m[0, 1] = (0.565838 + 7.23262j)/1.3205
 
 npts = 501
@@ -66,7 +66,8 @@ terms, E, H = fieldnlay(x, m, coord)
 Er = np.absolute(E)
 
 # |E|/|Eo|
-Eh = np.sqrt(Er[0, :, 0]**2 + Er[0, :, 1]**2 + Er[0, :, 2]**2)
+#Eh = np.sqrt(Er[0, :, 0]**2 + Er[0, :, 1]**2 + Er[0, :, 2]**2)
+Eh = Er[0, :, 0]**2 + Er[0, :, 1]**2 + Er[0, :, 2]**2
 
 result = np.vstack((coordX, coordY, coordZ, Eh)).transpose()
 
@@ -90,12 +91,13 @@ try:
     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(0, 2, 11)
+    scale_ticks = np.linspace(min_tick,min_tick*1.2, 11)
+    #scale_ticks = np.linspace(0, 2, 11)
 
     # 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, vmax = 2,
+                    #origin = 'lower', vmin = min_tick, vmax = max_tick,
+                    origin = 'lower', vmin = 0.25, vmax = 1,
                     extent = (min(scale_x), max(scale_x), min(scale_y), max(scale_y))
                     #,norm = LogNorm()
                     )