Add more polarized functions

Project.toml

@@ -1,11 +1,12 @@
 name = "PolarizedTypes"
 uuid = "d3c5d4cd-a8ee-40d6-aac7-e34df5a20044"
 authors = ["Paul Tiede <[email protected]> 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"
 

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")

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

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

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

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