Polarized Image and Instrumental Modeling

In this tutorial, we will analyze a simulated simple polarized dataset to demonstrate Comrade's polarized imaging capabilities.

Introduction to Polarized Imaging

The EHT is a polarized interferometer. However, like all VLBI interferometers, it does not directly measure the Stokes parameters (I, Q, U, V). Instead, it measures components related to the electric field at the telescope along two directions using feeds. There are two types of feeds at telescopes: circular, which measure $R / L$ components of the electric field, and linear feeds, which measure $X / Y$ components of the electric field. Most sites in the EHT use circular feeds, meaning they measure the right (R) and left electric field (L) at each telescope. Although note that ALMA actually uses linear feeds. Currently Comrade has the ability to fit natively mixed polarization data however, the publically released EHT data has been converted to circular polarization. For a VLBI array whose feeds are purely circluar the coherency matrices are given by,

C_{i j} = (\begin{matrix} R R^{*} & R L^{*} \\ L R^{*} & L L^{*} \end{matrix}) .

These coherency matrices are the fundamental object in interferometry and what the telescope observes. For a perfect interferometer, these circular coherency matrices are related to the usual Fourier transform of the stokes parameters by

(\begin{matrix} \tilde{I} \\ \tilde{Q} \\ \tilde{U} \\ \tilde{V} \end{matrix}) = \frac{1}{2} (\begin{matrix} R R^{*} + L L^{*} \\ R L^{*} + L R^{*} \\ i (L R^{*} - R L^{*}) \\ R R^{*} - L L^{*} \end{matrix}) .

Note

In this tutorial, we stick to circular feeds but Comrade has the capabilities to model linear (XX,XY, ...) and mixed basis coherencies (e.g., RX, RY, ...).

In reality, the measure coherencies are corrupted by both the atmosphere and the telescope itself. In Comrade we use the RIME formalism ^[1] to represent these corruptions, namely our measured coherency matrices $V_{i j}$ are given by

V_{i j} = J_{i} C_{i j} J_{j}^{†}

where $J$ is known as a Jones matrix and $i j$ denotes the baseline $i j$ with sites $i$ and $j$ .

Comrade is highly flexible with how the Jones matrices are formed and provides several convenience functions that parameterize standard Jones matrices. These matrices include:

JonesG which builds the set of complex gain Jones matrices

G = (\begin{matrix} g_{a} & 0 \\ 0 & g_{b} \end{matrix})

JonesD which builds the set of complex d-terms Jones matrices

D = (\begin{matrix} 1 & d_{a} \\ d_{b} & 1 \end{matrix})

JonesR is the basis transform matrix $T$ . This transformation is special and combines two things using the decomposition $T = F B$ . The first, $B$ , is the transformation from some reference basis to the observed coherency basis (this allows for mixed basis measurements). The second is the feed rotation, $F$ , that transforms from some reference axis to the axis of the telescope as the source moves in the sky. The feed rotation matrix F for circular feeds in terms of the per station feed rotation angle $φ$ is

F = (\begin{matrix} e^{- i φ} & 0 \\ 0 & e^{i φ} \end{matrix})

In the rest of the tutorial, we are going to solve for all of these instrument model terms on in addition to our image structure to reconstruct a polarized image of a synthetic dataset.

Load the Data

To get started we will load Comrade

julia

using Comrade

Load the Data

julia

using Pyehtim

For reproducibility we use a stable random number genreator

julia

using StableRNGs
rng = StableRNG(42)

StableRNGs.LehmerRNG(state=0x00000000000000000000000000000055)

Now we will load some synthetic polarized data.

julia

obs = Pyehtim.load_uvfits_and_array(
        joinpath(__DIR, "..", "..", "Data", "polarized_gaussian_all_corruptions.uvfits"),
        joinpath(__DIR, "..", "..", "Data", "array.txt"), polrep="circ")

Python: <ehtim.obsdata.Obsdata object at 0x7f55d45036a0>

Notice that, unlike other non-polarized tutorials, we need to include a second argument. This is the array file of the observation and is required to determine the feed rotation of the array.

Now we scan average the data since the data to boost the SNR and reduce the total data volume.

julia

obs = scan_average(obs)

Python: <ehtim.obsdata.Obsdata object at 0x7f55f6318a00>

Now we extract our observed/corrupted coherency matrices.

julia

dvis = extract_table(obs, Coherencies())

EHTObservationTable{Comrade.EHTCoherencyDatum{Float64}}
  source:      17.761122472222223:-28.992189444444445
  mjd:         51544
  bandwidth:   1.0e9
  sites:       [:AA, :AP, :AZ, :JC, :LM, :PV, :SM]
  nsamples:    315

##Building the Model/Posterior

To build the model, we first break it down into two parts:

The image or sky model. In Comrade, all polarized image models are written in terms of the Stokes parameters. In this tutorial, we will use a polarized image model based on Pesce (2021)^[2], and parameterizes each pixel in terms of the Poincare sphere. This parameterization ensures that we have all the physical properties of stokes parameters. Note that we also have a parameterization in terms of hyperbolic trig functions VLBISkyModels.PolExp2Map
The instrument model. The instrument model specifies the model that describes the impact of instrumental and atmospheric effects. We will be using the $J = G D R$ decomposition we described above. However, to parameterize the R/L complex gains, we will be using a gain product and ratio decomposition. The reason for this decomposition is that in realistic measurements, the gain ratios and products have different temporal characteristics. Namely, many of the EHT observations tend to demonstrate constant R/L gain ratios across an nights observations, compared to the gain products, which vary every scan. Additionally, the gain ratios tend to be smaller (i.e., closer to unity) than the gain products. Using this apriori knowledge, we can build this into our model and reduce the total number of parameters we need to model.

First we specify our sky model. As always Comrade requires this to be a two argument function where the first argument is typically a NamedTuple of parameters we will fit and the second are additional metadata required to build the model.

julia

using StatsFuns: logistic
function sky(θ, metadata)
    (;c, σ, p, p0, pσ, angparams) = θ
    (;ftot, grid) = metadata
    # Build the stokes I model
    rast = ftot*to_simplex(CenteredLR(), σ.*c.params)
    # The total polarization fraction is modeled in logit space so we transform it back
    pim = logistic.(p0 .+ pσ.*p.params)
    # Build our IntensityMap
    pmap = PoincareSphere2Map(rast, pim, angparams, grid)
    # Construct the actual image model which uses a third order B-spline pulse
    m = ContinuousImage(pmap, BSplinePulse{3}())
    # Finally find the image centroid and shift it to be at the center
    x0, y0 = centroid(pmap)
    ms = shifted(m, -x0, -y0)
    return ms
end

sky (generic function with 1 method)

Now, we define the model metadata required to build the model. We specify our image grid and cache model needed to define the polarimetric image model. Our image will be a 10x10 raster with a 60μas FOV.

julia

using Distributions, DistributionsAD
using VLBIImagePriors
fovx = μas2rad(60.0)
fovy = μas2rad(60.0)
nx = ny = 10
grid = imagepixels(fovx, fovy, nx, ny)

RectiGrid(
executor: Serial()
Dimensions: 
↓ X Sampled{Float64} LinRange{Float64}(-1.3089969389957472e-10, 1.3089969389957472e-10, 10) ForwardOrdered Regular Points,
→ Y Sampled{Float64} LinRange{Float64}(-1.3089969389957472e-10, 1.3089969389957472e-10, 10) ForwardOrdered Regular Points
)

For the image metadata we specify the grid and the total flux of the image, which is 1.0. Note that we specify the total flux out front since it is degenerate with an overall shift in the gain amplitudes.

julia

skymeta = (;ftot=1.0, grid)

(ftot = 1.0, grid = RectiGrid{Tuple{X{DimensionalData.Dimensions.Lookups.Sampled{Float64, LinRange{Float64, Int64}, DimensionalData.Dimensions.Lookups.ForwardOrdered, DimensionalData.Dimensions.Lookups.Regular{Float64}, DimensionalData.Dimensions.Lookups.Points, DimensionalData.Dimensions.Lookups.NoMetadata}}, Y{DimensionalData.Dimensions.Lookups.Sampled{Float64, LinRange{Float64, Int64}, DimensionalData.Dimensions.Lookups.ForwardOrdered, DimensionalData.Dimensions.Lookups.Regular{Float64}, DimensionalData.Dimensions.Lookups.Points, DimensionalData.Dimensions.Lookups.NoMetadata}}}, Serial, ComradeBase.NoHeader}((X{DimensionalData.Dimensions.Lookups.Sampled{Float64, LinRange{Float64, Int64}, DimensionalData.Dimensions.Lookups.ForwardOrdered, DimensionalData.Dimensions.Lookups.Regular{Float64}, DimensionalData.Dimensions.Lookups.Points, DimensionalData.Dimensions.Lookups.NoMetadata}}([-1.3089969389957472e-10, -1.0181087303300256e-10, -7.27220521664304e-11, -4.363323129985825e-11, -1.454441043328608e-11, 1.454441043328608e-11, 4.363323129985823e-11, 7.27220521664304e-11, 1.0181087303300256e-10, 1.3089969389957472e-10]), Y{DimensionalData.Dimensions.Lookups.Sampled{Float64, LinRange{Float64, Int64}, DimensionalData.Dimensions.Lookups.ForwardOrdered, DimensionalData.Dimensions.Lookups.Regular{Float64}, DimensionalData.Dimensions.Lookups.Points, DimensionalData.Dimensions.Lookups.NoMetadata}}([-1.3089969389957472e-10, -1.0181087303300256e-10, -7.27220521664304e-11, -4.363323129985825e-11, -1.454441043328608e-11, 1.454441043328608e-11, 4.363323129985823e-11, 7.27220521664304e-11, 1.0181087303300256e-10, 1.3089969389957472e-10])), Serial(), ComradeBase.NoHeader()))

For our image prior we will use a simpler prior than

We use again use a GMRF prior. For more information see the Imaging a Black Hole using only Closure Quantities tutorial.
For the total polarization fraction, p, we assume an uncorrelated uniform prior ImageUniform for each pixel.
To specify the orientation of the polarization, angparams, on the Poincare sphere, we use a uniform spherical distribution, ImageSphericalUniform.

julia

rat = beamsize(dvis)/step(grid.X)
cmarkov = ConditionalMarkov(GMRF, grid; order=1)
dρ = truncated(InverseGamma(1.0, -log(0.1)*rat); lower=2.0, upper=max(nx, ny))
cprior = HierarchicalPrior(cmarkov, dρ)

HierarchicalPrior(
	map: 
	ConditionalMarkov(
Random Field: GaussMarkovRandomField
Graph: MarkovRandomFieldGraph{1}(
dims: (10, 10)
)
)	hyper prior: 
	Truncated(Distributions.InverseGamma{Float64}(
invd: Distributions.Gamma{Float64}(α=1.0, θ=0.09157612316984907)
θ: 10.919876987424647
)
; lower=2.0, upper=10.0)

)

For all the calibration parameters, we use a helper function CalPrior which builds the prior given the named tuple of station priors and a JonesCache that specifies the segmentation scheme. For the gain products, we use the scancache, while for every other quantity, we use the trackcache.

julia

fwhmfac = 2.0*sqrt(2.0*log(2.0))
skyprior = (
    c = cprior,
    σ  = truncated(Normal(0.0, 1.0); lower=0.0),
    p  = cprior,
    p0 = Normal(-2.0, 2.0),
    pσ =  truncated(Normal(0.0, 1.0); lower=0.01),
    angparams = ImageSphericalUniform(nx, ny),
    )

skym = SkyModel(sky, skyprior, grid; metadata=skymeta)

 SkyModel
  with map: sky
   on grid: RectiGrid

Now we build the instrument model. Due to the complexity of VLBI the instrument model is critical to the success of imaging and getting reliable results. For this example we will use the standard instrument model used in polarized EHT analyses expressed in the RIME formalism. Our Jones decomposition will be given by GDR, where G are the complex gains, D are the d-terms, and R is what we call the ideal instrument response, which is how an ideal interferometer using the feed basis we observe relative to some reference basis.

Given the possible flexibility in different parameterizations of the individual Jones matrices each Jones matrix requires the user to specify a function that converts from parameters to specific parameterization f the jones matrices.

For the complex gain matrix, we used the JonesG jones matrix. The first argument is now a function that converts from the parameters to the complex gain matrix. In this case, we will use a amplitude and phase decomposition of the complex gain matrix. Note that since the gain matrix is a diagonal 2x2 matrix the function must return a 2-element tuple. The first element of the tuple is the gain for the first polarization feed (R) and the second is the gain for the second polarization feed (L).

julia

G = JonesG() do x
    gR = exp(x.lgR + 1im*x.gpR)
    gL = gR*exp(x.lgrat + 1im*x.gprat)
    return gR, gL
end

JonesG{Main.var"#1#2"}(Main.var"#1#2"())

Note that we are using the Julia do syntax here to define an anonymous function. This could've also been written as

julia

fgain(x) = (exp(x.lgR + 1im*x.gpR), exp(x.lgR + x.lgrat + 1im*(x.gpR + x.gprat)))
G = JonesG(fgain)

Similarly we provide a JonesD function for the leakage terms. Since we assume that we are in the small leakage limit, we will use the decomposition 1 d1 d2 1 Therefore, there are 2 free parameters for the JonesD our parameterization function must return a 2-element tuple. For d-terms we will use a re-im parameterization.

julia

D = JonesD() do x
    dR = complex(x.dRx, x.dRy)
    dL = complex(x.dLx, x.dLy)
    return dR, dL
end

JonesD{Main.var"#3#4"}(Main.var"#3#4"())

Finally we define our response Jones matrix. This matrix is a basis transform matrix plus the feed rotation angle for each station. These are typically set by the telescope so there are no free parameters, so no parameterization is necessary.

julia

R = JonesR(;add_fr=true)

JonesR{Nothing}(nothing, true)

Finally, we build our total Jones matrix by using the JonesSandwich function. The first argument is a function that specifies how to combine each Jones matrix. In this case, we are completely standard so we just need to multiply the different jones matrices. Note that if no function is provided, the default is to multiply the Jones matrices, so we could've removed the * argument in this case.

julia

J = JonesSandwich(splat(*), G, D, R)


intprior = (
    lgR  = ArrayPrior(IIDSitePrior(ScanSeg(), Normal(0.0, 0.1))),
    gpR  = ArrayPrior(IIDSitePrior(ScanSeg(), DiagonalVonMises(0.0, inv(π  ^2))); refant=SEFDReference(0.0), phase=false),
    lgrat= ArrayPrior(IIDSitePrior(ScanSeg(), Normal(0.0, 0.1)), phase=true),
    gprat= ArrayPrior(IIDSitePrior(ScanSeg(), Normal(0.0, 0.1)); refant = SingleReference(:AA, 0.0)),
    dRx  = ArrayPrior(IIDSitePrior(TrackSeg(), Normal(0.0, 0.2))),
    dRy  = ArrayPrior(IIDSitePrior(TrackSeg(), Normal(0.0, 0.2))),
    dLx  = ArrayPrior(IIDSitePrior(TrackSeg(), Normal(0.0, 0.2))),
    dLy  = ArrayPrior(IIDSitePrior(TrackSeg(), Normal(0.0, 0.2))),
)

intmodel = InstrumentModel(J, intprior)

InstrumentModel
  with Jones:  JonesSandwich
  with reference basis: CirBasis()

Putting it all together, we form our likelihood and posterior objects for optimization and sampling.

julia

post = VLBIPosterior(skym, intmodel, dvis)

VLBIPosterior
ObservedSkyModel
  with map: sky
   on grid: FourierDualDomainObservedInstrumentModel
  with Jones:  JonesSandwich
  with reference basis: CirBasis()Data Products: Comrade.EHTCoherencyDatum

Reconstructing the Image and Instrument Effects

To sample from this posterior, it is convenient to move from our constrained parameter space to an unconstrained one (i.e., the support of the transformed posterior is (-∞, ∞)). This transformation is done using the asflat function.

julia

tpost = asflat(post)

TransformedVLBIPosterior(
VLBIPosterior
ObservedSkyModel
  with map: sky
   on grid: FourierDualDomainObservedInstrumentModel
  with Jones:  JonesSandwich
  with reference basis: CirBasis()Data Products: Comrade.EHTCoherencyDatum
Transform: Params to ℝ^1329
)

We can also query the dimension of our posterior or the number of parameters we will sample.

Warning

This can often be different from what you would expect. This difference is especially true when using angular variables, where we often artificially increase the dimension of the parameter space to make sampling easier.

julia

ndim = dimension(tpost)

Now we optimize. Unlike other imaging examples, we move straight to gradient optimizers due to the higher dimension of the space. In addition the only AD package that can currently work with the polarized Comrade posterior is Zygote. Note that in the future we expect to shift entirely to Enzyme, and in fact large portions of Comrade's AD already uses Enzyme through custom rules.

julia

using Optimization
using OptimizationOptimisers
using Zygote
xopt, sol = comrade_opt(post, Optimisers.Adam(), Optimization.AutoZygote(); initial_params=prior_sample(rng, post), maxiters=20_000)

((sky = (c = (params = [0.5401395887749496 0.29767411969354224 … -0.44202524472491667 1.2848912786293005; 1.024304491101714 2.3047194901760846 … -0.6366268092217762 -0.039342471292512234; … ; 0.9589653802051481 -0.2256931168177961 … 0.2787567139000131 1.1719279875059738; -0.23142172521951512 -0.11168615273270131 … -0.38060192190595454 0.6800366096935025], hyperparams = 6.739499837007305), σ = 2.5498046691959524, p = (params = [0.23610273440709886 0.4015155853641987 … 0.14427969224493284 -0.21502959745303812; 0.5334285958819853 -0.13128015473102583 … 0.21499978340206696 0.03777380582529088; … ; 0.15297346552193392 0.24758928336322192 … 0.21789752580813157 -0.21024639209623136; 0.2567314170358642 0.34031551627289053 … 0.14774081848070256 0.03697616428128368], hyperparams = 3.872471458316115), p0 = -0.968934804845571, pσ = 3.906070802304991, angparams = ([0.7158720623283233 -0.3004319874607667 … -0.750142309602187 -0.23899802124777383; -0.10941528688639413 0.03260053024908511 … -0.3662554801167876 0.8699071045011207; … ; 0.022442224709763498 0.7074189827297981 … -0.48484645926342235 -0.3238320398356123; 0.6091464658884841 0.6173199971342858 … -0.22230357299809542 0.1926636236963871], [0.4294274316065278 0.48277054677041736 … 0.2639198295903918 0.9436505625257833; -0.9907946976570112 0.9160352329341656 … -0.8173659755896637 0.46779731275061803; … ; 0.9910949118082337 -0.3036581079416378 … -0.625725792372899 0.7444175097655561; -0.7728651178885186 -0.7326807218109247 … -0.36806929436897545 0.40048303335588464], [0.5505626861326638 -0.8226014953070325 … 0.6063273364228889 0.2289182423583175; -0.07971299825197736 0.39977075612246527 … 0.44471314938146206 0.15629236616603084; … ; -0.13125251364414275 0.6382398736798277 … -0.6110573980396854 -0.5839309729157673; -0.17782039434985025 0.2865222173318286 … -0.9028322745532555 0.8958203324876292])), instrument = (lgR = [-0.03721218215044625, -0.0372121821449676, -0.028634119762634868, -0.028634119555731474, -0.054077979316809396, -0.054077979271123996, -0.04852849845225605, -0.048528498412512534, 0.006759834803850662, 0.006759834632062441  …  -0.0266757897428605, -0.02667578974310832, -0.010497613603166964, -0.010497613603089387, -0.0056154166837877485, -0.005615416683592862, 0.005993835729025875, 0.005993835771549176, -0.010430116369704751, -0.010430116368565268], gpR = [0.0, 2.2075016234047182, 0.0, -0.7285975391903294, 0.0, -1.7176716200692645, 0.0, 1.165569135166229, 0.0, 0.9175746560422089  …  0.0, -0.7103533445897671, 0.0, 3.0785947217850245, 0.0, -1.7081935724313206, 0.0, -2.454591504866054, 0.0, 2.906697978389427], lgrat = [0.031948859775745096, 0.00441020154804069, 0.06184150205016922, -0.031516618178385004, 0.015621688581589951, 0.018541653136894785, 0.058965968830081254, -0.020515507004280682, 0.073543059844518, -0.03603939590325872  …  0.0605910559890326, -0.011320084412421745, 0.04077367954362536, 0.012046829360283997, 0.0037371628625889053, 0.042577431410459936, -0.010067517575235066, 0.05621159603298218, 0.009297607427958009, 0.03387922794214224], gprat = [-0.028217036984062922, -0.005968493996799712, -0.11597873107032731, -0.22984409551670193, 0.19339912150352545, -0.11331292590787252, 0.05088281652865524, 0.013684230518784332, -0.113685012122558, 0.02008458365210478  …  -0.1043715777671645, 0.077646715543949, 0.14472576463083622, 0.020410763557118038, -0.001802102467009324, -0.13551027603276766, -0.12458286221443525, -0.03225674749761245, 0.09974915918705138, 0.10359360176203264], dRx = [0.01018457389377081, -0.08144095033755594, 0.09181148173247386, -0.04296075815020972, 0.03090723285369758, -0.007080831343523114, 0.07969773378305804], dRy = [-0.020028066106080426, 0.07150478127400729, -0.10164187695454097, 0.04966839606257122, -0.020435890022575096, 0.01818816773573022, -0.07102867436934572], dLx = [0.029290407430036802, -0.05860839862244992, 0.0887893920968169, -0.05867764040763219, 0.009854691010116353, -0.02400130410261375, 0.05890125535125114], dLy = [-0.03922431765585126, 0.048769141436910114, -0.07852577106117238, 0.06877960455520886, -0.0002946001008739742, 0.03961433810604454, -0.05188049519249888])), retcode: Default
u: [0.5401395887749496, 1.024304491101714, 2.4779059202197664, 2.351063561251729, 2.285344394959323, 2.1710795640935685, -0.004556254346461523, 0.08868322890028957, 0.9589653802051481, -0.23142172521951512  …  0.009854691010116353, -0.02400130410261375, 0.05890125535125114, -0.03922431765585126, 0.048769141436910114, -0.07852577106117238, 0.06877960455520886, -0.0002946001008739742, 0.03961433810604454, -0.05188049519249888]
Final objective value:     1230.4087827022386
)

Warning

Fitting polarized images is generally much harder than Stokes I imaging. This difficulty means that optimization can take a long time, and starting from a good starting location is often required.

Now let's evaluate our fits by plotting the residuals

julia

using Plots
residual(post, xopt)

These look reasonable, although there may be some minor overfitting. Let's compare our results to the ground truth values we know in this example. First, we will load the polarized truth

julia

imgtrue = Comrade.load(joinpath(__DIR, "..", "..", "Data", "polarized_gaussian.fits"), IntensityMap{StokesParams})

╭─────────────────────────────────────────────────╮
│ 1024×1024 IntensityMap{StokesParams{Float64},2} │
├─────────────────────────────────────────────────┴────────────────────── dims ┐
  ↓ X Sampled{Float64} LinRange{Float64}(-4.843402302490755e-10, 4.843402302490755e-10, 1024) ForwardOrdered Regular Points,
  → Y Sampled{Float64} LinRange{Float64}(-4.843402302490772e-10, 4.843402302490772e-10, 1024) ForwardOrdered Regular Points
├──────────────────────────────────────────────────────────────────── metadata ┤
  ComradeBase.MinimalHeader{Float64}("SgrA", 266.4168370833333, -28.99218944444445, 51544.0, 2.3e11)
└──────────────────────────────────────────────────────────────────────────────┘
  ↓ →          …  4.8434e-10
 -4.8434e-10       [6.86545e-65, 4.08459e-78, 1.02982e-65, 1.02982e-66]
 -4.83393e-10      [8.99558e-65, 5.35191e-78, 1.34934e-65, 1.34934e-66]
 -4.82446e-10      [1.17804e-64, 7.00873e-78, 1.76706e-65, 1.76706e-66]
  ⋮            ⋱  
  4.82446e-10      [1.17804e-64, 7.00873e-78, 1.76706e-65, 1.76706e-66]
  4.83393e-10      [8.99558e-65, 5.35191e-78, 1.34934e-65, 1.34934e-66]
  4.8434e-10       [6.86545e-65, 4.08459e-78, 1.02982e-65, 1.02982e-66]

Select a reasonable zoom in of the image.

julia

imgtruesub = regrid(imgtrue, imagepixels(fovx, fovy, nx*4, ny*4))
img = intensitymap(Comrade.skymodel(post, xopt), axisdims(imgtruesub))

#Plotting the results gives
import CairoMakie as CM
using DisplayAs
fig = CM.Figure(;size=(450, 350));
axs = [CM.Axis(fig[1, i], xreversed=true, aspect=1) for i in 1:2]
polimage!(axs[1], imgtruesub,
                   nvec = 8,
                   length_norm=1/2, plot_total=true, pcolormap=:RdBu,
                   pcolorrange=(-0.25, 0.25),); axs[1].title="True"
polimage!(axs[2], img,
                   nvec = 8,
                   length_norm=1/2, plot_total=true, pcolormap=:RdBu,
                   pcolorrange=(-0.25, 0.25),);axs[2].title="Recon."
CM.Colorbar(fig[2,:], colormap=:RdBu, vertical=false, colorrange=(-0.25, 0.25), label="Signed Polarization Fraction sign(V)*|p|", flipaxis=false)
CM.colgap!(fig.layout, 3)
CM.rowgap!(fig.layout, 3)
CM.hidedecorations!.(fig.content[1:2])
fig |> DisplayAs.PNG |> DisplayAs.Text

Let's compare some image statics, like the total linear polarization fraction

julia

ftrue = flux(imgtruesub);
@info "Linear polarization true image: $(abs(linearpol(ftrue))/ftrue.I)"
frecon = flux(img);
@info "Linear polarization recon image: $(abs(linearpol(frecon))/frecon.I)"

[ Info: Linear polarization true image: 0.14999999999999963
[ Info: Linear polarization recon image: 0.15137370291731148

And the Circular polarization fraction

julia

@info "Circular polarization true image: $(ftrue.V/ftrue.I)"
@info "Circular polarization recon image: $(frecon.V/frecon.I)"

[ Info: Circular polarization true image: 0.014999999999999986
[ Info: Circular polarization recon image: 0.0350246630306509

Because we also fit the instrument model, we can inspect their parameters. To do this, Comrade provides a caltable function that converts the flattened gain parameters to a tabular format based on the time and its segmentation.

julia

dR = caltable(complex.(xopt.instrument.dRx, xopt.instrument.dRy))

───────────────────────────────────────────────────┬────────────────────────────
                                              time │                    AA     ⋯
───────────────────────────────────────────────────┼────────────────────────────
 IntegrationTime{Float64}(51544, 11.8333, 23.6667) │ 0.0101846-0.0200281im  -0 ⋯
───────────────────────────────────────────────────┴────────────────────────────
                                                               6 columns omitted

We can compare this to the ground truth d-terms

time	AA	AP	AZ	JC	LM	PV	SM
0.0	0.01-0.02im	-0.08+0.07im	0.09-0.10im	-0.04+0.05im	0.03-0.02im	-0.01+0.02im	0.08-0.07im

And same for the left-handed dterms

julia

dL = caltable(complex.(xopt.instrument.dLx, xopt.instrument.dLy))

───────────────────────────────────────────────────┬────────────────────────────
                                              time │                    AA     ⋯
───────────────────────────────────────────────────┼────────────────────────────
 IntegrationTime{Float64}(51544, 11.8333, 23.6667) │ 0.0292904-0.0392243im  -0 ⋯
───────────────────────────────────────────────────┴────────────────────────────
                                                               6 columns omitted

time	AA	AP	AZ	JC	LM	PV	SM
0.0	0.03-0.04im	-0.06+0.05im	0.09-0.08im	-0.06+0.07im	0.01-0.00im	-0.03+0.04im	0.06-0.05im

Looking at the gain phase ratio

julia

gphase_ratio = caltable(xopt.instrument.gprat)

───────────────────────────────────────────────────┬────────────────────────────
                                              time │      AA           AP      ⋯
───────────────────────────────────────────────────┼────────────────────────────
      IntegrationTime{Float64}(51544, 0.0, 0.0002) │ missing      missing    m ⋯
 IntegrationTime{Float64}(51544, 0.333333, 0.0002) │ missing      missing    m ⋯
 IntegrationTime{Float64}(51544, 0.666667, 0.0002) │ missing      missing    m ⋯
      IntegrationTime{Float64}(51544, 1.0, 0.0002) │ missing      missing    m ⋯
  IntegrationTime{Float64}(51544, 1.33333, 0.0002) │ missing      missing    m ⋯
  IntegrationTime{Float64}(51544, 1.66667, 0.0002) │ missing      missing    m ⋯
  IntegrationTime{Float64}(51544, 9.66667, 0.0002) │     0.0     0.217334    m ⋯
     IntegrationTime{Float64}(51544, 10.0, 0.0002) │     0.0   0.00102176    m ⋯
  IntegrationTime{Float64}(51544, 10.3333, 0.0002) │     0.0   -0.0776838    m ⋯
  IntegrationTime{Float64}(51544, 10.6667, 0.0002) │     0.0    0.0412084    m ⋯
     IntegrationTime{Float64}(51544, 11.0, 0.0002) │     0.0    -0.121324    m ⋯
  IntegrationTime{Float64}(51544, 11.3333, 0.0002) │     0.0    0.0600186    m ⋯
  IntegrationTime{Float64}(51544, 11.6667, 0.0002) │     0.0   -0.0913801    m ⋯
     IntegrationTime{Float64}(51544, 12.0, 0.0002) │     0.0    0.0359793    m ⋯
  IntegrationTime{Float64}(51544, 12.3333, 0.0002) │     0.0   -0.0163904    m ⋯
  IntegrationTime{Float64}(51544, 12.6667, 0.0002) │     0.0    0.0631079    m ⋯
                         ⋮                         │    ⋮          ⋮           ⋱
───────────────────────────────────────────────────┴────────────────────────────
                                                   5 columns and 33 rows omitted

we see that they are all very small. Which should be the case since this data doesn't have gain corruptions! Similarly our gain ratio amplitudes are also very close to unity as expected.

julia

gamp_ratio   = caltable(exp.(xopt.instrument.lgrat))

───────────────────────────────────────────────────┬────────────────────────────
                                              time │      AA       AP       AZ ⋯
───────────────────────────────────────────────────┼────────────────────────────
      IntegrationTime{Float64}(51544, 0.0, 0.0002) │ missing  missing  missing ⋯
 IntegrationTime{Float64}(51544, 0.333333, 0.0002) │ missing  missing  missing ⋯
 IntegrationTime{Float64}(51544, 0.666667, 0.0002) │ missing  missing  missing ⋯
      IntegrationTime{Float64}(51544, 1.0, 0.0002) │ missing  missing  missing ⋯
  IntegrationTime{Float64}(51544, 1.33333, 0.0002) │ missing  missing  missing ⋯
  IntegrationTime{Float64}(51544, 1.66667, 0.0002) │ missing  missing  missing ⋯
  IntegrationTime{Float64}(51544, 9.66667, 0.0002) │  1.0162  1.04847  missing ⋯
     IntegrationTime{Float64}(51544, 10.0, 0.0002) │ 1.01739  1.02809  missing ⋯
  IntegrationTime{Float64}(51544, 10.3333, 0.0002) │ 1.02107  1.02874  missing ⋯
  IntegrationTime{Float64}(51544, 10.6667, 0.0002) │ 1.00787   1.0351  missing ⋯
     IntegrationTime{Float64}(51544, 11.0, 0.0002) │ 1.02001  1.01093  missing ⋯
  IntegrationTime{Float64}(51544, 11.3333, 0.0002) │ 1.01468  1.04773  missing ⋯
  IntegrationTime{Float64}(51544, 11.6667, 0.0002) │ 1.01983  1.00455  missing ⋯
     IntegrationTime{Float64}(51544, 12.0, 0.0002) │  1.0194  1.02904  missing ⋯
  IntegrationTime{Float64}(51544, 12.3333, 0.0002) │ 1.02248  1.00859  missing ⋯
  IntegrationTime{Float64}(51544, 12.6667, 0.0002) │ 1.01814  1.00138  missing ⋯
                         ⋮                         │    ⋮        ⋮        ⋮    ⋱
───────────────────────────────────────────────────┴────────────────────────────
                                                   4 columns and 33 rows omitted

Plotting the gain phases, we see some offsets from zero. This is because the prior on the gain product phases is very broad, so we can't phase center the image. For realistic data this is always the case since the atmosphere effectively scrambles the phases.

julia

gphaseR = caltable(xopt.instrument.gpR)
p = Plots.plot(gphaseR, layout=(3,3), size=(650,500));
Plots.plot!(p, gphase_ratio, layout=(3,3), size=(650,500));
p |> DisplayAs.PNG |> DisplayAs.Text

Finally, the product gain amplitudes are all very close to unity as well, as expected since gain corruptions have not been added to the data.

julia

gampr = caltable(exp.(xopt.instrument.lgR))
p = Plots.plot(gampr, layout=(3,3), size=(650,500))
Plots.plot!(p, gamp_ratio, layout=(3,3), size=(650,500))
p |> DisplayAs.PNG |> DisplayAs.Text

This page was generated using Literate.jl.

Hamaker J.P, Bregman J.D., Sault R.J. (1996) [https://articles.adsabs.harvard.edu/pdf/1996A%26AS..117..137H] ↩︎
Pesce D. (2021) [https://ui.adsabs.harvard.edu/abs/2021AJ....161..178P/abstract] ↩︎

Polarized Image and Instrumental Modeling ​

Introduction to Polarized Imaging ​

Load the Data ​

Load the Data ​

Reconstructing the Image and Instrument Effects ​

Polarized Image and Instrumental Modeling

Introduction to Polarized Imaging

Load the Data

Load the Data

Reconstructing the Image and Instrument Effects