PSSFSS Examples
Symmetric Strip Grating
This example consists of a symmetric strip grating, i.e. a grating where the strip width is half the unit cell period $P$:
Only three of the infinite number of strips in the grating are shown, and they extend infinitely to the left and right. The grating lies in the $z=0$ plane with free space on both sides. The shaded areas represent metalization. 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:
"""
grating(kP, nterms=30) -> (Γ, T)
Compute the normal incidence reflecton and transmission coefficients of a symmetric grid of
zero-thickness conducting strips. The product of the period of the strips and the incident
electric field wavenumber is `kP` (dimensionless). The incident electric field is
perpendicular to the direction along the axis of the strips. The series have been
accelerated by applying a Kummer's transformation, using the first two terms in the Maclaurin
series for the inverse sin function. `kP` must be in the half-open interval [0,1). The
default number of summed terms `nterms` yields better than 10 digits of accuracy over the
interval [0.01,0.99].
"""
function grating(kP; nterms=30)
sum1 = 1.3862943611198906 # \sum_{n=1}^{\infty} 1/(n-1/2) - 1/n = log(4)
sum3 = 7.2123414189575710 # \sum_{n=1}^{\infty} (n-1/2)^{-3} - n^{-3} = 6 * \zeta(3)
x = kP/(4π)
θ = x*sum1 + x^3/6 * sum3
for n = 1:nterms
xonmhalf = x/(n - 0.5)
xon = x/n
term = asin(xonmhalf) - (xonmhalf + (xonmhalf)^3/6) -
(asin(xon) - (xon + xon^3/6))
θ += term
end
Γ = sin(θ) * cis(-π/2 - θ)
T = 1 + Γ
return (Γ, T)
endMain.gratingNote that using the extension of Babinet's Principle for electromagnetic fields 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.
Here is the PSSFSS code to analyze this structure using electric currents as the unknowns. We will scale the geometry so that the frequency in GHz is numerically equal to the period of the strips measured in wavelengths.
using Plots, PSSFSS
c = 11.802852677165355 # light speed [inch*GHz]
period = c # so the period/wavelength = freq in GHz
Py = period
Ly = period/2
Px = Lx = Ly/10 # Infinite in x direction so this can be anything
Ny = 60
Nx = round(Int, Ny*Lx/Ly)
sheet = rectstrip(;Px, Py, Lx, Ly, Nx, Ny, units=inch)
flist = 0.02:0.02:0.98
steering = (θ=0, ϕ=0)
strata = [Layer()
sheet
Layer()]
results_j = analyze(strata, flist, steering, showprogress=false,
resultfile=devnull, logfile=devnull);
p1 = plot(sheet)
p2 = plot(sheet, unitcell=true)
ptitle = plot(title = "Symmetric Strip Triangulation",
grid = false, showaxis = false, xtick=[], ytick=[],
bottom_margin = -50Plots.px)
plot(ptitle, p1, p2, layout = @layout([A{0.09h}; [B C]]))
Note that setting Lx = Px causes the strip to fully occupy the x-extent of the unit cell. PSSFSS automatically ensures that the triangle edges at these unit cell boundaries define basis functions that satisfy the Floquet (phase shift) boundary conditions, so that currents are free to flow across these unit cell boundaries.
We can also analyze the same structure using magnetic currents in the areas free of metalization as the unknowns:
sheet = rectstrip(;class='M', Px, Py, Lx, Ly, Nx, Ny, units=inch)
strata = [Layer()
sheet
Layer()]
results_m = analyze(strata, flist, steering, showprogress=false,
resultfile=devnull, logfile=devnull);Under Julia 1.8, the first 50-frequency run of analyze takes about 9 seconds for this geometry of 720 triangles on my machine, and the second run takes 4 seconds. The additional time for the first run is due to JIT (just-in-time) compilation. Under Julia 1.9 the time of the first run is 3.8 seconds and the second run requires 3.2 seconds, a much smaller difference. This is due to improvements in Julia 1.9's ability to save precompiled native code for later reuse. More detailed timing information for PSSFSS runs is available in the log file (which is omitted for generating this documentation).
We will compare the PSSFSS results to the analytic solution:
# Generate exact results:
rt = grating.(2π*flist)
rperp_exact = first.(rt)
tperp_exact = last.(rt)
rpar_exact = -tperp_exact
tpar_exact = -rperp_exact;Obtain PSSFSS results for electric and magnetic currents:
outrequest = @outputs s11(v,v) s21(v,v) s11(h,h) s21(h,h)
rperp_j, tperp_j, rpar_j, tpar_j = eachcol(extract_result(results_j, outrequest))
rperp_m, tperp_m, rpar_m, tpar_m = eachcol(extract_result(results_m, outrequest));Generate the comparison plots:
angdeg(z) = rad2deg(angle(z)) # Convenience function
p1 = plot(title = "Perpendicular Reflection Magnitude",
xlabel = "Period (wavelengths)",
ylabel = "Coefficient Magnitude",
legend=:topleft)
plot!(p1, flist, abs.(rperp_exact), ls=:dash, label="Exact")
plot!(p1, flist, abs.(rperp_j), label="PSSFSS J")
plot!(p1, flist, abs.(rperp_m), label="PSSFSS M")
p2 = plot(title = "Perpendicular Reflection Phase",
xlabel = "Period (wavelengths)",
ylabel = "Phase (deg)")
plot!(p2, flist, angdeg.(rperp_exact), ls=:dash, label="Exact")
plot!(p2, flist, angdeg.(rperp_j), label="PSSFSS J")
plot!(p2, flist, angdeg.(rperp_m), label="PSSFSS M")
p1 = plot(title = "Parallel Reflection Magnitude",
xlabel = "Period (wavelengths)",
ylabel = "Coefficient Magnitude")
plot!(p1, flist, abs.(rpar_exact), ls=:dash, label="Exact")
plot!(p1, flist, abs.(rpar_j), label="PSSFSS J")
plot!(p1, flist, abs.(rpar_m), label="PSSFSS M")
p2 = plot(title = "Parallel Reflection Phase",
xlabel = "Period (wavelengths)",
ylabel = "Phase (deg)")
plot!(p2, flist, angdeg.(rpar_exact), ls=:dash, label="Exact")
plot!(p2, flist, angdeg.(rpar_j), label="PSSFSS J")
plot!(p2, flist, angdeg.(rpar_m), label="PSSFSS M")
Now look at the transmission coefficients:
p1 = plot(title = "Perpendicular Transmission Magnitude",
xlabel = "Period (wavelengths)",
ylabel = "Coefficient Magnitude")
plot!(p1, flist, abs.(tperp_exact), ls=:dash, label="Exact")
plot!(p1, flist, abs.(tperp_j), label="PSSFSS J")
plot!(p1, flist, abs.(tperp_m), label="PSSFSS M")
p2 = plot(title = "Perpendicular Transmission Phase",
xlabel = "Period (wavelengths)",
ylabel = "Phase (deg)")
plot!(p2, flist, angdeg.(tperp_exact), ls=:dash, label="Exact")
plot!(p2, flist, angdeg.(tperp_j), label="PSSFSS J")
plot!(p2, flist, angdeg.(tperp_m), label="PSSFSS M")
p1 = plot(title = "Parallel Transmission Magnitude",
xlabel = "Period (wavelengths)",
ylabel = "Coefficient Magnitude", legend=:topleft)
plot!(p1, flist, abs.(tpar_exact), ls=:dash, label="Exact")
plot!(p1, flist, abs.(tpar_j), label="PSSFSS J")
plot!(p1, flist, abs.(tpar_m), label="PSSFSS M")
p2 = plot(title = "Parallel Transmission Phase",
xlabel = "Period (wavelengths)",
ylabel = "Phase (deg)")
plot!(p2, flist, angdeg.(tpar_exact), ls=:dash, label="Exact")
plot!(p2, flist, angdeg.(tpar_j), label="PSSFSS J")
plot!(p2, flist, angdeg.(tpar_m), label="PSSFSS M")
Conclusion
Although very good agreement is obtained, as expected the best agreement between all three results occurs for the lowest frequencies, where the triangles are smallest in terms of wavelength. This emphasizes the fact that it is necessary for the user to check that enough triangles have been requested for good convergence over the frequency band of interest. This example is an extremely demanding case in terms of bandwidth, as the ratio of maximum to minimum frequency here is $0.98/0.02 = 49:1$
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 is from Section 3.2, Figures 7 and 8. We will compare PSSFSS results to those digitized from the cited figures.
We start by defining a function that creates a patch of the desired sheet resistance:
using Plots, PSSFSS
patch(Z) = rectstrip(Nx=10, Ny=10, Px=1, Py=1, Lx=0.5, Ly=0.5, units=cm, Zsheet=Z)
plot(patch(0), unitcell=true)
The patches measure 0.5 cm on a side and lie in a square lattice of period 1 cm. Now we perform the analysis, looping over the desired values of sheet resistance.
steering = (ϕ=0, θ=0)
flist = 1:0.5:60
Rs = [0, 10, 30, 100]
calculated = zeros(length(flist), length(Rs)) # preallocate storage
outputs = @outputs s11mag(v,v)
for (i,R) in pairs(Rs)
strata = [Layer(), patch(R), Layer()]
results = analyze(strata, flist, steering, showprogress=false,
logfile=devnull, resultfile=devnull)
calculated[:,i] = extract_result(results, outputs)
endLooping over the four sheet resistance values, each evaluated at 119 frequencies required approximately 20 seconds total on my machine.
We plot the results, including those digitized from the paper for comparison:
using DelimitedFiles
markers = (:diamond, :utriangle, :square, :xcross)
colors = (:blue, :red, :green, :black)
p = plot(xlim=(-0.01,60.01), xtick = 0:10:60, ylim=(-0.01,1.01), ytick=0:0.1:1,
xlabel="Frequency (GHz)", ylabel="Reflection Coefficient Magnitude",
title = "Resistive Square Patch",
legend=:topright)
for (i,R) in pairs(Rs)
scatter!(p, flist, calculated[:,i], label="PSSFSS $R Ω", ms=2, shape=markers[i], color=colors[i])
data = readdlm("../src/assets/barlevy_patch_digitized_$(R)ohm.csv", ',')
plot!(p, data[:,1], data[:,2], label="Barlevy $R Ω", ls=:dash, color=colors[i])
end
p
Conclusion
PSSFSS results are indistinguishable from those reported in the cited paper.
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.
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:
using Plots, PSSFSS, DelimitedFiles
sheet = loadedcross(w=1.0, L1=0.6875, L2=0.0625, s1=[1.0,0.0],
s2=[0.0,1.0], ntri=600, units=cm)
plot(sheet, unitcell=true, linecolor=:red)
A few things to note. First, as of PSSFSS version 1.3, the mesh is structured. So there are redundant triangle face-pairs that PSSFSS can exploit to reduce execution time. Second, the number of triangle faces generated is only approximately equal to the requested value of 600. This can be verified by entering the Julia variable sheet at the REPL (i.e. the Julia prompt):
sheetRWGSheet: style=loadedcross, class=J, 405 nodes, 1044 edges, 640 faces, Zs=0.0 + 0.0im ΩAlternatively, the facecount function will return the number of triangle faces on the sheet:
facecount(sheet)640The cross FSS is etched on a dielectric sheet of thickness 3 mm. The dielectric constant is varied over the values 1, 2, and 4 to observe the effect on the resonant frequency. Following the reference, the list of analysis frequencies is varied slightly depending on the value of dielectric constant:
resultsstack = Any[]
steering = (ϕ=0, θ=0)
for eps in [1, 2, 4]
strata = [ Layer()
sheet
Layer(ϵᵣ=eps, width=3mm)
Layer()
]
if eps == 1
flist = 1:0.2:30
elseif eps == 2
flist = 1:0.2:26
else
flist = 1:0.2:20
end
results = analyze(strata, flist, steering, showprogress=false,
resultfile=devnull, logfile=devnull)
push!(resultsstack, results)
endThe above loop requires about 18 seconds of execution time on my machine. Compare PSSFSS results to those digitized from the dissertation figure:
col=[:red,:blue,:green]
p = plot(xlim=(0.,30), xtick = 0:5:30, ylim=(0,1), ytick=0:0.1:1,
xlabel="Frequency (GHz)", ylabel="Reflection Coefficient Magnitude",
legend=:topleft, lw=2)
for (i,eps) in enumerate([1,2,4])
data = extract_result(resultsstack[i], @outputs FGHz s11mag(v,v))
plot!(p, data[:,1], data[:,2], label="PSSFSS ϵᵣ = $eps", lc=col[i])
data = readdlm("../src/assets/barlevy_diss_eps$(eps).csv", ',')
plot!(p, data[:,1], data[:,2], label="Barlevy ϵᵣ = $eps", lc=col[i], ls=:dot)
end
p
Conclusion
PSSFSS results agree very well with those of the cited reference, especially when accounting for the fact that the reference results are obtained by digitizing a scanned figure.
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 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.
using Plots, PSSFSS, DelimitedFiles
D = 11 # Period of square lattice (mm)
r_outer = √2/2 * D/8 * [5,6,7] # radii of square outer vertices
thickness = D/16 * [1,2,3]
r_inner = r_outer - √2 * thickness # radii of square inner vertices
Rs = [15,40,70] # Sheet resistance (Ω/□)
labels = ["Thin", "Medium", "Thick"]
colors = [:green, :blue, :red]
p = plot(title="Costa Absorber", xlim=(0,25),ylim=(-35,0),xtick=0:5:25,ytick=-35:5:0,
xlabel="Frequency (GHz)", ylabel="Reflection Magnitude (dB)", legend=:bottomleft)
ps = []
for (ri, ro, label, color, R) in zip(r_inner, r_outer, labels, colors, Rs)
sheet = polyring(sides=4, s1=[D, 0], s2=[0, D], ntri=750, orient=45,
a=[ri], b=[ro], Zsheet=R, units=mm)
push!(ps, plot(sheet, unitcell=true, title=label, lc=color))
strata = [Layer()
sheet
Layer(width=5mm)
pecsheet() # Perfectly conducting ground plane
Layer()]
results = analyze(strata, 1:0.2:25, (ϕ=0, θ=0), showprogress=false,
resultfile=devnull, logfile=devnull)
data = extract_result(results, @outputs FGHz s11dB(h,h))
plot!(p, data[:,1], data[:,2], label="PSSFSS "*label, lc=color)
dat = readdlm("../src/assets/costa_2014_" * lowercase(label) * "_reflection.csv", ',')
plot!(p, dat[:,1], dat[:,2], label="Costa "*label, ls=:dash, lc=color)
end
plot(ps..., layout=(1,3), size=(600,220), margin=3Plots.mm)
This PSSFSS run of three geometries takes about 15 seconds on my machine.
p
It is useful to take a look at the log file created by PSSFSS for the last run above (from a previous run where the log file was not discarded):
Starting PSSFSS 1.2.2 analysis on 2023-03-04 at 19:44:10.800
Julia Version 1.8.5
Commit 17cfb8e65e (2023-01-08 06:45 UTC)
Platform Info:
OS: Windows (x86_64-w64-mingw32)
CPU: 8 × Intel(R) Core(TM) i7-9700 CPU @ 3.00GHz
WORD_SIZE: 64
LIBM: libopenlibm
LLVM: libLLVM-13.0.1 (ORCJIT, skylake)
Threads: 8 on 8 virtual cores
BLAS: LBTConfig([ILP64] libopenblas64_.dll)
Dielectric layer information...
Layer Width units epsr tandel mur mtandel modes beta1x beta1y beta2x beta2y
----- ------------- ------- ------ ------- ------ ----- ------- ------- ------- -------
1 0.0000 mm 1.00 0.0000 1.00 0.0000 2 571.2 -0.0 -0.0 571.2
================== Sheet 1 ======================== 571.2 -0.0 -0.0 571.2
2 5.0000 mm 1.00 0.0000 1.00 0.0000 42 571.2 -0.0 -0.0 571.2
================== Sheet 2 ======================== 0.0 0.0 0.0 0.0
3 0.0000 mm 1.00 0.0000 1.00 0.0000 2 571.2 -0.0 -0.0 571.2
PSS/FSS sheet information...
Sheet Loc Style Rot J/M Faces Edges Nodes Unknowns NUFP
----- --- ---------------- ----- --- ----- ----- ----- -------- ------
1 1 polyring 0.0 J 720 1152 432 1008 199676
2 2 NULL 0.0 E 0 0 0 0 0
⋮Note from the dielectric layer report that there are 42 modes defined in the region between the ground plane and the FSS sheet. This is the number of modes selected by the code to include in the generalized scattering matrix formulation to properly account for electromagnetic coupling between the two surfaces. If the 5 mm spacing were increased to, say, 7 mm then fewer modes would be needed. Also note in the FSS sheet information that NUFP (the number of unique face pairs) 199676, is less than the number of faces squared ($567009 = 753^2$), a consequence of the structured triangulation used for a 4-sided polyring.
Conclusion
PSSFSS results agree very well with those of the paper, except for the medium width loop, where the agreement is not quite as good. It was found empirically that using a slightly different value of Rs = 37 for this ring results in nearly perfect agreement with the digitized results.
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. 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 from the Figure 2(a) of the paper. The latter was obtained by the authors using CST Microwave Studio.
using PSSFSS, Plots, DelimitedFiles
p = 20 # Period of square unit cell (mm)
a, b, Rsquare = (0.5, 19, 120) # Parameters of squares
c, d, Rcircle = (1, 9, 400) # Parameters of circles
units = mm; Px = Py = p; Lx = Ly = p - 2a; Nx = Ny = 21
squares = rectstrip(; Px, Py, Lx, Ly, Nx, Ny, units, Zsheet=Rsquare)
s1 = [p,0]; s2 = [0, p]; ntri = 800; sides = 40
disks = polyring(; units, sides, a=[0], b=[d], s1, s2, ntri, Zsheet=Rcircle)
psq = plot(squares, unitcell=true, linecolor=:red)
pdi = plot(disks, unitcell=true, linecolor=:red)
psheet = plot(psq, pdi)
plot!(psheet, title=["Resistive Squares" "Resistive Disks"], size=(600,300))
silicone = Layer(width = 4mm, ϵᵣ = 2.9, tanδ = 0.1) # using unicode keywords
foam = Layer(width = 6mm, epsr = 1.05, tandel = 0.02) # using ASCII keywords
strata = [Layer(), silicone, disks, silicone, squares, foam, pecsheet(), Layer()]
FGHz = 0.5:0.1:20
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')
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)
This PSSFSS run takes about 44 seconds on my machine for 196 frequencies covering a 40:1 bandwidth.
Conclusion
PSSFSS results agree well with those of the paper.
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.
It consists of three concentric split rings with gaps sequentially rotated by 180° situated on a thin dielectric slab.
Here is a script that analyzes this geometry:
using PSSFSS, Plots, DelimitedFiles
r123 = [3.7, 4.25, 4.8]
d = 10.8
w = 0.3
g = 0.3
a = r123 .- w/2
b = a .+ w
gapwidth = [g, g, g]
gapcenter = [-90, 90, -90]
s1 = [d, 0]
s2 = [0, d]
sheet = splitring(;class='M', units=mm, sides=42, ntri=1500, a, b, s1, s2, gapwidth, gapcenter)
display(plot(sheet, unitcell=true))
strata = [
Layer()
sheet
Layer(ϵᵣ=10.2, tanδ=0.0023, width=0.13mm)
Layer()
]
freqs = 1:0.02:14)
steering = (θ = 0, ϕ = 0)
results = analyze(strata, freqs, steering)
s11dbvv = extract_result(results, @outputs s11db(v,v))
p = plot(xlabel="Frequency (GHz)", ylabel="S₁₁ Amplitude (dB)",
xtick=1:14, xlim=(1, 14), ytick = -40:5:0, ylim=(-30,0),
legend=:bottom)
plot!(p, freqs, s11dbvv, label="PSSFSS")
dat = readdlm("../src/assets/fabian2021_fig3_digitized.csv", ',')
plot!(p, dat[:,1], dat[:,2], label="Fabian (CST)")
This run of 651 frequencies requires about 40 seconds on my machine. Generally good agreement is seen between the PSSFSS predicted reflection amplitude and that digitized from the paper. (The latter was obtained from a CST frequency domain analysis, according to the paper's authors.) However, there is a small discrepancy in the predicted resonant frequencies that increases with frequency, likely because both results are less well converged at higher frequencies. Also, the reflection amplitudes of the higher-frequency peaks are less than unity for the CST results, possibly because the authors may have included the finite conductivity of the metal traces. This detail was not reported in the paper.
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.
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 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 desired number of triangles ntri:
function element(L2, ntri::Int)
L = 17 # Period of square lattice (mm)
k1 = 0.5
k2 = 0.125
L1 = k1 * L2
w = k2 * L2
a = [-1.0, (L2-2w)/√2]
b = [L1/√2, L2/√2]
return polyring(; sides=4, s1=[L, 0], s2=[0, L], ntri, orient=45, a, b, units=mm)
endelement (generic function with 1 method)Here is a plot of the two elements at the size extremes to be examined:
using PSSFSS, Plots
p1 = plot(element(2, 675), title="L2 = 2mm", unitcell=true, linecolor=:red)
p2 = plot(element(16.5, 3416), title="L2 = 16.5mm", unitcell=true, linecolor=:red)
p = plot(p1,p2)
Here is the script that generates the S-curve. To ensure accurate phases, a convergence check is performed for each distinct L2 value. The number of triangles is increased by a factor of 1.5 if the previous increase resulted in a reflection phase change of more than one degree. After computing the phases, comparison is made to the same data computed by Li, et al from Ansoft HFSS and CST Microwave Studio, and displayed in their Figure 6.
Note that the width for the first Layer is set to -7mm. This has the effect of referring reflection phases to the location of the ground plane, as was apparently done by Li, et al. Also note that the unwrap! function of the DSP package is used to "unwrap" phases to remove any possible discontinuous jumps of 360 degrees.
using PSSFSS
using Plots, DelimitedFiles
using DSP: unwrap!
p = plot(title="Reflection Phase Convergence", xlim=(2,18),ylim=(-350,150),xtick=2:2:18,ytick=-350:50:150,
xlabel="Large Square Dimension L2 (mm)", ylabel="Reflecton Phase (deg)", legend=:topright)
FGHz = 10.0
L2s = range(2.0, 16.5, length=30)
nfactor = 1.5 # Growth factor for ntri
phasetol = 1.0 # Phase convergence tolerance
record = Tuple{Float64, Int, Int}[] # Storage for (L2,ntri,facecount(sheet))
phases = zeros(length(L2s))
ntri = 300
for (il2, L2) in pairs(L2s)
print("L2 = $L2; ntri =")
if il2 > 1
global ntri = round(Int, ntri / (nfactor*nfactor))
end
oldphase = Inf
firstrun = true
while true
ntri = round(Int, nfactor * ntri)
if firstrun
print(" $ntri")
firstrun = false
else
print(", $ntri")
end
sheet = element(L2, ntri)
strata = [Layer(width=-7mm)
sheet
Layer(ϵᵣ=2.65, width=1mm)
Layer(ϵᵣ=1.07, width=6mm)
pecsheet() # ground plane
Layer()]
results = analyze(strata, FGHz, (ϕ=0, θ=0), showprogress=false, logfile="li2009.log", resultfile="li2009.res")
phase = extract_result(results, @outputs s11ang(h,h))[1,1]
phase < 0 && (phase += 360)
Δphase = abs(phase - oldphase)
if Δphase > phasetol
oldphase = phase
else
phases[il2] = phase
push!(record, (L2, ntri, facecount(sheet)))
println("; Δphase = $(round(Δphase, digits=2))°")
break
end
end
end
unwrap!(phases, range=360)
p = plot(title="Li 2009 Reflectarray", xlim=(2,18),ylim=(-350,150),xtick=2:2:18,ytick=-350:50:150,
xminorticks=2, yminorticks=2,
xlabel="Large Square Dimension L2 (mm)", ylabel="Reflecton Phase (deg)", legend=:topright)
plot!(p, L2s, phases, label="PSSFSS", color=:red, mscolor=:red, ms=3, shape=:circ)
dat = readdlm("li2009_ansoft_digitized.csv", ',')
scatter!(p, dat[:,1], dat[:,2], label="Li Ansoft", mc=:blue, msc=:blue, markershape=:square)
dat = readdlm("li2009_cst_digitized.csv", ',')
scatter!(p, dat[:,1], dat[:,2], label="Li CST", mc=:green, msc=:green, markershape=:circle)
display(p)
println()
record
It can be seen that the PSSFSS phases compare well with the HFSS phases, better than the CST and HFSS phases compare to each other. The authors of the paper do not discuss checking the convergence of their results or even any details of how they set up their HFSS and CST models.
The console output from the above script is shown below:
julia> include("li2009_convergence.jl")
L2 = 2.0; ntri = 450, 675; Δphase = 0.02°
L2 = 2.5; ntri = 450, 675; Δphase = 0.04°
L2 = 3.0; ntri = 450, 675; Δphase = 0.06°
L2 = 3.5; ntri = 450, 675; Δphase = 0.1°
L2 = 4.0; ntri = 450, 675; Δphase = 0.14°
L2 = 4.5; ntri = 450, 675; Δphase = 0.18°
L2 = 5.0; ntri = 450, 675; Δphase = 0.19°
L2 = 5.5; ntri = 450, 675; Δphase = 0.15°
L2 = 6.0; ntri = 450, 675; Δphase = 0.02°
L2 = 6.5; ntri = 450, 675; Δphase = 0.41°
L2 = 7.0; ntri = 450, 675, 1012; Δphase = 0.12°
L2 = 7.5; ntri = 675, 1012; Δphase = 0.23°
L2 = 8.0; ntri = 675, 1012; Δphase = 0.38°
L2 = 8.5; ntri = 675, 1012; Δphase = 0.55°
L2 = 9.0; ntri = 675, 1012; Δphase = 0.73°
L2 = 9.5; ntri = 675, 1012; Δphase = 0.93°
L2 = 10.0; ntri = 675, 1012, 1518, 2277, 3416; Δphase = 0.67°
L2 = 10.5; ntri = 2277, 3416; Δphase = 0.54°
L2 = 11.0; ntri = 2277, 3416; Δphase = 0.35°
L2 = 11.5; ntri = 2277, 3416; Δphase = 0.1°
L2 = 12.0; ntri = 2277, 3416; Δphase = 0.17°
L2 = 12.5; ntri = 2277, 3416; Δphase = 0.41°
L2 = 13.0; ntri = 2277, 3416; Δphase = 0.59°
L2 = 13.5; ntri = 2277, 3416; Δphase = 0.69°
L2 = 14.0; ntri = 2277, 3416; Δphase = 0.73°
L2 = 14.5; ntri = 2277, 3416; Δphase = 0.72°
L2 = 15.0; ntri = 2277, 3416; Δphase = 0.69°
L2 = 15.5; ntri = 2277, 3416; Δphase = 0.65°
L2 = 16.0; ntri = 2277, 3416; Δphase = 0.61°
L2 = 16.5; ntri = 2277, 3416; Δphase = 0.59°
30-element Vector{Tuple{Float64, Int64, Int64}}:
(2.0, 675, 722)
(2.5, 675, 722)
(3.0, 675, 722)
(3.5, 675, 722)
(4.0, 675, 722)
(4.5, 675, 722)
(5.0, 675, 722)
(5.5, 675, 722)
(6.0, 675, 722)
(6.5, 675, 722)
(7.0, 1012, 944)
(7.5, 1012, 944)
(8.0, 1012, 944)
(8.5, 1012, 944)
(9.0, 1012, 944)
(9.5, 1012, 944)
(10.0, 3416, 3314)
(10.5, 3416, 3314)
(11.0, 3416, 3314)
(11.5, 3416, 3314)
(12.0, 3416, 3314)
(12.5, 3416, 3314)
(13.0, 3416, 3314)
(13.5, 3416, 3314)
(14.0, 3416, 3314)
(14.5, 3416, 3314)
(15.0, 3416, 3314)
(15.5, 3416, 3314)
(16.0, 3416, 3314)
(16.5, 3416, 3314)The first list shows how the script increases the number of triangles requested for each geometry until the reflection phase is sufficiently converged. The printout of the record array shows both the final requested value of ntri and actual number of triangle faces generated by the mesher for each L2 value. The final cases with ntri=3416 required about 21 seconds each of execution time.
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 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 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.
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. The closely spaced sheets are a good test of the generalized scattering formulation implemented in PSSFSS. The sheet geometry is shown below. Remember that the entire sheet is metalized except for the region of the triangulation.
using Plots, PSSFSS
sheet = loadedcross(class='M', w=0.023, L1=0.8, L2=0.14,
s1=[0.861,0.0], s2=[0.0,0.861], ntri=2400, units=cm)
plot(sheet, linecolor=:red, unitcell=true)
steering = (ϕ=0, θ=0)
strata = [ Layer()
Layer(ϵᵣ=1.3, width=1.1cm)
sheet
Layer(ϵᵣ=1.9, width=0.6cm)
sheet
Layer(ϵᵣ=1.3, width=1.1cm)
Layer() ]
flist = 1:0.1:20
results = analyze(strata, flist, steering, resultfile=devnull,
logfile=devnull, showprogress=false)
data = extract_result(results, @outputs FGHz s21db(v,v) s11db(v,v))
using DelimitedFiles
dat = readdlm("../src/assets/hfss_loadedcross_bpf_full.csv", ',', skipstart=1)
p = plot(xlabel="Frequency (GHz)", ylabel="Reflection Coefficient (dB)",
legend=:left, title="Loaded Cross Band-Pass Filter", xtick=0:2:20, ytick=-30:5:0,
xlim=(-0.1,20.1), ylim=(-35,0.1))
plot!(p, data[:,1], data[:,3], label="PSSFSS", color=:red)
plot!(p, dat[:,1], dat[:,2], label="HFSS", color=:blue)
p2 = plot(xlabel="Frequency (GHz)", ylabel="Transmission Coefficient (dB)",
legend=:bottom, title="Loaded Cross Band-Pass Filter", xtick=0:2:20, ytick=-80:10:0,
xlim=(-0.1,20.1), ylim=(-80,0.1))
plot!(p2, data[:,1], data[:,2], label="PSSFSS", color=:red)
plot!(p2, dat[:,1], dat[:,4], label="HFSS", color=:blue)
This analysis takes about 173 seconds for 191 frequencies on my machine. The number of triangles requested is a compromise between accuracy and speed. Even better agreement with the HFSS result would be possible with more triangles, at the cost of execution time. Note that rather than including two separate invocations of the loadedcross function when defining the strata, I referenced the same sheet object in the two different locations. This allows PSSFSS to recognize that the triangulations are identical, and to exploit this fact in making the analysis more efficient. In fact, if both sheets had been embedded in similar dielectric claddings (in the same order), then the GSM (generalized scattering matrix) computed for the sheet in its first location could be reused without additional computation for its second location. In this case, though, only the spatial integrals are re-used. For an oblique incidence case, computing the spatial integrals is often the most expensive part of the analysis, so the savings from reusing the same sheet definition can be substantial.
Conclusion
Very good agreement is obtained versus HFSS over a large dynamic range of almost 80 dB, and over a 20:1 frequency bandwidth.
Meanderline-Based CPSS
A "CPSS" is a circular polarization selective structure, i.e., a structure that reacts differently 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.
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.
Here is the script that analyzes their design:
using PSSFSS
# Define convenience functions for sheets:
outer(rot) = meander(a=3.97, b=3.97, w1=0.13, w2=0.13, h=2.53+0.13, units=mm, ntri=600, rot=rot)
inner(rot) = meander(a=3.97*√2, b=3.97/√2, w1=0.1, w2=0.1, h=0.14+0.1, units=mm, ntri=600, rot=rot)
center(rot) = meander(a=3.97, b=3.97, w1=0.34, w2=0.34, h=2.51+0.34, units=mm, ntri=600, rot=rot)
# Note our definition of `h` differs from that of the reference by the width of the strip.
t1 = 4mm # Outer layers thickness
t2 = 2.45mm # Inner layers thickness
substrate = Layer(width=0.1mm, epsr=2.6)
foam(w) = Layer(width=w, epsr=1.05) # Foam layer convenience function
rot0 = 0 # rotation of first sheet
strata = [
Layer()
outer(rot0)
substrate
foam(t1)
inner(rot0 - 45)
substrate
foam(t2)
center(rot0 - 2*45)
substrate
foam(t2)
inner(rot0 - 3*45)
substrate
foam(t1)
outer(rot0 - 4*45)
substrate
Layer() ]
steering = (θ=0, ϕ=0)
flist = 10:0.1:20
#
results = analyze(strata, flist, steering, showprogress=false,
resultfile=devnull, logfile=devnull);Here are plots of the five meanderline sheets:
using Plots
plot(outer(rot0), unitcell=true, title="Sheet1")plot(inner(rot0-45), unitcell=true, title="Sheet2")plot(center(rot0-2*45), unitcell=true, title="Sheet3 (Center)")plot(inner(rot0-3*45), unitcell=true, title="Sheet4")plot(outer(rot0-4*45), unitcell=true, title="Sheet5")