From 26101ed8b1da69c65bc8bca5f514ec6f626dfc56 Mon Sep 17 00:00:00 2001
From: Ardavan Oskooi <ardavano@google.com>
Date: Wed, 6 Dec 2023 21:12:18 -0800
Subject: [PATCH] add documentation for bfast_scaled_k parameter of Simulation
 class constructor

---
 python/simulation.py              | 9 +++++++++
 python/tests/test_refl_angular.py | 4 ++--
 2 files changed, 11 insertions(+), 2 deletions(-)

diff --git a/python/simulation.py b/python/simulation.py
index 91ce950de..4ca15b0e5 100644
--- a/python/simulation.py
+++ b/python/simulation.py
@@ -1357,6 +1357,15 @@ def __init__(
           equivalent to taking MPB's `k_points` through its function
           `reciprocal->cartesian`.
 
+        + **`bfast_scaled_k` [`Vector3`]** — For pulsed sources at oblique incidence, a
+          fixed `k_point` results in the [various frequency components having different angles](Python_Tutorials/Basics.md##angular-reflectance-spectrum-of-a-planar-interface). To ensure all frequency
+          components of a source in a lossless medium with refractive index $n$ have the
+          same angle $\theta$ relative to e.g. the $z$ axis with rotation in the $xz$ plane
+          in a 3D cell with Bloch-periodic boundaries in $x$ and $y$ requires setting the 3-tuple
+          `bfast_scaled_k` to `(n\\sin(\theta), 0, 0)`. For stability, this also requires
+          setting the `Courant` parameter to be *less* than $(1 - \\sin(\theta)) / \\sqrt{D}$
+          where $D$ is the dimensionality of the cell.
+
         + **`kz_2d` [ `"complex"`, `"real/imag"`, or `"3d"` ]** — A 2d cell (i.e.,
           `dimensions=2`) combined with a `k_point` that has a *non-zero* component in $z$
           would normally result in a 3d simulation with complex fields. However, by
diff --git a/python/tests/test_refl_angular.py b/python/tests/test_refl_angular.py
index dd3221b44..63a730aec 100644
--- a/python/tests/test_refl_angular.py
+++ b/python/tests/test_refl_angular.py
@@ -47,9 +47,9 @@ def reflectance_angular(
         theta_rad = math.radians(theta_deg)
 
         if use_bfast:
-            bfast_scaled_k = (self.n1 * np.sin(theta_rad), 0, 0)
+            bfast_scaled_k = (self.n1 * math.sin(theta_rad), 0, 0)
 
-            Courant = (1 - bfast_scaled_k[0]) / 3**0.5
+            Courant = ((1 - math.sin(theta_rad)) / 3**0.5) * (1 / self.n1**2)
 
             k = mp.Vector3()
         else: