bug fix

examples/example-eval-force.cc

@@ -58,8 +58,8 @@ int main(int argc, char *argv[]) {
       multi_layer_mie.RunMieCalculation();
       double Qsca = multi_layer_mie.GetQsca();
       printf("Qsca = %g\n", Qsca);
-      double scale = 2.0*pi*(outer_width)/WL*1.001;  //Integration sphere radius.
-      //double scale = 2.0*pi*(110)/WL*2.001;  //Integration sphere radius.
+      //double scale = 2.0*pi*(outer_width)/WL*1.001;  //Integration sphere radius.
+      double scale = 2.0*pi*(110)/WL*2.001;  //Integration sphere radius.
       //double scale = 1.0001;  //Integration sphere radius.
       shell.Rescale(scale);
       // shell.RotateX(pi/2.0);
@@ -74,8 +74,8 @@ int main(int argc, char *argv[]) {
       shell.SetField(E,H);
       // auto F = shell.Integrate();
       // std::cout<<"F: " <<F[0]<<", "<< F[1] <<", "<<F[2] << std::endl<< std::endl;
-      auto F = shell.IntegrateByComp();
-      std::cout<<"F: " <<F[0]<<", "<< F[1] <<", "<<F[2] << std::endl;
+      auto F1 = shell.IntegrateByComp();
+      std::cout<<"F: " <<F1[0]<<", "<< F1[1] <<", "<<F1[2] << std::endl;
     }
   } catch( const std::invalid_argument& ia ) {
     // Will catch if  multi_layer_mie fails or other errors.

src/nmie-impl.hpp

@@ -1130,8 +1130,7 @@ namespace nmie {
         Phi = (Yp != 0.0) ? nmm::asin(Yp/nmm::sqrt(pow2(Xp) + pow2(Yp))) : 0.0;
       else
         Phi = nmm::acos(Xp/nmm::sqrt(pow2(Xp) + pow2(Yp)));
-      if (Yp != 0.0)
-        Phi *= Yp/nmm::sqrt(pow2(Yp));
+      if (Yp < 0.0) Phi *= -1;
       // Avoid convergence problems due to Rho too small
       if (Rho < 1e-5) Rho = 1e-5;
       // std::cout << "Xp: "<<Xp<< "  Yp: "<<Yp<< "  Zp: "<<Zp<<std::endl;

src/shell-generator.cc

@@ -303,9 +303,6 @@ namespace shell_generator {
       //auto H = 377.0*H_[i];
       auto H = H_[i];
       auto vert = vertices_[i];
-      std::cout <<"E "<<E[0]<<", "<< E[1] <<", "<<E[2] << std::endl;
-      std::cout <<"H "<<H[0]<<", "<< H[1] <<", "<<H[2] << std::endl;
-      std::cout <<"vert "<<vert[0]<<", "<< vert[1] <<", "<<vert[2] << std::endl<<std::endl;
       // Vector to unit product
       double r = norm(vert);
       std::vector<std::complex<double> > unit = { vert[0]/r, vert[1]/r, vert[2]/r};
@@ -319,6 +316,9 @@ namespace shell_generator {
                           +dot(H,vconj(H))
                           )*unit
               );
+      // auto    std::cout <<"E "<<E[0]<<", "<< E[1] <<", "<<E[2] << std::endl;
+      // std::cout <<"H "<<H[0]<<", "<< H[1] <<", "<<H[2] << std::endl;
+      // std::cout <<"vert "<<vert[0]<<", "<< vert[1] <<", "<<vert[2] << std::endl<<std::endl;
       integral = integral + per_vertice_area_*P;
     }
     return integral;