can store meep fields and read them back in

examples/dipole_ice/dipole_ice.cxx

@@ -18,6 +18,12 @@ unsigned int fact(unsigned arg) {
 int main(int argc, char* argv[]) {
 
   meep::initialize mpi(argc, argv);
+
+  if(argc < 2) {
+    throw std::runtime_error("Error: need to pass path to output directory!");
+  }
+
+  std::string wf_path = argv[1];
 
   auto eps = [](scalar_t r, scalar_t z) {
     return 1.78;
@@ -26,7 +32,7 @@ int main(int argc, char* argv[]) {
   auto impulse_response = [](scalar_t t) {
     unsigned int order = 6;
     double tp = 2.0;   
-    double retval = 1.0 / (tp * fact(order - 1)) * std::pow(t * (double)order / tp, (double)order) * std::exp(-t * (double)order / tp);    
+    double retval = 1.0 / (tp * fact(order - 1)) * std::pow(t * (double)order / tp, (double)order) * std::exp(-t * (double)order / tp);
     return retval;
   };
 
@@ -35,7 +41,7 @@ int main(int argc, char* argv[]) {
 
   scalar_t t_end = 25;
   WeightingFieldCalculator wfc(geom, dipole, t_end);
-  wfc.Calculate("./wf_dipole_ice");
+  wfc.Calculate(wf_path);
 
   return 0;
 }

src/Antenna.cxx

@@ -26,7 +26,7 @@ InfEDipoleAntenna::InfEDipoleAntenna(scalar_t start_time, scalar_t end_time, sca
   Antenna(start_time, end_time, z_pos, impulse_response_func) { };
 
 void InfEDipoleAntenna::AddToGeometry(meep::fields& f, Geometry& geom) const {
-  CoordVector antenna_pos = CoordUtils::MakeCoordVectorTRZ(0.0, 0.0, z_pos);
+  CoordVector antenna_pos = CoordUtils::MakeCoordVectorTRZ(0.0, 1, z_pos);
 
   // orientation hard-coded for now ... need to change
   f.add_point_source(meep::Ez, *this, geom.toMeepCoords(antenna_pos));

src/WeightingFieldCalculator.cxx

@@ -91,8 +91,8 @@ namespace meep {
     vec rshift(shift * (0.5*fc->gv.inva));  // shift into unit cell for PBC geometries
 
     // prepare the list of field components to fetch at each grid point
-    component components[] = {Ex, Ey, Ez};
-    chunkloop_field_components data(fc, cgrid, shift_phase, S, sn, 3, components);
+    component components[] = {Ez, Er};
+    chunkloop_field_components data(fc, cgrid, shift_phase, S, sn, 2, components);
 
     // loop over all grid points in chunk
     LOOP_OVER_IVECS(fc->gv, is, ie, idx) {
@@ -113,13 +113,12 @@ namespace meep {
 
       // fetch field components at child point
       data.update_values(idx);
-      double E_x_val = data.values[0].real();
-      double E_y_val = data.values[1].real();
-      double E_z_val = data.values[2].real();
-
-      double E_r_val = std::sqrt(std::pow(E_x_val, 2) + std::pow(E_y_val, 2));
+      double E_z_val = data.values[0].real();
+      double E_r_val = data.values[1].real();
+
+      // double E_r_val = std::sqrt(std::pow(E_x_val, 2) + std::pow(E_y_val, 2));
       double E_phi_val = 0.0; // TODO: to be able to compute this, need to extract x/y coordinates of this point!
-
+      
       IndexVector chunk_ind = global_ind - chunk_start_inds;
 
       chunk_buffer_E_r(chunk_ind) = E_r_val;
@@ -149,16 +148,18 @@ void WeightingFieldCalculator::Calculate(std::string tmpdir) {
   ChunkloopData cld(0, dwf);
 
   std::size_t stepcnt = 0;
-  while (m_f -> time() < m_t_end) {
+  for(double cur_t = 0.0; cur_t <= m_t_end; cur_t += 1) {
+
+    // Time-step the fields
+    while (m_f -> time() < cur_t) {
+      m_f -> step();
+    }
 
     if(meep::am_master()) {
       std::cout << "Simulation time: " << m_f -> time() << std::endl;
     }
 
-    // Time-step the fields
-    m_f -> step();
     cld.ind_t = stepcnt++;
-
     m_f -> loop_in_chunks(meep::saving_chunkloop, static_cast<void*>(&cld), m_f -> total_volume());
   }
 }