test Poynting integration over the sphere

examples/example-eval-force.cc

@@ -41,12 +41,12 @@ int main(int argc, char *argv[]) {
     // const std::complex<double> epsilon_Si(18.4631066585, 0.6259727805);
     //    const std::complex<double> epsilon_Ag(-8.5014154589, 0.7585845411);
     // const std::complex<double> index_Si = std::sqrt(epsilon_Si);
-    const std::complex<double> index_Si(1.1);
+    const std::complex<double> index_Si(1.1,0.2);
     //    const std::complex<double> index_Ag = std::sqrt(epsilon_Ag);
     double WL=545; //nm
     //double WL=400; //nm
     //double outer_width = 67.91; //nm  Si
-    double outer_width = 4; //nm  Si
+    double outer_width = 4*2*2; //nm  Si
     auto shift = 0.0;
     shell_generator::ShellGenerator shell;
     shell.Init();
@@ -77,8 +77,8 @@ int main(int argc, char *argv[]) {
         shell.SetField(E,H);
         //shell.SetFieldSph(Es,Hs);
         // auto F = shell.Integrate();
-        //auto F = shell.IntegrateByFaces();
-        auto F = shell.IntegrateByComp();
+        auto F = shell.IntegrateByFaces();
+        //auto F = shell.IntegrateByComp();
         std::cout << "integrate_R:\t" << scale*WL/(2.0*pi);
         std::cout<<"\tforce:\t" <<F[0]<<"\t"<< F[1] <<"\t"<<F[2] << std::endl;
         // auto F1 = shell.IntegrateByComp();

src/shell-generator.cc

@@ -355,23 +355,24 @@ namespace shell_generator {
       auto E = E_[i];
       //auto H = 377.0*H_[i];
       auto H = H_[i];
-      auto Es = Es_[i];
-      auto Hs = Hs_[i];
+      // auto Es = Es_[i];
+      // auto Hs = Hs_[i];
       
       auto vert = face_centers_[i];
       // Vector to unit product
       double r = norm(vert);
-      std::vector<std::complex<double> > unit = { vert[0]/r, vert[1]/r, vert[2]/r};
+      //std::vector<std::complex<double> > unit = { vert[0]/r, vert[1]/r, vert[2]/r};
       // std::cout << norm(unit) << std::endl;
       //const double pi = 3.1415926535897932384626433832795;
-      std::vector<double> P = (1/(2.0))
-        *real(
-              dot(unit,E)*vconj(E) +
-              dot(unit,H)*vconj(H) +
-              (-1.0/2.0)*(dot(E,vconj(E))
-                          +dot(H,vconj(H))
-                          )*unit
-              );
+      // std::vector<double> P = (1/(2.0))
+      //   *real(
+      //         dot(unit,E)*vconj(E) +
+      //         dot(unit,H)*vconj(H) +
+      //         (-1.0/2.0)*(dot(E,vconj(E))
+      //                     +dot(H,vconj(H))
+      //                     )*unit
+      //         );
+
       // std::vector<double> P = (1/(2.0))
       //   *real(
       //         Es[0]*vconj(E) +
@@ -393,7 +394,14 @@ namespace shell_generator {
       // 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_face_area_[i]*P;
+      //integral = integral + per_face_area_[i]*P;
+
+      // Test Poynting vector integration
+      std::vector<double> unit = { vert[0]/r, vert[1]/r, vert[2]/r};
+      std::vector<double> P = (1/(2.0))
+        *real(cross(E,vconj(H)));
+      integral[0] = integral[0] + per_face_area_[i]*dot(P,unit);
+
     }
     return integral;
   }