diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index cc0e7dc7..cd58d6b2 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -4,7 +4,9 @@ on: [push, pull_request] # needed to allow julia-actions/cache to delete old caches that it has created permissions: actions: write - contents: read + contents: write + pull-requests: read + statuses: write jobs: test: diff --git a/.gitignore b/.gitignore index 62c70a04..df131a34 100644 --- a/.gitignore +++ b/.gitignore @@ -13,3 +13,4 @@ Manifest.toml /docs/Manifest.toml +.vscode/settings.json diff --git a/Project.toml b/Project.toml index a7e7109b..2791de1f 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "PSSFSS" uuid = "6b20a5d4-3c6c-44cd-883b-1480592d72be" authors = ["Peter Simon "] -version = "1.10.0" +version = "1.12" [deps] Dates = "ade2ca70-3891-5945-98fb-dc099432e06a" @@ -39,7 +39,7 @@ FFTW = "1.8" FileIO = "1.16.2" GeoInterface = "1.3" JLD2 = "0.5" -LibGEOS = "0.9" +LibGEOS = "0.9.4" LoggingExtras = "1.1" MetalSurfaceImpedance = "1.0" NearestNeighbors = "0.4.8" @@ -53,7 +53,7 @@ RecipesBase = "1.3.4" Reexport = "1.2.2" SLEEFPirates = "0.6.42" StaticArrays = "1.9.2" -TicraUtilities = "1.1.0" +TicraUtilities = "1.1.2" Triangulate = "2.3.2" Unitful = "1.19.0" julia = "1" diff --git a/README.md b/README.md index 27af5802..d49556de 100644 --- a/README.md +++ b/README.md @@ -68,6 +68,8 @@ the user then invokes the `analyze` function to perform the analysis. Post-proc - Version 1.9: Added `orient` keyword to `rectstrip`, allowing geometry rotation within the stationary unit cell. - Version 1.10: New `res2tep` function allows saving results in the form of a TICRA-compatible TEP (tabulated electrical properties) file. +- Version 1.11: New `res2fresnel` function allows saving results in the form of an HFSS SBR+ Fresnel table. +- Version 1.12: Added `pixels` and `sympixels` elements. ## Installation You can obtain and install PSSFSS using Julia's Pkg REPL-mode (hitting `]` as the first character at the command prompt): diff --git a/docs/PSS_&_FSS_Element_Gallery/pixelated/demo_pixelated1.jl b/docs/PSS_&_FSS_Element_Gallery/pixelated/demo_pixelated1.jl new file mode 100644 index 00000000..3f95148b --- /dev/null +++ b/docs/PSS_&_FSS_Element_Gallery/pixelated/demo_pixelated1.jl @@ -0,0 +1,50 @@ +# --- +# title: Symmetric Pixelated Element +# cover: "assets/demo_pixelated1.png" +# description: "Pixelated element created by the sympixels function" +# --- + +# For a pixelated element created with the `sympixels` function, one should always specify +# `'M'` for the element class. This is discussed in the "checkerboard" usage example. + +using PSSFSS, Plots +P = 5.4 +units = mm +nrim = 1 +halfnint = 26 +pdiv = 1 +patternvec = Bool[ + 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 1, 1, 0, 1, 1, 0, 0, 1, 0, 0, 1, 0, 0, 0, 1, + 0, 1, 0, 1, 0, 1, 1, 1, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 1, 1, 0, 1, 0, + 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 1, 1, 0, 0, 1, 0, 1, 1, 1, 1, 0, 0, 1, + 0, 1, 0, 0, 0, 1, 0, 1, 1, 1, 1, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, + 1, 1, 1, 0, 0, 1, 0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 0, 1, 1, 0, 0, 1, + 1, 0, 1, 1, 1, 0, 1, 0, 1, 0, 1, 1, 1, 0, 0, 1, 1, 0, 1, 0, 1, + 0, 0, 0, 1, 0, 0, 1, 0, 0, 1, 1, 1, 0, 1, 0, 0, 1, 0, 0, 0, + 0, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 1, 0, 1, 1, 0, 0, 1, 0, + 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 1, 1, + 1, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 1, 0, 0, 1, 1, + 0, 1, 0, 1, 0, 0, 1, 0, 0, 1, 0, 1, 0, 0, 1, 1, + 1, 1, 0, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 0, 0, + 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, + 1, 1, 1, 0, 1, 0, 0, 1, 1, 1, 1, 0, 1, + 1, 1, 1, 1, 1, 0, 0, 1, 1, 0, 1, 1, + 0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 1, + 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, + 0, 1, 0, 0, 0, 1, 0, 0, 0, + 0, 0, 1, 1, 0, 0, 1, 1, + 0, 0, 0, 0, 1, 1, 0, + 1, 0, 0, 1, 0, 0, + 0, 0, 1, 1, 1, + 1, 0, 0, 0, + 0, 0, 1, + 0, 0, + 0] + +sheet = sympixels(; P, units, nrim, halfnint, patternvec, pdiv, class='M') +plot(sheet, axis=false, xlabel="", ylabel="", xtick=[], ytick=[], linecolor=:orange, lw=0.35, size=(400,400), rep=(3,3)) #src +savefig("assets/demo_pixelated1.png") #src +p1 = plot(sheet, linecolor = "red", unitcell=true) +p2 = plot(sheet, linecolor = "blue", lw=0.35, rep = (3,3)) +plot(p1, p2, layout = (1,2), size=(800,400)) + diff --git a/docs/PSS_&_FSS_Element_Gallery/pixelated/demo_pixelated2.jl b/docs/PSS_&_FSS_Element_Gallery/pixelated/demo_pixelated2.jl new file mode 100644 index 00000000..96346084 --- /dev/null +++ b/docs/PSS_&_FSS_Element_Gallery/pixelated/demo_pixelated2.jl @@ -0,0 +1,27 @@ +# --- +# title: Pixelated Square Patch +# cover: "assets/demo_pixelated2.png" +# description: "Pixelated element created by the pixels function" +# --- + +# Specifying `'J'` for the element class as in this example will generate a warning in the Julia REPL. +# In this case the warning can be ignored because there are no pixels that intersect with their neighbors +# at only a single point. See the discussion in the "checkerboard" usage example. + +using PSSFSS, Plots +units = cm +P = 1 +units = cm +pdiv = 10 +patternmat = Bool[1 0 0 1 + 0 0 0 0 + 0 0 0 0 + 1 0 0 1] + +sheet = pixels(; P, pdiv, patternmat, units, class='J') +p1 = plot(sheet, linecolor = :red, unitcell = true) +p2 = plot(sheet, linecolor = :blue, rep=(3,3)) +plot(sheet, axis=false, xlabel="", ylabel="", xtick=[], ytick=[], unitcell=true, linecolor=:green, size=(400,400), rep=(3,3)) #src +savefig("assets/demo_pixelated2.png") #src +plot(p1, p2, layout = (1,2), size=(800,400)) + diff --git a/docs/Project.toml b/docs/Project.toml index 7cc3a595..a95b270d 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -1,11 +1,14 @@ [deps] DemoCards = "311a05b2-6137-4a5a-b473-18580a3d38b5" Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" +DocumenterCitations = "daee34ce-89f3-4625-b898-19384cb65244" Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306" PSSFSS = "6b20a5d4-3c6c-44cd-883b-1480592d72be" Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" +PrettyTables = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d" [compat] DemoCards = "0.5.2" -Documenter = "1" -Literate = "2.15.0, 2.15.1" +Documenter = "1.8" +DocumenterCitations = "1.3.6" +Literate = "2.20" diff --git a/docs/literate/angular_ss_example.jl b/docs/literate/angular_ss_example.jl index aa80b731..5746207c 100644 --- a/docs/literate/angular_ss_example.jl +++ b/docs/literate/angular_ss_example.jl @@ -1,10 +1,6 @@ #nb # %% A slide [markdown] {"slideshow": {"slide_type": "Slide"}} # ## Angle Selective Surface -# This example is taken from Zhenting Chen, Chao Du, Jie Liu, Di Zhou, and Zhongxiang Shen, -# "Design Methodology of Dual-Polarized Angle-Selective Surface Based -# on Three-Layer Frequency-Selective Surfaces", IEEE Trans. Antennas Propagat., Vol. 71, No. 11, November 2023, -# pp. 8704--8713. -# +# This example is taken [chen2023design](@cite). # A three-sheet FSS is designed that is transparent to normally incident plane waves, but strongly # attenuates obliquely incident waves. All three sheets are swastika-shaped, with the outer two # sheets being identical. The shapes are generated in PSSFSS using the manji element. diff --git a/docs/literate/band_pass_filter.jl b/docs/literate/band_pass_filter.jl index 4c6c9e45..d80e47ee 100644 --- a/docs/literate/band_pass_filter.jl +++ b/docs/literate/band_pass_filter.jl @@ -1,19 +1,16 @@ #nb # %% A slide [markdown] {"slideshow": {"slide_type": "subslide"}} # ## Loaded Cross Band Pass Filter -# This example is originally from Fig. 7.9 of B. Munk, *Frequency Selective Surfaces, -# Theory and Design,* John Wiley and Sons, 2000. The same case was analyzed in L. Li, -# D. H. Werner et al, "A Model-Based Parameter Estimation Technique for -# Wide-Band Interpolation of Periodic Moment Method Impedance Matrices With Application to -# Genetic Algorithm Optimization of Frequency Selective Surfaces", *IEEE Trans. AP-S*, -# vol. 54, no. 3, March 2006, pp. 908-924, Fig. 6. Unfortunately, in neither reference +# This example is originally from [munk2000fss; Fig. 7.9](@cite). +# The same case was analyzed in [li2006model](@cite). Unfortunately, in neither reference # are the dimensions of the loaded cross stated, except for the square unit cell # period of 8.61 mm. I estimated the dimensions from the sketch in Fig. 6 of the second # reference. To provide a reliable comparison, I analyzed one-eighth of the structure -# in HFSS, a commercial finite element solver, using all three planes of symmetry -# (using symmetry in the z = constant centerline plane +# in HFSS, a commercial finite element solver, using all three planes of symmetry. +# Exploiting the symmetry in the z = constant centerline plane # required two analyses, once for an H-wall boundary condition, and once for an E-wall. Those -# results were then combined using the method of Reed and Wheeler (even/odd symmetry)). With -# a much reduced computational domain, it was then possible to drive HFSS well into convergence. +# results were then combined using the method of Reed and Wheeler [reed1956method](@cite), i.e. +# using even/odd symmetry. With a much reduced computational domain, it was then possible to +# drive HFSS well into convergence. # # Two identical loaded cross slot-type elements are separated by a 6 mm layer of dielectric # constant 1.9. Outboard of each sheet is a 1.1 cm layer of dielectric constant 1.3. diff --git a/docs/literate/checkerboard_example.jl b/docs/literate/checkerboard_example.jl new file mode 100644 index 00000000..65feac9e --- /dev/null +++ b/docs/literate/checkerboard_example.jl @@ -0,0 +1,172 @@ +#nb # %% A slide [markdown] {"slideshow": {"slide_type": "subslide"}} +# ## Checkerboard Metasurface +# Checkerboard-style metasurfaces have been exensively studied, for example in [nakata2013plane](@cite) +# and [kempa2010percolation](@cite). They exhibit some very unusual properties, which we will demonstrate here. +# +# We consider a series of checkerboad-like geometries. The square unit cell has dimension ``P = 5~\text{mm}``. +# The side length for a PEC square sheet rotated 45° and perfectly inscribed in the unit cell is ``L_0 = P / \sqrt{2}``. +# We will analyze a series of squares at 1 GHz with side lengths ``L = L_0 + δ``. The function `computed_and_plot` +# defined below will take a given value of ``δ`` as its input, then perform four distinct analyses based on two +# complementary triangulations. Let's run it first for ``δ = -0.5 \text{mm}`` and look +# at the resulting triangulations: + +using PSSFSS, Plots, PrettyTables +using Printf: @sprintf + +units = mm +P = 5 # period in x and y direction +L0 = P / √2 # Side length for self-complementary square +Nx = Ny = 10 + +function compute_and_plot(δ; P=P, units=units, L0=L0, Nx=Nx, Ny=Ny) + + Px = Py = P + Lx = Ly = L0 - abs(δ) + orient = 45 + islekwds = (; Lx, Ly, Px, Py, units, Nx, Ny, orient) + holekwds = (; units, s1=[P,0], s2=[0,P], sides=4, + a = iszero(δ) ? [P/2] : [(L0-abs(δ)) / √2], + b=[-2*Nx], ntri=2*Nx*Ny, structuredtri=false) + if δ ≤ 0 + redj = rectstrip(; class='J', islekwds...) + redm = rectstrip(; class='M', islekwds...) + bluej = polyring(; class='J', holekwds...) + bluem = polyring(; class='M', holekwds...) + else + redj = polyring(; class='J', holekwds...) + redm = polyring(; class='M', holekwds...) + bluej = rectstrip(; class='J', islekwds...) + bluem = rectstrip(; class='M', islekwds...) + end + + p1 = plot(redj, unitcell=true, faces=true, fillcolor=:red, title="Unit Cell") + p2 = plot(redj, faces=true, edges=false, fillcolor=:red, rep=(3,3), title="3×3 Array") + p3 = plot(bluej, unitcell=true, faces=true, fillcolor=:blue, title="Unit Cell") + p4 = plot(bluej, faces=true, edges=false, fillcolor=:blue, rep=(3,3), title="3×3 Array") + pl = plot(p1,p2,p3,p4, layout=(2,2), size=(550,600), plot_title=" Triangulations for δ = $δ") + + flist = 1 # Analysis frequency in GHz + steering = (θ=0, ϕ=0) # Normal incidence + logfile = resultfile = devnull # Suppress creation of output files + showprogress = false # Suppress screen output + + redjres = analyze([Layer(), redj, Layer()], flist, steering; logfile, resultfile, showprogress) + bluejres = analyze([Layer(), bluej, Layer()], flist, steering; logfile, resultfile, showprogress) + redmres = analyze([Layer(), redm, Layer()], flist, steering; logfile, resultfile, showprogress) + bluemres = analyze([Layer(), bluem, Layer()], flist, steering; logfile, resultfile, showprogress) + + + s11rj = extract_result(redjres, @outputs s11(v,v)) |> only + s21rj = extract_result(redjres, @outputs s21(v,v)) |> only + s11rm = extract_result(redmres, @outputs s11(v,v)) |> only + s21rm = extract_result(redmres, @outputs s21(v,v)) |> only + s11bj = extract_result(bluejres, @outputs s11(v,v)) |> only + s21bj = extract_result(bluejres, @outputs s21(v,v)) |> only + s11bm = extract_result(bluemres, @outputs s11(v,v)) |> only + s21bm = extract_result(bluemres, @outputs s21(v,v)) |> only + + return pl, (; s11rj, s21rj, s11rm, s21rm, s11bj, s21bj, s11bm, s21bm) +end + +pl, _ = compute_and_plot(-0.5) +pl + +# The function creates a "red" triangulation occupying the triangle of side length ``L_0 + δ``, and a "blue" triangulation, +# occupying the complementary portion of the unit cell. Since for this case the red square +# side length is 0.5 mm shorter than the critical length ``L_0``, it lies strictly inside the unit cell. So if we choose to +# use the red triangulation to model electric surface current, then we can consider the red regions to be "islands" of +# metal in otherwise empty space. We could also use the blue triangulation to model magnetic surface current, which again would +# lead to the conclusion that the small untriangulated squares are conducting patches or "islands" of metalization. +# Either of these two choices, when analyzed with PSSFSS, should yield the same values for computed reflection or transmission +# coefficients (within modeling accuracy). +# +# A different approach would be to choose the red triangulation for representing magnetic surface current, in which case +# the small red squares would represent "holes" in an otherwise solid metallic sheet. The same "hole" interpretation would +# follow from using the blue triangulation to represent electric surface current. In fact, for this case, the blue region +# in the full plane can be regarded as the union of an infinite number of metallic squares of dimension ``L_0 + δ``. +# So positive values of ``δ`` can be handled by using ``-δ`` and reversing the roles of the red and blue triangulations. +# This is exactly what is done in the function `compute_and_plot` above. +#- +# We'll now exercise the function for the set of ``δ`` values ``\{ -0.2, -0.05, 0, 0.05, 0.2 \}``, observing both the plotted +# triangulations and the resulting scattering parameter predictions for each of the four modeling choices outlined above. + +# ```julia +# δs = [-0.2, -0.05, 0, 0.05, 0.2] # Departure in mm from self-complementary square side length +# s11rj, s21rj, s11rm, s21rm, s11bj, s21bj, s11bm, s21bm = (zeros(ComplexF64, length(δs)) for _ in 1:8) +# for (i, δ) in pairs(δs) +# plt, r = compute_and_plot(δ) +# s11rj[i] = r.s11rj +# s21rj[i] = r.s21rj +# s11rm[i] = r.s11rm +# s21rm[i] = r.s21rm +# s11bj[i] = r.s11bj +# s21bj[i] = r.s21bj +# s11bm[i] = r.s11bm +# s21bm[i] = r.s21bm +# display(plt) +# end +# ``` + +# .png) + +# .png) + +# .png) + +# .png) + +# .png) + +# Let the letter "J" denote use of a triangulation to represent electric surface current and let "M" denote magnetic surface current. +# So, e.g. "Blue M" means that the blue triangulation is used to represent magnetic current. We can make the following observations about +# the above plots: +# +# 1. The red and blue triangulations alway occupy complementary regions of the unit cell. Then for a given choice +# of ``δ``, PSSFSS analysis of "Blue J" and "Red M" should result identical scattering parameters (apart from modeling errors). +# Likewise analysis of "Red J" and "Blue M" should result in the same scattering parameters. +# 2. It is well known from Babinet's principle [nakata2013plane](@cite) [tan2012babinet](@cite) that reflection coefficients +# of thin perforated screens are equal to the negative of the transmission coefficients of the complementary structure, provided the +# structures exhibit sufficient rotational symmetry as is the case here. "Red J" and "Blue J" form such a +# complementary pair, as do "Red M" and "Blue M". +# 3. All of the ``δ`` values are quite small compared to ``L_0 \approx 3.5355 \text{mm}``. In particular, the plots for +# ``δ = \pm 0.05`` are almost indistinguishable from the plot for ``δ = 0``. Therefore, we expect the scattering +# parameters of all of these cases to be nearly equal. +# 4. For the case of ``δ = 0``, the red and blue regions are both self-complementary. From our earlier considerations, we then +# expect equal reflection and transmission coefficients for any of the choices "Red J", "Blue J", "Red M", or "Blue M". In +# fact all four of these cases should yield the same values. +# +# The following code generates a summary table showing how well these expectations are satisfied: + +# ```julia +# mat = hcat(s11rj, s11bm, -s21rm, -s21bj) |> transpose +# row_labels = ["S₁₁: Red J ", "S₁₁: Blue M", "-S₂₁: Red M ", "-S₂₁: Blue J"] +# header = (["δ = $δ mm" for δ in δs], ["islands or holes" for δ in δs]) +# header[2][findfirst(iszero, δs)] = "SC (Self-Complementary)" +# formatters = (v,i,j) -> imag(v) > 0 ? @sprintf("%7.4f + %6.4fim", real(v), imag(v)) : +# @sprintf("%7.4f - %6.4fim", real(v), -imag(v)) +# highlighters = Highlighter((data, i, j) -> (j == findfirst(iszero,δs)), crayon"red bold") +# pretty_table(mat; header, row_labels, alignment=:c, formatters, highlighters) +# ``` +for line in eachline("./assets/checkerboard_prettytable.data"); println(line); end #hide +# For ``δ < 0`` the results seem reasonable. In this regime of electrically small metal islands (for Red J and Blue M), +# one wouldn't expect much reflection, and that is what is observed. For ``δ > 0``, the geometry (for Red J and Blue M) +# consists of a metal plate with small holes, so one would expect almost total reflection, as observed. Now consider +# how well our above observations are aligned with the computed results... +# +# From observations 1 and 2, the numbers in any one column should all be approximately equal, which they are +# except for the center ``δ = 0`` column. From observation 3 we expect all the entries +# in any one row to be nearly equal, but this is not true at all. In any one row there is a violent jump in +# amplitude at or near ``δ = 0``. Finally, from observation 4 we expect all of the numeric entries for ``δ = 0`` +# to be approximately equal--which they are most definitely not. What is going on here? +# +# As discussed in the previously cited references, the idealized model being analyzed for ``δ = 0`` is unphysical. +# As ``\delta \rightarrow 0`` the response of the structure does not approach a limit, since the violent jump in +# reflection coefficient occurs for arbitrarily small positive and negative excursions from ``L_0``. Any physically +# realizable structure must exhibit continuous dependence on its physical parameters. If one attempts +# to approximate this ideal surface in the real world one must use finite thickness and conductivity, so that neighboring +# unit cells will intersect all along the thickness of the metal, rather than at a single point, thus destroying the +# self-complementary property. Similarly, finite losses in any real metal preclude self-complementarity. +# +# The same problem observed here for an ideal model of a self-complementary surface arises when analyzing pixelated +# structures such as shown in the next example, where adjacent metal pixels +# can intersect at only a single point. In that case, using a "Blue M" approach is known to agree with measurements. \ No newline at end of file diff --git a/docs/literate/compile.jl b/docs/literate/compile.jl index acc352e8..3ede7141 100644 --- a/docs/literate/compile.jl +++ b/docs/literate/compile.jl @@ -10,7 +10,8 @@ examples_list = ["symmetric_strip.jl", "resistive_square_patch.jl", "cross_on_dielectric_substrate.jl", "square_loop_absorber.jl", "flexible_absorber.jl", "splitringexample.jl", "reflectarray_example.jl", "band_pass_filter.jl", "cpss1.jl", "cpss_optimization.jl", "cpss2.jl", "splitring_cpss.jl", - "angular_ss_example.jl", "tepfile_creation_example.jl"] + "angular_ss_example.jl", "checkerboard_example.jl", "sympixels_example.jl", + "tepfile_creation_example.jl", "fresneltable_creation_example.jl"] flist = ["manual.jl"] diff --git a/docs/literate/cpss1.jl b/docs/literate/cpss1.jl index da46bb61..9d084785 100644 --- a/docs/literate/cpss1.jl +++ b/docs/literate/cpss1.jl @@ -5,10 +5,7 @@ # We'll first look at analyzing a design presented in the literature, and then proceed to optimize another # design using PSSFSS as the analysis engine inside the optimization objective function. # ### Sjöberg and Ericsson Design -# This example comes from the paper D. Sjöberg and A. Ericsson, "A multi layer meander line circular -# polarization selective structure (MLML-CPSS)," The 8th European Conference on Antennas and Propagation -# (EuCAP 2014), 2014, pp. 464-468, doi: 10.1109/EuCAP.2014.6901792. -# +# This example comes from [sjoberg2014multi](@cite). # The authors describe an ingenious structure consisting of 5 progressively rotated meanderline sheets, which # acts as a circular polarization selective surface: it passes LHCP (almost) without attenuation or # reflection, and reflects RHCP (without changing its sense!) almost without attenuation or transmission. @@ -208,7 +205,7 @@ plot!(p, foptlimits, ARgoal, color=:black, lw=4, label="Goal") # the entire unit cell, and the unit cell for sheets 2 and 4 are not square. Since the periodicity of # the sheets in the structure varies from sheet to sheet, higher order Floquet modes common to neighboring # sheets cannot be defined, so we are forced to use only the dominant (0,0) modes which are independent of -# the periodicity. This limitation is removed in a later example. +# the periodicity. This limitation is removed in the next example. # Meanwhile, it is of interest to note that their high-accuracy runs # required 10 hours for CST and 19 hours for COMSOL on large engineering workstations versus about 22 # seconds for PSSFSS on my desktop machine. diff --git a/docs/literate/cpss2.jl b/docs/literate/cpss2.jl index 49415fb7..8d35df46 100644 --- a/docs/literate/cpss2.jl +++ b/docs/literate/cpss2.jl @@ -1,9 +1,6 @@ #nb # %% A slide [markdown] {"slideshow": {"slide_type": "Slide"}} # ## Meanderline/Strip-Based CPSS -# This example comes from the same authors as the previous example. The paper is -# A. Ericsson and D. Sjöberg, "Design and Analysis of a Multilayer Meander Line -# Circular Polarization Selective Structure", IEEE Trans. Antennas Propagat., -# Vol. 65, No. 8, Aug 2017, pp. 4089-4101. +# This example comes from [ericsson2017design](@cite), by the same authors as the previous example. # The design is similar to that of the previous example except that here, the two ``\pm 45^\circ`` # rotated meanderlines are replaced with rectangular strips. # This allows us to employ the `diagstrip` element and the `orient` keyword for the @@ -133,7 +130,7 @@ plot(ps..., layout=5) # ... # ``` #- -# Layers 3 and 9 were assigned 10 modes each. Layers 5 and 7, being thinner were assigned +# Layers 3 and 9 were assigned 10 modes each. Layers 5 and 7, being thinner, were assigned # 18 modes each. The numbers of modes are determined automatically by PSSFSS to ensure # accurate cascading. #- diff --git a/docs/literate/cpss_optimization.jl b/docs/literate/cpss_optimization.jl index 084d417b..e2c55a72 100644 --- a/docs/literate/cpss_optimization.jl +++ b/docs/literate/cpss_optimization.jl @@ -159,4 +159,4 @@ #  # As hoped for, the performance meets the more stringent design goals over a broader bandwidth than the -# Sjöberg and Ericsson design, presumably because of the greater design flexibility allowed here. +# Sjöberg and Ericsson design of [sjoberg2014multi](@cite), presumably because of the greater design flexibility allowed here. diff --git a/docs/literate/cross_on_dielectric_substrate.jl b/docs/literate/cross_on_dielectric_substrate.jl index ad9d6106..d332ab0b 100644 --- a/docs/literate/cross_on_dielectric_substrate.jl +++ b/docs/literate/cross_on_dielectric_substrate.jl @@ -1,11 +1,9 @@ #nb # %% A slide [markdown] {"slideshow": {"slide_type": "subslide"}} # ## Cross on Dielectric Substrate -# This example is also taken from the paper by Alon S. Barlevy and -# Yahya Rahmat-Samii, "Fundamental Constraints on the Electrical Characteristics -# of Frequency Selective Surfaces", **Electromagnetics**, vol. 17, 1997, pp. 41-68. -# This particular example is from Section 3.2, Figures 7 and 8. It also appeared at -# higher resolution in Barlevy's PhD dissertation from which the comparison curves -# were digitized. +# This example is also taken from the paper [barlevy1997fundamental](@cite) by Barlevy and +# Rahmat-Samii. This particular example is from Section 3.2, Figures 7 and 8. It also appeared at +# higher resolution in Barlevy's PhD dissertation [barlevy1998dissertation](@cite) +# from which the comparison curves were digitized. # # We use the `loadedcross` element where we choose `w > L2/2`, so that the Cross # is "unloaded", i.e. the center section is filled in with metalization: @@ -57,7 +55,7 @@ for eps in [1, 2, 4] end # The above loop requires about 18 seconds of execution time on my machine. -# Compare PSSFSS results to those digitized from the dissertation figure: +# Compare PSSFSS results to those digitized from the dissertation: col=[:red,:blue,:green] p = plot(xlim=(0.,30), xtick = 0:5:30, ylim=(0,1), ytick=0:0.1:1, diff --git a/docs/literate/flexible_absorber.jl b/docs/literate/flexible_absorber.jl index 31f7a7d2..08a7a7cf 100644 --- a/docs/literate/flexible_absorber.jl +++ b/docs/literate/flexible_absorber.jl @@ -1,8 +1,6 @@ #nb # %% A slide [markdown] {"slideshow": {"slide_type": "subslide"}} # ## Flexible Absorber -# This example is from Yize Li, et al., "Ultra-wideband, polarization-insensitive flexible metamaterial -# absorber base on laser printed graphene using equivalent circuit design method," Carbon, Vol 212, 2023, -# available for free download from [here](https://doi.org/10.1016/j.carbon.2023.118166). +# This example is from [li2023ultra](@cite). # It uses square and circular resistive FSS elements sandwiched between layers of flexible dielectrics to # realize a reflective absorber (i.e. a "rabsorber"). # We compare the reflection coefficient magnitude computed by PSSFSS to that digitized @@ -38,13 +36,13 @@ scan = (θ=0, ϕ=0) results = analyze(strata, FGHz, scan, showprogress=false, resultfile=devnull, logfile=devnull) s11db = extract_result(results, @outputs s11db(te,te)) -yize = readdlm("../src/assets/yize2023_fig2a_s11db_digitized.csv", ',', Float64, '\n') +li = readdlm("../src/assets/li2023_fig2a_s11db_digitized.csv", ',', Float64, '\n') ps11 = plot(title="Normal Incidence Reflection Magnitude", xlabel="Frequency (GHz)", ylabel="20log₁₀|s₁₁|", xlim=(0,20), ylim=(-18,0), xtick=0:2:20, ytick=-20:2:0, framestyle=:box) plot!(ps11, FGHz, s11db, color=:blue, label="PSSFSS", lw=2) -plot!(ps11, yize[:,1], yize[:,2], color=:red, label="Yize et al.", lw=2) +plot!(ps11, li[:,1], li[:,2], color=:red, label="Li et al.", lw=2) #md savefig(ps11, "flexibleabsorbers11.png"); nothing # hide #- #md #  diff --git a/docs/literate/fresneltable_creation_example.jl b/docs/literate/fresneltable_creation_example.jl new file mode 100644 index 00000000..60d344c9 --- /dev/null +++ b/docs/literate/fresneltable_creation_example.jl @@ -0,0 +1,28 @@ +#nb # %% A slide [markdown] {"slideshow": {"slide_type": "Slide"}} +# ## Fresnel Table Creation +# Here we show how to create an HFSS-compatible Fresnel table using the [`res2fresnel`](@ref) +# function included with PSSFSS. The geometry for this example is a double square loop +# in a 5 cm square unit cell. The code for analyzing this geometry and creating the +# Fresnel table is shown below: +# ```Julia +# using PSSFSS +# dwidth = 3mm +# duroid = Layer(epsr=2.2, tandel=0.0009, width=dwidth) +# a = √2 * [1, 2.125] +# b = √2 * [1.5, 3.125] +# units = mm +# sides = 4 +# orient = 45 +# sheet = polyring(;a, b, units, sides, orient, s1=[8,0], s2=[0,8], ntri=2000) +# strata = [Layer(), duroid, sheet, duroid, Layer(width=-2*dwidth)] +# steering = (θ=0:5:45, ϕ=0) +# FGHz = 10:2:20 +# results = analyze(strata, FGHz, steering; resultfile="double_square_loop.res") +# res2fresnel(results, "double_square_loop.rttbl") +# # Alternatively: res2fresnel("double_square_loop.res", "double_square_loop.rttbl") +# ``` +# Note that the final layer specifies a width that is the negative of the sum of the remaining layers' widths. +# This has the effect of moving the output phase reference plane to coincide with the input phase reference plane. +# Please see the documentation for [`res2fresnel`](@ref) for discussion of the setup requirements and restrictions +# necessary for creating a valid Fresnel table. + diff --git a/docs/literate/manual.jl b/docs/literate/manual.jl index 4bbed9c2..48e97b6f 100644 --- a/docs/literate/manual.jl +++ b/docs/literate/manual.jl @@ -25,9 +25,7 @@ #nb # %% A slide [markdown] {"slideshow": {"slide_type": "slide"}} # ### A Quick Example # Here is an example run of PSSFSS for a 4-sheet meanderline polarizer design -# from one of the first papers on the subject: T. L. Blackney, J. R. Burnett, and -# S. B. Cohn, “A design method for meander-line circular polarizers” -# presented at 22nd Annual Antenna Symposium, Oct. 1972. +# from [blbc:72](@cite), one of the first papers on the subject. # Detailed explanations of the code are omitted for now. #- #nb # %% A slide [markdown] {"slideshow": {"slide_type": "subslide"}} @@ -75,7 +73,7 @@ # #  # -# The PSSFSS run took about 24 seconds on my machine for 4 meanderline PSS sheets +# The PSSFSS run took about 14 seconds on my machine for 4 meanderline PSS sheets # analyzed at 61 frequencies. # # I hope this example whetted your appetite to learn @@ -130,7 +128,7 @@ # the [Pluto Home Page](https://plutojl.org/) for details of using Pluto. #nb # %% A slide [markdown] {"slideshow": {"slide_type": "fragment"}} -# It is almost essential that PSSFSS users also install the [Plots](https://github.com/JuliaPlots/Plots.jl) +# It is essential that PSSFSS users also install the [Plots](https://github.com/JuliaPlots/Plots.jl) # package. This will allow easy visualization of the FSS/PSS element triangulations produced by PSSFSS, in # addition to providing a convenient means to plot analysis results. @@ -167,8 +165,12 @@ #nb # %% A slide [markdown] {"slideshow": {"slide_type": "slide"}} # ## Strata -# ### Layer -# Dielectric layers are created with the [`Layer`](@ref) function: +# The geometry to be analyzed is passed as the first argument of the [`analyze`](@ref) function in the form +# of a Julia `Vector` (i.e., a one-dimensional array). This vector must contain at least two `Layer` +# objects representing dielectric layers, and zero or more `RWGSheet` instances, representing the zero-thickness +# PSS or FSS sheets located between the dielectric layers. These are described below. +# ### Dielectric Layers +# Dielectric layers are created with the [`Layer`](@ref) constructor function: using PSSFSS # Brings PSSFSS functions and types into scope Layer() # Defaults to zero-thickness vacuum layer #- @@ -191,7 +193,7 @@ foam = Layer(width=0.25inch, tandel=0.001, epsr=1.05) # You can stick to ASCII i patch = rectstrip(Nx=10, Ny=10, Px=1, Py=1, Lx=0.5, Ly=0.5, units=cm) #nb # %% A slide [markdown] {"slideshow": {"slide_type": "fragment"}} -# The call to [`rectstrip`](@ref) above creates a `RWGSheet` object for a rectangular strip +# The call to [`rectstrip`](@ref) above creates an `RWGSheet` object for a rectangular strip # of dimensions 0.5 cm in the x and y directions, lying in a square unit cell of dimension # 1 cm. The triangulation uses 10 edges in the x and y directions (`Nx` and `Ny`). # @@ -234,8 +236,7 @@ patch = rectstrip(Nx=10, Ny=10, Px=1, Py=1, Lx=0.5, Ly=0.5, units=cm) # default value, suitable for the so-called "oxide side" of a planar conductor, or `:rayleigh` (suitable for the # "foil side" of the conductor, i.e. the side bonded to the dielectric substrate). Together, the conductivity, # surface roughness, and roughness distribution type are used by PSSFSS internally to compute a frequency-dependent surface -# impedance, using the so-called Gradient Model, as described in D. N. Grujić, “Simple and Accurate Approximation of -# Rough Conductor Surface Impedance,” IEEE Trans. Microwave Theory Tech., vol. 70, no. 4, pp. 2053-2059, April 2022. +# impedance, using the so-called Gradient Model, as described in [grujic2021simple](@cite). # Obviously, only one of `Zsheet` and `sigma` may be specified as keyword arguments for a given `RWGSheet`. # #### Perfectly Conducting Walls @@ -441,7 +442,7 @@ flist = union(7:0.5:10, 20:0.5:25) # Two frequency bands #nb # %% A slide [markdown] {"slideshow": {"slide_type": "fragment"}} # It is often desired to use a set of polarization basis vectors other than TE/TM to define field # coefficients. PSSFSS supports in addition to TE/TM the use of H/V for horizontal/vertical components -# (in the [Ludwig 3](https://ieeexplore.ieee.org/document/1140406) sense), and L/R for left-hand circular +# (in the Ludwig 3 [ludw:73](@cite) sense), and L/R for left-hand circular # and right-hand circular polarization (in the # [IEEE sense](https://en.wikipedia.org/wiki/Circular_polarization#Uses_of_the_two_conventions)). @@ -542,7 +543,8 @@ flist = union(7:0.5:10, 20:0.5:25) # Two frequency bands # #### Using the `analyze` Function Return Value # The variable returned from [`analyze`](@ref) can be used in a similar way to the result file to obtain any quantity -# available to the [`@outputs`](@ref) macro, but for this purpose we use the [`extract_result`](@ref) function: +# available to the [`@outputs`](@ref) macro. For this purpose we pass the return value of the `analyze` function +# as the first argument to the [`extract_result`](@ref) function: # # ```julia # julia> results = analyze(strata, flist, steering, outlist=outputs); @@ -565,7 +567,28 @@ flist = union(7:0.5:10, 20:0.5:25) # Two frequency bands # with each column corresponding to a parameter of the [`@outputs`](@ref) macro. #nb # %% A slide [markdown] {"slideshow": {"slide_type": "subslide"}} -# ### Saving Analysis Results As TEP Files +# ## Exporting Analysis Results +# Besides the ability to export results to CSV files as described [here](@ref "The `outlist` Keyword Argument to `analyze`"), +# one can export results to TICRA-compatible TEP files or HFSS-compatible Fresnel tables, as described in the following +# subsections. +#nb # %% A slide [markdown] {"slideshow": {"slide_type": "subslide"}} +# ### Exporting Analysis Results to TEP Files # The vector of `Result` objects returned by the [`analyze`](@ref) function (or the corresponding # result file) can be converted to a TICRA-compatible TEP (tabulated electrical properties) file -# using the function [`res2tep`](@ref). The code needed to do so is shown in [this example](@ref "TEP File Creation"). +# using the function [`res2tep`](@ref). TEP files contain the equivalent of the full +# 4×4 scattering matrix computed by PSSFSS. Therefore, there is no requirement/limitation for unit cell geometry +# to exhibit any particular symmetry when results are to be exported to a TEP file. There are, however, requirements +# on the choice of steering angles to be analyzed. For details, see the documentation of [`res2tep`](@ref). +# Sample code for creating a TEP file is shown in [this example](@ref "TEP File Creation"). + +#nb # %% A slide [markdown] {"slideshow": {"slide_type": "subslide"}} +# ### Exporting Analysis Results to Fresnel Tables +# The vector of `Result` objects returned by the [`analyze`](@ref) function (or the corresponding +# result file) can be converted to an HFSS SBR+-compatible Fresnel table +# using the function [`res2fresnel`](@ref). Fresnel tables contain reflection +# and transmission coefficients for only "front" surface incidence, and only for co-polarized +# (TE → TE and TM → TM) scattering. In addition, coefficients for only one azimuth angle ϕ are +# retained in such a table. These limitations impose restrictions on the type of structure that +# can be analyzed for use in generating a valid Fresnel table. For details, see the documentation +# of [`res2fresnel`](@ref). Sample code for creating a Fresnel table is shown in +# [this example](@ref "Fresnel Table Creation"). diff --git a/docs/literate/reflectarray_example.jl b/docs/literate/reflectarray_example.jl index 5ce36b31..35a86b31 100644 --- a/docs/literate/reflectarray_example.jl +++ b/docs/literate/reflectarray_example.jl @@ -1,13 +1,9 @@ #nb # %% A slide [markdown] {"slideshow": {"slide_type": "Slide"}} # ## Reflectarray Element -# This example is taken from Figure 6 of -# Li, Jiao, and Zhao: "A Novel Microstrip Rectangular-Patch/Ring- -# Combination Reflectarray Element and Its Application", **IEEE Antennas and Wireless Propagation Letters**, -# VOL. 8, 2009, pp. 1119-1112. -# +# This example is taken from [li2009novel; Figure 6](@cite). # It generates the so-called "S-curve" for reflection phase of a reflectarray element. The element # consists of a square patch in a square ring, separated from a ground plane by two dielectric layers. -# Reflection phase is plotted versus the `L2` parameter as defined by Li, et al., which characterizes the +# Reflection phase is plotted versus the `L2` parameter as defined in [li2009novel](@cite), which characterizes the # overall size of the element. # diff --git a/docs/literate/resistive_square_patch.jl b/docs/literate/resistive_square_patch.jl index daaa7956..b96aafad 100644 --- a/docs/literate/resistive_square_patch.jl +++ b/docs/literate/resistive_square_patch.jl @@ -2,10 +2,8 @@ # ## Resistive Square Patch # This example will demonstrate the ability of PSSFSS to accurately model finite # conductivity of FSS metalization. It consists of a square finitely conducting -# patch in a square lattice. It is taken from a paper by Alon S. Barlevy and -# Yahya Rahmat-Samii, -# "Fundamental Constraints on the Electrical Characteristics of Frequency Selective -# Surfaces", **Electromagnetics**, vol. 17, 1997, pp. 41-68. This particular example +# patch in a square lattice. It is taken from a paper by Barlevy and +# Rahmat-Samii [barlevy1997fundamental](@cite). This particular example # is from Section 3.2, Figures 7 and 8. We will compare PSSFSS results to those digitized # from the cited figures. # diff --git a/docs/literate/splitring_cpss.jl b/docs/literate/splitring_cpss.jl index d7c42f28..f6e94fc8 100644 --- a/docs/literate/splitring_cpss.jl +++ b/docs/literate/splitring_cpss.jl @@ -1,10 +1,7 @@ #nb # %% A slide [markdown] {"slideshow": {"slide_type": "Slide"}} # ## Split Ring-Based CPSS # This circular polarization selective surface (CPSS) example comes from the paper -# L.-X. Wu, K. Chen, T. Jiang, J. Zhao and Y. Feng, "Circular-Polarization-Selective -# Metasurface and Its Applications to Transmit-Reflect-Array Antenna and Bidirectional -# Antenna," in IEEE Trans. Antennas and Propag., vol. 70, no. 11, pp. 10207-10217, -# Nov. 2022, doi: 10.1109/TAP.2022.3191213. +# [wu2022circular](@cite) by Wu et al. # The design consists of three sequentially rotated split rings separated by dielectric # layers. Since the unit cells for all three rings are identical, PSSFSS can rigorously # account for multiple scattering between the individual sheets using multiple diff --git a/docs/literate/splitringexample.jl b/docs/literate/splitringexample.jl index 7191101b..acb652d3 100644 --- a/docs/literate/splitringexample.jl +++ b/docs/literate/splitringexample.jl @@ -1,10 +1,6 @@ #nb # %% A slide [markdown] {"slideshow": {"slide_type": "Slide"}} # ## Split-Ring Resonator -# This example is taken from Figure 3 of -# Fabian‐Gongora, et al, "Independently Tunable Closely Spaced Triband -# Frequency Selective Surface Unit Cell Using the Third Resonant Mode of Split Ring Slots", -# IEEE Access, 8/3/2021, Digital Object Identifier 10.1109/ACCESS.2021.3100325. -# +# This example is taken from [fabian2021independently; Figure 3](@cite). # It consists of three concentric split rings with gaps sequentially rotated by 180° situated # on a thin dielectric slab. # diff --git a/docs/literate/square_loop_absorber.jl b/docs/literate/square_loop_absorber.jl index 70de327c..13203da0 100644 --- a/docs/literate/square_loop_absorber.jl +++ b/docs/literate/square_loop_absorber.jl @@ -1,13 +1,12 @@ #nb # %% A slide [markdown] {"slideshow": {"slide_type": "subslide"}} # ## Square Loop Absorber -# This example is from Figure 7 of Costa and Monorchio: "A frequency selective -# radome with wideband absorbing properties", *IEEE Trans. AP-S*, -# Vol. 60, no. 6, June 2012, pp. 2740--2747. It shows how one can use the `polyring` +# This example, from Figure 7 of [costa2012frequency](@cite), +# shows how one can use the [`polyring`](@ref) # function to model square loop elements. Three different designs are examined # that employ different loop thicknesses and different values of sheet resistance. # We compare the reflection coefficient magnitudes computed by PSSFSS with those digitized # from the cited figure when the sheet is suspended -# 5 mm above a ground plane, hence we will also make use of the `pecsheet` function. +# 5 mm above a ground plane, hence we will also make use of the [`pecsheet`](@ref) function. using Plots, PSSFSS, DelimitedFiles D = 11 # Period of square lattice (mm) diff --git a/docs/literate/symmetric_strip.jl b/docs/literate/symmetric_strip.jl index 0a388724..391a01b8 100644 --- a/docs/literate/symmetric_strip.jl +++ b/docs/literate/symmetric_strip.jl @@ -10,10 +10,10 @@ # The dashed lines show two possible choices for the unit cell location: "J" for a formulation in terms of electric # surface currents, and "M" for magnetic surface currents. # -# For normal incidence there is a closed-form solution due to Weinstein, -# but for a more recent reference one can find the solution in Problem 10.6 of R. E. Collin, -# *Field Theory of Guided Waves, Second Ed.*, -# IEEE Press, 1991. Here is the code for computing the exact solution: +# For normal incidence there is a closed-form solution due to Weinstein [weinstein1969theory](@cite), +# but there are more recent (and accessible) references in which one can also find the +# solution. For example [collin1991field; Problem 10.6](@cite) and [daniele1990diffraction](@cite). +# Here is the code for computing the exact solution, based on [collin1991field](@cite): """ grating(kP, nterms=30) -> (Γ, T) @@ -43,8 +43,8 @@ function grating(kP; nterms=30) return (Γ, T) end -# Note that using the extension of -# [Babinet's Principle for electromagnetic fields](http://kirkmcd.princeton.edu/examples/babinet.pdf) +# Note that using the extension of +# Babinet's Principle for electromagnetic fields [tan2012babinet](@cite), [nakata2013plane](@cite), # this also provides the solution (upon appropriate interchange and sign change of the coefficients) for # the case where the incident wave polarization is parallel to the direction of the strips. diff --git a/docs/literate/sympixels_example.jl b/docs/literate/sympixels_example.jl new file mode 100644 index 00000000..37f3c935 --- /dev/null +++ b/docs/literate/sympixels_example.jl @@ -0,0 +1,88 @@ +#nb # %% A slide [markdown] {"slideshow": {"slide_type": "Slide"}} +# ## Pixelated Band Pass Filter +# This example of from [hong2021design; Fig. 4a](@cite). It consists of a single-layer +# band-pass filter consisting of many square metallic "pixels", optimized using a special binary +# optimizer described in the cited paper. +# +# The geometry for the pixelated filter is created using the [`sympixels`](@ref) function as +# show below. +using PSSFSS, Plots +P = 5.4 +units = mm +nrim = 1 +halfnint = 26 +patternvec = Bool[ + 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 1, 1, 0, 1, 1, 0, 0, 1, 0, 0, 1, 0, 0, 0, 1, + 0, 1, 0, 1, 0, 1, 1, 1, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 1, 1, 0, 1, 0, + 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 1, 1, 0, 0, 1, 0, 1, 1, 1, 1, 0, 0, 1, + 0, 1, 0, 0, 0, 1, 0, 1, 1, 1, 1, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, + 1, 1, 1, 0, 0, 1, 0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 0, 1, 1, 0, 0, 1, + 1, 0, 1, 1, 1, 0, 1, 0, 1, 0, 1, 1, 1, 0, 0, 1, 1, 0, 1, 0, 1, + 0, 0, 0, 1, 0, 0, 1, 0, 0, 1, 1, 1, 0, 1, 0, 0, 1, 0, 0, 0, + 0, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 1, 0, 1, 1, 0, 0, 1, 0, + 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 1, 1, + 1, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 1, 0, 0, 1, 1, + 0, 1, 0, 1, 0, 0, 1, 0, 0, 1, 0, 1, 0, 0, 1, 1, + 1, 1, 0, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 0, 0, + 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, + 1, 1, 1, 0, 1, 0, 0, 1, 1, 1, 1, 0, 1, + 1, 1, 1, 1, 1, 0, 0, 1, 1, 0, 1, 1, + 0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 1, + 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, + 0, 1, 0, 0, 0, 1, 0, 0, 0, + 0, 0, 1, 1, 0, 0, 1, 1, + 0, 0, 0, 0, 1, 1, 0, + 1, 0, 0, 1, 0, 0, + 0, 0, 1, 1, 1, + 1, 0, 0, 0, + 0, 0, 1, + 0, 0, + 0] + + +sheet1 = sympixels(; P, units, nrim, halfnint, patternvec, pdiv=1, class='M') +plot(sheet1, size=(450,450), unitcell=true) +#- +# As described in the documentation for [`sympixels`](@ref), the `1` entries in `patternvec` above +# denote the locations of the metallized pixels. However, by setting `class = 'M'` above we choose to +# triangulate the empty pixel locations, i.e. those that are not occupied by metal. As discussed in +# [Checkerboard Metasurface](@ref), this choice is necessary to avoid spurious results when analyzing +# a structure of this type. +# +# Note that the above triangulation is created by forming a square for each triangulated pixel and then +# adding a single diagonal. A finer triangulation is required for an accurate analysis result. This can +# be accomplished by increasing `pdiv` to a value greater than `1`. Here is the triangulation that results +# from `pdiv = 2`: +sheet2 = sympixels(; P, units, nrim, halfnint, patternvec, pdiv=2, class='M') +plot(sheet2, lw=0.5, size=(450,450), unitcell=true) +#- +# Now each pixel has been divided into an array of 2×2 squares, each of which receives a diagonal to +# form triangles. The code for analyzing the `pdiv = 1` case would look like this: +# ```Julia +# strata = [Layer(), sheet1, Layer(epsr=3.28, tandel=0.007, width=0.05mm), Layer()] +# steering = (θ=0, ϕ=0) +# flist = range(start=20, stop=36, length=401) +# logfile = "bp_filter_pdiv1.log" +# resultfile = "bp_filter_pdiv1.res" +# results = analyze(strata, flist, steering; logfile, resultfile) +# s21db = extract_result(results, @outputs s21db(te,te)) +# ``` +#- +# Accurately analyzing a structure like this consisting of a huge number of pixels can be become +# very expensive (in terms of memory usage and CPU time) very quickly. For example, the `pdiv = 1` +# case generated 3056 basis functions and took about 30 seconds to analyze on a 32 GByte machine. +# But the `pdiv = 2` case generated 15,088 unknowns and required 683 seconds of execution time. It +# also failed on the 32 GByte machine due to lack of memory, but was run successfully on a 64 GByte +# machine. +#- +# Because the analysis in [hong2021design](@cite) used nonzero metal thickness, I ran the zero-thickness +# geometry in HFSS to observe the significance of the thickness. Here is a plot comparing the digitized +# results from [hong2021design; Fig. 4a](@cite), with the zero-thickness HFSS analysis and the two +# PSSFSS analyses: +#- +#  +#- +# As seen above, the finite metallization thickness does not have a large effect on the transmission trace. +# Using `pdiv = 2` produces much better agreement with HFSS, but it appears that an even larger value would +# be required to achieve full convergence. + diff --git a/docs/literate/tepfile_creation_example.jl b/docs/literate/tepfile_creation_example.jl index 93270cbe..d9e081f0 100644 --- a/docs/literate/tepfile_creation_example.jl +++ b/docs/literate/tepfile_creation_example.jl @@ -1,8 +1,8 @@ #nb # %% A slide [markdown] {"slideshow": {"slide_type": "Slide"}} # ## TEP File Creation -# Here we show how to create a TICRA-compatible TEP (tabulated electrical properties) file using PSSFSS. -# The geometry for this example is a rectangular copper strip measuring 4 cm × 0.2 cm in a 5 cm square unit cell. -# The code for analyzing this geometry and creating the TEP file is shown below: +# Here we show how to create a TICRA-compatible TEP (tabulated electrical properties) file using the [`res2tep`](@ref) +# function included with PSSFSS. The geometry for this example is a rectangular copper strip measuring +# 4 cm × 0.2 cm in a 5 cm square unit cell. The code for analyzing this geometry and creating the TEP file is shown below: # ```Julia # using PSSFSS # FGHz = 3.0 @@ -14,11 +14,14 @@ # sheet = rectstrip(; Px, Py, Lx, Ly, Nx, Ny, units=cm, sigma=5.7e8) # steering = (θ=0:5:70, ϕ=0:15:345) # strata = [Layer(), sheet, Layer()] -# results = analyze(strata, FGHz, steering) +# resultfile = "strip.res" +# results = analyze(strata, FGHz, steering; resultfile) # res2tep(results, "dipole_pssfss.tep"; name = "dipole", class = "pssfss") +# # Alternatively: res2tep(resultfile, "dipole_pssfss.tep"; name = "dipole", class = "pssfss") # ``` # Note that a very fine discretization has been specified, resulting in 1560 triangles. There are also a large # number of steering angles requested (15 × 24 = 360). This analysis required about 155 seconds on my machine. +# Please see the documentation for [`res2tep`](@ref) for requirements on setting up analysis scan angles. # # For comparison, the same geometry was analyzed using the QUPES program of Ticra Tools Student Edition # 2024 (hence the specification for `sigma` above, which is the default conductivity used by QUPES). For the diff --git a/docs/make.jl b/docs/make.jl index 8c1b7fa7..141c5a16 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -1,6 +1,15 @@ -using PSSFSS -using Documenter -using DemoCards +import Pkg +Pkg.activate(@__DIR__) + +using Documenter, DocumenterCitations, DemoCards, PSSFSS + +olddir = pwd() +cd(@__DIR__) + +bib = CitationBibliography( + joinpath(@__DIR__, "src", "refs.bib"); + style=:numeric +) #= using Literate @@ -10,9 +19,10 @@ end =# demopage, postprocess_cb, demo_assets = makedemos("PSS_&_FSS_Element_Gallery") -@show demopage + assets = String[] isnothing(demo_assets) || (push!(assets, demo_assets)) +push!(assets, "assets/citations.css") makedocs(; clean=false, @@ -21,6 +31,7 @@ makedocs(; modules=[PSSFSS], authors="Peter Simon and contributors", sitename="PSSFSS.jl", + plugins = [bib,], format=Documenter.HTML(; repolink="https://github.com/simonp0420/PSSFSS.jl/blob/{commit}{path}#L{line}", prettyurls = true, @@ -29,12 +40,14 @@ makedocs(; size_threshold = nothing, ), pages=[ + "Contents" => "contents.md", "Home" => "index.md", "User Manual" => "manual.md", demopage, "Usage Examples" => "examples.md", - "Function Reference" => "reference.md", - "Index" => "function_index.md" + "API Documentation" => "reference.md", + "API Index" => "function_index.md", + "References" => "references.md" ], ) postprocess_cb() # For DemoCards @@ -43,3 +56,5 @@ deploydocs(; repo="github.com/simonp0420/PSSFSS.jl", devbranch = "main" ) + +cd(olddir) \ No newline at end of file diff --git a/docs/notebooks/angular_ss_example.ipynb b/docs/notebooks/angular_ss_example.ipynb index 9df2e194..cadb5d88 100644 --- a/docs/notebooks/angular_ss_example.ipynb +++ b/docs/notebooks/angular_ss_example.ipynb @@ -4,11 +4,7 @@ "cell_type": "markdown", "source": [ "## Angle Selective Surface\n", - "This example is taken from Zhenting Chen, Chao Du, Jie Liu, Di Zhou, and Zhongxiang Shen,\n", - "\"Design Methodology of Dual-Polarized Angle-Selective Surface Based\n", - "on Three-Layer Frequency-Selective Surfaces\", IEEE Trans. Antennas Propagat., Vol. 71, No. 11, November 2023,\n", - "pp. 8704--8713.\n", - "\n", + "This example is taken [chen2023design](@cite).\n", "A three-sheet FSS is designed that is transparent to normally incident plane waves, but strongly\n", "attenuates obliquely incident waves. All three sheets are swastika-shaped, with the outer two\n", "sheets being identical. The shapes are generated in PSSFSS using the manji element." @@ -162,11 +158,11 @@ "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", - "version": "1.11.2" + "version": "1.11.3" }, "kernelspec": { "name": "julia-1.11", - "display_name": "Julia 1.11.2", + "display_name": "Julia 1.11.3", "language": "julia" } }, diff --git a/docs/notebooks/band_pass_filter.ipynb b/docs/notebooks/band_pass_filter.ipynb index df9a2d1b..a03f6537 100644 --- a/docs/notebooks/band_pass_filter.ipynb +++ b/docs/notebooks/band_pass_filter.ipynb @@ -4,20 +4,17 @@ "cell_type": "markdown", "source": [ "## Loaded Cross Band Pass Filter\n", - "This example is originally from Fig. 7.9 of B. Munk, *Frequency Selective Surfaces,\n", - "Theory and Design,* John Wiley and Sons, 2000. The same case was analyzed in L. Li,\n", - "D. H. Werner et al, \"A Model-Based Parameter Estimation Technique for\n", - "Wide-Band Interpolation of Periodic Moment Method Impedance Matrices With Application to\n", - "Genetic Algorithm Optimization of Frequency Selective Surfaces\", *IEEE Trans. AP-S*,\n", - "vol. 54, no. 3, March 2006, pp. 908-924, Fig. 6. Unfortunately, in neither reference\n", + "This example is originally from [munk2000fss; Fig. 7.9](@cite).\n", + "The same case was analyzed in [li2006model](@cite). Unfortunately, in neither reference\n", "are the dimensions of the loaded cross stated, except for the square unit cell\n", "period of 8.61 mm. I estimated the dimensions from the sketch in Fig. 6 of the second\n", "reference. To provide a reliable comparison, I analyzed one-eighth of the structure\n", - "in HFSS, a commercial finite element solver, using all three planes of symmetry\n", - "(using symmetry in the z = constant centerline plane\n", + "in HFSS, a commercial finite element solver, using all three planes of symmetry.\n", + "Exploiting the symmetry in the z = constant centerline plane\n", "required two analyses, once for an H-wall boundary condition, and once for an E-wall. Those\n", - "results were then combined using the method of Reed and Wheeler (even/odd symmetry)). With\n", - "a much reduced computational domain, it was then possible to drive HFSS well into convergence.\n", + "results were then combined using the method of Reed and Wheeler [reed1956method](@cite), i.e.\n", + "using even/odd symmetry. With a much reduced computational domain, it was then possible to\n", + "drive HFSS well into convergence.\n", "\n", "Two identical loaded cross slot-type elements are separated by a 6 mm layer of dielectric\n", "constant 1.9. Outboard of each sheet is a 1.1 cm layer of dielectric constant 1.3.\n", @@ -127,11 +124,11 @@ "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", - "version": "1.11.2" + "version": "1.11.3" }, "kernelspec": { "name": "julia-1.11", - "display_name": "Julia 1.11.2", + "display_name": "Julia 1.11.3", "language": "julia" } }, diff --git a/docs/notebooks/checkerboard_example.ipynb b/docs/notebooks/checkerboard_example.ipynb new file mode 100644 index 00000000..21923fb4 --- /dev/null +++ b/docs/notebooks/checkerboard_example.ipynb @@ -0,0 +1,281 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "source": [ + "## Checkerboard Metasurface\n", + "Checkerboard-style metasurfaces have been exensively studied, for example in [nakata2013plane](@cite)\n", + "and [kempa2010percolation](@cite). They exhibit some very unusual properties, which we will demonstrate here.\n", + "\n", + "We consider a series of checkerboad-like geometries. The square unit cell has dimension $P = 5~\\text{mm}$.\n", + "The side length for a PEC square sheet rotated 45° and perfectly inscribed in the unit cell is $L_0 = P / \\sqrt{2}$.\n", + "We will analyze a series of squares at 1 GHz with side lengths $L = L_0 + δ$. The function `computed_and_plot`\n", + "defined below will take a given value of $δ$ as its input, then perform four distinct analyses based on two\n", + "complementary triangulations. Let's run it first for $δ = -0.5 \\text{mm}$ and look\n", + "at the resulting triangulations:" + ], + "metadata": { + "name": "A slide ", + "slideshow": { + "slide_type": "subslide" + } + } + }, + { + "outputs": [], + "cell_type": "code", + "source": [ + "using PSSFSS, Plots, PrettyTables\n", + "using Printf: @sprintf\n", + "\n", + "units = mm\n", + "P = 5 # period in x and y direction\n", + "L0 = P / √2 # Side length for self-complementary square\n", + "Nx = Ny = 10\n", + "\n", + "function compute_and_plot(δ; P=P, units=units, L0=L0, Nx=Nx, Ny=Ny)\n", + "\n", + " Px = Py = P\n", + " Lx = Ly = L0 - abs(δ)\n", + " orient = 45\n", + " islekwds = (; Lx, Ly, Px, Py, units, Nx, Ny, orient)\n", + " holekwds = (; units, s1=[P,0], s2=[0,P], sides=4,\n", + " a = iszero(δ) ? [P/2] : [(L0-abs(δ)) / √2],\n", + " b=[-2*Nx], ntri=2*Nx*Ny, structuredtri=false)\n", + " if δ ≤ 0\n", + " redj = rectstrip(; class='J', islekwds...)\n", + " redm = rectstrip(; class='M', islekwds...)\n", + " bluej = polyring(; class='J', holekwds...)\n", + " bluem = polyring(; class='M', holekwds...)\n", + " else\n", + " redj = polyring(; class='J', holekwds...)\n", + " redm = polyring(; class='M', holekwds...)\n", + " bluej = rectstrip(; class='J', islekwds...)\n", + " bluem = rectstrip(; class='M', islekwds...)\n", + " end\n", + "\n", + " p1 = plot(redj, unitcell=true, faces=true, fillcolor=:red, title=\"Unit Cell\")\n", + " p2 = plot(redj, faces=true, edges=false, fillcolor=:red, rep=(3,3), title=\"3×3 Array\")\n", + " p3 = plot(bluej, unitcell=true, faces=true, fillcolor=:blue, title=\"Unit Cell\")\n", + " p4 = plot(bluej, faces=true, edges=false, fillcolor=:blue, rep=(3,3), title=\"3×3 Array\")\n", + " pl = plot(p1,p2,p3,p4, layout=(2,2), size=(550,600), plot_title=\" Triangulations for δ = $δ\")\n", + "\n", + " flist = 1 # Analysis frequency in GHz\n", + " steering = (θ=0, ϕ=0) # Normal incidence\n", + " logfile = resultfile = devnull # Suppress creation of output files\n", + " showprogress = false # Suppress screen output\n", + "\n", + " redjres = analyze([Layer(), redj, Layer()], flist, steering; logfile, resultfile, showprogress)\n", + " bluejres = analyze([Layer(), bluej, Layer()], flist, steering; logfile, resultfile, showprogress)\n", + " redmres = analyze([Layer(), redm, Layer()], flist, steering; logfile, resultfile, showprogress)\n", + " bluemres = analyze([Layer(), bluem, Layer()], flist, steering; logfile, resultfile, showprogress)\n", + "\n", + "\n", + " s11rj = extract_result(redjres, @outputs s11(v,v)) |> only\n", + " s21rj = extract_result(redjres, @outputs s21(v,v)) |> only\n", + " s11rm = extract_result(redmres, @outputs s11(v,v)) |> only\n", + " s21rm = extract_result(redmres, @outputs s21(v,v)) |> only\n", + " s11bj = extract_result(bluejres, @outputs s11(v,v)) |> only\n", + " s21bj = extract_result(bluejres, @outputs s21(v,v)) |> only\n", + " s11bm = extract_result(bluemres, @outputs s11(v,v)) |> only\n", + " s21bm = extract_result(bluemres, @outputs s21(v,v)) |> only\n", + "\n", + " return pl, (; s11rj, s21rj, s11rm, s21rm, s11bj, s21bj, s11bm, s21bm)\n", + "end\n", + "\n", + "pl, _ = compute_and_plot(-0.5)\n", + "pl" + ], + "metadata": {}, + "execution_count": null + }, + { + "cell_type": "markdown", + "source": [ + "The function creates a \"red\" triangulation occupying the triangle of side length $L_0 + δ$, and a \"blue\" triangulation,\n", + "occupying the complementary portion of the unit cell. Since for this case the red square\n", + "side length is 0.5 mm shorter than the critical length $L_0$, it lies strictly inside the unit cell. So if we choose to\n", + "use the red triangulation to model electric surface current, then we can consider the red regions to be \"islands\" of\n", + "metal in otherwise empty space. We could also use the blue triangulation to model magnetic surface current, which again would\n", + "lead to the conclusion that the small untriangulated squares are conducting patches or \"islands\" of metalization.\n", + "Either of these two choices, when analyzed with PSSFSS, should yield the same values for computed reflection or transmission\n", + "coefficients (within modeling accuracy).\n", + "\n", + "A different approach would be to choose the red triangulation for representing magnetic surface current, in which case\n", + "the small red squares would represent \"holes\" in an otherwise solid metallic sheet. The same \"hole\" interpretation would\n", + "follow from using the blue triangulation to represent electric surface current. In fact, for this case, the blue region\n", + "in the full plane can be regarded as the union of an infinite number of metallic squares of dimension $L_0 + δ$.\n", + "So positive values of $δ$ can be handled by using $-δ$ and reversing the roles of the red and blue triangulations.\n", + "This is exactly what is done in the function `compute_and_plot` above." + ], + "metadata": {} + }, + { + "cell_type": "markdown", + "source": [ + "We'll now exercise the function for the set of $δ$ values $\\{ -0.2, -0.05, 0, 0.05, 0.2 \\}$, observing both the plotted\n", + "triangulations and the resulting scattering parameter predictions for each of the four modeling choices outlined above." + ], + "metadata": {} + }, + { + "cell_type": "markdown", + "source": [ + "```julia\n", + "δs = [-0.2, -0.05, 0, 0.05, 0.2] # Departure in mm from self-complementary square side length\n", + "s11rj, s21rj, s11rm, s21rm, s11bj, s21bj, s11bm, s21bm = (zeros(ComplexF64, length(δs)) for _ in 1:8)\n", + "for (i, δ) in pairs(δs)\n", + " plt, r = compute_and_plot(δ)\n", + " s11rj[i] = r.s11rj\n", + " s21rj[i] = r.s21rj\n", + " s11rm[i] = r.s11rm\n", + " s21rm[i] = r.s21rm\n", + " s11bj[i] = r.s11bj\n", + " s21bj[i] = r.s21bj\n", + " s11bm[i] = r.s11bm\n", + " s21bm[i] = r.s21bm\n", + " display(plt)\n", + "end\n", + "```" + ], + "metadata": {} + }, + { + "cell_type": "markdown", + "source": [ + ".png)" + ], + "metadata": {} + }, + { + "cell_type": "markdown", + "source": [ + ".png)" + ], + "metadata": {} + }, + { + "cell_type": "markdown", + "source": [ + ".png)" + ], + "metadata": {} + }, + { + "cell_type": "markdown", + "source": [ + ".png)" + ], + "metadata": {} + }, + { + "cell_type": "markdown", + "source": [ + ".png)" + ], + "metadata": {} + }, + { + "cell_type": "markdown", + "source": [ + "Let the letter \"J\" denote use of a triangulation to represent electric surface current and let \"M\" denote magnetic surface current.\n", + "So, e.g. \"Blue M\" means that the blue triangulation is used to represent magnetic current. We can make the following observations about\n", + "the above plots:\n", + "\n", + "1. The red and blue triangulations alway occupy complementary regions of the unit cell. Then for a given choice\n", + " of $δ$, PSSFSS analysis of \"Blue J\" and \"Red M\" should result identical scattering parameters (apart from modeling errors).\n", + " Likewise analysis of \"Red J\" and \"Blue M\" should result in the same scattering parameters.\n", + "2. It is well known from Babinet's principle [nakata2013plane](@cite) [tan2012babinet](@cite) that reflection coefficients\n", + " of thin perforated screens are equal to the negative of the transmission coefficients of the complementary structure, provided the\n", + " structures exhibit sufficient rotational symmetry as is the case here. \"Red J\" and \"Blue J\" form such a\n", + " complementary pair, as do \"Red M\" and \"Blue M\".\n", + "3. All of the $δ$ values are quite small compared to $L_0 \\approx 3.5355 \\text{mm}$. In particular, the plots for\n", + " $δ = \\pm 0.05$ are almost indistinguishable from the plot for $δ = 0$. Therefore, we expect the scattering\n", + " parameters of all of these cases to be nearly equal.\n", + "4. For the case of $δ = 0$, the red and blue regions are both self-complementary. From our earlier considerations, we then\n", + " expect equal reflection and transmission coefficients for any of the choices \"Red J\", \"Blue J\", \"Red M\", or \"Blue M\". In\n", + " fact all four of these cases should yield the same values.\n", + "\n", + "The following code generates a summary table showing how well these expectations are satisfied:" + ], + "metadata": {} + }, + { + "cell_type": "markdown", + "source": [ + "```julia\n", + "mat = hcat(s11rj, s11bm, -s21rm, -s21bj) |> transpose\n", + "row_labels = [\"S₁₁: Red J \", \"S₁₁: Blue M\", \"-S₂₁: Red M \", \"-S₂₁: Blue J\"]\n", + "header = ([\"δ = $δ mm\" for δ in δs], [\"islands or holes\" for δ in δs])\n", + "header[2][findfirst(iszero, δs)] = \"SC (Self-Complementary)\"\n", + "formatters = (v,i,j) -> imag(v) > 0 ? @sprintf(\"%7.4f + %6.4fim\", real(v), imag(v)) :\n", + " @sprintf(\"%7.4f - %6.4fim\", real(v), -imag(v))\n", + "highlighters = Highlighter((data, i, j) -> (j == findfirst(iszero,δs)), crayon\"red bold\")\n", + "pretty_table(mat; header, row_labels, alignment=:c, formatters, highlighters)\n", + "```" + ], + "metadata": {} + }, + { + "outputs": [], + "cell_type": "code", + "source": [ + "for line in eachline(\"./assets/checkerboard_prettytable.data\"); println(line); end #hide" + ], + "metadata": {}, + "execution_count": null + }, + { + "cell_type": "markdown", + "source": [ + "For $δ < 0$ the results seem reasonable. In this regime of electrically small metal islands (for Red J and Blue M),\n", + "one wouldn't expect much reflection, and that is what is observed. For $δ > 0$, the geometry (for Red J and Blue M)\n", + "consists of a metal plate with small holes, so one would expect almost total reflection, as observed. Now consider\n", + "how well our above observations are aligned with the computed results...\n", + "\n", + "From observations 1 and 2, the numbers in any one column should all be approximately equal, which they are\n", + "except for the center $δ = 0$ column. From observation 3 we expect all the entries\n", + "in any one row to be nearly equal, but this is not true at all. In any one row there is a violent jump in\n", + "amplitude at or near $δ = 0$. Finally, from observation 4 we expect all of the numeric entries for $δ = 0$\n", + "to be approximately equal--which they are most definitely not. What is going on here?\n", + "\n", + "As discussed in the previously cited references, the idealized model being analyzed for $δ = 0$ is unphysical.\n", + "As $\\delta \\rightarrow 0$ the response of the structure does not approach a limit, since the violent jump in\n", + "reflection coefficient occurs for arbitrarily small positive and negative excursions from $L_0$. Any physically\n", + "realizable structure must exhibit continuous dependence on its physical parameters. If one attempts\n", + "to approximate this ideal surface in the real world one must use finite thickness and conductivity, so that neighboring\n", + "unit cells will intersect all along the thickness of the metal, rather than at a single point, thus destroying the\n", + "self-complementary property. Similarly, finite losses in any real metal preclude self-complementarity.\n", + "\n", + "The same problem observed here for an ideal model of a self-complementary surface arises when analyzing pixelated\n", + "structures such as shown in the next example, where adjacent metal pixels\n", + "can intersect at only a single point. In that case, using a \"Blue M\" approach is known to agree with measurements." + ], + "metadata": {} + }, + { + "cell_type": "markdown", + "source": [ + "---\n", + "\n", + "*This notebook was generated using [Literate.jl](https://github.com/fredrikekre/Literate.jl).*" + ], + "metadata": {} + } + ], + "nbformat_minor": 3, + "metadata": { + "language_info": { + "file_extension": ".jl", + "mimetype": "application/julia", + "name": "julia", + "version": "1.11.3" + }, + "kernelspec": { + "name": "julia-1.11", + "display_name": "Julia 1.11.3", + "language": "julia" + } + }, + "nbformat": 4 +} diff --git a/docs/notebooks/cpss1.ipynb b/docs/notebooks/cpss1.ipynb index 82769026..9cafb0e6 100644 --- a/docs/notebooks/cpss1.ipynb +++ b/docs/notebooks/cpss1.ipynb @@ -9,10 +9,7 @@ "We'll first look at analyzing a design presented in the literature, and then proceed to optimize another\n", "design using PSSFSS as the analysis engine inside the optimization objective function.\n", "### Sjöberg and Ericsson Design\n", - "This example comes from the paper D. Sjöberg and A. Ericsson, \"A multi layer meander line circular\n", - "polarization selective structure (MLML-CPSS),\" The 8th European Conference on Antennas and Propagation\n", - "(EuCAP 2014), 2014, pp. 464-468, doi: 10.1109/EuCAP.2014.6901792.\n", - "\n", + "This example comes from [sjoberg2014multi](@cite).\n", "The authors describe an ingenious structure consisting of 5 progressively rotated meanderline sheets, which\n", "acts as a circular polarization selective surface: it passes LHCP (almost) without attenuation or\n", "reflection, and reflects RHCP (without changing its sense!) almost without attenuation or transmission." @@ -297,7 +294,7 @@ "the entire unit cell, and the unit cell for sheets 2 and 4 are not square. Since the periodicity of\n", "the sheets in the structure varies from sheet to sheet, higher order Floquet modes common to neighboring\n", "sheets cannot be defined, so we are forced to use only the dominant (0,0) modes which are independent of\n", - "the periodicity. This limitation is removed in a later example.\n", + "the periodicity. This limitation is removed in the next example.\n", "Meanwhile, it is of interest to note that their high-accuracy runs\n", "required 10 hours for CST and 19 hours for COMSOL on large engineering workstations versus about 22\n", "seconds for PSSFSS on my desktop machine." @@ -320,11 +317,11 @@ "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", - "version": "1.11.2" + "version": "1.11.3" }, "kernelspec": { "name": "julia-1.11", - "display_name": "Julia 1.11.2", + "display_name": "Julia 1.11.3", "language": "julia" } }, diff --git a/docs/notebooks/cpss2.ipynb b/docs/notebooks/cpss2.ipynb index cd7e49a5..0a2e24dd 100644 --- a/docs/notebooks/cpss2.ipynb +++ b/docs/notebooks/cpss2.ipynb @@ -4,10 +4,7 @@ "cell_type": "markdown", "source": [ "## Meanderline/Strip-Based CPSS\n", - "This example comes from the same authors as the previous example. The paper is\n", - "A. Ericsson and D. Sjöberg, \"Design and Analysis of a Multilayer Meander Line\n", - "Circular Polarization Selective Structure\", IEEE Trans. Antennas Propagat.,\n", - "Vol. 65, No. 8, Aug 2017, pp. 4089-4101.\n", + "This example comes from [ericsson2017design](@cite), by the same authors as the previous example.\n", "The design is similar to that of the previous example except that here, the two $\\pm 45^\\circ$\n", "rotated meanderlines are replaced with rectangular strips.\n", "This allows us to employ the `diagstrip` element and the `orient` keyword for the\n", @@ -192,7 +189,7 @@ { "cell_type": "markdown", "source": [ - "Layers 3 and 9 were assigned 10 modes each. Layers 5 and 7, being thinner were assigned\n", + "Layers 3 and 9 were assigned 10 modes each. Layers 5 and 7, being thinner, were assigned\n", "18 modes each. The numbers of modes are determined automatically by PSSFSS to ensure\n", "accurate cascading." ], @@ -287,11 +284,11 @@ "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", - "version": "1.11.2" + "version": "1.11.3" }, "kernelspec": { "name": "julia-1.11", - "display_name": "Julia 1.11.2", + "display_name": "Julia 1.11.3", "language": "julia" } }, diff --git a/docs/notebooks/cpss_optimization.ipynb b/docs/notebooks/cpss_optimization.ipynb index ca20d303..43a66ac7 100644 --- a/docs/notebooks/cpss_optimization.ipynb +++ b/docs/notebooks/cpss_optimization.ipynb @@ -229,7 +229,7 @@ "cell_type": "markdown", "source": [ "As hoped for, the performance meets the more stringent design goals over a broader bandwidth than the\n", - "Sjöberg and Ericsson design, presumably because of the greater design flexibility allowed here." + "Sjöberg and Ericsson design of [sjoberg2014multi](@cite), presumably because of the greater design flexibility allowed here." ], "metadata": {} }, @@ -249,11 +249,11 @@ "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", - "version": "1.11.2" + "version": "1.11.3" }, "kernelspec": { "name": "julia-1.11", - "display_name": "Julia 1.11.2", + "display_name": "Julia 1.11.3", "language": "julia" } }, diff --git a/docs/notebooks/cross_on_dielectric_substrate.ipynb b/docs/notebooks/cross_on_dielectric_substrate.ipynb index 87bb3038..aa967863 100644 --- a/docs/notebooks/cross_on_dielectric_substrate.ipynb +++ b/docs/notebooks/cross_on_dielectric_substrate.ipynb @@ -4,12 +4,10 @@ "cell_type": "markdown", "source": [ "## Cross on Dielectric Substrate\n", - "This example is also taken from the paper by Alon S. Barlevy and\n", - "Yahya Rahmat-Samii, \"Fundamental Constraints on the Electrical Characteristics\n", - "of Frequency Selective Surfaces\", **Electromagnetics**, vol. 17, 1997, pp. 41-68.\n", - "This particular example is from Section 3.2, Figures 7 and 8. It also appeared at\n", - "higher resolution in Barlevy's PhD dissertation from which the comparison curves\n", - "were digitized.\n", + "This example is also taken from the paper [barlevy1997fundamental](@cite) by Barlevy and\n", + "Rahmat-Samii. This particular example is from Section 3.2, Figures 7 and 8. It also appeared at\n", + "higher resolution in Barlevy's PhD dissertation [barlevy1998dissertation](@cite)\n", + "from which the comparison curves were digitized.\n", "\n", "We use the `loadedcross` element where we choose `w > L2/2`, so that the Cross\n", "is \"unloaded\", i.e. the center section is filled in with metalization:" @@ -111,7 +109,7 @@ "cell_type": "markdown", "source": [ "The above loop requires about 18 seconds of execution time on my machine.\n", - "Compare PSSFSS results to those digitized from the dissertation figure:" + "Compare PSSFSS results to those digitized from the dissertation:" ], "metadata": {} }, @@ -160,11 +158,11 @@ "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", - "version": "1.11.2" + "version": "1.11.3" }, "kernelspec": { "name": "julia-1.11", - "display_name": "Julia 1.11.2", + "display_name": "Julia 1.11.3", "language": "julia" } }, diff --git a/docs/notebooks/flexible_absorber.ipynb b/docs/notebooks/flexible_absorber.ipynb index ba588c0e..0ce1c357 100644 --- a/docs/notebooks/flexible_absorber.ipynb +++ b/docs/notebooks/flexible_absorber.ipynb @@ -4,9 +4,7 @@ "cell_type": "markdown", "source": [ "## Flexible Absorber\n", - "This example is from Yize Li, et al., \"Ultra-wideband, polarization-insensitive flexible metamaterial\n", - "absorber base on laser printed graphene using equivalent circuit design method,\" Carbon, Vol 212, 2023,\n", - "available for free download from [here](https://doi.org/10.1016/j.carbon.2023.118166).\n", + "This example is from [li2023ultra](@cite).\n", "It uses square and circular resistive FSS elements sandwiched between layers of flexible dielectrics to\n", "realize a reflective absorber (i.e. a \"rabsorber\").\n", "We compare the reflection coefficient magnitude computed by PSSFSS to that digitized\n", @@ -57,13 +55,13 @@ "results = analyze(strata, FGHz, scan, showprogress=false, resultfile=devnull, logfile=devnull)\n", "\n", "s11db = extract_result(results, @outputs s11db(te,te))\n", - "yize = readdlm(\"../src/assets/yize2023_fig2a_s11db_digitized.csv\", ',', Float64, '\\n')\n", + "li = readdlm(\"../src/assets/li2023_fig2a_s11db_digitized.csv\", ',', Float64, '\\n')\n", "\n", "ps11 = plot(title=\"Normal Incidence Reflection Magnitude\",\n", " xlabel=\"Frequency (GHz)\", ylabel=\"20log₁₀|s₁₁|\", xlim=(0,20),\n", " ylim=(-18,0), xtick=0:2:20, ytick=-20:2:0, framestyle=:box)\n", "plot!(ps11, FGHz, s11db, color=:blue, label=\"PSSFSS\", lw=2)\n", - "plot!(ps11, yize[:,1], yize[:,2], color=:red, label=\"Yize et al.\", lw=2)" + "plot!(ps11, li[:,1], li[:,2], color=:red, label=\"Li et al.\", lw=2)" ], "metadata": {}, "execution_count": null @@ -99,11 +97,11 @@ "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", - "version": "1.11.2" + "version": "1.11.3" }, "kernelspec": { "name": "julia-1.11", - "display_name": "Julia 1.11.2", + "display_name": "Julia 1.11.3", "language": "julia" } }, diff --git a/docs/notebooks/fresneltable_creation_example.ipynb b/docs/notebooks/fresneltable_creation_example.ipynb new file mode 100644 index 00000000..96b359ef --- /dev/null +++ b/docs/notebooks/fresneltable_creation_example.ipynb @@ -0,0 +1,65 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "source": [ + "## Fresnel Table Creation\n", + "Here we show how to create an HFSS-compatible Fresnel table using the `res2fresnel`\n", + "function included with PSSFSS. The geometry for this example is a double square loop\n", + "in a 5 cm square unit cell. The code for analyzing this geometry and creating the\n", + "Fresnel table is shown below:\n", + "```Julia\n", + "using PSSFSS\n", + "dwidth = 3mm\n", + "duroid = Layer(epsr=2.2, tandel=0.0009, width=dwidth)\n", + "a = √2 * [1, 2.125]\n", + "b = √2 * [1.5, 3.125]\n", + "units = mm\n", + "sides = 4\n", + "orient = 45\n", + "sheet = polyring(;a, b, units, sides, orient, s1=[8,0], s2=[0,8], ntri=2000)\n", + "strata = [Layer(), duroid, sheet, duroid, Layer(width=-2*dwidth)]\n", + "steering = (θ=0:5:45, ϕ=0)\n", + "FGHz = 10:2:20\n", + "results = analyze(strata, FGHz, steering; resultfile=\"double_square_loop.res\")\n", + "res2fresnel(results, \"double_square_loop.rttbl\")\n", + "# Alternatively: res2fresnel(\"double_square_loop.res\", \"double_square_loop.rttbl\")\n", + "```\n", + "Note that the final layer specifies a width that is the negative of the sum of the remaining layers' widths.\n", + "This has the effect of moving the output phase reference plane to coincide with the input phase reference plane.\n", + "Please see the documentation for `res2fresnel` for discussion of the setup requirements and restrictions\n", + "necessary for creating a valid Fresnel table." + ], + "metadata": { + "name": "A slide ", + "slideshow": { + "slide_type": "Slide" + } + } + }, + { + "cell_type": "markdown", + "source": [ + "---\n", + "\n", + "*This notebook was generated using [Literate.jl](https://github.com/fredrikekre/Literate.jl).*" + ], + "metadata": {} + } + ], + "nbformat_minor": 3, + "metadata": { + "language_info": { + "file_extension": ".jl", + "mimetype": "application/julia", + "name": "julia", + "version": "1.11.3" + }, + "kernelspec": { + "name": "julia-1.11", + "display_name": "Julia 1.11.3", + "language": "julia" + } + }, + "nbformat": 4 +} diff --git a/docs/notebooks/manual.ipynb b/docs/notebooks/manual.ipynb index 35703c72..e61fbde8 100644 --- a/docs/notebooks/manual.ipynb +++ b/docs/notebooks/manual.ipynb @@ -61,9 +61,7 @@ "source": [ "### A Quick Example\n", "Here is an example run of PSSFSS for a 4-sheet meanderline polarizer design\n", - "from one of the first papers on the subject: T. L. Blackney, J. R. Burnett, and\n", - "S. B. Cohn, “A design method for meander-line circular polarizers”\n", - "presented at 22nd Annual Antenna Symposium, Oct. 1972.\n", + "from [blbc:72](@cite), one of the first papers on the subject.\n", "Detailed explanations of the code are omitted for now." ], "metadata": { @@ -138,7 +136,7 @@ "\n", "\n", "\n", - "The PSSFSS run took about 24 seconds on my machine for 4 meanderline PSS sheets\n", + "The PSSFSS run took about 14 seconds on my machine for 4 meanderline PSS sheets\n", "analyzed at 61 frequencies.\n", "\n", "I hope this example whetted your appetite to learn\n", @@ -242,7 +240,7 @@ { "cell_type": "markdown", "source": [ - "It is almost essential that PSSFSS users also install the [Plots](https://github.com/JuliaPlots/Plots.jl)\n", + "It is essential that PSSFSS users also install the [Plots](https://github.com/JuliaPlots/Plots.jl)\n", "package. This will allow easy visualization of the FSS/PSS element triangulations produced by PSSFSS, in\n", "addition to providing a convenient means to plot analysis results." ], @@ -364,8 +362,12 @@ "cell_type": "markdown", "source": [ "## Strata\n", - "### Layer\n", - "Dielectric layers are created with the `Layer` function:" + "The geometry to be analyzed is passed as the first argument of the `analyze` function in the form\n", + "of a Julia `Vector` (i.e., a one-dimensional array). This vector must contain at least two `Layer`\n", + "objects representing dielectric layers, and zero or more `RWGSheet` instances, representing the zero-thickness\n", + "PSS or FSS sheets located between the dielectric layers. These are described below.\n", + "### Dielectric Layers\n", + "Dielectric layers are created with the `Layer` constructor function:" ], "metadata": { "name": "A slide ", @@ -459,7 +461,7 @@ { "cell_type": "markdown", "source": [ - "The call to `rectstrip` above creates a `RWGSheet` object for a rectangular strip\n", + "The call to `rectstrip` above creates an `RWGSheet` object for a rectangular strip\n", "of dimensions 0.5 cm in the x and y directions, lying in a square unit cell of dimension\n", "1 cm. The triangulation uses 10 edges in the x and y directions (`Nx` and `Ny`).\n", "\n", @@ -519,8 +521,7 @@ "default value, suitable for the so-called \"oxide side\" of a planar conductor, or `:rayleigh` (suitable for the\n", "\"foil side\" of the conductor, i.e. the side bonded to the dielectric substrate). Together, the conductivity,\n", "surface roughness, and roughness distribution type are used by PSSFSS internally to compute a frequency-dependent surface\n", - "impedance, using the so-called Gradient Model, as described in D. N. Grujić, “Simple and Accurate Approximation of\n", - "Rough Conductor Surface Impedance,” IEEE Trans. Microwave Theory Tech., vol. 70, no. 4, pp. 2053-2059, April 2022.\n", + "impedance, using the so-called Gradient Model, as described in [grujic2021simple](@cite).\n", "Obviously, only one of `Zsheet` and `sigma` may be specified as keyword arguments for a given `RWGSheet`." ], "metadata": { @@ -954,7 +955,7 @@ "source": [ "It is often desired to use a set of polarization basis vectors other than TE/TM to define field\n", "coefficients. PSSFSS supports in addition to TE/TM the use of H/V for horizontal/vertical components\n", - "(in the [Ludwig 3](https://ieeexplore.ieee.org/document/1140406) sense), and L/R for left-hand circular\n", + "(in the Ludwig 3 [ludw:73](@cite) sense), and L/R for left-hand circular\n", "and right-hand circular polarization (in the\n", "[IEEE sense](https://en.wikipedia.org/wiki/Circular_polarization#Uses_of_the_two_conventions))." ], @@ -1114,7 +1115,8 @@ "source": [ "#### Using the `analyze` Function Return Value\n", "The variable returned from `analyze` can be used in a similar way to the result file to obtain any quantity\n", - "available to the `@outputs` macro, but for this purpose we use the `extract_result` function:\n", + "available to the `@outputs` macro. For this purpose we pass the return value of the `analyze` function\n", + "as the first argument to the `extract_result` function:\n", "\n", "```julia\n", "julia> results = analyze(strata, flist, steering, outlist=outputs);\n", @@ -1141,10 +1143,40 @@ { "cell_type": "markdown", "source": [ - "### Saving Analysis Results As TEP Files\n", + "## Exporting Analysis Results\n", + "Besides the ability to export results to CSV files as described here,\n", + "one can export results to TICRA-compatible TEP files or HFSS-compatible Fresnel tables, as described in the following\n", + "subsections.\n", + "%% A slide [markdown] {\"slideshow\": {\"slide_type\": \"subslide\"}}\n", + "### Exporting Analysis Results to TEP Files\n", "The vector of `Result` objects returned by the `analyze` function (or the corresponding\n", "result file) can be converted to a TICRA-compatible TEP (tabulated electrical properties) file\n", - "using the function `res2tep`. The code needed to do so is shown in this example." + "using the function `res2tep`. TEP files contain the equivalent of the full\n", + "4×4 scattering matrix computed by PSSFSS. Therefore, there is no requirement/limitation for unit cell geometry\n", + "to exhibit any particular symmetry when results are to be exported to a TEP file. There are, however, requirements\n", + "on the choice of steering angles to be analyzed. For details, see the documentation of `res2tep`.\n", + "Sample code for creating a TEP file is shown in this example." + ], + "metadata": { + "name": "A slide ", + "slideshow": { + "slide_type": "subslide" + } + } + }, + { + "cell_type": "markdown", + "source": [ + "### Exporting Analysis Results to Fresnel Tables\n", + "The vector of `Result` objects returned by the `analyze` function (or the corresponding\n", + "result file) can be converted to an HFSS SBR+-compatible Fresnel table\n", + "using the function `res2fresnel`. Fresnel tables contain reflection\n", + "and transmission coefficients for only \"front\" surface incidence, and only for co-polarized\n", + "(TE → TE and TM → TM) scattering. In addition, coefficients for only one azimuth angle ϕ are\n", + "retained in such a table. These limitations impose restrictions on the type of structure that\n", + "can be analyzed for use in generating a valid Fresnel table. For details, see the documentation\n", + "of `res2fresnel`. Sample code for creating a Fresnel table is shown in\n", + "this example." ], "metadata": { "name": "A slide ", @@ -1169,11 +1201,11 @@ "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", - "version": "1.11.2" + "version": "1.11.3" }, "kernelspec": { "name": "julia-1.11", - "display_name": "Julia 1.11.2", + "display_name": "Julia 1.11.3", "language": "julia" } }, diff --git a/docs/notebooks/reflectarray_example.ipynb b/docs/notebooks/reflectarray_example.ipynb index 0657dfa4..bc2274bf 100644 --- a/docs/notebooks/reflectarray_example.ipynb +++ b/docs/notebooks/reflectarray_example.ipynb @@ -4,14 +4,10 @@ "cell_type": "markdown", "source": [ "## Reflectarray Element\n", - "This example is taken from Figure 6 of\n", - "Li, Jiao, and Zhao: \"A Novel Microstrip Rectangular-Patch/Ring-\n", - "Combination Reflectarray Element and Its Application\", **IEEE Antennas and Wireless Propagation Letters**,\n", - "VOL. 8, 2009, pp. 1119-1112.\n", - "\n", + "This example is taken from [li2009novel; Figure 6](@cite).\n", "It generates the so-called \"S-curve\" for reflection phase of a reflectarray element. The element\n", "consists of a square patch in a square ring, separated from a ground plane by two dielectric layers.\n", - "Reflection phase is plotted versus the `L2` parameter as defined by Li, et al., which characterizes the\n", + "Reflection phase is plotted versus the `L2` parameter as defined in [li2009novel](@cite), which characterizes the\n", "overall size of the element." ], "metadata": { @@ -273,11 +269,11 @@ "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", - "version": "1.11.2" + "version": "1.11.3" }, "kernelspec": { "name": "julia-1.11", - "display_name": "Julia 1.11.2", + "display_name": "Julia 1.11.3", "language": "julia" } }, diff --git a/docs/notebooks/resistive_square_patch.ipynb b/docs/notebooks/resistive_square_patch.ipynb index f2f47ca6..993cf1d1 100644 --- a/docs/notebooks/resistive_square_patch.ipynb +++ b/docs/notebooks/resistive_square_patch.ipynb @@ -6,10 +6,8 @@ "## Resistive Square Patch\n", "This example will demonstrate the ability of PSSFSS to accurately model finite\n", "conductivity of FSS metalization. It consists of a square finitely conducting\n", - "patch in a square lattice. It is taken from a paper by Alon S. Barlevy and\n", - "Yahya Rahmat-Samii,\n", - "\"Fundamental Constraints on the Electrical Characteristics of Frequency Selective\n", - "Surfaces\", **Electromagnetics**, vol. 17, 1997, pp. 41-68. This particular example\n", + "patch in a square lattice. It is taken from a paper by Barlevy and\n", + "Rahmat-Samii [barlevy1997fundamental](@cite). This particular example\n", "is from Section 3.2, Figures 7 and 8. We will compare PSSFSS results to those digitized\n", "from the cited figures.\n", "\n", @@ -115,11 +113,11 @@ "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", - "version": "1.11.2" + "version": "1.11.3" }, "kernelspec": { "name": "julia-1.11", - "display_name": "Julia 1.11.2", + "display_name": "Julia 1.11.3", "language": "julia" } }, diff --git a/docs/notebooks/splitring_cpss.ipynb b/docs/notebooks/splitring_cpss.ipynb index f0fa51f2..04da8f1d 100644 --- a/docs/notebooks/splitring_cpss.ipynb +++ b/docs/notebooks/splitring_cpss.ipynb @@ -5,10 +5,7 @@ "source": [ "## Split Ring-Based CPSS\n", "This circular polarization selective surface (CPSS) example comes from the paper\n", - "L.-X. Wu, K. Chen, T. Jiang, J. Zhao and Y. Feng, \"Circular-Polarization-Selective\n", - "Metasurface and Its Applications to Transmit-Reflect-Array Antenna and Bidirectional\n", - "Antenna,\" in IEEE Trans. Antennas and Propag., vol. 70, no. 11, pp. 10207-10217,\n", - "Nov. 2022, doi: 10.1109/TAP.2022.3191213.\n", + "[wu2022circular](@cite) by Wu et al.\n", "The design consists of three sequentially rotated split rings separated by dielectric\n", "layers. Since the unit cells for all three rings are identical, PSSFSS can rigorously\n", "account for multiple scattering between the individual sheets using multiple\n", @@ -226,11 +223,11 @@ "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", - "version": "1.11.2" + "version": "1.11.3" }, "kernelspec": { "name": "julia-1.11", - "display_name": "Julia 1.11.2", + "display_name": "Julia 1.11.3", "language": "julia" } }, diff --git a/docs/notebooks/splitringexample.ipynb b/docs/notebooks/splitringexample.ipynb index 761695ab..3b48e9a5 100644 --- a/docs/notebooks/splitringexample.ipynb +++ b/docs/notebooks/splitringexample.ipynb @@ -4,11 +4,7 @@ "cell_type": "markdown", "source": [ "## Split-Ring Resonator\n", - "This example is taken from Figure 3 of\n", - "Fabian‐Gongora, et al, \"Independently Tunable Closely Spaced Triband\n", - "Frequency Selective Surface Unit Cell Using the Third Resonant Mode of Split Ring Slots\",\n", - "IEEE Access, 8/3/2021, Digital Object Identifier 10.1109/ACCESS.2021.3100325.\n", - "\n", + "This example is taken from [fabian2021independently; Figure 3](@cite).\n", "It consists of three concentric split rings with gaps sequentially rotated by 180° situated\n", "on a thin dielectric slab." ], @@ -111,11 +107,11 @@ "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", - "version": "1.11.2" + "version": "1.11.3" }, "kernelspec": { "name": "julia-1.11", - "display_name": "Julia 1.11.2", + "display_name": "Julia 1.11.3", "language": "julia" } }, diff --git a/docs/notebooks/square_loop_absorber.ipynb b/docs/notebooks/square_loop_absorber.ipynb index 564f221b..6d66381a 100644 --- a/docs/notebooks/square_loop_absorber.ipynb +++ b/docs/notebooks/square_loop_absorber.ipynb @@ -4,9 +4,8 @@ "cell_type": "markdown", "source": [ "## Square Loop Absorber\n", - "This example is from Figure 7 of Costa and Monorchio: \"A frequency selective\n", - "radome with wideband absorbing properties\", *IEEE Trans. AP-S*,\n", - "Vol. 60, no. 6, June 2012, pp. 2740--2747. It shows how one can use the `polyring`\n", + "This example, from Figure 7 of [costa2012frequency](@cite),\n", + "shows how one can use the `polyring`\n", "function to model square loop elements. Three different designs are examined\n", "that employ different loop thicknesses and different values of sheet resistance.\n", "We compare the reflection coefficient magnitudes computed by PSSFSS with those digitized\n", @@ -155,11 +154,11 @@ "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", - "version": "1.11.2" + "version": "1.11.3" }, "kernelspec": { "name": "julia-1.11", - "display_name": "Julia 1.11.2", + "display_name": "Julia 1.11.3", "language": "julia" } }, diff --git a/docs/notebooks/symmetric_strip.ipynb b/docs/notebooks/symmetric_strip.ipynb index 2d4b3b0b..d51c3b3e 100644 --- a/docs/notebooks/symmetric_strip.ipynb +++ b/docs/notebooks/symmetric_strip.ipynb @@ -14,10 +14,10 @@ "The dashed lines show two possible choices for the unit cell location: \"J\" for a formulation in terms of electric\n", "surface currents, and \"M\" for magnetic surface currents.\n", "\n", - "For normal incidence there is a closed-form solution due to Weinstein,\n", - "but for a more recent reference one can find the solution in Problem 10.6 of R. E. Collin,\n", - "*Field Theory of Guided Waves, Second Ed.*,\n", - "IEEE Press, 1991. Here is the code for computing the exact solution:" + "For normal incidence there is a closed-form solution due to Weinstein [weinstein1969theory](@cite),\n", + "but there are more recent (and accessible) references in which one can also find the\n", + "solution. For example [collin1991field; Problem 10.6](@cite) and [daniele1990diffraction](@cite).\n", + "Here is the code for computing the exact solution, based on [collin1991field](@cite):" ], "metadata": { "name": "A slide ", @@ -65,8 +65,8 @@ { "cell_type": "markdown", "source": [ - " Note that using the extension of\n", - "[Babinet's Principle for electromagnetic fields](http://kirkmcd.princeton.edu/examples/babinet.pdf)\n", + "Note that using the extension of\n", + "Babinet's Principle for electromagnetic fields [tan2012babinet](@cite), [nakata2013plane](@cite),\n", "this also provides the solution (upon appropriate interchange and sign change of the coefficients) for\n", "the case where the incident wave polarization is parallel to the direction of the strips." ], @@ -356,11 +356,11 @@ "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", - "version": "1.11.2" + "version": "1.11.3" }, "kernelspec": { "name": "julia-1.11", - "display_name": "Julia 1.11.2", + "display_name": "Julia 1.11.3", "language": "julia" } }, diff --git a/docs/notebooks/sympixels_example.ipynb b/docs/notebooks/sympixels_example.ipynb new file mode 100644 index 00000000..bcd6c53b --- /dev/null +++ b/docs/notebooks/sympixels_example.ipynb @@ -0,0 +1,171 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "source": [ + "## Pixelated Band Pass Filter\n", + "This example of from [hong2021design; Fig. 4a](@cite). It consists of a single-layer\n", + "band-pass filter consisting of many square metallic \"pixels\", optimized using a special binary\n", + "optimizer described in the cited paper.\n", + "\n", + "The geometry for the pixelated filter is created using the `sympixels` function as\n", + "show below." + ], + "metadata": { + "name": "A slide ", + "slideshow": { + "slide_type": "Slide" + } + } + }, + { + "outputs": [], + "cell_type": "code", + "source": [ + "using PSSFSS, Plots\n", + "P = 5.4\n", + "units = mm\n", + "nrim = 1\n", + "halfnint = 26\n", + "patternvec = Bool[\n", + " 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 1, 1, 0, 1, 1, 0, 0, 1, 0, 0, 1, 0, 0, 0, 1,\n", + " 0, 1, 0, 1, 0, 1, 1, 1, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 1, 1, 0, 1, 0,\n", + " 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 1, 1, 0, 0, 1, 0, 1, 1, 1, 1, 0, 0, 1,\n", + " 0, 1, 0, 0, 0, 1, 0, 1, 1, 1, 1, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0,\n", + " 1, 1, 1, 0, 0, 1, 0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 0, 1, 1, 0, 0, 1,\n", + " 1, 0, 1, 1, 1, 0, 1, 0, 1, 0, 1, 1, 1, 0, 0, 1, 1, 0, 1, 0, 1,\n", + " 0, 0, 0, 1, 0, 0, 1, 0, 0, 1, 1, 1, 0, 1, 0, 0, 1, 0, 0, 0,\n", + " 0, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 1, 0, 1, 1, 0, 0, 1, 0,\n", + " 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 1, 1,\n", + " 1, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 1, 0, 0, 1, 1,\n", + " 0, 1, 0, 1, 0, 0, 1, 0, 0, 1, 0, 1, 0, 0, 1, 1,\n", + " 1, 1, 0, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 0, 0,\n", + " 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0,\n", + " 1, 1, 1, 0, 1, 0, 0, 1, 1, 1, 1, 0, 1,\n", + " 1, 1, 1, 1, 1, 0, 0, 1, 1, 0, 1, 1,\n", + " 0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 1,\n", + " 1, 0, 1, 1, 1, 1, 0, 1, 1, 1,\n", + " 0, 1, 0, 0, 0, 1, 0, 0, 0,\n", + " 0, 0, 1, 1, 0, 0, 1, 1,\n", + " 0, 0, 0, 0, 1, 1, 0,\n", + " 1, 0, 0, 1, 0, 0,\n", + " 0, 0, 1, 1, 1,\n", + " 1, 0, 0, 0,\n", + " 0, 0, 1,\n", + " 0, 0,\n", + " 0]\n", + "\n", + "\n", + "sheet1 = sympixels(; P, units, nrim, halfnint, patternvec, pdiv=1, class='M')\n", + "plot(sheet1, size=(450,450), unitcell=true)" + ], + "metadata": {}, + "execution_count": null + }, + { + "cell_type": "markdown", + "source": [ + "As described in the documentation for `sympixels`, the `1` entries in `patternvec` above\n", + "denote the locations of the metallized pixels. However, by setting `class = 'M'` above we choose to\n", + "triangulate the empty pixel locations, i.e. those that are not occupied by metal. As discussed in\n", + "Checkerboard Metasurface, this choice is necessary to avoid spurious results when analyzing\n", + "a structure of this type.\n", + "\n", + "Note that the above triangulation is created by forming a square for each triangulated pixel and then\n", + "adding a single diagonal. A finer triangulation is required for an accurate analysis result. This can\n", + "be accomplished by increasing `pdiv` to a value greater than `1`. Here is the triangulation that results\n", + "from `pdiv = 2`:" + ], + "metadata": {} + }, + { + "outputs": [], + "cell_type": "code", + "source": [ + "sheet2 = sympixels(; P, units, nrim, halfnint, patternvec, pdiv=2, class='M')\n", + "plot(sheet2, lw=0.5, size=(450,450), unitcell=true)" + ], + "metadata": {}, + "execution_count": null + }, + { + "cell_type": "markdown", + "source": [ + "Now each pixel has been divided into an array of 2×2 squares, each of which receives a diagonal to\n", + "form triangles. The code for analyzing the `pdiv = 1` case would look like this:\n", + "```Julia\n", + "strata = [Layer(), sheet1, Layer(epsr=3.28, tandel=0.007, width=0.05mm), Layer()]\n", + "steering = (θ=0, ϕ=0)\n", + "flist = range(start=20, stop=36, length=401)\n", + "logfile = \"bp_filter_pdiv1.log\"\n", + "resultfile = \"bp_filter_pdiv1.res\"\n", + "results = analyze(strata, flist, steering; logfile, resultfile)\n", + "s21db = extract_result(results, @outputs s21db(te,te))\n", + "```" + ], + "metadata": {} + }, + { + "cell_type": "markdown", + "source": [ + "Accurately analyzing a structure like this consisting of a huge number of pixels can be become\n", + "very expensive (in terms of memory usage and CPU time) very quickly. For example, the `pdiv = 1`\n", + "case generated 3056 basis functions and took about 30 seconds to analyze on a 32 GByte machine.\n", + "But the `pdiv = 2` case generated 15,088 unknowns and required 683 seconds of execution time. It\n", + "also failed on the 32 GByte machine due to lack of memory, but was run successfully on a 64 GByte\n", + "machine." + ], + "metadata": {} + }, + { + "cell_type": "markdown", + "source": [ + "Because the analysis in [hong2021design](@cite) used nonzero metal thickness, I ran the zero-thickness\n", + "geometry in HFSS to observe the significance of the thickness. Here is a plot comparing the digitized\n", + "results from [hong2021design; Fig. 4a](@cite), with the zero-thickness HFSS analysis and the two\n", + "PSSFSS analyses:" + ], + "metadata": {} + }, + { + "cell_type": "markdown", + "source": [ + "" + ], + "metadata": {} + }, + { + "cell_type": "markdown", + "source": [ + "As seen above, the finite metallization thickness does not have a large effect on the transmission trace.\n", + "Using `pdiv = 2` produces much better agreement with HFSS, but it appears that an even larger value would\n", + "be required to achieve full convergence." + ], + "metadata": {} + }, + { + "cell_type": "markdown", + "source": [ + "---\n", + "\n", + "*This notebook was generated using [Literate.jl](https://github.com/fredrikekre/Literate.jl).*" + ], + "metadata": {} + } + ], + "nbformat_minor": 3, + "metadata": { + "language_info": { + "file_extension": ".jl", + "mimetype": "application/julia", + "name": "julia", + "version": "1.11.3" + }, + "kernelspec": { + "name": "julia-1.11", + "display_name": "Julia 1.11.3", + "language": "julia" + } + }, + "nbformat": 4 +} diff --git a/docs/notebooks/tepfile_creation_example.ipynb b/docs/notebooks/tepfile_creation_example.ipynb index d77c34be..bc373c09 100644 --- a/docs/notebooks/tepfile_creation_example.ipynb +++ b/docs/notebooks/tepfile_creation_example.ipynb @@ -4,9 +4,9 @@ "cell_type": "markdown", "source": [ "## TEP File Creation\n", - "Here we show how to create a TICRA-compatible TEP (tabulated electrical properties) file using PSSFSS.\n", - "The geometry for this example is a rectangular copper strip measuring 4 cm × 0.2 cm in a 5 cm square unit cell.\n", - "The code for analyzing this geometry and creating the TEP file is shown below:\n", + "Here we show how to create a TICRA-compatible TEP (tabulated electrical properties) file using the `res2tep`\n", + "function included with PSSFSS. The geometry for this example is a rectangular copper strip measuring\n", + "4 cm × 0.2 cm in a 5 cm square unit cell. The code for analyzing this geometry and creating the TEP file is shown below:\n", "```Julia\n", "using PSSFSS\n", "FGHz = 3.0\n", @@ -18,11 +18,14 @@ "sheet = rectstrip(; Px, Py, Lx, Ly, Nx, Ny, units=cm, sigma=5.7e8)\n", "steering = (θ=0:5:70, ϕ=0:15:345)\n", "strata = [Layer(), sheet, Layer()]\n", - "results = analyze(strata, FGHz, steering)\n", + "resultfile = \"strip.res\"\n", + "results = analyze(strata, FGHz, steering; resultfile)\n", "res2tep(results, \"dipole_pssfss.tep\"; name = \"dipole\", class = \"pssfss\")\n", + "# Alternatively: res2tep(resultfile, \"dipole_pssfss.tep\"; name = \"dipole\", class = \"pssfss\")\n", "```\n", "Note that a very fine discretization has been specified, resulting in 1560 triangles. There are also a large\n", "number of steering angles requested (15 × 24 = 360). This analysis required about 155 seconds on my machine.\n", + "Please see the documentation for `res2tep` for requirements on setting up analysis scan angles.\n", "\n", "For comparison, the same geometry was analyzed using the QUPES program of Ticra Tools Student Edition\n", "2024 (hence the specification for `sigma` above, which is the default conductivity used by QUPES). For the\n", @@ -55,11 +58,11 @@ "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", - "version": "1.11.2" + "version": "1.11.3" }, "kernelspec": { "name": "julia-1.11", - "display_name": "Julia 1.11.2", + "display_name": "Julia 1.11.3", "language": "julia" } }, diff --git a/docs/src/assets/checkerboard_delta(-0.05).png b/docs/src/assets/checkerboard_delta(-0.05).png new file mode 100644 index 00000000..2ab63334 Binary files /dev/null and b/docs/src/assets/checkerboard_delta(-0.05).png differ diff --git a/docs/src/assets/checkerboard_delta(-0.2).png b/docs/src/assets/checkerboard_delta(-0.2).png new file mode 100644 index 00000000..6d196af4 Binary files /dev/null and b/docs/src/assets/checkerboard_delta(-0.2).png differ diff --git a/docs/src/assets/checkerboard_delta(0.0).png b/docs/src/assets/checkerboard_delta(0.0).png new file mode 100644 index 00000000..b44a05e8 Binary files /dev/null and b/docs/src/assets/checkerboard_delta(0.0).png differ diff --git a/docs/src/assets/checkerboard_delta(0.05).png b/docs/src/assets/checkerboard_delta(0.05).png new file mode 100644 index 00000000..578e1760 Binary files /dev/null and b/docs/src/assets/checkerboard_delta(0.05).png differ diff --git a/docs/src/assets/checkerboard_delta(0.2).png b/docs/src/assets/checkerboard_delta(0.2).png new file mode 100644 index 00000000..1a473638 Binary files /dev/null and b/docs/src/assets/checkerboard_delta(0.2).png differ diff --git a/docs/src/assets/checkerboard_prettytable.data b/docs/src/assets/checkerboard_prettytable.data new file mode 100644 index 00000000..f9d9c70b --- /dev/null +++ b/docs/src/assets/checkerboard_prettytable.data @@ -0,0 +1,9 @@ +┌─────────────────────────────────────┬────────────────────┬────────────────────┬─────────────────────────┬────────────────────┬────────────────────┐ +│ │ δ = -0.2 mm │ δ = -0.05 mm │ δ = 0.0 mm │ δ = 0.05 mm │ δ = 0.2 mm │ +│ │ islands or holes │ islands or holes │ SC (Self-Complementary) │ islands or holes │ islands or holes │ +├─────────────────────────────────────┼────────────────────┼────────────────────┼─────────────────────────┼────────────────────┼────────────────────┤ +│ S₁₁: Red J (islands → SC → holes) │ -0.0003 - 0.0178im │ -0.0005 - 0.0213im │ -0.0005 - 0.0227im │ -0.9993 + 0.0268im │ -0.9995 + 0.0223im │ +│ S₁₁: Blue M (islands → SC → holes) │ -0.0005 - 0.0223im │ -0.0007 - 0.0268im │ -0.9994 + 0.0240im │ -0.9995 + 0.0218im │ -0.9997 + 0.0181im │ +│ -S₂₁: Red M (holes → SC → islands) │ -0.0003 - 0.0181im │ -0.0005 - 0.0218im │ -0.0005 - 0.0234im │ -0.9993 + 0.0268im │ -0.9995 + 0.0223im │ +│ -S₂₁: Blue J (holes → SC → islands) │ -0.0005 - 0.0223im │ -0.0007 - 0.0268im │ -0.9994 + 0.0240im │ -0.9995 + 0.0213im │ -0.9997 + 0.0178im │ +└─────────────────────────────────────┴────────────────────┴────────────────────┴─────────────────────────┴────────────────────┴────────────────────┘ diff --git a/docs/src/assets/citations.css b/docs/src/assets/citations.css new file mode 100644 index 00000000..53f1937c --- /dev/null +++ b/docs/src/assets/citations.css @@ -0,0 +1,17 @@ +.citation dl { + display: grid; + grid-template-columns: max-content auto; } + .citation dt { + grid-column-start: 1; } + .citation dd { + grid-column-start: 2; + margin-bottom: 0.75em; } + .citation ul { + padding: 0 0 2.25em 0; + margin: 0; + list-style: none !important;} + .citation ul li { + text-indent: -2.25em; + margin: 0.33em 0.5em 0.5em 2.25em;} + .citation ol li { + padding-left:0.75em;} diff --git a/docs/src/assets/yize2023_fig2a_s11db_digitized.csv b/docs/src/assets/li2023_fig2a_s11db_digitized.csv similarity index 100% rename from docs/src/assets/yize2023_fig2a_s11db_digitized.csv rename to docs/src/assets/li2023_fig2a_s11db_digitized.csv diff --git a/docs/src/assets/sympixel_with_irzone_numbering.png b/docs/src/assets/sympixel_with_irzone_numbering.png new file mode 100644 index 00000000..9451c9c0 Binary files /dev/null and b/docs/src/assets/sympixel_with_irzone_numbering.png differ diff --git a/docs/src/assets/sympixels_fig4a_comparison.svg b/docs/src/assets/sympixels_fig4a_comparison.svg new file mode 100644 index 00000000..66defb82 --- /dev/null +++ b/docs/src/assets/sympixels_fig4a_comparison.svg @@ -0,0 +1,198 @@ + + diff --git a/docs/src/contents.md b/docs/src/contents.md new file mode 100644 index 00000000..9ed71816 --- /dev/null +++ b/docs/src/contents.md @@ -0,0 +1,8 @@ + +# Contents + +```@contents +Depth=5 +Pages = ["index.md", "manual.md", "examples.md", "reference.md"] +``` + diff --git a/docs/src/examples.md b/docs/src/examples.md index 66332083..fb0e24b6 100644 --- a/docs/src/examples.md +++ b/docs/src/examples.md @@ -19,10 +19,10 @@ The grating lies in the ``z=0`` plane with free space on both sides. The shaded The dashed lines show two possible choices for the unit cell location: "J" for a formulation in terms of electric surface currents, and "M" for magnetic surface currents. -For normal incidence there is a closed-form solution due to Weinstein, -but for a more recent reference one can find the solution in Problem 10.6 of R. E. Collin, -*Field Theory of Guided Waves, Second Ed.*, -IEEE Press, 1991. Here is the code for computing the exact solution: +For normal incidence there is a closed-form solution due to Weinstein [weinstein1969theory](@cite), +but there are more recent (and accessible) references in which one can also find the +solution. For example [collin1991field; Problem 10.6](@cite) and [daniele1990diffraction](@cite). +Here is the code for computing the exact solution, based on [collin1991field](@cite): ````@example symmetric_strip """ @@ -55,8 +55,8 @@ function grating(kP; nterms=30) end ```` - Note that using the extension of -[Babinet's Principle for electromagnetic fields](http://kirkmcd.princeton.edu/examples/babinet.pdf) +Note that using the extension of +Babinet's Principle for electromagnetic fields [tan2012babinet](@cite), [nakata2013plane](@cite), this also provides the solution (upon appropriate interchange and sign change of the coefficients) for the case where the incident wave polarization is parallel to the direction of the strips. @@ -260,10 +260,8 @@ EditURL = "../literate/resistive_square_patch.jl" ## Resistive Square Patch This example will demonstrate the ability of PSSFSS to accurately model finite conductivity of FSS metalization. It consists of a square finitely conducting -patch in a square lattice. It is taken from a paper by Alon S. Barlevy and -Yahya Rahmat-Samii, -"Fundamental Constraints on the Electrical Characteristics of Frequency Selective -Surfaces", **Electromagnetics**, vol. 17, 1997, pp. 41-68. This particular example +patch in a square lattice. It is taken from a paper by Barlevy and +Rahmat-Samii [barlevy1997fundamental](@cite). This particular example is from Section 3.2, Figures 7 and 8. We will compare PSSFSS results to those digitized from the cited figures. @@ -327,12 +325,10 @@ EditURL = "../literate/cross_on_dielectric_substrate.jl" ``` ## Cross on Dielectric Substrate -This example is also taken from the paper by Alon S. Barlevy and -Yahya Rahmat-Samii, "Fundamental Constraints on the Electrical Characteristics -of Frequency Selective Surfaces", **Electromagnetics**, vol. 17, 1997, pp. 41-68. -This particular example is from Section 3.2, Figures 7 and 8. It also appeared at -higher resolution in Barlevy's PhD dissertation from which the comparison curves -were digitized. +This example is also taken from the paper [barlevy1997fundamental](@cite) by Barlevy and +Rahmat-Samii. This particular example is from Section 3.2, Figures 7 and 8. It also appeared at +higher resolution in Barlevy's PhD dissertation [barlevy1998dissertation](@cite) +from which the comparison curves were digitized. We use the `loadedcross` element where we choose `w > L2/2`, so that the Cross is "unloaded", i.e. the center section is filled in with metalization: @@ -392,7 +388,7 @@ end ```` The above loop requires about 18 seconds of execution time on my machine. -Compare PSSFSS results to those digitized from the dissertation figure: +Compare PSSFSS results to those digitized from the dissertation: ````@example cross_on_dielectric_substrate col=[:red,:blue,:green] @@ -421,14 +417,13 @@ EditURL = "../literate/square_loop_absorber.jl" ``` ## Square Loop Absorber -This example is from Figure 7 of Costa and Monorchio: "A frequency selective -radome with wideband absorbing properties", *IEEE Trans. AP-S*, -Vol. 60, no. 6, June 2012, pp. 2740--2747. It shows how one can use the `polyring` +This example, from Figure 7 of [costa2012frequency](@cite), +shows how one can use the [`polyring`](@ref) function to model square loop elements. Three different designs are examined that employ different loop thicknesses and different values of sheet resistance. We compare the reflection coefficient magnitudes computed by PSSFSS with those digitized from the cited figure when the sheet is suspended -5 mm above a ground plane, hence we will also make use of the `pecsheet` function. +5 mm above a ground plane, hence we will also make use of the [`pecsheet`](@ref) function. ````@example square_loop_absorber using Plots, PSSFSS, DelimitedFiles @@ -530,9 +525,7 @@ EditURL = "../literate/flexible_absorber.jl" ``` ## Flexible Absorber -This example is from Yize Li, et al., "Ultra-wideband, polarization-insensitive flexible metamaterial -absorber base on laser printed graphene using equivalent circuit design method," Carbon, Vol 212, 2023, -available for free download from [here](https://doi.org/10.1016/j.carbon.2023.118166). +This example is from [li2023ultra](@cite). It uses square and circular resistive FSS elements sandwiched between layers of flexible dielectrics to realize a reflective absorber (i.e. a "rabsorber"). We compare the reflection coefficient magnitude computed by PSSFSS to that digitized @@ -571,13 +564,13 @@ scan = (θ=0, ϕ=0) results = analyze(strata, FGHz, scan, showprogress=false, resultfile=devnull, logfile=devnull) s11db = extract_result(results, @outputs s11db(te,te)) -yize = readdlm("../src/assets/yize2023_fig2a_s11db_digitized.csv", ',', Float64, '\n') +li = readdlm("../src/assets/li2023_fig2a_s11db_digitized.csv", ',', Float64, '\n') ps11 = plot(title="Normal Incidence Reflection Magnitude", xlabel="Frequency (GHz)", ylabel="20log₁₀|s₁₁|", xlim=(0,20), ylim=(-18,0), xtick=0:2:20, ytick=-20:2:0, framestyle=:box) plot!(ps11, FGHz, s11db, color=:blue, label="PSSFSS", lw=2) -plot!(ps11, yize[:,1], yize[:,2], color=:red, label="Yize et al.", lw=2) +plot!(ps11, li[:,1], li[:,2], color=:red, label="Li et al.", lw=2) savefig(ps11, "flexibleabsorbers11.png"); nothing # hide ```` @@ -593,11 +586,7 @@ EditURL = "../literate/splitringexample.jl" ``` ## Split-Ring Resonator -This example is taken from Figure 3 of -Fabian‐Gongora, et al, "Independently Tunable Closely Spaced Triband -Frequency Selective Surface Unit Cell Using the Third Resonant Mode of Split Ring Slots", -IEEE Access, 8/3/2021, Digital Object Identifier 10.1109/ACCESS.2021.3100325. - +This example is taken from [fabian2021independently; Figure 3](@cite). It consists of three concentric split rings with gaps sequentially rotated by 180° situated on a thin dielectric slab. @@ -657,14 +646,10 @@ EditURL = "../literate/reflectarray_example.jl" ``` ## Reflectarray Element -This example is taken from Figure 6 of -Li, Jiao, and Zhao: "A Novel Microstrip Rectangular-Patch/Ring- -Combination Reflectarray Element and Its Application", **IEEE Antennas and Wireless Propagation Letters**, -VOL. 8, 2009, pp. 1119-1112. - +This example is taken from [li2009novel; Figure 6](@cite). It generates the so-called "S-curve" for reflection phase of a reflectarray element. The element consists of a square patch in a square ring, separated from a ground plane by two dielectric layers. -Reflection phase is plotted versus the `L2` parameter as defined by Li, et al., which characterizes the +Reflection phase is plotted versus the `L2` parameter as defined in [li2009novel](@cite), which characterizes the overall size of the element. We start by defining a convenience function to generate a `RWGSheet` for a given value of `L2` in mm and @@ -856,20 +841,17 @@ EditURL = "../literate/band_pass_filter.jl" ``` ## Loaded Cross Band Pass Filter -This example is originally from Fig. 7.9 of B. Munk, *Frequency Selective Surfaces, -Theory and Design,* John Wiley and Sons, 2000. The same case was analyzed in L. Li, -D. H. Werner et al, "A Model-Based Parameter Estimation Technique for -Wide-Band Interpolation of Periodic Moment Method Impedance Matrices With Application to -Genetic Algorithm Optimization of Frequency Selective Surfaces", *IEEE Trans. AP-S*, -vol. 54, no. 3, March 2006, pp. 908-924, Fig. 6. Unfortunately, in neither reference +This example is originally from [munk2000fss; Fig. 7.9](@cite). +The same case was analyzed in [li2006model](@cite). Unfortunately, in neither reference are the dimensions of the loaded cross stated, except for the square unit cell period of 8.61 mm. I estimated the dimensions from the sketch in Fig. 6 of the second reference. To provide a reliable comparison, I analyzed one-eighth of the structure -in HFSS, a commercial finite element solver, using all three planes of symmetry -(using symmetry in the z = constant centerline plane +in HFSS, a commercial finite element solver, using all three planes of symmetry. +Exploiting the symmetry in the z = constant centerline plane required two analyses, once for an H-wall boundary condition, and once for an E-wall. Those -results were then combined using the method of Reed and Wheeler (even/odd symmetry)). With -a much reduced computational domain, it was then possible to drive HFSS well into convergence. +results were then combined using the method of Reed and Wheeler [reed1956method](@cite), i.e. +using even/odd symmetry. With a much reduced computational domain, it was then possible to +drive HFSS well into convergence. Two identical loaded cross slot-type elements are separated by a 6 mm layer of dielectric constant 1.9. Outboard of each sheet is a 1.1 cm layer of dielectric constant 1.3. @@ -950,10 +932,7 @@ to the two senses of circular polarization. We'll first look at analyzing a design presented in the literature, and then proceed to optimize another design using PSSFSS as the analysis engine inside the optimization objective function. ### Sjöberg and Ericsson Design -This example comes from the paper D. Sjöberg and A. Ericsson, "A multi layer meander line circular -polarization selective structure (MLML-CPSS)," The 8th European Conference on Antennas and Propagation -(EuCAP 2014), 2014, pp. 464-468, doi: 10.1109/EuCAP.2014.6901792. - +This example comes from [sjoberg2014multi](@cite). The authors describe an ingenious structure consisting of 5 progressively rotated meanderline sheets, which acts as a circular polarization selective surface: it passes LHCP (almost) without attenuation or reflection, and reflects RHCP (without changing its sense!) almost without attenuation or transmission. @@ -1171,7 +1150,7 @@ As previosly discussed, this is necessary because the rotated meanderlines are a the entire unit cell, and the unit cell for sheets 2 and 4 are not square. Since the periodicity of the sheets in the structure varies from sheet to sheet, higher order Floquet modes common to neighboring sheets cannot be defined, so we are forced to use only the dominant (0,0) modes which are independent of -the periodicity. This limitation is removed in a later example. +the periodicity. This limitation is removed in the next example. Meanwhile, it is of interest to note that their high-accuracy runs required 10 hours for CST and 19 hours for COMSOL on large engineering workstations versus about 22 seconds for PSSFSS on my desktop machine. @@ -1341,17 +1320,14 @@ The final sheet geometries and performance of this design are shown below:  As hoped for, the performance meets the more stringent design goals over a broader bandwidth than the -Sjöberg and Ericsson design, presumably because of the greater design flexibility allowed here. +Sjöberg and Ericsson design of [sjoberg2014multi](@cite), presumably because of the greater design flexibility allowed here. ```@meta EditURL = "../literate/cpss2.jl" ``` ## Meanderline/Strip-Based CPSS -This example comes from the same authors as the previous example. The paper is -A. Ericsson and D. Sjöberg, "Design and Analysis of a Multilayer Meander Line -Circular Polarization Selective Structure", IEEE Trans. Antennas Propagat., -Vol. 65, No. 8, Aug 2017, pp. 4089-4101. +This example comes from [ericsson2017design](@cite), by the same authors as the previous example. The design is similar to that of the previous example except that here, the two ``\pm 45^\circ`` rotated meanderlines are replaced with rectangular strips. This allows us to employ the `diagstrip` element and the `orient` keyword for the @@ -1485,7 +1461,7 @@ Dielectric layer information... ... ``` -Layers 3 and 9 were assigned 10 modes each. Layers 5 and 7, being thinner were assigned +Layers 3 and 9 were assigned 10 modes each. Layers 5 and 7, being thinner, were assigned 18 modes each. The numbers of modes are determined automatically by PSSFSS to ensure accurate cascading. @@ -1550,10 +1526,7 @@ EditURL = "../literate/splitring_cpss.jl" ## Split Ring-Based CPSS This circular polarization selective surface (CPSS) example comes from the paper -L.-X. Wu, K. Chen, T. Jiang, J. Zhao and Y. Feng, "Circular-Polarization-Selective -Metasurface and Its Applications to Transmit-Reflect-Array Antenna and Bidirectional -Antenna," in IEEE Trans. Antennas and Propag., vol. 70, no. 11, pp. 10207-10217, -Nov. 2022, doi: 10.1109/TAP.2022.3191213. +[wu2022circular](@cite) by Wu et al. The design consists of three sequentially rotated split rings separated by dielectric layers. Since the unit cells for all three rings are identical, PSSFSS can rigorously account for multiple scattering between the individual sheets using multiple @@ -1703,11 +1676,7 @@ EditURL = "../literate/angular_ss_example.jl" ``` ## Angle Selective Surface -This example is taken from Zhenting Chen, Chao Du, Jie Liu, Di Zhou, and Zhongxiang Shen, -"Design Methodology of Dual-Polarized Angle-Selective Surface Based -on Three-Layer Frequency-Selective Surfaces", IEEE Trans. Antennas Propagat., Vol. 71, No. 11, November 2023, -pp. 8704--8713. - +This example is taken [chen2023design](@cite). A three-sheet FSS is designed that is transparent to normally incident plane waves, but strongly attenuates obliquely incident waves. All three sheets are swastika-shaped, with the outer two sheets being identical. The shapes are generated in PSSFSS using the manji element. @@ -1809,14 +1778,293 @@ about 13.5 minutes to complete, much longer than for normal incidence. This is d the incremental phase shift ψ₁ is not held constant during the analysis, requiring the spatial integrals to be recomputed for each new incidence angle. +```@meta +EditURL = "../literate/checkerboard_example.jl" +``` + +## Checkerboard Metasurface +Checkerboard-style metasurfaces have been exensively studied, for example in [nakata2013plane](@cite) +and [kempa2010percolation](@cite). They exhibit some very unusual properties, which we will demonstrate here. + +We consider a series of checkerboad-like geometries. The square unit cell has dimension ``P = 5~\text{mm}``. +The side length for a PEC square sheet rotated 45° and perfectly inscribed in the unit cell is ``L_0 = P / \sqrt{2}``. +We will analyze a series of squares at 1 GHz with side lengths ``L = L_0 + δ``. The function `computed_and_plot` +defined below will take a given value of ``δ`` as its input, then perform four distinct analyses based on two +complementary triangulations. Let's run it first for ``δ = -0.5 \text{mm}`` and look +at the resulting triangulations: + +````@example checkerboard_example +using PSSFSS, Plots, PrettyTables +using Printf: @sprintf + +units = mm +P = 5 # period in x and y direction +L0 = P / √2 # Side length for self-complementary square +Nx = Ny = 10 + +function compute_and_plot(δ; P=P, units=units, L0=L0, Nx=Nx, Ny=Ny) + + Px = Py = P + Lx = Ly = L0 - abs(δ) + orient = 45 + islekwds = (; Lx, Ly, Px, Py, units, Nx, Ny, orient) + holekwds = (; units, s1=[P,0], s2=[0,P], sides=4, + a = iszero(δ) ? [P/2] : [(L0-abs(δ)) / √2], + b=[-2*Nx], ntri=2*Nx*Ny, structuredtri=false) + if δ ≤ 0 + redj = rectstrip(; class='J', islekwds...) + redm = rectstrip(; class='M', islekwds...) + bluej = polyring(; class='J', holekwds...) + bluem = polyring(; class='M', holekwds...) + else + redj = polyring(; class='J', holekwds...) + redm = polyring(; class='M', holekwds...) + bluej = rectstrip(; class='J', islekwds...) + bluem = rectstrip(; class='M', islekwds...) + end + + p1 = plot(redj, unitcell=true, faces=true, fillcolor=:red, title="Unit Cell") + p2 = plot(redj, faces=true, edges=false, fillcolor=:red, rep=(3,3), title="3×3 Array") + p3 = plot(bluej, unitcell=true, faces=true, fillcolor=:blue, title="Unit Cell") + p4 = plot(bluej, faces=true, edges=false, fillcolor=:blue, rep=(3,3), title="3×3 Array") + pl = plot(p1,p2,p3,p4, layout=(2,2), size=(550,600), plot_title=" Triangulations for δ = $δ") + + flist = 1 # Analysis frequency in GHz + steering = (θ=0, ϕ=0) # Normal incidence + logfile = resultfile = devnull # Suppress creation of output files + showprogress = false # Suppress screen output + + redjres = analyze([Layer(), redj, Layer()], flist, steering; logfile, resultfile, showprogress) + bluejres = analyze([Layer(), bluej, Layer()], flist, steering; logfile, resultfile, showprogress) + redmres = analyze([Layer(), redm, Layer()], flist, steering; logfile, resultfile, showprogress) + bluemres = analyze([Layer(), bluem, Layer()], flist, steering; logfile, resultfile, showprogress) + + + s11rj = extract_result(redjres, @outputs s11(v,v)) |> only + s21rj = extract_result(redjres, @outputs s21(v,v)) |> only + s11rm = extract_result(redmres, @outputs s11(v,v)) |> only + s21rm = extract_result(redmres, @outputs s21(v,v)) |> only + s11bj = extract_result(bluejres, @outputs s11(v,v)) |> only + s21bj = extract_result(bluejres, @outputs s21(v,v)) |> only + s11bm = extract_result(bluemres, @outputs s11(v,v)) |> only + s21bm = extract_result(bluemres, @outputs s21(v,v)) |> only + + return pl, (; s11rj, s21rj, s11rm, s21rm, s11bj, s21bj, s11bm, s21bm) +end + +pl, _ = compute_and_plot(-0.5) +pl +```` + +The function creates a "red" triangulation occupying the triangle of side length ``L_0 + δ``, and a "blue" triangulation, +occupying the complementary portion of the unit cell. Since for this case the red square +side length is 0.5 mm shorter than the critical length ``L_0``, it lies strictly inside the unit cell. So if we choose to +use the red triangulation to model electric surface current, then we can consider the red regions to be "islands" of +metal in otherwise empty space. We could also use the blue triangulation to model magnetic surface current, which again would +lead to the conclusion that the small untriangulated squares are conducting patches or "islands" of metalization. +Either of these two choices, when analyzed with PSSFSS, should yield the same values for computed reflection or transmission +coefficients (within modeling accuracy). + +A different approach would be to choose the red triangulation for representing magnetic surface current, in which case +the small red squares would represent "holes" in an otherwise solid metallic sheet. The same "hole" interpretation would +follow from using the blue triangulation to represent electric surface current. In fact, for this case, the blue region +in the full plane can be regarded as the union of an infinite number of metallic squares of dimension ``L_0 + δ``. +So positive values of ``δ`` can be handled by using ``-δ`` and reversing the roles of the red and blue triangulations. +This is exactly what is done in the function `compute_and_plot` above. + +We'll now exercise the function for the set of ``δ`` values ``\{ -0.2, -0.05, 0, 0.05, 0.2 \}``, observing both the plotted +triangulations and the resulting scattering parameter predictions for each of the four modeling choices outlined above. + +```julia +δs = [-0.2, -0.05, 0, 0.05, 0.2] # Departure in mm from self-complementary square side length +s11rj, s21rj, s11rm, s21rm, s11bj, s21bj, s11bm, s21bm = (zeros(ComplexF64, length(δs)) for _ in 1:8) +for (i, δ) in pairs(δs) + plt, r = compute_and_plot(δ) + s11rj[i] = r.s11rj + s21rj[i] = r.s21rj + s11rm[i] = r.s11rm + s21rm[i] = r.s21rm + s11bj[i] = r.s11bj + s21bj[i] = r.s21bj + s11bm[i] = r.s11bm + s21bm[i] = r.s21bm + display(plt) +end +``` + +.png) + +.png) + +.png) + +.png) + +.png) + +Let the letter "J" denote use of a triangulation to represent electric surface current and let "M" denote magnetic surface current. +So, e.g. "Blue M" means that the blue triangulation is used to represent magnetic current. We can make the following observations about +the above plots: + +1. The red and blue triangulations alway occupy complementary regions of the unit cell. Then for a given choice + of ``δ``, PSSFSS analysis of "Blue J" and "Red M" should result identical scattering parameters (apart from modeling errors). + Likewise analysis of "Red J" and "Blue M" should result in the same scattering parameters. +2. It is well known from Babinet's principle [nakata2013plane](@cite) [tan2012babinet](@cite) that reflection coefficients + of thin perforated screens are equal to the negative of the transmission coefficients of the complementary structure, provided the + structures exhibit sufficient rotational symmetry as is the case here. "Red J" and "Blue J" form such a + complementary pair, as do "Red M" and "Blue M". +3. All of the ``δ`` values are quite small compared to ``L_0 \approx 3.5355 \text{mm}``. In particular, the plots for + ``δ = \pm 0.05`` are almost indistinguishable from the plot for ``δ = 0``. Therefore, we expect the scattering + parameters of all of these cases to be nearly equal. +4. For the case of ``δ = 0``, the red and blue regions are both self-complementary. From our earlier considerations, we then + expect equal reflection and transmission coefficients for any of the choices "Red J", "Blue J", "Red M", or "Blue M". In + fact all four of these cases should yield the same values. + +The following code generates a summary table showing how well these expectations are satisfied: + +```julia +mat = hcat(s11rj, s11bm, -s21rm, -s21bj) |> transpose +row_labels = ["S₁₁: Red J ", "S₁₁: Blue M", "-S₂₁: Red M ", "-S₂₁: Blue J"] +header = (["δ = $δ mm" for δ in δs], ["islands or holes" for δ in δs]) +header[2][findfirst(iszero, δs)] = "SC (Self-Complementary)" +formatters = (v,i,j) -> imag(v) > 0 ? @sprintf("%7.4f + %6.4fim", real(v), imag(v)) : + @sprintf("%7.4f - %6.4fim", real(v), -imag(v)) +highlighters = Highlighter((data, i, j) -> (j == findfirst(iszero,δs)), crayon"red bold") +pretty_table(mat; header, row_labels, alignment=:c, formatters, highlighters) +``` + +````@example checkerboard_example +for line in eachline("./assets/checkerboard_prettytable.data"); println(line); end #hide +```` + +For ``δ < 0`` the results seem reasonable. In this regime of electrically small metal islands (for Red J and Blue M), +one wouldn't expect much reflection, and that is what is observed. For ``δ > 0``, the geometry (for Red J and Blue M) +consists of a metal plate with small holes, so one would expect almost total reflection, as observed. Now consider +how well our above observations are aligned with the computed results... + +From observations 1 and 2, the numbers in any one column should all be approximately equal, which they are +except for the center ``δ = 0`` column. From observation 3 we expect all the entries +in any one row to be nearly equal, but this is not true at all. In any one row there is a violent jump in +amplitude at or near ``δ = 0``. Finally, from observation 4 we expect all of the numeric entries for ``δ = 0`` +to be approximately equal--which they are most definitely not. What is going on here? + +As discussed in the previously cited references, the idealized model being analyzed for ``δ = 0`` is unphysical. +As ``\delta \rightarrow 0`` the response of the structure does not approach a limit, since the violent jump in +reflection coefficient occurs for arbitrarily small positive and negative excursions from ``L_0``. Any physically +realizable structure must exhibit continuous dependence on its physical parameters. If one attempts +to approximate this ideal surface in the real world one must use finite thickness and conductivity, so that neighboring +unit cells will intersect all along the thickness of the metal, rather than at a single point, thus destroying the +self-complementary property. Similarly, finite losses in any real metal preclude self-complementarity. + +The same problem observed here for an ideal model of a self-complementary surface arises when analyzing pixelated +structures such as shown in the next example, where adjacent metal pixels +can intersect at only a single point. In that case, using a "Blue M" approach is known to agree with measurements. + +```@meta +EditURL = "../literate/sympixels_example.jl" +``` + +## Pixelated Band Pass Filter +This example of from [hong2021design; Fig. 4a](@cite). It consists of a single-layer +band-pass filter consisting of many square metallic "pixels", optimized using a special binary +optimizer described in the cited paper. + +The geometry for the pixelated filter is created using the [`sympixels`](@ref) function as +show below. + +````@example sympixels_example +using PSSFSS, Plots +P = 5.4 +units = mm +nrim = 1 +halfnint = 26 +patternvec = Bool[ + 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 1, 1, 0, 1, 1, 0, 0, 1, 0, 0, 1, 0, 0, 0, 1, + 0, 1, 0, 1, 0, 1, 1, 1, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 1, 1, 0, 1, 0, + 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 1, 1, 0, 0, 1, 0, 1, 1, 1, 1, 0, 0, 1, + 0, 1, 0, 0, 0, 1, 0, 1, 1, 1, 1, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, + 1, 1, 1, 0, 0, 1, 0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 0, 1, 1, 0, 0, 1, + 1, 0, 1, 1, 1, 0, 1, 0, 1, 0, 1, 1, 1, 0, 0, 1, 1, 0, 1, 0, 1, + 0, 0, 0, 1, 0, 0, 1, 0, 0, 1, 1, 1, 0, 1, 0, 0, 1, 0, 0, 0, + 0, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 1, 0, 1, 1, 0, 0, 1, 0, + 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 1, 1, + 1, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 1, 0, 0, 1, 1, + 0, 1, 0, 1, 0, 0, 1, 0, 0, 1, 0, 1, 0, 0, 1, 1, + 1, 1, 0, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 0, 0, + 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, + 1, 1, 1, 0, 1, 0, 0, 1, 1, 1, 1, 0, 1, + 1, 1, 1, 1, 1, 0, 0, 1, 1, 0, 1, 1, + 0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 1, + 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, + 0, 1, 0, 0, 0, 1, 0, 0, 0, + 0, 0, 1, 1, 0, 0, 1, 1, + 0, 0, 0, 0, 1, 1, 0, + 1, 0, 0, 1, 0, 0, + 0, 0, 1, 1, 1, + 1, 0, 0, 0, + 0, 0, 1, + 0, 0, + 0] + + +sheet1 = sympixels(; P, units, nrim, halfnint, patternvec, pdiv=1, class='M') +plot(sheet1, size=(450,450), unitcell=true) +```` + +As described in the documentation for [`sympixels`](@ref), the `1` entries in `patternvec` above +denote the locations of the metallized pixels. However, by setting `class = 'M'` above we choose to +triangulate the empty pixel locations, i.e. those that are not occupied by metal. As discussed in +[Checkerboard Metasurface](@ref), this choice is necessary to avoid spurious results when analyzing +a structure of this type. + +Note that the above triangulation is created by forming a square for each triangulated pixel and then +adding a single diagonal. A finer triangulation is required for an accurate analysis result. This can +be accomplished by increasing `pdiv` to a value greater than `1`. Here is the triangulation that results +from `pdiv = 2`: + +````@example sympixels_example +sheet2 = sympixels(; P, units, nrim, halfnint, patternvec, pdiv=2, class='M') +plot(sheet2, lw=0.5, size=(450,450), unitcell=true) +```` + +Now each pixel has been divided into an array of 2×2 squares, each of which receives a diagonal to +form triangles. The code for analyzing the `pdiv = 1` case would look like this: +```Julia +strata = [Layer(), sheet1, Layer(epsr=3.28, tandel=0.007, width=0.05mm), Layer()] +steering = (θ=0, ϕ=0) +flist = range(start=20, stop=36, length=401) +logfile = "bp_filter_pdiv1.log" +resultfile = "bp_filter_pdiv1.res" +results = analyze(strata, flist, steering; logfile, resultfile) +s21db = extract_result(results, @outputs s21db(te,te)) +``` + +Accurately analyzing a structure like this consisting of a huge number of pixels can be become +very expensive (in terms of memory usage and CPU time) very quickly. For example, the `pdiv = 1` +case generated 3056 basis functions and took about 30 seconds to analyze on a 32 GByte machine. +But the `pdiv = 2` case generated 15,088 unknowns and required 683 seconds of execution time. It +also failed on the 32 GByte machine due to lack of memory, but was run successfully on a 64 GByte +machine. + +Because the analysis in [hong2021design](@cite) used nonzero metal thickness, I ran the zero-thickness +geometry in HFSS to observe the significance of the thickness. Here is a plot comparing the digitized +results from [hong2021design; Fig. 4a](@cite), with the zero-thickness HFSS analysis and the two +PSSFSS analyses: + + + +As seen above, the finite metallization thickness does not have a large effect on the transmission trace. +Using `pdiv = 2` produces much better agreement with HFSS, but it appears that an even larger value would +be required to achieve full convergence. + ```@meta EditURL = "../literate/tepfile_creation_example.jl" ``` ## TEP File Creation -Here we show how to create a TICRA-compatible TEP (tabulated electrical properties) file using PSSFSS. -The geometry for this example is a rectangular copper strip measuring 4 cm × 0.2 cm in a 5 cm square unit cell. -The code for analyzing this geometry and creating the TEP file is shown below: +Here we show how to create a TICRA-compatible TEP (tabulated electrical properties) file using the [`res2tep`](@ref) +function included with PSSFSS. The geometry for this example is a rectangular copper strip measuring +4 cm × 0.2 cm in a 5 cm square unit cell. The code for analyzing this geometry and creating the TEP file is shown below: ```Julia using PSSFSS FGHz = 3.0 @@ -1828,11 +2076,14 @@ Ny = 6 sheet = rectstrip(; Px, Py, Lx, Ly, Nx, Ny, units=cm, sigma=5.7e8) steering = (θ=0:5:70, ϕ=0:15:345) strata = [Layer(), sheet, Layer()] -results = analyze(strata, FGHz, steering) +resultfile = "strip.res" +results = analyze(strata, FGHz, steering; resultfile) res2tep(results, "dipole_pssfss.tep"; name = "dipole", class = "pssfss") +# Alternatively: res2tep(resultfile, "dipole_pssfss.tep"; name = "dipole", class = "pssfss") ``` Note that a very fine discretization has been specified, resulting in 1560 triangles. There are also a large number of steering angles requested (15 × 24 = 360). This analysis required about 155 seconds on my machine. +Please see the documentation for [`res2tep`](@ref) for requirements on setting up analysis scan angles. For comparison, the same geometry was analyzed using the QUPES program of Ticra Tools Student Edition 2024 (hence the specification for `sigma` above, which is the default conductivity used by QUPES). For the @@ -1842,3 +2093,34 @@ for any of these 360 steering angles was approximately 0.0021. Convergence stud complex scattering coefficients were still varying slightly in the third decimal place for the settings used in this example. +```@meta +EditURL = "../literate/fresneltable_creation_example.jl" +``` + +## Fresnel Table Creation +Here we show how to create an HFSS-compatible Fresnel table using the [`res2fresnel`](@ref) +function included with PSSFSS. The geometry for this example is a double square loop +in a 5 cm square unit cell. The code for analyzing this geometry and creating the +Fresnel table is shown below: +```Julia +using PSSFSS +dwidth = 3mm +duroid = Layer(epsr=2.2, tandel=0.0009, width=dwidth) +a = √2 * [1, 2.125] +b = √2 * [1.5, 3.125] +units = mm +sides = 4 +orient = 45 +sheet = polyring(;a, b, units, sides, orient, s1=[8,0], s2=[0,8], ntri=2000) +strata = [Layer(), duroid, sheet, duroid, Layer(width=-2*dwidth)] +steering = (θ=0:5:45, ϕ=0) +FGHz = 10:2:20 +results = analyze(strata, FGHz, steering; resultfile="double_square_loop.res") +res2fresnel(results, "double_square_loop.rttbl") +# Alternatively: res2fresnel("double_square_loop.res", "double_square_loop.rttbl") +``` +Note that the final layer specifies a width that is the negative of the sum of the remaining layers' widths. +This has the effect of moving the output phase reference plane to coincide with the input phase reference plane. +Please see the documentation for [`res2fresnel`](@ref) for discussion of the setup requirements and restrictions +necessary for creating a valid Fresnel table. + diff --git a/docs/src/function_index.md b/docs/src/function_index.md index 6e2d728d..8300b78b 100644 --- a/docs/src/function_index.md +++ b/docs/src/function_index.md @@ -1,5 +1,5 @@ -# Index +# Public API Index ```@index ``` diff --git a/docs/src/manual.md b/docs/src/manual.md index 074b7e2c..56a34268 100644 --- a/docs/src/manual.md +++ b/docs/src/manual.md @@ -25,9 +25,7 @@ of arbitrary polarization incident from either Region ``1`` or Region ``N``. ### A Quick Example Here is an example run of PSSFSS for a 4-sheet meanderline polarizer design -from one of the first papers on the subject: T. L. Blackney, J. R. Burnett, and -S. B. Cohn, “A design method for meander-line circular polarizers” -presented at 22nd Annual Antenna Symposium, Oct. 1972. +from [blbc:72](@cite), one of the first papers on the subject. Detailed explanations of the code are omitted for now. ```julia @@ -72,7 +70,7 @@ presented in the paper. Here is the comparison plot:  -The PSSFSS run took about 24 seconds on my machine for 4 meanderline PSS sheets +The PSSFSS run took about 14 seconds on my machine for 4 meanderline PSS sheets analyzed at 61 frequencies. I hope this example whetted your appetite to learn @@ -121,7 +119,7 @@ An alternative to using a text editor is to run scripts from a Jupyter or Pluto [IJulia](https://github.com/JuliaLang/IJulia.jl) package for further details on using Jupyter, and the [Pluto Home Page](https://plutojl.org/) for details of using Pluto. -It is almost essential that PSSFSS users also install the [Plots](https://github.com/JuliaPlots/Plots.jl) +It is essential that PSSFSS users also install the [Plots](https://github.com/JuliaPlots/Plots.jl) package. This will allow easy visualization of the FSS/PSS element triangulations produced by PSSFSS, in addition to providing a convenient means to plot analysis results. @@ -148,8 +146,12 @@ Here are the steps in the analysis process. Each of these steps is examined in detail below... ## Strata -### Layer -Dielectric layers are created with the [`Layer`](@ref) function: +The geometry to be analyzed is passed as the first argument of the [`analyze`](@ref) function in the form +of a Julia `Vector` (i.e., a one-dimensional array). This vector must contain at least two `Layer` +objects representing dielectric layers, and zero or more `RWGSheet` instances, representing the zero-thickness +PSS or FSS sheets located between the dielectric layers. These are described below. +### Dielectric Layers +Dielectric layers are created with the [`Layer`](@ref) constructor function: ```@repl manual using PSSFSS # Brings PSSFSS functions and types into scope @@ -178,7 +180,7 @@ created by calling a constructor function for a particular style of sheet: patch = rectstrip(Nx=10, Ny=10, Px=1, Py=1, Lx=0.5, Ly=0.5, units=cm) ``` -The call to [`rectstrip`](@ref) above creates a `RWGSheet` object for a rectangular strip +The call to [`rectstrip`](@ref) above creates an `RWGSheet` object for a rectangular strip of dimensions 0.5 cm in the x and y directions, lying in a square unit cell of dimension 1 cm. The triangulation uses 10 edges in the x and y directions (`Nx` and `Ny`). @@ -218,8 +220,7 @@ distribution function for roughness is specified by the keyword `disttype` which default value, suitable for the so-called "oxide side" of a planar conductor, or `:rayleigh` (suitable for the "foil side" of the conductor, i.e. the side bonded to the dielectric substrate). Together, the conductivity, surface roughness, and roughness distribution type are used by PSSFSS internally to compute a frequency-dependent surface -impedance, using the so-called Gradient Model, as described in D. N. Grujić, “Simple and Accurate Approximation of -Rough Conductor Surface Impedance,” IEEE Trans. Microwave Theory Tech., vol. 70, no. 4, pp. 2053-2059, April 2022. +impedance, using the so-called Gradient Model, as described in [grujic2021simple](@cite). Obviously, only one of `Zsheet` and `sigma` may be specified as keyword arguments for a given `RWGSheet`. #### Perfectly Conducting Walls @@ -419,7 +420,7 @@ e.g. `S21ang(te,tm)`. Again, case is not significant, and "1" and "te" can be f It is often desired to use a set of polarization basis vectors other than TE/TM to define field coefficients. PSSFSS supports in addition to TE/TM the use of H/V for horizontal/vertical components -(in the [Ludwig 3](https://ieeexplore.ieee.org/document/1140406) sense), and L/R for left-hand circular +(in the Ludwig 3 [ludw:73](@cite) sense), and L/R for left-hand circular and right-hand circular polarization (in the [IEEE sense](https://en.wikipedia.org/wiki/Circular_polarization#Uses_of_the_two_conventions)). @@ -514,7 +515,8 @@ parameter of the [`@outputs`](@ref) macro. #### Using the `analyze` Function Return Value The variable returned from [`analyze`](@ref) can be used in a similar way to the result file to obtain any quantity -available to the [`@outputs`](@ref) macro, but for this purpose we use the [`extract_result`](@ref) function: +available to the [`@outputs`](@ref) macro. For this purpose we pass the return value of the `analyze` function +as the first argument to the [`extract_result`](@ref) function: ```julia julia> results = analyze(strata, flist, steering, outlist=outputs); @@ -536,8 +538,27 @@ julia> dat = extract_result(results, @outputs FGHz ar21db(h)) As in the previous example, the value returned by [`extract_result`](@ref) is a two-dimensional array (a `Matrix`), with each column corresponding to a parameter of the [`@outputs`](@ref) macro. -### Saving Analysis Results As TEP Files +## Exporting Analysis Results +Besides the ability to export results to CSV files as described [here](@ref "The `outlist` Keyword Argument to `analyze`"), +one can export results to TICRA-compatible TEP files or HFSS-compatible Fresnel tables, as described in the following +subsections. +### Exporting Analysis Results to TEP Files The vector of `Result` objects returned by the [`analyze`](@ref) function (or the corresponding result file) can be converted to a TICRA-compatible TEP (tabulated electrical properties) file -using the function [`res2tep`](@ref). The code needed to do so is shown in [this example](@ref "TEP File Creation"). +using the function [`res2tep`](@ref). TEP files contain the equivalent of the full +4×4 scattering matrix computed by PSSFSS. Therefore, there is no requirement/limitation for unit cell geometry +to exhibit any particular symmetry when results are to be exported to a TEP file. There are, however, requirements +on the choice of steering angles to be analyzed. For details, see the documentation of [`res2tep`](@ref). +Sample code for creating a TEP file is shown in [this example](@ref "TEP File Creation"). + +### Exporting Analysis Results to Fresnel Tables +The vector of `Result` objects returned by the [`analyze`](@ref) function (or the corresponding +result file) can be converted to an HFSS SBR+-compatible Fresnel table +using the function [`res2fresnel`](@ref). Fresnel tables contain reflection +and transmission coefficients for only "front" surface incidence, and only for co-polarized +(TE → TE and TM → TM) scattering. In addition, coefficients for only one azimuth angle ϕ are +retained in such a table. These limitations impose restrictions on the type of structure that +can be analyzed for use in generating a valid Fresnel table. For details, see the documentation +of [`res2fresnel`](@ref). Sample code for creating a Fresnel table is shown in +[this example](@ref "Fresnel Table Creation"). diff --git a/docs/src/reference.md b/docs/src/reference.md index f6cc880c..6eb6440f 100644 --- a/docs/src/reference.md +++ b/docs/src/reference.md @@ -1,5 +1,5 @@ -# [PSSFSS](https://github.com/simonp0420/PSSFSS) Function Reference +# [PSSFSS](https://github.com/simonp0420/PSSFSS) Public API Documentation ```@docs analyze @@ -16,12 +16,16 @@ meander nodecount @outputs pecsheet +pixels pmcsheet polyring read_sheet_data rectstrip +res2fresnel res2tep +RWGSheet sinuous splitring +sympixels ``` diff --git a/docs/src/references.md b/docs/src/references.md new file mode 100644 index 00000000..78f29bc4 --- /dev/null +++ b/docs/src/references.md @@ -0,0 +1,4 @@ +# References + +```@bibliography +``` diff --git a/docs/src/refs.bib b/docs/src/refs.bib new file mode 100644 index 00000000..30314472 --- /dev/null +++ b/docs/src/refs.bib @@ -0,0 +1,234 @@ +@inproceedings{blbc:72, + title = "A design method for meander-line circular polarizers", + author = "T. L. Blackney and J. R. Burnett and S. B. Cohn", + booktitle = "22cd Annual Antenna Symposium", + year = "1972", + month = oct, + kwds = "meanderline polarizer", + address = "University of Illinois" +} + +@article{grujic2021simple, + title={Simple and accurate approximation of rough conductor surface impedance}, + author={Dušan N. Grujić}, + journal={IEEE Trans. Microwave Theory Tech.}, + volume={70}, + number={4}, + pages={2053-2059}, + year={2021}, + publisher={IEEE} +} + +@article{ludw:73, + title = "The definition of cross polarization", + author = "Arthur C. Ludwig", + journal={IEEE Trans. Antennas Prop.}, + year = "1973", + volume = "AP-21", + number = "1", + pages = "116-119", + month = "January", +} + +@article{chen2023design, + title={Design methodology of dual-polarized angle-selective surface based on three-layer frequency-selective surfaces}, + author={Chen, Zhenting and Du, Chao and Liu, Jie and Zhou, Di and Shen, Zhongxiang}, + journal={IEEE Trans. Antennas Prop.}, + volume={71}, + number={11}, + pages={8704--8713}, + year={2023}, + publisher={IEEE} +} + +@book{munk2000fss, + title = "Frequency Selective Surfaces: Theory and Design", + booktitle = "Frequency Selective Surfaces: Theory and Design", + author = "Ben A. Munk", + year = 2000, + kwds = "FSS, Array, Frequency Selective Surface", + publisher = "John Wiley \& Sons, Inc." +} + +@article{li2006model, + title={A model-based parameter estimation technique for wide-band interpolation of periodic moment method impedance matrices with application to genetic algorithm optimization of frequency selective surfaces}, + author={Li, Ling and Werner, Douglas H and Bossard, Jeremy A and Mayer, Theresa S}, + journal={IEEE Trans. Antennas Prop.}, + volume={54}, + number={3}, + pages={908--924}, + year={2006}, + publisher={IEEE} +} + +@article{reed1956method, + title={A method of analysis of symmetrical four-port networks}, + author={Reed, John and Wheeler, GJ}, + journal={IRE Trans. Microwave Theory Tech.}, + volume={4}, + number={4}, + pages={246--252}, + year={1956}, + publisher={IEEE} +} + +@inproceedings{sjoberg2014multi, + title={A multi layer meander line circular polarization selective structure (MLML-CPSS)}, + author={Sjöberg, Daniel and Ericsson, Andreas}, + booktitle={The 8th European Conference on Antennas and Propagation (EuCAP 2014)}, + pages={464--468}, + year={2014}, + organization={IEEE} +} + +@article{ericsson2017design, + title={Design and analysis of a multilayer meander line circular polarization selective structure}, + author={Ericsson, Andreas and Sjöberg, Daniel}, + journal={IEEE Trans. Antennas Prop.}, + volume={65}, + number={8}, + pages={4089--4101}, + year={2017}, + publisher={IEEE} +} + +@article{barlevy1997fundamental, + title={Fundamental constraints on the electrical characteristics of frequency selective surfaces}, + author={Barlevy, Alon S and Rahmet-Samii, Yahya}, + journal={Electromagnetics}, + volume={17}, + number={1}, + pages={41--68}, + year={1997}, + publisher={Taylor \& Francis} +} + +@phdthesis{barlevy1998dissertation, + title={Properties of electromagnetic scattering from periodic metallic structures: Application to frequency-selective and bandgap structures}, + author={Barlevy, Alon S}, + school={University of California, Los Angeles}, + type={PhD thesis}, + year={1998}, +} + +@article{li2023ultra, + title={Ultra-wideband, polarization-insensitive flexible metamaterial absorber base on laser printed graphene using equivalent circuit design method}, + author={Li, Yize and Fang, Yixian and Huang, Yihe and Pan, Kewen and Xiao, Xiaoyu and Liu, Xuzhao and Li, Lin and Hu, Zhirun}, + journal={Carbon}, + volume={212}, + pages={118166}, + year={2023}, + publisher={Elsevier}, + note = {Available at \url{https://doi.org/10.1016/j.carbon.2023.118166}}, +} + +@article{li2009novel, + title={A novel microstrip rectangular-patch/ring-combination reflectarray element and its application}, + author={Li, Qin-Yi and Jiao, Yong-Chang and Zhao, Gang}, + journal={IEEE Antennas Wireless Prop. Letters}, + volume={8}, + pages={1119--1122}, + year={2009}, + publisher={IEEE} +} + +@article{wu2022circular, + title={Circular-polarization-selective metasurface and its applications to transmit-reflect-array antenna and bidirectional antenna}, + author={Wu, Lin-Xiao and Chen, Ke and Jiang, Tian and Zhao, Junming and Feng, Yijun}, + journal={IEEE Trans. Antennas Prop.}, + volume={70}, + number={11}, + pages={10207--10217}, + year={2022}, + publisher={IEEE} +} + +@article{fabian2021independently, + title={Independently tunable closely spaced triband frequency selective surface unit cell using the third resonant mode of split ring slots}, + author={Fabian-Gongora, Henry and Martynyuk, Alexander E and Rodriguez-Cuevas, Jorge and Martinez-Lopez, Lourdes and Martinez-Lopez, Rosalba and Martinez-Lopez, Jose I}, + journal={IEEE Access}, + volume={9}, + pages={105564--105576}, + year={2021}, + publisher={IEEE} +} + +@article{costa2012frequency, + title={A frequency selective radome with wideband absorbing properties}, + author={Costa, Filippo and Monorchio, Agostino}, + journal={IEEE Trans. Antennas Prop.}, + volume={60}, + number={6}, + pages={2740--2747}, + year={2012}, + publisher={IEEE} +} + +@book{collin1991field, + title={Field Theory of Guided Waves}, + author={Collin, Robert E.}, + edition = {2cd}, + year={1991}, + publisher={IEEE Press} +} + +@book{weinstein1969theory, + author = {L. A. Weinstein}, + title = {The Theory of Diffraction and the Factorization Method}, + publisher = {The Golem Press}, + address = {Boulder, Colorado}, + year = {1969}, +} + +@article{daniele1990diffraction, + title={Diffraction of a plane wave by a strip grating}, + author={Daniele, Vito and Gilli, Marco and Viterbo, Emanuele}, + journal={Electromagnetics}, + volume={10}, + number={3}, + pages={245--269}, + year={1990}, + publisher={Taylor \& Francis} +} + +@article{tan2012babinet, + title={Babinet’s principle for electromagnetic fields}, + author={Tan, Zhong Ming and McDonald, Kirk T}, + journal={Joseph Henry Laboratories, Princeton University, Princeton, NJ}, + volume={8544}, + year={2012}, + note = {Available at \url{http://kirkmcd.princeton.edu/examples/babinet.pdf}}, +} + +@article{nakata2013plane, + title={Plane-wave scattering by self-complementary metasurfaces in terms of electromagnetic duality and Babinet's principle}, + author={Nakata, Yosuke and Urade, Yoshiro and Nakanishi, Toshihiro and Kitano, Masao}, + journal={Physical Review B—Condensed Matter and Materials Physics}, + volume={88}, + number={20}, + pages={205138}, + year={2013}, + publisher={APS}, + note = {Available at \url{https://arxiv.org/pdf/1311.3839}}, +} + +@article{kempa2010percolation, + title={Percolation effects in the checkerboard Babinet series of metamaterial structures}, + author={Kempa, K}, + journal={physica status solidi (RRL)--Rapid Research Letters}, + volume={4}, + number={8-9}, + pages={218--220}, + year={2010}, + publisher={Wiley Online Library} +} + +@article{hong2021design, + title={Design of single-layer metasurface filter by conformational space annealing algorithm for 5G mm-wave communications}, + author={Hong, Young-Pyo and Hwang, In-June and Yun, Dal-Jae and Lee, Dong-Joon and Lee, In-Ho}, + journal={IEEE Access}, + volume={9}, + pages={29764--29774}, + year={2021}, + publisher={IEEE} +} \ No newline at end of file diff --git a/docs/test/Project.toml b/docs/test/Project.toml new file mode 100644 index 00000000..81648c0b --- /dev/null +++ b/docs/test/Project.toml @@ -0,0 +1 @@ +[deps] diff --git a/src/Elements.jl b/src/Elements.jl index 72be295c..1924cdd7 100644 --- a/src/Elements.jl +++ b/src/Elements.jl @@ -1,11 +1,12 @@ module Elements -export diagstrip, jerusalemcross, loadedcross, manji, meander, pecsheet, pmcsheet, polyring, rectstrip, sinuous, splitring +export diagstrip, jerusalemcross, loadedcross, manji, meander, pecsheet, + pixels, pmcsheet, polyring, rectstrip, sinuous, splitring, sympixels using ..PSSFSSLen: mm, cm, inch, mil, PSSFSSLength -using ..Sheets: RWGSheet, rotate!, translate!, combine, recttri, SV2, orient! +using ..Sheets: RWGSheet, rotate!, translate!, combine, recttri, SV2, orient!, test_fefv using ..Meshsub: meshsub -using StaticArrays: SA +using StaticArrays: SA, @SVector using LinearAlgebra: norm, ⋅, × using Printf: @sprintf import LibGEOS # difference, readgeom, Polygon, MultiPolygon @@ -201,8 +202,8 @@ const optional_kwargs = """ with those of pair (C,D) by a simple translation. If there are many such equivalent pairs, a significant decrease in matrix fill time ensues by exploiting the equivalence. The tradeoff is the time needed to identify them. The default value is `true` for the `strip`, `diagstrip`, - `meander`, `manji`, `loadedcross`, `jerusalemcross`, and 4-sided `polyring` styles (those employing - structured meshes) and `false` for the remaining styles (those employing unstructured meshes). + `meander`, `manji`, `loadedcross`, `jerusalemcross`, `pixels`, 4-sided `polyring` styles (those + employing structured meshes), and `sympixels`, and `false` for the remaining styles (those employing unstructured meshes). - `save::String=""` Specifies a file name to which the sheet triangulation and unit cell data is to be written, typically to be plotted later. @@ -613,7 +614,7 @@ function jerusalemcross_unstructured(; P::Real, L1::Real, L2::Real, A::Real, B:: end sheet.style = "jerusalemcross" - sheet.ξη_check = false + sheet.ξη_check = true sheet.units = units sheet.s₁ = s1 sheet.s₂ = s2 @@ -862,7 +863,7 @@ function loadedcross_unstructured(; s1::Vector{<:Real}, s2::Vector{<:Real}, L1:: end sheet.style = "loadedcross" - sheet.ξη_check = false + sheet.ξη_check = true sheet.units = units sheet.s₁ = s1 sheet.s₂ = s2 @@ -1965,4 +1966,6 @@ end include("structuredtri.jl") # Code for structured meshes for loadedcross, jerusalemcross +include("pixels.jl") # Code for pixels and sympixels elements + end # module diff --git a/src/FastSweep.jl b/src/FastSweep.jl index d7699713..72bea408 100644 --- a/src/FastSweep.jl +++ b/src/FastSweep.jl @@ -25,7 +25,7 @@ Rational function interpolation using a Path II Neville lozenge, as defined in t ## Return Values - `Sinterp`: The rational function interpolation at `x0`. - - `errest`: An error estimate for `norm(Sinterp - S(x0), Inf)`. + - `errest`: An error estimate for `norm(Sinterp - S(x0))`. ## Reference Ma, X., Wan, G. and Wan, W., 2012. "A Multi-Dimensional Adaptive Sampling Method for Analysis @@ -86,7 +86,7 @@ function interp_path2( Sjk[j] += (x0j[j+k] - x0j[j]) * den1 .* den2 ./ bigden end end - return (Sjk[0], norm(Sjk[0] - Sjkm1[0]), Inf) + return (Sjk[0], norm(Sjk[0] - Sjkm1[0])) end fixbigden(x::T) where {T <: Number} = ifelse(iszero(x), one(T), x) diff --git a/src/Meshsub.jl b/src/Meshsub.jl index 3437030d..fa64d132 100644 --- a/src/Meshsub.jl +++ b/src/Meshsub.jl @@ -61,15 +61,16 @@ function meshsub(; points::Matrix{<:Real}, seglist::Matrix{<:Integer}, iter = 1 end + # Set up for call to triangulate: # Set up for call to triangulate: triin = TriangulateIO() triin.pointlist = Matrix{Cdouble}(points) triin.segmentlist = Matrix{Cint}(seglist) isempty(segmarkers) || (triin.segmentmarkerlist = Vector{Cint}(segmarkers)) isempty(holes) || (triin.holelist = Matrix{Cdouble}(holes)) + triout = TriangulateIO() # Establish scope areaold = area - triout = deepcopy(triin) # Establish scope for k in 1:iter # Perform triangulation: (triout, vorout) = triangulate(switches, triin) diff --git a/src/Outputs.jl b/src/Outputs.jl index 80a01455..f7f7239b 100644 --- a/src/Outputs.jl +++ b/src/Outputs.jl @@ -1,5 +1,5 @@ module Outputs -export @outputs, Result, append_result_data, read_result_file, extract_result_file, extract_result, res2tep +export @outputs, Result, append_result_data, read_result_file, extract_result_file, extract_result, res2fresnel, res2tep using LinearAlgebra: ⋅, norm, × using ..UnitVectors: ẑ @@ -12,6 +12,8 @@ using StaticArrays: @SVector, @SMatrix, SMatrix using JLD2: JLD2, jldopen using FileIO: load using TicraUtilities: TicraUtilities +using Printf: @printf +using Dates: now @enum HorV H = 1 V = 2 @enum RorL R = 1 L = 2 @@ -492,7 +494,7 @@ function _check_results_for_tep!(results::Vector{Result}) phis == phi_range || error("ϕ values must be equivalent to 0:Δϕ:(360-Δϕ)") dtheta = thetas[2] - thetas[1] theta_range = 0:dtheta:last(thetas) - thetas == theta_range || error("θ values must be equivalent to 0:Δθ:(360-Δθ)") + thetas == theta_range || error("θ values must be equivalent to 0:Δθ:θmax") freqs_vec = [f*u"GHz" for f in freqs] # Verify that full cartesian product of angles and frequencies is present: i = 0 @@ -519,6 +521,15 @@ TEP (tabulated electrical properties) file. If the first positional argument is assumed to be the name of a PSSFSS results file, from which the vector of results will be read. The keyword arguments are used to provide values for the same-named fields in the TEP structure. +## Requirements for TEP File Compatibility +Because a TEP file contains all of the information of the full 4×4 scattering matrix computed by PSSFSS, there are no +limitations on the type of unit cell geometry that can be used for creating TEP files. + +TEP files use the concept of "front" and "rear" incidence. When converting a PSSFSS analysis result to TEP format, Region 1 +(the first layer in the `strata` vector) is taken as the "front" incidence region, and Region `n` (the last layer) is +taken to be the "rear" region. Both of these layers should have zero width and assume vacuum electrical parameters. I.e., +they should be specified as `Layer()` in the `strata` stackup. + `results` (or `resultfile`) must contain the results of a PSSFSS analysis sweep over θ and ϕ (and possibly frequency) such that 1. Incidence angles θ and ϕ rather than incremental phasings ψ₁ and ψ₂ were used. 1. If more than one ϕ value is used, then all ϕ values in the range `0:Δϕ:(360-Δϕ)` must be present. @@ -562,5 +573,181 @@ function res2tep(results::Vector{Result}; name="tep", class="res2tep") end # function +""" + _prepare_results_for_fresnel(results::Vector{Result}) -> Vector{Result} -> newresults::Vector{Result} + +Verify that the input vector of `Result` objects is suitable for creation of an HFSS-compatible Fresnel table. + +The following checks are performed: +1. Ensure that incidence angles rather than incremental phasings are used. +2. Ensure that the θ angles begin at 0 and are uniformly spaced up to the maximum θ value present. +3. Ensure that the increment in θ values divides evenly into 90. +4. Ensure that if multiple frequencies are present, then they have a uniform spacing. + +The output vector is modified in the following ways: +1. Results for all ϕ angles other than the minimum ϕ in absolute value are discarded. +2. The remaining results are sorted into order of increasing θ, then increasing frequency. +3. The vector is sorted into the appropriate order for creating a Fresnel table: frequency + varies most rapidly, then θ. +4. If the input θ values do not extend all the way to 90°, then the results for the maximum supplied θ are copied and appended + a sufficient number of times to simulate data all the way to 90°. +""" +function _prepare_results_for_fresnel(results::Vector{Result}) + kys = keys(results[1].steering) + kys != (:θ, :ϕ) && error("results do not use θ and ϕ for steering") + + ϕmin = minimum(abs, (r.steering.ϕ for r in results)) + results = filter(r -> r.steering.ϕ == ϕmin, results) + + sort!(results, by = r -> (r.steering.θ, r.FGHz)) + + freqs = unique!([r.FGHz for r in results]) + thetas = unique!([r.steering.θ for r in results]) + iszero(first(thetas)) || error("First θ value is not equal to zero") + nf = length(freqs) + if nf > 1 + df = freqs[2] - freqs[1] + frange = first(freqs):df:last(freqs) + freqs == frange || error("frequencies are not uniformly spaced") + end + nt = length(thetas) + nt ≤ 1 && error("Too few θ value") + dtheta = thetas[2] - thetas[1] + isinteger(90 / dtheta) || error("θ increment does not divide evenly into 90") + theta_range = 0:dtheta:last(thetas) + thetas == theta_range || error("θ values are not uniformly spaced") + + # Verify that full cartesian product of angles and frequencies is present: + i = 0 + for t in thetas, f in freqs + i += 1 + r = results[i] + if (r.steering.θ != t) || (r.FGHz != f) + error("Missing case (FGHz, ϕ, θ) = ($f, $phimin, $t) in input vector") + end + end + + # Fill in any missing theta values + if last(thetas) < 90 + newthetas = float((dtheta + last(thetas)):dtheta:90) + irng = (length(results) - nf + 1):length(results) # Indices in results for all frequencies of last θ input value + for newtheta in newthetas, i in irng + r = results[i] + steer = (θ = newtheta, ϕ = r.steering.ϕ) + rnew = Result(r.gsm, steer, r.β⃗₀₀, r.FGHz, r.ϵᵣin, r.μᵣin, r.β₁in, r.β₂in, r.ϵᵣout, r.μᵣout, r.β₁out, r.β₂out) + push!(results, rnew) + end + end + + return results, union(thetas, newthetas), freqs +end # function + + + +""" + res2fresnel(results::Vector{Result}, fresnelfile::AbstractString) + res2fresnel(resultfile::AbstractString, fresnelfile::AbstractString) + +Create an HFSS-compatible "Fresnel table" file from `results`, the vector of `Result` objects returned by +the `analyze` function. If the first positional argument is an `AbstractString`, it is +assumed to be the name of a PSSFSS results file, from which the vector of results will be read. + +Since Fresnel tables contain data for only a single ϕ value, if the input `result` vector contains data for multiple +ϕ values, only the value with minimum magnitude will be used. + +Fresnel tables may be formatted to contain only reflection coefficients (for a so-called "opaque" structure), or they +may contain both reflection and transmission coefficients (a "non-opaque" structure). +An opaque structure is one for which the s21 partition of the generalized scattering matrix is identically zero +for all frequencies and scan angles. The correct format to be written will be selected automatically by `res2fresnel`. + +## Requirements for Fresnel Table Compatibility +The data in `results` must satisfy the following requirements: +1. Incidence angles rather than incremental phasings must be used. +2. θ angles must begin at 0 and be uniformly spaced up to the maximum θ value present. +3. The increment in θ values must divide evenly into 90. +4. If multiple frequencies are present, then they must have a uniform spacing. + +A Fresnel table must contain θ values equally spaced between 0 and 90, inclusive. +If the `results` vector provided as input does not contain θ values all the way to 90, then the scattering matrix values +corresponding to the maximum provided θ value will be copied into the remaining angular "slots" as necessary to provide +a complete Fresnel table. + +There are some limitations on the type of unit cell geometry that should be used for creating Fresnel tables. First, a Fresnel +table contains data for only a single ϕ value. This means that the geometry being analyzed must be such that the scattering +matrix of the structure is essentially independent of ϕ. As a counterexample, a strip grid is not a suitable structure, since +its scattering properties are strongly dependent on ϕ. Second, a Fresnel table records only co-polarized +(TE → TE and TM → TM) transmission and reflection coefficients. This means that the structure being analyzed must not +generate cross-polarized (TE → TM or TM → TE) transmission or reflection coefficients of significant amplitude. + +Fresnel tables consider only incidence from a single "front" region. When creating the Fresnel table, the front region is taken +to be Region 1 of the PSSFSS model (i.e. the first layer present in the PSSFSS `strata` vector). + +### Additional Requirements for Non-Opaque Structures +When used in an HFSS SBR+ model, the scattering properties read from the Fresnel table are applied to a zero-thickness surface, +so that the transmitted ray is launched from the same "hit" point of the surface that was encountered by the incident +ray. Because of this, the phase reference plane for both input and output ports of the PSSFSS model should be located +at this front surface (i.e. the first interface plane in the `strata` vector). This is accomplished by specifying zero +width for the first `Layer` object (i.e. using `Layer()` for the first layer), and then specifying the final layer's width +to be the negative of the sum of all the other layer widths in the `strata` vector. The negative width value shifts +the output port reference plane to coincide with that of the input port. As an example: +```julia +strata = [Layer(), Layer(width=2mm, ϵᵣ=2.2) Layer(width=3.3mm, ϵᵣ=3.0), Layer(width=2mm, ϵᵣ=2.2), Layer(width=-7.3mm)] +``` +""" +function res2fresnel(results::Vector{Result}, fresnelfile::AbstractString) + + results, thetas, freqs = _prepare_results_for_fresnel(results) + nt = length(thetas) + nf = length(freqs) + ronly = all(r -> iszero(r.gsm[2,1]), results) # reflection only + open(fresnelfile, "w") do fid + date, clock = split(string(now()), 'T') + if ronly + println(fid, "# HFSS-compatible Fresnel reflection table created by PSSFSS") + println(fid, "# Created on ", date, " at ", clock) + println(fid, "ReflTab1e") + else + println(fid, "# HFSS-compatible Fresnel reflection/transmission table created by PSSFSS") + println(fid, "# Created on ", date, " at ", clock) + println(fid, "RTTable") + end + + println(fid, "# = - 1") + println(fid, nt - 1) + + if nf == 1 + println(fid, "# Mono freq. table for ", only(freqs), " GHz") + println(fid, "MonoFreq") + println(fid, "# Data section follows.") + else + println(fid, "# MultiFreq ") + println(fid, "MultiFreq ", first(freqs), " ", last(freqs), " ", nf-1) + println(fid, "# Data section follows. Frequency loops within theta") + end + + if ronly + println(fid, "# ") + else + println(fid, "# ") + end + + for res in results + r = res.gsm[1,1]; rte = r[1,1]; rtm = -r[2,2] + t = res.gsm[2,1]; tte = t[1,1]; ttm = t[2,2] + @printf(fid, "%8.5f %8.5f %8.5f %8.5f", real(rte), imag(rte), real(rtm), imag(rtm)) + if ronly + println(fid) + else + @printf(fid, " %8.5f %8.5f %8.5f %8.5f\n", real(tte), imag(tte), real(ttm), imag(ttm)) + end + end + end +end # function + +function res2fresnel(resultfile::AbstractString, fresnelfile::AbstractString) + results = read_result_file(resultfile) + res2fresnel(results, fresnelfile) +end + end # module \ No newline at end of file diff --git a/src/PSSFSS.jl b/src/PSSFSS.jl index b65edd06..27da6025 100644 --- a/src/PSSFSS.jl +++ b/src/PSSFSS.jl @@ -63,9 +63,10 @@ using .Constants: twopi, c₀, tdigits, dbmin using .Log: pssfss_logger, @logfile @reexport using .PSSFSSLen @reexport using .Layers: Layer -@reexport using .Elements: rectstrip, diagstrip, polyring, manji, meander, loadedcross, - jerusalemcross, pecsheet, pmcsheet, sinuous, splitring -@reexport using .Outputs: @outputs, extract_result_file, extract_result, res2tep +@reexport using .Elements: diagstrip, jerusalemcross, loadedcross, manji, meander, + pixels, pecsheet, pmcsheet, polyring, rectstrip, sinuous, + splitring, sympixels +@reexport using .Outputs: @outputs, extract_result_file, extract_result, res2fresnel, res2tep using .Outputs: Result, append_result_data using .FastSweep: interpolate_band @@ -75,14 +76,13 @@ Base.isfile(f::Base.DevNull) = false Base.open(f::Base.DevNull, ::AbstractString) = f """ - result = analyze(strata::Vector, flist, steering; outlist=[], logfile="pssfss.log", resultfile="pssfss.res", - showprogress::Bool=true, fastsweep=true) + result = analyze(strata, flist, steering; outlist=[], logfile="pssfss.log", resultfile="pssfss.res", showprogress::Bool=true, fastsweep=true) Analyze a full FSS/PSS structure over a range of frequencies and steering angles/phasings. Generate output files as specified in `outlist`. ## Positional Arguments -- `strata`: A vector of `Layer` and `Sheet` objects. The first and last entries must be of type `Layer`. +- `strata`: A vector of `Layer` and `RWGSheet` objects. The first and last entries must be of type `Layer`. - `flist`: An iterable containing the analysis frequencies in GHz. @@ -93,7 +93,8 @@ Generate output files as specified in `outlist`. - one of {`:psi1` ,`:ψ₁`} and one of {`:psi2`, `:ψ₂`}. - All steering parameters are input in degrees. + All steering parameters are input in degrees. Examples of valid `steering` tuples: + `(θ=0, ϕ=0)`, `(theta=0:10:60, phi=[0, 45])`, `(theta=40, ϕ=90)`, `(psi1=0, psi2=90)`, `(ψ₁=0, psi2=34.3)`. The program will analyze while iterating over a triple loop over the two steering parameters and frequency, with frequency in the innermost loop (i.e. varying the fastest). @@ -127,8 +128,8 @@ Generate output files as specified in `outlist`. ## Return Value - `result`: A vector of `Result` objects, one for each scan angle/frequency combination. This -vector can be passed as an input to the [`extract_result`](@refs) function to obtain any desired -performance parameters that are supported by the [`@outputs`](@refs) macro. + vector can be passed as an input to the [`extract_result`](@refs) function to obtain any desired + performance parameters that are supported by the [`@outputs`](@refs) macro. """ function analyze(strata::Vector, flist, steering; outlist=[], logfile="pssfss.log", resultfile="pssfss.res", showprogress::Bool=true, fastsweep::Bool=length(flist)≥6) diff --git a/src/RWG.jl b/src/RWG.jl index 949f10fc..a6f0017f 100644 --- a/src/RWG.jl +++ b/src/RWG.jl @@ -305,17 +305,18 @@ function setup_rwg(sheet::RWGSheet, leafsize::Int=9)::RWGData rnm1[1], rnm1[2], rnm2[1], rnm2[2], rnm3[1], rnm3[2]) end - kdtree = KDTree(data, leafsize=leafsize) + kdtree = KDTree(data, leafsize=leafsize, reorder=false) ufp2fp = Vector{Int}[] mn = 0 # Initialize face/pair index. nufp = 0 found = fill(false, nface^2) r = 1e-6 * norm(centroid[1] - centroid[2]) + idxs = Int[] for n in 1:nface, m in 1:nface mn += 1 # Bump global face/pair counter. found[mn] && continue nufp += 1 - idxs = inrange(kdtree, data[mn], r, true) + idxs = inrange(kdtree, data[mn], r, false) found[idxs] .= true push!(ufp2fp, idxs) for i in idxs diff --git a/src/Sheets.jl b/src/Sheets.jl index 01a54d92..72720693 100644 --- a/src/Sheets.jl +++ b/src/Sheets.jl @@ -17,7 +17,13 @@ const SV2 = SVector{2,Float64} abstract type Sheet end - +""" + RWGSheet +A type that represents a zero-thickness sheet of periodically patterned metalization. Particular instances +are created by calling a constructor function for a specific type of sheet geometry. These include: +`diagstrip`, `jerusalemcross`, `loadedcross`, `manji`, `meander`, `pecsheet`, `pmcsheet`, `polyring`, +`rectstrip`, `sinuous`, and `splitring`. +""" mutable struct RWGSheet <: Sheet style::String units::PSSFSSLength # Length unit @@ -351,24 +357,33 @@ function combine(sh1::RWGSheet, sh2::RWGSheet, dup_coor::Char, dup_coor_value::R # Count number of vertices located at the duplicate coordinate. # Save vertex indices of matching points in vcen1 and vcen2 # + Nv1 = length(sh1.ρ) + Ne1 = length(sh1.e1) + Nf1 = size(sh1.fe, 2) + Nv2 = length(sh2.ρ) + Ne2 = length(sh2.e1) + Nf2 = size(sh2.fe, 2) + + # Find coincident vertices Nvcen = 0 - Necen = 0 - vcen1 = Int[] - vcen2 = Int[] + vcen1 = Int[] # coincident vertices in sh1 + vcen2 = Int[] # coincident vertices in sh2 + if dup_coor == ' ' + elseif dup_coor == 'x' + ρind = 1 + elseif dup_coor == 'y' + ρind = 2 + else + throw(ArgumentError("Illegal value for dup_coor: '$dup_coor'")) + end if dup_coor ≠ ' ' tol = 0.5e-4 * norm(sh1.ρ[sh1.e1[1]] - sh1.ρ[sh1.e2[1]]) - for i in 1:length(sh2.ρ) - if dup_coor == 'x' - test = sh2.ρ[i][1] - elseif dup_coor == 'y' - test = sh2.ρ[i][2] - else - test = typemax(typeof(test)) - end + for i in 1:Nv2 + test = sh2.ρ[i][ρind] if abs(test - dup_coor_value) < tol # Test to see if there is a sh1 vertex at same coordinate: n1match = 0 - for n1 in 1:length(sh1.ρ) + for n1 in 1:Nv1 if norm(sh1.ρ[n1] - sh2.ρ[i]) < tol n1match = n1 break @@ -381,23 +396,25 @@ function combine(sh1::RWGSheet, sh2::RWGSheet, dup_coor::Char, dup_coor_value::R end end end - Necen = Nvcen - 1 # Number of shared edges. - ecen1 = zeros(Int, Necen) - ecen2 = zeros(Int, Necen) + # Locate and save indices in sh1 and sh2 of edges along the center line: - Necen = 0 - for i in 1:length(sh2.e1) # Loop over sh2 edges + Necen = 0 # Initialize count of shared center edges + ecen1 = Int[] # Coincident edges in sh1 + ecen2 = Int[] # Coincident edges in sh1 + for i in 1:Ne2 # Loop over sh2 edges if (sh2.e1[i] in vcen2) && (sh2.e2[i] in vcen2) # Edge i of sh2 is a shared edge. Necen += 1 - ecen2[Necen] = i - i1 = findfirst(x -> x == sh2.e1[i], vcen2) - i2 = findfirst(x -> x == sh2.e2[i], vcen2) + push!(ecen2, i) + i1 = findfirst(==(sh2.e1[i]), vcen2) + i2 = findfirst(==(sh2.e2[i]), vcen2) + # vcen1[k] and vcen2[k] refer to nodes in sh1 and sh2, resp, but at the same coordinates + vcen1nodes = extrema((vcen1[i1], vcen1[i2])) # Same node locations but wrt sh1 # Now find matching sh1 edge: - for j1 in 1:length(sh1.e1) - if ((sh1.e1[j1] == vcen1[i1]) && (sh1.e2[j1] == vcen1[i2])) || - ((sh1.e2[j1] == vcen1[i1]) && (sh1.e1[j1] == vcen1[i2])) - ecen1[Necen] = j1 + for j1 in 1:Ne1 + if extrema((sh1.e1[j1], sh1.e2[j1])) == vcen1nodes + push!(ecen1, j1) + # ecen1[k] and ecen2[k] refer to edges in sh1 and sh2, resp, but at the same locations @goto sh2edges end end @@ -411,96 +428,104 @@ function combine(sh1::RWGSheet, sh2::RWGSheet, dup_coor::Char, dup_coor_value::R # Allocate triangulation arrays in new sheet: sh3 = RWGSheet() - sh3.e1 = zeros(Int, length(sh1.e1) + length(sh2.e1) - Necen) - sh3.e2 = zeros(Int, length(sh1.e2) + length(sh2.e2) - Necen) - sh3.ρ = Vector{SV2}(undef, length(sh1.ρ) + length(sh2.ρ) - Nvcen) - sh3.fe = zeros(Int, 3, size(sh1.fe, 2) + size(sh2.fe, 2)) - sh3.fv = zeros(Int, 3, size(sh1.fv, 2) + size(sh2.fv, 2)) + sh3.e1 = zeros(Int, Ne1 + Ne2) # Will remove Necen of these later + sh3.e2 = zeros(Int, Ne1 + Ne2) # Will remove Necen of these later + sh3.ρ = Vector{SV2}(undef, Nv1 + Nv2 - Nvcen) + sh3.fe = zeros(Int, 3, Nf1 + Nf2) + sh3.fv = zeros(Int, 3, Nf1 + Nf2) # Copy vertex locations: - sh3.ρ[1:length(sh1.ρ)] = sh1.ρ + sh3.ρ[1:Nv1] = sh1.ρ if dup_coor == ' ' - sh3.ρ[(1+length(sh1.ρ)):end] = sh2.ρ + sh3.ρ[(1+Nv1):end] = sh2.ρ else - i2 = 1 + length(sh1.ρ) # Node counter - for i in 1:length(sh2.ρ) - if !(i in vcen2[1:Nvcen]) + i2 = 1 + Nv1 # Node counter + for i in 1:Nv2 + if !(i in vcen2) sh3.ρ[i2] = sh2.ρ[i] i2 += 1 end end end + # Copy edge-node matrices - eoffset = length(sh1.e1) - voffset = length(sh1.ρ) - sh3.e1[1:length(sh1.e1)] = sh1.e1 - sh3.e2[1:length(sh1.e2)] = sh1.e2 - if dup_coor == ' ' - sh3.e1[(eoffset+1):end] = sh2.e1 .+ voffset - sh3.e2[(eoffset+1):end] = sh2.e2 .+ voffset - else - for i in 1:length(sh2.e1) - if i in ecen2 - # Edge i of sh2 is located on duplication line. - # Do not include this duplicate edge in sh3, but - # decrement the edge index offset: - eoffset -= 1 - else - if sh2.e1[i] in vcen2 - # The initial point of edge i of sh2 is a duplicate vertex. - i2 = findfirst(x -> x == sh2.e1[i], vcen2) - sh3.e1[i+eoffset] = vcen1[i2] + eoffset = Ne1 + voffset = Nv1 + sh3.e1[1:Ne1] = sh1.e1 + sh3.e2[1:Ne1] = sh1.e2 + sh3.e1[Ne1+1:end] .= sh2.e1 .+ voffset # Offset will be corrected later if necessary + sh3.e2[Ne1+1:end] .= sh2.e2 .+ voffset # Offset will be corrected later if necessary + + # Correct if there are duplicate edges + sh3e1new = @view sh3.e1[Ne1+1:end] + sh3e2new = @view sh3.e2[Ne1+1:end] + if dup_coor ≠ ' ' + for i in eachindex(sh2.e1, sh2.e2, sh3e1new, sh3e2new) + if i ∉ ecen2 # Edges in ecen2 will be removed later + i2 = findfirst(==(sh2.e1[i]), vcen2) + if isnothing(i2) + # ordinary point not on shared border: + sh3e1new[i] -= count(<(sh2.e1[i]), vcen2) else - # Ordinary point - sh3.e1[i+eoffset] = sh2.e1[i] + voffset + # initial point of edge i of sh2 is a duplicate vertex: + sh3e1new[i] = vcen1[i2] end # - if sh2.e2[i] in vcen2 - # The terminal point of edge i of sh2 is on the duplication edge. - i2 = findfirst(x -> x == sh2.e2[i], vcen2) - sh3.e2[i+eoffset] = vcen1[i2] + i2 = findfirst(==(sh2.e2[i]), vcen2) + if isnothing(i2) + # ordinary point not on shared border: + sh3e2new[i] -= count(<(sh2.e2[i]), vcen2) else - # Ordinary point: - sh3.e2[i+eoffset] = sh2.e2[i] + voffset + # initial point of edge i of sh2 is a duplicate vertex: + sh3e2new[i] = vcen1[i2] end end end - # Correct vertex indices: - for i in Nvcen:-1:1 - sh3.e1[sh3.e1.>(length(sh1.ρ)+vcen2[i])] .-= 1 - sh3.e2[sh3.e2.>(length(sh1.ρ)+vcen2[i])] .-= 1 - end + # Remove duplicate edges + delindices = Ne1 .+ ecen2 + deleteat!(sh3.e1, delindices) + deleteat!(sh3.e2, delindices) end + # Copy face/vertex matrix - sh3.fv[:, 1:size(sh1.fv, 2)] = sh1.fv - sh3.fv[:, 1+size(sh1.fv, 2):end] = sh2.fv .+ voffset # offset will be corrected later - # Correct duplicate vertices from sh2: - for n2 in 1:Nvcen # Examine each duplicate vertex in sh2 - n1 = vcen1[n2] # Initialize index of matching vertex in sh1 - # Replace all references to vertex vcen2[n2] with ref to vcen1[n1] - sh3.fv[sh3.fv.==vcen2[n2]+voffset] .= n1 - end - # Correct vertex indices: - for i in Nvcen:-1:1 - sh3.fv[sh3.fv.>voffset+vcen2[i]] .-= 1 + sh3.fv[:, 1:Nf1] = sh1.fv + sh3.fv[:, 1+Nf1:end] = sh2.fv .+ voffset # offset will be corrected later + if dup_coor ≠ ' ' + # Correct duplicate vertices from sh2: + sh3fvnew = @view sh3.fv[:, 1+Nf1:end] + for n in eachindex(sh3fvnew, sh2.fv) + v2 = sh2.fv[n] + ncen = findfirst(==(v2), vcen2) + if isnothing(ncen) + # ordinary, non-duplicate point + sh3fvnew[n] -= count(<(v2), vcen2) # correct offset + else + sh3fvnew[n] = vcen1[ncen] # Replace duplicate node + end + end end + # Copy face/edge matrix: - foffset = size(sh1.fv, 2) eoffset = length(sh1.e1) - sh3.fe[:, 1:size(sh1.fe, 2)] = sh1.fe - sh3.fe[:, 1+size(sh1.fe, 2):end] = sh2.fe .+ eoffset # offset will be corrected later + sh3.fe[:, 1:Nf1] = sh1.fe + sh3.fe[:, 1+Nf1:end] .= sh2.fe .+ eoffset # offset will be corrected later # Correct duplicate edges from sh2: if dup_coor ≠ ' ' - for n2 in 1:length(ecen2) # Examine each duplicate edge in sh2 - e2 = ecen2[n2] # Index of duplicate edge in sh2. - e1 = ecen1[n2] # Index of duplicate edge in sh1. - sh3.fe[sh3.fe.==(e2+eoffset)] .= e1 - end - # Correct edge indices: - for i in length(ecen2):-1:1 - sh3.fe[sh3.fe.>eoffset+ecen2[i]] .-= 1 + sh3fenew = @view sh3.fe[:, 1+Nf1:end] + for n in eachindex(sh2.fe, sh3fenew) + e2 = sh2.fe[n] + ncen = findfirst(==(e2), ecen2) + if isnothing(ncen) + # ordinary (non-duplicate) edge + sh3fenew[n] -= count(<(e2), ecen2) # correct the offset + else + sh3fenew[n] = ecen1[ncen] + end end end + + test_fefv(sh3) + return sh3 end @@ -771,4 +796,31 @@ end end end +function test_fefv(sheet::RWGSheet) + ip1 = (2, 3, 1) + ip2 = (3, 1, 2) + for iface in axes(sheet.fv, 2) + for ivlocal in axes(sheet.fv, 1) # local node index within current face + ivglobal = sheet.fv[ivlocal, iface] # Global node index + oppedge= sheet.fe[ivlocal, iface] # Should be global index of edge opposite this vertex + edgen1n2 = Set((sheet.e1[oppedge], sheet.e2[oppedge])) # global node indices of edge vertices + facen1n2 = Set((sheet.fv[ip1[ivlocal], iface], sheet.fv[ip2[ivlocal], iface])) # global node indices of other face vertices + if edgen1n2 ≠ facen1n2 + @show edgen1n2 + @show facen1n2 + println("Bad triangle fe/fv matching for face #", iface, " local vertex ", ivlocal) + println("face edges: " ) + for e in @view sheet.fe[:,iface] + println(" edge ", e, " from node ", sheet.e1[e], " to node ", sheet.e2[e]) + end + print("face nodes: ") + foreach(x -> print(x," "), @view sheet.fv[:,iface]) + println() + error("bad") + end + end + end # iface loop +end + + end # module diff --git a/src/pixels.jl b/src/pixels.jl new file mode 100644 index 00000000..a6982353 --- /dev/null +++ b/src/pixels.jl @@ -0,0 +1,328 @@ +# Note: this file is included by Elements.jl + +""" + SymPixels + +A struct containing the pixel parameters of a symmetric pixel mesh element. Fields: + - `halfnint`: Half the number of pixels across the interior (non-rim) portion of + the unit cell. Must be a positive integer. + - `nrim`: Width of the rim in pixels. Must be a nonnegative integer. +""" +struct SymPixels{T1<:Integer, T2<:Integer} + halfnint::T1 + nrim::T2 # Width of the rim in pixels ≥ 0 + function SymPixels(halfint::T1, nrim::T2) where {T1<:Integer, T2<:Integer} + halfint < 1 && throw(ArgumentError("SymPixels first argument must be ≥ 1")) + nrim < 0 && throw(ArgumentError("SymPixels second argument must be ≥ 0")) + return new{T1, T2}(halfint, nrim) + end +end + +SymPixels(; halfnint::Integer, nrim::Integer=0) = SymPixels(halfnint, nrim) # keyword constructor + +ucsidelen(s::SymPixels) = 2 * (s.halfnint + s.nrim) # Total number of pixels along side of unit cell + +""" + patternveclen(s::SymPixels) -> vlen::Integer + + +Return the length of the 1/0 pattern vector needed to define the unit cell metallization pattern +for a symmetric pixel mesh. This is the number of pixels in the irreducible zone for the mesh. +""" +patternveclen(s::SymPixels) = s.halfnint * (1 + s.halfnint) ÷ 2 + +""" + sympixmat(s::SymPixels, v::AbstractVector) -> mat::Matrix{Bool} + +Compute the metalization pixel pattern matrix for a symmetric pixel mesh. + +The matrix will be square of dimension `ucsidelen(s)`. `true` elements correspond +to metalized pixels in the unit cell. Element (i,j) of the matrix corresponds +to to the pixel in the unit cell located at the i'th x location and j'th y location, +counting from the lower-left corner of the unit cell. + +## Input Arguments +- `s::SymPixels`: Defines the pixel parameters of the mesh. +- `v::AbstractVector`: This is the metallization pattern vector. It must contain only ones and + zeros, and its length should be equal to the number of pixels in the + irreducible zone of the mesh. A greater length is also allowed, in which case a warning is + issued and the extra elements in `v` are ignored. + +""" +function sympixmat(s::SymPixels, v::AbstractVector) + niz = patternveclen(s) + length(v) ≥ niz || throw(ArgumentError("v is not long enough")) + all(x -> isone(x) || iszero(x), v) || throw(ArgumentError("v must consist only of 1s and 0s")) + length(v) > niz && @warn "v is longer than necessary. Ignoring extra entries" + + totsidelen = ucsidelen(s) + mat = falses(totsidelen, totsidelen) + + # Fill the rim + nrim = s.nrim + if nrim > 0 + mat[1:nrim, :] .= true + mat[:, 1:nrim] .= true + mat[:, end-nrim+1:end] .= true + mat[end-nrim+1:end, :] .= true + end + + nir = patternveclen(s) # Number of pixels in irreducible zone + halfnint = s.halfnint + nint = 2 * halfnint + + # Work on the interior + mint = @view mat[1+nrim:end-nrim, 1+nrim:end-nrim] # interior matrix + + # Fill upper right quadrant + mintur = @view mint[1:halfnint, halfnint+1:end] # upper-right quadrant of interior matrix + for (i,j,k) in zip(1:nir, ulindices(halfnint), lrindices(halfnint)) + mintur[j] = mintur[k] = v[i] + end + # Copy it into the upper left quadrant + mintul = @view mint[1:halfnint, halfnint:-1:1] # Reverse the columns + for col in 1:s.halfnint + mintul[:, col] .= @view mintur[:, col] + end + # Copy into the bottom half: + mintbot = @view mint[nint:-1:halfnint+1, :] # reverse the rows + for row in 1:halfnint + mintbot[row, :] .= @view mint[row, :] + end + + return mat +end + +""" + ulindices(n::Integer) + +Return a CartesianIndex generator that iterates over the indices in the upper left half +of a square matrix of order n. The upper left half is bounded by (and includes) the +antidiagonal of the matrix. The elements are iterated in column major order. +""" +function ulindices(n::Integer) + n > 0 || throw(ArgumentError("uladindices requires a a positive argument")) + return (CartesianIndex(i,j) for (j,k) in zip(1:n, n:-1:1) for i in 1:n if i ≤ k) +end + +""" + lrindices(n::Integer) + +Return a CartesianIndex generator that iterates over the indices in the lower right half +of a square matrix of order n. The lower right half is bounded by (and includes) the +antidiagonal of the matrix. The elements are iterated in row major order, but with the +elements of each row reversed. This places the elements in 1:1 correspondence with their +mirror images in the upper left half. +""" +function lrindices(n::Integer) + n > 0 || throw(ArgumentError("lradindices requires a a positive argument")) + return (CartesianIndex(n-j+1, n-i+1) for (j,k) in zip(1:n, n:-1:1) for i in 1:n if i ≤ k) +end + +const keepstart = first(findfirst("- `dx", optional_kwargs)) - 1 +const pix_optional_kwargs = +""" +- `class::Char='M'` Specify the class, either `'J'` or `'M'`. If `'J'`, the unknowns are electric surface + currents in the areas corresponding to `1` values of the pixels. If `'M'`, the unknowns are magnetic surface + currents in the areas corresponding to `0` values of the pixels. It is known that using `'J'` can result in + grossly incorrect results for some geometries where adjacent metallic pixels intersect at only a single point. + Therefore, use of only `'M'` is strongly recommended for most `pixels` elements and all `sympixels` elements. +""" * +optional_kwargs[keepstart:end] + + +""" + function sympixels(; P, nrim, halfnint, patternvec, units, kwargs...) + +# Description: + +Create a variable of type `RWGSheet` that contains the triangulation for a symmetrically pixelated +square unit cell. The pattern of metallic pixels has C4 symmetry, as well as left-right and up-down mirror symmetry, +as well as each quadrant exhibiting antidiagonal mirror symmetry. The returned value has fields `s₁`, `s₂`, `β₁`, +`β₂`, `ρ`, `e1`, `e2`, `fv`, `fe`, and `fr` properly initialized. The pixels included in the triangulation are +determined by the `patternvec` input vector as described below. + + +# Arguments: + +All arguments are keyword arguments which can be entered in any order. + +## Required arguments: +- `P::Real > 0`: The side length of the square unit cell, specified in units defined by the `units` keyword argument. +- `nrim::Integer`: The width of the solid metallic rim placed just inside the unit cell boundary, in units of pixels. + A value of `0` implies that there is no rim. +- `halfnint`: Half the number of pixels in the side length for the interior (strictly inside the rim) square pixelated + region of the unit cell. +- `patternvec::AbstractVector{<:Integer}`: A vector of length `halfnint*(halfnint+1)÷2`, consisting solely of 1's and 0's. + The elements of this vector are mapped to pixels in the irreducible zone of the unit cell as shown in the + diagram at ![https://simonp0420.github.io/PSSFSS.jl/stable/assets/sympixel_with_irzone_numbering.png]\ + (https://simonp0420.github.io/PSSFSS.jl/stable/assets/sympixel_with_irzone_numbering.png). + Within the irreducible zone, pixels corresponding to a value of `1` (or `true`) + are taken to be areas of metallization, while `0` or `false` values are metal-free (void) areas. This holds for either + `J` or `M` as the `class` value (see the `class` argument below for important limitations). +- `units`: Length units for `P` (either `mm`, `cm`, `inch`, or `mil`). + +## Optional arguments: +- `pdiv::Int = 1`: The number of "chops" or subdivisions applied to each square pixel side when forming the triangulation. + A value of `1` (the default) means that the pixels included in the triangulation (`1` or `true` values for `class='J'`, + `0` or `false` values for `class='M'`) are not subdivided any further, except for a single diagonal across each square + pixel to form triangles. A value of `n>1` means that each square pixel is first divided into `n×n` square + subpixels, after which a single diagonal edge (if `quad=false`) is added to each subpixel to form triangles. +- `quad::Bool=false`: If `true`, each subpixel (or pixel, if `pdiv` is 1) is divided into four triangles by adding + two diagonals. If `false` (the default), only a single diagonal is added to each square to produce two triangles. +- `sym::Bool = false`: If true, the diagonals added to the squares will exhibit the same left-right and up-down + mirror symmetry as the collection of `true` (`false`) pixel locations. If `false` and `quad=false`, then all added diagonals + across the unit cell will have the same orientation. `sym` has no effect (i.e. is redundant) if `quad=true` since + in that case the two added diagonals already ensure mirror symmetry of the triangulation. +$(pix_optional_kwargs) +""" +function sympixels(; P::Real, nrim::Integer, halfnint::Integer, class::Char='M', + patternvec::AbstractVector{<:Integer}, units, pdiv::Integer=1, quad::Bool=false, + kwargs...)::RWGSheet + + @testpos(P) + @testpos(halfnint) + @testpos(pdiv) + nrim ≥ 0 || throw(ArgumentError("nrim must be nonnegative")) + length(patternvec) ≥ halfnint * (halfnint + 1) ÷ 2 || throw(ArgumentError("patternvec is not long enough")) + all(x -> isone(x) || iszero(x), patternvec) || throw(ArgumentError("patternvec must consist of ones and zeros")) + + s = SymPixels(halfnint, nrim) + patternmat = sympixmat(s, patternvec) + + sheet = pixels(; P, patternmat, units, pdiv, quad, class, kwargs...) + return sheet + +end + +""" + pixels(; P, patternmat, units, kwargs...) + +# Description: + +Create a variable of type `RWGSheet` that contains the triangulation for an arbitrarily pixelated +square unit cell. The returned value has fields `s₁`, `s₂`, `β₁`, `β₂`, `ρ`, `e1`, +`e2`, `fv`, `fe`, and `fr` properly initialized. + + +# Arguments: + +All arguments are keyword arguments which can be entered in any order. + +## Required arguments: +- `P::Real > 0`: The side length of the square unit cell, specified in units defined by the `units` keyword argument. +- `patternmat::AbstractMatrix{<:Integer}`: A square matrix, consisting solely of 1's and 0's. The matrix entries control + the metallization pattern in the unit cell, with a 1 or `true` value denoting a metallized pixel, and a 0 or `false` value + indicating no metallization. The `(i,j)` entry corresponds to the pixel centered at `(x,y) = ((j-1/2)d, P-(i-1/2)d)`, where + `d = P / size(patternmat, 1)`. +- `units`: Length units for `P` (either `mm`, `cm`, `inch`, or `mil`). + +## Optional arguments: +- `pdiv::Int = 1`: The number of "chops" or subdivisions applied to each square pixel side when forming the geometry triangulation. + A value of `1` (the default) means that the pixels are not subdivided any further, except for a single diagonal across each + pixel to form triangles. A value of `2` would subdivide each pixel into 4 squares. A diagonal edge would be added to + each of the resulting squares to form triangles. +- `quad::Bool=false`: If `true`, each subpixel (or pixel, if `pdiv` is 1) is divided into four triangles by adding + two diagonals. If `false` (the default), only a single diagonal is added to produce two triangles. +- `sym::Bool = false`: A `true` value states that the pixel matrix has vertical and horizontal mirror symmetry, and + consists of an even number of rows (and columns). If `true`, the diagonals added to the squares to form + triangles will exhibit the same left-right and up-down mirror symmetry as the collection of `true` (`false`) pixel + locations. If `sym=false` and `quad=false`, then all added diagonals across the unit cell will have the same + orientation. `sym` has no effect (i.e. is redundant) if `quad=true` since + in that case the two added diagonals already ensure mirror symmetry of the triangulation. +$(pix_optional_kwargs) +""" +function pixels(; P::Real, patternmat::AbstractMatrix{<:Integer}, pdiv::Integer=1, class::Char='M', + sym::Bool=false, quad::Bool=false, units, kwarg...) + + kwargs = Dict{Symbol,Any}(kwarg) + haskey(kwargs, :fufp) || (kwargs[:fufp] = true) + kwargs[:class] = class + check_optional_kw_arguments!(kwargs) + @testpos(P) + @testpos(pdiv) + all(x -> isone(x) || iszero(x), patternmat) || throw(ArgumentError("patternmat must consist of ones and zeros")) + size(patternmat,1) == size(patternmat,2) || throw(ArgumentError("patternmat must be a square matrix")) + + if kwargs[:class] == 'M' + patternmat = (!).(patternmat) # Triangulate the non-metal regions + else + @warn """ + + `class='J'` detected for `sympixels` or `pixels` element. + This is is known to produce incorrect results for certain geometries containing metallic islands that intersect in a single point. + It is strongly recommended to use `class='M'` for most `sympixels` and `pixels` elements. + """ + end + + s1 = [P, 0.0] + s2 = [0.0, P] + + d = P / size(patternmat, 1) # Pixel side length + npside = round(Int, P / d) # pixels on side of unit cell + + xrequired = range(start = 0.0, stop = P, length = 1 + pdiv * npside) |> collect + yrequired = copy(xrequired) + + function is_inside(x::Real, y::Real) + # predicate to determine if a point is within the region to be triangulated + (x < 0 || x > P || y < 0 || y > P) && return false + # column index into patternmat: + if iszero(x) + j = 1 + elseif x == P + j = npside + else + j = 1 + trunc(Int, x / d) + end + # row index into patternmat: + if iszero(y) + i = npside + elseif y == P + i = 1 + else + i = 1 + trunc(Int, (P - y) / d) + end + return isone(patternmat[i,j]) ? true : false + end + + # Triangulate + npixtri = count(isone, patternmat) # Number of triangulated pixels + areat = npixtri * d^2 # Total area to triangulate + ntri = npixtri * pdiv^2 * 2 # Number of triangles expected + ntri_requested = ntri ÷ 4 # Much smaller to ensure no extra triangles are added by make_plaid_mesh + if sym && !quad + xrequired .-= xrequired[(length(xrequired)+1) ÷ 2] + yrequired .-= yrequired[(length(yrequired)+1) ÷ 2] + #xrequired = xrequired[((length(xrequired)+1) ÷ 2):end] + #yrequired = yrequired[((length(xrequired)+1) ÷ 2):end] + is_inside_sym(x,y) = is_inside(x + P/2, y + P/2) + sheet = make_sym_plaid_mesh(xrequired, yrequired, areat, ntri_requested, is_inside_sym) + else + sheet = make_plaid_mesh(xrequired, yrequired, areat, ntri_requested, is_inside, quad) + end + + sheet.Zs = kwargs[:Zsheet] + sheet.σ = kwargs[:σ] + sheet.Rq = kwargs[:Rq] + sheet.disttype = kwargs[:disttype] + + # Handle remaining optional arguments + sheet.fufp = kwargs[:fufp] + sheet.class = kwargs[:class] + rotate!(sheet, kwargs[:rot]) + dxdy = SV2([kwargs[:dx], kwargs[:dy]]) + if dxdy ≠ [0.0, 0.0] + sheet.ρ .= (dxdy + xy for xy in sheet.ρ) + end + + sheet.style = "pixels" + sheet.ξη_check = true + sheet.units = units + sheet.s₁ = SV2(s1) + sheet.s₂ = SV2(s2) + sheet.β₁, sheet.β₂ = s₁s₂2β₁β₂(sheet.s₁, sheet.s₂) + + return sheet + +end # function pixels diff --git a/src/structuredtri.jl b/src/structuredtri.jl index 89ef7bde..a7d16770 100644 --- a/src/structuredtri.jl +++ b/src/structuredtri.jl @@ -1,3 +1,5 @@ +# Note: this file is included by Elements.jl. Thus the definitions herein are part of the Elements module. + function loadedcross_structured(; s1::Vector{<:Real}, s2::Vector{<:Real}, L1::Real, L2::Real, w::Real, ntri::Int, orient::Real=0.0, units::PSSFSSLength, kwarg...) kwargs = Dict{Symbol,Any}(kwarg) @@ -60,7 +62,7 @@ function loadedcross_structured(; s1::Vector{<:Real}, s2::Vector{<:Real}, L1::Re end sheet.style = "loadedcross" - sheet.ξη_check = false + sheet.ξη_check = true sheet.units = units sheet.s₁ = s1 sheet.s₂ = s2 @@ -173,7 +175,7 @@ function jerusalemcross_structured(; P::Real, L1::Real, L2::Real, A::Real, B::Re end sheet.style = "jerusalemcross" - sheet.ξη_check = false + sheet.ξη_check = true sheet.units = units sheet.s₁ = SV2([P, 0]) sheet.s₂ = SV2([0, P]) @@ -471,16 +473,20 @@ function manji(; L1::Real, L2::Real, L3::Real, L4::Real=0.0, a::Real=0.0, w::Rea end # function manji """ - make_plaid_mesh(xr, yr, area, ntri, is_inside) -> sheet::RWGSheet + make_plaid_mesh(xr, yr, area, ntri, is_inside, quad=false) -> sheet::RWGSheet Generate a structured, plaid triangular mesh from list of required coordinates and predicate function -# Input Arguments +# Positional Input Arguments - `xr`, `yr`: Vectors of required x and y coordinates for vertices of the geometry to be meshed. - `area`: The area of the geometry to be meshed. - `ntri`: The desired number of triangles for the area to be meshed. - `is_inside`: A predicate function where `is_inside(x,y) -> tf::Bool` determines whether a point (x,y) is within the region to be meshed. +- `quad::Bool=false`: If `true`, each subpixel (or pixel, if `pdiv` is 1) is divided into four triangles by adding + two diagonals. If `false` (the default), only a single diagonal is added to produce two triangles. + + # Return value: - `sheet`: A variable of type RWGSheet with fields ρ, e1, e2, fe, and fv properly initialized. The @@ -488,14 +494,13 @@ Generate a structured, plaid triangular mesh from list of required coordinates a product of `xr` and `yr`, the latter supplemented with additional points to refine the mesh, and then converted to a triangular tesselation by adding a diagonal to each rectangle. """ -function make_plaid_mesh(xr::AbstractVector, yr::AbstractVector, area, ntri, is_inside)::RWGSheet +function make_plaid_mesh(xr::AbstractVector, yr::AbstractVector, area, ntri, is_inside, quad::Bool=false)::RWGSheet xr, yr = sort.((xr, yr)) bigarea = (xr[end] - xr[begin]) * (yr[end] - yr[begin]) # area of circumscribing rectangle bignsq = ceil(Int, bigarea / area * ntri/2) # desired number of squares to form in circumscribing rectangle s = sqrt(bigarea / bignsq) # ideal side length for squares used to tesselate the big area - facevs = Tuple{Tuple{Int,Int}, Tuple{Int,Int}, Tuple{Int,Int}}[] - edgevs = Tuple{Tuple{Int,Int}, Tuple{Int,Int}}[] + # Add new vertex locations as needed to generate at least desired number of triangles: xn = xr[begin:begin] yn = yr[begin:begin] for (tr, tn) in ((xr, xn), (yr, yn)), i in eachindex(tr)[begin+1:end] @@ -506,6 +511,14 @@ function make_plaid_mesh(xr::AbstractVector, yr::AbstractVector, area, ntri, is_ end # xn and yn now contain the plaid vertex coordinates + if quad + # Add centerpoint values + xnc = [0.5(xn[i] + xn[i+1]) for i in 1:(length(xn)-1)] + xn = sort!(append!(xn, xnc)) + ync = [0.5(yn[i] + yn[i+1]) for i in 1:(length(yn)-1)] + yn = sort!(append!(yn, ync)) + end + # Initialize vectors of face and edge indices into xn and yn: facevs = Tuple{Tuple{Int,Int}, Tuple{Int,Int}, Tuple{Int,Int}}[] edgevs = Tuple{Tuple{Int,Int}, Tuple{Int,Int}}[] @@ -513,19 +526,37 @@ function make_plaid_mesh(xr::AbstractVector, yr::AbstractVector, area, ntri, is_ # Add triangular faces and edges within the desired geometry. Assumes that # if the center of the rectangle is inside, then so are both triangles # partitioning that rectangle. - for i in eachindex(xn)[begin+1:end] - xcen = 0.5 * (xn[i] + xn[i-1]) - for j in eachindex(yn)[begin+1:end] - ycen = 0.5 * (yn[j] + yn[j-1]) + irange = quad ? (3:2:length(xn)) : (2:length(xn)) + jrange = quad ? (3:2:length(yn)) : (2:length(yn)) + for i in irange + xcen = quad ? xn[i-1] : 0.5 * (xn[i] + xn[i-1]) + for j in jrange + ycen = quad ? yn[j-1] : 0.5 * (yn[j] + yn[j-1]) is_inside(xcen, ycen) || continue - topface = ((i,j-1), (i,j), (i-1,j)) - botface = ((i,j-1), (i-1,j), (i-1,j-1)) - push!(facevs, topface, botface) - push!(edgevs, ((i-1,j-1), (i,j-1)), - ((i,j-1), (i,j)), - ((i,j), (i-1,j)), - ((i-1,j),(i-1,j-1)), - ((i,j-1), (i-1,j))) + if quad + leftface = ((i-1, j-1), (i-2, j), (i-2, j-2)) + rightface = ((i-1, j-1), (i, j-2), (i, j)) + topface = ((i-1, j-1), (i,j), (i-2,j)) + botface = ((i-1, j-1), (i-2, j-2), (i, j-2)) + push!(facevs, leftface, rightface, topface, botface) + push!(edgevs, ((i-1, j-1), (i-2, j)), + ((i-2, j), (i-2, j-2)), + ((i-2, j-2), (i-1, j-1)), + ((i-1, j-1), (i, j-2)), + ((i, j-2), (i, j)), + ((i, j), (i-1, j-1)), + ((i-2, j), (i, j)), + ((i-2, j-2), (i, j-2))) + else + topface = ((i,j-1), (i,j), (i-1,j)) + botface = ((i,j-1), (i-1,j), (i-1,j-1)) + push!(facevs, topface, botface) + push!(edgevs, ((i-1,j-1), (i,j-1)), + ((i,j-1), (i,j)), + ((i,j), (i-1,j)), + ((i-1,j),(i-1,j-1)), + ((i,j-1), (i-1,j))) + end end end @@ -556,6 +587,8 @@ function make_plaid_mesh(xr::AbstractVector, yr::AbstractVector, area, ntri, is_ return sh end + + _plain_rim_area(P, w) = 4 * (P * w - w^2) function _fancy_rim_area(P, w, c) qarea = 3 * _plain_rim_area(c, w) # corner squares @@ -633,3 +666,155 @@ function _squarerim(P::Real, w::Real, c::Real, ntri::Integer) return sheet end + +""" + make_sym_plaid_mesh(xr, yr, area, ntri, is_inside) -> sheet::RWGSheet + +Generate a structured, symmetrical, plaid triangular mesh from list of required coordinates and predicate function. + + +# Input Arguments +- `xr`, `yr`: Vectors of required x and y coordinates for vertices of the geometry to be meshed. + Both of these must contain points distributed symmetrically about the zero, and must also contain + zero. Therefore, each must contain an odd number of elements. +- `area`: The area of the geometry to be meshed. +- `ntri`: The desired number of triangles for the area to be meshed. +- `is_inside`: A predicate function where `is_inside(x,y) -> tf::Bool` determines whether a point (x,y) + is within the region to be meshed. + +# Return value: +- `sheet`: A variable of type RWGSheet with fields ρ, e1, e2, fe, and fv properly initialized. The + mesh results from a plaid rectangular tesselation containing at least the vertices in the Cartesian + product of `xr` and `yr`, the latter supplemented with additional points to refine the mesh, and then + converted to a triangular tesselation by adding a diagonal to each rectangle. +""" +function make_sym_plaid_mesh(xr::AbstractVector, yr::AbstractVector, area, ntri, is_inside)::RWGSheet + nx = length(xr) + ny = length(yr) + isodd(nx) || error("xr must contain an odd number of points") + isodd(ny) || error("yr must contain an odd number of points") + xr, yr = sort.((xr, yr)) + nxlh = (nx - 1) ÷ 2 + nylh = (ny - 1) ÷ 2 + for i in 1:nxlh + xr[i] ≈ -xr[nx - i + 1] || error("xr must by symmetrically disposed about 0") + end + for i in 1:nylh + yr[i] ≈ -yr[nx - i + 1] || error("yr must by symmetrically disposed about 0") + end + iszero(xr[nxlh+1]) || error("xr must contain 0") + iszero(yr[nylh+1]) || error("yr must contain 0") + + + bigarea = (xr[end] - xr[begin]) * (yr[end] - yr[begin]) # area of circumscribing rectangle + bignsq = ceil(Int, bigarea / area * ntri/2) # desired number of squares to form in circumscribing rectangle + s = sqrt(bigarea / bignsq) # ideal side length for squares used to tesselate the big area + + # Restrict consideration to 1st quadrant for now: + xr = xr[nxlh+1:end] + yr = yr[nxlh+1:end] + # Add new vertex locations in 1st quadrant as needed to generate at least desired number of triangles: + xn = xr[begin:begin] + yn = yr[begin:begin] + for (tr, tn) in ((xr, xn), (yr, yn)), i in eachindex(tr)[begin+1:end] + dt = tr[i] - tr[i-1] + nt = max(1, round(Int, dt / s)) + append!(tn, tr[i-1] .+ collect((1:nt) * (dt / nt))) + tn[end] = tr[i] # correct rounding errors + end + + # xn and yn now contain the plaid vertex coordinates + # Initialize vectors of face and edge indices into xn and yn: + facevs = Tuple{Tuple{Int,Int}, Tuple{Int,Int}, Tuple{Int,Int}}[] + edgevs = Tuple{Tuple{Int,Int}, Tuple{Int,Int}}[] + + # Add triangular faces and edges within the desired geometry. Assumes that + # if the center of the rectangle is inside, then so are both triangles + # partitioning that rectangle. + for i in eachindex(xn)[begin+1:end] + xcen = 0.5 * (xn[i] + xn[i-1]) + for j in eachindex(yn)[begin+1:end] + ycen = 0.5 * (yn[j] + yn[j-1]) + is_inside(xcen, ycen) || continue + topface = ((i,j-1), (i,j), (i-1,j)) + botface = ((i,j-1), (i-1,j), (i-1,j-1)) + push!(facevs, topface, botface) + push!(edgevs, ((i-1,j-1), (i,j-1)), + ((i,j-1), (i,j)), + ((i,j), (i-1,j)), + ((i-1,j),(i-1,j-1)), + ((i,j-1), (i-1,j))) + end + end + + nodes = unique!([s[i] for s in facevs for i in 1:3]) # list of node (ix,iy) values + nface = length(facevs) + inodes = Dict(t => i for (i,t) in pairs(nodes)) # returns linear node index given (ix,iy) + edgenodes = unique!([extrema((inodes[e[1]], inodes[e[2]])) for e in edgevs]) # List of edge (m,n) node indices with m < n + iedges = Dict(t => i for (i,t) in pairs(edgenodes)) # returns linear edge index given (m,n) node indices with m < n + e1 = first.(edgenodes) + e2 = last.(edgenodes) + fv = [inodes[f[k]] for k in 1:3, f in facevs] # face-vertex matrix + fe = zeros(Int, 3, nface) + previ = (3, 1, 2) + nexti = (2, 3, 1) + for (jf, f) in pairs(facevs), i in 1:3 + edgenodes = (inodes[f[nexti[i]]], inodes[f[previ[i]]]) + fe[i, jf] = iedges[extrema(edgenodes)] + end + + ρ = [SV2([xn[ix], yn[iy]]) for (ix, iy) in nodes] + + sh = RWGSheet() + sh.ρ = ρ + sh.e1 = e1 + sh.e2 = e2 + sh.fv = fv + sh.fe = fe + + test_fefv(sh) + + # Mirror in x: + sh2 = deepcopy(sh) + for i in eachindex(sh2.ρ) + sh2.ρ[i] = sh2.ρ[i] .* @SVector [-1, 1] + end + # Reverse order of face edges and nodes to stay CCW: + for i in axes(sh2.fv, 2) + sh2.fv[1, i] = sh.fv[3, i] + sh2.fv[3, i] = sh.fv[1, i] + end + for i in axes(sh2.fe, 2) + sh2.fe[1, i] = sh.fe[3, i] + sh2.fe[3, i] = sh.fe[1, i] + end + test_fefv(sh2) + + sh = combine(sh, sh2, 'x', 0.0) + test_fefv(sh) + + # Now mirror in y + sh2 = deepcopy(sh) + for i in eachindex(sh2.ρ) + sh2.ρ[i] = sh2.ρ[i] .* @SVector [1, -1] + end + # Reverse order of face edges and nodes to stay CCW: + for i in axes(sh2.fv, 2) + sh2.fv[1, i] = sh.fv[3, i] + sh2.fv[3, i] = sh.fv[1, i] + end + for i in axes(sh2.fe, 2) + sh2.fe[1, i] = sh.fe[3, i] + sh2.fe[3, i] = sh.fe[1, i] + end + test_fefv(sh2) + + sh = combine(sh, sh2, 'y', 0.0) + test_fefv(sh) + + # Restore original coordinates: + sh.ρ .+= Ref([xr[end], yr[end]]) + return sh + +end + diff --git a/test/FastSweep_test.jl b/test/FastSweep_test.jl index 0f660f73..ae577dca 100644 --- a/test/FastSweep_test.jl +++ b/test/FastSweep_test.jl @@ -18,7 +18,7 @@ using Test for (s4x4interp, dr) in zip(smat4x4s, direct_results) gsm = dr.gsm s4x4direct = MArray{Tuple{4,4}}([gsm[1,1] gsm[1,2]; gsm[2,1] gsm[2,2]]) - @test norm(s4x4interp - s4x4direct, Inf) ≤ err_limit + @test norm(s4x4interp - s4x4direct) ≤ err_limit end end nothing diff --git a/test/RWGData_test.jl b/test/RWGData_test.jl index d9f96aad..81be2c4f 100644 --- a/test/RWGData_test.jl +++ b/test/RWGData_test.jl @@ -17,8 +17,8 @@ using Test 2 6 10 12 3 7 1 5 4 8 2 6] - @test rwgdata.ufp2fp == [[1, 11], [2, 12], [3], [4], [5, 15], [6, 16], [7], - [8], [9], [10], [13], [14]] + @test Set.(rwgdata.ufp2fp) == Set.([[1, 11], [2, 12], [3], [4], [5, 15], [6, 16], [7], + [8], [9], [10], [13], [14]]) @test rwgdata.nufp == 12 end diff --git a/test/RWGSheet_test.jl b/test/RWGSheet_test.jl index d4cb2beb..5395a758 100644 --- a/test/RWGSheet_test.jl +++ b/test/RWGSheet_test.jl @@ -106,7 +106,7 @@ end positions = [[Float32(v[1]), Float32(v[2]), 0.0f0] for v in vec(sheet.ρ[sheet.fv])] @test mshb.position == positions @test msha.position ≈ mshb.position - @test all(==([0,0,1]), msha.normals) - @test all(==([0,0,1]), mshb.normals) + @test all(==([0,0,1]), msha.normal) + @test all(==([0,0,1]), mshb.normal) end \ No newline at end of file diff --git a/test/fresnel_test.jl b/test/fresnel_test.jl new file mode 100644 index 00000000..164c9438 --- /dev/null +++ b/test/fresnel_test.jl @@ -0,0 +1,21 @@ +using PSSFSS +using Test +@testset "TEP Creation" begin +dwidth = 3mm +duroid = Layer(epsr=2.2, tandel=0.0009, width=dwidth) +strata = [Layer(), duroid, Layer(width=-dwidth)] +steering = (θ=0:10:80, ϕ=0) +FGHz = 1:2 +rttbl = tempname() +resultfile = tempname() +logfile = devnull +showprogress = false +analyze(strata, FGHz, steering; logfile, resultfile, showprogress) +testfile_fresnel = tempname() +res2fresnel(resultfile, testfile_fresnel) +testdat = readlines(testfile_fresnel) +testdat = testdat[10:end] +gooddat = readlines(joinpath(@__DIR__, "pssfss_diel_fresnel_table.rttbl")) +gooddat = gooddat[10:end] +@test all(x -> isequal(x[1],x[2]), zip(testdat, gooddat)) +end diff --git a/test/pssfss_diel_fresnel_table.rttbl b/test/pssfss_diel_fresnel_table.rttbl new file mode 100644 index 00000000..988e44fb --- /dev/null +++ b/test/pssfss_diel_fresnel_table.rttbl @@ -0,0 +1,29 @@ +# HFSS-compatible Fresnel reflection/transmission table created by PSSFSS +# Created on 2025-01-10 at 14:10:43.349 +RTTable +# = - 1 +9 +# MultiFreq +MultiFreq 1.0 2.0 1 +# Data section follows. Frequency loops within theta +# +-0.00384 -0.03745 0.00384 0.03745 0.99852 -0.03764 0.99852 -0.03764 +-0.01504 -0.07327 0.01504 0.07327 0.99427 -0.07477 0.99427 -0.07477 +-0.00389 -0.03802 0.00369 0.03636 0.99848 -0.03822 0.99856 -0.03759 +-0.01522 -0.07441 0.01447 0.07118 0.99409 -0.07591 0.99443 -0.07472 +-0.00403 -0.03985 0.00326 0.03309 0.99833 -0.04004 0.99868 -0.03752 +-0.01579 -0.07801 0.01279 0.06486 0.99352 -0.07951 0.99490 -0.07464 +-0.00432 -0.04325 0.00259 0.02755 0.99804 -0.04344 0.99885 -0.03757 +-0.01691 -0.08467 0.01017 0.05412 0.99239 -0.08618 0.99556 -0.07486 +-0.00485 -0.04889 0.00173 0.01955 0.99750 -0.04908 0.99903 -0.03812 +-0.01898 -0.09568 0.00679 0.03849 0.99031 -0.09722 0.99623 -0.07606 +-0.00588 -0.05823 0.00072 0.00855 0.99648 -0.05843 0.99912 -0.03989 +-0.02294 -0.11378 0.00279 0.01688 0.98632 -0.11539 0.99657 -0.07969 +-0.00813 -0.07472 -0.00051 -0.00683 0.99423 -0.07495 0.99893 -0.04457 +-0.03153 -0.14524 -0.00205 -0.01351 0.97766 -0.14703 0.99583 -0.08902 +-0.01451 -0.10856 -0.00247 -0.03123 0.98782 -0.10885 0.99782 -0.05709 +-0.05539 -0.20736 -0.00979 -0.06162 0.95357 -0.20962 0.99145 -0.11366 +-0.04757 -0.20666 -0.01007 -0.08822 0.95469 -0.20717 0.99079 -0.10143 +-0.16628 -0.36070 -0.03903 -0.17068 0.84163 -0.36428 0.96432 -0.19776 +-0.04757 -0.20666 -0.01007 -0.08822 0.95469 -0.20717 0.99079 -0.10143 +-0.16628 -0.36070 -0.03903 -0.17068 0.84163 -0.36428 0.96432 -0.19776 diff --git a/test/runtests.jl b/test/runtests.jl index ec623958..5d9fb728 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -37,6 +37,13 @@ end @safetestset "TEP test" begin include("tep_test.jl") end +@safetestset "Fresnel test" begin + include("fresnel_test.jl") +end + +@safetestset "sympixels test" begin + include("sympixels_test.jl") +end if get(ENV, "BENCHMARK", "false") == "true" include("benchmark.jl") diff --git a/test/sympixels_test.jl b/test/sympixels_test.jl new file mode 100644 index 00000000..0c5e103b --- /dev/null +++ b/test/sympixels_test.jl @@ -0,0 +1,21 @@ +using PSSFSS, Test +using LinearAlgebra: I +@testset "sympixels" begin + P = 5 + units = mm + nrim = 0 + halfnint = 2 + pdiv = 2 + + veclen = halfnint * (halfnint + 1) ÷ 2 + patternvec = ones(Bool, veclen) + + sheet = sympixels(; P, nrim, halfnint, patternvec, units, pdiv, class='J') + flist = 1 + steering = (θ=0, ϕ=0) + logfile = resultfile = devnull + showprogress = false + res = analyze([Layer(), sheet, Layer()], flist, steering; showprogress, logfile, resultfile) |> only + gsm4x4 = [res.gsm[1,1] res.gsm[1,2]; res.gsm[2,1] res.gsm[2,2]] + @test maximum(abs, gsm4x4 + I) < 1e-5 +end \ No newline at end of file