Merge branch 'main' into compathelper/new_version/2023-07-11-19-07-14…

…-916-01957477830

Project.toml

@@ -1,16 +1,25 @@
 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]
+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]
+DocStringExtensions = "0.9"
 HypergeometricFunctions = "0.3"
+NonlinearSolve = "1"
+SpecialFunctions = "2"
 julia = "1.9"
 
 [extras]

src/ScatteringOptics.jl

@@ -1,17 +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 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

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

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

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

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