From 481fbabbf553d5f0924ef1e3af529f0be079d869 Mon Sep 17 00:00:00 2001
From: Philipp Windischhofer <windischhofer@uchicago.edu>
Date: Wed, 31 Jan 2024 13:10:13 -0600
Subject: [PATCH] can store meep fields and read them back in

---
 examples/dipole_ice/dipole_ice.cxx | 10 ++++++++--
 src/Antenna.cxx                    |  2 +-
 src/WeightingFieldCalculator.cxx   | 25 +++++++++++++------------
 3 files changed, 22 insertions(+), 15 deletions(-)

diff --git a/examples/dipole_ice/dipole_ice.cxx b/examples/dipole_ice/dipole_ice.cxx
index 858f6b617..0083aadad 100644
--- a/examples/dipole_ice/dipole_ice.cxx
+++ b/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;
 }
diff --git a/src/Antenna.cxx b/src/Antenna.cxx
index 0d306f9e4..9d2023a32 100644
--- a/src/Antenna.cxx
+++ b/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));
diff --git a/src/WeightingFieldCalculator.cxx b/src/WeightingFieldCalculator.cxx
index 93a265b11..e4e0bda5e 100644
--- a/src/WeightingFieldCalculator.cxx
+++ b/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());
   }
 }