repair Phi bug

examples/example-eval-force.cc

@@ -44,46 +44,39 @@ int main(int argc, char *argv[]) {
     const std::complex<double> index_Si(2.0);
     //    const std::complex<double> index_Ag = std::sqrt(epsilon_Ag);
     double WL=500; //nm
-    // double core_width = 5.27; //nm Si
-    // double inner_width = 8.22; //nm Ag
     double outer_width = 67.91; //nm  Si
-    // core_width = 5.27; //nm Si
-    // inner_width = 8.22; //nm Ag
-    //outer_width = WL/(2.0*pi); //nm  Si
-    outer_width = 50; //nm  Si
-    // multi_layer_mie.AddTargetLayer(core_width, index_Si);
-    // multi_layer_mie.AddTargetLayer(inner_width, index_Ag);
-    multi_layer_mie.AddTargetLayer(outer_width, index_Si);
-    multi_layer_mie.SetWavelength(WL);
-    multi_layer_mie.RunMieCalculation();
-    double Qsca = multi_layer_mie.GetQsca();
-    printf("Qsca = %g\n", Qsca);
-
     shell_generator::ShellGenerator shell;
     shell.Init();
-    // shell.Refine();
-    // shell.Refine();
-    //shell.Refine();
+    shell.Refine();
+    shell.Refine();
+    shell.Refine();
 
-    //double scale = 2.0*pi*(core_width + inner_width +  outer_width)/WL*8.0001;  //Integration sphere radius.
-    double scale = 2.0*pi*(outer_width)/WL*1.4001;  //Integration sphere radius.
-    //double scale = 1.0001;  //Integration sphere radius.
-    
-    shell.Rescale(scale);
-    // shell.RotateX(pi/2.0);
-    // shell.RotateY(pi/2.0);
-    // shell.RotateZ(pi/2.0);
-    //shell.PrintVerts();
-    auto points = shell.GetVerticesT();
-    multi_layer_mie.SetFieldPointsSP(points);
-    multi_layer_mie.RunFieldCalculation();
-    auto E = nmie::ConvertComplexVectorVector<double>(multi_layer_mie.GetFieldE());
-    auto H = nmie::ConvertComplexVectorVector<double>(multi_layer_mie.GetFieldH());
-    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;
+    for (int i=0; i<10; ++i) {
+      outer_width = 10+10*i; //nm  Si
+      multi_layer_mie.AddTargetLayer(outer_width, index_Si);
+      multi_layer_mie.SetWavelength(WL);
+      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 = 1.0001;  //Integration sphere radius.
+      shell.Rescale(scale);
+      // shell.RotateX(pi/2.0);
+      // shell.RotateY(pi/2.0);
+      // shell.RotateZ(pi/2.0);
+      //shell.PrintVerts();
+      auto points = shell.GetVerticesT();
+      multi_layer_mie.SetFieldPointsSP(points);
+      multi_layer_mie.RunFieldCalculation();
+      auto E = nmie::ConvertComplexVectorVector<double>(multi_layer_mie.GetFieldE());
+      auto H = nmie::ConvertComplexVectorVector<double>(multi_layer_mie.GetFieldH());
+      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;
+    }
   } catch( const std::invalid_argument& ia ) {
     // Will catch if  multi_layer_mie fails or other errors.
     std::cerr << "Invalid argument: " << ia.what() << std::endl;

src/nmie-impl.hpp

@@ -570,12 +570,15 @@ namespace nmie {
     Mo1n[0] = c_zero;
     Mo1n[1] = cos(Phi)*Pi*rn/Rho;
     Mo1n[2] = -sin(Phi)*Tau*rn/Rho;
+
     Me1n[0] = c_zero;
     Me1n[1] = -sin(Phi)*Pi*rn/Rho;
     Me1n[2] = -cos(Phi)*Tau*rn/Rho;
+
     No1n[0] = sin(Phi)*(n*n + n)*sin(Theta)*Pi*rn/Rho/Rho;
     No1n[1] = sin(Phi)*Tau*Dn*rn/Rho;
     No1n[2] = cos(Phi)*Pi*Dn*rn/Rho;
+
     Ne1n[0] = cos(Phi)*(n*n + n)*sin(Theta)*Pi*rn/Rho/Rho;
     Ne1n[1] = cos(Phi)*Tau*Dn*rn/Rho;
     Ne1n[2] = -sin(Phi)*Pi*Dn*rn/Rho;
@@ -1127,9 +1130,12 @@ 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));
       // 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;
+      // std::cout << "  Rho: "<<Rho<<" Theta: "<<Theta<<"  Phi:"<<Phi<<std::endl<<std::endl;
 
       //*******************************************************//
       // external scattering field = incident + scattered      //

src/shell-generator.cc

@@ -257,9 +257,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;
       // Vector to unit product
       double r = norm(vert);
       std::vector<double> n = { vert[0]/r, vert[1]/r, vert[2]/r};
@@ -286,7 +283,10 @@ namespace shell_generator {
           F[i] += T[i][j]*n[j];
         }
       }
-      std::cout<<"F: " <<F[0]<<", "<< F[1] <<", "<<F[2] << std::endl<< std::endl;
+      // 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::cout<<"F: " <<F[0]<<", "<< F[1] <<", "<<F[2] << std::endl<< std::endl;
       integral = integral + per_vertice_area_*F;
     }
     return integral;
@@ -488,10 +488,12 @@ namespace shell_generator {
 
     //std::vector< std::vector<double> > points_debug = {{1,0,0},{-1,0,0}};
     //std::vector< std::vector<double> > points_debug = {{0,1,0},{0,-1,0}};
-    std::vector< std::vector<double> > points_debug = {{1,1,0},{1,-1,0},{-1,-1,0},{-1,1,0}};
+    //std::vector< std::vector<double> > points_debug = {{1,1,0},{1,-1,0},{-1,-1,0},{-1,1,0}};
+    //std::vector< std::vector<double> > points_debug = {{0,1,1},{0,1,-1},{0,-1,-1},{0,-1,1}};
     //std::vector< std::vector<double> > points_debug = {};
     // std::vector< std::vector<double> > points_debug = {{0,0,1},{0,0,-1}};
-    vertices_ = std::move(points_debug);
+
+    //vertices_ = std::move(points_debug);
 
     // for (auto v : vertices_) {
     //   for (auto p : v)