feat: proper theoretical SI references

Project.toml

@@ -15,6 +15,7 @@ FileIO = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549"
 FourierTools = "b18b359b-aebc-45ac-a139-9c0ccbb2871e"
 H5Zblosc = "c8ec2601-a99c-407f-b158-e79c03c2f5f7"
 HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f"
+HypergeometricFunctions = "34004b35-14d8-5ef3-9330-4cdb6864b03a"
 Optim = "429524aa-4258-5aef-a3af-852621145aeb"
 OrderedCollections = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
 PlotlyKaleido = "f2990250-8cf9-495f-b13a-cce12b45703c"

dev/fix_scintillation_index.jl

@@ -75,14 +75,13 @@ function analyze_batch_scintillation_index(dirname; xscale="linear", yscale="lin
     # choose column of xval if multiple parameters have been swept simultaneously
     xvals = xval_df[!, names(xval_df)[dfcol]]
 
-    σ² = zeros(Nf)
-    σ²ᵣ_pw = zeros(Nf)
-    σ²ᵣ_sw = zeros(Nf)
+    SI_est = zeros(Nf)
     SI_pw = zeros(Nf)
+    SI_sw = zeros(Nf)
     SI_gaussian = zeros(Nf)
 
     for (k, file) ∈ enumerate(h5files)
-        σ²[k] = scintillation_index(file)
+        SI_est[k] = scintillation_index(file)
 
         metadata = h5open(file, "r") do h5
             Dict(attrs(h5))
@@ -92,25 +91,20 @@ function analyze_batch_scintillation_index(dirname; xscale="linear", yscale="lin
         λ = metadata["lambda"]
         Dz = metadata["Dz"]
         w0 = metadata["w0"]
-        σ²ᵣ_pw[k] = rytov_pw(2π / λ, Dz, Cn²)
-        σ²ᵣ_sw[k] = rytov_sw(2π / λ, Dz, Cn²)
-        SI_pw[k] = exp(0.49σ²ᵣ_pw[k] / (1 + 1.11σ²ᵣ_pw[k]^(6 // 5))^(7 // 6) + 0.51σ²ᵣ_pw[k] / (1 + 0.69σ²ᵣ_pw[k]^(6 // 5))^(5 // 6)) - 1
-        # Gaussian, according to Andrews 2019 p.54
-        Λ0 = λ * Dz / (π * w0^2) # inverse Rayleigh length
-        Θ0 = 1 # 1 - Dz / F0 # F0 is the radius of curvature at z=0, we consider a collimated beam, i.e., F0 = ∞
-        Θ = Θ0 / (Θ0^2 + Λ0^2)
-        Λ = Λ0 / (Θ0^2 + Λ0^2)
-        σ²B = 3.86 * σ²ᵣ_pw[k] * ((0.4 * (1 + 2Θ)^2 + 4Λ^2)^(5 // 12) * cos(5 / 6 * atan(1 + 2Θ, 2Λ)) - 11 / 16 * Λ^(5 / 6))
-        SI_gaussian[k] = exp(0.49σ²B / (1 + 0.56 * (1 + Θ) * σ²B^(6 // 5))^(7 // 6) + 0.51σ²B / (1 + 0.69σ²B^(6 // 5))^(5 // 6)) - 1
+        
+        l0 = metadata["l0"]
+        L0 = metadata["L0"]
+
+        SI_pw[k] = scintillation_index_pw(2π / λ, Dz, Cn², l0, L0)
+        SI_sw[k] = scintillation_index_sw(2π / λ, Dz, Cn², l0, L0)
+        SI_gaussian[k] = scintillation_index_gaussian(2π / λ, Dz, Cn², w0, l0, L0)
     end
 
     # visualize results
-    plt = plot.scatter(x=xvals, y=σ², name="Estimated SI")
-    plt.scatter(x=xvals, y=σ²ᵣ_pw, line=Config(color="black"), name="Rytov variance - PW")
-    plt.scatter(x=xvals, y=SI_pw, line=Config(color="gray"), name="SI - PW")
-    plt.scatter(x=xvals, y=σ²ᵣ_sw, line=Config(color="black"), name="Rytov variance - SW")
+    plt = plot.scatter(x=xvals, y=SI_est, name="Estimated SI")
     plt.scatter(x=xvals, y=SI_gaussian, name="SI - Gaussian")
-    # plt = plot.scatter(x=xvals, y=σ² ./ SI_pw, name="Estimated SI vs. PW-SI")
+    plt.scatter(x=xvals, y=SI_pw, line=Config(color="black"), name="SI - PW")
+    plt.scatter(x=xvals, y=SI_sw, line=Config(color="gray"), name="SI - SW")
     # set lower yaxis limit to 0
     yl = plt.layout.yaxis.range
     yl[1] = 0

src/optical_parameters.jl

+using HypergeometricFunctions
 using LinearAlgebra
 using Optim
 
@@ -38,11 +39,135 @@ fried_sw_layered(k0, Δz::AbstractVector, Cn²::AbstractVector) = begin
 end
 
 # Rytov variance
-# for weak fluctuations, the scintillation index leads to the Rytov variance
+# for weak fluctuations, the scintillation index corresponds to the Rytov variance
 rytov_pw(k0, z, Cn²) = 1.23Cn² * k0^(7 // 6) * z^(11 // 6)
 rytov_sw(k0, z, Cn²) = 0.5Cn² * k0^(7 // 6) * z^(11 // 6)
 
 
+# Scintillation index
+# plane wave - compare Andrews 2019 - Field Guide ... p. 50/51 : the inner and outer scale are considered assuming the Andrews spectrum
+scintillation_index_pw(k0, Δz, Cn², l0=0, L0=Inf; extended=true) = begin
+
+    σ²_R = rytov_pw(k0, Δz, Cn²)
+
+    κl = 3.3 / l0
+    Qₗ = κl^2 * Δz / k0
+    Q̂₀ = 64π^2 * Δz / (k0 * L0^2)
+
+    # p. 51
+    c1(Qₗ) = 2.61 / (2.61 + Qₗ + 0.45σ²_R * Qₗ^(7 // 6))
+    σ²_lnX(Qₗ) = 0.16σ²_R * (Qₗ * c1(Qₗ))^(7 // 6) * (1 + 1.75 * sqrt(c1(Qₗ)) - 0.25 * c1(Qₗ)^(7 // 12))
+
+    c2(Qₗ, Q̂₀) = 2.61Q̂₀ / (2.61 * (Q̂₀ + Qₗ) + Q̂₀ * Qₗ * (1 + 0.45σ²_R * Qₗ^(1 // 6)))
+    σ²_lnX(Qₗ, Q̂₀) = 0.16σ²_R * (Qₗ * c2(Qₗ, Q̂₀))^(7 // 6) * (1 + 1.75 * sqrt(c2(Qₗ, Q̂₀)) - 0.25 * c2(Qₗ, Q̂₀)^(7 // 12))
+
+    if iszero(l0)
+        (L0 < Inf) && @warn "Outer scale L0=$L0 ignored."
+        return exp(0.49σ²_R / (1 + 1.11σ²_R^(6 // 5))^(7 // 6) + 0.51σ²_R / (1 + 0.69σ²_R^(6 // 5))^(5 // 6)) - 1
+    else
+        # bottom of p. 50
+        c3 = atan(Qₗ)
+        σ²_pl = 3.86σ²_R * ((1 + 1 / Qₗ^2)^(11 // 12) * (sin(11 // 6 * c3) + 1.507 / (1 + Qₗ^2)^(1 // 4) * sin(4 // 3 * c3) - 0.273 / (1 + Qₗ^2)^(7 // 24) * sin(5 // 4 * c3)) - 3.50 / Qₗ^(5 // 6))
+
+        if extended
+            return exp(σ²_lnX(Qₗ) - σ²_lnX(Qₗ, Q̂₀) + (0.51 * σ²_pl) / (1 + 0.69 * σ²_pl^(6 // 5))^(5 // 6)) - 1
+        else
+            return σ²_pl
+        end
+    end
+end
+
+# the inner and outer scale are considered assuming the Andrews spectrum
+scintillation_index_sw(k0, Δz, Cn², l0=0, L0=Inf; extended=true) = begin
+
+    β₀² = rytov_sw(k0, Δz, Cn²)
+
+    κl = 3.3 / l0
+    Qₗ = κl^2 * Δz / k0
+    Q̂₀ = 64π^2 * Δz / (k0 * L0^2)
+
+    # p. 53
+    c1(Qₗ) = 8.56 / (8.56 + Qₗ + 0.20β₀² * Qₗ^(7 // 6))
+    σ²_lnX(Qₗ) = 0.04β₀² * (c1(Qₗ) * Qₗ)^(7 // 6) * (1 + 1.75 * sqrt(c1(Qₗ)) - 0.25 * c1(Qₗ)^(7 // 12))
+
+    c2(Qₗ, Q̂₀) = 8.56Q̂₀ / (8.56 * (Q̂₀ + Qₗ) + Q̂₀ * Qₗ * (1 + 0.20β₀² * Qₗ^(1 // 6)))
+    σ²_lnX(Qₗ, Q̂₀) = 0.04β₀² * (Qₗ * c2(Qₗ, Q̂₀))^(7 // 6) * (1 + 1.75 * sqrt(c2(Qₗ, Q̂₀)) - 0.25 * c2(Qₗ, Q̂₀)^(7 // 12))
+
+    if iszero(l0)
+        # (L0 < Inf) && @warn "Outer scale L0=$L0 ignored."
+        return exp(0.49β₀² / (1 + 0.56β₀²^(6 // 5))^(7 // 6) + 0.51β₀² / (1 + 0.69β₀²^(6 // 5))^(5 // 6)) - 1
+    else
+        # bottom of p. 52
+        c3 = atan(Qₗ, 3)
+        σ²_sph = 9.65β₀² * (0.40 * (1 + 9 / Qₗ^2)^(11 // 12) * (sin(11 // 6 * c3) + 2.610 / (9 + Qₗ^2)^(1 // 4) * sin(4 // 3 * c3) - 0.518 / (9 + Qₗ^2)^(7 // 24) * sin(5 // 4 * c3)) - 3.50 / Qₗ^(5 // 6))
+
+        # make sure we do not end up with negative variances ... they might be an artifact from the (propbably) approximative nature of the expression for σ²_sph
+        if σ²_sph < 0
+            return zero(σ²_sph)
+        end
+        if extended
+            return exp(σ²_lnX(Qₗ) - σ²_lnX(Qₗ, Q̂₀) + (0.51 * σ²_sph) / (1 + 0.69 * σ²_sph^(6 / 5))^(5 / 6)) - 1
+        else
+            return σ²_sph
+        end
+    end
+end
+
+# on-axis SI for Gaussian beam: the inner and outer scale are considered assuming the Andrews spectrum
+scintillation_index_gaussian(k0, Δz, Cn², w0, l0=0, L0=Inf; extended=true) = begin
+
+    σ²_R = rytov_pw(k0, Δz, Cn²)
+
+    κl = 3.3 / l0
+    Qₗ = κl^2 * Δz / k0
+    Q̂₀ = 64π^2 * Δz / (k0 * L0^2)
+
+    # beam parameters in Andrew's weird notation
+    Λ0 = 2Δz / (k0 * w0^2) # inverse Rayleigh length
+    Θ0 = 1 # 1 - Dz / F0 # F0 is the radius of curvature at z=0, we consider a collimated beam, i.e., F0 = ∞
+    Θ = Θ0 / (Θ0^2 + Λ0^2)
+    Λ = Λ0 / (Θ0^2 + Λ0^2)
+    Θ̄ = 1 - Θ # compare p. 19
+    W = w0 * sqrt(Θ0^2 + Λ0^2) # compare p. 18
+
+    # p. 56
+    η_X = inv(0.38 / (1 - 3.21Θ̄ + 5.29Θ̄^2) + 0.47σ²_R * Qₗ^(1 // 6) * ((1 // 3 - Θ̄ / 2 + Θ̄^2 / 5) / (1 + 2.20Θ̄))^(6 // 7))
+    η_X0 = η_X * Q̂₀ / (η_X + Q̂₀)
+
+    σ²_lnX(Qₗ, η_X0) = 0.49σ²_R * (1 // 3 - Θ̄ / 2 + Θ̄^2 / 5) * (η_X0 * Qₗ / (η_X0 + Qₗ))^(7 // 6) * (1 + 1.75 * sqrt(η_X0 / (η_X0 + Qₗ)) - 0.25 * (η_X0 / (η_X0 + Qₗ))^(7 // 12))
+    σ²_lnX(Qₗ) = σ²_lnX(Qₗ, η_X)
+
+    # compute scintillation index for Gaussian beam assuming Kolmogorov spectrum and weak fluctuations
+    # exact solution, holds only for r < W!
+    σ²_I(r) = 3.86σ²_R * real(cispi(5 // 12) * _₂F₁(-5 // 6, 11 // 6, 17 // 6, Θ̄ + 1im * Λ)) - 2.65σ²_R * Λ^(5 // 6) * _₁F₁(-5 // 6, 1, 2r^2 / W^2)
+    # approximation
+    # σ²_I(r) = 3.86σ²_R * (0.4 * ((1 + 2Θ)^2 + 4Λ^2)^(5 // 12) * cos(5 // 6 * atan(1 + 2Θ, 2Λ)) - 11 // 16 * Λ^(5 // 6)) + 4.42σ²_R * Λ^(5 // 6) * r^2 / W^2
+    σ²_B = σ²_I(0)
+
+    # the scintillation index of the Gaussian beam is composed of two parts, a radial and a longitudinal component. However, for r=0, i.e., on axis, the radial component vanishes, compare equation at the top of p. 56
+    if iszero(l0)
+        (L0 < Inf) && @warn "Outer scale L0=$L0 ignored."
+        # extend validity into strong turbulence regime
+        return exp(0.49σ²_B / (1 + 0.56 * (1 + Θ) * σ²_B^(6 // 5))^(7 // 6) + 0.51σ²_B / (1 + 0.69σ²_B^(6 // 5))^(5 // 6)) - 1
+    else
+        # bottom of p. 55
+        φ₁ = atan(2Λ, 1 + 2Θ)
+        φ₂ = atan((1 + 2Θ) * Qₗ, 3 + 2Λ * Qₗ)
+        # upper half of p. 55
+
+        σ²_G = 3.86σ²_R * (0.40 * ((1 + 2Θ)^2 + (2Λ + 3 / Qₗ)^2)^(11 // 12) / sqrt((1 + 2Θ)^2 + 4Λ^2) *
+                           (sin(11 / 6 * φ₂ + φ₁) + 2.610 / ((1 + 2Θ)^2 * Qₗ^2 + (3 + 2Λ * Qₗ)^2)^(1 / 4) * sin(4 / 3 * φ₂ + φ₁)
+                            - 0.518 / ((1 + 2Θ)^2 * Qₗ^2 + (3 + 2Λ * Qₗ)^2)^(7 / 24) * sin(5 / 4 * φ₂ + φ₁))
+                           - 13.401Λ / (Qₗ^(11 / 6) * ((1 + 2Θ)^2 + 4Λ^2)) - 11 / 6 * Qₗ^(-5 / 6) * ((1 + 0.31Λ * Qₗ)^(5 / 6) + 1.096 * (1 + 0.27Λ * Qₗ)^(1 / 3) - 0.186 * (1 + 0.24Λ * Qₗ)^(1 / 4)))
+        if extended
+            exp(σ²_lnX(Qₗ) - σ²_lnX(Qₗ, η_X0) + (0.51 * σ²_G) / (1 + 0.69 * σ²_G^(6 / 5))^(5 / 6)) - 1
+        else
+            σ²_G
+        end
+    end
+end
+
+
 # log-amplitude variances (evaluate integrals (9.63) and (9.64) with wolframalpha, where gamma(17/6) / gamma(7/3) ≈ 1.448411512633859)
 sigma_chi_pw(k0, Δz, Cn²) = 0.563 * k0^(7 / 6) * Δz^(5 / 6) * Cn² * (6Δz / 11)
 sigma_chi_sw(k0, Δz, Cn²) = 0.563 * k0^(7 / 6) * Cn² * (3 * sqrt(π) / 22 / 2^(2 / 3) * Δz^(11 / 6) * 1.448411512633859)