module PSSFSS

if isdefined(Base, :Experimental) && isdefined(Base.Experimental, Symbol("@optlevel"))

@eval Base.Experimental.@optlevel 3

using Reexport

using PkgVersion

using InteractiveUtils: versioninfo

using Dates: now

using DelimitedFiles: writedlm

using Printf: @sprintf

using LinearAlgebra: ×, norm, ⋅, factorize, lu!, ldiv!, BLAS

using StaticArrays: StaticArrays, SVector, SArray, @SVector, MArray

using Unitful: ustrip, @u_str

using Logging: with_logger

using ProgressMeter

using PrecompileTools

include("Constants.jl")

include("ZhatCross.jl")

include("Log.jl")

include("PSSFSSLen.jl")

include("Rings.jl")

include("Layers.jl")

include("Sheets.jl")

include("Meshsub.jl")

include("Elements.jl")

include("RWG.jl")

include("PGF.jl")

include("Zint.jl")

include("FillZY.jl")

include("GSMs.jl")

include("Modes.jl")

include("Outputs.jl")

include("FastSweep.jl")

using .ZhatCross: ẑ

using .Rings

@reexport using .Sheets: Sheet, RWGSheet, read_sheet_data, nodecount, facecount, edgecount,

export_sheet, STL_ASCII, STL_BINARY

using .RWG: setup_rwg, rwgbfft!, RWGData

using .GSMs: GSM, cascade, cascade!, gsm_electric_gblock, gsm_magnetic_gblock,

gsm_slab_interface, translate_gsm!, choose_gblocks, Gblock, pecgsm, pmcgsm

using .FillZY: fillz, filly

using .Modes: choose_layer_modes!, setup_modes!

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

using .Outputs: Result, append_result_data

using .FastSweep: interpolate_band

export analyze

Base.isfile(f::Base.DevNull) = false

Base.open(f::Base.DevNull, ::AbstractString) = f

function analyze(strata::Vector, flist, steering; outlist=[], logfile="pssfss.log",

resultfile="pssfss.res", showprogress::Bool=true, fastsweep::Bool=true)

tstart = time()

layers = Layer[deepcopy(s) for s in strata if s isa Layer]

sheets = RWGSheet[s for s in strata if s isa Sheet]

islayer = map(x -> x isa Layer, strata)

issheet = map(x -> x isa RWGSheet, strata)

nl = length(layers)

nj = nl - 1

ns = length(sheets)

sint = cumsum(islayer)[issheet] # sint[k] contains dielectric interface number of k'th sheet

junc = zeros(Int, nj)

junc[sint] = 1:ns # junc[i] is the sheet number present at interface i, or 0 if no sheet is there

freqstemp = float.(collect(flist))

if length(freqstemp) < 2

freqs = Float64[freqstemp]

else

freqs::Vector{Float64} = freqstemp

end

stkeys::Tuple{Symbol,Symbol} = keys(steering)

stvaluestemp = [float.(collect(s)) for s in steering]

stvalues = Vector{Float64}[]

for (i, s) in pairs(stvaluestemp)

if length(s) < 2

push!(stvalues, Float64[s])

else

push!(stvalues, s)

end

end

with_logger(pssfss_logger(logfile)) do

_analyze(layers, sheets, junc, freqs, stkeys, stvalues;

outlist, resultfile, showprogress, tstart, fastsweep)

end

end # function

function _analyze(layers, sheets, junc, freqs, stkeys, stvalues;

outlist=[], resultfile="pssfss.res", showprogress::Bool=true, tstart=time(), fastsweep=true)

showprogress && println("Beginning PSSFSS Analysis")

ncount = 0 # Number of analyses performed

ntotal = length(freqs) * length(stvalues[1]) * length(stvalues[2])

showprogress && (progress = Progress(ntotal; dt = 1))

showprogress && update!(progress, 0)

isfile(resultfile) && rm(resultfile)

date, clock = split(string(now()), 'T')

pssfssv = PkgVersion.Version(PSSFSS)

ss = "Environment"

io = IOBuffer()

versioninfo(io)

juliainfo = String(take!(io)) * ss

i = findfirst(ss, juliainfo)

juliainfo = juliainfo[1:first(i)-1]

juliainfo *= " BLAS: $(BLAS.get_config())\n"

if VERSION < v"1.8"

juliainfo = juliainfo * " Threads.nthreads() = $(Threads.nthreads())\n"

end

@logfile "\n\nStarting PSSFSS $(pssfssv) analysis on $(date) at $(clock)\n$(juliainfo)\n\n"

check_inputs(layers, sheets, junc, freqs, stkeys, stvalues, outlist)

k0min, k0max = twopi * 1e9 / c₀ .* extrema(freqs)

gbls = choose_gblocks(layers, sheets, junc, k0min)

gsm_save = Vector{GSM}(undef, length(gbls)) # Storage for reusable GSMs

choose_layer_modes!(layers, sheets, junc, gbls, k0max, dbmin)

gbldup = get_gbldup(gbls, layers, sheets, junc)

usi = unique_indices(sheets)

rwgdat = RWGData[setup_rwg(sheet) for sheet in sheets]

report_layers_sheets(layers, sheets, junc, rwgdat, usi)

uvec::Vector{Float64} = map(sheets) do sh # Green's function smoothing factors

sh.style == "NULL" && (return 0.0)

ufactor = 0.5 * ustrip(Float64, sh.units, 1u"m")

ufactor * max(norm(sh.β₁), norm(sh.β₂))

end

# Begin analysis loops over steering angles and frequency

firstoutput = true

results = Result[]

for stout in stvalues[1], stin in stvalues[2]

steer = getsttuple(stkeys, stout, stin)

if keys(steer)[1] == :ψ₁

ψ₁, ψ₂ = deg2rad.(steer) # radians

upm::Float64 = ustrip(Float64, sheets[1].units, 1u"m")

β₁, β₂ = sheets[1].β₁ * upm, sheets[1].β₂ * upm

β⃗₀₀k1 = (ψ₁ * β₁ + ψ₂ * β₂) / twopi # Eq. (2.13b)

else

θ, ϕ = steer # degrees

st = sind(θ)

sp, cp = sincosd(ϕ)

β⃗₀₀k1 = @SVector([st * cp, st * sp])

end

@logfile "Beginning $(steer)"

if fastsweep

smat4x4s = interpolate_band(freqs; showprogress, xlabel="GHz") do fghz

(_, result) = compute_next_freq(fghz, β⃗₀₀k1, steer, layers, sheets, usi,

rwgdat, uvec, junc, gbls, gbldup, gsm_save)

smat = MArray{Tuple{4,4}}([result.gsm[1,1] result.gsm[1,2]; result.gsm[2,1] result.gsm[2,2]])

return smat

end

for (fghz, s4x4) in zip(freqs, smat4x4s)

k0 = pi * fghz * 2e9 / c₀

if keys(steer)[1] == :θ

k1 = k0 * sqrt(real(layers[1].ϵᵣ * layers[1].μᵣ))

β⃗₀₀ = k1 * β⃗₀₀k1

else

β⃗₀₀ = copy(β⃗₀₀k1)

end

gsm = GSM(s4x4[1:2,1:2], s4x4[1:2,3:4], s4x4[3:4,1:2], s4x4[3:4,3:4])

result = Result(gsm, steer, β⃗₀₀, fghz, layers[1].ϵᵣ, layers[1].μᵣ,

layers[1].β₁, layers[1].β₂, layers[end].ϵᵣ, layers[end].μᵣ,

layers[end].β₁, layers[end].β₂)

push!(results, result)

end

end

for fghz in freqs

ncount += 1

if fastsweep

result = results[ncount]

else

(β⃗₀₀, result) = compute_next_freq(fghz, β⃗₀₀k1, steer, layers, sheets, usi,

rwgdat, uvec, junc, gbls, gbldup, gsm_save)

push!(results, result)

end

# Write to output files

append_result_data(resultfile, string(ncount), result)

for row in eachrow(outlist)

if firstoutput

open(row[1], "w") do io

writedlm(io, permutedims([r.label for r in row[2]]), ',')

end

end

open(row[1], "a") do io

writedlm(io, permutedims([r(result) for r in row[2]]), ',')

end

end

firstoutput = false

showprogress && next!(progress) # Bump progress meter

end # Frequency loop

end # steering angle loop

date, clock = split(string(now()), 'T')

telapsed = round(time() - tstart, digits=1)

@logfile "\n\n PSSFSS analysis exiting on $(date) at $(clock) ($(telapsed) seconds elapsed time)\n\n"

showprogress && println("")

return results

end # function

function compute_next_freq(fghz, β⃗₀₀k1, steer, layers, sheets, usi, rwgdat, uvec, junc, gbls, gbldup, gsm_save)

@logfile " $(fghz) GHz"

t_freq = time()

k0 = pi * fghz * 2e9 / c₀

if keys(steer)[1] == :θ

k1 = k0 * sqrt(real(layers[1].ϵᵣ * layers[1].μᵣ))

β⃗₀₀ = k1 * β⃗₀₀k1

else

β⃗₀₀ = copy(β⃗₀₀k1)

end

t_cascade = 0.0

setup_modes!.(layers, k0, Ref(β⃗₀₀))

if !(angle(layers[begin].γ[1]) ≈ angle(layers[end].γ[1]) ≈ π / 2)

@logfile " Skipping $(fghz) GHz due to cutoff principal modes in ambient medium"

return nothing

end

# Initialize overall GSM and propagate it through layer 1's width:

n1 = length(layers[1].P)

gsma = GSM(n1, n1)

gsmc = GSM(n1, n1)

t_cascade1 = time()

cascade!(gsma, layers[1])

t_cascade += time() - t_cascade1

for (ig, gbl) in pairs(gbls) # Walk through the Gblocks

i1 = first(gbl.rng) # Index of layer to left of Gblock

i2 = 1 + last(gbl.rng) # Index of layer to right of Gblock

i_junc = gbl.j # junction where FSS is located, or 0 if no sheet

i_sheet = i_junc == 0 ? 0 : junc[i_junc]

if i_sheet ≠ 0

if gbldup[ig] > 0

gsmb::GSM = deepcopy(gsm_save[gbldup[ig]]) # Use previously calculated GSM

else

region = @view layers[i1:i2]

sheet = sheets[i_sheet]

s = gbl.j - i1 + 1 # sheet interface location within `region`

if sheet.class == 'J'

gsmb = calculate_jtype_gsm(region, sheet, uvec[i_sheet],

rwgdat[i_sheet], s, k0, β⃗₀₀, usi[i_sheet])

elseif sheet.class == 'M'

gsmb = calculate_mtype_gsm(region, sheet, uvec[i_sheet],

rwgdat[i_sheet], s, k0, β⃗₀₀, usi[i_sheet])

elseif sheet.class == 'E'

gsmb = pecgsm(length(region[1].P), length(region[end].P))

elseif sheet.class == 'H'

gsmb = pmcgsm(length(region[1].P), length(region[end].P))

else

error("Illegal sheet class: $(sheet.class)")

end

gbldup[ig] < 0 && (gsm_save[ig] = deepcopy(gsmb))

end

# Apply translations if requested:

sheet = sheets[i_sheet]

if sheet.dx ≠ 0 || sheet.dy ≠ 0

upm = ustrip(Float64, sheet.units, 1u"m")

dx = sheet.dx / upm

dy = sheet.dy / upm

translate_gsm!(gsmb, dx, dy, first(region), last(region))

end

else # no sheet

@assert i2 - i1 == 1

gsmb = gsm_slab_interface(layers[i1], layers[i2], k0)

end

t_cascade1 = time()

gsmc = cascade(gsma, gsmb)

cascade!(gsmc, layers[i2])

t_cascade += time() - t_cascade1

gsma = gsmc

end # Gblock loop

t_freq = round(time() - t_freq, digits=tdigits)

t_cascade = round(t_cascade, digits=tdigits)

@logfile " $(t_cascade) seconds for cascading at $(fghz) GHz"

@logfile " $(t_freq) seconds total at $(fghz) GHz"

result = Result(gsmc, steer, β⃗₀₀, fghz, layers[1].ϵᵣ, layers[1].μᵣ,

layers[1].β₁, layers[1].β₂, layers[end].ϵᵣ, layers[end].μᵣ,

layers[end].β₁, layers[end].β₂)

return (β⃗₀₀, result)

end # function

function calculate_jtype_gsm(layers::AbstractVector{Layer}, sheet::RWGSheet, u::Float64,

rwgdat::RWGData, s::Int, k0::Float64, k⃗inc::SVector{2,Float64}, is::Int)

one_meter = ustrip(Float64, sheet.units, 1u"m")

area = norm(sheet.s₁ × sheet.s₂) / one_meter^2 # Unit cell area (m^2).

nmodesmax = max(length(layers[begin].P), length(layers[end].P))

nbf = size(rwgdat.bfe, 2) # Number of basis functions

if isempty(rwgdat.bfftstore)

rwgdat.bfftstore = zeros(SArray{Tuple{2},ComplexF64,1,2}, (nbf, 2, nmodesmax))

end

bfftstore = rwgdat.bfftstore

# Compute area correction factors for the mode normalization constants of

# the two end regions:

tempvec = zeros(2)

for (i, l) in enumerate(@view layers[[begin, end]])

area_i = twopi * twopi / norm(l.β₁ × l.β₂)

tempvec[i] = √(area_i / area)

end

acf = SVector(tempvec[1], tempvec[2])

# Set up the partial GSM due to incident field

(gsm, tlgfvi, vincs) = gsm_electric_gblock(layers, s, k0)

if sheet.style == "NULL"

return gsm

end

# Calculate the scattered field partial scattering matrix of the

# FSS sheet at this junction. Then add it to GSM already computed

# for the dielectric discontinuity...

# Fill the interaction matrix for the current sheet:

t_temp = time()

ψ₁ = k⃗inc ⋅ sheet.s₁ / one_meter

ψ₂ = k⃗inc ⋅ sheet.s₂ / one_meter

@logfile " Beginning matrix fill for sheet $(is)"

zmat = fillz(k0, u, layers, s, ψ₁, ψ₂, sheet, rwgdat)

t_fill = round(time() - t_temp, digits=tdigits)

@logfile " $(t_fill) seconds total matrix fill time for sheet $(is)"

# Factor the matrix:

t_temp = time()

zmatf = lu!(zmat)

t_factor = round(time() - t_temp, digits=tdigits)

@logfile " $(t_factor) seconds to factor matrix for sheet $(is)"

t_temp = time()

# Compute and store the basis function Fourier transforms:

i_ft = 0

for (sr, l) in enumerate(@view layers[[begin, end]]) # loop over possible source regions

for qp = 1:length(l.P)

kvec = l.β[qp]

# If desired F.T. has already been computed, then copy it.

if qp > 1 && kvec ≈ l.β[qp-1]

@inbounds for i in firstindex(bfftstore, 1):lastindex(bfftstore,1)

bfftstore[i, sr, qp] = bfftstore[i, sr, qp-1]

end

continue

end

if sr == 2 && length(layers[1].β) ≥ qp && kvec ≈ layers[1].β[qp]

@inbounds for i in firstindex(bfftstore, 1):lastindex(bfftstore,1)

bfftstore[i, 2, qp] = bfftstore[i, 1, qp]

end

continue

end

bfft = @view bfftstore[:, sr, qp]

rwgbfft!(bfft, rwgdat, sheet, kvec, ψ₁, ψ₂) # Otherwise, compute from scratch

i_ft += 1

end

end

t_fft = round(time() - t_temp, digits=tdigits)

#@logfile " $(t_fft) seconds for basis function Fourier transforms at $(i_ft) points"

nsolve = 0

t_extract = 0.0

i_extract = 0

t_solve = 0.0

imat = rwgdat.rhs

for (sr, ls) in enumerate(@view layers[[begin, end]]) # Loop over source regions

for qp in 1:length(ls.P) # Loop over srce reg modes

# Incident field for source layer in absence of the FSS sheet:

sourcevec = vincs[qp, sr] * ls.c[qp] * acf[sr] * ls.tvec[qp]

# Compute generalized voltage vector:

@inbounds for i in firstindex(bfftstore, 1):lastindex(bfftstore, 1)

imat[i] = bfftstore[i, sr, qp] ⋅ sourcevec # Eq. (7.39)

end

# Solve the matrix equation

t_solve1 = time()

ldiv!(zmatf, imat)

t_solve2 = time()

t_solve += t_solve2 - t_solve1

nsolve += 1

t_extract1 = time()

for (or, lo) in enumerate(@view layers[[begin, end]]) # Loop over obs. regions

smat = gsm[or, sr]

for q in 1:length(lo.P) # Loop obs. regn modes

# Extract partial scattering parameter due to scattered fields...

FTJ = sum((imat[n] * bfftstore[n, or, q] for n in 1:nbf)) # FT of total current

smat[q, qp] -= (lo.tvec[q] ⋅ FTJ) * (tlgfvi[q, or] /

(lo.c[q] * acf[or] * area)) # Eq (6.18)

i_extract += 1

end

end

t_extract2 = time()

t_extract += (t_extract2 - t_extract1)

end

end

t_extract = round(t_extract, digits=tdigits)

@logfile " $(t_extract) seconds to extract $(i_extract) GSM entries"

return gsm

function calculate_mtype_gsm(layers::AbstractVector{Layer}, sheet::RWGSheet, u::Float64,

rwgdat::RWGData, s::Int, k0::Float64, k⃗inc::SVector{2,Float64}, is::Int)

one_meter = ustrip(Float64, sheet.units, 1u"m")

area = norm(sheet.s₁ × sheet.s₂) / one_meter^2 # Unit cell area (m^2).

nmodesmax = max(length(layers[begin].P), length(layers[end].P))

nbf = size(rwgdat.bfe, 2) # Number of basis functions

if isempty(rwgdat.bfftstore)

rwgdat.bfftstore = zeros(SArray{Tuple{2},ComplexF64,1,2}, (nbf, 2, nmodesmax))

end

bfftstore = rwgdat.bfftstore

# Compute area correction factors for the mode normalization constants of

# the two end regions:

tempvec = zeros(2)

for (i, l) in enumerate(@view layers[[begin, end]])

area_i = twopi * twopi / norm(l.β₁ × l.β₂)

tempvec[i] = √(area_i / area)

end

acf = SVector(tempvec[1], tempvec[2])

# Set up the partial GSM due to incident field

(gsm, tlgfiv, iincs) = gsm_magnetic_gblock(layers, s, k0)

if sheet.style == "NULL"

return gsm

end

# Calculate the scattered field partial scattering matrix of the

# FSS sheet at this junction. Then add it to GSM already computed

# for the dielectric discontinuity...

# Fill the interaction matrix for the current sheet:

t_temp = time()

ψ₁ = k⃗inc ⋅ sheet.s₁ / one_meter

ψ₂ = k⃗inc ⋅ sheet.s₂ / one_meter

@logfile " Beginning matrix fill for sheet $(is)"

ymat = filly(k0, u, layers, s, ψ₁, ψ₂, sheet, rwgdat)

t_fill = round(time() - t_temp, digits=tdigits)

@logfile " $(t_fill) seconds total matrix fill time for sheet $(is)"

# Factor the matrix:

t_temp = time()

ymatf = lu!(ymat)

t_factor = round(time() - t_temp, digits=tdigits)

@logfile " $(t_factor) seconds to factor matrix for sheet $(is)"

t_temp = time()

# Compute and store the basis function Fourier transforms:

i_ft = 0

for (sr, l) in enumerate(@view layers[[begin, end]]) # loop over possible source regions

for qp = 1:length(l.P)

kvec = l.β[qp]

# If desired F.T. has already been computed, then copy it.

if qp > 1 && kvec ≈ l.β[qp-1]

@inbounds for i in firstindex(bfftstore, 1):lastindex(bfftstore,1)

bfftstore[i, sr, qp] = bfftstore[i, sr, qp-1]

end

continue

end

if sr == 2 && length(layers[1].β) ≥ qp && kvec ≈ layers[1].β[qp]

@inbounds for i in firstindex(bfftstore, 1):lastindex(bfftstore,1)

bfftstore[i, 2, qp] = bfftstore[i, 1, qp]

end

continue

end

bfft = @view bfftstore[:, sr, qp]

rwgbfft!(bfft, rwgdat, sheet, kvec, ψ₁, ψ₂) # Otherwise, compute from scratch

i_ft += 1

end

end

t_fft = round(time() - t_temp, digits=tdigits)

@logfile " $(t_fft) seconds for basis function Fourier transforms at $(i_ft) points"

nsolve = 0

t_extract = 0.0

t_solve = 0.0

vmat = rwgdat.rhs

i_extract = 0

σ = -1

for (sr, ls) in enumerate(@view layers[[begin, end]]) # Loop over source regions

σ *= -1 # 1 for sr == 1, and -1 for sr == 2

for qp in 1:length(ls.P) # Loop over srce reg modes

# Incident field for Region sr (Eq. (7.64))

sourcevec = iincs[qp, sr] * ls.c[qp] * ls.Y[qp] * (ẑ × ls.tvec[qp])

# Compute generalized current vector:

@inbounds for i in firstindex(bfftstore, 1):lastindex(bfftstore,1)

vmat[i] = bfftstore[i, sr, qp] ⋅ sourcevec # Eq. (7.64)

end

# Solve the matrix equation

t_solve1 = time()

ldiv!(ymatf, vmat)

t_solve2 = time()

t_solve += t_solve2 - t_solve1

nsolve += 1

t_extract1 = time()

for (or, lo) in enumerate(@view layers[[begin, end]]) # Loop over obs. regions

smat = gsm[or, sr]

for q in 1:length(lo.P) # Loop obs. regn. modes

# Extract partial scattering parameter due to scattered fields...

FTM = sum((vmat[n] * bfftstore[n, or, q] for n in 1:nbf)) # FT of total mag. current

smat[q, qp] += ((ẑ × lo.tvec[q]) ⋅ FTM) *

(σ * tlgfiv[q, or] * lo.c[q]) # Eq. (6.37)

i_extract += 1

end

end

t_extract2 = time()

t_extract += t_extract2 - t_extract1

end

end

t_extract = round(t_extract, digits=tdigits)

@logfile " $(t_extract) seconds to extract $(i_extract) GSM entries"

return gsm

function getsttuple(stkeys::Tuple{Symbol,Symbol}, stout::Float64, stin::Float64)

if stkeys[1] ∈ (:phi, :ϕ, :Phi, :PHI, :Φ)

return (θ=stin, ϕ=stout)

elseif stkeys[2] ∈ (:phi, :ϕ, :Phi, :PHI, :Φ)

return (θ=stout, ϕ=stin)

elseif stkeys[1] ∈ (:psi1, :Psi1, :PSI1, :ψ₁, :Ψ₁, :ψ1, :Ψ1)

return (ψ₁=stout, ψ₂=stin)

else

return (ψ₁=stin, ψ₂=stout)

end

function check_inputs(layers, sheets, junc, freqs, stkeys, stvalues, outlist)

# Check that input and output media are lossless and the same

imag(first(layers).ϵᵣ) == imag(last(layers).ϵᵣ) ==

imag(first(layers).μᵣ) == imag(last(layers).μᵣ) == 0 ||

error("First and last layers must be lossless")

(first(layers).ϵᵣ == last(layers).ϵᵣ) &&

(first(layers).μᵣ == last(layers).μᵣ) ||

error("First and last layers must have identical electrical parameters")

# todo:

# Check that ψ₁ and ψ₂ are not specified when there are no nonnull sheets

return

function unique_indices(v::Vector)

n = length(v)

nexti = 0

ui = zeros(Int, n)

for io in 1:n

iszero(ui[io]) || continue

nexti += 1

ui[io] = nexti

for it in io+1:n

iszero(ui[it]) || continue

v[io] === v[it] && (ui[it] = nexti)

end

end

return ui

function get_gbldup(gbls::Vector{Gblock}, layers::Vector{Layer}, sheets::Vector{RWGSheet}, junc::Vector{Int})

gbldup = zeros(Int, length(gbls))

for (g1, gbl1) in pairs(gbls)

(gbldup[g1] ≠ 0 || gbl1.j == 0) && continue

j1, rng1 = gbl1.j, gbl1.rng

n1 = length(layers[first(rng1)].P) # modes at left side

n2 = length(layers[last(rng1)+1].P) # modes at right side

# Examine succeeding Gblocks to see if they match:

for g2 in g1+1:length(gbls)

gbl2 = gbls[g2]

rng2 = gbl2.rng

j2 = gbl2.j

(gbldup[g2] ≠ 0 || j2 == 0) && continue

sheets[junc[j1]] === sheets[junc[j2]] || continue

last(rng1) - j1 ≠ last(rng2) - j2 && continue

# Check that layers within blocks gbl1 and gbl2 are identical:

length(rng1) ≠ length(rng2) && continue

n1 ≠ length(layers[first(rng2)].P) && continue

n2 ≠ length(layers[1+last(rng2)].P) && continue

# Check interior layers:

rng1test = 1+first(rng1):last(rng1)

rng2test = 1+first(rng2):last(rng2)

all(zip(rng1test, rng2test)) do (i1, i2)

layers[i1] == layers[i2]

end || continue

# Check layers bounding the Gblocks:

rng1test = (first(rng1), 1 + last(rng1))

rng2test = (first(rng2), 1 + last(rng2))

all(zip(rng1test, rng2test)) do (i1, i2)

l1 = layers[i1]

l2 = layers[i2]

(l1.ϵᵣ == l2.ϵᵣ) && (l1.μᵣ == l2.μᵣ) &&

(l1.P == l2.P) && (l1.β₁ == l2.β₁) && (l1.β₂ == l2.β₂)

end || continue

# If we made it to here, the two Gblocks are identical:

gbldup[g1] = -1 # Indicate that GSM of Gblock g1 is to be saved

gbldup[g2] = g1 # GSM of Gblock g2 is obtained from saved GSM of block g1

end

end

return gbldup

end # function

function report_layers_sheets(layers, sheets, junc, rwgdat, usi)

@logfile "Dielectric layer information... \n"

@logfile " Layer Width units epsr tandel mur mtandel modes beta1x beta1y beta2x beta2y"

@logfile " ----- ------------- ------- ------ ------- ------ ----- ------- ------- ------- -------"

for (jl, l) in pairs(layers)

eps = real(l.ϵᵣ)

tandel = -imag(l.ϵᵣ) / eps

mu = real(l.μᵣ)

mtandel = -imag(l.μᵣ) / mu

nmode = length(l.P)

units = string(unit(l.user_width))

units == "inch" && (units = "in")

uw_unitless = ustrip(l.user_width)

str = @sprintf(" %5i %9.4f %3s %7.2f %6.4f %7.2f %6.4f %5i %7.1f %7.1f %7.1f %7.1f",

jl, uw_unitless, units, eps, tandel, mu, mtandel, nmode, l.β₁[1], l.β₁[2],

l.β₂[1], l.β₂[2])

@logfile "$str"

if jl < length(layers) && junc[jl] ≠ 0

js = junc[jl]

s = sheets[js]

om = ustrip(Float64, s.units, 1.0u"m")

str = @sprintf(

" ================== Sheet %3i ======================== %7.1f %7.1f %7.1f %7.1f",

usi[js], s.β₁[1] * om, s.β₁[2] * om, s.β₂[1] * om, s.β₂[2] * om)

@logfile "$str"

end

end

sint = findall(junc .≠ 0) # sint[k] contains dielectric interface number of k'th sheet

@logfile "\n\n\nPSS/FSS sheet information...\n"

@logfile "Sheet Loc Style Rot J/M Faces Edges Nodes Unknowns NUFP"

@logfile "----- --- ---------------- ----- --- ----- ----- ----- -------- ------"

for (js, s) in pairs(sheets)

str = @sprintf("%4i %3i %16s %5.1f %1s %5i %5i %5i %6i %7i",

usi[js], sint[js], s.style, s.rot, s.class, size(s.fe, 2), length(s.e1),

length(s.ρ), size(rwgdat[js].bfe, 2), length(rwgdat[js].ufp2fp))

@logfile "$str"

end

@logfile "\n\n"

nothing

@setup_workload begin

# Putting some things in `setup` can reduce the size of the

# precompile file and potentially make loading faster.

outer(rot) = meander(a=3.97, b=3.97, w1=0.13, w2=0.13, h=2.53+0.13, units=mm, ntri=30, 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=30, rot=rot, class='M')

#center(rot) = meander(a=3.97, b=3.97, w1=0.34, w2=0.34, h=2.51+0.34, units=mm, ntri=300, rot=rot)

t1 = 4

t2 = 2.45

foam(w) = Layer(width=w, epsr=1.05)

#(r::Int, uρ⃗₀₀::SV2, us₁::SV2, us₂::SV2, ψ₁::Float64, ψ₂::Float64, tid::Int) = (1, [1.7453292519943293, 28.099800957108705], [3.141592653589793, 0.0], [0.0, 62.831853071799586], 0.0, 0.0, 4)

substrate = Layer(width=0.1mm, epsr=2.6)

strata = [

Layer()

outer(0)

substrate

foam(t1*mm)

inner(0)

substrate

Layer(epsr = 2.3, width = 5mm)

Layer() ]

steering = (θ=0:1, ϕ=0)

flist = 10

@compile_workload begin

resultfile = tempname()

logfile = tempname()

results = redirect_stdout(devnull) do

results = analyze(strata, flist, steering; resultfile, logfile)

end

RL11rr = -extract_result(results, @outputs s11db(r,r))

AR11r = extract_result(results, @outputs ar11db(r))

IL21L = -extract_result(results, @outputs s21db(L,L))

AR21L = extract_result(results, @outputs ar21db(L))

RL11rr = -extract_result(resultfile, @outputs s11db(r,r))

AR11r = extract_result(resultfile, @outputs ar11db(r))

IL21L = -extract_result(resultfile, @outputs s21db(L,L))

AR21L = extract_result(resultfile, @outputs ar21db(L))

end

end # module