diff --git a/Project.toml b/Project.toml index 2357b96..b8ff540 100644 --- a/Project.toml +++ b/Project.toml @@ -1,11 +1,12 @@ name = "PolarizedTypes" uuid = "d3c5d4cd-a8ee-40d6-aac7-e34df5a20044" authors = ["Paul Tiede and contributors"] -version = "0.1.1" +version = "0.1.2" [deps] ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" +LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" diff --git a/src/PolarizedTypes.jl b/src/PolarizedTypes.jl index bc254c9..78ad510 100644 --- a/src/PolarizedTypes.jl +++ b/src/PolarizedTypes.jl @@ -4,12 +4,14 @@ using ChainRulesCore using DocStringExtensions using PrecompileTools using StaticArrays +using LinearAlgebra # Write your package code here. -export StokesParams, CoherencyMatrix, SingleStokes, RPol, LPol, XPol, YPol, +export StokesParams, CoherencyMatrix, RPol, LPol, XPol, YPol, PolBasis, CirBasis, LinBasis, basis_components, basis_transform, coherencymatrix, stokesparams, - linearpol, mbreve, m̆, evpa + linearpol, mbreve, m̆, evpa, polarization, fracpolarization, mpol, + polellipse include("types.jl") include("basis_transforms.jl") diff --git a/src/basis_transforms.jl b/src/basis_transforms.jl index 9ad910b..391ce58 100644 --- a/src/basis_transforms.jl +++ b/src/basis_transforms.jl @@ -97,3 +97,9 @@ function basis_transform end @inline basis_transform(::Type{T}, p::Pair{B1,B2}) where {T, B1<:PolBasis, B2<:PolBasis} = basis_transform(T, B1(), B2()) @inline basis_transform(::Pair{B1,B2}) where {B1<:PolBasis, B2<:PolBasis} = basis_transform(Float64, B1(), B2()) + +@inline function basis_transform(c::CoherencyMatrix{B1, B2, Complex{T}}, e1::PolBasis, e2::PolBasis) where {B1, B2, T} + t1 = basis_transform(T, B1()=>e1) + t2 = basis_transform(T, e2=>B2()) + return CoherencyMatrix(t1*c*t2, e1, e2) +end diff --git a/src/functions.jl b/src/functions.jl index c612636..923ec3f 100644 --- a/src/functions.jl +++ b/src/functions.jl @@ -9,18 +9,99 @@ end """ $(SIGNATURES) -Compute the fractional linear polarization of a stokes vector -or coherency matrix + +Returns the (Q, U, V) polarization vector as a 3-element static vector. +""" +function polarization(s::StokesParams) + return SVector(s.Q, s.U, s.V) +end + +""" + $(SIGNATURES) + +Returns the polarization ellipse of the Stokes parameters `s`. The results +is a named tuple with elements + - `a` : The semi-major axis of the polarization ellipse + - `b` : The semi-minor axis of the polarization ellipse + - `evpa`: The electric vector position angle of `s` or the PA of the ellipse. + - `sn` : The sign of the Stokes `V`. + + +## Notes + +The semi-major and semi-minor axes are defined as +``` + a = 1/2(Iₚ + |L|) + b = 1/2(Iₚ - |L|) +``` +where `Iₚ = √(Q² + U² + V²)` and `|L| = √(Q² + U²)`. + +In general the area of the ellipse is given by +``` + πab = π/4|V|² +``` + +For sources with zero linear polarization `a = b` so we +have a circle with radius `|V|`. For purely linear polarization `b = 0` giving a line +with length `|L|`. + +""" +function polellipse(s::StokesParams{<:Real}) + l = linearpol(s) + p = norm(polarization(s)) + a = (p + abs(l))/2 + b = (p - abs(l))/2 + ev = evpa(s) + sn = sign(s.V) + return (;a, b, evpa=ev, sn) +end + + + +""" + $(SIGNATURES) + +Returns the (Q/I, U/I, V/I) fractional polarization vector as a 3-element static vector. +""" +fracpolarization(s::StokesParams{T}) where {T} = polarization(s)*inv(s.I + eps(T)) +fracpolarization(s::StokesParams{Complex{T}}) where {T} = polarization(s)*inv(s.I + eps(T)) + +""" + mpol(m::StokesParameters{<:Real}) + + +Compute the complex fractional linear polarization of a Stokes Parameter `m` +""" +mpol(m::StokesParams{T}) where {T<:Real} = (m.Q + 1im*m.U)/(m.I + eps(T)) + +""" + m̆(m::Union{StokesParameters{<:Complex}, CoherencyMatrix) + +Computes the complex fractional linear polarization of the complex or visibility quantities. +Note that this function can also be called used [`mbreve`](@ref) """ -m̆(m::StokesParams{T}) where {T} = (m.Q + 1im*m.U)/(m.I + eps(T)) m̆(m::StokesParams{Complex{T}}) where {T} = (m.Q + 1im*m.U)/(m.I + eps(T)) -m̆(m::CoherencyMatrix{CirBasis,CirBasis}) = 2*m.e12/(m.e11+m.e22) -mbreve(m::Union{StokesParams, CoherencyMatrix}) = m̆(m) + """ $(SIGNATURES) -Compute the evpa of a stokes vector or cohereny matrix. + +Computes the complex fractional linear polarization of the complex or visibility quantities. +Note that this function can also be called used [`m̆`](@ref) +""" +mbreve(m::StokesParams) = m̆(m) + +""" + evpa(m::Union{StokesParams, CoherencyMatrix}) +Compute the evpa of a stokes vect or cohereny matrix. """ evpa(m::StokesParams) = 1/2*atan(m.U, m.Q) evpa(m::StokesParams{<:Complex}) = 1/2*angle(m.U/m.Q) -evpa(m::CoherencyMatrix{CirBasis, CirBasis}) = angle(m.e12) + + +# Now define the functions for Coherency matrix +for f in (:linearpol, :polarization, :fracpolarization, :m̆, :mbreve, :evpa) + @eval begin + $(f)(c::CoherencyMatrix) = $(f)(StokesParams(c)) + end +end diff --git a/src/types.jl b/src/types.jl index e668d92..e91537f 100644 --- a/src/types.jl +++ b/src/types.jl @@ -87,6 +87,8 @@ struct StokesParams{T} <: FieldVector{4,T} V::T end +StokesParams(::SArray) = throw(ArgumentError("argument does not have a basis please wrap it in a `CoherencyMatrix`")) + StaticArraysCore.similar_type(::Type{StokesParams}, ::Type{T}, s::Size{(4,)}) where {T} = StokesParams{T} @@ -283,6 +285,10 @@ end return CoherencyMatrix(RR, LR, RL, LL, CirBasis(), CirBasis()) end +@inline function CoherencyMatrix{B1, B2}(s::StokesParams) where {B1, B2} + return CoherencyMatrix(s, B1(), B2()) +end + @inline function CoherencyMatrix{LinBasis, LinBasis}(s::StokesParams) (;I,Q,U,V) = s diff --git a/test/runtests.jl b/test/runtests.jl index 473c5e7..ed9dd26 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -92,12 +92,31 @@ using Test @test s ≈ StokesParams(CoherencyMatrix(s, PolBasis{YPol,XPol}(), PolBasis{LPol,RPol}())) end + @testset "Mixed Pol" begin + I = 2.0 + 0.5im + Q = rand(ComplexF64) - 0.5 + U = rand(ComplexF64) - 0.5 + V = rand(ComplexF64) - 0.5 + s = StokesParams(I, Q, U, V) + + c1 = CoherencyMatrix(s, CirBasis(), LinBasis()) + c2 = CoherencyMatrix{CirBasis, LinBasis}(s) + c3 = basis_transform(CoherencyMatrix(s, CirBasis()), CirBasis(), LinBasis()) + @test c1 ≈ c2 + @test c1 ≈ c3 + + @test StokesParams(c1) ≈ s + @test StokesParams(c2) ≈ s + @test StokesParams(c3) ≈ s + end + @testset "Conversion Consistency" begin s = StokesParams(1.0 .+ 0.0im, 0.2 + 0.2im, 0.2 - 0.2im, 0.1+0.05im) c_lin1 = CoherencyMatrix(s, LinBasis()) c_lin2 = CoherencyMatrix(s, PolBasis{XPol,YPol}()) c_lin3 = CoherencyMatrix(s, PolBasis{XPol,YPol}(), PolBasis{XPol,YPol}()) + @test c_lin1 ≈ c_lin2 ≈ c_lin3 c_cir1 = CoherencyMatrix(s, CirBasis()) @@ -112,6 +131,8 @@ using Test @test t2*c_cir1*t1 ≈ c_lin1 @test t1*c_lin1*t2 ≈ c_cir1 + @test_throws ArgumentError StokesParams(t1*c_lin1*t2) + # Test the mixed basis @test c_cir1*t1 ≈ t1*c_lin1 @test c_lin1*t2 ≈ t2*c_cir1 @@ -125,4 +146,48 @@ using Test @test_opt StokesParams(CoherencyMatrix(s, LinBasis(), CirBasis())) @test_opt StokesParams(CoherencyMatrix(s, LinBasis(), CirBasis(), LinBasis())) end + + @testset "Polarized Functions" begin + s = StokesParams(2.0, 0.25, -0.25, 0.25) + @test linearpol(s) == complex(0.25, -0.25) + @test evpa(s) ≈ atan(-0.5, 0.5)/2 + @test s[2:end] ≈ polarization(s) + + + slin = StokesParams(2.0, 0.2, 0.2, 0.0) + fp = fracpolarization(slin) + @test complex(fp[1], fp[2]) ≈ linearpol(slin)/slin.I + @test fp[end] ≈ 0 + + @testset "ellipse" begin + p = polellipse(s) + @test p.a*p.b ≈ s.V^2/4 + @test p.evpa ≈ evpa(s) + plin = polellipse(slin) + @test plin.a ≈ (abs(linearpol(slin))) + @test isapprox(plin.b, 0, atol=1e-8) + @test plin.evpa ≈ evpa(slin) + @test p.sn ≈ sign(s.V) + end + + + @test mpol(s) ≈ complex(0.25, -0.25)/2 + + @testset "Complex Vis" begin + I = 2.0 + 0.5im + Q = rand(ComplexF64) - 0.5 + U = rand(ComplexF64) - 0.5 + V = rand(ComplexF64) - 0.5 + s = StokesParams(I, Q, U, V) + @test mbreve(s) == m̆(s) + c = CoherencyMatrix{CirBasis, LinBasis}(s) + @test linearpol(c) ≈ linearpol(s) + @test polarization(c) ≈ polarization(s) + @test fracpolarization(c) ≈ fracpolarization(s) + @test m̆(c) ≈ m̆(s) + @test mbreve(c) ≈ mbreve(s) + @test evpa(c) ≈ evpa(s) + end + + end end