feat: improved source gen & optional single-prec wave propagation

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/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²)