Skip to content
1 change: 1 addition & 0 deletions .gitignore
Original file line number Diff line number Diff line change
Expand Up @@ -4,3 +4,4 @@
deps/deps.jl
*.mem
*.h5
.vscode/settings.json
20 changes: 16 additions & 4 deletions src/LinearOps.jl
Original file line number Diff line number Diff line change
Expand Up @@ -3,23 +3,32 @@ import FFTW
import Luna: Modes, Grid, PhysData
import Luna.PhysData: wlfreq

# limit α so that we do not get overflow in exp(α*dz)
function αlim!(α)
# magic number: this is 1300 dB/m it could probably be bigger
α[α .> 300] .= 300
end

function make_const_linop(grid::Grid.RealGrid, βfun, αfun, frame_vel)
β = similar(grid.ω)
βfun(β, grid.ω, 0)
α = αfun(grid.ω, 0)
α = zeros(length(grid.ω))
α[2:end] .= αfun(grid.ω[2:end], 0)
αlim!(α)
β1 = 1/frame_vel(0)
linop = @. -im*(β-β1*grid.ω) - α/2
linop[1] = 1
linop[1] = 0
return linop
end

function make_const_linop(grid::Grid.EnvGrid, βfun, αfun, frame_vel, β0ref)
β = similar(grid.ω)
βfun(β, grid.ω, 0)
α = αfun(grid.ω, 0)
αlim!(α)
β1 = 1/frame_vel(0)
linop = -im.*(β .- β1.*(grid.ω .- grid.ω0) .- β0ref) .- α./2
linop[.!grid.sidx] .= 1
linop[.!grid.sidx] .= 0
return linop
end

Expand Down Expand Up @@ -62,7 +71,9 @@ function make_const_linop(grid::Grid.RealGrid, modes, λ0; ref_mode=1)
βconst = zero(grid.ω)
βconst[2:end] = Modes.β.(modes[i], grid.ω[2:end])
βconst[1] = 1
α = Modes.α.(modes[i], grid.ω)
α = zeros(length(grid.ω))
α[2:end] .= Modes.α.(modes[i], grid.ω[2:end])
αlim!(α)
linops[:,i] = im.*(-βconst .+ grid.ω./vel) .- α./2
end
linops
Expand All @@ -82,6 +93,7 @@ function make_const_linop(grid::Grid.EnvGrid, modes, λ0; ref_mode=1, thg=false)
βconst[grid.sidx] = Modes.β.(modes[i], grid.ω[grid.sidx])
βconst[.!grid.sidx] .= 1
α = Modes.α.(modes[i], grid.ω)
αlim!(α)
linops[:,i] = -im.*(βconst .- (grid.ω .- grid.ω0)./vel .- βref) .- α./2
end
linops
Expand Down
Loading