dev/fix_scintillation_index.jl

@@ -8,7 +8,7 @@ include(joinpath("..", "examples", "analyze_turbulent_channel_model.jl"))
 
     xvals has to be a vector/range/any iterable
 """
-function batch_scintillation_index(args=Dict(:Nz => 2:11, :dz => 1e3 ./ (1:10)); iters=100, λ=1.55e-6, showprogress=true, samplingcheck=false, kwargs...)
+function batch_scintillation_index(args=Dict(:Nz => 4:16, :dz => 1e3 ./ (3:15)); iters=100, λ=1.55e-6, showprogress=true, samplingcheck=false, kwargs...)
 
     # make sure, that all value arrays of the argument dict are of the same size
     valuelen = length.(values(args))

dev/pointsrc.jl

-include(joinpath("..", "examples", "repeated_wave_propagation.jl"))
-
-function pointsrc(; Nx=256, dx1=1e-3, dz=1e1, λ=1.55e-6, dx2=λ * dz / (Nx * dx1), Nplt=256)
-
-    # The bandwidth of the point source is 1/dx1, i.e., the largest spatial frequency is 1/2dx1 which corresponds to  cos(α) = kx / k0 = λ/2dx1. Accordingly, we find tan(α) = kz / kx. Hence, after traversing a distance dz, the spatial extent of the point source has increased by 2dz / tan(α) = 2dz / tan(arccos(λ / 2dx1)) = 2dz / √((2dx1 / λ)^2 - 1) ≈ dz * λ / dx1
-    D1 = dx1
-    D2 = 0.75dz * λ / dx1
-    @info "Observation plane aperture size D2=$D2"
-    
-    # Check the angular spectrum method sampling
-    check_sampling(dz, Nx, dx1, D1, dx2, D2; λ, R=Inf, debug=true)
-
-    # Initialize the field with a point source at the origin
-    u = zeros(ComplexF64, Nx, Nx)
-    u[1] = 1
-
-    # Propagate the field using the two-step angular spectrum method
-    twostep_prop_ang!(u, dx1, dx2, 2π / λ, dz)
-    us = fftshift(u)
-
-    # Plot the result
-    x0 = (Nx > Nplt) ? -Nplt / 2 * dx2 : -Nx / 2 * dx2
-    uplt = (Nx > Nplt) ? center_extract(us, (Nplt, Nplt)) : us
-    Plot(complex_plot(; x0=x0, y0=x0, dx=dx2, dy=dx2, z=uplt))
-end
\ No newline at end of file

examples/repeated_wave_propagation.jl

 include(joinpath("..", "src", "sampling_requirements.jl"))
 include(joinpath("..", "src", "turbulent_channel.jl"))
 include(joinpath("..", "src", "monte_carlo.jl"))
+include(joinpath("..", "src", "source.jl"))
 
 include("plot_results.jl")
 
@@ -54,7 +55,7 @@ function repeated_wave_propagation(utx::AbstractArray{Complex{T}}, midx, domain:
     fftplan = plan_fft!(utx[:, :, 1])
 
     for (j, chunk) ∈ enumerate(Iterators.partition(1:iters, chunksize))
-        urx = maybe_cuda(zeros(ComplexF64, size(utx)..., length(chunk)))
+        urx = maybe_cuda(zeros(Complex{T}, size(utx)..., length(chunk)))
         # monte carlo loop
         mc_iterations!(urx, utx, r0, domain, channel, chunksize, j, nchunks; showprogress, λ, uniform, fftplan)
         # receive power (within aperture)
@@ -110,39 +111,34 @@ function repeated_wave_propagation(utx::AbstractArray{Complex{T}}, midx, domain:
     nothing
 end
 
-function repeated_wave_propagation(iters; Nx::Int=128, λ=1.55e-6, w0=1e-2, dx1=1e-3, dxn=1e-3, dz=5e1, Nz::Int=21, m=0, n=0, kwargs...)
-    beam = Beam(m, n, λ, w0)
+function repeated_wave_propagation(iters, T::DataType=Float64; Nx::Int=128, λ=1.55e-6, w0=1e-2, dx1=1e-3, dxn=1e-3, dz=5e1, Nz::Int=21, m=0, n=0, kwargs...)
+    beam = Beam{T}(m, n, λ, w0)
 
     domain = if dx1 == dxn
-        Domain(dx1, dz, Nx, Nz)
+        Domain{T}(dx1, dz, Nx, Nz)
     else
-        Domain(dx1, dxn, dz, Nx, Nz)
+        Domain{T}(dx1, dxn, dz, Nx, Nz)
     end
 
     midx = [(m, n)]
+    utx = gaussian_source(domain, w0, λ, midx...)
 
-    utx = maybe_cuda(broadcast(beam.fun, domain.x(1), domain.y(1) |> transpose, 0))
     repeated_wave_propagation(utx, midx, domain, iters; beam.λ, beam.w0, kwargs...)
 end
 
-function repeated_wave_propagation_multi(iters; Nx::Int=128, λ=1.55e-6, w0=1e-2, dx1=1e-3, dxn=1e-3, dz=5e1, Nz::Int=21, ng=3, kwargs...)
+function repeated_wave_propagation_multi(iters, T::DataType=Float64; Nx::Int=128, λ=1.55e-6, w0=1e-2, dx1=1e-3, dxn=1e-3, dz=5e1, Nz::Int=21, ng=3, kwargs...)
     # simulation domain
     domain = if dx1 == dxn
-        Domain(dx1, dz, Nx, Nz)
+        Domain{T}(dx1, dz, Nx, Nz)
     else
-        Domain(dx1, dxn, dz, Nx, Nz)
+        Domain{T}(dx1, dxn, dz, Nx, Nz)
     end
 
     # we consider the modes from the first ng mode groups
     midx = sort!([(m, n) for m ∈ range(0, ng - 1) for n ∈ range(0, ng - m - 1)], by=sum)
-    Nmodes = length(midx)
 
     # set up transmit modes
-    lgfuns = Dict(m => lg_fun_cart(m[1], m[2], w0, λ) for m ∈ midx)
-    utx = maybe_cuda(zeros(ComplexF64, domain.Nx, domain.Ny, Nmodes))
-    for (k, m) ∈ enumerate(midx)
-        broadcast!(lgfuns[m], view(utx, :, :, k), domain.x(1), domain.y(1) |> transpose, 0)
-    end
+    utx = gaussian_source(domain, w0, λ, midx)
 
     repeated_wave_propagation(utx, midx, domain, iters; λ, w0, kwargs...)
 end
\ No newline at end of file

src/optical_parameters.jl

@@ -89,18 +89,14 @@ end
 """
 To properly represent the atmosphere, we have to find optimal r0s for each partial propagation such that the first few moments match.
 """
-function partial_r0_pw(k0, z, Cn²; debug=false)
+partial_r0_pw(k0, z, Cn²; debug=false) = partial_r0_pw!(zero(z), k0, z, Cn²; debug)
+function partial_r0_pw!(r0s::AbstractVector{T}, k0, z::AbstractVector{T}, Cn²; debug=false) where {T}
     Nz = length(z)
 
     # define target values
     r0 = fried_pw(k0, last(z), Cn²)
     σ²χ = sigma_chi_pw(k0, last(z), Cn²)
 
-    # single propagation does not require any optimization
-    if isone(Nz)
-        return [r0]
-    end
-
     # define A-matrix for plane wave, compare Schmidt (9.75)
     A = zeros(2, Nz)
     @. A[1, :] = 1
@@ -129,17 +125,22 @@ function partial_r0_pw(k0, z, Cn²; debug=false)
     X = Optim.minimizer(result)
 
     # Convert optimization result to screen r0s
-    r0s = X .^ (-3 / 5)
+    @. r0s = X^(-3 / 5)
 
     # check result
     bp = A * X
     r0res = bp[1]^(-3 / 5)
     σ²χres = bp[2] * 1.33 * (z[end] / k0)^(5 / 6)
-    
-    abs((r0res - r0) / r0) > 1e-3 && @warn "Total r0: $(r0res). Target r0: $(r0)." 
-    abs((σ²χres - σ²χ) / σ²χ) > 1e-3 && @warn "Total σ²χ: $(σ²χres). Target σ²χ: $(σ²χ)." 
+
+    abs((r0res - r0) / r0) > 1e-3 && @warn "Total r0: $(r0res). Target r0: $(r0)."
+    abs((σ²χres - σ²χ) / σ²χ) > 1e-3 && @warn "Total σ²χ: $(σ²χres). Target σ²χ: $(σ²χ)."
 
     if debug
+        @info "Coefficient matrix:"
+        display(A)
+        @info "Lower limit for r0:"
+        display(x2 .^ (-3 / 5))
+
         # Check resulting r0sw & rytov
         bp = A * X
         r0res = bp[1]^(-3 / 5)
@@ -151,25 +152,20 @@ function partial_r0_pw(k0, z, Cn²; debug=false)
     return r0s
 end
 
-
-function partial_r0_sw(k0, z, Cn²; debug=false)
+partial_r0_sw(k0, z, Cn²; debug=false) = partial_r0_sw!(zero(z), k0, z, Cn²; debug)
+function partial_r0_sw!(r0s::AbstractVector{T}, k0, z::AbstractVector{T}, Cn²; debug=false) where {T}
     Nz = length(z)
 
     # define target values
     r0 = fried_sw(k0, last(z), Cn²)
     σ²χ = sigma_chi_sw(k0, last(z), Cn²)
 
-    # single propagation does not require any optimization
-    if isone(Nz)
-        return [r0]
-    end
-
     # define A-matrix for plane wave, compare Schmidt (9.75)
     A = zeros(2, Nz)
     @. A[1, :] = (z / z[end])^(5 / 3)
     @views @. A[2, 1:end-1] = (1 - z[1:end-1] / z[end])^(5 / 6) * (z[1:end-1] / z[end])^(5 / 6)
     A[2, end] = 0
-    # define target vector for plane wave
+    # define target vector for spherical wave
     b = [r0^(-5 / 3); σ²χ / 1.33 * (k0 / z[end])^(5 / 6)]
 
     # Objective function
@@ -195,17 +191,22 @@ function partial_r0_sw(k0, z, Cn²; debug=false)
     X = Optim.minimizer(result)
 
     # Convert optimization result to screen r0s
-    r0s = X .^ (-3 / 5)
+    @. r0s = X^(-3 / 5)
 
     # check result
     bp = A * X
     r0res = bp[1]^(-3 / 5)
     σ²χres = bp[2] * 1.33 * (z[end] / k0)^(5 / 6)
-    
-    abs((r0res - r0) / r0) > 1e-3 && @warn "Total r0: $(r0res). Target r0: $(r0)." 
-    abs((σ²χres - σ²χ) / σ²χ) > 1e-3 && @warn "Total σ²χ: $(σ²χres). Target σ²χ: $(σ²χ)." 
+
+    abs((r0res - r0) / r0) > 1e-3 && @warn "Total r0: $(r0res). Target r0: $(r0)."
+    abs((σ²χres - σ²χ) / σ²χ) > 1e-3 && @warn "Total σ²χ: $(σ²χres). Target σ²χ: $(σ²χ)."
 
     if debug
+        @info "Coefficient matrix:"
+        display(A)
+        @info "Lower limit for r0:"
+        display(x2 .^ (-3 / 5))
+
         # Check resulting r0sw & rytov
         bp = A * X
         r0res = bp[1]^(-3 / 5)

src/source.jl

+include("helpers.jl")
+include("generate_lg.jl")
+include("newtypes.jl")
+
+function gaussian_source(domain::Domain{T}, w0, λ, midx=(0, 0)) where {T}
+    utx = zeros(Complex{T}, domain.Nx, domain.Ny)
+    broadcast!(lg_fun_cart(midx..., T(w0), T(λ)), utx, domain.x(1), transpose(domain.y(1)), zero(T))
+
+    maybe_cuda(utx)
+end
+function gaussian_source(domain::Domain{T}, w0, λ, midx::AbstractVector=[(0, 0)]) where {T}
+    utx = zeros(Complex{T}, domain.Nx, domain.Ny, length(midx))
+
+    for (k, m) ∈ enumerate(midx)
+        broadcast!(lg_fun_cart(m..., T(w0), T(λ)), view(utx, :, :, k), domain.x(1), transpose(domain.y(1)), zero(T))
+    end
+
+    maybe_cuda(utx)
+end
+"""
+    point_source(domain[, λ, ROI])
+
+Combination of sinc and Gaussian point-source models according to Schmidt 2010 p. 111, i.e., listing 6.9
+"""
+function point_source(domain::Domain{T}, λ=1.55e-6, ROI=0.94domain.Dx[end]) where {T}
+    utx = maybe_cuda(zeros(Complex{T}, domain.Nx, domain.Ny))
+    x = domain.x(1)
+    yt = domain.y(1) |> transpose
+
+    arg = ROI / (λ * domain.Dz)
+
+    @. utx = λ * domain.Dz * cispi((x^2 + yt^2) / (λ * domain.Dz)) * arg^2 * sinc(arg * x) * sinc(arg * yt) * exp(-(arg / 4)^2 * (x^2 + yt^2))
+end
+
+function super_gaussian_source(domain::Domain{T}, w=0.25, N=8) where {T}
+    utx = maybe_cuda(zeros(Complex{T}, domain.Nx, domain.Ny))
+    x = domain.x(1)
+    yt = domain.y(1) |> transpose
+    @. utx = exp(-((x^2 + yt^2) / T(w)^2)^N)
+end
+
+function plane_source(domain::Domain{T}) where {T}
+    maybe_cuda(fill(one(Complex{T}), domain.Nx, domain.Ny))
+end
\ No newline at end of file

src/turbulent_channel.jl

@@ -136,7 +136,7 @@ end
 function turbulent_channel(T::DataType=Float64; domain::Domain=Domain{T}(0.3 / 256, 1e2, 256, 11), channel::Channel=Channel{T}(1e-14, 1e-2, 1e2, andrews), beam::Beam=Beam{T}(0, 0, 1.55e-6, 1e-2), kwargs...)
 
     # generate Gaussian beam
-    u = broadcast(beam.fun, domain.x(1), domain.y(1) |> transpose, zero(T))
+    u = gaussian_source(domain, beam.w0, beam.λ, (beam.l, beam.m))
 
     # estimate proper r0
     r0 = partial_r0_pw(beam.k0, domain.z, channel.Cn²)