Fix normalinc betahat (#56)

* Correct computation of β̂ when theta is zero and phi is not

* Update compats and bump patch version

* Correct calculation of β̂

* Correct betahat unit vector calc

Project.toml

@@ -1,7 +1,7 @@
 name = "PSSFSS"
 uuid = "6b20a5d4-3c6c-44cd-883b-1480592d72be"
 authors = ["Peter Simon <psimon0420@gmail.com>"]
-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"

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

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])

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