examples/analyze_turbulent_channel_model.jl

@@ -156,21 +156,33 @@ function scintillation_index(fn; basedir=pwd(), debug=false, showplot=false)
     λ = metadata["lambda"]
     Dz = metadata["Dz"]
     Cn² = metadata["Cn2"]
+    w0 = metadata["w0"]
+    l0 = metadata["l0"]
+    L0 = metadata["L0"]
 
     # scintillation index along highest dimension, note that the dumped fields are not fftshifted!
     SI = squeeze(scv(urx; dims=ndims(urx)))
     SI0 = CUDA.@allowscalar first(SI) # without the fftshift, the first element is the on-axis one
 
     if debug
+        # compute theoretical reference values
+        SI_pw = scintillation_index_pw(2π / λ, Dz, Cn², l0, L0)
+        SI_sw = scintillation_index_sw(2π / λ, Dz, Cn², l0, L0)
+        SI_gaussian = scintillation_index_gaussian(2π / λ, Dz, Cn², w0, l0, L0)
+        
         @info "Estimated on-axis scintillation index: " SI0
+        @info "Theoretical result (plane wave): " SI_pw
+        @info "Theoretical result (spherical wave): " SI_sw
+        @info "Theoretical result (gaussian beam): " SI_gaussian
 
-        # for weak fluctuations, the scintillation index leads to the Rytov variance
-        σ²ᵣ_pw = rytov_pw(2π / λ, Dz, Cn²)
-        @info "Rytov variance (plane wave): " σ²ᵣ_pw
-        σ²ᵣ_sw = rytov_sw(2π / λ, Dz, Cn²)
-        @info "Rytov variance (spherical wave): " σ²ᵣ_sw
+        # check convergence on axis
+        SI_evo = squeeze(cumscv(urx[1, 1, :]))
+        y = Adapt.adapt(Array, SI_evo)
+        p = plot.scatter(y=y)
+        p.layout.xaxis.title = "Number of iterations"
+        p.layout.yaxis.title = "Estimated scintillation index"
 
-        @info "Theoretical scintillation index (PW, Kolmogorov): " exp(0.49σ²ᵣ_pw / (1 + 1.11σ²ᵣ_pw^(6 // 5))^(7 // 6) + 0.51σ²ᵣ_pw / (1 + 0.69σ²ᵣ_pw^(6 // 5))^(5 // 6)) - 1
+        display(p)
     end
     
     if showplot

examples/batch_analysis.jl

@@ -107,17 +107,16 @@ function batch_scintillation_index(dirname)
 
     # gather all hdf5 files in the given directory `dirname`
     files = readdir(dirname, join=true)
-    csvidx = findfirst(contains("xvals.csv"), files)
     h5files = filter(x -> splitext(x)[2] == ".hdf5", files)
 
     Nf = length(h5files)
 
-    # gather xvals from corresponding csv if possible - otherwise 1:N
-    df = if isnothing(csvidx)
-        DataFrame(:File => range(1, Nf))
-    else
-        CSV.read(files[csvidx], DataFrame)
+    # gather xvals from corresponding csv if possible - otherwise generate xvals.csv and fill with 1:N
+    xvals_fn = joinpath(dirname, "xvals.csv")
+    !isfile(xvals_fn) && open(xvals_fn, "w") do fid
+        CSV.write(fid, DataFrame(:file => range(1, Nf)))
     end
+    df = CSV.read(xvals_fn, DataFrame)
 
     SI_est = zeros(Nf)
     SI_pw = zeros(Nf)
@@ -133,7 +132,9 @@ function batch_scintillation_index(dirname)
 
         # make sure the order the files are read matches the xvals in the DataFrame
         xname = first(names(df))
-        @assert df[k, xname] == metadata[string(xname)]
+        if xname != "file"
+            @assert df[k, xname] == metadata[string(xname)]
+        end
 
         Cn² = metadata["Cn2"]
         λ = metadata["lambda"]

src/estimate_statistics.jl

@@ -232,6 +232,12 @@ function scv(u::AbstractArray{T,N}; dims=:) where {T,N}
     mean(abs2, I; dims) ./ mean(I; dims) .^ 2 .- 1 # it might be advantageous to not use `var` to avoid computing the mean twice
 end
 
+function cumscv(u::AbstractArray{T,N}; dims=N) where {T,N}
+    I = abs2.(u)
+    I² = I .^ 2
+    cumsum(I²; dims) .* (1:size(I, dims))./ cumsum(I; dims) .^ 2 .- 1
+end
+
 
 function estimate_ccdf(v)
     # local ccdf