diff --git a/examples/README.md b/examples/README.md index a49f173..9836d31 100644 --- a/examples/README.md +++ b/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. diff --git a/examples/patch_antenna_pattern.csv b/examples/patch_antenna_pattern.csv new file mode 100644 index 0000000..8708216 --- /dev/null +++ b/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 diff --git a/include/vectorem/post/ntf.hpp b/include/vectorem/post/ntf.hpp new file mode 100644 index 0000000..114ae25 --- /dev/null +++ b/include/vectorem/post/ntf.hpp @@ -0,0 +1,40 @@ +#pragma once + +#include +#include +#include + +#include + +namespace vectorem { + +struct FFPoint2D { + double theta_deg; ///< Polar angle in degrees + std::complex e_theta; ///< \f$E_\theta\f$ far-field component + std::complex 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 stratton_chu_2d( + const std::vector &r, + const std::vector &n, + const std::vector &E, + const std::vector &H, + const std::vector &area, + const std::vector &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 &pat); + +} // namespace vectorem + diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 9529085..78b9124 100644 --- a/src/CMakeLists.txt +++ b/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 ) diff --git a/src/post/ntf.cpp b/src/post/ntf.cpp new file mode 100644 index 0000000..f6e884e --- /dev/null +++ b/src/post/ntf.cpp @@ -0,0 +1,61 @@ +#include "vectorem/post/ntf.hpp" + +#include +#include + +namespace vectorem { + +static constexpr double Z0 = 376.730313668; // free-space impedance [ohm] + +std::vector stratton_chu_2d( + const std::vector &r, + const std::vector &n, + const std::vector &E, + const std::vector &H, + const std::vector &area, + const std::vector &theta_rad, + double phi_rad, double k0) { + std::vector 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 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 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 &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 + diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 992830b..fb6d623 100644 --- a/tests/CMakeLists.txt +++ b/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}) diff --git a/tests/test_ntf.cpp b/tests/test_ntf.cpp new file mode 100644 index 0000000..7a62094 --- /dev/null +++ b/tests/test_ntf.cpp @@ -0,0 +1,52 @@ +#include +#include +#include + +#include "vectorem/post/ntf.hpp" + +using namespace vectorem; + +int main() { + // Zero-current sanity check + { + std::vector r = {Eigen::Vector3d::Zero()}; + std::vector n = {Eigen::Vector3d::UnitZ()}; + std::vector E = {Eigen::Vector3cd::Zero()}; + std::vector H = {Eigen::Vector3cd::Zero()}; + std::vector area = {1.0}; + std::vector 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 r = {Eigen::Vector3d(0.5, 0.0, 0.0), + Eigen::Vector3d(-0.5, 0.0, 0.0)}; + std::vector n(2, Eigen::Vector3d::UnitZ()); + std::vector E(2, Eigen::Vector3cd::UnitX()); + std::vector H(2, Eigen::Vector3cd::UnitY()); + std::vector area = {1.0, 1.0}; + std::vector 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; +} diff --git a/tools/plot_pattern.py b/tools/plot_pattern.py new file mode 100755 index 0000000..fa2f339 --- /dev/null +++ b/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])