integrate by components

examples/example-eval-force.cc

@@ -41,64 +41,49 @@ 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(2.0);
+    const std::complex<double> index_Si(1.1);
     //    const std::complex<double> index_Ag = std::sqrt(epsilon_Ag);
-    double WL=500; //nm
-    double outer_width = 67.91; //nm  Si
+    double WL=545; //nm
+    //double WL=400; //nm
+    //double outer_width = 67.91; //nm  Si
+    double outer_width = 4; //nm  Si
     auto shift = 0.0;
     shell_generator::ShellGenerator shell;
     shell.Init();
-    for (int refines=0; refines<4; ++refines) {
+    shell.Refine();
+    shell.Refine();
+    for (int refines=0; refines<1; ++refines) {
       shell.Refine();
-    
-    for (int i=0; i<1; ++i) {
-      auto integration_radius = 10+20*i;
-      outer_width = 10; //+10*i; //nm  Si
-      multi_layer_mie.ClearAllDesign();
-      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*(integration_radius)/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.PrintVerts();
-      shell.Rescale(scale);
-      //shell.PrintVerts();
-      std::cout << "rescale with scale factor: " << scale << std::endl;
-
-      // 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 F1 = shell.IntegrateByComp();
-      // std::cout<<"F: " <<F1[0]<<", "<< F1[1] <<", "<<F1[2] << std::endl;
-      auto charge = 2.54;
-      {
-        shift = 0.0 * scale;
-        auto F1 = shell.IntegrateGaussSimple(charge, shift);
-        std::cout<<"charge: "<< charge << "   integral_R: " << scale << "   \tshift_centerX: " << shift << "  \tcharge result: " <<F1 << std::endl;
-      }
-      {
-        shift = 0.5 * scale;
-        auto F1 = shell.IntegrateGaussSimple(charge, shift);
-        std::cout<<"charge: "<< charge << "   integral_R: " << scale << "   \tshift_centerX: " << shift << "  \tresult: " <<F1 << std::endl;
-      }
-      {
-        shift = 0.9 * scale;
-        auto F1 = shell.IntegrateGaussSimple(charge, shift);
-        std::cout<<"charge: "<< charge << "   integral_R: " << scale << "   \tshift_centerX: " << shift << "  \tresult: " <<F1 << std::endl;
+      for (int i=0; i<170; ++i) {
+        //outer_width = 40 + 5*i;
+        auto integration_radius = outer_width  + 5*i ;
+        //outer_width = 10; //+10*i; //nm  Si
+        multi_layer_mie.ClearAllDesign();
+        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\t", Qsca);
+        double scale = 2.0*pi*(integration_radius)/WL*1.00001;  //Integration sphere radius.
+        shell.Rescale(scale);
+        // auto points = shell.GetVerticesT();
+        auto points = shell.GetFaceCentersT();
+        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());
+        // auto Es = nmie::ConvertComplexVectorVector<double>(multi_layer_mie.GetFieldEs());
+        // auto Hs = nmie::ConvertComplexVectorVector<double>(multi_layer_mie.GetFieldHs());
+        shell.SetField(E,H);
+        //shell.SetFieldSph(Es,Hs);
+        // auto F = shell.Integrate();
+        //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();
+        // std::cout<<"F: " <<F1[0]<<", "<< F1[1] <<", "<<F1[2] << std::endl;        
       }
-    }
 
     }  // end for refines
   } catch( const std::invalid_argument& ia ) {

src/nmie-impl.hpp

@@ -1111,8 +1111,12 @@ namespace nmie {
     long total_points = coords_[0].size();
     E_.resize(total_points);
     H_.resize(total_points);
+    Es_.resize(total_points);
+    Hs_.resize(total_points);
     for (auto& f : E_) f.resize(3);
     for (auto& f : H_) f.resize(3);
+    for (auto& f : Es_) f.resize(3);
+    for (auto& f : Hs_) f.resize(3);
 
     for (int point = 0; point < total_points; point++) {
       const FloatType& Xp = coords_[0][point];
@@ -1147,7 +1151,10 @@ namespace nmie {
 
       // Do the actual calculation of electric and magnetic field
       calcField(Rho, Theta, Phi, Es, Hs);
-
+      for (int sph_coord = 0; sph_coord<3; ++sph_coord) {
+        Es_[point][sph_coord] = Es[sph_coord];
+        Hs_[point][sph_coord] = Hs[sph_coord];
+      }
       { //Now, convert the fields back to cartesian coordinates
         using nmm::sin;
         using nmm::cos;

src/nmie.hpp

@@ -108,6 +108,9 @@ namespace nmie {
 
     std::vector<std::vector< std::complex<FloatType> > > GetFieldE(){return E_;};   // {X[], Y[], Z[]}
     std::vector<std::vector< std::complex<FloatType> > > GetFieldH(){return H_;};
+    // Get fields in spherical coordinates.
+    std::vector<std::vector< std::complex<FloatType> > > GetFieldEs(){return E_;};   // {rho[], teha[], phi[]}
+    std::vector<std::vector< std::complex<FloatType> > > GetFieldHs(){return H_;};
 
   protected:
     // Size parameter for all layers
@@ -167,6 +170,7 @@ namespace nmie {
     /// Store result
     FloatType Qsca_ = 0.0, Qext_ = 0.0, Qabs_ = 0.0, Qbk_ = 0.0, Qpr_ = 0.0, asymmetry_factor_ = 0.0, albedo_ = 0.0;
     std::vector<std::vector< std::complex<FloatType> > > E_, H_;  // {X[], Y[], Z[]}
+    std::vector<std::vector< std::complex<FloatType> > > Es_, Hs_;  // {X[], Y[], Z[]}
     std::vector<std::complex<FloatType> > S1_, S2_;
 
 

src/shell-generator.cc

@@ -230,7 +230,6 @@ namespace shell_generator {
   // ********************************************************************** //
   // ********************************************************************** //
   double ShellGenerator::IntegrateGaussSimple(double charge, double shift) {
-    EvalFaces();
     double integral = 0.0;
     // return field at point p from the charge, located at (shift, 0,0)
     auto field = [](double charge, double shift, std::vector<double> p){
@@ -307,7 +306,8 @@ namespace shell_generator {
       auto E = E_[i];
       //auto H = 377.0*H_[i];
       auto H = H_[i];
-      auto vert = vertices_[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};
@@ -323,22 +323,77 @@ 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<double> >  T = {{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] = (1/(2.0/* *4.0*pi */))
-            * real(E[i]*std::conj(E[j]) + H[i]*std::conj(H[j])
-                   -1.0/2.0*S*d[i][j]);
-          F[i] += T[i][j]*n[j];
+          T[i][j] =  E[i]*std::conj(E[j]) + H[i]*std::conj(H[j])
+            -1.0/2.0*S*d[i][j];
+          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_vertice_area_*F;
+      integral = integral + per_face_area_[i]*F;      
+    }
+    return integral;
+  }
+  // ********************************************************************** //
+  // ********************************************************************** //
+  // ********************************************************************** //
+  
+  std::vector<double> ShellGenerator::IntegrateByFaces() {
+    std::vector<double> integral = {0.0, 0.0, 0.0};
+    //simple 
+    for (long unsigned int i=0; i<E_.size(); ++i) {
+      //std::cout << i << " ";
+      auto E = E_[i];
+      //auto H = 377.0*H_[i];
+      auto H = H_[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::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(
+      //         Es[0]*vconj(E) +
+      //         Hs[0]*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)) +
+      //     real(dot(unit,H)*vconj(H))) +
+      //   (-1.0/4.0)*(dot(E,vconj(E))*unit
+      //               +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_face_area_[i]*P;
     }
     return integral;
   }
@@ -390,6 +445,20 @@ namespace shell_generator {
   // ********************************************************************** //
   // ********************************************************************** //
   // ********************************************************************** //
+  std::vector< std::vector<double> > ShellGenerator::GetFaceCentersT() {
+    EvalFaces();
+    std::vector< std::vector<double> > vertices_t;
+    vertices_t.resize(3); 
+    for(const auto vert : face_centers_){
+      vertices_t[0].push_back(vert[0]);
+      vertices_t[1].push_back(vert[1]);
+      vertices_t[2].push_back(vert[2]);
+    }
+    return vertices_t;
+  }
+  // ********************************************************************** //
+  // ********************************************************************** //
+  // ********************************************************************** //
   double ShellGenerator::norm(std::vector<double> a){
     double norm_value = 0;
     for (auto coord:a)
@@ -404,7 +473,7 @@ namespace shell_generator {
        double factor = norm(vert);
        //std::cout<< factor <<std::endl;
        for (auto &coord:vert) {
-         coord*=scale/factor;
+         coord = coord*scale/factor;
        }
        //std::cout << " " << norm(vert) << " ";
      }

src/shell-generator.hpp

@@ -41,6 +41,7 @@ namespace shell_generator {
     double norm(std::vector<double> a);
     std::vector< std::vector<double> > GetVertices(){return vertices_;};
     std::vector< std::vector<double> > GetVerticesT();
+    std::vector< std::vector<double> > GetFaceCentersT();
     std::vector< std::vector<long unsigned int> > GetEdges(){return edges_;};
     std::vector< std::vector<long unsigned int> > GetFaces(){return faces_;};
     void EvalEdgeLength();
@@ -49,6 +50,7 @@ namespace shell_generator {
     void GenerateFacesInit();
     void Init();
     std::vector<double> Integrate();
+    std::vector<double> IntegrateByFaces();
     std::vector<double> IntegrateByComp();
     double IntegrateGauss(double charge, double dist);
     double IntegrateGaussSimple(double charge, double dist);
@@ -64,10 +66,12 @@ namespace shell_generator {
 
     void SetField(std::vector<std::vector< std::complex<double> > > &E,
                   std::vector<std::vector< std::complex<double> > > &H) {E_ = E; H_=H;};
+    void SetFieldSph(std::vector<std::vector< std::complex<double> > > &E,
+                  std::vector<std::vector< std::complex<double> > > &H) {Es_ = E; Hs_=H;};
     void SetInitialVerticesIcosahedron();
     void SetInitialVerticesTetrahedron();
   private:
-    std::vector<std::vector< std::complex<double> > > E_, H_;
+    std::vector<std::vector< std::complex<double> > > E_, H_, Es_, Hs_;
     double edge_length_ = 0.0;
     double face_area_ = 0.0;
     std::vector<double> per_face_area_;