diff --git a/Project.toml b/Project.toml index 04111d57..2cdb499b 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "PSSFSS" uuid = "6b20a5d4-3c6c-44cd-883b-1480592d72be" authors = ["Peter Simon "] -version = "1.8.2" +version = "1.8.3" [deps] Dates = "ade2ca70-3891-5945-98fb-dc099432e06a" @@ -36,13 +36,13 @@ Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d" FFTW = "1.8" FileIO = "1.16.2" GeoInterface = "1.3.3" -JLD2 = "0.4.46" +JLD2 = "0.5.8" LibGEOS = "0.9" LoggingExtras = "1.0.3" MetalSurfaceImpedance = "1.0" NearestNeighbors = "0.4.16" OffsetArrays = "1.13" -OhMyThreads = "0.5" +OhMyThreads = "0.7" PkgVersion = "0.3.3" PolygonOps = "0.1.2" PrecompileTools = "1.2" diff --git a/src/Modes.jl b/src/Modes.jl index ff1c4825..f00014fc 100644 --- a/src/Modes.jl +++ b/src/Modes.jl @@ -276,7 +276,7 @@ function choose_layer_modes!(layers::Vector{Layer}, sheets::Vector{RWGSheet}, ju end """ - setup_modes!(layer::Layer, k0::Real, kvec::AbstractVector) + setup_modes!(layer::Layer, k0::Real, kvec::AbstractVector, β̂default::AbstractVector) Fill the modal layer fields. Needed for layers not contained in a Gblock. The arrays are assumed to have been already allocated, and the index arrays @@ -290,13 +290,14 @@ The arrays are assumed to have been already allocated, and the index arrays - `kvec`: A real-valued 2-vector containing the x and y components of the incident plane wave unit vector that defines the unit cell incremental phase shifts. +- `β̂default`: Region 1 incident unit tangential wave vector that depends correctly on ϕ in all cases. ### Outputs There is no explicit output, but the fields of `layer` will be modified, including `β`, `tvec`, `γ`, `c`, and `Y`. """ -function setup_modes!(layer::Layer, k0::Real, kvec::AbstractVector) +function setup_modes!(layer::Layer, k0::Real, kvec::AbstractVector, β̂default::AbstractVector) β₀₀ = SVector(kvec[1], kvec[2]) β₁, β₂ = layer.β₁, layer.β₂ area = twopi^2 / norm(β₁ × β₂) @@ -305,7 +306,11 @@ function setup_modes!(layer::Layer, k0::Real, kvec::AbstractVector) m, n, p = layer.M[mode], layer.N[mode], layer.P[mode] β = β₀₀ + m * β₁ + n * β₂ β² = β ⋅ β - β̂ = β² * area < 1e-14 ? SVector(1.0, 0.0) : β / norm(β) + if β² * area < 1e-14 + β̂ = copy(β̂default) + else + β̂ = β / norm(β) + end layer.β[mode] = β layer.γ[mode] = γ = mysqrt(β² - ksq) if p == TE diff --git a/src/Outputs.jl b/src/Outputs.jl index ceb703e8..21b94da0 100644 --- a/src/Outputs.jl +++ b/src/Outputs.jl @@ -110,7 +110,7 @@ function sourcemat(j::Int, n::HorV, o::Result) ĥ = @view ĥ3[1:2] # Only need x and y components due to dot product later v̂ = @view v̂3[1:2] # Only need x and y components due to dot product later β₀₀ = norm(o.β⃗₀₀) - β̂₀₀ = (β₀₀ == 0) ? @SVector([1.0, 0.0]) : o.β⃗₀₀ / β₀₀ + β̂₀₀ = (β₀₀ == 0) ? @SVector([cosd(ϕ1inc), sind(ϕ1inc)]) : o.β⃗₀₀ / β₀₀ t̂₁ = ẑ × β̂₀₀ t̂₂ = β̂₀₀ ct = cosd(θ) @@ -134,7 +134,7 @@ function sourcemat(j::Int, n::RorL, o::Result) L̂ = view((ĥ + sgn * im * v̂) / √2, 1:2) # Only need x and y components due to dot product later R̂ = view((ĥ - sgn * im * v̂) / √2, 1:2) β₀₀ = norm(o.β⃗₀₀) - β̂₀₀ = (β₀₀ == 0) ? @SVector([1.0, 0.0]) : o.β⃗₀₀ / β₀₀ + β̂₀₀ = (β₀₀ == 0) ? @SVector([cosd(ϕ1inc), sind(ϕ1inc)]) : o.β⃗₀₀ / β₀₀ t̂₁ = ẑ × β̂₀₀ t̂₂ = β̂₀₀ ct = cosd(θ) @@ -165,7 +165,7 @@ function obsmat(i::Int, n::HorV, o::Result) end (ĥ, v̂) = ĥv̂(θ, ϕ) β₀₀ = norm(o.β⃗₀₀) - β̂₀₀ = (β₀₀ == 0) ? @SVector([1.0, 0.0]) : o.β⃗₀₀ / β₀₀ + β̂₀₀ = (β₀₀ == 0) ? @SVector([cosd(ϕ1inc), sind(ϕ1inc)]) : o.β⃗₀₀ / β₀₀ t̂₁2 = ẑ × β̂₀₀ t̂₁ = @SVector([t̂₁2[1], t̂₁2[2], 0.0]) t̂₂ = @SVector([β̂₀₀[1], β̂₀₀[2], sgn * tand(θ)]) # term from Eqs. (8.20) @@ -190,7 +190,7 @@ function obsmat(i::Int, n::RorL, o::Result) L̂ = (ĥ - sgn * im * v̂) / √2 R̂ = (ĥ + sgn * im * v̂) / √2 β₀₀ = norm(o.β⃗₀₀) - β̂₀₀ = (β₀₀ == 0) ? @SVector([1.0, 0.0]) : o.β⃗₀₀ / β₀₀ + β̂₀₀ = (β₀₀ == 0) ? @SVector([cosd(ϕ1inc), sind(ϕ1inc)]) : o.β⃗₀₀ / β₀₀ t̂₁2 = ẑ × β̂₀₀ t̂₁ = @SVector([t̂₁2[1], t̂₁2[2], 0.0]) t̂₂ = @SVector([β̂₀₀[1], β̂₀₀[2], sgn * tand(θ)]) # term from Eqs. (8.20) @@ -336,6 +336,7 @@ function θϕ(o::Result) β₀₀² == 0 && return (0.0, get(o.steering, :ϕ, 0.0)) k² = (twopi * o.FGHz * 1e9 / c₀)^2 * real(o.ϵᵣin * o.μᵣin) β₀₀² > k² && error("Cut-off dominant mode") + haskey(o.steering, :θ) && return (o.steering.θ, o.steering.ϕ) kz = √(k² - β₀₀²) # for out-going wave vector in Layer 1 θ = acosd(kz / sqrt(k²)) ϕ = atand(o.β⃗₀₀[2], o.β⃗₀₀[1]) diff --git a/src/PSSFSS.jl b/src/PSSFSS.jl index e6b5140a..2b5ef7fc 100644 --- a/src/PSSFSS.jl +++ b/src/PSSFSS.jl @@ -263,7 +263,7 @@ function _analyze(layers, sheets, junc, freqs, stkeys, stvalues; for stout in stvalues[1], stin in stvalues[2] steer = getsttuple(stkeys, stout, stin) if keys(steer)[1] == :ψ₁ - ψ₁, ψ₂ = deg2rad.(steer) # radians + ψ₁, ψ₂ = (deg2rad(x) for x in steer) # radians upm::Float64 = ustrip(Float64, sheets[1].units, 1u"m") β₁, β₂ = sheets[1].β₁ * upm, sheets[1].β₂ * upm β⃗₀₀k1 = (ψ₁ * β₁ + ψ₂ * β₂) / twopi # Eq. (2.13b) @@ -337,12 +337,15 @@ function compute_next_freq(fghz, β⃗₀₀k1, steer, layers, sheets, usi, rwgd if keys(steer)[1] == :θ k1 = k0 * sqrt(real(layers[1].ϵᵣ * layers[1].μᵣ)) β⃗₀₀ = k1 * β⃗₀₀k1 + sϕdefault, cϕdefault = sincosd(last(steer)) + β̂default = SVector(cϕdefault, sϕdefault) # Needed in case iszero(θ) else β⃗₀₀ = copy(β⃗₀₀k1) + β̂default = SVector(1.0, 0.0) # Needed in case iszero(β⃗₀₀) end t_cascade = 0.0 - setup_modes!.(layers, k0, Ref(β⃗₀₀)) + setup_modes!.(layers, k0, Ref(β⃗₀₀), Ref(β̂default)) if !(angle(layers[begin].γ[1]) ≈ angle(layers[end].γ[1]) ≈ π / 2) @logfile " Skipping $(fghz) GHz due to cutoff principal modes in ambient medium" return nothing