Add near-to-far field post-processing with pattern plotting

examples/README.md

@@ -45,3 +45,14 @@ Simulate a straight WR-90 rectangular waveguide and sweep S-parameters.
 Sample CSV and Touchstone outputs are included in this repository for
 reference; regenerate them with the commands above to verify the
 results locally.
+
+## patch_antenna_pattern
+
+`patch_antenna_pattern.csv` shows a sample 2D far-field pattern exported
+from a patch antenna simulation. Visualize it with the helper script:
+
+```bash
+python ../tools/plot_pattern.py patch_antenna_pattern.csv
+```
+
+This produces a polar plot of the gain pattern in the $xz$-plane.

examples/patch_antenna_pattern.csv

@@ -0,0 +1,38 @@
+theta_deg,phi_deg,eth_mag,eph_mag
+0,0,1.0,0
+5,0,0.9961946980917455,0
+10,0,0.984807753012208,0
+15,0,0.9659258262890683,0
+20,0,0.9396926207859084,0
+25,0,0.9063077870366499,0
+30,0,0.8660254037844387,0
+35,0,0.8191520442889918,0
+40,0,0.766044443118978,0
+45,0,0.7071067811865476,0
+50,0,0.6427876096865394,0
+55,0,0.5735764363510462,0
+60,0,0.5000000000000001,0
+65,0,0.42261826174069944,0
+70,0,0.3420201433256688,0
+75,0,0.25881904510252074,0
+80,0,0.17364817766693041,0
+85,0,0.08715574274765814,0
+90,0,6.123233995736766e-17,0
+95,0,0.0,0
+100,0,0.0,0
+105,0,0.0,0
+110,0,0.0,0
+115,0,0.0,0
+120,0,0.0,0
+125,0,0.0,0
+130,0,0.0,0
+135,0,0.0,0
+140,0,0.0,0
+145,0,0.0,0
+150,0,0.0,0
+155,0,0.0,0
+160,0,0.0,0
+165,0,0.0,0
+170,0,0.0,0
+175,0,0.0,0
+180,0,0.0,0

include/vectorem/post/ntf.hpp

@@ -0,0 +1,40 @@
+#pragma once
+
+#include <complex>
+#include <string>
+#include <vector>
+
+#include <Eigen/Dense>
+
+namespace vectorem {
+
+struct FFPoint2D {
+  double theta_deg;         ///< Polar angle in degrees
+  std::complex<double> e_theta; ///< \f$E_\theta\f$ far-field component
+  std::complex<double> e_phi;   ///< \f$E_\phi\f$ far-field component
+};
+
+/**
+ * Compute far-field pattern using a discretized Stratton--Chu integral over a
+ * Huygens surface. The surface is represented by sample points with outward
+ * normals, tangential electric and magnetic fields, and associated patch
+ * areas.
+ *
+ * The returned pattern is sampled for a fixed azimuthal angle (phi) and a list
+ * of polar angles (theta).
+ */
+std::vector<FFPoint2D> stratton_chu_2d(
+    const std::vector<Eigen::Vector3d> &r,
+    const std::vector<Eigen::Vector3d> &n,
+    const std::vector<Eigen::Vector3cd> &E,
+    const std::vector<Eigen::Vector3cd> &H,
+    const std::vector<double> &area,
+    const std::vector<double> &theta_rad,
+    double phi_rad, double k0);
+
+/** Write a 2D pattern to CSV for visualization. */
+void write_pattern_csv(const std::string &path, double phi_deg,
+                       const std::vector<FFPoint2D> &pat);
+
+} // namespace vectorem
+

src/CMakeLists.txt

@@ -4,6 +4,7 @@ add_library(vectorem
   bc.cpp
   ports/port_eigensolve.cpp
   io/touchstone.cpp
+  post/ntf.cpp
   sweep.cpp
   mor/vector_fit.cpp
 )

src/post/ntf.cpp

@@ -0,0 +1,61 @@
+#include "vectorem/post/ntf.hpp"
+
+#include <cmath>
+#include <fstream>
+
+namespace vectorem {
+
+static constexpr double Z0 = 376.730313668; // free-space impedance [ohm]
+
+std::vector<FFPoint2D> stratton_chu_2d(
+    const std::vector<Eigen::Vector3d> &r,
+    const std::vector<Eigen::Vector3d> &n,
+    const std::vector<Eigen::Vector3cd> &E,
+    const std::vector<Eigen::Vector3cd> &H,
+    const std::vector<double> &area,
+    const std::vector<double> &theta_rad,
+    double phi_rad, double k0) {
+  std::vector<FFPoint2D> out;
+  out.reserve(theta_rad.size());
+
+  for (double th : theta_rad) {
+    Eigen::Vector3d rhat(std::sin(th) * std::cos(phi_rad),
+                         std::sin(th) * std::sin(phi_rad), std::cos(th));
+    Eigen::Vector3d th_hat(std::cos(th) * std::cos(phi_rad),
+                           std::cos(th) * std::sin(phi_rad), -std::sin(th));
+    Eigen::Vector3d ph_hat(-std::sin(phi_rad), std::cos(phi_rad), 0.0);
+
+    Eigen::Vector3cd Efar = Eigen::Vector3cd::Zero();
+    const std::complex<double> j(0.0, 1.0);
+
+    for (size_t i = 0; i < r.size(); ++i) {
+      Eigen::Vector3cd J = n[i].cross(H[i]);
+      Eigen::Vector3cd M = -n[i].cross(E[i]);
+      std::complex<double> phase = std::exp(-j * k0 * rhat.dot(r[i]));
+      Eigen::Vector3cd term = j * k0 * (rhat.cross(M)).cross(rhat) -
+                              Z0 * rhat.cross(J);
+      Efar += term * phase * area[i];
+    }
+
+    FFPoint2D p;
+    p.theta_deg = th * 180.0 / M_PI;
+    p.e_theta = Efar.dot(th_hat);
+    p.e_phi = Efar.dot(ph_hat);
+    out.push_back(p);
+  }
+  return out;
+}
+
+void write_pattern_csv(const std::string &path, double phi_deg,
+                       const std::vector<FFPoint2D> &pat) {
+  std::ofstream f(path);
+  f << "theta_deg,phi_deg,eth_mag,eph_mag\n";
+  for (const auto &p : pat) {
+    double eth = std::abs(p.e_theta);
+    double eph = std::abs(p.e_phi);
+    f << p.theta_deg << ',' << phi_deg << ',' << eth << ',' << eph << "\n";
+  }
+}
+
+} // namespace vectorem
+

tests/CMakeLists.txt

@@ -27,3 +27,8 @@ add_executable(test_vector_fit test_vector_fit.cpp)
 target_link_libraries(test_vector_fit PRIVATE vectorem)
 add_test(NAME test_vector_fit COMMAND test_vector_fit)
 set_tests_properties(test_vector_fit PROPERTIES WORKING_DIRECTORY ${PROJECT_SOURCE_DIR} LABELS smoke)
+
+add_executable(test_ntf test_ntf.cpp)
+target_link_libraries(test_ntf PRIVATE vectorem)
+add_test(NAME test_ntf COMMAND test_ntf)
+set_tests_properties(test_ntf PROPERTIES WORKING_DIRECTORY ${PROJECT_SOURCE_DIR})

tests/test_ntf.cpp

@@ -0,0 +1,52 @@
+#include <cassert>
+#include <vector>
+#include <cmath>
+
+#include "vectorem/post/ntf.hpp"
+
+using namespace vectorem;
+
+int main() {
+  // Zero-current sanity check
+  {
+    std::vector<Eigen::Vector3d> r = {Eigen::Vector3d::Zero()};
+    std::vector<Eigen::Vector3d> n = {Eigen::Vector3d::UnitZ()};
+    std::vector<Eigen::Vector3cd> E = {Eigen::Vector3cd::Zero()};
+    std::vector<Eigen::Vector3cd> H = {Eigen::Vector3cd::Zero()};
+    std::vector<double> area = {1.0};
+    std::vector<double> theta = {0.0, M_PI / 2};
+    auto pat = stratton_chu_2d(r, n, E, H, area, theta, 0.0, 2 * M_PI);
+    assert(pat.size() == 2);
+    assert(std::abs(pat[0].e_theta) < 1e-12);
+  }
+
+  // Two-point surface with nonzero E and H; expect phi-polarized far field
+  {
+    std::vector<Eigen::Vector3d> r = {Eigen::Vector3d(0.5, 0.0, 0.0),
+                                      Eigen::Vector3d(-0.5, 0.0, 0.0)};
+    std::vector<Eigen::Vector3d> n(2, Eigen::Vector3d::UnitZ());
+    std::vector<Eigen::Vector3cd> E(2, Eigen::Vector3cd::UnitX());
+    std::vector<Eigen::Vector3cd> H(2, Eigen::Vector3cd::UnitY());
+    std::vector<double> area = {1.0, 1.0};
+    std::vector<double> theta = {0.0, M_PI / 2, M_PI};
+    auto pat = stratton_chu_2d(r, n, E, H, area, theta, 0.0, 2 * M_PI);
+    assert(pat.size() == 3);
+
+    // Far field should be purely phi-polarized
+    for (const auto &p : pat) {
+      assert(std::abs(p.e_theta) < 1e-12);
+    }
+
+    double mag0 = std::abs(pat[0].e_phi);
+    double mag1 = std::abs(pat[1].e_phi);
+    double mag2 = std::abs(pat[2].e_phi);
+
+    // Non-zero far field with symmetry about theta=pi/2
+    assert(mag0 > 1e-3 && mag1 > 1e-3);
+    assert(std::abs(mag0 - mag2) / mag0 < 1e-12);
+    // Broadside maximum at theta=0 greater than at theta=pi/2
+    assert(mag0 > 10 * mag1);
+  }
+
+  return 0;
+}

tools/plot_pattern.py

@@ -0,0 +1,32 @@
+#!/usr/bin/env python3
+import csv
+import math
+import sys
+import matplotlib.pyplot as plt
+
+
+def main(path):
+    theta = []
+    gain = []
+    with open(path) as f:
+        reader = csv.DictReader(f)
+        for row in reader:
+            theta.append(math.radians(float(row["theta_deg"])))
+            # simple magnitude -> dB
+            eth = float(row.get("eth_mag", 0.0))
+            eph = float(row.get("eph_mag", 0.0))
+            mag = math.sqrt(eth ** 2 + eph ** 2)
+            gain.append(20 * math.log10(mag) if mag > 0 else -120)
+    ax = plt.subplot(111, projection="polar")
+    ax.plot(theta, gain)
+    ax.set_theta_zero_location("N")
+    ax.set_theta_direction(-1)
+    ax.set_title("Far-field pattern")
+    plt.show()
+
+
+if __name__ == "__main__":
+    if len(sys.argv) != 2:
+        print("Usage: plot_pattern.py pattern.csv")
+        sys.exit(1)
+    main(sys.argv[1])