From 513745bc7c3f9c5534b821dae744ece7061e1f8c Mon Sep 17 00:00:00 2001 From: CompatHelper Julia Date: Tue, 11 Jul 2023 19:07:00 +0000 Subject: [PATCH 01/19] CompatHelper: add new compat entry for DocStringExtensions at version 0.9, (keep existing compat) --- Project.toml | 1 + 1 file changed, 1 insertion(+) diff --git a/Project.toml b/Project.toml index e8386bf..36cb836 100644 --- a/Project.toml +++ b/Project.toml @@ -10,6 +10,7 @@ NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec" SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b" [compat] +DocStringExtensions = "0.9" julia = "1.9" [extras] From fd910abe4c112b639245d747c5f312322acadfd2 Mon Sep 17 00:00:00 2001 From: CompatHelper Julia Date: Tue, 11 Jul 2023 19:07:09 +0000 Subject: [PATCH 02/19] CompatHelper: add new compat entry for NonlinearSolve at version 1, (keep existing compat) --- Project.toml | 1 + 1 file changed, 1 insertion(+) diff --git a/Project.toml b/Project.toml index e8386bf..d571dce 100644 --- a/Project.toml +++ b/Project.toml @@ -10,6 +10,7 @@ NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec" SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b" [compat] +NonlinearSolve = "1" julia = "1.9" [extras] From e38d4ff883f72789d2ca1236f97a65d02c2a198a Mon Sep 17 00:00:00 2001 From: CompatHelper Julia Date: Tue, 11 Jul 2023 19:07:12 +0000 Subject: [PATCH 03/19] CompatHelper: add new compat entry for SpecialFunctions at version 2, (keep existing compat) --- Project.toml | 1 + 1 file changed, 1 insertion(+) diff --git a/Project.toml b/Project.toml index e8386bf..60d5285 100644 --- a/Project.toml +++ b/Project.toml @@ -10,6 +10,7 @@ NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec" SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b" [compat] +SpecialFunctions = "2" julia = "1.9" [extras] From ed0c2e4073c63584fd9f802677bfbe7dd1ce65a0 Mon Sep 17 00:00:00 2001 From: Anna Tartaglia Date: Tue, 11 Jul 2023 16:19:51 -0400 Subject: [PATCH 04/19] math corrections --- src/scatteringkernels/abstractscatteringkernel.jl | 2 +- src/scatteringkernels/phasestructure.jl | 6 +++--- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/src/scatteringkernels/abstractscatteringkernel.jl b/src/scatteringkernels/abstractscatteringkernel.jl index eb57379..a4c0c23 100644 --- a/src/scatteringkernels/abstractscatteringkernel.jl +++ b/src/scatteringkernels/abstractscatteringkernel.jl @@ -81,6 +81,6 @@ end Maps visibility kernel for a given observing wavelength λ and fourier space coordinates u, v """ function visibility_point(ps::AbstractScatteringKernel, λ::Number, u::Number, v::Number) - b = (u/λ, v/λ)./(1+ps.M) + b = (u*λ, v*λ)./(1+ps.M) return exp(-.5 * phase_structure_point(ps, b...)) end diff --git a/src/scatteringkernels/phasestructure.jl b/src/scatteringkernels/phasestructure.jl index 5c60470..b152e81 100644 --- a/src/scatteringkernels/phasestructure.jl +++ b/src/scatteringkernels/phasestructure.jl @@ -25,13 +25,13 @@ struct PhaseStructure <: AbstractScatteringKernel hyp1 = _₂F₁(1/2., (2+α)/2., 1., -k) hyp2 = _₂F₁(1/2., -α/2., 1., -k) num = 2^(4. - α)/(α^2.)/gamma(α/2.) - Bmaj = num * hyp1 / √(1+k) / (1+ζ) - Bmin = num * hyp2/(1-ζ) + Bmaj = num / hyp1 / √(1+k) / (1+ζ) + Bmin = num / hyp2/(1-ζ) bm = √(8*log(2)/2*π) * λ0/θmaj l1 = (1. +M)*inscale*(2. /(α * Bmaj))^(1/(α-2.)) l2 = 4. ^(1/(2. -α))*Bmaj^(2/(α-2.))*α^(α/(α-2.))/(1+ζ) C = l2/(((bm/l1)^2. +1)^(α/2.) - 1) - return new(α, inscale, θmaj, θmin, ϕpt, λ0, M, A, ζ, k, Bmaj, Bmin, C) + return new(α, inscale, θmaj, θmin, ϕpt, λ0, M, ζ, A, k, Bmaj, Bmin, C) end end From 71b8b5016a68281a29a7cd394f6e20e4c792b220 Mon Sep 17 00:00:00 2001 From: Anna Tartaglia Date: Tue, 11 Jul 2023 16:24:59 -0400 Subject: [PATCH 05/19] testing math --- src/scatteringkernels/phasestructure.jl | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/src/scatteringkernels/phasestructure.jl b/src/scatteringkernels/phasestructure.jl index b152e81..29cc289 100644 --- a/src/scatteringkernels/phasestructure.jl +++ b/src/scatteringkernels/phasestructure.jl @@ -28,10 +28,7 @@ struct PhaseStructure <: AbstractScatteringKernel Bmaj = num / hyp1 / √(1+k) / (1+ζ) Bmin = num / hyp2/(1-ζ) - bm = √(8*log(2)/2*π) * λ0/θmaj - l1 = (1. +M)*inscale*(2. /(α * Bmaj))^(1/(α-2.)) - l2 = 4. ^(1/(2. -α))*Bmaj^(2/(α-2.))*α^(α/(α-2.))/(1+ζ) - C = l2/(((bm/l1)^2. +1)^(α/2.) - 1) + C = ((1+M)*ps.inscale / (2*π*ps.λ0 ))^2. * (ps.θmaj^2 + ps.θmin^2)/(8 * log(2.)) return new(α, inscale, θmaj, θmin, ϕpt, λ0, M, ζ, A, k, Bmaj, Bmin, C) end end From 1d8afbd50ed605c2d76acb030119447ed32ef068 Mon Sep 17 00:00:00 2001 From: Anna Tartaglia Date: Tue, 11 Jul 2023 16:26:43 -0400 Subject: [PATCH 06/19] math testing --- src/scatteringkernels/phasestructure.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/scatteringkernels/phasestructure.jl b/src/scatteringkernels/phasestructure.jl index 29cc289..d6191ca 100644 --- a/src/scatteringkernels/phasestructure.jl +++ b/src/scatteringkernels/phasestructure.jl @@ -28,7 +28,7 @@ struct PhaseStructure <: AbstractScatteringKernel Bmaj = num / hyp1 / √(1+k) / (1+ζ) Bmin = num / hyp2/(1-ζ) - C = ((1+M)*ps.inscale / (2*π*ps.λ0 ))^2. * (ps.θmaj^2 + ps.θmin^2)/(8 * log(2.)) + C = ((1+M)*inscale / (2*π*λ0 ))^2. * (θmaj^2 + θmin^2)/(8 * log(2.)) return new(α, inscale, θmaj, θmin, ϕpt, λ0, M, ζ, A, k, Bmaj, Bmin, C) end end From 155542a6228f080c96aabe60f4cfe7be8fc09f2e Mon Sep 17 00:00:00 2001 From: Anna Tartaglia Date: Wed, 12 Jul 2023 10:30:18 -0400 Subject: [PATCH 07/19] fix C calc --- src/scatteringkernels/phasestructure.jl | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/src/scatteringkernels/phasestructure.jl b/src/scatteringkernels/phasestructure.jl index d6191ca..bba884c 100644 --- a/src/scatteringkernels/phasestructure.jl +++ b/src/scatteringkernels/phasestructure.jl @@ -18,17 +18,21 @@ struct PhaseStructure <: AbstractScatteringKernel Bmin::Number C::Number function PhaseStructure(α, inscale, θmaj, θmin, ϕpt, λ0, M) + # solve for k from dipole model in Psaltis 2018 A = θmaj/θmin ζ = (A^2. - 1)/(A^2. +1) k = findk(KFinder(α, A)) + # B parameters hyp1 = _₂F₁(1/2., (2+α)/2., 1., -k) hyp2 = _₂F₁(1/2., -α/2., 1., -k) num = 2^(4. - α)/(α^2.)/gamma(α/2.) Bmaj = num / hyp1 / √(1+k) / (1+ζ) Bmin = num / hyp2/(1-ζ) - - C = ((1+M)*inscale / (2*π*λ0 ))^2. * (θmaj^2 + θmin^2)/(8 * log(2.)) + + # C parameter + Qbar = 2.0/gamma((2-α)/2.) * (ps.inscale^2*(1.0 + M)/((2.0 * log(2.0))^0.5/π*(ps.λ0/(2.0*π))^2) )^2 * ( (ps.θmaj^2 + ps.θmin^2)) + C = (ps.λ0/(2.0*π))^2 * Qbar*gamma(1.0 - α/2.0)/(8.0*π^2*ps.inscale^2) return new(α, inscale, θmaj, θmin, ϕpt, λ0, M, ζ, A, k, Bmaj, Bmin, C) end end From 762fc4ac7ab3b75a1edfe45c34762d37cc948476 Mon Sep 17 00:00:00 2001 From: Anna Tartaglia Date: Wed, 12 Jul 2023 10:35:39 -0400 Subject: [PATCH 08/19] phasestructure --- src/scatteringkernels/phasestructure.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/scatteringkernels/phasestructure.jl b/src/scatteringkernels/phasestructure.jl index bba884c..b5e2239 100644 --- a/src/scatteringkernels/phasestructure.jl +++ b/src/scatteringkernels/phasestructure.jl @@ -31,8 +31,8 @@ struct PhaseStructure <: AbstractScatteringKernel Bmin = num / hyp2/(1-ζ) # C parameter - Qbar = 2.0/gamma((2-α)/2.) * (ps.inscale^2*(1.0 + M)/((2.0 * log(2.0))^0.5/π*(ps.λ0/(2.0*π))^2) )^2 * ( (ps.θmaj^2 + ps.θmin^2)) - C = (ps.λ0/(2.0*π))^2 * Qbar*gamma(1.0 - α/2.0)/(8.0*π^2*ps.inscale^2) + Qbar = 2.0/gamma((2-α)/2.) * (inscale^2*(1.0 + M)/((2.0 * log(2.0))^0.5/π*(λ0/(2.0*π))^2) )^2 * ( (θmaj^2 + θmin^2)) + C = (λ0/(2.0*π))^2 * Qbar*gamma(1.0 - α/2.0)/(8.0*π^2*inscale^2) return new(α, inscale, θmaj, θmin, ϕpt, λ0, M, ζ, A, k, Bmaj, Bmin, C) end end From 6caefb6adb10abfb2a85af26e6b114fed1c9aff6 Mon Sep 17 00:00:00 2001 From: Anna Tartaglia Date: Wed, 12 Jul 2023 13:29:59 -0400 Subject: [PATCH 09/19] math corrections in scattering kernel --- .../abstractscatteringkernel.jl | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/src/scatteringkernels/abstractscatteringkernel.jl b/src/scatteringkernels/abstractscatteringkernel.jl index a4c0c23..a1080ba 100644 --- a/src/scatteringkernels/abstractscatteringkernel.jl +++ b/src/scatteringkernels/abstractscatteringkernel.jl @@ -27,7 +27,7 @@ abstract type AbstractScatteringKernel end Maps D_maj(r) for given r Equation 33 of Psaltis et al. 2018 """ -@inline function calc_dmaj(ps::AbstractScatteringKernel, r::Number) +@inline function calc_dmaj(ps::AbstractScatteringKernel, λ::Number, r::Number) Bmaj = ps.Bmaj C = ps.C α = ps.α @@ -36,7 +36,7 @@ Equation 33 of Psaltis et al. 2018 d1 = C*(1. +ζ)/2. * Bmaj*(2. /(α*Bmaj))^(-α/(2. -α)) d2 = (2. /(α*Bmaj))^(2. /(2. -α)) - return d1*((1+d2*(r/rin)^2.)^(α/2.)-1) + return d1*((1+d2*(r/rin)^2.)^(α/2.)-1) * (λ/ps.λ0)^2. end """ @@ -45,7 +45,7 @@ end Maps D_min(r) for given r Equation 34 of Psaltis et al. 2018 """ -@inline function calc_dmin(ps::AbstractScatteringKernel, r::Number) +@inline function calc_dmin(ps::AbstractScatteringKernel, λ::Number, r::Number) Bmin = ps.Bmin C = ps.C α = ps.α @@ -54,7 +54,7 @@ Equation 34 of Psaltis et al. 2018 d1 = C*(1. -ζ)/2. * Bmin*(2. /(α*Bmin))^(-α/(2. -α)) d2 = (2. /(α*Bmin))^(2. /(2. -α)) - return d1*((1+d2*(r/rin)^2.)^(α/2.)-1) + return d1*((1+d2*(r/rin)^2.)^(α/2.)-1) * (λ/ps.λ0)^2. end """ @@ -63,12 +63,12 @@ end Maps phase structure function D_ϕ(r, ϕ), first converting x, y into polar coordinates Equation 35 of Psaltis et al. 2018 """ -function phase_structure_point(ps::AbstractScatteringKernel, x::Number, y::Number) +function phase_structure_point(ps::AbstractScatteringKernel, λ::Number, x::Number, y::Number) ϕ0 = 0 r = √(x^2. + y^2.) ϕ = atan(y,x) - dmaj = calc_dmaj(ps, r) - dmin = calc_dmin(ps, r) + dmaj = calc_dmaj(ps, λ, r) + dmin = calc_dmin(ps, λ, r) add = (dmaj + dmin)/2. sub = (dmaj - dmin)/2. return add + sub * cos(2. *(ϕ-ϕ0)) @@ -81,6 +81,6 @@ end Maps visibility kernel for a given observing wavelength λ and fourier space coordinates u, v """ function visibility_point(ps::AbstractScatteringKernel, λ::Number, u::Number, v::Number) - b = (u*λ, v*λ)./(1+ps.M) - return exp(-.5 * phase_structure_point(ps, b...)) + b = (u, v).* (λ/(1+ps.M)) + return exp(-.5 * phase_structure_point(ps, λ, b...)) end From 8d2b0eec0fdc388bcece968bbb93d4d0edee365f Mon Sep 17 00:00:00 2001 From: Anna Tartaglia Date: Wed, 12 Jul 2023 16:12:44 -0400 Subject: [PATCH 10/19] fixing scattering kernel equations --- src/scatteringkernels/abstractscatteringkernel.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/scatteringkernels/abstractscatteringkernel.jl b/src/scatteringkernels/abstractscatteringkernel.jl index a1080ba..9463b6d 100644 --- a/src/scatteringkernels/abstractscatteringkernel.jl +++ b/src/scatteringkernels/abstractscatteringkernel.jl @@ -81,6 +81,6 @@ end Maps visibility kernel for a given observing wavelength λ and fourier space coordinates u, v """ function visibility_point(ps::AbstractScatteringKernel, λ::Number, u::Number, v::Number) - b = (u, v).* (λ/(1+ps.M)) + b = (u*λ, v*λ)./(1+ps.M) return exp(-.5 * phase_structure_point(ps, λ, b...)) -end +end \ No newline at end of file From 4eff68a17b69825a7cbba71637a96fdc27d9b11b Mon Sep 17 00:00:00 2001 From: Anna Tartaglia Date: Thu, 13 Jul 2023 13:58:16 -0400 Subject: [PATCH 11/19] unit conversions --- Manifest.toml | 5 +++++ Project.toml | 1 + src/ScatteringOptics.jl | 1 + src/scatteringkernels/phasestructure.jl | 6 +++++- 4 files changed, 12 insertions(+), 1 deletion(-) diff --git a/Manifest.toml b/Manifest.toml index 13b01d9..85bb4c5 100644 --- a/Manifest.toml +++ b/Manifest.toml @@ -56,6 +56,11 @@ version = "0.1.29" [[Artifacts]] uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33" +[[AstroAngles]] +git-tree-sha1 = "41621fa5ed5f7614b75eea8e0b3cfd967b284c87" +uuid = "5c4adb95-c1fc-4c53-b4ea-2a94080c53d2" +version = "0.1.3" + [[Base64]] uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f" diff --git a/Project.toml b/Project.toml index e8386bf..f0acb63 100644 --- a/Project.toml +++ b/Project.toml @@ -4,6 +4,7 @@ authors = ["Anna Tartaglia", "Kazunori Akiyama"] version = "1.0.0-DEV" [deps] +AstroAngles = "5c4adb95-c1fc-4c53-b4ea-2a94080c53d2" DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" HypergeometricFunctions = "34004b35-14d8-5ef3-9330-4cdb6864b03a" NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec" diff --git a/src/ScatteringOptics.jl b/src/ScatteringOptics.jl index f1043bc..f075046 100644 --- a/src/ScatteringOptics.jl +++ b/src/ScatteringOptics.jl @@ -5,6 +5,7 @@ using DocStringExtensions using HypergeometricFunctions using SpecialFunctions using NonlinearSolve +using AstroAngles #k finders include("./kfinders/abstractkfinder.jl") diff --git a/src/scatteringkernels/phasestructure.jl b/src/scatteringkernels/phasestructure.jl index b5e2239..77df0f7 100644 --- a/src/scatteringkernels/phasestructure.jl +++ b/src/scatteringkernels/phasestructure.jl @@ -23,6 +23,10 @@ struct PhaseStructure <: AbstractScatteringKernel ζ = (A^2. - 1)/(A^2. +1) k = findk(KFinder(α, A)) + #milliarcseconds to radians + θmaj_rad = dms2rad((0,0,θmaj*10^-3)) + θmin_rad = dms2rad((0,0,θmin*10^-3)) + # B parameters hyp1 = _₂F₁(1/2., (2+α)/2., 1., -k) hyp2 = _₂F₁(1/2., -α/2., 1., -k) @@ -31,7 +35,7 @@ struct PhaseStructure <: AbstractScatteringKernel Bmin = num / hyp2/(1-ζ) # C parameter - Qbar = 2.0/gamma((2-α)/2.) * (inscale^2*(1.0 + M)/((2.0 * log(2.0))^0.5/π*(λ0/(2.0*π))^2) )^2 * ( (θmaj^2 + θmin^2)) + Qbar = 2.0/gamma((2-α)/2.) * (inscale^2*(1.0 + M)/((2.0 * log(2.0))^0.5/π*(λ0/(2.0*π))^2) )^2 * ( (θmaj_rad^2 + θmin_rad^2)) C = (λ0/(2.0*π))^2 * Qbar*gamma(1.0 - α/2.0)/(8.0*π^2*inscale^2) return new(α, inscale, θmaj, θmin, ϕpt, λ0, M, ζ, A, k, Bmaj, Bmin, C) end From e276a4166dffee63685400863606b7b2836dd225 Mon Sep 17 00:00:00 2001 From: Anna Tartaglia Date: Thu, 13 Jul 2023 15:17:50 -0400 Subject: [PATCH 12/19] fix unit conversions, add Sag A* defaults for ps --- .../abstractscatteringkernel.jl | 20 +++++++++---------- src/scatteringkernels/phasestructure.jl | 2 +- 2 files changed, 11 insertions(+), 11 deletions(-) diff --git a/src/scatteringkernels/abstractscatteringkernel.jl b/src/scatteringkernels/abstractscatteringkernel.jl index 9463b6d..712d372 100644 --- a/src/scatteringkernels/abstractscatteringkernel.jl +++ b/src/scatteringkernels/abstractscatteringkernel.jl @@ -27,7 +27,7 @@ abstract type AbstractScatteringKernel end Maps D_maj(r) for given r Equation 33 of Psaltis et al. 2018 """ -@inline function calc_dmaj(ps::AbstractScatteringKernel, λ::Number, r::Number) +@inline function calc_dmaj(ps::AbstractScatteringKernel, r::Number) Bmaj = ps.Bmaj C = ps.C α = ps.α @@ -36,7 +36,7 @@ Equation 33 of Psaltis et al. 2018 d1 = C*(1. +ζ)/2. * Bmaj*(2. /(α*Bmaj))^(-α/(2. -α)) d2 = (2. /(α*Bmaj))^(2. /(2. -α)) - return d1*((1+d2*(r/rin)^2.)^(α/2.)-1) * (λ/ps.λ0)^2. + return d1*((1+d2*(r/rin)^2.)^(α/2.)-1) end """ @@ -45,7 +45,7 @@ end Maps D_min(r) for given r Equation 34 of Psaltis et al. 2018 """ -@inline function calc_dmin(ps::AbstractScatteringKernel, λ::Number, r::Number) +@inline function calc_dmin(ps::AbstractScatteringKernel, r::Number) Bmin = ps.Bmin C = ps.C α = ps.α @@ -54,7 +54,7 @@ Equation 34 of Psaltis et al. 2018 d1 = C*(1. -ζ)/2. * Bmin*(2. /(α*Bmin))^(-α/(2. -α)) d2 = (2. /(α*Bmin))^(2. /(2. -α)) - return d1*((1+d2*(r/rin)^2.)^(α/2.)-1) * (λ/ps.λ0)^2. + return d1*((1+d2*(r/rin)^2.)^(α/2.)-1) end """ @@ -63,12 +63,12 @@ end Maps phase structure function D_ϕ(r, ϕ), first converting x, y into polar coordinates Equation 35 of Psaltis et al. 2018 """ -function phase_structure_point(ps::AbstractScatteringKernel, λ::Number, x::Number, y::Number) - ϕ0 = 0 +function phase_structure_point(ps::AbstractScatteringKernel, x::Number, y::Number) + ϕ0 = deg2rad(90 - ps.ϕpt) r = √(x^2. + y^2.) ϕ = atan(y,x) - dmaj = calc_dmaj(ps, λ, r) - dmin = calc_dmin(ps, λ, r) + dmaj = calc_dmaj(ps, r) + dmin = calc_dmin(ps, r) add = (dmaj + dmin)/2. sub = (dmaj - dmin)/2. return add + sub * cos(2. *(ϕ-ϕ0)) @@ -80,7 +80,7 @@ end Maps visibility kernel for a given observing wavelength λ and fourier space coordinates u, v """ -function visibility_point(ps::AbstractScatteringKernel, λ::Number, u::Number, v::Number) - b = (u*λ, v*λ)./(1+ps.M) +function visibility_point(ps::AbstractScatteringKernel, u::Number, v::Number) + b = (u, v).* (λ/(1+ps.M)) return exp(-.5 * phase_structure_point(ps, λ, b...)) end \ No newline at end of file diff --git a/src/scatteringkernels/phasestructure.jl b/src/scatteringkernels/phasestructure.jl index 77df0f7..678f5c5 100644 --- a/src/scatteringkernels/phasestructure.jl +++ b/src/scatteringkernels/phasestructure.jl @@ -17,7 +17,7 @@ struct PhaseStructure <: AbstractScatteringKernel Bmaj::Number Bmin::Number C::Number - function PhaseStructure(α, inscale, θmaj, θmin, ϕpt, λ0, M) + function PhaseStructure(α=1.38, inscale=800e5, θmaj=1.380, θmin=.703, ϕpt=81.9, λ0=1, M=.53) # solve for k from dipole model in Psaltis 2018 A = θmaj/θmin ζ = (A^2. - 1)/(A^2. +1) From 896306b25fef35ea10809efff8bbe6c4ae112f19 Mon Sep 17 00:00:00 2001 From: Anna Tartaglia Date: Thu, 13 Jul 2023 15:20:36 -0400 Subject: [PATCH 13/19] fix units, add default ps values --- src/scatteringkernels/abstractscatteringkernel.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/scatteringkernels/abstractscatteringkernel.jl b/src/scatteringkernels/abstractscatteringkernel.jl index 712d372..b40f13b 100644 --- a/src/scatteringkernels/abstractscatteringkernel.jl +++ b/src/scatteringkernels/abstractscatteringkernel.jl @@ -80,7 +80,7 @@ end Maps visibility kernel for a given observing wavelength λ and fourier space coordinates u, v """ -function visibility_point(ps::AbstractScatteringKernel, u::Number, v::Number) +function visibility_point(ps::AbstractScatteringKernel, λ::Number, u::Number, v::Number) b = (u, v).* (λ/(1+ps.M)) return exp(-.5 * phase_structure_point(ps, λ, b...)) end \ No newline at end of file From db08e566f9ee9063419e688a82bda07ee1479af1 Mon Sep 17 00:00:00 2001 From: Anna Tartaglia Date: Fri, 14 Jul 2023 11:25:16 -0400 Subject: [PATCH 14/19] more math modifications --- src/scatteringkernels/abstractscatteringkernel.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/scatteringkernels/abstractscatteringkernel.jl b/src/scatteringkernels/abstractscatteringkernel.jl index b40f13b..4969d12 100644 --- a/src/scatteringkernels/abstractscatteringkernel.jl +++ b/src/scatteringkernels/abstractscatteringkernel.jl @@ -82,5 +82,5 @@ Maps visibility kernel for a given observing wavelength λ and fourier space coo """ function visibility_point(ps::AbstractScatteringKernel, λ::Number, u::Number, v::Number) b = (u, v).* (λ/(1+ps.M)) - return exp(-.5 * phase_structure_point(ps, λ, b...)) + return exp(-.5 * phase_structure_point(ps, b...)) end \ No newline at end of file From 13f16efb3d29e83e8f04961b7ed68279adf939d7 Mon Sep 17 00:00:00 2001 From: Anna Tartaglia Date: Mon, 17 Jul 2023 16:43:02 -0400 Subject: [PATCH 15/19] dphi exact functions --- Manifest.toml | 6 +++ Project.toml | 1 + src/ScatteringOptics.jl | 1 + .../abstractscatteringkernel.jl | 47 ++++++++++++++----- src/scatteringkernels/phasestructure.jl | 5 +- 5 files changed, 48 insertions(+), 12 deletions(-) diff --git a/Manifest.toml b/Manifest.toml index 85bb4c5..dd6a9e1 100644 --- a/Manifest.toml +++ b/Manifest.toml @@ -567,6 +567,12 @@ version = "1.4.0" deps = ["Unicode"] uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7" +[[QuadGK]] +deps = ["DataStructures", "LinearAlgebra"] +git-tree-sha1 = "6ec7ac8412e83d57e313393220879ede1740f9ee" +uuid = "1fd47b50-473d-5c70-9696-f719f8f3bcdc" +version = "2.8.2" + [[REPL]] deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"] uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb" diff --git a/Project.toml b/Project.toml index f0acb63..7cf3fcc 100644 --- a/Project.toml +++ b/Project.toml @@ -8,6 +8,7 @@ AstroAngles = "5c4adb95-c1fc-4c53-b4ea-2a94080c53d2" DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" HypergeometricFunctions = "34004b35-14d8-5ef3-9330-4cdb6864b03a" NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec" +QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc" SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b" [compat] diff --git a/src/ScatteringOptics.jl b/src/ScatteringOptics.jl index f075046..94b8989 100644 --- a/src/ScatteringOptics.jl +++ b/src/ScatteringOptics.jl @@ -6,6 +6,7 @@ using HypergeometricFunctions using SpecialFunctions using NonlinearSolve using AstroAngles +using QuadGK #k finders include("./kfinders/abstractkfinder.jl") diff --git a/src/scatteringkernels/abstractscatteringkernel.jl b/src/scatteringkernels/abstractscatteringkernel.jl index 4969d12..2424b97 100644 --- a/src/scatteringkernels/abstractscatteringkernel.jl +++ b/src/scatteringkernels/abstractscatteringkernel.jl @@ -27,7 +27,7 @@ abstract type AbstractScatteringKernel end Maps D_maj(r) for given r Equation 33 of Psaltis et al. 2018 """ -@inline function calc_dmaj(ps::AbstractScatteringKernel, r::Number) +@inline function calc_dmaj(ps::AbstractScatteringKernel, λ::Number, r::Number) Bmaj = ps.Bmaj C = ps.C α = ps.α @@ -36,7 +36,7 @@ Equation 33 of Psaltis et al. 2018 d1 = C*(1. +ζ)/2. * Bmaj*(2. /(α*Bmaj))^(-α/(2. -α)) d2 = (2. /(α*Bmaj))^(2. /(2. -α)) - return d1*((1+d2*(r/rin)^2.)^(α/2.)-1) + return d1*((1+d2*(r/rin)^2.)^(α/2.)-1) * (λ/ps.λ0)^2. end """ @@ -45,7 +45,7 @@ end Maps D_min(r) for given r Equation 34 of Psaltis et al. 2018 """ -@inline function calc_dmin(ps::AbstractScatteringKernel, r::Number) +@inline function calc_dmin(ps::AbstractScatteringKernel, λ::Number, r::Number) Bmin = ps.Bmin C = ps.C α = ps.α @@ -54,33 +54,58 @@ Equation 34 of Psaltis et al. 2018 d1 = C*(1. -ζ)/2. * Bmin*(2. /(α*Bmin))^(-α/(2. -α)) d2 = (2. /(α*Bmin))^(2. /(2. -α)) - return d1*((1+d2*(r/rin)^2.)^(α/2.)-1) + return d1*((1+d2*(r/rin)^2.)^(α/2.)-1) * (λ/ps.λ0)^2. end """ - phase_structure_point(ps::AbstractScatteringKernel, x::Number, y::Number) + Dϕ_approx(ps::AbstractScatteringKernel, λ::Number, x::Number, y::Number) -Maps phase structure function D_ϕ(r, ϕ), first converting x, y into polar coordinates +Maps approximate phase structure function D_ϕ(r, ϕ) at observing wavelength λ, first converting +x and y into polar coordinates Equation 35 of Psaltis et al. 2018 """ -function phase_structure_point(ps::AbstractScatteringKernel, x::Number, y::Number) +function Dϕ_approx(ps::AbstractScatteringKernel, λ::Number, x::Number, y::Number) ϕ0 = deg2rad(90 - ps.ϕpt) r = √(x^2. + y^2.) ϕ = atan(y,x) - dmaj = calc_dmaj(ps, r) - dmin = calc_dmin(ps, r) + dmaj = calc_dmaj(ps, λ, r) + dmin = calc_dmin(ps, λ, r) add = (dmaj + dmin)/2. sub = (dmaj - dmin)/2. return add + sub * cos(2. *(ϕ-ϕ0)) end +P_ϕ(ps::AbstractScatteringKernel, ϕ) = ps.P_ϕ0*(1. +ps.k*sin(ϕ-deg2rad(90-ps.ϕpt))^2.)^(-(ps.α+2.)/2.) + +function dDϕ_dz(ps::AbstractScatteringKernel, λ::Number, r::Number, ϕ::Number, ϕ_q) +"""differential contribution to the phase structure function +""" + return 4.0 * (λ/ps.λ0)^2. * ps.C/ps.α * (_₁F₁(-ps.α/2.0, 0.5, -r^2. /(4. *ps.inscale^2.)*cos(ϕ - ϕ_q)^2) - 1.) +end + +""" + Dϕ_exact(ps::AbstractScatteringKernel, λ::Number, x::Number, y::Number) + +Maps exact phase structure function D_ϕ(r, ϕ)at observing wavelength λ, first converting +x and y into polar coordinates +""" +function Dϕ_exact(ps::AbstractScatteringKernel, λ::Number, x::Number, y::Number) + r = (x^2 + y^2)^0.5 + ϕ = atan(y, x) + f(ϕ_q) = dDϕ_dz(ps, λ, r, ϕ, ϕ_q)*P_ϕ(ps, ϕ_q) + return quadgk(f, 0, 2.0*π)[1] +end """ visibility_point(ps::AbstractScatteringKernel, λ::Number, u::Number, v::Number) Maps visibility kernel for a given observing wavelength λ and fourier space coordinates u, v """ -function visibility_point(ps::AbstractScatteringKernel, λ::Number, u::Number, v::Number) +function visibility_point(ps::AbstractScatteringKernel, λ::Number, u::Number, v::Number, use_approximate_form=true) b = (u, v).* (λ/(1+ps.M)) - return exp(-.5 * phase_structure_point(ps, b...)) + if use_approximate_form == true + return exp(-.5 * Dϕ_approx(ps, λ, b...)) + else + return exp(-.5 * Dϕ_exact(ps, λ, b...)) + end end \ No newline at end of file diff --git a/src/scatteringkernels/phasestructure.jl b/src/scatteringkernels/phasestructure.jl index 678f5c5..665a092 100644 --- a/src/scatteringkernels/phasestructure.jl +++ b/src/scatteringkernels/phasestructure.jl @@ -17,6 +17,7 @@ struct PhaseStructure <: AbstractScatteringKernel Bmaj::Number Bmin::Number C::Number + P_ϕ0::Number function PhaseStructure(α=1.38, inscale=800e5, θmaj=1.380, θmin=.703, ϕpt=81.9, λ0=1, M=.53) # solve for k from dipole model in Psaltis 2018 A = θmaj/θmin @@ -37,6 +38,8 @@ struct PhaseStructure <: AbstractScatteringKernel # C parameter Qbar = 2.0/gamma((2-α)/2.) * (inscale^2*(1.0 + M)/((2.0 * log(2.0))^0.5/π*(λ0/(2.0*π))^2) )^2 * ( (θmaj_rad^2 + θmin_rad^2)) C = (λ0/(2.0*π))^2 * Qbar*gamma(1.0 - α/2.0)/(8.0*π^2*inscale^2) - return new(α, inscale, θmaj, θmin, ϕpt, λ0, M, ζ, A, k, Bmaj, Bmin, C) + + P_ϕ0 = 1.0/(2.0*π* _₂F₁((α + 2.0)/2.0, 0.5, 1.0, -k)) + return new(α, inscale, θmaj, θmin, ϕpt, λ0, M, ζ, A, k, Bmaj, Bmin, C, P_ϕ0) end end From 653aacb331adfe98bbddfff628252177cb0fee26 Mon Sep 17 00:00:00 2001 From: Anna Tartaglia Date: Thu, 20 Jul 2023 10:56:24 -0400 Subject: [PATCH 16/19] Change equations for consistency --- .../abstractscatteringkernel.jl | 24 ++++++------- src/scatteringkernels/phasestructure.jl | 36 ++++++++++++------- 2 files changed, 36 insertions(+), 24 deletions(-) diff --git a/src/scatteringkernels/abstractscatteringkernel.jl b/src/scatteringkernels/abstractscatteringkernel.jl index 2424b97..ca5b31d 100644 --- a/src/scatteringkernels/abstractscatteringkernel.jl +++ b/src/scatteringkernels/abstractscatteringkernel.jl @@ -1,9 +1,13 @@ export AbstractScatteringKernel export calc_dmaj export calc_dmin -export phase_structure_point +export Dϕ_approx +export P_ϕ +export dDϕ_dz +export Dϕ_exact export visibility_point + """ AbstractScatteringKernel @@ -29,13 +33,12 @@ Equation 33 of Psaltis et al. 2018 """ @inline function calc_dmaj(ps::AbstractScatteringKernel, λ::Number, r::Number) Bmaj = ps.Bmaj - C = ps.C + Amaj = ps.Amaj α = ps.α - ζ = ps.ζ rin = ps.inscale - d1 = C*(1. +ζ)/2. * Bmaj*(2. /(α*Bmaj))^(-α/(2. -α)) - d2 = (2. /(α*Bmaj))^(2. /(2. -α)) + d1 = Bmaj*(2. *Amaj/(α*Bmaj))^(-α/(2. -α)) + d2 = (2. *Amaj/(α*Bmaj))^(2. /(2. -α)) return d1*((1+d2*(r/rin)^2.)^(α/2.)-1) * (λ/ps.λ0)^2. end @@ -47,13 +50,12 @@ Equation 34 of Psaltis et al. 2018 """ @inline function calc_dmin(ps::AbstractScatteringKernel, λ::Number, r::Number) Bmin = ps.Bmin - C = ps.C + Amin = ps.Amin α = ps.α - ζ = ps.ζ rin = ps.inscale - d1 = C*(1. -ζ)/2. * Bmin*(2. /(α*Bmin))^(-α/(2. -α)) - d2 = (2. /(α*Bmin))^(2. /(2. -α)) + d1 = Bmin*(2. *Amin/(α*Bmin))^(-α/(2. -α)) + d2 = (2. *Amin/(α*Bmin))^(2. /(2. -α)) return d1*((1+d2*(r/rin)^2.)^(α/2.)-1) * (λ/ps.λ0)^2. end @@ -65,17 +67,15 @@ x and y into polar coordinates Equation 35 of Psaltis et al. 2018 """ function Dϕ_approx(ps::AbstractScatteringKernel, λ::Number, x::Number, y::Number) - ϕ0 = deg2rad(90 - ps.ϕpt) r = √(x^2. + y^2.) ϕ = atan(y,x) dmaj = calc_dmaj(ps, λ, r) dmin = calc_dmin(ps, λ, r) add = (dmaj + dmin)/2. sub = (dmaj - dmin)/2. - return add + sub * cos(2. *(ϕ-ϕ0)) + return add + sub * cos(2. *(ϕ-ps.ϕ0)) end -P_ϕ(ps::AbstractScatteringKernel, ϕ) = ps.P_ϕ0*(1. +ps.k*sin(ϕ-deg2rad(90-ps.ϕpt))^2.)^(-(ps.α+2.)/2.) function dDϕ_dz(ps::AbstractScatteringKernel, λ::Number, r::Number, ϕ::Number, ϕ_q) """differential contribution to the phase structure function diff --git a/src/scatteringkernels/phasestructure.jl b/src/scatteringkernels/phasestructure.jl index 665a092..aa20fbb 100644 --- a/src/scatteringkernels/phasestructure.jl +++ b/src/scatteringkernels/phasestructure.jl @@ -3,8 +3,8 @@ export PhaseStructure $(TYPEDEF) """ -struct PhaseStructure <: AbstractScatteringKernel - α::Number +struct PhaseStructure{T<:Number, F<:Function} <: AbstractScatteringKernel + α::T inscale::Number θmaj::Number θmin::Number @@ -16,9 +16,13 @@ struct PhaseStructure <: AbstractScatteringKernel k::Number Bmaj::Number Bmin::Number + Qbar::Number C::Number - P_ϕ0::Number - function PhaseStructure(α=1.38, inscale=800e5, θmaj=1.380, θmin=.703, ϕpt=81.9, λ0=1, M=.53) + Amaj::Number + Amin::Number + ϕ0::Number + P_ϕ::F + function PhaseStructure(α=1.38, inscale=800e5, θmaj=1.380, θmin=.703, ϕpt=81.9, λ0=1., M=.53) # solve for k from dipole model in Psaltis 2018 A = θmaj/θmin ζ = (A^2. - 1)/(A^2. +1) @@ -28,18 +32,26 @@ struct PhaseStructure <: AbstractScatteringKernel θmaj_rad = dms2rad((0,0,θmaj*10^-3)) θmin_rad = dms2rad((0,0,θmin*10^-3)) + ϕ0 = deg2rad(90-ϕpt) + P_ϕ0 = 1.0/(2.0*π* _₂F₁((α + 2.0)/2.0, 0.5, 1.0, -k)) + P_ϕ(ϕ) = P_ϕ0*(1. +k*sin(ϕ-ϕ0)^2.)^(-(α+2.)/2.) + # B parameters - hyp1 = _₂F₁(1/2., (2+α)/2., 1., -k) - hyp2 = _₂F₁(1/2., -α/2., 1., -k) - num = 2^(4. - α)/(α^2.)/gamma(α/2.) - Bmaj = num / hyp1 / √(1+k) / (1+ζ) - Bmin = num / hyp2/(1-ζ) + fmaj(ϕ) = abs(cos(ϕ0-ϕ))^α*P_ϕ(ϕ) + fmin(ϕ) = abs(sin(ϕ0-ϕ))^α*P_ϕ(ϕ) + B_prefac = C*2.0^(2.0-α)*π^0.5/(α*gamma((α+1.0)/2.0)) + Bmaj = B_prefac*quadgk(fmaj, 0, 2*π)[1] + Bmin = B_prefac*quadgk(fmin, 0, 2*π)[1] # C parameter - Qbar = 2.0/gamma((2-α)/2.) * (inscale^2*(1.0 + M)/((2.0 * log(2.0))^0.5/π*(λ0/(2.0*π))^2) )^2 * ( (θmaj_rad^2 + θmin_rad^2)) + FWHM = (2.0 * log(2.0))^0.5/π + Qbar = 2.0/gamma((2-α)/2.) * (inscale^2*(1.0 + M)/(FWHM*(λ0/(2.0*π))^2) )^2 * ( (θmaj_rad^2 + θmin_rad^2)) C = (λ0/(2.0*π))^2 * Qbar*gamma(1.0 - α/2.0)/(8.0*π^2*inscale^2) - P_ϕ0 = 1.0/(2.0*π* _₂F₁((α + 2.0)/2.0, 0.5, 1.0, -k)) - return new(α, inscale, θmaj, θmin, ϕpt, λ0, M, ζ, A, k, Bmaj, Bmin, C, P_ϕ0) + Amaj = (inscale*(1.0 + M) * θmaj_rad/(FWHM * (λ0/(2.0*π)) * 2.0*π ))^2 + Amin = (inscale*(1.0 + M) * θmin_rad/(FWHM * (λ0/(2.0*π)) * 2.0*π ))^2 + + + return new(α, inscale, θmaj, θmin, ϕpt, λ0, M, ζ, A, k, Bmaj, Bmin, Qbar, C, Amaj, Amin, ϕ0, P_ϕ) end end From be453ef5b679f255b7f6b2521c3513319b6346f1 Mon Sep 17 00:00:00 2001 From: Anna Tartaglia Date: Thu, 20 Jul 2023 16:24:33 -0400 Subject: [PATCH 17/19] change kernel equations for ehtim consistency --- .../abstractscatteringkernel.jl | 3 +-- src/scatteringkernels/phasestructure.jl | 18 ++++++++---------- 2 files changed, 9 insertions(+), 12 deletions(-) diff --git a/src/scatteringkernels/abstractscatteringkernel.jl b/src/scatteringkernels/abstractscatteringkernel.jl index ca5b31d..c6f0be3 100644 --- a/src/scatteringkernels/abstractscatteringkernel.jl +++ b/src/scatteringkernels/abstractscatteringkernel.jl @@ -92,8 +92,7 @@ x and y into polar coordinates function Dϕ_exact(ps::AbstractScatteringKernel, λ::Number, x::Number, y::Number) r = (x^2 + y^2)^0.5 ϕ = atan(y, x) - f(ϕ_q) = dDϕ_dz(ps, λ, r, ϕ, ϕ_q)*P_ϕ(ps, ϕ_q) - return quadgk(f, 0, 2.0*π)[1] + return quadgk(ϕ_q -> dDϕ_dz(ps, λ, r, ϕ, ϕ_q)*ps.P_ϕ(ϕ_q), 0, 2.0*π)[1] end """ diff --git a/src/scatteringkernels/phasestructure.jl b/src/scatteringkernels/phasestructure.jl index aa20fbb..9832a59 100644 --- a/src/scatteringkernels/phasestructure.jl +++ b/src/scatteringkernels/phasestructure.jl @@ -22,7 +22,7 @@ struct PhaseStructure{T<:Number, F<:Function} <: AbstractScatteringKernel Amin::Number ϕ0::Number P_ϕ::F - function PhaseStructure(α=1.38, inscale=800e5, θmaj=1.380, θmin=.703, ϕpt=81.9, λ0=1., M=.53) + function PhaseStructure(α=1.38, inscale=800e5, θmaj=1.380, θmin=.703, ϕpt=81.9, λ0=1., M=2.82/5.53) # solve for k from dipole model in Psaltis 2018 A = θmaj/θmin ζ = (A^2. - 1)/(A^2. +1) @@ -36,22 +36,20 @@ struct PhaseStructure{T<:Number, F<:Function} <: AbstractScatteringKernel P_ϕ0 = 1.0/(2.0*π* _₂F₁((α + 2.0)/2.0, 0.5, 1.0, -k)) P_ϕ(ϕ) = P_ϕ0*(1. +k*sin(ϕ-ϕ0)^2.)^(-(α+2.)/2.) - # B parameters - fmaj(ϕ) = abs(cos(ϕ0-ϕ))^α*P_ϕ(ϕ) - fmin(ϕ) = abs(sin(ϕ0-ϕ))^α*P_ϕ(ϕ) - B_prefac = C*2.0^(2.0-α)*π^0.5/(α*gamma((α+1.0)/2.0)) - Bmaj = B_prefac*quadgk(fmaj, 0, 2*π)[1] - Bmin = B_prefac*quadgk(fmin, 0, 2*π)[1] - # C parameter FWHM = (2.0 * log(2.0))^0.5/π Qbar = 2.0/gamma((2-α)/2.) * (inscale^2*(1.0 + M)/(FWHM*(λ0/(2.0*π))^2) )^2 * ( (θmaj_rad^2 + θmin_rad^2)) C = (λ0/(2.0*π))^2 * Qbar*gamma(1.0 - α/2.0)/(8.0*π^2*inscale^2) + # B parameters + B_prefac = C*2.0^(2.0-α)*π^0.5/(α*gamma((α+1.0)/2.0)) + Bmaj = B_prefac*quadgk(ϕ -> abs(cos(ϕ0-ϕ))^α*P_ϕ(ϕ), 0, 2*π)[1] + Bmin = B_prefac*quadgk(ϕ -> abs(sin(ϕ0-ϕ))^α*P_ϕ(ϕ), 0, 2*π)[1] + + Amaj = (inscale*(1.0 + M) * θmaj_rad/(FWHM * (λ0/(2.0*π)) * 2.0*π ))^2 Amin = (inscale*(1.0 + M) * θmin_rad/(FWHM * (λ0/(2.0*π)) * 2.0*π ))^2 - - return new(α, inscale, θmaj, θmin, ϕpt, λ0, M, ζ, A, k, Bmaj, Bmin, Qbar, C, Amaj, Amin, ϕ0, P_ϕ) + return new{typeof(α), Function}(α, inscale, θmaj, θmin, ϕpt, λ0, M, ζ, A, k, Bmaj, Bmin, Qbar, C, Amaj, Amin, ϕ0, P_ϕ) end end From 7c930cce467b3b99b5d37f5e49ec892e21a22608 Mon Sep 17 00:00:00 2001 From: Kazunori Akiyama Date: Wed, 26 Jul 2023 23:04:21 -0400 Subject: [PATCH 18/19] The entire library was clearned up. The major change includes remaining of source codes, types and methods. Now the library includes - three field wonder models - abstract type, data types and methods for kzeta finding - abstract type, data types and methods for scattering models - abstract type, data types and methods for scattering kernels --- Manifest.toml | 833 +++++++++++++++++- Project.toml | 8 +- src/ScatteringOptics.jl | 34 +- src/kfinders/abstractkfinder.jl | 38 - src/kfinders/kfinder.jl | 12 - src/kzetafinders/abstractkzetafinder.jl | 33 + src/kzetafinders/dipole.jl | 24 + src/kzetafinders/periodicboxcar.jl | 20 + src/kzetafinders/vonmises.jl | 21 + .../abstractscatteringkernel.jl | 110 --- src/scatteringkernels/phasestructure.jl | 55 -- .../abstractscatteringmodel.jl | 148 ++++ src/scatteringmodels/commonfunctions.jl | 29 + src/scatteringmodels/kernels/abstract.jl | 17 + src/scatteringmodels/kernels/approx.jl | 33 + src/scatteringmodels/kernels/exact.jl | 29 + src/scatteringmodels/kernels/kernelmodel.jl | 23 + src/scatteringmodels/models/dipole.jl | 95 ++ src/scatteringmodels/models/periodicboxcar.jl | 92 ++ src/scatteringmodels/models/vonmises.jl | 85 ++ 20 files changed, 1467 insertions(+), 272 deletions(-) delete mode 100644 src/kfinders/abstractkfinder.jl delete mode 100644 src/kfinders/kfinder.jl create mode 100644 src/kzetafinders/abstractkzetafinder.jl create mode 100644 src/kzetafinders/dipole.jl create mode 100644 src/kzetafinders/periodicboxcar.jl create mode 100644 src/kzetafinders/vonmises.jl delete mode 100644 src/scatteringkernels/abstractscatteringkernel.jl delete mode 100644 src/scatteringkernels/phasestructure.jl create mode 100644 src/scatteringmodels/abstractscatteringmodel.jl create mode 100644 src/scatteringmodels/commonfunctions.jl create mode 100644 src/scatteringmodels/kernels/abstract.jl create mode 100644 src/scatteringmodels/kernels/approx.jl create mode 100644 src/scatteringmodels/kernels/exact.jl create mode 100644 src/scatteringmodels/kernels/kernelmodel.jl create mode 100644 src/scatteringmodels/models/dipole.jl create mode 100644 src/scatteringmodels/models/periodicboxcar.jl create mode 100644 src/scatteringmodels/models/vonmises.jl diff --git a/Manifest.toml b/Manifest.toml index dd6a9e1..4ac0d51 100644 --- a/Manifest.toml +++ b/Manifest.toml @@ -1,9 +1,38 @@ # This file is machine-generated - editing it directly is not advised [[ADTypes]] -git-tree-sha1 = "e58c18d2312749847a74f5be80bb0fa53da102bd" +git-tree-sha1 = "f5c25e8a5b29b5e941b7408bc8cc79fea4d9ef9a" uuid = "47edcb42-4c32-4615-8424-f2b9edc5f35b" -version = "0.1.5" +version = "0.1.6" + +[[AbstractFFTs]] +deps = ["LinearAlgebra"] +git-tree-sha1 = "cad4c758c0038eea30394b1b671526921ca85b21" +uuid = "621f4979-c628-5d54-868e-fcf4e3e8185c" +version = "1.4.0" +weakdeps = ["ChainRulesCore"] + + [AbstractFFTs.extensions] + AbstractFFTsChainRulesCoreExt = "ChainRulesCore" + +[[AbstractNFFTs]] +deps = ["LinearAlgebra", "Printf"] +git-tree-sha1 = "292e21e99dedb8621c15f185b8fdb4260bb3c429" +uuid = "7f219486-4aa7-41d6-80a7-e08ef20ceed7" +version = "0.8.2" + +[[Accessors]] +deps = ["CompositionsBase", "ConstructionBase", "Dates", "InverseFunctions", "LinearAlgebra", "MacroTools", "Requires", "Test"] +git-tree-sha1 = "954634616d5846d8e216df1298be2298d55280b2" +uuid = "7d9f7c33-5ae7-4f3b-8dc6-eff91059b697" +version = "0.1.32" +weakdeps = ["AxisKeys", "IntervalSets", "StaticArrays", "StructArrays"] + + [Accessors.extensions] + AccessorsAxisKeysExt = "AxisKeys" + AccessorsIntervalSetsExt = "IntervalSets" + AccessorsStaticArraysExt = "StaticArrays" + AccessorsStructArraysExt = "StructArrays" [[Adapt]] deps = ["LinearAlgebra", "Requires"] @@ -15,6 +44,11 @@ weakdeps = ["StaticArrays"] [Adapt.extensions] AdaptStaticArraysExt = "StaticArrays" +[[ArgCheck]] +git-tree-sha1 = "a3a402a35a2f7e0b87828ccabbd5ebfbebe356b4" +uuid = "dce04be8-c92d-5529-be00-80e4d2c0e197" +version = "2.3.0" + [[ArgTools]] uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f" version = "1.1.1" @@ -56,20 +90,86 @@ version = "0.1.29" [[Artifacts]] uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33" -[[AstroAngles]] -git-tree-sha1 = "41621fa5ed5f7614b75eea8e0b3cfd967b284c87" -uuid = "5c4adb95-c1fc-4c53-b4ea-2a94080c53d2" -version = "0.1.3" +[[AxisKeys]] +deps = ["AbstractFFTs", "ChainRulesCore", "CovarianceEstimation", "IntervalSets", "InvertedIndices", "LazyStack", "LinearAlgebra", "NamedDims", "OffsetArrays", "Statistics", "StatsBase", "Tables"] +git-tree-sha1 = "dba0fdaa3a95e591aa9cbe0df9aba41e295a2011" +uuid = "94b1ba4f-4ee9-5380-92f1-94cde586c3c5" +version = "0.2.13" + +[[BangBang]] +deps = ["Compat", "ConstructionBase", "InitialValues", "LinearAlgebra", "Requires", "Setfield", "Tables"] +git-tree-sha1 = "e28912ce94077686443433c2800104b061a827ed" +uuid = "198e06fe-97b7-11e9-32a5-e1d131e6ad66" +version = "0.3.39" + + [BangBang.extensions] + BangBangChainRulesCoreExt = "ChainRulesCore" + BangBangDataFramesExt = "DataFrames" + BangBangStaticArraysExt = "StaticArrays" + BangBangStructArraysExt = "StructArrays" + BangBangTypedTablesExt = "TypedTables" + + [BangBang.weakdeps] + ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" + DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" + StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" + StructArrays = "09ab397b-f2b6-538f-b94a-2f83cf4a842a" + TypedTables = "9d95f2ec-7b3d-5a63-8d20-e2491e220bb9" [[Base64]] uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f" +[[Baselet]] +git-tree-sha1 = "aebf55e6d7795e02ca500a689d326ac979aaf89e" +uuid = "9718e550-a3fa-408a-8086-8db961cd8217" +version = "0.1.1" + +[[BasicInterpolators]] +deps = ["LinearAlgebra", "Memoize", "Random"] +git-tree-sha1 = "3f7be532673fc4a22825e7884e9e0e876236b12a" +uuid = "26cce99e-4866-4b6d-ab74-862489e035e0" +version = "0.7.1" + +[[Bessels]] +git-tree-sha1 = "4435559dc39793d53a9e3d278e185e920b4619ef" +uuid = "0e736298-9ec6-45e8-9647-e4fc86a2fe38" +version = "0.2.8" + [[BitTwiddlingConvenienceFunctions]] deps = ["Static"] git-tree-sha1 = "0c5f81f47bbbcf4aea7b2959135713459170798b" uuid = "62783981-4cbd-42fc-bca8-16325de8dc4b" version = "0.1.5" +[[Bzip2_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "19a35467a82e236ff51bc17a3a44b69ef35185a2" +uuid = "6e34b625-4abd-537c-b88f-471c36dfa7a0" +version = "1.0.8+0" + +[[CEnum]] +git-tree-sha1 = "eb4cb44a499229b3b8426dcfb5dd85333951ff90" +uuid = "fa961155-64e5-5f13-b03f-caf6b980ea82" +version = "0.4.2" + +[[CFITSIO]] +deps = ["CFITSIO_jll"] +git-tree-sha1 = "8425c47db102577eefb93cb37b4480e750116b0d" +uuid = "3b1b4be9-1499-4b22-8d78-7db3344d1961" +version = "1.4.1" + +[[CFITSIO_jll]] +deps = ["Artifacts", "JLLWrappers", "LibCURL_jll", "Libdl", "Pkg", "Zlib_jll"] +git-tree-sha1 = "9c91a9358de42043c3101e3a29e60883345b0b39" +uuid = "b3e40c51-02ae-5482-8a39-3ace5868dcf4" +version = "4.0.0+0" + +[[CFTime]] +deps = ["Dates", "Printf"] +git-tree-sha1 = "ed2e76c1c3c43fd9d0cb9248674620b29d71f2d1" +uuid = "179af706-886a-5703-950a-314cd64e0468" +version = "0.1.2" + [[CPUSummary]] deps = ["CpuId", "IfElse", "PrecompileTools", "Static"] git-tree-sha1 = "89e0654ed8c7aebad6d5ad235d6242c2d737a928" @@ -94,6 +194,24 @@ git-tree-sha1 = "70232f82ffaab9dc52585e0dd043b5e0c6b714f1" uuid = "fb6a15b2-703c-40df-9091-08a04967cfa9" version = "0.1.12" +[[ColorTypes]] +deps = ["FixedPointNumbers", "Random"] +git-tree-sha1 = "eb7f0f8307f71fac7c606984ea5fb2817275d6e4" +uuid = "3da002f7-5984-5a60-b8a6-cbb66c0b333f" +version = "0.11.4" + +[[Colors]] +deps = ["ColorTypes", "FixedPointNumbers", "Reexport"] +git-tree-sha1 = "fc08e5930ee9a4e03f84bfb5211cb54e7769758a" +uuid = "5ae59095-9a9b-59fe-a467-6f913c188581" +version = "0.12.10" + +[[CommonDataModel]] +deps = ["CFTime", "DataStructures", "Dates", "Preferences", "Printf"] +git-tree-sha1 = "2678b3fc170d582655a14d22867b031b6e43c2d4" +uuid = "1fbeeb36-5f17-413c-809b-666fb144f157" +version = "0.2.4" + [[CommonSolve]] git-tree-sha1 = "0eee5eb66b1cf62cd6ad1b460238e60e4b09400c" uuid = "38540f10-b2f7-11e9-35d8-d573e4eb0ff2" @@ -107,9 +225,9 @@ version = "0.3.0" [[Compat]] deps = ["UUIDs"] -git-tree-sha1 = "4e88377ae7ebeaf29a047aa1ee40826e0b708a5d" +git-tree-sha1 = "5ce999a19f4ca23ea484e92a1774a61b8ca4cf8e" uuid = "34da2185-b29b-5c13-b0c7-acf172513d20" -version = "4.7.0" +version = "4.8.0" weakdeps = ["Dates", "LinearAlgebra"] [Compat.extensions] @@ -118,21 +236,51 @@ weakdeps = ["Dates", "LinearAlgebra"] [[CompilerSupportLibraries_jll]] deps = ["Artifacts", "Libdl"] uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae" -version = "1.0.2+0" +version = "1.0.5+0" + +[[CompositionsBase]] +git-tree-sha1 = "802bb88cd69dfd1509f6670416bd4434015693ad" +uuid = "a33af91c-f02d-484b-be07-31d278c5ca2b" +version = "0.1.2" +weakdeps = ["InverseFunctions"] + + [CompositionsBase.extensions] + CompositionsBaseInverseFunctionsExt = "InverseFunctions" + +[[ComradeBase]] +deps = ["AxisKeys", "ChainRulesCore", "DocStringExtensions", "FITSIO", "FillArrays", "LinearAlgebra", "PolarizedTypes", "PrecompileTools", "RectiGrids", "Reexport", "StaticArrays", "Statistics", "StructArrays"] +git-tree-sha1 = "408c20277b7c21a3a05ba0f07241d81c8a5d480d" +uuid = "6d8c423b-a35f-4ef1-850c-862fe21f82c4" +version = "0.5.2" + +[[CondaPkg]] +deps = ["JSON3", "Markdown", "MicroMamba", "Pidfile", "Pkg", "TOML"] +git-tree-sha1 = "741146cf2ced5859faae76a84b541aa9af1a78bb" +uuid = "992eb4ea-22a4-4c89-a5bb-47a3300528ab" +version = "0.2.18" [[ConstructionBase]] deps = ["LinearAlgebra"] -git-tree-sha1 = "738fec4d684a9a6ee9598a8bfee305b26831f28c" +git-tree-sha1 = "fe2838a593b5f776e1597e086dcd47560d94e816" uuid = "187b0558-2788-49d3-abe0-74a17ed4e7c9" -version = "1.5.2" +version = "1.5.3" +weakdeps = ["IntervalSets", "StaticArrays"] [ConstructionBase.extensions] ConstructionBaseIntervalSetsExt = "IntervalSets" ConstructionBaseStaticArraysExt = "StaticArrays" - [ConstructionBase.weakdeps] - IntervalSets = "8197267c-284f-5f27-9208-e0e47529a953" - StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" +[[ContextVariablesX]] +deps = ["Compat", "Logging", "UUIDs"] +git-tree-sha1 = "25cc3803f1030ab855e383129dcd3dc294e322cc" +uuid = "6add18c4-b38d-439d-96f6-d6bc489c04c5" +version = "0.1.3" + +[[CovarianceEstimation]] +deps = ["LinearAlgebra", "Statistics", "StatsBase"] +git-tree-sha1 = "6711ad240bb8861dda376bad332d3f89e2ac5f30" +uuid = "587fd27a-f159-11e8-2dae-1979310e6154" +version = "0.2.9" [[CpuId]] deps = ["Markdown"] @@ -140,11 +288,22 @@ git-tree-sha1 = "fcbb72b032692610bfbdb15018ac16a36cf2e406" uuid = "adafc99b-e345-5852-983c-f28acb93d879" version = "0.3.1" +[[Crayons]] +git-tree-sha1 = "249fe38abf76d48563e2f4556bebd215aa317e15" +uuid = "a8cc5b0e-0ffa-5ad4-8c14-923d3ee1735f" +version = "4.1.1" + [[DataAPI]] git-tree-sha1 = "8da84edb865b0b5b0100c0666a9bc9a0b71c553c" uuid = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a" version = "1.15.0" +[[DataFrames]] +deps = ["Compat", "DataAPI", "DataStructures", "Future", "InlineStrings", "InvertedIndices", "IteratorInterfaceExtensions", "LinearAlgebra", "Markdown", "Missings", "PooledArrays", "PrecompileTools", "PrettyTables", "Printf", "REPL", "Random", "Reexport", "SentinelArrays", "SortingAlgorithms", "Statistics", "TableTraits", "Tables", "Unicode"] +git-tree-sha1 = "04c738083f29f86e62c8afc341f0967d8717bdb8" +uuid = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" +version = "1.6.1" + [[DataStructures]] deps = ["Compat", "InteractiveUtils", "OrderedCollections"] git-tree-sha1 = "cf25ccb972fec4e4817764d01c82386ae94f77b4" @@ -160,11 +319,22 @@ version = "1.0.0" deps = ["Printf"] uuid = "ade2ca70-3891-5945-98fb-dc099432e06a" +[[DefineSingletons]] +git-tree-sha1 = "0fba8b706d0178b4dc7fd44a96a92382c9065c2c" +uuid = "244e2a9f-e319-4986-a169-4d1fe445cd52" +version = "0.1.2" + +[[DelimitedFiles]] +deps = ["Mmap"] +git-tree-sha1 = "9e2f36d3c96a820c678f2f1f1782582fcf685bae" +uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab" +version = "1.9.1" + [[DiffEqBase]] deps = ["ArrayInterface", "ChainRulesCore", "DataStructures", "DocStringExtensions", "EnumX", "FastBroadcast", "ForwardDiff", "FunctionWrappers", "FunctionWrappersWrappers", "LinearAlgebra", "Logging", "Markdown", "MuladdMacro", "Parameters", "PreallocationTools", "Printf", "RecursiveArrayTools", "Reexport", "Requires", "SciMLBase", "SciMLOperators", "Setfield", "SparseArrays", "Static", "StaticArraysCore", "Statistics", "Tricks", "TruncatedStacktraces", "ZygoteRules"] -git-tree-sha1 = "62c41421bd0facc43dfe4e9776135fe21fd1e1b9" +git-tree-sha1 = "c5692436e7f8279503466db216c74165d1b301e4" uuid = "2b5f629d-d688-5b77-993f-72d75c75574e" -version = "6.126.0" +version = "6.127.0" [DiffEqBase.extensions] DiffEqBaseDistributionsExt = "Distributions" @@ -200,6 +370,12 @@ git-tree-sha1 = "23163d55f885173722d1e4cf0f6110cdbaf7e272" uuid = "b552c78f-8df3-52c6-915a-8e097449b14b" version = "1.15.1" +[[DimensionalData]] +deps = ["Adapt", "ArrayInterfaceCore", "ConstructionBase", "Dates", "Extents", "IntervalSets", "IteratorInterfaceExtensions", "LinearAlgebra", "Random", "RecipesBase", "SnoopPrecompile", "SparseArrays", "Statistics", "TableTraits", "Tables"] +git-tree-sha1 = "dda58a378971eabba69c526ca159ee9ca6715f4f" +uuid = "0703355e-b756-11e9-17c0-8b28908087d0" +version = "0.24.12" + [[Distributed]] deps = ["Random", "Serialization", "Sockets"] uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b" @@ -221,21 +397,110 @@ git-tree-sha1 = "5837a837389fccf076445fce071c8ddaea35a566" uuid = "fa6b7ba4-c1ee-5f82-b5fc-ecf0adba8f74" version = "0.6.8" +[[EHTDimensionalData]] +deps = ["DimensionalData", "EHTUtils", "OrderedCollections", "Reexport"] +git-tree-sha1 = "55c00f1238fa0a371c1c7daa61123fe5855254a5" +uuid = "1ac9d6a4-0fea-4ca3-8a21-bd292c035e25" +version = "0.1.2" + +[[EHTImages]] +deps = ["Dates", "DocStringExtensions", "EHTDimensionalData", "EHTModels", "EHTNCDBase", "EHTUVData", "EHTUtils", "FFTW", "FITSIO", "FLoops", "Formatting", "Logging", "Missings", "NCDatasets", "OrderedCollections", "Parameters", "PythonCall", "PythonPlot", "Unitful", "UnitfulAngles", "UnitfulAstro"] +git-tree-sha1 = "c7ffc09cc002fe824aa23f1658d13e5430bf533d" +uuid = "a1be09f2-4a9d-4f8b-a1e3-812015e25e30" +version = "0.2.2" + +[[EHTModels]] +deps = ["Accessors", "DocStringExtensions", "EHTUtils", "NFFT"] +git-tree-sha1 = "89a98ac9d2369582d4b393ea61fc0f16e37f641a" +uuid = "aba06f68-8285-4329-b7bf-b3157c6bc00b" +version = "0.1.2" + +[[EHTNCDBase]] +deps = ["NCDatasets"] +git-tree-sha1 = "cd499d4eb498cb8208c32e5ff65fddde66a61d1c" +uuid = "0d1a6071-c6ae-40ce-b68e-43bf1c0beb66" +version = "0.1.0" + +[[EHTUVData]] +deps = ["CondaPkg", "DataFrames", "Dates", "EHTDimensionalData", "EHTUtils", "FLoops", "Formatting", "OrderedCollections", "PythonCall", "Statistics", "Unitful", "UnitfulAngles", "UnitfulAstro"] +git-tree-sha1 = "6d6237de8c3ad1e3e80044bbb62d9c13c51ecf5e" +uuid = "56783127-7e15-4e7d-9691-d34b50bc826a" +version = "0.1.4" + +[[EHTUtils]] +deps = ["Dates", "Logging", "PhysicalConstants", "Unitful", "UnitfulAngles", "UnitfulAstro"] +git-tree-sha1 = "69fc67cfb95f5d9e753ffea647069a0f34deac9a" +uuid = "9d0fa6a6-ae25-4c2e-8152-6c4b7f2016aa" +version = "0.1.2" + [[EnumX]] git-tree-sha1 = "bdb1942cd4c45e3c678fd11569d5cccd80976237" uuid = "4e289a0a-7415-4d19-859d-a7e5c4648b56" version = "1.0.4" +[[Enzyme]] +deps = ["CEnum", "EnzymeCore", "Enzyme_jll", "GPUCompiler", "LLVM", "Libdl", "LinearAlgebra", "ObjectFile", "Printf", "Random"] +git-tree-sha1 = "5eb8de9b18c31683be5c7b5a5034aaf5e24dd068" +uuid = "7da242da-08ed-463a-9acd-ee780be4f1d9" +version = "0.11.6" + +[[EnzymeCore]] +deps = ["Adapt"] +git-tree-sha1 = "209c7b307d1a8971e85587a970381bf1e53fda30" +uuid = "f151be2c-9106-41f4-ab19-57ee4f262869" +version = "0.5.1" + +[[Enzyme_jll]] +deps = ["Artifacts", "JLLWrappers", "LazyArtifacts", "Libdl", "TOML"] +git-tree-sha1 = "dca536dc737737d14228b485bb46abeeecfa9695" +uuid = "7cc45869-7501-5eee-bdea-0790c847d4ef" +version = "0.0.78+0" + [[ExprTools]] git-tree-sha1 = "c1d06d129da9f55715c6c212866f5b1bddc5fa00" uuid = "e2ba6199-217a-4e67-a87a-7c52f15ade04" version = "0.1.9" +[[Extents]] +git-tree-sha1 = "5e1e4c53fa39afe63a7d356e30452249365fba99" +uuid = "411431e0-e8b7-467b-b5e0-f676ba4f2910" +version = "0.1.1" + +[[FFTW]] +deps = ["AbstractFFTs", "FFTW_jll", "LinearAlgebra", "MKL_jll", "Preferences", "Reexport"] +git-tree-sha1 = "b4fbdd20c889804969571cc589900803edda16b7" +uuid = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341" +version = "1.7.1" + +[[FFTW_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "c6033cc3892d0ef5bb9cd29b7f2f0331ea5184ea" +uuid = "f5851436-0d7a-5f13-b9de-f02708fd171a" +version = "3.3.10+0" + +[[FITSIO]] +deps = ["CFITSIO", "Printf", "Reexport", "Tables"] +git-tree-sha1 = "a8924c203d66d4c5d72980572c6810213422a59d" +uuid = "525bcba6-941b-5504-bd06-fd0dc1a4d2eb" +version = "0.17.1" + +[[FLoops]] +deps = ["BangBang", "Compat", "FLoopsBase", "InitialValues", "JuliaVariables", "MLStyle", "Serialization", "Setfield", "Transducers"] +git-tree-sha1 = "ffb97765602e3cbe59a0589d237bf07f245a8576" +uuid = "cc61a311-1640-44b5-9fba-1b764f453329" +version = "0.2.1" + +[[FLoopsBase]] +deps = ["ContextVariablesX"] +git-tree-sha1 = "656f7a6859be8673bf1f35da5670246b923964f7" +uuid = "b9860ae5-e623-471e-878b-f6a53c775ea6" +version = "0.1.1" + [[FastBroadcast]] deps = ["ArrayInterface", "LinearAlgebra", "Polyester", "Static", "StaticArrayInterface", "StrideArraysCore"] -git-tree-sha1 = "d1248fceea0b26493fd33e8e9e8c553270da03bd" +git-tree-sha1 = "aa9925a229d45fe3018715238956766fa21804d1" uuid = "7034ab61-46d4-4ed7-9d0f-46aef9175898" -version = "0.2.5" +version = "0.2.6" [[FastLapackInterface]] deps = ["LinearAlgebra"] @@ -246,6 +511,12 @@ version = "2.0.0" [[FileWatching]] uuid = "7b1f6079-737a-58dc-b8bc-7a2ca5c1b5ee" +[[FillArrays]] +deps = ["LinearAlgebra", "Random", "SparseArrays", "Statistics"] +git-tree-sha1 = "f372472e8672b1d993e93dada09e23139b509f9e" +uuid = "1a297f60-69ca-5386-bcde-b61e274b549b" +version = "1.5.0" + [[FiniteDiff]] deps = ["ArrayInterface", "LinearAlgebra", "Requires", "Setfield", "SparseArrays"] git-tree-sha1 = "c6e4a1fbe73b31a3dea94b1da449503b8830c306" @@ -262,6 +533,18 @@ version = "2.21.1" BlockBandedMatrices = "ffab5731-97b5-5995-9138-79e8c1846df0" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" +[[FixedPointNumbers]] +deps = ["Statistics"] +git-tree-sha1 = "335bfdceacc84c5cdf16aadc768aa5ddfc5383cc" +uuid = "53c48c17-4a7d-5ca2-90c5-79b7896eea93" +version = "0.8.4" + +[[Formatting]] +deps = ["Printf"] +git-tree-sha1 = "8339d61043228fdd3eb658d86c926cb282ae72a8" +uuid = "59287772-0a20-5a39-b81b-1366585eb4c0" +version = "0.4.2" + [[ForwardDiff]] deps = ["CommonSubexpressions", "DiffResults", "DiffRules", "LinearAlgebra", "LogExpFunctions", "NaNMath", "Preferences", "Printf", "Random", "SpecialFunctions"] git-tree-sha1 = "00e252f4d706b3d55a8863432e742bf5717b498d" @@ -293,12 +576,24 @@ git-tree-sha1 = "2d6ca471a6c7b536127afccfa7564b5b39227fe0" uuid = "46192b85-c4d5-4398-a991-12ede77f4527" version = "0.1.5" +[[GPUCompiler]] +deps = ["ExprTools", "InteractiveUtils", "LLVM", "Libdl", "Logging", "Scratch", "TimerOutputs", "UUIDs"] +git-tree-sha1 = "72b2e3c2ba583d1a7aa35129e56cf92e07c083e3" +uuid = "61eb1bfa-7361-4325-ad38-22787b887f55" +version = "0.21.4" + [[Graphs]] deps = ["ArnoldiMethod", "Compat", "DataStructures", "Distributed", "Inflate", "LinearAlgebra", "Random", "SharedArrays", "SimpleTraits", "SparseArrays", "Statistics"] git-tree-sha1 = "1cf1d7dcb4bc32d7b4a5add4232db3750c27ecb4" uuid = "86223c79-3864-5bf0-83f7-82e725a168b6" version = "1.8.0" +[[HDF5_jll]] +deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "LLVMOpenMP_jll", "LazyArtifacts", "LibCURL_jll", "Libdl", "MPICH_jll", "MPIPreferences", "MPItrampoline_jll", "MicrosoftMPI_jll", "OpenMPI_jll", "OpenSSL_jll", "TOML", "Zlib_jll", "libaec_jll"] +git-tree-sha1 = "3b20c3ce9c14aedd0adca2bc8c882927844bd53d" +uuid = "0234f1f7-429e-5d53-9886-15a909be8d59" +version = "1.14.0+0" + [[HostCPUFeatures]] deps = ["BitTwiddlingConvenienceFunctions", "IfElse", "Libdl", "Static"] git-tree-sha1 = "d38bd0d9759e3c6cfa19bdccc314eccf8ce596cc" @@ -307,9 +602,9 @@ version = "0.1.15" [[HypergeometricFunctions]] deps = ["DualNumbers", "LinearAlgebra", "OpenLibm_jll", "SpecialFunctions"] -git-tree-sha1 = "ce7ea9cc5db29563b1fe20196b6d23ab3b111384" +git-tree-sha1 = "83e95aaab9dc184a6dcd9c4c52aa0dc26cd14a1d" uuid = "34004b35-14d8-5ef3-9330-4cdb6864b03a" -version = "0.3.18" +version = "0.3.21" [[IfElse]] git-tree-sha1 = "debdd00ffef04665ccbb3e150747a77560e8fad1" @@ -321,10 +616,44 @@ git-tree-sha1 = "5cd07aab533df5170988219191dfad0519391428" uuid = "d25df0c9-e2be-5dd7-82c8-3ad0b3e990b9" version = "0.1.3" +[[InitialValues]] +git-tree-sha1 = "4da0f88e9a39111c2fa3add390ab15f3a44f3ca3" +uuid = "22cec73e-a1b8-11e9-2c92-598750a2cf9c" +version = "0.3.1" + +[[InlineStrings]] +deps = ["Parsers"] +git-tree-sha1 = "9cc2baf75c6d09f9da536ddf58eb2f29dedaf461" +uuid = "842dd82b-1e85-43dc-bf29-5d0ee9dffc48" +version = "1.4.0" + +[[IntelOpenMP_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "0cb9352ef2e01574eeebdb102948a58740dcaf83" +uuid = "1d5cc7b8-4909-519e-a0f8-d0f5ad9712d0" +version = "2023.1.0+0" + [[InteractiveUtils]] deps = ["Markdown"] uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240" +[[IntervalSets]] +deps = ["Dates", "Random", "Statistics"] +git-tree-sha1 = "16c0cc91853084cb5f58a78bd209513900206ce6" +uuid = "8197267c-284f-5f27-9208-e0e47529a953" +version = "0.7.4" + +[[InverseFunctions]] +deps = ["Test"] +git-tree-sha1 = "eabe3125edba5c9c10b60a160b1779a000dc8b29" +uuid = "3587e190-3f89-42d0-90ee-14403ec27112" +version = "0.1.11" + +[[InvertedIndices]] +git-tree-sha1 = "0dc7b50b8d436461be01300fd8cd45aa0274b038" +uuid = "41ab1584-1d38-5bbf-9106-f11c6c58b48f" +version = "1.3.0" + [[IrrationalConstants]] git-tree-sha1 = "630b497eafcc20001bba38a4651b327dcfc491d2" uuid = "92d709cd-6900-40b7-9082-c6be49f344b6" @@ -341,6 +670,18 @@ git-tree-sha1 = "abc9885a7ca2052a736a600f7fa66209f96506e1" uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210" version = "1.4.1" +[[JSON3]] +deps = ["Dates", "Mmap", "Parsers", "PrecompileTools", "StructTypes", "UUIDs"] +git-tree-sha1 = "5b62d93f2582b09e469b3099d839c2d2ebf5066d" +uuid = "0f8b85d8-7281-11e9-16c2-39a750bddbf1" +version = "1.13.1" + +[[JuliaVariables]] +deps = ["MLStyle", "NameResolution"] +git-tree-sha1 = "49fb3cb53362ddadb4415e9b73926d6b40709e70" +uuid = "b14d175d-62b4-44ba-8fb7-3064adc8c3ec" +version = "0.2.4" + [[KLU]] deps = ["LinearAlgebra", "SparseArrays", "SuiteSparse_jll"] git-tree-sha1 = "764164ed65c30738750965d55652db9c94c59bfe" @@ -349,9 +690,32 @@ version = "0.4.0" [[Krylov]] deps = ["LinearAlgebra", "Printf", "SparseArrays"] -git-tree-sha1 = "0356a64062656b0cbb43c504ad5de338251f4bda" +git-tree-sha1 = "6dc4ad9cd74ad4ca0a8e219e945dbd22039f2125" uuid = "ba0b0d4f-ebba-5204-a429-3ac8c609bfb7" -version = "0.9.1" +version = "0.9.2" + +[[LLVM]] +deps = ["CEnum", "LLVMExtra_jll", "Libdl", "Printf", "Unicode"] +git-tree-sha1 = "8695a49bfe05a2dc0feeefd06b4ca6361a018729" +uuid = "929cbde3-209d-540e-8aea-75f648917ca0" +version = "6.1.0" + +[[LLVMExtra_jll]] +deps = ["Artifacts", "JLLWrappers", "LazyArtifacts", "Libdl", "TOML"] +git-tree-sha1 = "c35203c1e1002747da220ffc3c0762ce7754b08c" +uuid = "dad2f222-ce93-54a1-a47d-0025e8a3acab" +version = "0.0.23+0" + +[[LLVMOpenMP_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "f689897ccbe049adb19a065c495e75f372ecd42b" +uuid = "1d63c593-3942-5779-bab2-d838dc0a180e" +version = "15.0.4+0" + +[[LaTeXStrings]] +git-tree-sha1 = "f2355693d6778a178ade15952b7ac47a4ff97996" +uuid = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f" +version = "1.3.0" [[LayoutPointers]] deps = ["ArrayInterface", "LinearAlgebra", "ManualMemory", "SIMDTypes", "Static", "StaticArrayInterface"] @@ -365,6 +729,16 @@ git-tree-sha1 = "1370f8202dac30758f3c345f9909b97f53d87d3f" uuid = "50d2b5c4-7a5e-59d5-8109-a42b560f39c0" version = "0.15.1" +[[LazyArtifacts]] +deps = ["Artifacts", "Pkg"] +uuid = "4af54fe1-eca0-43a8-85a7-787d91b784e3" + +[[LazyStack]] +deps = ["ChainRulesCore", "LinearAlgebra", "NamedDims", "OffsetArrays"] +git-tree-sha1 = "2eb4a5bf2eb0519ebf40c797ba5637d327863637" +uuid = "1fad7336-0346-5a1a-a56f-a06ba010965b" +version = "0.0.8" + [[LibCURL]] deps = ["LibCURL_jll", "MozillaCACerts_jll"] uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21" @@ -387,15 +761,21 @@ version = "1.10.2+0" [[Libdl]] uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb" +[[Libiconv_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "c7cb1f5d892775ba13767a87c7ada0b980ea0a71" +uuid = "94ce4f54-9a6c-5748-9c1c-f9c7231a4531" +version = "1.16.1+2" + [[LinearAlgebra]] deps = ["Libdl", "OpenBLAS_jll", "libblastrampoline_jll"] uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" [[LinearSolve]] deps = ["ArrayInterface", "DocStringExtensions", "EnumX", "FastLapackInterface", "GPUArraysCore", "InteractiveUtils", "KLU", "Krylov", "LinearAlgebra", "PrecompileTools", "Preferences", "RecursiveFactorization", "Reexport", "Requires", "SciMLBase", "SciMLOperators", "Setfield", "SparseArrays", "Sparspak", "SuiteSparse", "UnPack"] -git-tree-sha1 = "93f3a0d88c8f5498bed3ad37ddd844b939fdf899" +git-tree-sha1 = "1b55771f2c211583ad52af5a5ca6475be374c961" uuid = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae" -version = "2.3.0" +version = "2.4.1" [LinearSolve.extensions] LinearSolveCUDAExt = "CUDA" @@ -432,15 +812,44 @@ uuid = "56ddb016-857b-54e1-b83d-db4d58db5568" [[LoopVectorization]] deps = ["ArrayInterface", "ArrayInterfaceCore", "CPUSummary", "CloseOpenIntervals", "DocStringExtensions", "HostCPUFeatures", "IfElse", "LayoutPointers", "LinearAlgebra", "OffsetArrays", "PolyesterWeave", "PrecompileTools", "SIMDTypes", "SLEEFPirates", "Static", "StaticArrayInterface", "ThreadingUtilities", "UnPack", "VectorizationBase"] -git-tree-sha1 = "24e6c5697a6c93b5e10af2acf95f0b2e15303332" +git-tree-sha1 = "c88a4afe1703d731b1c4fdf4e3c7e77e3b176ea2" uuid = "bdcacae8-1622-11e9-2a5c-532679323890" -version = "0.12.163" +version = "0.12.165" weakdeps = ["ChainRulesCore", "ForwardDiff", "SpecialFunctions"] [LoopVectorization.extensions] ForwardDiffExt = ["ChainRulesCore", "ForwardDiff"] SpecialFunctionsExt = "SpecialFunctions" +[[MKL_jll]] +deps = ["Artifacts", "IntelOpenMP_jll", "JLLWrappers", "LazyArtifacts", "Libdl", "Pkg"] +git-tree-sha1 = "154d7aaa82d24db6d8f7e4ffcfe596f40bff214b" +uuid = "856f044c-d86e-5d09-b602-aeab76dc8ba7" +version = "2023.1.0+0" + +[[MLStyle]] +git-tree-sha1 = "bc38dff0548128765760c79eb7388a4b37fae2c8" +uuid = "d8e11817-5142-5d16-987a-aa16d5891078" +version = "0.4.17" + +[[MPICH_jll]] +deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "LazyArtifacts", "Libdl", "MPIPreferences", "TOML"] +git-tree-sha1 = "8a5b4d2220377d1ece13f49438d71ad20cf1ba83" +uuid = "7cb0a576-ebde-5e09-9194-50597f1243b4" +version = "4.1.2+0" + +[[MPIPreferences]] +deps = ["Libdl", "Preferences"] +git-tree-sha1 = "781916a2ebf2841467cda03b6f1af43e23839d85" +uuid = "3da0fdf6-3ccc-4f1b-acd9-58baa6c99267" +version = "0.1.9" + +[[MPItrampoline_jll]] +deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "LazyArtifacts", "Libdl", "MPIPreferences", "TOML"] +git-tree-sha1 = "6979eccb6a9edbbb62681e158443e79ecc0d056a" +uuid = "f1f71cc9-e9ae-5b93-9b94-4fe0e1ad3748" +version = "5.3.1+0" + [[MacroTools]] deps = ["Markdown", "Random"] git-tree-sha1 = "42324d08725e200c23d4dfb549e0d5d89dede2d2" @@ -461,6 +870,42 @@ deps = ["Artifacts", "Libdl"] uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1" version = "2.28.2+0" +[[Measurements]] +deps = ["Calculus", "LinearAlgebra", "Printf", "RecipesBase", "Requires"] +git-tree-sha1 = "51d946d38d62709d6a2d37ea9bcc30c80c686801" +uuid = "eff96d63-e80a-5855-80a2-b1b0885c5ab7" +version = "2.9.0" + +[[Memoize]] +deps = ["MacroTools"] +git-tree-sha1 = "2b1dfcba103de714d31c033b5dacc2e4a12c7caa" +uuid = "c03570c3-d221-55d1-a50c-7939bbd78826" +version = "0.4.4" + +[[MicroCollections]] +deps = ["BangBang", "InitialValues", "Setfield"] +git-tree-sha1 = "629afd7d10dbc6935ec59b32daeb33bc4460a42e" +uuid = "128add7d-3638-4c79-886c-908ea0c25c34" +version = "0.1.4" + +[[MicroMamba]] +deps = ["Pkg", "Scratch", "micromamba_jll"] +git-tree-sha1 = "011cab361eae7bcd7d278f0a7a00ff9c69000c51" +uuid = "0b3b1443-0f03-428d-bdfb-f27f9c1191ea" +version = "0.1.14" + +[[MicrosoftMPI_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "a8027af3d1743b3bfae34e54872359fdebb31422" +uuid = "9237b28f-5490-5468-be7b-bb81f5f5e6cf" +version = "10.1.3+4" + +[[Missings]] +deps = ["DataAPI"] +git-tree-sha1 = "f66bdc5de519e8f8ae43bdc598782d35a25b1272" +uuid = "e1d29d7a-bbdc-5cf2-9ac0-f12de2c33e28" +version = "1.1.0" + [[Mmap]] uuid = "a63ad114-7e13-5084-954f-fe012c677804" @@ -473,12 +918,47 @@ git-tree-sha1 = "cac9cc5499c25554cba55cd3c30543cff5ca4fab" uuid = "46d2c3a1-f734-5fdb-9937-b9b9aeba4221" version = "0.2.4" +[[NCDatasets]] +deps = ["CFTime", "CommonDataModel", "DataStructures", "Dates", "NetCDF_jll", "NetworkOptions", "Printf"] +git-tree-sha1 = "4263c4220f22e20729329838bf7e94a49d1ac32f" +uuid = "85f8d34a-cbdd-5861-8df4-14fed0d494ab" +version = "0.12.17" + +[[NFFT]] +deps = ["AbstractNFFTs", "BasicInterpolators", "Distributed", "FFTW", "FLoops", "LinearAlgebra", "Printf", "Random", "Reexport", "SnoopPrecompile", "SparseArrays", "SpecialFunctions"] +git-tree-sha1 = "93a5f32dd6cf09456b0b81afcb8fc29f06535ffd" +uuid = "efe261a4-0d2b-5849-be55-fc731d526b0d" +version = "0.13.3" + [[NaNMath]] deps = ["OpenLibm_jll"] git-tree-sha1 = "0877504529a3e5c3343c6f8b4c0381e57e4387e4" uuid = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3" version = "1.0.2" +[[NameResolution]] +deps = ["PrettyPrint"] +git-tree-sha1 = "1a0fa0e9613f46c9b8c11eee38ebb4f590013c5e" +uuid = "71a1bf82-56d0-4bbc-8a3c-48b961074391" +version = "0.1.5" + +[[NamedDims]] +deps = ["AbstractFFTs", "ChainRulesCore", "CovarianceEstimation", "LinearAlgebra", "Pkg", "Requires", "Statistics"] +git-tree-sha1 = "dc9144f80a79b302b48c282ad29b1dc2f10a9792" +uuid = "356022a1-0364-5f58-8944-0da4b18d706f" +version = "1.2.1" + +[[NamedTupleTools]] +git-tree-sha1 = "90914795fc59df44120fe3fff6742bb0d7adb1d0" +uuid = "d9ec5142-1e00-5aa0-9d6a-321866360f50" +version = "0.14.3" + +[[NetCDF_jll]] +deps = ["Artifacts", "Bzip2_jll", "HDF5_jll", "JLLWrappers", "LibCURL_jll", "Libdl", "XML2_jll", "Zlib_jll", "Zstd_jll"] +git-tree-sha1 = "10c612c81eaffdd6b7c28a45a554cdd9d2f40ff1" +uuid = "7243133f-43d8-5620-bbf4-c2c921802cf3" +version = "400.902.208+0" + [[NetworkOptions]] uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908" version = "1.2.0" @@ -489,6 +969,12 @@ git-tree-sha1 = "2a7f28c62eb2c16b9c375c38f664cdcf22313cf5" uuid = "8913a72c-1f9b-4ce2-8d82-65094dcecaec" version = "1.8.0" +[[ObjectFile]] +deps = ["Reexport", "StructIO"] +git-tree-sha1 = "69607899b46e1f8ead70396bc51a4c361478d8f6" +uuid = "d8793406-e978-5875-9003-1fc021f44a92" +version = "0.4.0" + [[OffsetArrays]] deps = ["Adapt"] git-tree-sha1 = "2ac17d29c523ce1cd38e27785a7d23024853a4bb" @@ -505,6 +991,18 @@ deps = ["Artifacts", "Libdl"] uuid = "05823500-19ac-5b8b-9628-191a04bc5112" version = "0.8.1+0" +[[OpenMPI_jll]] +deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "LazyArtifacts", "Libdl", "MPIPreferences", "TOML"] +git-tree-sha1 = "f3080f4212a8ba2ceb10a34b938601b862094314" +uuid = "fe0851c0-eecd-5654-98d4-656369965a5c" +version = "4.1.5+0" + +[[OpenSSL_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl"] +git-tree-sha1 = "cae3153c7f6cf3f069a853883fd1919a6e5bab5b" +uuid = "458c3c95-2e84-50aa-8efc-19380b2a3a95" +version = "3.0.9+0" + [[OpenSpecFun_jll]] deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Pkg"] git-tree-sha1 = "13652491f6856acfd2db29360e1bbcd4565d04f1" @@ -512,9 +1010,21 @@ uuid = "efe28fd5-8261-553b-a9e1-b2916fc3738e" version = "0.5.5+0" [[OrderedCollections]] -git-tree-sha1 = "d321bf2de576bf25ec4d3e4360faca399afca282" +git-tree-sha1 = "2e73fe17cac3c62ad1aebe70d44c963c3cfdc3e3" uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d" -version = "1.6.0" +version = "1.6.2" + +[[PackageExtensionCompat]] +git-tree-sha1 = "32f3d52212a8d1c5d589a58851b1f04c97339110" +uuid = "65ce6f38-6b18-4e1d-a461-8949797d7930" +version = "1.0.0" +weakdeps = ["Requires", "TOML"] + +[[PaddedViews]] +deps = ["OffsetArrays"] +git-tree-sha1 = "0fac6313486baae819364c52b4f483450a9d793f" +uuid = "5432bcbf-9aad-5242-b902-cca2824c8663" +version = "0.5.12" [[Parameters]] deps = ["OrderedCollections", "UnPack"] @@ -522,16 +1032,40 @@ git-tree-sha1 = "34c0e9ad262e5f7fc75b10a9952ca7692cfc5fbe" uuid = "d96e819e-fc66-5662-9728-84c9c7592b0a" version = "0.12.3" +[[Parsers]] +deps = ["Dates", "PrecompileTools", "UUIDs"] +git-tree-sha1 = "4b2e829ee66d4218e0cef22c0a64ee37cf258c29" +uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0" +version = "2.7.1" + +[[PhysicalConstants]] +deps = ["Measurements", "Roots", "Unitful"] +git-tree-sha1 = "cd4da9d1890bc2204b08fe95ebafa55e9366ae4e" +uuid = "5ad8b20f-a522-5ce9-bfc9-ddf1d5bda6ab" +version = "0.2.3" + +[[Pidfile]] +deps = ["FileWatching", "Test"] +git-tree-sha1 = "2d8aaf8ee10df53d0dfb9b8ee44ae7c04ced2b03" +uuid = "fa939f87-e72e-5be4-a000-7fc836dbe307" +version = "1.3.0" + [[Pkg]] deps = ["Artifacts", "Dates", "Downloads", "FileWatching", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"] uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" -version = "1.9.0" +version = "1.9.2" + +[[PolarizedTypes]] +deps = ["ChainRulesCore", "DocStringExtensions", "PrecompileTools", "StaticArrays"] +git-tree-sha1 = "4f9da54a2a92e24104f1fd1770897e4d366ec99d" +uuid = "d3c5d4cd-a8ee-40d6-aac7-e34df5a20044" +version = "0.1.0" [[Polyester]] deps = ["ArrayInterface", "BitTwiddlingConvenienceFunctions", "CPUSummary", "IfElse", "ManualMemory", "PolyesterWeave", "Requires", "Static", "StaticArrayInterface", "StrideArraysCore", "ThreadingUtilities"] -git-tree-sha1 = "0fe4e7c4d8ff4c70bfa507f0dd96fa161b115777" +git-tree-sha1 = "3d811babe092a6e7b130beee84998fe7663348b6" uuid = "f517fe37-dbe3-4b94-8317-1923a5111588" -version = "0.7.3" +version = "0.7.5" [[PolyesterWeave]] deps = ["BitTwiddlingConvenienceFunctions", "CPUSummary", "IfElse", "Static", "ThreadingUtilities"] @@ -539,6 +1073,12 @@ git-tree-sha1 = "240d7170f5ffdb285f9427b92333c3463bf65bf6" uuid = "1d0040c9-8b98-4ee7-8388-3f51789ca0ad" version = "0.2.1" +[[PooledArrays]] +deps = ["DataAPI", "Future"] +git-tree-sha1 = "a6062fe4063cdafe78f4a0a81cfffb89721b30e7" +uuid = "2dfb63ee-cc39-5dd5-95bd-886bf059d720" +version = "1.4.2" + [[PreallocationTools]] deps = ["Adapt", "ArrayInterface", "ForwardDiff", "Requires"] git-tree-sha1 = "f739b1b3cc7b9949af3b35089931f2b58c289163" @@ -563,10 +1103,33 @@ git-tree-sha1 = "7eb1686b4f04b82f96ed7a4ea5890a4f0c7a09f1" uuid = "21216c6a-2e73-6563-6e65-726566657250" version = "1.4.0" +[[PrettyPrint]] +git-tree-sha1 = "632eb4abab3449ab30c5e1afaa874f0b98b586e4" +uuid = "8162dcfd-2161-5ef2-ae6c-7681170c5f98" +version = "0.2.0" + +[[PrettyTables]] +deps = ["Crayons", "LaTeXStrings", "Markdown", "Printf", "Reexport", "StringManipulation", "Tables"] +git-tree-sha1 = "ee094908d720185ddbdc58dbe0c1cbe35453ec7a" +uuid = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d" +version = "2.2.7" + [[Printf]] deps = ["Unicode"] uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7" +[[PythonCall]] +deps = ["CondaPkg", "Dates", "Libdl", "MacroTools", "Markdown", "Pkg", "REPL", "Requires", "Serialization", "Tables", "UnsafePointers"] +git-tree-sha1 = "70af6bdbde63d7d0a4ea99f3e890ebdb55e9d464" +uuid = "6099a3de-0909-46bc-b1f4-468b9a2dfc0d" +version = "0.9.14" + +[[PythonPlot]] +deps = ["Colors", "CondaPkg", "LaTeXStrings", "PythonCall", "Sockets", "Test", "VersionParsing"] +git-tree-sha1 = "63a18f73425ce301fef12c3f1fb201b720c05015" +uuid = "274fc56d-3b97-40fa-a1cd-1b4a50311bf9" +version = "1.0.2" + [[QuadGK]] deps = ["DataStructures", "LinearAlgebra"] git-tree-sha1 = "6ec7ac8412e83d57e313393220879ede1740f9ee" @@ -587,11 +1150,17 @@ git-tree-sha1 = "5c3d09cc4f31f5fc6af001c250bf1278733100ff" uuid = "3cdcf5f2-1ef4-517c-9805-6587b60abb01" version = "1.3.4" +[[RectiGrids]] +deps = ["AxisKeys", "ConstructionBase", "Random"] +git-tree-sha1 = "ce3a813e8ddbfc0c7516e440a8c18464f6b7a5d0" +uuid = "8ac6971d-971d-971d-971d-971d5ab1a71a" +version = "0.1.18" + [[RecursiveArrayTools]] deps = ["Adapt", "ArrayInterface", "DocStringExtensions", "GPUArraysCore", "IteratorInterfaceExtensions", "LinearAlgebra", "RecipesBase", "Requires", "StaticArraysCore", "Statistics", "SymbolicIndexingInterface", "Tables"] -git-tree-sha1 = "02ef02926f30d53b94be443bfaea010c47f6b556" +git-tree-sha1 = "7ed35fb5f831aaf09c2d7c8736d44667a1afdcb0" uuid = "731186ca-8d62-57ce-b412-fbd966d074cd" -version = "2.38.5" +version = "2.38.7" [RecursiveArrayTools.extensions] RecursiveArrayToolsMeasurementsExt = "Measurements" @@ -620,6 +1189,22 @@ git-tree-sha1 = "838a3a4188e2ded87a4f9f184b4b0d78a1e91cb7" uuid = "ae029012-a4dd-5104-9daa-d747884805df" version = "1.3.0" +[[Roots]] +deps = ["ChainRulesCore", "CommonSolve", "Printf", "Setfield"] +git-tree-sha1 = "de432823e8aab4dd1a985be4be768f95acf152d4" +uuid = "f2b01f46-fcfa-551c-844a-d8ac1e96c665" +version = "2.0.17" + + [Roots.extensions] + RootsForwardDiffExt = "ForwardDiff" + RootsIntervalRootFindingExt = "IntervalRootFinding" + RootsSymPyExt = "SymPy" + + [Roots.weakdeps] + ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" + IntervalRootFinding = "d2bf35a9-74e0-55ec-b149-d360ff49b807" + SymPy = "24249f21-da20-56a4-8eb1-6a02cf4ae2e6" + [[RuntimeGeneratedFunctions]] deps = ["ExprTools", "SHA", "Serialization"] git-tree-sha1 = "237edc1563bbf078629b4f8d194bd334e97907cf" @@ -643,15 +1228,27 @@ version = "0.6.39" [[SciMLBase]] deps = ["ADTypes", "ArrayInterface", "CommonSolve", "ConstructionBase", "Distributed", "DocStringExtensions", "EnumX", "FunctionWrappersWrappers", "IteratorInterfaceExtensions", "LinearAlgebra", "Logging", "Markdown", "PrecompileTools", "Preferences", "RecipesBase", "RecursiveArrayTools", "Reexport", "RuntimeGeneratedFunctions", "SciMLOperators", "StaticArraysCore", "Statistics", "SymbolicIndexingInterface", "Tables", "TruncatedStacktraces"] -git-tree-sha1 = "a22a12db91f6a921e28a7ae39a9546eed93fd92e" +git-tree-sha1 = "04370090604cd399db5bebddb636d80ab9d338e9" uuid = "0bca4576-84f4-4d90-8ffe-ffa030f20462" -version = "1.93.0" +version = "1.94.0" [[SciMLOperators]] deps = ["ArrayInterface", "DocStringExtensions", "Lazy", "LinearAlgebra", "Setfield", "SparseArrays", "StaticArraysCore", "Tricks"] -git-tree-sha1 = "b1fe33c9984c6789b58419e62e7a2b92f9aa813e" +git-tree-sha1 = "745755a5b932c9a664d7e9e4beb60c692b211d4b" uuid = "c0aeaf25-5076-4817-a8d5-81caf7dfa961" -version = "0.3.3" +version = "0.3.5" + +[[Scratch]] +deps = ["Dates"] +git-tree-sha1 = "30449ee12237627992a99d5e30ae63e4d78cd24a" +uuid = "6c6a2e73-6563-6170-7368-637461726353" +version = "1.2.0" + +[[SentinelArrays]] +deps = ["Dates", "Random"] +git-tree-sha1 = "04bdff0b09c65ff3e06a05e3eb7b120223da3d39" +uuid = "91c51154-3ec4-41a3-a24f-3f23e20d615c" +version = "1.4.0" [[Serialization]] uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b" @@ -667,13 +1264,13 @@ deps = ["Distributed", "Mmap", "Random", "Serialization"] uuid = "1a1011a3-84de-559e-8e89-a11a2f7dc383" [[SimpleNonlinearSolve]] -deps = ["ArrayInterface", "DiffEqBase", "FiniteDiff", "ForwardDiff", "LinearAlgebra", "PrecompileTools", "Reexport", "Requires", "SciMLBase", "StaticArraysCore"] -git-tree-sha1 = "56aa73a93cdca493af5155a0338a864ed314222b" +deps = ["ArrayInterface", "DiffEqBase", "FiniteDiff", "ForwardDiff", "LinearAlgebra", "PackageExtensionCompat", "PrecompileTools", "Reexport", "SciMLBase", "StaticArraysCore"] +git-tree-sha1 = "91fcc402c4ab978ad5759489db9a9c5a71732f2d" uuid = "727e6d20-b764-4bd8-a329-72de5adea6c7" -version = "0.1.16" +version = "0.1.18" [SimpleNonlinearSolve.extensions] - SimpleBatchedNonlinearSolveExt = "NNlib" + SimpleNonlinearSolveNNlibExt = "NNlib" [SimpleNonlinearSolve.weakdeps] NNlib = "872c559c-99b0-510c-b3b7-b6c96a88d5cd" @@ -693,6 +1290,12 @@ version = "1.0.3" [[Sockets]] uuid = "6462fe0b-24de-5631-8697-dd941f90decc" +[[SortingAlgorithms]] +deps = ["DataStructures"] +git-tree-sha1 = "c60ec5c62180f27efea3ba2908480f8055e17cee" +uuid = "a2af1166-a08f-5f64-846c-94a0d3cef48c" +version = "1.1.1" + [[SparseArrays]] deps = ["Libdl", "LinearAlgebra", "Random", "Serialization", "SuiteSparse_jll"] uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" @@ -725,11 +1328,17 @@ weakdeps = ["ChainRulesCore"] [SpecialFunctions.extensions] SpecialFunctionsChainRulesCoreExt = "ChainRulesCore" +[[SplittablesBase]] +deps = ["Setfield", "Test"] +git-tree-sha1 = "e08a62abc517eb79667d0a29dc08a3b589516bb5" +uuid = "171d559e-b47b-412a-8079-5efa626c420e" +version = "0.1.15" + [[Static]] deps = ["IfElse"] -git-tree-sha1 = "dbde6766fc677423598138a5951269432b0fcc90" +git-tree-sha1 = "f295e0a1da4ca425659c57441bcb59abb035a4bc" uuid = "aedffcd0-7271-4cad-89d0-dc628f76c6d3" -version = "0.8.7" +version = "0.8.8" [[StaticArrayInterface]] deps = ["ArrayInterface", "Compat", "IfElse", "LinearAlgebra", "Requires", "SnoopPrecompile", "SparseArrays", "Static", "SuiteSparse"] @@ -744,30 +1353,65 @@ weakdeps = ["OffsetArrays", "StaticArrays"] [[StaticArrays]] deps = ["LinearAlgebra", "Random", "StaticArraysCore"] -git-tree-sha1 = "0da7e6b70d1bb40b1ace3b576da9ea2992f76318" +git-tree-sha1 = "9cabadf6e7cd2349b6cf49f1915ad2028d65e881" uuid = "90137ffa-7385-5640-81b9-e52037218182" -version = "1.6.0" +version = "1.6.2" weakdeps = ["Statistics"] [StaticArrays.extensions] StaticArraysStatisticsExt = "Statistics" [[StaticArraysCore]] -git-tree-sha1 = "6b7ba252635a5eff6a0b0664a41ee140a1c9e72a" +git-tree-sha1 = "36b3d696ce6366023a0ea192b4cd442268995a0d" uuid = "1e83bf80-4336-4d27-bf5d-d5a4f845583c" -version = "1.4.0" +version = "1.4.2" [[Statistics]] deps = ["LinearAlgebra", "SparseArrays"] uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" version = "1.9.0" +[[StatsAPI]] +deps = ["LinearAlgebra"] +git-tree-sha1 = "45a7769a04a3cf80da1c1c7c60caf932e6f4c9f7" +uuid = "82ae8749-77ed-4fe6-ae5f-f523153014b0" +version = "1.6.0" + +[[StatsBase]] +deps = ["DataAPI", "DataStructures", "LinearAlgebra", "LogExpFunctions", "Missings", "Printf", "Random", "SortingAlgorithms", "SparseArrays", "Statistics", "StatsAPI"] +git-tree-sha1 = "75ebe04c5bed70b91614d684259b661c9e6274a4" +uuid = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" +version = "0.34.0" + [[StrideArraysCore]] deps = ["ArrayInterface", "CloseOpenIntervals", "IfElse", "LayoutPointers", "ManualMemory", "SIMDTypes", "Static", "StaticArrayInterface", "ThreadingUtilities"] git-tree-sha1 = "f02eb61eb5c97b48c153861c72fbbfdddc607e06" uuid = "7792a7ef-975c-4747-a70f-980b88e8d1da" version = "0.4.17" +[[StringManipulation]] +git-tree-sha1 = "46da2434b41f41ac3594ee9816ce5541c6096123" +uuid = "892a3eda-7b42-436c-8928-eab12a02cf0e" +version = "0.3.0" + +[[StructArrays]] +deps = ["Adapt", "DataAPI", "GPUArraysCore", "StaticArraysCore", "Tables"] +git-tree-sha1 = "521a0e828e98bb69042fec1809c1b5a680eb7389" +uuid = "09ab397b-f2b6-538f-b94a-2f83cf4a842a" +version = "0.6.15" + +[[StructIO]] +deps = ["Test"] +git-tree-sha1 = "010dc73c7146869c042b49adcdb6bf528c12e859" +uuid = "53d494c1-5632-5724-8f4c-31dff12d585f" +version = "0.3.0" + +[[StructTypes]] +deps = ["Dates", "UUIDs"] +git-tree-sha1 = "ca4bccb03acf9faaf4137a9abc1881ed1841aa70" +uuid = "856f2bd8-1eba-4b0a-8007-ebc267875bd4" +version = "1.10.0" + [[SuiteSparse]] deps = ["Libdl", "LinearAlgebra", "Serialization", "SparseArrays"] uuid = "4607b0f0-06f3-5cda-b6b1-a6196a1729e9" @@ -815,6 +1459,32 @@ git-tree-sha1 = "eda08f7e9818eb53661b3deb74e3159460dfbc27" uuid = "8290d209-cae3-49c0-8002-c8c24d57dab5" version = "0.5.2" +[[TimerOutputs]] +deps = ["ExprTools", "Printf"] +git-tree-sha1 = "f548a9e9c490030e545f72074a41edfd0e5bcdd7" +uuid = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f" +version = "0.5.23" + +[[Transducers]] +deps = ["Adapt", "ArgCheck", "BangBang", "Baselet", "CompositionsBase", "ConstructionBase", "DefineSingletons", "Distributed", "InitialValues", "Logging", "Markdown", "MicroCollections", "Requires", "Setfield", "SplittablesBase", "Tables"] +git-tree-sha1 = "53bd5978b182fa7c57577bdb452c35e5b4fb73a5" +uuid = "28d57a85-8fef-5791-bfe6-a80928e7c999" +version = "0.4.78" + + [Transducers.extensions] + TransducersBlockArraysExt = "BlockArrays" + TransducersDataFramesExt = "DataFrames" + TransducersLazyArraysExt = "LazyArrays" + TransducersOnlineStatsBaseExt = "OnlineStatsBase" + TransducersReferenceablesExt = "Referenceables" + + [Transducers.weakdeps] + BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e" + DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" + LazyArrays = "5078a376-72f3-5289-bfd5-ec5146d43c02" + OnlineStatsBase = "925886fa-5bf2-5e8e-b522-a9147a512338" + Referenceables = "42d2dcc6-99eb-4e98-b66c-637b7d73030e" + [[TriangularSolve]] deps = ["CloseOpenIntervals", "IfElse", "LayoutPointers", "LinearAlgebra", "LoopVectorization", "Polyester", "Static", "VectorizationBase"] git-tree-sha1 = "31eedbc0b6d07c08a700e26d31298ac27ef330eb" @@ -844,34 +1514,103 @@ version = "1.0.2" [[Unicode]] uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5" +[[Unitful]] +deps = ["Dates", "LinearAlgebra", "Random"] +git-tree-sha1 = "c4d2a349259c8eba66a00a540d550f122a3ab228" +uuid = "1986cc42-f94f-5a68-af5c-568840ba703d" +version = "1.15.0" +weakdeps = ["ConstructionBase", "InverseFunctions"] + + [Unitful.extensions] + ConstructionBaseUnitfulExt = "ConstructionBase" + InverseFunctionsUnitfulExt = "InverseFunctions" + +[[UnitfulAngles]] +deps = ["Dates", "Unitful"] +git-tree-sha1 = "d6cfdb6ddeb388af1aea38d2b9905fa014d92d98" +uuid = "6fb2a4bd-7999-5318-a3b2-8ad61056cd98" +version = "0.6.2" + +[[UnitfulAstro]] +deps = ["Unitful", "UnitfulAngles"] +git-tree-sha1 = "05adf5e3a3bd1038dd50ff6760cddd42380a7260" +uuid = "6112ee07-acf9-5e0f-b108-d242c714bf9f" +version = "1.2.0" + +[[UnsafePointers]] +git-tree-sha1 = "c81331b3b2e60a982be57c046ec91f599ede674a" +uuid = "e17b2a0c-0bdf-430a-bd0c-3a23cae4ff39" +version = "1.0.0" + +[[VLBISkyModels]] +deps = ["AbstractFFTs", "Accessors", "ArgCheck", "AxisKeys", "BasicInterpolators", "ChainRulesCore", "ComradeBase", "DelimitedFiles", "DocStringExtensions", "Enzyme", "EnzymeCore", "FFTW", "FastBroadcast", "FillArrays", "ForwardDiff", "LinearAlgebra", "NFFT", "NamedTupleTools", "PaddedViews", "PolarizedTypes", "Printf", "RecipesBase", "Reexport", "Requires", "SpecialFunctions", "Static", "StaticArrays", "StructArrays"] +git-tree-sha1 = "592aa5fbd89ca01afd5546b143fdcd479816745a" +uuid = "d6343c73-7174-4e0f-bb64-562643efbeca" +version = "0.2.1" + + [VLBISkyModels.extensions] + VLBISkyModelsMakieExt = "Makie" + + [VLBISkyModels.weakdeps] + Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a" + [[VectorizationBase]] deps = ["ArrayInterface", "CPUSummary", "HostCPUFeatures", "IfElse", "LayoutPointers", "Libdl", "LinearAlgebra", "SIMDTypes", "Static", "StaticArrayInterface"] git-tree-sha1 = "b182207d4af54ac64cbc71797765068fdeff475d" uuid = "3d5dd08c-fd9d-11e8-17fa-ed2836048c2f" version = "0.21.64" +[[VersionParsing]] +git-tree-sha1 = "58d6e80b4ee071f5efd07fda82cb9fbe17200868" +uuid = "81def892-9a0e-5fdd-b105-ffc91e053289" +version = "1.3.0" + [[VertexSafeGraphs]] deps = ["Graphs"] git-tree-sha1 = "8351f8d73d7e880bfc042a8b6922684ebeafb35c" uuid = "19fa3120-7c27-5ec5-8db8-b0b0aa330d6f" version = "0.2.0" +[[XML2_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Libiconv_jll", "Pkg", "Zlib_jll"] +git-tree-sha1 = "93c41695bc1c08c46c5899f4fe06d6ead504bb73" +uuid = "02c8fc9c-b97f-50b9-bbe4-9be30ff0a78a" +version = "2.10.3+0" + [[Zlib_jll]] deps = ["Libdl"] uuid = "83775a58-1f1d-513f-b197-d71354ab007a" version = "1.2.13+0" +[[Zstd_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl"] +git-tree-sha1 = "49ce682769cd5de6c72dcf1b94ed7790cd08974c" +uuid = "3161d3a3-bdf6-5164-811a-617609db77b4" +version = "1.5.5+0" + [[ZygoteRules]] deps = ["ChainRulesCore", "MacroTools"] git-tree-sha1 = "977aed5d006b840e2e40c0b48984f7463109046d" uuid = "700de1a5-db45-46bc-99cf-38207098b444" version = "0.2.3" +[[libaec_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl"] +git-tree-sha1 = "eddd19a8dea6b139ea97bdc8a0e2667d4b661720" +uuid = "477f73a3-ac25-53e9-8cc3-50b2fa2566f0" +version = "1.0.6+1" + [[libblastrampoline_jll]] deps = ["Artifacts", "Libdl"] uuid = "8e850b90-86db-534c-a0d3-1478176c7d93" version = "5.8.0+0" +[[micromamba_jll]] +deps = ["Artifacts", "JLLWrappers", "LazyArtifacts", "Libdl"] +git-tree-sha1 = "ed38e87f1a2f42427603a2a188b4ec5d9d28fb44" +uuid = "f8abcde7-e9b7-5caa-b8af-a437887ae8e4" +version = "1.4.7+0" + [[nghttp2_jll]] deps = ["Artifacts", "Libdl"] uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d" diff --git a/Project.toml b/Project.toml index 7cf3fcc..e0d55a6 100644 --- a/Project.toml +++ b/Project.toml @@ -1,15 +1,19 @@ name = "ScatteringOptics" uuid = "e317c96f-4a6f-4ae9-8016-294fb9063ff5" -authors = ["Anna Tartaglia", "Kazunori Akiyama"] +authors = ["Anna Tartaglia", "Kazunori Akiyama", "Paul Tiede"] version = "1.0.0-DEV" [deps] -AstroAngles = "5c4adb95-c1fc-4c53-b4ea-2a94080c53d2" +Bessels = "0e736298-9ec6-45e8-9647-e4fc86a2fe38" +ComradeBase = "6d8c423b-a35f-4ef1-850c-862fe21f82c4" DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" +EHTImages = "a1be09f2-4a9d-4f8b-a1e3-812015e25e30" +EHTUtils = "9d0fa6a6-ae25-4c2e-8152-6c4b7f2016aa" HypergeometricFunctions = "34004b35-14d8-5ef3-9330-4cdb6864b03a" NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec" QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc" SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b" +VLBISkyModels = "d6343c73-7174-4e0f-bb64-562643efbeca" [compat] julia = "1.9" diff --git a/src/ScatteringOptics.jl b/src/ScatteringOptics.jl index 94b8989..5a90105 100644 --- a/src/ScatteringOptics.jl +++ b/src/ScatteringOptics.jl @@ -1,19 +1,37 @@ module ScatteringOptics # Import Modules +using Bessels +import ComradeBase: + AbstractModel, IsPrimitive, IsAnalytic, NotAnalytic, + IsPolarized, NotPolarized, + visanalytic, imanalytic, isprimitive, ispolarized using DocStringExtensions +using EHTUtils: mas2rad using HypergeometricFunctions -using SpecialFunctions using NonlinearSolve -using AstroAngles using QuadGK +using SpecialFunctions + +# kζ finders +include("./kzetafinders/abstractkzetafinder.jl") +include("./kzetafinders/dipole.jl") +include("./kzetafinders/periodicboxcar.jl") +include("./kzetafinders/vonmises.jl") + +# scattering models: general defenitions +include("./scatteringmodels/commonfunctions.jl") +include("./scatteringmodels/abstractscatteringmodel.jl") -#k finders -include("./kfinders/abstractkfinder.jl") -include("./kfinders/kfinder.jl") +# scattering kernels +include("./scatteringmodels/kernels/abstract.jl") +include("./scatteringmodels/kernels/approx.jl") +include("./scatteringmodels/kernels/exact.jl") +include("./scatteringmodels/kernels/kernelmodel.jl") -#scattering kernels -include("./scatteringkernels/abstractscatteringkernel.jl") -include("./scatteringkernels/phasestructure.jl") +# implementation of specific models +include("./scatteringmodels/models/dipole.jl") +include("./scatteringmodels/models/periodicboxcar.jl") +include("./scatteringmodels/models/vonmises.jl") end diff --git a/src/kfinders/abstractkfinder.jl b/src/kfinders/abstractkfinder.jl deleted file mode 100644 index f79a55c..0000000 --- a/src/kfinders/abstractkfinder.jl +++ /dev/null @@ -1,38 +0,0 @@ -export AbstractKFinder -export Psaltis18Eq43 -export findk - -""" - AbstractKFinder - -This is an abstract data type for the hypergeometic computation of the concentration parameter k_ζ,2 for -use in the phase structure function - -**Mandatory fields** -- `α::Number`: spectrum index -- `A::Number`: defined as θmax/θmin (major over minor axis) -""" -abstract type AbstractKFinder end - - -""" - Psaltis18Eq43(k, constants::KFinder) -Equation 43 of Psaltis et al. 2018 -""" -function Psaltis18Eq43(k, constants::AbstractKFinder) - α = constants.α - A = constants.A - ᾱ = (α+2.)/2 - return _₂F₁.(ᾱ, 1/2., 2, .-k) ./ _₂F₁.(ᾱ, 3/2., 2, .-k) .- A^2. -end - -""" - findk(finder::KFinder) - -Solves Psaltis18Eq43 for k given A and α parameters in finder -""" -function findk(finder::AbstractKFinder) - probN = NonlinearProblem(Psaltis18Eq43, [0., 0.], finder) - return solve(probN, NewtonRaphson())[1] # 13 micro s, 10 -end - diff --git a/src/kfinders/kfinder.jl b/src/kfinders/kfinder.jl deleted file mode 100644 index 7ca849c..0000000 --- a/src/kfinders/kfinder.jl +++ /dev/null @@ -1,12 +0,0 @@ -export KFinder - -""" - $(TYPEDEF) - -Data type for a Kfinder object that holds parameters α and A -""" - -struct KFinder{T <: Number, Y <: Number} <: AbstractKFinder - α::T - A::Y -end diff --git a/src/kzetafinders/abstractkzetafinder.jl b/src/kzetafinders/abstractkzetafinder.jl new file mode 100644 index 0000000..456f371 --- /dev/null +++ b/src/kzetafinders/abstractkzetafinder.jl @@ -0,0 +1,33 @@ +""" + AbstractKzetaFinder + +This is an abstract data type to set up equations and provide a solver for the concentration parameter k_ζ +for an anistropic interstellar scattering model. See Psaltis et al. 2018, arxiv::1805.01242v1 for details. + +**Mandatory methods** +- `kzetafinder_equation(kzeta, finder::AbstractKzetaFinder)`: should privide a equation where k will be derived. + +**Methods provided** +- `findkzeta_exact(finder::AbstractKzetaFinder)`: solves the equation for k_ζ,2 given the set of parameters in the finder. +""" +abstract type AbstractKzetaFinder end + + +""" + findkzeta_exact(finder::AbstractKzetaFinder; kwargs...) + +Solves the equation for the concentration parameter k_ζ,2 given the set of parameters in the finder. +It uses NonlinearSolve.jl. +""" +@inline function findkzeta_exact(finder::AbstractKzetaFinder; solver=NewtonRaphson()) + probN = NonlinearProblem(kzetafinder_equation, [0.0, 0.0], finder) + return solve(probN, solver)[1] +end + + +""" + kzetafinder_equation(kzeta, finder::AbstractKzetaFinder) + +This equation privide a root-finding function `f(kzeta, finder)` to find kzeta from the equation `f(kzeta, finder)=0`. +""" +function kzetafinder_equation(kzeta, ::AbstractKzetaFinder) end diff --git a/src/kzetafinders/dipole.jl b/src/kzetafinders/dipole.jl new file mode 100644 index 0000000..ee1e287 --- /dev/null +++ b/src/kzetafinders/dipole.jl @@ -0,0 +1,24 @@ +""" + $(TYPEDEF) + +The finder of the concentration parameter k_ζ,2 for dipole anistropic scattering models. +See Psaltis et al. 2018, arxiv::1805.01242v1 for details. The equation for k_ζ,2 is +given by the equation 43 of Psaltis et al. 2018. + +**Mandatory fields** +- `α::Number`: The power-law index of the phase fluctuations (Kolmogorov is 5/3). +- `A::Number`: The anisotropy parameter of the angular broaderning defined by θmaj/θmin. +""" +struct Dipole_KzetaFinder{T<:Number} <: AbstractKzetaFinder + α::T + A::T +end + + +# Equation 43 of Psaltis et al. 2018. +@inline function kzetafinder_equation(kzeta, finder::Dipole_KzetaFinder) + α = finder.α + A = finder.A + ᾱ = (α + 2) / 2 + return _₂F₁.(ᾱ, 0.5, 2, .-kzeta) ./ _₂F₁.(ᾱ, 1.5, 2, .-kzeta) .- A^2 +end diff --git a/src/kzetafinders/periodicboxcar.jl b/src/kzetafinders/periodicboxcar.jl new file mode 100644 index 0000000..cc68f2f --- /dev/null +++ b/src/kzetafinders/periodicboxcar.jl @@ -0,0 +1,20 @@ +""" + $(TYPEDEF) + +The finder of the concentration parameter k_ζ,3 for the Periodic Box Car anistropic scattering model. +See Psaltis et al. 2018, arxiv::1805.01242v1 for details. The equation for k_ζ,3 is +given by the equation 47 of Psaltis et al. 2018. + +**Mandatory fields** +- `A::Number`: The anisotropy parameter of the angular broaderning defined by θmaj/θmin. +""" +struct PeriodicBoxCar_KzetaFinder{T<:Number} <: AbstractKzetaFinder + ζ0::T +end + + +# Equation 47 of Psaltis et al. 2018. +@inline function kzetafinder_equation(kzeta, finder::PeriodicBoxCar_KzetaFinder) + k̃zeta .= π ./ (1 .+ kzeta) + return sin(k̃zeta) ./ k̃zeta .- finder.ζ0 +end diff --git a/src/kzetafinders/vonmises.jl b/src/kzetafinders/vonmises.jl new file mode 100644 index 0000000..8b6d8d4 --- /dev/null +++ b/src/kzetafinders/vonmises.jl @@ -0,0 +1,21 @@ +""" + $(TYPEDEF) + +The finder of the concentration parameter k_ζ,1 for the von Mises anistropic scattering model. +See Psaltis et al. 2018, arxiv::1805.01242v1 for details. The equation for k_ζ,1 is originally +given by the equation 37 of Psaltis et al. 2018, but this is different in the implementation of +Johnson et al. 2018 in eht-imaging library. We follow eht-imaging's implementation. + +**Mandatory fields** +- `A::Number`: The anisotropy parameter of the angular broaderning defined by θmaj/θmin. +""" +struct vonMises_KzetaFinder{T<:Number} <: AbstractKzetaFinder + A::T +end + + +# Originally Equation 37 of Psaltis et al. 2018, but apparely different equation is used in eht-imaging library. +# We here use the version used in eht-imaging library. +@inline function kzetafinder_equation(kzeta, finder::vonMises_KzetaFinder) + return besseli0.(kzeta) ./ besseli1.(kzeta) .- finder.A^2 +end diff --git a/src/scatteringkernels/abstractscatteringkernel.jl b/src/scatteringkernels/abstractscatteringkernel.jl deleted file mode 100644 index c6f0be3..0000000 --- a/src/scatteringkernels/abstractscatteringkernel.jl +++ /dev/null @@ -1,110 +0,0 @@ -export AbstractScatteringKernel -export calc_dmaj -export calc_dmin -export Dϕ_approx -export P_ϕ -export dDϕ_dz -export Dϕ_exact -export visibility_point - - -""" - AbstractScatteringKernel - -This is an abstract data type for a scattering screen phase structure function. - -**Mandatory fields** -- `α::Number`: spectrum index -- `inscale::Number` inner scale r_in -- `θmaj::Number` reference major axis size -- `θmin::Number` reference minor axis size -- `ϕpt::Number`: position angle of major axis -- `λ0::Number`: reference wavelength -- `M::Number`: magnification D/R -""" - -abstract type AbstractScatteringKernel end - -""" - dmaj(r, ps::AbstractScatteringKernel) - -Maps D_maj(r) for given r -Equation 33 of Psaltis et al. 2018 -""" -@inline function calc_dmaj(ps::AbstractScatteringKernel, λ::Number, r::Number) - Bmaj = ps.Bmaj - Amaj = ps.Amaj - α = ps.α - rin = ps.inscale - - d1 = Bmaj*(2. *Amaj/(α*Bmaj))^(-α/(2. -α)) - d2 = (2. *Amaj/(α*Bmaj))^(2. /(2. -α)) - return d1*((1+d2*(r/rin)^2.)^(α/2.)-1) * (λ/ps.λ0)^2. -end - -""" - dmin(r, ps::AbstractScatteringKernel) - -Maps D_min(r) for given r -Equation 34 of Psaltis et al. 2018 -""" -@inline function calc_dmin(ps::AbstractScatteringKernel, λ::Number, r::Number) - Bmin = ps.Bmin - Amin = ps.Amin - α = ps.α - rin = ps.inscale - - d1 = Bmin*(2. *Amin/(α*Bmin))^(-α/(2. -α)) - d2 = (2. *Amin/(α*Bmin))^(2. /(2. -α)) - return d1*((1+d2*(r/rin)^2.)^(α/2.)-1) * (λ/ps.λ0)^2. -end - -""" - Dϕ_approx(ps::AbstractScatteringKernel, λ::Number, x::Number, y::Number) - -Maps approximate phase structure function D_ϕ(r, ϕ) at observing wavelength λ, first converting -x and y into polar coordinates -Equation 35 of Psaltis et al. 2018 -""" -function Dϕ_approx(ps::AbstractScatteringKernel, λ::Number, x::Number, y::Number) - r = √(x^2. + y^2.) - ϕ = atan(y,x) - dmaj = calc_dmaj(ps, λ, r) - dmin = calc_dmin(ps, λ, r) - add = (dmaj + dmin)/2. - sub = (dmaj - dmin)/2. - return add + sub * cos(2. *(ϕ-ps.ϕ0)) -end - - -function dDϕ_dz(ps::AbstractScatteringKernel, λ::Number, r::Number, ϕ::Number, ϕ_q) -"""differential contribution to the phase structure function -""" - return 4.0 * (λ/ps.λ0)^2. * ps.C/ps.α * (_₁F₁(-ps.α/2.0, 0.5, -r^2. /(4. *ps.inscale^2.)*cos(ϕ - ϕ_q)^2) - 1.) -end - -""" - Dϕ_exact(ps::AbstractScatteringKernel, λ::Number, x::Number, y::Number) - -Maps exact phase structure function D_ϕ(r, ϕ)at observing wavelength λ, first converting -x and y into polar coordinates -""" -function Dϕ_exact(ps::AbstractScatteringKernel, λ::Number, x::Number, y::Number) - r = (x^2 + y^2)^0.5 - ϕ = atan(y, x) - return quadgk(ϕ_q -> dDϕ_dz(ps, λ, r, ϕ, ϕ_q)*ps.P_ϕ(ϕ_q), 0, 2.0*π)[1] -end - -""" - visibility_point(ps::AbstractScatteringKernel, λ::Number, u::Number, v::Number) - -Maps visibility kernel for a given observing wavelength λ and fourier space coordinates u, v -""" -function visibility_point(ps::AbstractScatteringKernel, λ::Number, u::Number, v::Number, use_approximate_form=true) - b = (u, v).* (λ/(1+ps.M)) - if use_approximate_form == true - return exp(-.5 * Dϕ_approx(ps, λ, b...)) - else - return exp(-.5 * Dϕ_exact(ps, λ, b...)) - end -end \ No newline at end of file diff --git a/src/scatteringkernels/phasestructure.jl b/src/scatteringkernels/phasestructure.jl deleted file mode 100644 index 9832a59..0000000 --- a/src/scatteringkernels/phasestructure.jl +++ /dev/null @@ -1,55 +0,0 @@ -export PhaseStructure -""" - $(TYPEDEF) -""" - -struct PhaseStructure{T<:Number, F<:Function} <: AbstractScatteringKernel - α::T - inscale::Number - θmaj::Number - θmin::Number - ϕpt::Number - λ0::Number - M::Number - ζ::Number - A::Number - k::Number - Bmaj::Number - Bmin::Number - Qbar::Number - C::Number - Amaj::Number - Amin::Number - ϕ0::Number - P_ϕ::F - function PhaseStructure(α=1.38, inscale=800e5, θmaj=1.380, θmin=.703, ϕpt=81.9, λ0=1., M=2.82/5.53) - # solve for k from dipole model in Psaltis 2018 - A = θmaj/θmin - ζ = (A^2. - 1)/(A^2. +1) - k = findk(KFinder(α, A)) - - #milliarcseconds to radians - θmaj_rad = dms2rad((0,0,θmaj*10^-3)) - θmin_rad = dms2rad((0,0,θmin*10^-3)) - - ϕ0 = deg2rad(90-ϕpt) - P_ϕ0 = 1.0/(2.0*π* _₂F₁((α + 2.0)/2.0, 0.5, 1.0, -k)) - P_ϕ(ϕ) = P_ϕ0*(1. +k*sin(ϕ-ϕ0)^2.)^(-(α+2.)/2.) - - # C parameter - FWHM = (2.0 * log(2.0))^0.5/π - Qbar = 2.0/gamma((2-α)/2.) * (inscale^2*(1.0 + M)/(FWHM*(λ0/(2.0*π))^2) )^2 * ( (θmaj_rad^2 + θmin_rad^2)) - C = (λ0/(2.0*π))^2 * Qbar*gamma(1.0 - α/2.0)/(8.0*π^2*inscale^2) - - # B parameters - B_prefac = C*2.0^(2.0-α)*π^0.5/(α*gamma((α+1.0)/2.0)) - Bmaj = B_prefac*quadgk(ϕ -> abs(cos(ϕ0-ϕ))^α*P_ϕ(ϕ), 0, 2*π)[1] - Bmin = B_prefac*quadgk(ϕ -> abs(sin(ϕ0-ϕ))^α*P_ϕ(ϕ), 0, 2*π)[1] - - - Amaj = (inscale*(1.0 + M) * θmaj_rad/(FWHM * (λ0/(2.0*π)) * 2.0*π ))^2 - Amin = (inscale*(1.0 + M) * θmin_rad/(FWHM * (λ0/(2.0*π)) * 2.0*π ))^2 - - return new{typeof(α), Function}(α, inscale, θmaj, θmin, ϕpt, λ0, M, ζ, A, k, Bmaj, Bmin, Qbar, C, Amaj, Amin, ϕ0, P_ϕ) - end -end diff --git a/src/scatteringmodels/abstractscatteringmodel.jl b/src/scatteringmodels/abstractscatteringmodel.jl new file mode 100644 index 0000000..8108920 --- /dev/null +++ b/src/scatteringmodels/abstractscatteringmodel.jl @@ -0,0 +1,148 @@ +export AbstractScatteringModel +#export calc_Dmaj +#export calc_Dmin +#export Dϕ_approx +#export Pϕ +#export dDϕ_dz +#export Dϕ_exact +export visibility_point_apporoximate +export visibility_point_exact + + +""" + $(TYPEDEF) + +An abstract anistropic scattering model based on a thin-screen approximation. +In this package, we provide a reference implementation of +the dipole (`DipoleScatteringModel`), von Mises (`vonMisesScatteringModel`) and +periodic Box Car models (`PeriodicBoxCarScatteringModel`) all introduced in Psaltis et al. 2018. + +**Mandatory fields** +The scattering model will be fundamentally governed by the following parameters. +Ideally, a subtype of this abstract model should have a constructor only with these arguments. +- `α::Number`: The power-law index of the phase fluctuations (Kolmogorov is 5/3). +- `rin::Number`: The inner scale of the scattering screen in cm. +- `θmaj::Number`: FWHM in mas of the major axis angular broadening at the specified reference wavelength. +- `θmin::Number`: FWHM in mas of the minor axis angular broadening at the specified reference wavelength. +- `ϕ::Number`: The position angle of the major axis of the scattering in degree. +- `λ0::Number`: The reference wavelength for the scattering model in cm. +- `D::Number`: The distance from the observer to the scattering screen in pc. +- `R::Number`: The distance from the source to the scattering screen in pc. +Furthermore the following parameters need to be precomputed. +- `M::Number`: +- `ζ0::Number`: +- `A::Number`: +- `kζ::Number` +- `Bmaj::Number`: +- `Bmin::Number`: +- `Qbar::Number`: +- `C::Number`: +- `Amaj::Number`: +- `Amin::Number`: +- `ϕ0::Number`: + +** Methods need to be defined ** +- `Pϕ(::AbstractScatteringModel, ϕ::Number)` and `Pϕ(::Type{<:AbstractScatteringModel}, ϕ::Number, args...)`: + Normalized probability distribution describing the distribution of the field wonder. + The function should depend on the field wonder model. +""" +abstract type AbstractScatteringModel end + +function Pϕ(::AbstractScatteringModel, ϕ::Number) end + +""" + Dmaj(r, sm::AbstractScatteringModel) + +Masm D_maj(r) for given r. Based on Equation 33 of Psaltis et al. 2018 +""" +@inline function calc_Dmaj(sm::AbstractScatteringModel, λ::Number, r::Number) + Bmaj = sm.Bmaj + Amaj = sm.Amaj + α = sm.α + rin = sm.rin + + d1 = Bmaj * (2 * Amaj / (α * Bmaj))^(-α / (2 - α)) + d2 = (2 * Amaj / (α * Bmaj))^(2 / (2 - α)) + return d1 * ((1 + d2 * (r / rin)^2)^(α / 2) - 1) * (λ / sm.λ0)^2 +end + + +""" + calc_Dmin(r, sm::AbstractScatteringModel) + +Masm D_min(r) for given r. Based on Equation 34 of Psaltis et al. 2018 +""" +@inline function calc_Dmin(sm::AbstractScatteringModel, λ::Number, r::Number) + Bmin = sm.Bmin + Amin = sm.Amin + α = sm.α + rin = sm.rin + + d1 = Bmin * (2 * Amin / (α * Bmin))^(-α / (2 - α)) + d2 = (2 * Amin / (α * Bmin))^(2 / (2 - α)) + return d1 * ((1 + d2 * (r / rin)^2)^(α / 2) - 1) * (λ / sm.λ0)^2 +end + + +""" + Dϕ_approx(sm::AbstractScatteringModel, λ::Number, x::Number, y::Number) + +Masm approximate phase structure function Dϕ(r, ϕ) at observing wavelength λ, first converting +x and y into polar coordinates. Based on Equation 35 of Psaltis et al. 2018. +""" +@inline function Dϕ_approx(sm::AbstractScatteringModel, λ::Number, x::Number, y::Number) + r = √(x^2 + y^2) + ϕ = atan(y, x) + Dmaj = calc_Dmaj(sm, λ, r) + Dmin = calc_Dmin(sm, λ, r) + add = (Dmaj + Dmin) / 2 + sub = (Dmaj - Dmin) / 2 + return add + sub * cos(2 * (ϕ - sm.ϕ0)) +end + +""" + dDϕ_dz(sm::AbstractScatteringModel, λ::Number, r::Number, ϕ::Number, ϕq) + +Differential contribution to the phase structure function. +""" +@inline function dDϕ_dz(sm::AbstractScatteringModel, λ::Number, r::Number, ϕ::Number, ϕq) + + return 4 * (λ / sm.λ0)^2 * sm.C / sm.α * (_₁F₁(-sm.α / 2, 0.5, -r^2 / (4 * sm.rin^2) * cos(ϕ - ϕq)^2) - 1) +end + + +""" + Dϕ_exact(sm::AbstractScatteringModel, λ::Number, x::Number, y::Number) + +Masm exact phase structure function Dϕ(r, ϕ) at observing wavelength `λ`, first converting +`x` and `y` into the polar coordinates +""" +@inline function Dϕ_exact(sm::AbstractScatteringModel, λ::Number, x::Number, y::Number) + r = √(x^2 + y^2) + ϕ = atan(y, x) + return quadgk(ϕq -> dDϕ_dz(sm, λ, r, ϕ, ϕq) * sm.Pϕ(ϕq), 0, 2π)[1] +end + + +""" + visibility_point_approx(sm::AbstractScatteringModel, λ::Number, u::Number, v::Number) + +Compute the diffractive kernel for a given observing wavelength `λ` and fourier space coordinates `u`, `v` +using the approximated formula of the phase structure function. +""" +@inline function visibility_point_approx(sm::AbstractScatteringModel, λ::Number, u::Number, v::Number) + b = (u, v) .* (λ / (1 + sm.M)) + return exp(-0.5 * Dϕ_approx(sm, λ, b...)) +end + + +""" + visibility_point_exact(sm::AbstractScatteringModel, λ::Number, u::Number, v::Number) + +Compute the diffractive kernel for a given observing wavelength λ and fourier space coordinates `u`, `v` +using the exact formula of the phase structure function. +""" +@inline function visibility_point_exact(sm::AbstractScatteringModel, λ::Number, u::Number, v::Number) + b = (u, v) .* (λ / (1 + sm.M)) + return exp(-0.5 * Dϕ_exact(sm, λ, b...)) +end diff --git a/src/scatteringmodels/commonfunctions.jl b/src/scatteringmodels/commonfunctions.jl new file mode 100644 index 0000000..18019ce --- /dev/null +++ b/src/scatteringmodels/commonfunctions.jl @@ -0,0 +1,29 @@ +# constants +const FWHMfac = √(2 * log(2)) / π +const c_cgs = 29979245800.0 + +""" +Best-fit parameters of the dipole scattering model derived in Johnson et al. 2018 +""" +const Params_Johnson2018 = (α=1.38, rin_cm=800e5, θmaj_mas=1.380, θmin_mas=0.703, ϕpa_deg=81.9, λ0_cm=1.0, D_pc=2.82, R_pc=5.53) + +# common functions to precomute key constants +@inline calc_A(θmaj, θmin) = θmaj / θmin +@inline calc_ζ0(A) = (A^2 - 1) / (A^2 + 1) +@inline calc_M(D, R) = D / R + +@inline calc_θrad(θmas) = θmas * mas2rad # milliarcseconds to radians +@inline calc_Amaj(rin_cm, λ0_cm, M, θrad) = (rin_cm * (1 + M) * θrad / (FWHMfac * λ0_cm))^2 +@inline calc_Amin = calc_Amaj + +@inline calc_Qbar(α, rin_cm, λ0_cm, M, θmaj_rad, θmin_rad) = 2 / gamma((2 - α) / 2) * (rin_cm^2 * (1 + M) / (FWHMfac * (λ0_cm / 2π)^2))^2 * (θmaj_rad^2 + θmin_rad^2) +@inline calc_C(α, rin_cm, λ0_cm, Qbar) = (λ0_cm / 2π)^2 * Qbar * gamma(1 - α / 2) / (8π^2 * rin_cm^2) + +@inline calc_ϕ0(ϕpa_deg) = deg2rad(90 - ϕpa_deg) + +@inline calc_B_prefac(α, C) = C * 2^(2 - α) * √π / (α * gamma((α + 1) / 2)) +@inline calc_Bmaj(α, ϕ0, Pϕ, B_prefac) = B_prefac * quadgk(ϕ -> abs(cos(ϕ0 - ϕ))^α * Pϕ(ϕ), 0, 2π)[1] +@inline calc_Bmin(α, ϕ0, Pϕ, B_prefac) = B_prefac * quadgk(ϕ -> abs(sin(ϕ0 - ϕ))^α * Pϕ(ϕ), 0, 2π)[1] + +@inline λcm2ν(λcm) = c_cgs / λcm +@inline ν2λcm(ν) = c_cgs / ν diff --git a/src/scatteringmodels/kernels/abstract.jl b/src/scatteringmodels/kernels/abstract.jl new file mode 100644 index 0000000..ccfcca1 --- /dev/null +++ b/src/scatteringmodels/kernels/abstract.jl @@ -0,0 +1,17 @@ +export AbstractScatteringKernel + +""" + $(TYPEDEF) + +An abstract Comrade Sky Model for the diffractive scattering kernel. These are usually +primitive models, and are usually analytic in Fourier but not analytic in the image domain. +As a result a user only needs to implement the following methods. + +- `visibility_point` +- `radialextent` +""" +abstract type AbstractScatteringKernel{T} <: AbstractModel end +@inline flux(::AbstractScatteringKernel{T}) where {T} = one(T) +@inline isprimitive(::Type{<:AbstractScatteringKernel}) = IsPrimitive() +@inline visanalytic(::Type{<:AbstractScatteringKernel}) = IsAnalytic() +@inline imanalytic(::Type{<:AbstractScatteringKernel}) = NotAnalytic() diff --git a/src/scatteringmodels/kernels/approx.jl b/src/scatteringmodels/kernels/approx.jl new file mode 100644 index 0000000..2ade551 --- /dev/null +++ b/src/scatteringmodels/kernels/approx.jl @@ -0,0 +1,33 @@ +export ScatteringKernel +export ApproximatedScatteringKernel + +""" + $(TYPEDEF) + +A Comrade VLBI Sky Model for the scattering kernel based on a Scattering Model +`sm <: AbstractScatteringModel` using the fast approximation formula in Psaltis et al. (2018). + +By default if `T` isn't given, `T` defaults to `Float64` +""" +struct ApproximatedScatteringKernel{T,S,N} <: AbstractScatteringKernel{T} + """Scattering Model""" + sm::S + """Reference Frequency for radial extent. Pick the lowest frequency of your skymodel""" + νref::N + function ApproximatedScatteringKernel{T}(sm::S, νref::N) where {T,S,N} + return new{T,S,N}(sm, νref) + end +end + +ApproximatedScatteringKernel(sm::S, νref::N) where {S,N} = ApproximatedScatteringKernel{Float64}(sm, νref) + +function radialextent(skm::ApproximatedScatteringKernel{T,S,N}) where {T,S,N} + return convert(T, 5) * calc_θrad(skm.sm.θmaj) * (ν2λcm(νref) / skm.sm.λ0)^2 +end + +@inline function visibility_point(skm::ApproximatedScatteringKernel{T,S,N}, u, v, time, freq) where {T,S,N} + return visibility_point_approx(skm.sm, ν2λcm(freq), u, v) + zero(T)im +end + +# use the approximated kernel as the default +ScatteringKernel = ApproximatedScatteringKernel diff --git a/src/scatteringmodels/kernels/exact.jl b/src/scatteringmodels/kernels/exact.jl new file mode 100644 index 0000000..c6887e2 --- /dev/null +++ b/src/scatteringmodels/kernels/exact.jl @@ -0,0 +1,29 @@ +export ExactScatteringKernel + +""" + $(TYPEDEF) + +A Comrade VLBI Sky Model for the scattering kernel based on a Scattering Model +`sm <: AbstractScatteringModel` using the exact formula in Psaltis et al. (2018). + +By default if `T` isn't given, `T` defaults to `Float64` +""" +struct ExactScatteringKernel{T,S,N} <: AbstractScatteringKernel{T} + """Scattering Model""" + sm::S + """Reference frequency in Hz for radial extent. Pick the lowest frequency of your skymodel.""" + νref::N + function ExactScatteringKernel{T}(sm::S, νref::N) where {T,S,N} + return new{T,S,N}(sm, νref) + end +end + +ExactScatteringKernel(sm::S, νref::N) where {S,N} = ExactScatteringKernel{Float64}(sm, νref) + +function radialextent(skm::ExactScatteringKernel{T,S,N}) where {T,S,N} + return convert(T, 5) * calc_θrad(skm.sm.θmaj) * (ν2λcm(νref) / skm.sm.λ0)^2 +end + +@inline function visibility_point(sm::ExactScatteringKernel{T,S,N}, u, v, time, freq) where {T,S,N} + return visibility_point_exact(sm, ν2λcm(freq), u, v) + zero(T)im +end diff --git a/src/scatteringmodels/kernels/kernelmodel.jl b/src/scatteringmodels/kernels/kernelmodel.jl new file mode 100644 index 0000000..45a2839 --- /dev/null +++ b/src/scatteringmodels/kernels/kernelmodel.jl @@ -0,0 +1,23 @@ +export kernelmodel + +""" + kernelmodel(sm::AbstractScatteringModel; νref::Number=c_cgs, use_approx::Bool == true) + +Return a Comrade Sky Model for the diffractive scattering kernel of the input scattering model. + +** Keyword Argurments ** +- `νref::Number`: + the reference frequency in Hz to give a radial extent of the kernel, + which is ideally the lowest frequency of your data sets as the kenerl size is roughly scale with λ^2. + `νref` defaults to the light speed in cgs unit, providing the wavelength of 1 cm. +- `use_approx::Bool==true`: + if `true`, returns a model using the approximated fourmula to compute visibiltiies (`ScatteringKernel`). + Otherwise, return a model instead using the exact formula (`ExactScatteringKernel`). +""" +@inline function kernelmodel(sm::AbstractScatteringModel; νref::Number=c_cgs, use_approx::Bool=true) + if use_approx + return ScatteringKernel(sm, νref) + else + return ExactScatteringKernel(sm, νref) + end +end diff --git a/src/scatteringmodels/models/dipole.jl b/src/scatteringmodels/models/dipole.jl new file mode 100644 index 0000000..86df8eb --- /dev/null +++ b/src/scatteringmodels/models/dipole.jl @@ -0,0 +1,95 @@ +export ScatteringModel +export DipoleScatteringModel + +""" + $(TYPEDEF) + +An anistropic scattering model based on a thin-screen approximation. +This scattering model adopts the dipole field wonder described in Psaltis et al. 2018. + +** Keywords for the constructor ** +The default numbers are based on the best-fit parameters presented in Johnson et al. 2018. +- `α::Number`: The power-law index of the phase fluctuations (Kolmogorov is 5/3). +- `rin_cm::Number`: The inner scale of the scattering screen in cm. +- `θmaj_mas::Number`: FWHM in mas of the major axis angular broadening at the specified reference wavelength. +- `θmin_nas::Number`: FWHM in mas of the minor axis angular broadening at the specified reference wavelength. +- `ϕ_deg::Number`: The position angle of the major axis of the scattering in degree. +- `λ0_cm::Number`: The reference wavelength for the scattering model in cm. +- `D_pc::Number`: The distance from the observer to the scattering screen in pc. +- `R_pc::Number`: The distance from the source to the scattering screen in pc. +""" +struct DipoleScatteringModel{T<:Number,F<:Function} <: AbstractScatteringModel + # Mandatory fields for AbstractScatteringModel + # fundamental parameters + α::T + rin::T + θmaj::T + θmin::T + ϕpa::T + λ0::T + R::T + D::T + + # precomputed constants + M::T + ζ0::T + A::T + kζ::T + Bmaj::T + Bmin::T + Qbar::T + C::T + Amaj::T + Amin::T + ϕ0::T + + # Constructor + function DipoleScatteringModel(; α=1.38, rin_cm=800e5, θmaj_mas=1.380, θmin_mas=0.703, ϕpa_deg=81.9, λ0_cm=1.0, D_pc=2.82, R_pc=5.53) + # compute asymmetry parameters and magnification parameter + A = calc_A(θmaj_mas, θmin_mas) + ζ0 = calc_ζ0(A) + M = calc_M(D_pc, R_pc) + + # position angle (measured from Dec axis in CCW) + # to a more tranditional angle measured from RA axis in CW + ϕ0 = calc_ϕ0(ϕpa_deg) + + # Parameters for the approximate phase structure function + θmaj_rad = calc_θrad(θmaj_mas) # milliarcseconds to radians + θmin_rad = calc_θrad(θmin_mas) # milliarcseconds to radians + Amaj = calc_Amaj(rin_cm, λ0_cm, M, θmaj_rad) + Amin = calc_Amin(rin_cm, λ0_cm, M, θmin_rad) + + # C parameters that scale the powerspectrum of the phase screen + Qbar = calc_Qbar(α, rin_cm, λ0_cm, M, θmaj_rad, θmin_rad) + C = calc_C(α, rin_cm, λ0_cm, Qbar) + + # find kζ + # note: this is depending on the type of the scattering model + kζ = findkzeta_exact(Dipole_KzetaFinder(α, A)) + + # precomputing factor for Pϕ + # note: both lines are depending on the type of the scattering model + Pϕ0 = 1 / (2π * _₂F₁((α + 2) / 2, 0.5, 1, -kζ)) + Pϕfunc(ϕ) = Pϕ(DipoleScatteringModel, ϕ, α, ϕ0, kζ, Pϕ0) + + # B parameters + B_prefac = calc_B_prefac(α, C) + Bmaj = calc_Bmaj(α, ϕ0, Pϕfunc, B_prefac) + Bmin = calc_Bmin(α, ϕ0, Pϕfunc, B_prefac) + + return new{typeof(α),Function}(α, rin_cm, θmaj_mas, θmin_mas, ϕpa_deg, λ0_cm, R_pc, D_pc, M, ζ0, A, kζ, Bmaj, Bmin, Qbar, C, Amaj, Amin, ϕ0) + end +end + + +# the dipole model is the default scattering model +ScatteringModel = DipoleScatteringModel + + +@inline function Pϕ(::Type{<:DipoleScatteringModel}, ϕ, α, ϕ0, kζ, Pϕ0) + return Pϕ0 * (1 + kζ * sin(ϕ - ϕ0)^2)^(-(α + 2) / 2) +end + + +@inline Pϕ(sm::DipoleScatteringModel, ϕ) = Pϕ(DipoleScatteringModel, ϕ, sm.α, sm.ϕ0, sm.kζ, sm.Pϕ0) diff --git a/src/scatteringmodels/models/periodicboxcar.jl b/src/scatteringmodels/models/periodicboxcar.jl new file mode 100644 index 0000000..164499b --- /dev/null +++ b/src/scatteringmodels/models/periodicboxcar.jl @@ -0,0 +1,92 @@ +export PeriodicBoxCarScatteringModel + +""" + $(TYPEDEF) + +An anistropic scattering model based on a thin-screen approximation. +This scattering adopts the periodic boxcar field wonder described in Psaltis et al. 2018. + +** Keywords for the constructor ** +The default numbers are based on the best-fit parameters presented in Johnson et al. 2018. +- `α::Number`: The power-law index of the phase fluctuations (Kolmogorov is 5/3). +- `rin_cm::Number`: The inner scale of the scattering screen in cm. +- `θmaj_mas::Number`: FWHM in mas of the major axis angular broadening at the specified reference wavelength. +- `θmin_nas::Number`: FWHM in mas of the minor axis angular broadening at the specified reference wavelength. +- `ϕ_deg::Number`: The position angle of the major axis of the scattering in degree. +- `λ0_cm::Number`: The reference wavelength for the scattering model in cm. +- `D_pc::Number`: The distance from the observer to the scattering screen in pc. +- `R_pc::Number`: The distance from the source to the scattering screen in pc. +""" +struct PeriodicBoxCarScatteringModel{T<:Number,F<:Function} <: AbstractScatteringModel + # Mandatory fields for AbstractScatteringModel + # fundamental parameters + α::T + rin::T + θmaj::T + θmin::T + ϕpa::T + λ0::T + R::T + D::T + # precomputed constants + M::T + ζ0::T + A::T + kζ::T + Bmaj::T + Bmin::T + Qbar::T + C::T + Amaj::T + Amin::T + ϕ0::T + + function PeriodicBoxCarScatteringModel(; α=1.38, rin_cm=800e5, θmaj_mas=1.380, θmin_mas=0.703, ϕpa_deg=81.9, λ0_cm=1.0, D_pc=2.82, R_pc=5.53) + # compute asymmetry parameters and magnification parameter + A = calc_A(θmaj_mas, θmin_mas) + ζ0 = calc_ζ0(A) + M = calc_M(D_pc, R_pc) + + # Parameters for the approximate phase structure function + θmaj_rad = calc_θrad(θmaj_mas) # milliarcseconds to radians + θmin_rad = calc_θrad(θmin_mas) # milliarcseconds to radians + Amaj = calc_Amaj(rin_cm, λ0_cm, M, θmaj_rad) + Amin = calc_Amin(rin_cm, λ0_cm, M, θmin_rad) + + # C parameters that scale the powerspectrum of the phase screen + Qbar = calc_Qbar(α, rin_cm, λ0_cm, M, θmaj_rad, θmin_rad) + C = calc_C(α, rin_cm, λ0_cm, Qbar) + + # position angle (measured from Dec axis in CCW) to a more tranditional angle measured from RA axis in CW + ϕ0 = calc_ϕ0(ϕpa_deg) + + # find kζ + # note: this is depending on the type of the scattering model + kζ = findkzeta_exact(PeriodicBoxCar_KzetaFinder(ζ0)) + + # precomputing factor for Pϕ + # note: both lines are depending on the type of the scattering model + Pϕ0 = (1 + kζ) / 2π + Pϕfunc(ϕ) = Pϕ(PeriodicBoxCarScatteringModel, ϕ, ϕ0, kζ, Pϕ0) + + # B parameters + B_prefac = calc_B_prefac(α, C) + Bmaj = calc_Bmaj(α, ϕ0, Pϕfunc, B_prefac) + Bmin = calc_Bmin(α, ϕ0, Pϕfunc, B_prefac) + + return new{typeof(α),Function}(α, rin_cm, θmaj_mas, θmin_mas, ϕpa_deg, λ0_cm, M, R, D, ζ0, A, kζ, Bmaj, Bmin, Qbar, C, Amaj, Amin, ϕ0) + end +end + +@inline function Pϕ(::Type{<:PeriodicBoxCarScatteringModel}, ϕ, ϕ0, kζ, Pϕ0) + # set up a range for the box car function + Δϕ̃0 = (ϕ - ϕ0) % π + # orig: (np.pi/(2.0*(1.0 + self.kzeta)) < (phi - self.phi0) % np.pi) + condition1 = π / (2 * (1 + kζ)) < Δϕ̃0 + # orig: (phi - self.phi0) % np.pi < np.pi - np.pi/(2.0*(1.0 + self.kzeta)) + condition2 = Δϕ̃0 < π * (1 - 0.5 / (1 + kζ)) + + return Pϕ0 * (1 - condition1 & condition2) +end + +@inline Pϕ(sm::PeriodicBoxCarScatteringModel, ϕ) = Pϕ(PeriodicBoxCarScatteringModel, ϕ, sm.ϕ0, sm.kζ, sm.Pϕ0) diff --git a/src/scatteringmodels/models/vonmises.jl b/src/scatteringmodels/models/vonmises.jl new file mode 100644 index 0000000..60a01c6 --- /dev/null +++ b/src/scatteringmodels/models/vonmises.jl @@ -0,0 +1,85 @@ +export vonMisesScatteringModel + +""" + $(TYPEDEF) + +An anistropic scattering model based on a thin-screen approximation. +This scattering adopts the von Mises field wonder described in Psaltis et al. 2018. + +** Keywords for the constructor ** +The default numbers are based on the best-fit parameters presented in Johnson et al. 2018. +- `α::Number`: The power-law index of the phase fluctuations (Kolmogorov is 5/3). +- `rin_cm::Number`: The inner scale of the scattering screen in cm. +- `θmaj_mas::Number`: FWHM in mas of the major axis angular broadening at the specified reference wavelength. +- `θmin_nas::Number`: FWHM in mas of the minor axis angular broadening at the specified reference wavelength. +- `ϕ_deg::Number`: The position angle of the major axis of the scattering in degree. +- `λ0_cm::Number`: The reference wavelength for the scattering model in cm. +- `D_pc::Number`: The distance from the observer to the scattering screen in pc. +- `R_pc::Number`: The distance from the source to the scattering screen in pc. +""" +struct vonMisesScatteringModel{T<:Number,F<:Function} <: AbstractScatteringModel + # Mandatory fields for AbstractScatteringModel + # fundamental parameters + α::T + rin::T + θmaj::T + θmin::T + ϕpa::T + λ0::T + R::T + D::T + # precomputed constants + M::T + ζ0::T + A::T + kζ::T + Bmaj::T + Bmin::T + Qbar::T + C::T + Amaj::T + Amin::T + ϕ0::T + + function vonMisesScatteringModel(; α=1.38, rin_cm=800e5, θmaj_mas=1.380, θmin_mas=0.703, ϕpa_deg=81.9, λ0_cm=1.0, D_pc=2.82, R_pc=5.53) + # compute asymmetry parameters and magnification parameter + A = calc_A(θmaj_mas, θmin_mas) + ζ0 = calc_ζ0(A) + M = calc_M(D_pc, R_pc) + + # Parameters for the approximate phase structure function + θmaj_rad = calc_θrad(θmaj_mas) # milliarcseconds to radians + θmin_rad = calc_θrad(θmin_mas) # milliarcseconds to radians + Amaj = calc_Amaj(rin_cm, λ0_cm, M, θmaj_rad) + Amin = calc_Amin(rin_cm, λ0_cm, M, θmin_rad) + + # C parameters that scale the powerspectrum of the phase screen + Qbar = calc_Qbar(α, rin_cm, λ0_cm, M, θmaj_rad, θmin_rad) + C = calc_C(α, rin_cm, λ0_cm, Qbar) + + # position angle (measured from Dec axis in CCW) to a more tranditional angle measured from RA axis in CW + ϕ0 = calc_ϕ0(ϕpa_deg) + + # find kζ + # note: this is depending on the type of the scattering model + kζ = findkzeta_exact(vonMises_KzetaFinder(A)) + + # precomputing factor for Pϕ + # note: both lines are depending on the type of the scattering model + Pϕ0 = 1 / (2π * besseli0(kζ)) + Pϕfunc(ϕ) = Pϕ(vonMisesScatteringModel, ϕ, ϕ0, kζ, Pϕ0) + + # B parameters + B_prefac = calc_B_prefac(α, C) + Bmaj = calc_Bmaj(α, ϕ0, Pϕfunc, B_prefac) + Bmin = calc_Bmin(α, ϕ0, Pϕfunc, B_prefac) + + return new{typeof(α),Function}(α, rin_cm, θmaj_mas, θmin_mas, ϕpa_deg, λ0_cm, M, R, D, ζ0, A, kζ, Bmaj, Bmin, Qbar, C, Amaj, Amin, ϕ0) + end +end + +@inline function Pϕ(::Type{<:vonMisesScatteringModel}, ϕ, ϕ0, kζ, Pϕ0) + return Pϕ0 * cosh(kζ * cos(ϕ - ϕ0)) +end + +@inline Pϕ(sm::vonMisesScatteringModel, ϕ) = Pϕ(vonMisesScatteringModel, ϕ, sm.ϕ0, sm.kζ, sm.Pϕ0) From cdfcefb681ebae05d815db0e0055f3f9eb2aeede Mon Sep 17 00:00:00 2001 From: Kazunori Akiyama Date: Thu, 27 Jul 2023 04:05:50 -0400 Subject: [PATCH 19/19] typo correction --- src/scatteringmodels/kernels/exact.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/scatteringmodels/kernels/exact.jl b/src/scatteringmodels/kernels/exact.jl index c6887e2..1d12e1e 100644 --- a/src/scatteringmodels/kernels/exact.jl +++ b/src/scatteringmodels/kernels/exact.jl @@ -24,6 +24,6 @@ function radialextent(skm::ExactScatteringKernel{T,S,N}) where {T,S,N} return convert(T, 5) * calc_θrad(skm.sm.θmaj) * (ν2λcm(νref) / skm.sm.λ0)^2 end -@inline function visibility_point(sm::ExactScatteringKernel{T,S,N}, u, v, time, freq) where {T,S,N} - return visibility_point_exact(sm, ν2λcm(freq), u, v) + zero(T)im +@inline function visibility_point(skm::ExactScatteringKernel{T,S,N}, u, v, time, freq) where {T,S,N} + return visibility_point_exact(skm.sm, ν2λcm(freq), u, v) + zero(T)im end