Merge pull request #71 from EHTJulia/ptiede-polupdates

Loop for pol transform

Author: ptiede

ext/VLBISkyModelsMakieExt.jl

@@ -105,7 +105,7 @@ end
                 colormap = :grayC,
                 colorrange = Makie.automatic,
                 pcolorrange=Makie.automatic,
-                pcolormap=:coolwarm,
+                pcolormap=Makie.Reverse(:RdBu),
                 nvec = 30,
                 min_frac = 0.1,
                 min_pol_frac=0.2,
@@ -337,7 +337,7 @@ function Makie.plot!(plot::PolImage)
     pcm = lift(plot.pcolormap, plot.plot_total) do pc, pt
         (pc != Makie.automatic) && return pc
         if pt
-            return :coolwarm
+            return Makie.Reverse(:RdBu)
         else
             return :rainbow1
         end
@@ -408,7 +408,7 @@ function _imgviz!(fig, ax, img::IntensityMap{<:Real}; scale_length=fieldofview(i
     dkwargs = Dict(kwargs)
     crange = get(dkwargs, :colorrange, colorrange_default)
     delete!(dkwargs, :colorrange)
-    cmap = get(dkwargs, :colormap, :magma)
+    cmap = get(dkwargs, :colormap, :inferno)
     delete!(dkwargs, :colormap)
 
     hm = heatmap!(ax, img; colorrange=crange, colormap=cmap, dkwargs...)

src/models/polarized.jl

@@ -255,12 +255,22 @@ The arguments are:
 
 """
 function PoincareSphere2Map(I, p, X, grid)
-    pimgI = I .* p
-    stokesI = I
-    stokesQ = pimgI .* X[1]
-    stokesU = pimgI .* X[2]
-    stokesV = pimgI .* X[3]
-    return stokes_intensitymap(stokesI, stokesQ, stokesU, stokesV, grid)
+    Q = similar(I)
+    U = similar(I)
+    V = similar(I)
+
+    t1 = X[1]
+    t2 = X[2]
+    t3 = X[3]
+
+    @inbounds @simd for i in eachindex(I, Q, U, V)
+        pI = I[i] * p[i]
+        Q[i] = pI * t1[i]
+        U[i] = pI * t2[i]
+        V[i] = pI * t3[i]
+    end
+
+    return stokes_intensitymap(I, Q, U, V, grid)
 end
 
 # function PoincareSphere2Map(I, p, X, grid)
@@ -291,20 +301,39 @@ Each Stokes parameter is parameterized as
 
 where `a,b,c,d` are real numbers with no conditions, and `p=√(a² + b² + c²)`.
 """
-function PolExp2Map(a::AbstractArray, b::AbstractArray, c::AbstractArray, d::AbstractArray,
-                    grid::AbstractRectiGrid)
-    p = sqrt.(b .^ 2 .+ c .^ 2 .+ d .^ 2)
-
-    tmp = exp.(a) .* sinh.(p) .* inv.(p)
-    pimgI = exp.(a) .* cosh.(p)
-    pimgQ = tmp .* b
-    pimgU = tmp .* c
-    pimgV = tmp .* d
-    stokesI = pimgI
-    stokesQ = pimgQ
-    stokesU = pimgU
-    stokesV = pimgV
-    return stokes_intensitymap(stokesI, stokesQ, stokesU, stokesV, grid)
+@fastmath function PolExp2Map(a::AbstractArray,
+                              b::AbstractArray,
+                              c::AbstractArray,
+                              d::AbstractArray,
+                              grid::AbstractRectiGrid)
+    pimgI = similar(a)
+    pimgQ = similar(b)
+    pimgU = similar(c)
+    pimgV = similar(d)
+
+    # This is just faster because it is a 1 pass algorithm
+    @inbounds @simd for i in eachindex(pimgI, pimgQ, pimgU, pimgV)
+        p = sqrt(b[i]^2 + c[i]^2 + d[i]^2)
+        ea = exp(a[i])
+        tmp = ea * sinh(p) / p
+        pimgI[i] = ea * cosh(p)
+        pimgQ[i] = tmp * b[i]
+        pimgU[i] = tmp * c[i]
+        pimgV[i] = tmp * d[i]
+    end
+
+    # p = @.. sqrt(b ^ 2 + c ^ 2 + d ^ 2)
+    # ea = @.. exp(a)
+    # tmp = @.. ea * sinh(p) * inv(p)
+    # pimgI = @.. ea * cosh(p)
+    # pimgQ = @.. tmp * b
+    # pimgU = @.. tmp * c
+    # pimgV = @.. tmp * d
+    # stokesI = pimgI
+    # stokesQ = pimgQ
+    # stokesU = pimgU
+    # stokesV = pimgV
+    return stokes_intensitymap(pimgI, pimgQ, pimgU, pimgV, grid)
 end
 
 """