two improvements to BFAST API

python/examples/refl_angular_bfast.ipynb

@@ -80,8 +80,7 @@
     "        k_point=mp.Vector3(),\n",
     "        dimensions=dimensions,\n",
     "        resolution=resolution,\n",
-    "        need_bfast = True,\n",
-    "        bfast_k_bar = (np.sin(theta_r),0,0), #sets the angle of the incident wave\n",
+    "        bfast_scaled_k = (np.sin(theta_r),0,0), #sets the angle of the incident wave\n",
     "        Courant = Courant \n",
     "    )\n",
     "\n",
@@ -116,8 +115,7 @@
     "        k_point=mp.Vector3(),\n",
     "        dimensions=dimensions,\n",
     "        resolution=resolution,\n",
-    "        need_bfast = True,\n",
-    "        bfast_k_bar = (np.sin(theta_r),0,0),\n",
+    "        bfast_scaled_k = (np.sin(theta_r),0,0),\n",
     "        Courant = Courant  \n",
     "    )\n",
     "\n",

python/simulation.py

@@ -1230,8 +1230,7 @@ def __init__(
         force_complex_fields: bool = False,
         default_material: Medium = mp.Medium(),
         m: float = 0,
-        need_bfast: bool = False,
-        bfast_k_bar: Vector3Type = (0, 0, 0),
+        bfast_scaled_k: Optional[Vector3Type] = None,
         k_point: Union[Vector3Type, bool] = False,
         kz_2d: str = "complex",
         extra_materials: Optional[List[Medium]] = None,
@@ -1529,8 +1528,7 @@ def __init__(
         self.last_eps_filename = ""
         self.output_h5_hook = lambda fname: False
         self.interactive = False
-        self.need_bfast = need_bfast
-        self.bfast_k_bar = bfast_k_bar
+        self.bfast_scaled_k = (0, 0, 0) if bfast_scaled_k is None else bfast_scaled_k
         self.is_cylindrical = False
         self.material_function = material_function
         self.epsilon_func = epsilon_func
@@ -2479,8 +2477,7 @@ def init_sim(self):
             not self.accurate_fields_near_cylorigin,
             self.loop_tile_base_db,
             self.loop_tile_base_eh,
-            self.need_bfast,
-            self.bfast_k_bar,
+            self.bfast_scaled_k,
         )
 
         if self.force_all_components and self.dimensions != 1:

python/tests/test_refl_angular.py

@@ -47,13 +47,13 @@ def reflectance_angular(
         theta_rad = math.radians(theta_deg)
 
         if need_bfast:
-            bfast_k_bar = (self.n1 * np.sin(theta_rad), 0, 0)
+            bfast_scaled_k = (self.n1 * np.sin(theta_rad), 0, 0)
 
-            Courant = (1 - bfast_k_bar[0]) / 3**0.5
+            Courant = (1 - bfast_scaled_k[0]) / 3**0.5
 
             k = mp.Vector3()
         else:
-            bfast_k_bar = (0, 0, 0)
+            bfast_scaled_k = (0, 0, 0)
 
             Courant = 0.5
 
@@ -89,8 +89,7 @@ def reflectance_angular(
             sources=sources,
             boundary_layers=pml_layers,
             k_point=k,
-            need_bfast=need_bfast,
-            bfast_k_bar=bfast_k_bar,
+            bfast_scaled_k=bfast_scaled_k,
             Courant=Courant,
         )
 
@@ -126,8 +125,7 @@ def reflectance_angular(
             sources=sources,
             boundary_layers=pml_layers,
             k_point=k,
-            need_bfast=need_bfast,
-            bfast_k_bar=bfast_k_bar,
+            bfast_scaled_k=bfast_scaled_k,
             Courant=Courant,
             geometry=geometry,
         )

src/fields.cpp

@@ -30,11 +30,10 @@ using namespace std;
 namespace meep {
 
 fields::fields(structure *s, double m, double beta, bool zero_fields_near_cylorigin,
-               int loop_tile_base_db, int loop_tile_base_eh, bool need_bfast,
-               std::vector<double> bfast_k_bar)
+               int loop_tile_base_db, int loop_tile_base_eh, std::vector<double> bfast_scaled_k)
     : S(s->S), gv(s->gv), user_volume(s->user_volume), v(s->v), m(m), beta(beta),
       loop_tile_base_db(loop_tile_base_db), loop_tile_base_eh(loop_tile_base_eh),
-      working_on(&times_spent), need_bfast(need_bfast), bfast_k_bar(bfast_k_bar) {
+      working_on(&times_spent), bfast_scaled_k(bfast_scaled_k) {
   shared_chunks = s->shared_chunks;
   components_allocated = false;
   synchronized_magnetic_fields = 0;
@@ -61,7 +60,7 @@ fields::fields(structure *s, double m, double beta, bool zero_fields_near_cylori
   chunks = new fields_chunk_ptr[num_chunks];
   for (int i = 0; i < num_chunks; i++)
     chunks[i] = new fields_chunk(s->chunks[i], outdir, m, beta, zero_fields_near_cylorigin, i,
-                                 loop_tile_base_db, need_bfast, bfast_k_bar);
+                                 loop_tile_base_db, bfast_scaled_k);
   FOR_FIELD_TYPES(ft) {
     typedef realnum *realnum_ptr;
     comm_blocks[ft] = new realnum_ptr[num_chunks * num_chunks];
@@ -94,8 +93,7 @@ fields::fields(const fields &thef)
   outdir = new char[strlen(thef.outdir) + 1];
   strcpy(outdir, thef.outdir);
   m = thef.m;
-  need_bfast = thef.need_bfast;
-  bfast_k_bar = thef.bfast_k_bar;
+  bfast_scaled_k = thef.bfast_scaled_k;
   beta = thef.beta;
   phasein_time = thef.phasein_time;
   for (int d = 0; d < 5; d++) {
@@ -243,9 +241,9 @@ void check_tiles(grid_volume gv, const std::vector<grid_volume> &gvs) {
 
 fields_chunk::fields_chunk(structure_chunk *the_s, const char *od, double m, double beta,
                            bool zero_fields_near_cylorigin, int chunkidx, int loop_tile_base_db,
-                           bool need_bfast, std::vector<double> bfast_k_bar)
+                           std::vector<double> bfast_scaled_k)
     : gv(the_s->gv), v(the_s->v), m(m), zero_fields_near_cylorigin(zero_fields_near_cylorigin),
-      beta(beta), need_bfast(need_bfast), bfast_k_bar(bfast_k_bar) {
+      beta(beta), bfast_scaled_k(bfast_scaled_k) {
   s = the_s;
   chunk_idx = chunkidx;
   s->refcount++;
@@ -307,8 +305,7 @@ fields_chunk::fields_chunk(const fields_chunk &thef, int chunkidx) : gv(thef.gv)
   s->refcount++;
   outdir = thef.outdir;
   m = thef.m;
-  need_bfast = thef.need_bfast;
-  bfast_k_bar = thef.bfast_k_bar;
+  bfast_scaled_k = thef.bfast_scaled_k;
   zero_fields_near_cylorigin = thef.zero_fields_near_cylorigin;
   beta = thef.beta;
   new_s = thef.new_s;

src/meep.hpp

@@ -1502,8 +1502,7 @@ class fields_chunk {
   volume v;
   double m;                        // angular dependence in cyl. coords
   bool zero_fields_near_cylorigin; // fields=0 m pixels near r=0 for stability
-  bool need_bfast;
-  std::vector<double> bfast_k_bar;
+  std::vector<double> bfast_scaled_k;
   double beta;
   int is_real;
   std::vector<src_vol> sources[NUM_FIELD_TYPES];
@@ -1514,7 +1513,7 @@ class fields_chunk {
 
   fields_chunk(structure_chunk *, const char *outdir, double m, double beta,
                bool zero_fields_near_cylorigin, int chunkidx, int loop_tile_base_db,
-               bool need_bfast, std::vector<double> bfast_k_bar);
+               std::vector<double> bfast_scaled_k);
 
   fields_chunk(const fields_chunk &, int chunkidx);
   ~fields_chunk();
@@ -1742,8 +1741,7 @@ class fields {
   grid_volume gv, user_volume;
   volume v;
   double m;
-  bool need_bfast;
-  std::vector<double> bfast_k_bar;
+  std::vector<double> bfast_scaled_k;
   double beta;
   int t, phasein_time, is_real;
   std::complex<double> k[5], eikna[5];
@@ -1755,8 +1753,8 @@ class fields {
 
   // fields.cpp methods:
   fields(structure *, double m = 0, double beta = 0, bool zero_fields_near_cylorigin = true,
-         int loop_tile_base_db = 0, int loop_tile_base_eh = 0, bool need_bfast = false,
-         std::vector<double> bfast_k_bar = {0, 0, 0});
+         int loop_tile_base_db = 0, int loop_tile_base_eh = 0,
+         std::vector<double> bfast_scaled_k = {0, 0, 0});
   fields(const fields &);
   ~fields();
   bool equal_layout(const fields &f) const;

src/step_db.cpp

@@ -62,6 +62,7 @@ bool fields_chunk::step_db(field_type ft) {
         realnum *f_p = have_p ? f[c_p][cmp] : NULL;
         realnum *f_m = have_m ? f[c_m][cmp] : NULL;
         realnum *the_f = f[cc][cmp];
+        bool need_bfast = (bfast_scaled_k[0] || bfast_scaled_k[1] || bfast_scaled_k[2]) ? 1 : 0;
 
         if (dsig != NO_DIRECTION && s->conductivity[cc][d_c] && !f_cond[cc][cmp]) {
           f_cond[cc][cmp] = new realnum[gv.ntot()];
@@ -126,8 +127,10 @@ bool fields_chunk::step_db(field_type ft) {
                   s->conductivity[cc][d_c], s->condinv[cc][d_c], f_cond[cc][cmp]);
 
         if (need_bfast) {
-          realnum k1 = have_m ? bfast_k_bar[component_index(c_m)] : 0; // puts k1 in direction of g2
-          realnum k2 = have_p ? bfast_k_bar[component_index(c_p)] : 0; // puts k2 in direction of g1
+          realnum k1 =
+              have_m ? bfast_scaled_k[component_index(c_m)] : 0; // puts k1 in direction of g2
+          realnum k2 =
+              have_p ? bfast_scaled_k[component_index(c_p)] : 0; // puts k2 in direction of g1
           if (ft == D_stuff) {
             k1 = -k1;
             k2 = -k2;