test Maxwell stress tensor integration using pre-component notation

examples/example-eval-force.cc

@@ -41,7 +41,7 @@ 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,0.2);
+    const std::complex<double> index_Si(1.1,0.0);
     //    const std::complex<double> index_Ag = std::sqrt(epsilon_Ag);
     double WL=545; //nm
     //double WL=400; //nm
@@ -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

@@ -301,20 +301,12 @@ namespace shell_generator {
   
   std::vector<double> ShellGenerator::IntegrateByComp() {
     std::vector<double> integral = {0.0, 0.0, 0.0};
-    //simple 
     for (unsigned int i=0; i<E_.size(); ++i) {
       auto E = E_[i];
-      //auto H = 377.0*H_[i];
       auto H = H_[i];
-      //auto vert = vertices_[i];
       auto vert = face_centers_[i];
-      // Vector to unit product
       double r = norm(vert);
       std::vector<double> n = { vert[0]/r, vert[1]/r, vert[2]/r};
-      //std::cout<<"norm: " <<n[0]<<", "<< n[1] <<", "<<n[2] << std::endl;
-          
-      // std::cout << norm(unit) << std::endl;
-      //const double pi = 3.1415926535897932384626433832795;
       const std::vector< std::vector<double> > d = {{1.0, 0.0, 0.0},
                                                     {0.0, 1.0, 0.0},
                                                     {0.0, 0.0, 1.0}};
@@ -323,9 +315,10 @@ namespace shell_generator {
       std::complex<double> S(0.0);
       for (int ii = 0; ii < 3; ++ii)
         S += E[ii]*std::conj(E[ii]) + H[ii]*std::conj(H[ii]);
-      std::vector< std::vector<std::complex<double> > >  T = {{0.0, 0.0, 0.0},
-                               {0.0, 0.0, 0.0},
-                               {0.0, 0.0, 0.0}};
+      std::vector< std::vector<std::complex<double> > >
+        T = {{0.0, 0.0, 0.0},
+             {0.0, 0.0, 0.0},
+             {0.0, 0.0, 0.0}};
       for (int i = 0; i < 3; ++i) {
         for (int j = 0; j < 3; ++j) {
           T[i][j] =  E[i]*std::conj(E[j]) + H[i]*std::conj(H[j])
@@ -333,12 +326,6 @@ namespace shell_generator {
           F[i] += (1/(2.0/* *4.0*pi */))*real(T[i][j]*n[j]);
         }
       }
-      // 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;
       integral = integral + per_face_area_[i]*F;      
     }
     return integral;