-
Notifications
You must be signed in to change notification settings - Fork 0
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
Add refractive phase screen + power law
- Loading branch information
Anna Tartaglia
committed
Aug 6, 2023
1 parent
e7acef8
commit d967678
Showing
10 changed files
with
155 additions
and
24 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,51 @@ | ||
| export PhaseScreenPowerLaw | ||
| @doc raw""" | ||
| $(TYPEDEF) | ||
| # Fields | ||
| $(FIELDS) | ||
| """ | ||
| struct PhaseScreenPowerLaw{N,S<:AbstractScatteringModel, T<:Number} <: AbstractPowerSpectrumModel{2} | ||
| sm::S | ||
| dx::T | ||
| dy::T | ||
| Vx_km_per_s::T | ||
| Vy_km_per_s::T | ||
|
|
||
| function PhaseScreenPowerLaw(sm::S, dx, dy, Vx_km_per_s=0.0::T, Vy_km_per_s=0.0::T) where {S, T} | ||
| new{2,typeof(sm), typeof(dx)}(sm, dx, dy, Vx_km_per_s, Vy_km_per_s) | ||
| end | ||
| end | ||
|
|
||
| @inline fourieranalytic(::PhaseScreenPowerLaw) = IsAnalytic() | ||
|
|
||
| function StationaryRandomFields.power_point(model::PhaseScreenPowerLaw, q...) | ||
| @assert length(q) == 2 | ||
|
|
||
| q_r = √sum(abs2, q) + 1e-12/model.sm.rin | ||
| q_ϕ = atan(q[2], q[1]) | ||
| return model.sm.Qbar * (q_r*model.sm.rin)^(-(model.sm.α + 2.0)) * exp(-(q_r * model.sm.rin)^2) * Pϕ(model.sm, q_ϕ) | ||
| end | ||
|
|
||
| function StationaryRandomFields.amplitude_point(model::PhaseScreenPowerLaw, t_hr, N, FOV, dq, q...) | ||
| x_offset_pixels = (model.Vx_km_per_s * 1.e5) * (t_hr*3600.) / (FOV/float(N)) | ||
| y_offset_pixels = (model.Vy_km_per_s * 1.e5) * (t_hr*3600.) / (FOV/float(N)) | ||
| q = (q[1]/model.dx, q[2]/model.dy) | ||
|
|
||
| return √(power_point(model, dq.*q...)) * exp(2.0*π*1im*(q[1]*x_offset_pixels + q[2]*y_offset_pixels)/float(N)) | ||
| end | ||
|
|
||
| function StationaryRandomFields.amplitude_map(model::PhaseScreenPowerLaw, noisesignal::Union{AbstractNoiseSignal,AbstractContinuousNoiseSignal}, t_hr = 0.) | ||
| N = noisesignal.dims[1] | ||
| FOV = N * model.sm.D | ||
| dq = 2*π/FOV | ||
| gridofq = rfftfreq(noisesignal) | ||
|
|
||
| prod = collect(Iterators.product(gridofq...)) | ||
| amp = zeros(size(prod)...) | ||
| for i in eachindex(prod)[2:end] | ||
| amp[i] = amplitude_point(model, t_hr, N, FOV, dq, prod[i]...) | ||
| end | ||
| return amp | ||
| end |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,48 @@ | ||
| export AbstractPhaseScreen | ||
| export RefractivePhaseScreen | ||
| export phase_screen | ||
| export wrapped_grad | ||
|
|
||
| abstract type AbstractPhaseScreen end | ||
|
|
||
| struct RefractivePhaseScreen{S, T <: Number, N<:AbstractNoiseSignal, P <: AbstractPowerSpectrumModel} <: AbstractPhaseScreen | ||
| sm::S | ||
| dx::T | ||
| dy::T | ||
| signal::N | ||
| Q::P | ||
| function RefractivePhaseScreen(sm::S, Nx::Number, Ny::Number, dx::T, dy::T, Vx_km_per_s=0.0::T, Vy_km_per_s=0.0::T) where {S, T} | ||
| dx = dx*sm.D | ||
| dy = dy*sm.D | ||
| Q = PhaseScreenPowerLaw(sm, dx, dy, Vx_km_per_s, Vy_km_per_s) | ||
| signal = NoiseSignal((Nx, Ny)) | ||
| return new{S, T, NoiseSignal, typeof(Q)}(sm, dx, dy, signal, Q) | ||
| end | ||
| end | ||
|
|
||
| StationaryRandomFields.generate_gaussian_noise(psm::AbstractPhaseScreen; rng = Random.default_rng()) = generate_gaussian_noise(psm.signal; rng = rng) | ||
|
|
||
| function phase_screen(psm::AbstractPhaseScreen, noise_screen=nothing) | ||
| if isnothing(noise_screen) | ||
| noise_screen = generate_gaussian_noise(psm) | ||
| end | ||
| cns = ContinuousNoiseSignal(psm.signal) | ||
| screengen = PSNoiseGenerator(psm.Q, cns) | ||
| return generate_signal_noise(screengen) | ||
| end | ||
|
|
||
| function wrapped_grad(ϕ) | ||
| nx, ny = size(ϕ) | ||
| gradϕ_x = zeros(Float64, nx, ny) | ||
| gradϕ_y = zeros(Float64, nx, ny) | ||
| for i=1:nx, j=1:ny | ||
| i0 = i == 1 ? nx : i - 1 | ||
| j0 = j == 1 ? ny : j - 1 | ||
| i1 = i == nx ? 1 : i + 1 | ||
| j1 = j == ny ? 1 : j + 1 | ||
|
|
||
| @inbounds gradϕ_x[i, j] = 0.5 * (ϕ[i1, j] - ϕ[i0, j]) | ||
| @inbounds gradϕ_y[i, j] = 0.5 * (ϕ[i, j1] - ϕ[i, j0]) | ||
| end | ||
| return (gradϕ_x, gradϕ_y) | ||
| end |