Stokes I Simultaneous Image and Instrument Modeling

In this tutorial, we will create a preliminary reconstruction of the 2017 M87 data on April 6 by simultaneously creating an image and model for the instrument. By instrument model, we mean something akin to self-calibration in traditional VLBI imaging terminology. However, unlike traditional self-cal, we will solve for the gains each time we update the image self-consistently. This allows us to model the correlations between gains and the image.

To get started we load Comrade.

julia

using Comrade


using Pyehtim
using LinearAlgebra

    CondaPkg Found dependencies: /home/runner/.julia/packages/DimensionalData/hv9KC/CondaPkg.toml
    CondaPkg Found dependencies: /home/runner/.julia/packages/CondaPkg/0UqYV/CondaPkg.toml
    CondaPkg Found dependencies: /home/runner/.julia/packages/PythonCall/wkBj7/CondaPkg.toml
    CondaPkg Found dependencies: /home/runner/.julia/packages/Pyehtim/bQtHC/CondaPkg.toml
    CondaPkg Resolving changes
             + ehtim (pip)
             + libstdcxx
             + libstdcxx-ng
             + numpy
             + numpy (pip)
             + openssl
             + pandas
             + python
             + setuptools (pip)
             + uv
             + xarray
    CondaPkg Initialising pixi
             │ /home/runner/.julia/artifacts/cefba4912c2b400756d043a2563ef77a0088866b/bin/pixi
             │ init
             │ --format pixi
             └ /home/runner/work/Comrade.jl/Comrade.jl/examples/intermediate/StokesIImaging/.CondaPkg
✔ Created /home/runner/work/Comrade.jl/Comrade.jl/examples/intermediate/StokesIImaging/.CondaPkg/pixi.toml
    CondaPkg Wrote /home/runner/work/Comrade.jl/Comrade.jl/examples/intermediate/StokesIImaging/.CondaPkg/pixi.toml
             │ [dependencies]
             │ openssl = ">=3, <3.6, >=3, <3.6"
             │ libstdcxx = ">=3.4,<15.0"
             │ uv = ">=0.4"
             │ libstdcxx-ng = ">=3.4,<15.0"
             │ pandas = "<2"
             │ xarray = "*"
             │ numpy = ">=1.24, <2.0"
             │ 
             │     [dependencies.python]
             │     channel = "conda-forge"
             │     build = "*cp*"
             │     version = ">=3.10,<3.14, >=3.6,<=3.12"
             │ 
             │ [project]
             │ name = ".CondaPkg"
             │ platforms = ["linux-64"]
             │ channels = ["conda-forge"]
             │ channel-priority = "strict"
             │ description = "automatically generated by CondaPkg.jl"
             │ 
             │ [pypi-dependencies]
             │ ehtim = ">=1.2.10, <2.0"
             │ numpy = ">=1.24, <2.0"
             └ setuptools = "*"
    CondaPkg Installing packages
             │ /home/runner/.julia/artifacts/cefba4912c2b400756d043a2563ef77a0088866b/bin/pixi
             │ install
             └ --manifest-path /home/runner/work/Comrade.jl/Comrade.jl/examples/intermediate/StokesIImaging/.CondaPkg/pixi.toml
✔ The default environment has been installed.
/home/runner/work/Comrade.jl/Comrade.jl/examples/intermediate/StokesIImaging/.CondaPkg/.pixi/envs/default/lib/python3.11/site-packages/ehtim/__init__.py:58: UserWarning: pkg_resources is deprecated as an API. See https://setuptools.pypa.io/en/latest/pkg_resources.html. The pkg_resources package is slated for removal as early as 2025-11-30. Refrain from using this package or pin to Setuptools<81.
  import pkg_resources

For reproducibility we use a stable random number genreator

julia

using StableRNGs
rng = StableRNG(42)

StableRNGs.LehmerRNG(state=0x00000000000000000000000000000055)

Load the Data

To download the data visit https://doi.org/10.25739/g85n-f134 First we will load our data:

julia

obs = ehtim.obsdata.load_uvfits(joinpath(__DIR, "..", "..", "Data", "SR1_M87_2017_096_lo_hops_netcal_StokesI.uvfits"))

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

Now we do some minor preprocessing:

Scan average the data since the data have been preprocessed so that the gain phases coherent.
Add 1% systematic noise to deal with calibration issues that cause 1% non-closing errors.

julia

obs = scan_average(obs).add_fractional_noise(0.02)

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

Now we extract our complex visibilities.

julia

dvis = extract_table(obs, Visibilities())

EHTObservationTable{Comrade.EHTVisibilityDatum{:I}}
  source:      M87
  mjd:         57849
  bandwidth:   1.856e9
  sites:       [:AA, :AP, :AZ, :JC, :LM, :PV, :SM]
  nsamples:    274

##Building the Model/Posterior

Now, we must build our intensity/visibility model. That is, the model that takes in a named tuple of parameters and perhaps some metadata required to construct the model. For our model, we will use a raster or ContinuousImage for our image model. The model is given below:

The model construction is very similar to Imaging a Black Hole using only Closure Quantities, except we include a large scale gaussian since we want to model the zero baselines. For more information about the image model please read the closure-only example.

julia

function sky(θ, metadata)
    (; fg, c, σimg) = θ
    (; ftot, mimg) = metadata
    # Apply the GMRF fluctuations to the image
    rast = apply_fluctuations(CenteredLR(), mimg, σimg .* c.params)
    pimg = parent(rast)
    @. pimg = (ftot * (1 - fg)) * pimg
    m = ContinuousImage(rast, BSplinePulse{3}())
    x0, y0 = centroid(m)
    # Add a large-scale gaussian to deal with the over-resolved mas flux
    g = modify(Gaussian(), Stretch(μas2rad(500.0), μas2rad(500.0)), Renormalize(ftot * fg))
    return shifted(m, -x0, -y0) + g
end

sky (generic function with 1 method)

Now, let's set up our image model. The EHT's nominal resolution is 20-25 μas. Additionally, the EHT is not very sensitive to a larger field of view. Typically 60-80 μas is enough to describe the compact flux of M87. Given this, we only need to use a small number of pixels to describe our image.

julia

npix = 48
fovx = μas2rad(200.0)
fovy = μas2rad(200.0)

9.69627362219072e-10

Now let's form our cache's. First, we have our usual image cache which is needed to numerically compute the visibilities.

julia

grid = imagepixels(fovx, fovy, npix, npix)

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

Now we need to specify our image prior. For this work we will use a Gaussian Markov Random field prior Since we are using a Gaussian Markov random field prior we need to first specify our mean image. This behaves somewhat similary to a entropy regularizer in that it will start with an initial guess for the image structure. For this tutorial we will use a a symmetric Gaussian with a FWHM of 50 μas

julia

using VLBIImagePriors
using Distributions
fwhmfac = 2 * sqrt(2 * log(2))
mpr = modify(Gaussian(), Stretch(μas2rad(50.0) ./ fwhmfac))
mimg = intensitymap(mpr, grid)

┌ 48×48 IntensityMap{Float64, 2} ┐
├────────────────────────────────┴─────────────────────────────────────── dims ┐
  ↓ X Sampled{Float64} LinRange{Float64}(-4.747133960864206e-10, 4.747133960864206e-10, 48) ForwardOrdered Regular Points,
  → Y Sampled{Float64} LinRange{Float64}(-4.747133960864206e-10, 4.747133960864206e-10, 48) ForwardOrdered Regular Points
└──────────────────────────────────────────────────────────────────────────────┘
  ↓ →          -4.74713e-10  -4.54513e-10  -4.34312e-10  -4.14112e-10  -3.93911e-10  -3.73711e-10  -3.5351e-10   -3.33309e-10  -3.13109e-10  -2.92908e-10  -2.72708e-10  -2.52507e-10  -2.32307e-10  -2.12106e-10  -1.91905e-10  -1.71705e-10  -1.51504e-10  -1.31304e-10  -1.11103e-10  -9.09026e-11  -7.0702e-11   -5.05014e-11  -3.03009e-11  -1.01003e-11  1.01003e-11  3.03009e-11  5.05014e-11  7.0702e-11   9.09026e-11  1.11103e-10  1.31304e-10  1.51504e-10  1.71705e-10  1.91905e-10  2.12106e-10  2.32307e-10  2.52507e-10  2.72708e-10  2.92908e-10  3.13109e-10  3.33309e-10  3.5351e-10   3.73711e-10  3.93911e-10  4.14112e-10  4.34312e-10  4.54513e-10  4.74713e-10
 -4.74713e-10   3.56128e-12   8.63485e-12   2.01456e-11   4.52254e-11   9.76922e-11   2.03055e-10   4.0611e-10    7.81538e-10   1.44721e-9    2.57864e-9    4.42105e-9    7.29349e-9    1.15777e-8    1.76842e-8    2.59911e-8    3.67569e-8    5.00184e-8    6.54933e-8    8.25164e-8    1.00037e-7    1.16696e-7    1.30987e-7    1.41473e-7    1.47028e-7   1.47028e-7   1.41473e-7   1.30987e-7   1.16696e-7   1.00037e-7   8.25164e-8   6.54933e-8   5.00184e-8   3.67569e-8   2.59911e-8   1.76842e-8   1.15777e-8   7.29349e-9   4.42105e-9   2.57864e-9   1.44721e-9   7.81538e-10  4.0611e-10   2.03055e-10  9.76922e-11  4.52254e-11  2.01456e-11  8.63485e-12  3.56128e-12
 -4.54513e-10   8.63485e-12   2.09365e-11   4.88461e-11   1.09656e-10   2.3687e-10    4.92338e-10   9.84676e-10   1.89496e-9    3.50899e-9    6.2523e-9     1.07195e-8    1.76842e-8    2.80719e-8    4.2878e-8     6.30193e-8    8.91227e-8    1.21277e-7    1.58799e-7    2.00074e-7    2.42555e-7    2.82947e-7    3.17597e-7    3.43024e-7    3.56491e-7   3.56491e-7   3.43024e-7   3.17597e-7   2.82947e-7   2.42555e-7   2.00074e-7   1.58799e-7   1.21277e-7   8.91227e-8   6.30193e-8   4.2878e-8    2.80719e-8   1.76842e-8   1.07195e-8   6.2523e-9    3.50899e-9   1.89496e-9   9.84676e-10  4.92338e-10  2.3687e-10   1.09656e-10  4.88461e-11  2.09365e-11  8.63485e-12
 -4.34312e-10   2.01456e-11   4.88461e-11   1.13961e-10   2.55833e-10   5.52631e-10   1.14865e-9    2.29731e-9    4.42105e-9    8.18667e-9    1.4587e-8     2.50092e-8    4.12582e-8    6.54933e-8    1.00037e-7    1.47028e-7    2.07928e-7    2.82947e-7    3.70486e-7    4.66784e-7    5.65894e-7    6.60132e-7    7.40973e-7    8.00295e-7    8.31714e-7   8.31714e-7   8.00295e-7   7.40973e-7   6.60132e-7   5.65894e-7   4.66784e-7   3.70486e-7   2.82947e-7   2.07928e-7   1.47028e-7   1.00037e-7   6.54933e-8   4.12582e-8   2.50092e-8   1.4587e-8    8.18667e-9   4.42105e-9   2.29731e-9   1.14865e-9   5.52631e-10  2.55833e-10  1.13961e-10  4.88461e-11  2.01456e-11
 -4.14112e-10   4.52254e-11   1.09656e-10   2.55833e-10   5.74327e-10   1.24061e-9    2.57864e-9    5.15728e-9    9.92491e-9    1.83784e-8    3.27467e-8    5.61438e-8    9.26216e-8    1.47028e-7    2.24575e-7    3.30066e-7    4.66784e-7    6.35194e-7    8.31714e-7    1.04789e-6    1.27039e-6    1.48195e-6    1.66343e-6    1.7966e-6     1.86713e-6   1.86713e-6   1.7966e-6    1.66343e-6   1.48195e-6   1.27039e-6   1.04789e-6   8.31714e-7   6.35194e-7   4.66784e-7   3.30066e-7   2.24575e-7   1.47028e-7   9.26216e-8   5.61438e-8   3.27467e-8   1.83784e-8   9.92491e-9   5.15728e-9   2.57864e-9   1.24061e-9   5.74327e-10  2.55833e-10  1.09656e-10  4.52254e-11
 -3.93911e-10   9.76922e-11   2.3687e-10    5.52631e-10   1.24061e-9    2.67987e-9    5.57017e-9    1.11403e-8    2.1439e-8     3.96996e-8    7.07367e-8    1.21277e-7    2.00074e-7    3.17597e-7    4.85109e-7    7.12982e-7    1.00831e-6    1.3721e-6     1.7966e-6     2.26358e-6    2.74419e-6    3.20118e-6    3.5932e-6     3.88087e-6    4.03323e-6   4.03323e-6   3.88087e-6   3.5932e-6    3.20118e-6   2.74419e-6   2.26358e-6   1.7966e-6    1.3721e-6    1.00831e-6   7.12982e-7   4.85109e-7   3.17597e-7   2.00074e-7   1.21277e-7   7.07367e-8   3.96996e-8   2.1439e-8    1.11403e-8   5.57017e-9   2.67987e-9   1.24061e-9   5.52631e-10  2.3687e-10   9.76922e-11
 -3.73711e-10   2.03055e-10   4.92338e-10   1.14865e-9    2.57864e-9    5.57017e-9    1.15777e-8    2.31554e-8    4.45613e-8    8.25164e-8    1.47028e-7    2.52077e-7    4.15857e-7    6.60132e-7    1.00831e-6    1.48195e-6    2.09579e-6    2.85193e-6    3.73427e-6    4.70488e-6    5.70385e-6    6.65371e-6    7.46854e-6    8.06647e-6    8.38315e-6   8.38315e-6   8.06647e-6   7.46854e-6   6.65371e-6   5.70385e-6   4.70488e-6   3.73427e-6   2.85193e-6   2.09579e-6   1.48195e-6   1.00831e-6   6.60132e-7   4.15857e-7   2.52077e-7   1.47028e-7   8.25164e-8   4.45613e-8   2.31554e-8   1.15777e-8   5.57017e-9   2.57864e-9   1.14865e-9   4.92338e-10  2.03055e-10
 -3.5351e-10    4.0611e-10    9.84676e-10   2.29731e-9    5.15728e-9    1.11403e-8    2.31554e-8    4.63108e-8    8.91227e-8    1.65033e-7    2.94055e-7    5.04154e-7    8.31714e-7    1.32026e-6    2.01662e-6    2.96389e-6    4.19157e-6    5.70385e-6    7.46854e-6    9.40977e-6    1.14077e-5    1.33074e-5    1.49371e-5    1.61329e-5    1.67663e-5   1.67663e-5   1.61329e-5   1.49371e-5   1.33074e-5   1.14077e-5   9.40977e-6   7.46854e-6   5.70385e-6   4.19157e-6   2.96389e-6   2.01662e-6   1.32026e-6   8.31714e-7   5.04154e-7   2.94055e-7   1.65033e-7   8.91227e-8   4.63108e-8   2.31554e-8   1.11403e-8   5.15728e-9   2.29731e-9   9.84676e-10  4.0611e-10
 -3.33309e-10   7.81538e-10   1.89496e-9    4.42105e-9    9.92491e-9    2.1439e-8     4.45613e-8    8.91227e-8    1.71512e-7    3.17597e-7    5.65894e-7    9.70218e-7    1.60059e-6    2.54078e-6    3.88087e-6    5.70385e-6    8.06647e-6    1.09768e-5    1.43728e-5    1.81086e-5    2.19535e-5    2.56094e-5    2.87456e-5    3.1047e-5     3.22659e-5   3.22659e-5   3.1047e-5    2.87456e-5   2.56094e-5   2.19535e-5   1.81086e-5   1.43728e-5   1.09768e-5   8.06647e-6   5.70385e-6   3.88087e-6   2.54078e-6   1.60059e-6   9.70218e-7   5.65894e-7   3.17597e-7   1.71512e-7   8.91227e-8   4.45613e-8   2.1439e-8    9.92491e-9   4.42105e-9   1.89496e-9   7.81538e-10
 -3.13109e-10   1.44721e-9    3.50899e-9    8.18667e-9    1.83784e-8    3.96996e-8    8.25164e-8    1.65033e-7    3.17597e-7    5.8811e-7     1.04789e-6    1.7966e-6     2.96389e-6    4.70488e-6    7.1864e-6     1.05621e-5    1.49371e-5    2.03262e-5    2.66148e-5    3.35326e-5    4.06524e-5    4.74222e-5    5.32297e-5    5.74912e-5    5.97483e-5   5.97483e-5   5.74912e-5   5.32297e-5   4.74222e-5   4.06524e-5   3.35326e-5   2.66148e-5   2.03262e-5   1.49371e-5   1.05621e-5   7.1864e-6    4.70488e-6   2.96389e-6   1.7966e-6    1.04789e-6   5.8811e-7    3.17597e-7   1.65033e-7   8.25164e-8   3.96996e-8   1.83784e-8   8.18667e-9   3.50899e-9   1.44721e-9
 -2.92908e-10   2.57864e-9    6.2523e-9     1.4587e-8     3.27467e-8    7.07367e-8    1.47028e-7    2.94055e-7    5.65894e-7    1.04789e-6    1.86713e-6    3.20118e-6    5.28105e-6    8.38315e-6    1.28047e-5    1.88195e-5    2.66148e-5    3.62172e-5    4.74222e-5    5.97483e-5    7.24344e-5    8.44968e-5    9.48445e-5    0.000102438   0.000106459  0.000106459  0.000102438  9.48445e-5   8.44968e-5   7.24344e-5   5.97483e-5   4.74222e-5   3.62172e-5   2.66148e-5   1.88195e-5   1.28047e-5   8.38315e-6   5.28105e-6   3.20118e-6   1.86713e-6   1.04789e-6   5.65894e-7   2.94055e-7   1.47028e-7   7.07367e-8   3.27467e-8   1.4587e-8    6.2523e-9    2.57864e-9
 -2.72708e-10   4.42105e-9    1.07195e-8    2.50092e-8    5.61438e-8    1.21277e-7    2.52077e-7    5.04154e-7    9.70218e-7    1.7966e-6     3.20118e-6    5.48838e-6    9.0543e-6     1.43728e-5    2.19535e-5    3.22659e-5    4.56308e-5    6.2094e-5     8.13049e-5    0.000102438   0.000124188   0.000144869   0.00016261    0.000175628   0.000182523  0.000182523  0.000175628  0.00016261   0.000144869  0.000124188  0.000102438  8.13049e-5   6.2094e-5    4.56308e-5   3.22659e-5   2.19535e-5   1.43728e-5   9.0543e-6    5.48838e-6   3.20118e-6   1.7966e-6    9.70218e-7   5.04154e-7   2.52077e-7   1.21277e-7   5.61438e-8   2.50092e-8   1.07195e-8   4.42105e-9
 -2.52507e-10   7.29349e-9    1.76842e-8    4.12582e-8    9.26216e-8    2.00074e-7    4.15857e-7    8.31714e-7    1.60059e-6    2.96389e-6    5.28105e-6    9.0543e-6     1.49371e-5    2.37111e-5    3.62172e-5    5.32297e-5    7.52781e-5    0.000102438   0.00013413    0.000168994   0.000204875   0.000238993   0.000268261   0.000289738   0.000301113  0.000301113  0.000289738  0.000268261  0.000238993  0.000204875  0.000168994  0.00013413   0.000102438  7.52781e-5   5.32297e-5   3.62172e-5   2.37111e-5   1.49371e-5   9.0543e-6    5.28105e-6   2.96389e-6   1.60059e-6   8.31714e-7   4.15857e-7   2.00074e-7   9.26216e-8   4.12582e-8   1.76842e-8   7.29349e-9
 -2.32307e-10   1.15777e-8    2.80719e-8    6.54933e-8    1.47028e-7    3.17597e-7    6.60132e-7    1.32026e-6    2.54078e-6    4.70488e-6    8.38315e-6    1.43728e-5    2.37111e-5    3.76391e-5    5.74912e-5    8.44968e-5    0.000119497   0.00016261    0.000212919   0.000268261   0.00032522    0.000379378   0.000425837   0.00045993    0.000477986  0.000477986  0.00045993   0.000425837  0.000379378  0.00032522   0.000268261  0.000212919  0.00016261   0.000119497  8.44968e-5   5.74912e-5   3.76391e-5   2.37111e-5   1.43728e-5   8.38315e-6   4.70488e-6   2.54078e-6   1.32026e-6   6.60132e-7   3.17597e-7   1.47028e-7   6.54933e-8   2.80719e-8   1.15777e-8
 -2.12106e-10   1.76842e-8    4.2878e-8     1.00037e-7    2.24575e-7    4.85109e-7    1.00831e-6    2.01662e-6    3.88087e-6    7.1864e-6     1.28047e-5    2.19535e-5    3.62172e-5    5.74912e-5    8.78141e-5    0.000129063   0.000182523   0.000248376   0.00032522    0.000409751   0.000496752   0.000579475   0.000650439   0.000702513   0.000730093  0.000730093  0.000702513  0.000650439  0.000579475  0.000496752  0.000409751  0.00032522   0.000248376  0.000182523  0.000129063  8.78141e-5   5.74912e-5   3.62172e-5   2.19535e-5   1.28047e-5   7.1864e-6    3.88087e-6   2.01662e-6   1.00831e-6   4.85109e-7   2.24575e-7   1.00037e-7   4.2878e-8    1.76842e-8
 -1.91905e-10   2.59911e-8    6.30193e-8    1.47028e-7    3.30066e-7    7.12982e-7    1.48195e-6    2.96389e-6    5.70385e-6    1.05621e-5    1.88195e-5    3.22659e-5    5.32297e-5    8.44968e-5    0.000129063   0.000189689   0.000268261   0.000365047   0.000477986   0.000602225   0.000730093   0.000851675   0.000955973   0.00103251    0.00107304   0.00107304   0.00103251   0.000955973  0.000851675  0.000730093  0.000602225  0.000477986  0.000365047  0.000268261  0.000189689  0.000129063  8.44968e-5   5.32297e-5   3.22659e-5   1.88195e-5   1.05621e-5   5.70385e-6   2.96389e-6   1.48195e-6   7.12982e-7   3.30066e-7   1.47028e-7   6.30193e-8   2.59911e-8
 -1.71705e-10   3.67569e-8    8.91227e-8    2.07928e-7    4.66784e-7    1.00831e-6    2.09579e-6    4.19157e-6    8.06647e-6    1.49371e-5    2.66148e-5    4.56308e-5    7.52781e-5    0.000119497   0.000182523   0.000268261   0.000379378   0.000516254   0.000675975   0.000851675   0.00103251    0.00120445    0.00135195    0.00146019    0.00151751   0.00151751   0.00146019   0.00135195   0.00120445   0.00103251   0.000851675  0.000675975  0.000516254  0.000379378  0.000268261  0.000182523  0.000119497  7.52781e-5   4.56308e-5   2.66148e-5   1.49371e-5   8.06647e-6   4.19157e-6   2.09579e-6   1.00831e-6   4.66784e-7   2.07928e-7   8.91227e-8   3.67569e-8
 -1.51504e-10   5.00184e-8    1.21277e-7    2.82947e-7    6.35194e-7    1.3721e-6     2.85193e-6    5.70385e-6    1.09768e-5    2.03262e-5    3.62172e-5    6.2094e-5     0.000102438   0.00016261    0.000248376   0.000365047   0.000516254   0.000702513   0.00091986    0.00115895    0.00140503    0.001639      0.00183972    0.00198701    0.00206502   0.00206502   0.00198701   0.00183972   0.001639     0.00140503   0.00115895   0.00091986   0.000702513  0.000516254  0.000365047  0.000248376  0.00016261   0.000102438  6.2094e-5    3.62172e-5   2.03262e-5   1.09768e-5   5.70385e-6   2.85193e-6   1.3721e-6    6.35194e-7   2.82947e-7   1.21277e-7   5.00184e-8
 -1.31304e-10   6.54933e-8    1.58799e-7    3.70486e-7    8.31714e-7    1.7966e-6     3.73427e-6    7.46854e-6    1.43728e-5    2.66148e-5    4.74222e-5    8.13049e-5    0.00013413    0.000212919   0.00032522    0.000477986   0.000675975   0.00091986    0.00120445    0.00151751    0.00183972    0.00214609    0.0024089     0.00260176    0.0027039    0.0027039    0.00260176   0.0024089    0.00214609   0.00183972   0.00151751   0.00120445   0.00091986   0.000675975  0.000477986  0.00032522   0.000212919  0.00013413   8.13049e-5   4.74222e-5   2.66148e-5   1.43728e-5   7.46854e-6   3.73427e-6   1.7966e-6    8.31714e-7   3.70486e-7   1.58799e-7   6.54933e-8
 -1.11103e-10   8.25164e-8    2.00074e-7    4.66784e-7    1.04789e-6    2.26358e-6    4.70488e-6    9.40977e-6    1.81086e-5    3.35326e-5    5.97483e-5    0.000102438   0.000168994   0.000268261   0.000409751   0.000602225   0.000851675   0.00115895    0.00151751    0.00191195    0.0023179     0.0027039     0.00303502    0.00327801    0.0034067    0.0034067    0.00327801   0.00303502   0.0027039    0.0023179    0.00191195   0.00151751   0.00115895   0.000851675  0.000602225  0.000409751  0.000268261  0.000168994  0.000102438  5.97483e-5   3.35326e-5   1.81086e-5   9.40977e-6   4.70488e-6   2.26358e-6   1.04789e-6   4.66784e-7   2.00074e-7   8.25164e-8
 -9.09026e-11   1.00037e-7    2.42555e-7    5.65894e-7    1.27039e-6    2.74419e-6    5.70385e-6    1.14077e-5    2.19535e-5    4.06524e-5    7.24344e-5    0.000124188   0.000204875   0.00032522    0.000496752   0.000730093   0.00103251    0.00140503    0.00183972    0.0023179     0.00281005    0.00327801    0.00367944    0.00397401    0.00413003   0.00413003   0.00397401   0.00367944   0.00327801   0.00281005   0.0023179    0.00183972   0.00140503   0.00103251   0.000730093  0.000496752  0.00032522   0.000204875  0.000124188  7.24344e-5   4.06524e-5   2.19535e-5   1.14077e-5   5.70385e-6   2.74419e-6   1.27039e-6   5.65894e-7   2.42555e-7   1.00037e-7
 -7.0702e-11    1.16696e-7    2.82947e-7    6.60132e-7    1.48195e-6    3.20118e-6    6.65371e-6    1.33074e-5    2.56094e-5    4.74222e-5    8.44968e-5    0.000144869   0.000238993   0.000379378   0.000579475   0.000851675   0.00120445    0.001639      0.00214609    0.0027039     0.00327801    0.00382389    0.00429217    0.0046358     0.0048178    0.0048178    0.0046358    0.00429217   0.00382389   0.00327801   0.0027039    0.00214609   0.001639     0.00120445   0.000851675  0.000579475  0.000379378  0.000238993  0.000144869  8.44968e-5   4.74222e-5   2.56094e-5   1.33074e-5   6.65371e-6   3.20118e-6   1.48195e-6   6.60132e-7   2.82947e-7   1.16696e-7
 -5.05014e-11   1.30987e-7    3.17597e-7    7.40973e-7    1.66343e-6    3.5932e-6     7.46854e-6    1.49371e-5    2.87456e-5    5.32297e-5    9.48445e-5    0.00016261    0.000268261   0.000425837   0.000650439   0.000955973   0.00135195    0.00183972    0.0024089     0.00303502    0.00367944    0.00429217    0.0048178     0.00520351    0.0054078    0.0054078    0.00520351   0.0048178    0.00429217   0.00367944   0.00303502   0.0024089    0.00183972   0.00135195   0.000955973  0.000650439  0.000425837  0.000268261  0.00016261   9.48445e-5   5.32297e-5   2.87456e-5   1.49371e-5   7.46854e-6   3.5932e-6    1.66343e-6   7.40973e-7   3.17597e-7   1.30987e-7
 -3.03009e-11   1.41473e-7    3.43024e-7    8.00295e-7    1.7966e-6     3.88087e-6    8.06647e-6    1.61329e-5    3.1047e-5     5.74912e-5    0.000102438   0.000175628   0.000289738   0.00045993    0.000702513   0.00103251    0.00146019    0.00198701    0.00260176    0.00327801    0.00397401    0.0046358     0.00520351    0.0056201     0.00584074   0.00584074   0.0056201    0.00520351   0.0046358    0.00397401   0.00327801   0.00260176   0.00198701   0.00146019   0.00103251   0.000702513  0.00045993   0.000289738  0.000175628  0.000102438  5.74912e-5   3.1047e-5    1.61329e-5   8.06647e-6   3.88087e-6   1.7966e-6    8.00295e-7   3.43024e-7   1.41473e-7
 -1.01003e-11   1.47028e-7    3.56491e-7    8.31714e-7    1.86713e-6    4.03323e-6    8.38315e-6    1.67663e-5    3.22659e-5    5.97483e-5    0.000106459   0.000182523   0.000301113   0.000477986   0.000730093   0.00107304    0.00151751    0.00206502    0.0027039     0.0034067     0.00413003    0.0048178     0.0054078     0.00584074    0.00607005   0.00607005   0.00584074   0.0054078    0.0048178    0.00413003   0.0034067    0.0027039    0.00206502   0.00151751   0.00107304   0.000730093  0.000477986  0.000301113  0.000182523  0.000106459  5.97483e-5   3.22659e-5   1.67663e-5   8.38315e-6   4.03323e-6   1.86713e-6   8.31714e-7   3.56491e-7   1.47028e-7
  1.01003e-11   1.47028e-7    3.56491e-7    8.31714e-7    1.86713e-6    4.03323e-6    8.38315e-6    1.67663e-5    3.22659e-5    5.97483e-5    0.000106459   0.000182523   0.000301113   0.000477986   0.000730093   0.00107304    0.00151751    0.00206502    0.0027039     0.0034067     0.00413003    0.0048178     0.0054078     0.00584074    0.00607005   0.00607005   0.00584074   0.0054078    0.0048178    0.00413003   0.0034067    0.0027039    0.00206502   0.00151751   0.00107304   0.000730093  0.000477986  0.000301113  0.000182523  0.000106459  5.97483e-5   3.22659e-5   1.67663e-5   8.38315e-6   4.03323e-6   1.86713e-6   8.31714e-7   3.56491e-7   1.47028e-7
  3.03009e-11   1.41473e-7    3.43024e-7    8.00295e-7    1.7966e-6     3.88087e-6    8.06647e-6    1.61329e-5    3.1047e-5     5.74912e-5    0.000102438   0.000175628   0.000289738   0.00045993    0.000702513   0.00103251    0.00146019    0.00198701    0.00260176    0.00327801    0.00397401    0.0046358     0.00520351    0.0056201     0.00584074   0.00584074   0.0056201    0.00520351   0.0046358    0.00397401   0.00327801   0.00260176   0.00198701   0.00146019   0.00103251   0.000702513  0.00045993   0.000289738  0.000175628  0.000102438  5.74912e-5   3.1047e-5    1.61329e-5   8.06647e-6   3.88087e-6   1.7966e-6    8.00295e-7   3.43024e-7   1.41473e-7
  5.05014e-11   1.30987e-7    3.17597e-7    7.40973e-7    1.66343e-6    3.5932e-6     7.46854e-6    1.49371e-5    2.87456e-5    5.32297e-5    9.48445e-5    0.00016261    0.000268261   0.000425837   0.000650439   0.000955973   0.00135195    0.00183972    0.0024089     0.00303502    0.00367944    0.00429217    0.0048178     0.00520351    0.0054078    0.0054078    0.00520351   0.0048178    0.00429217   0.00367944   0.00303502   0.0024089    0.00183972   0.00135195   0.000955973  0.000650439  0.000425837  0.000268261  0.00016261   9.48445e-5   5.32297e-5   2.87456e-5   1.49371e-5   7.46854e-6   3.5932e-6    1.66343e-6   7.40973e-7   3.17597e-7   1.30987e-7
  7.0702e-11    1.16696e-7    2.82947e-7    6.60132e-7    1.48195e-6    3.20118e-6    6.65371e-6    1.33074e-5    2.56094e-5    4.74222e-5    8.44968e-5    0.000144869   0.000238993   0.000379378   0.000579475   0.000851675   0.00120445    0.001639      0.00214609    0.0027039     0.00327801    0.00382389    0.00429217    0.0046358     0.0048178    0.0048178    0.0046358    0.00429217   0.00382389   0.00327801   0.0027039    0.00214609   0.001639     0.00120445   0.000851675  0.000579475  0.000379378  0.000238993  0.000144869  8.44968e-5   4.74222e-5   2.56094e-5   1.33074e-5   6.65371e-6   3.20118e-6   1.48195e-6   6.60132e-7   2.82947e-7   1.16696e-7
  9.09026e-11   1.00037e-7    2.42555e-7    5.65894e-7    1.27039e-6    2.74419e-6    5.70385e-6    1.14077e-5    2.19535e-5    4.06524e-5    7.24344e-5    0.000124188   0.000204875   0.00032522    0.000496752   0.000730093   0.00103251    0.00140503    0.00183972    0.0023179     0.00281005    0.00327801    0.00367944    0.00397401    0.00413003   0.00413003   0.00397401   0.00367944   0.00327801   0.00281005   0.0023179    0.00183972   0.00140503   0.00103251   0.000730093  0.000496752  0.00032522   0.000204875  0.000124188  7.24344e-5   4.06524e-5   2.19535e-5   1.14077e-5   5.70385e-6   2.74419e-6   1.27039e-6   5.65894e-7   2.42555e-7   1.00037e-7
  1.11103e-10   8.25164e-8    2.00074e-7    4.66784e-7    1.04789e-6    2.26358e-6    4.70488e-6    9.40977e-6    1.81086e-5    3.35326e-5    5.97483e-5    0.000102438   0.000168994   0.000268261   0.000409751   0.000602225   0.000851675   0.00115895    0.00151751    0.00191195    0.0023179     0.0027039     0.00303502    0.00327801    0.0034067    0.0034067    0.00327801   0.00303502   0.0027039    0.0023179    0.00191195   0.00151751   0.00115895   0.000851675  0.000602225  0.000409751  0.000268261  0.000168994  0.000102438  5.97483e-5   3.35326e-5   1.81086e-5   9.40977e-6   4.70488e-6   2.26358e-6   1.04789e-6   4.66784e-7   2.00074e-7   8.25164e-8
  1.31304e-10   6.54933e-8    1.58799e-7    3.70486e-7    8.31714e-7    1.7966e-6     3.73427e-6    7.46854e-6    1.43728e-5    2.66148e-5    4.74222e-5    8.13049e-5    0.00013413    0.000212919   0.00032522    0.000477986   0.000675975   0.00091986    0.00120445    0.00151751    0.00183972    0.00214609    0.0024089     0.00260176    0.0027039    0.0027039    0.00260176   0.0024089    0.00214609   0.00183972   0.00151751   0.00120445   0.00091986   0.000675975  0.000477986  0.00032522   0.000212919  0.00013413   8.13049e-5   4.74222e-5   2.66148e-5   1.43728e-5   7.46854e-6   3.73427e-6   1.7966e-6    8.31714e-7   3.70486e-7   1.58799e-7   6.54933e-8
  1.51504e-10   5.00184e-8    1.21277e-7    2.82947e-7    6.35194e-7    1.3721e-6     2.85193e-6    5.70385e-6    1.09768e-5    2.03262e-5    3.62172e-5    6.2094e-5     0.000102438   0.00016261    0.000248376   0.000365047   0.000516254   0.000702513   0.00091986    0.00115895    0.00140503    0.001639      0.00183972    0.00198701    0.00206502   0.00206502   0.00198701   0.00183972   0.001639     0.00140503   0.00115895   0.00091986   0.000702513  0.000516254  0.000365047  0.000248376  0.00016261   0.000102438  6.2094e-5    3.62172e-5   2.03262e-5   1.09768e-5   5.70385e-6   2.85193e-6   1.3721e-6    6.35194e-7   2.82947e-7   1.21277e-7   5.00184e-8
  1.71705e-10   3.67569e-8    8.91227e-8    2.07928e-7    4.66784e-7    1.00831e-6    2.09579e-6    4.19157e-6    8.06647e-6    1.49371e-5    2.66148e-5    4.56308e-5    7.52781e-5    0.000119497   0.000182523   0.000268261   0.000379378   0.000516254   0.000675975   0.000851675   0.00103251    0.00120445    0.00135195    0.00146019    0.00151751   0.00151751   0.00146019   0.00135195   0.00120445   0.00103251   0.000851675  0.000675975  0.000516254  0.000379378  0.000268261  0.000182523  0.000119497  7.52781e-5   4.56308e-5   2.66148e-5   1.49371e-5   8.06647e-6   4.19157e-6   2.09579e-6   1.00831e-6   4.66784e-7   2.07928e-7   8.91227e-8   3.67569e-8
  1.91905e-10   2.59911e-8    6.30193e-8    1.47028e-7    3.30066e-7    7.12982e-7    1.48195e-6    2.96389e-6    5.70385e-6    1.05621e-5    1.88195e-5    3.22659e-5    5.32297e-5    8.44968e-5    0.000129063   0.000189689   0.000268261   0.000365047   0.000477986   0.000602225   0.000730093   0.000851675   0.000955973   0.00103251    0.00107304   0.00107304   0.00103251   0.000955973  0.000851675  0.000730093  0.000602225  0.000477986  0.000365047  0.000268261  0.000189689  0.000129063  8.44968e-5   5.32297e-5   3.22659e-5   1.88195e-5   1.05621e-5   5.70385e-6   2.96389e-6   1.48195e-6   7.12982e-7   3.30066e-7   1.47028e-7   6.30193e-8   2.59911e-8
  2.12106e-10   1.76842e-8    4.2878e-8     1.00037e-7    2.24575e-7    4.85109e-7    1.00831e-6    2.01662e-6    3.88087e-6    7.1864e-6     1.28047e-5    2.19535e-5    3.62172e-5    5.74912e-5    8.78141e-5    0.000129063   0.000182523   0.000248376   0.00032522    0.000409751   0.000496752   0.000579475   0.000650439   0.000702513   0.000730093  0.000730093  0.000702513  0.000650439  0.000579475  0.000496752  0.000409751  0.00032522   0.000248376  0.000182523  0.000129063  8.78141e-5   5.74912e-5   3.62172e-5   2.19535e-5   1.28047e-5   7.1864e-6    3.88087e-6   2.01662e-6   1.00831e-6   4.85109e-7   2.24575e-7   1.00037e-7   4.2878e-8    1.76842e-8
  2.32307e-10   1.15777e-8    2.80719e-8    6.54933e-8    1.47028e-7    3.17597e-7    6.60132e-7    1.32026e-6    2.54078e-6    4.70488e-6    8.38315e-6    1.43728e-5    2.37111e-5    3.76391e-5    5.74912e-5    8.44968e-5    0.000119497   0.00016261    0.000212919   0.000268261   0.00032522    0.000379378   0.000425837   0.00045993    0.000477986  0.000477986  0.00045993   0.000425837  0.000379378  0.00032522   0.000268261  0.000212919  0.00016261   0.000119497  8.44968e-5   5.74912e-5   3.76391e-5   2.37111e-5   1.43728e-5   8.38315e-6   4.70488e-6   2.54078e-6   1.32026e-6   6.60132e-7   3.17597e-7   1.47028e-7   6.54933e-8   2.80719e-8   1.15777e-8
  2.52507e-10   7.29349e-9    1.76842e-8    4.12582e-8    9.26216e-8    2.00074e-7    4.15857e-7    8.31714e-7    1.60059e-6    2.96389e-6    5.28105e-6    9.0543e-6     1.49371e-5    2.37111e-5    3.62172e-5    5.32297e-5    7.52781e-5    0.000102438   0.00013413    0.000168994   0.000204875   0.000238993   0.000268261   0.000289738   0.000301113  0.000301113  0.000289738  0.000268261  0.000238993  0.000204875  0.000168994  0.00013413   0.000102438  7.52781e-5   5.32297e-5   3.62172e-5   2.37111e-5   1.49371e-5   9.0543e-6    5.28105e-6   2.96389e-6   1.60059e-6   8.31714e-7   4.15857e-7   2.00074e-7   9.26216e-8   4.12582e-8   1.76842e-8   7.29349e-9
  2.72708e-10   4.42105e-9    1.07195e-8    2.50092e-8    5.61438e-8    1.21277e-7    2.52077e-7    5.04154e-7    9.70218e-7    1.7966e-6     3.20118e-6    5.48838e-6    9.0543e-6     1.43728e-5    2.19535e-5    3.22659e-5    4.56308e-5    6.2094e-5     8.13049e-5    0.000102438   0.000124188   0.000144869   0.00016261    0.000175628   0.000182523  0.000182523  0.000175628  0.00016261   0.000144869  0.000124188  0.000102438  8.13049e-5   6.2094e-5    4.56308e-5   3.22659e-5   2.19535e-5   1.43728e-5   9.0543e-6    5.48838e-6   3.20118e-6   1.7966e-6    9.70218e-7   5.04154e-7   2.52077e-7   1.21277e-7   5.61438e-8   2.50092e-8   1.07195e-8   4.42105e-9
  2.92908e-10   2.57864e-9    6.2523e-9     1.4587e-8     3.27467e-8    7.07367e-8    1.47028e-7    2.94055e-7    5.65894e-7    1.04789e-6    1.86713e-6    3.20118e-6    5.28105e-6    8.38315e-6    1.28047e-5    1.88195e-5    2.66148e-5    3.62172e-5    4.74222e-5    5.97483e-5    7.24344e-5    8.44968e-5    9.48445e-5    0.000102438   0.000106459  0.000106459  0.000102438  9.48445e-5   8.44968e-5   7.24344e-5   5.97483e-5   4.74222e-5   3.62172e-5   2.66148e-5   1.88195e-5   1.28047e-5   8.38315e-6   5.28105e-6   3.20118e-6   1.86713e-6   1.04789e-6   5.65894e-7   2.94055e-7   1.47028e-7   7.07367e-8   3.27467e-8   1.4587e-8    6.2523e-9    2.57864e-9
  3.13109e-10   1.44721e-9    3.50899e-9    8.18667e-9    1.83784e-8    3.96996e-8    8.25164e-8    1.65033e-7    3.17597e-7    5.8811e-7     1.04789e-6    1.7966e-6     2.96389e-6    4.70488e-6    7.1864e-6     1.05621e-5    1.49371e-5    2.03262e-5    2.66148e-5    3.35326e-5    4.06524e-5    4.74222e-5    5.32297e-5    5.74912e-5    5.97483e-5   5.97483e-5   5.74912e-5   5.32297e-5   4.74222e-5   4.06524e-5   3.35326e-5   2.66148e-5   2.03262e-5   1.49371e-5   1.05621e-5   7.1864e-6    4.70488e-6   2.96389e-6   1.7966e-6    1.04789e-6   5.8811e-7    3.17597e-7   1.65033e-7   8.25164e-8   3.96996e-8   1.83784e-8   8.18667e-9   3.50899e-9   1.44721e-9
  3.33309e-10   7.81538e-10   1.89496e-9    4.42105e-9    9.92491e-9    2.1439e-8     4.45613e-8    8.91227e-8    1.71512e-7    3.17597e-7    5.65894e-7    9.70218e-7    1.60059e-6    2.54078e-6    3.88087e-6    5.70385e-6    8.06647e-6    1.09768e-5    1.43728e-5    1.81086e-5    2.19535e-5    2.56094e-5    2.87456e-5    3.1047e-5     3.22659e-5   3.22659e-5   3.1047e-5    2.87456e-5   2.56094e-5   2.19535e-5   1.81086e-5   1.43728e-5   1.09768e-5   8.06647e-6   5.70385e-6   3.88087e-6   2.54078e-6   1.60059e-6   9.70218e-7   5.65894e-7   3.17597e-7   1.71512e-7   8.91227e-8   4.45613e-8   2.1439e-8    9.92491e-9   4.42105e-9   1.89496e-9   7.81538e-10
  3.5351e-10    4.0611e-10    9.84676e-10   2.29731e-9    5.15728e-9    1.11403e-8    2.31554e-8    4.63108e-8    8.91227e-8    1.65033e-7    2.94055e-7    5.04154e-7    8.31714e-7    1.32026e-6    2.01662e-6    2.96389e-6    4.19157e-6    5.70385e-6    7.46854e-6    9.40977e-6    1.14077e-5    1.33074e-5    1.49371e-5    1.61329e-5    1.67663e-5   1.67663e-5   1.61329e-5   1.49371e-5   1.33074e-5   1.14077e-5   9.40977e-6   7.46854e-6   5.70385e-6   4.19157e-6   2.96389e-6   2.01662e-6   1.32026e-6   8.31714e-7   5.04154e-7   2.94055e-7   1.65033e-7   8.91227e-8   4.63108e-8   2.31554e-8   1.11403e-8   5.15728e-9   2.29731e-9   9.84676e-10  4.0611e-10
  3.73711e-10   2.03055e-10   4.92338e-10   1.14865e-9    2.57864e-9    5.57017e-9    1.15777e-8    2.31554e-8    4.45613e-8    8.25164e-8    1.47028e-7    2.52077e-7    4.15857e-7    6.60132e-7    1.00831e-6    1.48195e-6    2.09579e-6    2.85193e-6    3.73427e-6    4.70488e-6    5.70385e-6    6.65371e-6    7.46854e-6    8.06647e-6    8.38315e-6   8.38315e-6   8.06647e-6   7.46854e-6   6.65371e-6   5.70385e-6   4.70488e-6   3.73427e-6   2.85193e-6   2.09579e-6   1.48195e-6   1.00831e-6   6.60132e-7   4.15857e-7   2.52077e-7   1.47028e-7   8.25164e-8   4.45613e-8   2.31554e-8   1.15777e-8   5.57017e-9   2.57864e-9   1.14865e-9   4.92338e-10  2.03055e-10
  3.93911e-10   9.76922e-11   2.3687e-10    5.52631e-10   1.24061e-9    2.67987e-9    5.57017e-9    1.11403e-8    2.1439e-8     3.96996e-8    7.07367e-8    1.21277e-7    2.00074e-7    3.17597e-7    4.85109e-7    7.12982e-7    1.00831e-6    1.3721e-6     1.7966e-6     2.26358e-6    2.74419e-6    3.20118e-6    3.5932e-6     3.88087e-6    4.03323e-6   4.03323e-6   3.88087e-6   3.5932e-6    3.20118e-6   2.74419e-6   2.26358e-6   1.7966e-6    1.3721e-6    1.00831e-6   7.12982e-7   4.85109e-7   3.17597e-7   2.00074e-7   1.21277e-7   7.07367e-8   3.96996e-8   2.1439e-8    1.11403e-8   5.57017e-9   2.67987e-9   1.24061e-9   5.52631e-10  2.3687e-10   9.76922e-11
  4.14112e-10   4.52254e-11   1.09656e-10   2.55833e-10   5.74327e-10   1.24061e-9    2.57864e-9    5.15728e-9    9.92491e-9    1.83784e-8    3.27467e-8    5.61438e-8    9.26216e-8    1.47028e-7    2.24575e-7    3.30066e-7    4.66784e-7    6.35194e-7    8.31714e-7    1.04789e-6    1.27039e-6    1.48195e-6    1.66343e-6    1.7966e-6     1.86713e-6   1.86713e-6   1.7966e-6    1.66343e-6   1.48195e-6   1.27039e-6   1.04789e-6   8.31714e-7   6.35194e-7   4.66784e-7   3.30066e-7   2.24575e-7   1.47028e-7   9.26216e-8   5.61438e-8   3.27467e-8   1.83784e-8   9.92491e-9   5.15728e-9   2.57864e-9   1.24061e-9   5.74327e-10  2.55833e-10  1.09656e-10  4.52254e-11
  4.34312e-10   2.01456e-11   4.88461e-11   1.13961e-10   2.55833e-10   5.52631e-10   1.14865e-9    2.29731e-9    4.42105e-9    8.18667e-9    1.4587e-8     2.50092e-8    4.12582e-8    6.54933e-8    1.00037e-7    1.47028e-7    2.07928e-7    2.82947e-7    3.70486e-7    4.66784e-7    5.65894e-7    6.60132e-7    7.40973e-7    8.00295e-7    8.31714e-7   8.31714e-7   8.00295e-7   7.40973e-7   6.60132e-7   5.65894e-7   4.66784e-7   3.70486e-7   2.82947e-7   2.07928e-7   1.47028e-7   1.00037e-7   6.54933e-8   4.12582e-8   2.50092e-8   1.4587e-8    8.18667e-9   4.42105e-9   2.29731e-9   1.14865e-9   5.52631e-10  2.55833e-10  1.13961e-10  4.88461e-11  2.01456e-11
  4.54513e-10   8.63485e-12   2.09365e-11   4.88461e-11   1.09656e-10   2.3687e-10    4.92338e-10   9.84676e-10   1.89496e-9    3.50899e-9    6.2523e-9     1.07195e-8    1.76842e-8    2.80719e-8    4.2878e-8     6.30193e-8    8.91227e-8    1.21277e-7    1.58799e-7    2.00074e-7    2.42555e-7    2.82947e-7    3.17597e-7    3.43024e-7    3.56491e-7   3.56491e-7   3.43024e-7   3.17597e-7   2.82947e-7   2.42555e-7   2.00074e-7   1.58799e-7   1.21277e-7   8.91227e-8   6.30193e-8   4.2878e-8    2.80719e-8   1.76842e-8   1.07195e-8   6.2523e-9    3.50899e-9   1.89496e-9   9.84676e-10  4.92338e-10  2.3687e-10   1.09656e-10  4.88461e-11  2.09365e-11  8.63485e-12
  4.74713e-10   3.56128e-12   8.63485e-12   2.01456e-11   4.52254e-11   9.76922e-11   2.03055e-10   4.0611e-10    7.81538e-10   1.44721e-9    2.57864e-9    4.42105e-9    7.29349e-9    1.15777e-8    1.76842e-8    2.59911e-8    3.67569e-8    5.00184e-8    6.54933e-8    8.25164e-8    1.00037e-7    1.16696e-7    1.30987e-7    1.41473e-7    1.47028e-7   1.47028e-7   1.41473e-7   1.30987e-7   1.16696e-7   1.00037e-7   8.25164e-8   6.54933e-8   5.00184e-8   3.67569e-8   2.59911e-8   1.76842e-8   1.15777e-8   7.29349e-9   4.42105e-9   2.57864e-9   1.44721e-9   7.81538e-10  4.0611e-10   2.03055e-10  9.76922e-11  4.52254e-11  2.01456e-11  8.63485e-12  3.56128e-12

Now we can form our metadata we need to fully define our model. We will also fix the total flux to be the observed value 1.1. This is because total flux is degenerate with a global shift in the gain amplitudes making the problem degenerate. To fix this we use the observed total flux as our value.

julia

skymeta = (; ftot = 1.1, mimg = mimg ./ flux(mimg))

(ftot = 1.1, mimg = [3.5612926136987482e-12 8.634895734275082e-12 2.0145713255487537e-11 4.522559713081575e-11 9.769269325671449e-11 2.0305606556587994e-10 4.061121311317591e-10 7.815415460537175e-10 1.447219108186102e-9 2.5786512967024023e-9 4.4210666159488215e-9 7.293527272855018e-9 1.15777528654888e-8 1.7684266463795293e-8 2.599117639243255e-8 3.67570741562095e-8 5.001865894743796e-8 6.549366049671296e-8 8.25168414944769e-8 1.0003731789487592e-7 1.166964363656801e-7 1.3098732099342573e-7 1.4147413171036213e-7 1.4702829662483802e-7 1.4702829662483802e-7 1.4147413171036213e-7 1.3098732099342573e-7 1.166964363656801e-7 1.0003731789487592e-7 8.25168414944769e-8 6.549366049671285e-8 5.001865894743796e-8 3.67570741562095e-8 2.599117639243255e-8 1.768426646379523e-8 1.15777528654888e-8 7.293527272855018e-9 4.421066615948815e-9 2.5786512967024106e-9 1.447219108186102e-9 7.815415460537175e-10 4.061121311317591e-10 2.0305606556587994e-10 9.769269325671449e-11 4.522559713081575e-11 2.0145713255487537e-11 8.634895734275082e-12 3.5612926136987482e-12; 8.634895734275082e-12 2.0936618365757505e-11 4.884634662835752e-11 1.0965634057779147e-10 2.368707971446701e-10 4.923403226202322e-10 9.846806452404627e-10 1.8949663771573584e-9 3.509002898489335e-9 6.252332368429731e-9 1.0719548603267813e-8 1.7684266463795293e-8 2.8072023187914533e-8 4.287819441307118e-8 6.301956129538964e-8 8.912311827874256e-8 1.212778481380712e-7 1.587993436624031e-7 2.0007463578975186e-7 2.4255569627614246e-7 2.8294826342072426e-7 3.175986873248062e-7 3.430255553045695e-7 3.564924731149702e-7 3.564924731149702e-7 3.430255553045695e-7 3.175986873248062e-7 2.8294826342072426e-7 2.4255569627614246e-7 2.0007463578975186e-7 1.587993436624031e-7 1.212778481380712e-7 8.912311827874256e-8 6.301956129538964e-8 4.287819441307118e-8 2.807202318791463e-8 1.7684266463795293e-8 1.0719548603267791e-8 6.252332368429754e-9 3.509002898489335e-9 1.8949663771573584e-9 9.846806452404627e-10 4.923403226202322e-10 2.368707971446701e-10 1.0965634057779147e-10 4.884634662835752e-11 2.0936618365757505e-11 8.634895734275082e-12; 2.0145713255487537e-11 4.884634662835752e-11 1.1396136363835936e-10 2.5583461131528484e-10 5.526333269935996e-10 1.1486585673815487e-9 2.297317134763089e-9 4.421066615948807e-9 8.186707562089107e-9 1.4587054545709983e-8 2.5009329473718933e-8 4.125842074723837e-8 6.549366049671273e-8 1.0003731789487575e-7 1.4702829662483776e-7 2.079294111394597e-7 2.8294826342072373e-7 3.704880916956411e-7 4.667857454627196e-7 5.658965268414476e-7 6.60134731955813e-7 7.40976183391281e-7 8.002985431590045e-7 8.317176445578389e-7 8.317176445578389e-7 8.002985431590045e-7 7.40976183391281e-7 6.60134731955813e-7 5.658965268414476e-7 4.667857454627196e-7 3.704880916956404e-7 2.8294826342072373e-7 2.079294111394597e-7 1.4702829662483776e-7 1.0003731789487539e-7 6.549366049671273e-8 4.125842074723837e-8 2.5009329473718887e-8 1.4587054545710036e-8 8.186707562089107e-9 4.421066615948807e-9 2.297317134763089e-9 1.1486585673815487e-9 5.526333269935996e-10 2.5583461131528484e-10 1.1396136363835936e-10 4.884634662835752e-11 2.0145713255487537e-11; 4.522559713081575e-11 1.0965634057779147e-10 2.5583461131528484e-10 5.743292836907731e-10 1.2406198723625213e-9 2.5786512967024023e-9 5.1573025934047955e-9 9.924958978900158e-9 1.837853707810478e-8 3.274683024835631e-8 5.6144046375829066e-8 9.262202292391026e-8 1.4702829662483752e-7 2.2457618550331592e-7 3.3006736597790707e-7 4.667857454627196e-7 6.351973746496115e-7 8.317176445578389e-7 1.0478985679474042e-6 1.2703947492992233e-6 1.481952366782562e-6 1.6634352891156783e-6 1.7966094840265274e-6 1.8671429818508787e-6 1.8671429818508787e-6 1.7966094840265274e-6 1.6634352891156783e-6 1.481952366782562e-6 1.2703947492992233e-6 1.0478985679474042e-6 8.317176445578389e-7 6.351973746496115e-7 4.667857454627196e-7 3.3006736597790707e-7 2.2457618550331592e-7 1.4702829662483802e-7 9.262202292391026e-8 5.614404637582897e-8 3.274683024835643e-8 1.837853707810478e-8 9.924958978900158e-9 5.1573025934047955e-9 2.578651296702407e-9 1.2406198723625213e-9 5.743292836907731e-10 2.5583461131528484e-10 1.0965634057779147e-10 4.522559713081575e-11; 9.769269325671449e-11 2.368707971446701e-10 5.526333269935996e-10 1.2406198723625213e-9 2.6798871508169387e-9 5.5701948924213965e-9 1.1140389784842775e-8 2.143909720653547e-8 3.9699835915600776e-8 7.073706585518068e-8 1.2127784813807075e-7 2.0007463578975112e-7 3.1759868732480456e-7 4.851113925522823e-7 7.129849462299381e-7 1.0083129807262293e-6 1.3721022212182733e-6 1.7966094840265248e-6 2.263586107365787e-6 2.7442044424365453e-6 3.2011941726360163e-6 3.5932189680530496e-6 3.88089114041826e-6 4.033251922904918e-6 4.033251922904918e-6 3.88089114041826e-6 3.5932189680530496e-6 3.2011941726360163e-6 2.7442044424365453e-6 2.263586107365787e-6 1.7966094840265248e-6 1.3721022212182733e-6 1.0083129807262293e-6 7.129849462299381e-7 4.851113925522823e-7 3.1759868732480573e-7 2.0007463578975112e-7 1.2127784813807057e-7 7.073706585518093e-8 3.9699835915600776e-8 2.143909720653547e-8 1.1140389784842775e-8 5.570194892421407e-9 2.6798871508169387e-9 1.2406198723625213e-9 5.526333269935996e-10 2.368707971446701e-10 9.769269325671449e-11; 2.0305606556587994e-10 4.923403226202322e-10 1.1486585673815487e-9 2.5786512967024023e-9 5.5701948924213965e-9 1.157775286548882e-8 2.315550573097756e-8 4.456155913937119e-8 8.251684149447703e-8 1.4702829662483776e-7 2.5207824518155823e-7 4.158588222789203e-7 6.601347319558142e-7 1.0083129807262331e-6 1.4819523667825646e-6 2.095797135894812e-6 2.8519397849197576e-6 3.734285963701764e-6 4.704905491994815e-6 5.703879569839515e-6 6.653741156462725e-6 7.468571927403523e-6 8.066503845809857e-6 8.383188543579245e-6 8.383188543579245e-6 8.066503845809857e-6 7.468571927403523e-6 6.653741156462725e-6 5.703879569839515e-6 4.704905491994815e-6 3.734285963701764e-6 2.85193978491976e-6 2.095797135894812e-6 1.4819523667825646e-6 1.0083129807262312e-6 6.601347319558153e-7 4.158588222789203e-7 2.5207824518155776e-7 1.4702829662483828e-7 8.251684149447703e-8 4.456155913937119e-8 2.315550573097756e-8 1.1577752865488842e-8 5.5701948924213965e-9 2.5786512967024023e-9 1.1486585673815487e-9 4.923403226202322e-10 2.0305606556587994e-10; 4.061121311317591e-10 9.846806452404627e-10 2.297317134763089e-9 5.1573025934047955e-9 1.1140389784842775e-8 2.315550573097756e-8 4.631101146195504e-8 8.912311827874205e-8 1.650336829889538e-7 2.9405659324967445e-7 5.041564903631155e-7 8.317176445578389e-7 1.320269463911624e-6 2.016625961452459e-6 2.963904733565124e-6 4.191594271789614e-6 5.703879569839505e-6 7.46857192740351e-6 9.409810983989606e-6 1.1407759139679001e-5 1.330748231292542e-5 1.4937143854807011e-5 1.6133007691619673e-5 1.6766377087158443e-5 1.6766377087158443e-5 1.6133007691619673e-5 1.4937143854807011e-5 1.330748231292542e-5 1.1407759139679001e-5 9.409810983989606e-6 7.46857192740351e-6 5.703879569839505e-6 4.191594271789614e-6 2.963904733565124e-6 2.016625961452459e-6 1.3202694639116285e-6 8.317176445578389e-7 5.041564903631147e-7 2.9405659324967556e-7 1.650336829889538e-7 8.912311827874205e-8 4.631101146195504e-8 2.31555057309776e-8 1.1140389784842775e-8 5.1573025934047955e-9 2.297317134763089e-9 9.846806452404627e-10 4.061121311317591e-10; 7.815415460537175e-10 1.8949663771573584e-9 4.421066615948807e-9 9.924958978900158e-9 2.143909720653547e-8 4.456155913937119e-8 8.912311827874205e-8 1.7151277765228381e-7 3.1759868732480573e-7 5.658965268414456e-7 9.702227851045645e-7 1.6005970863180064e-6 2.5407894985984374e-6 3.88089114041826e-6 5.7038795698395e-6 8.066503845809836e-6 1.0976817769746173e-5 1.43728758722122e-5 1.810868885892628e-5 2.1953635539492345e-5 2.5609553381088117e-5 2.8745751744424377e-5 3.104712912334606e-5 3.2266015383239325e-5 3.2266015383239325e-5 3.104712912334606e-5 2.8745751744424377e-5 2.5609553381088117e-5 2.1953635539492345e-5 1.810868885892628e-5 1.43728758722122e-5 1.0976817769746181e-5 8.066503845809836e-6 5.7038795698395e-6 3.880891140418253e-6 2.5407894985984416e-6 1.6005970863180064e-6 9.702227851045628e-7 5.658965268414476e-7 3.1759868732480573e-7 1.7151277765228381e-7 8.912311827874205e-8 4.4561559139371265e-8 2.143909720653547e-8 9.924958978900158e-9 4.421066615948807e-9 1.8949663771573584e-9 7.815415460537175e-10; 1.447219108186102e-9 3.509002898489335e-9 8.186707562089107e-9 1.837853707810478e-8 3.9699835915600776e-8 8.251684149447703e-8 1.650336829889538e-7 3.1759868732480573e-7 5.881131864993532e-7 1.0478985679474042e-6 1.796609484026531e-6 2.963904733565132e-6 4.70490549199481e-6 7.186437936106125e-6 1.056215571129304e-5 1.493714385480706e-5 2.0326315988787573e-5 2.6614964625850903e-5 3.353275417431698e-5 4.065263197757515e-5 4.742247573704209e-5 5.322992925170177e-5 5.749150348884896e-5 5.97485754192282e-5 5.97485754192282e-5 5.749150348884896e-5 5.322992925170177e-5 4.742247573704209e-5 4.065263197757515e-5 3.353275417431698e-5 2.6614964625850903e-5 2.032631598878759e-5 1.493714385480706e-5 1.056215571129304e-5 7.186437936106112e-6 4.704905491994819e-6 2.963904733565132e-6 1.796609484026531e-6 1.047898567947408e-6 5.881131864993532e-7 3.1759868732480573e-7 1.650336829889538e-7 8.251684149447703e-8 3.9699835915600776e-8 1.837853707810478e-8 8.186707562089107e-9 3.509002898489335e-9 1.447219108186102e-9; 2.5786512967024023e-9 6.252332368429731e-9 1.4587054545709983e-8 3.274683024835631e-8 7.073706585518068e-8 1.4702829662483776e-7 2.9405659324967445e-7 5.658965268414456e-7 1.0478985679474042e-6 1.8671429818508722e-6 3.201194172636014e-6 5.2810778556465e-6 8.383188543579206e-6 1.2804776690544057e-5 1.8819621967979195e-5 2.661496462585082e-5 3.6217377717852536e-5 4.742247573704192e-5 5.9748575419228004e-5 7.243475543570508e-5 8.449724569034398e-5 9.484495147408377e-5 0.0001024382135243524 0.00010645985850340317 0.00010645985850340317 0.0001024382135243524 9.484495147408377e-5 8.449724569034398e-5 7.243475543570508e-5 5.9748575419228004e-5 4.742247573704192e-5 3.621737771785257e-5 2.661496462585082e-5 1.8819621967979195e-5 1.2804776690544035e-5 8.383188543579222e-6 5.2810778556465e-6 3.2011941726360103e-6 1.8671429818508787e-6 1.0478985679474042e-6 5.658965268414456e-7 2.9405659324967445e-7 1.4702829662483802e-7 7.073706585518068e-8 3.274683024835631e-8 1.4587054545709983e-8 6.252332368429731e-9 2.5786512967024023e-9; 4.4210666159488215e-9 1.0719548603267813e-8 2.5009329473718933e-8 5.6144046375829066e-8 1.2127784813807075e-7 2.5207824518155823e-7 5.041564903631155e-7 9.702227851045645e-7 1.796609484026531e-6 3.201194172636014e-6 5.488408884873099e-6 9.054344429463164e-6 1.43728758722122e-5 2.1953635539492382e-5 3.2266015383239386e-5 4.563103655871605e-5 6.209425824669222e-5 8.130526395515017e-5 0.00010243821352435257 0.0001241885164933844 0.0001448695108714104 0.0001626105279103003 0.00017562908431593906 0.00018252414623486416 0.00018252414623486416 0.00017562908431593906 0.0001626105279103003 0.0001448695108714104 0.00012418851649338445 0.00010243821352435257 8.130526395515017e-5 6.209425824669228e-5 4.563103655871605e-5 3.2266015383239386e-5 2.1953635539492345e-5 1.4372875872212226e-5 9.054344429463164e-6 5.488408884873096e-6 3.201194172636025e-6 1.796609484026531e-6 9.702227851045645e-7 5.041564903631155e-7 2.5207824518155866e-7 1.2127784813807075e-7 5.6144046375829066e-8 2.5009329473718933e-8 1.0719548603267813e-8 4.4210666159488215e-9; 7.293527272855018e-9 1.7684266463795293e-8 4.125842074723837e-8 9.262202292391026e-8 2.0007463578975112e-7 4.158588222789203e-7 8.317176445578389e-7 1.6005970863180064e-6 2.963904733565132e-6 5.2810778556465e-6 9.054344429463164e-6 1.4937143854807035e-5 2.371123786852098e-5 3.621737771785263e-5 5.322992925170173e-5 7.527848787191693e-5 0.00010243821352435266 0.00013413101669726773 0.00016899449138068818 0.00020487642704870515 0.00023899430167691242 0.0002682620333945352 0.0002897390217428208 0.0003011139514876675 0.0003011139514876675 0.0002897390217428208 0.0002682620333945352 0.00023899430167691231 0.00020487642704870529 0.00016899449138068818 0.00013413101669726773 0.00010243821352435266 7.527848787191693e-5 5.322992925170173e-5 3.621737771785257e-5 2.371123786852102e-5 1.4937143854807035e-5 9.054344429463147e-6 5.28107785564652e-6 2.963904733565132e-6 1.6005970863180064e-6 8.317176445578389e-7 4.1585882227892094e-7 2.0007463578975112e-7 9.262202292391026e-8 4.125842074723837e-8 1.7684266463795293e-8 7.293527272855018e-9; 1.15777528654888e-8 2.8072023187914533e-8 6.549366049671273e-8 1.4702829662483752e-7 3.1759868732480456e-7 6.601347319558142e-7 1.320269463911624e-6 2.5407894985984374e-6 4.70490549199481e-6 8.383188543579206e-6 1.43728758722122e-5 2.371123786852098e-5 3.763924393595837e-5 5.749150348884881e-5 8.449724569034405e-5 0.00011949715083845609 0.00016261052791030013 0.00021291971700680657 0.0002682620333945348 0.0003252210558206002 0.0003793798058963355 0.0004258394340136128 0.00045993202791079024 0.00047798860335382403 0.00047798860335382403 0.00045993202791079024 0.0004258394340136128 0.0003793798058963353 0.00032522105582060036 0.0002682620333945348 0.00021291971700680657 0.00016261052791030018 0.00011949715083845609 8.449724569034405e-5 5.749150348884871e-5 3.7639243935958425e-5 2.371123786852098e-5 1.4372875872212187e-5 8.383188543579237e-6 4.70490549199481e-6 2.5407894985984374e-6 1.320269463911624e-6 6.601347319558153e-7 3.1759868732480456e-7 1.4702829662483752e-7 6.549366049671273e-8 2.8072023187914533e-8 1.15777528654888e-8; 1.7684266463795293e-8 4.287819441307118e-8 1.0003731789487575e-7 2.2457618550331592e-7 4.851113925522823e-7 1.0083129807262331e-6 2.016625961452459e-6 3.88089114041826e-6 7.186437936106125e-6 1.2804776690544057e-5 2.1953635539492382e-5 3.621737771785263e-5 5.749150348884881e-5 8.781454215796957e-5 0.00012906406153295754 0.00018252414623486427 0.0002483770329867688 0.00032522105582060074 0.00040975285409741024 0.0004967540659735375 0.0005794780434856418 0.000650442111641201 0.0007025163372637559 0.0007300965849394566 0.0007300965849394566 0.0007025163372637559 0.000650442111641201 0.0005794780434856414 0.0004967540659735377 0.00040975285409741024 0.00032522105582060074 0.00024837703298676895 0.00018252414623486427 0.00012906406153295754 8.78145421579694e-5 5.7491503488848916e-5 3.621737771785263e-5 2.1953635539492362e-5 1.2804776690544102e-5 7.186437936106125e-6 3.88089114041826e-6 2.016625961452459e-6 1.0083129807262346e-6 4.851113925522823e-7 2.2457618550331592e-7 1.0003731789487575e-7 4.287819441307118e-8 1.7684266463795293e-8; 2.599117639243255e-8 6.301956129538964e-8 1.4702829662483776e-7 3.3006736597790707e-7 7.129849462299381e-7 1.4819523667825646e-6 2.963904733565124e-6 5.7038795698395e-6 1.056215571129304e-5 1.8819621967979195e-5 3.2266015383239386e-5 5.322992925170173e-5 8.449724569034405e-5 0.00012906406153295754 0.00018968990294816794 0.0002682620333945353 0.0003650482924697282 0.00047798860335382485 0.0006022279029753348 0.0007300965849394566 0.0008516788680272268 0.0009559772067076492 0.0010325124922636595 0.0010730481335781406 0.0010730481335781406 0.0010325124922636595 0.0009559772067076492 0.0008516788680272265 0.0007300965849394566 0.0006022279029753348 0.00047798860335382463 0.00036504829246972854 0.0002682620333945353 0.00018968990294816794 0.0001290640615329573 8.449724569034421e-5 5.322992925170173e-5 3.226601538323935e-5 1.8819621967979263e-5 1.056215571129304e-5 5.7038795698395e-6 2.963904733565124e-6 1.4819523667825673e-6 7.129849462299381e-7 3.3006736597790707e-7 1.4702829662483776e-7 6.301956129538964e-8 2.599117639243255e-8; 3.67570741562095e-8 8.912311827874256e-8 2.079294111394597e-7 4.667857454627196e-7 1.0083129807262293e-6 2.095797135894812e-6 4.191594271789614e-6 8.066503845809836e-6 1.493714385480706e-5 2.661496462585082e-5 4.563103655871605e-5 7.527848787191693e-5 0.00011949715083845609 0.00018252414623486427 0.0002682620333945353 0.00037937980589633605 0.0005162562461318299 0.000675977965522753 0.0008516788680272268 0.0010325124922636595 0.0012044558059506696 0.0013519559310455054 0.0014601931698789128 0.0015175192235853429 0.0015175192235853429 0.0014601931698789128 0.0013519559310455054 0.0012044558059506696 0.0010325124922636597 0.0008516788680272268 0.0006759779655227528 0.0005162562461318303 0.00037937980589633605 0.0002682620333945353 0.00018252414623486392 0.0001194971508384563 7.527848787191693e-5 4.563103655871602e-5 2.6614964625850903e-5 1.493714385480706e-5 8.066503845809836e-6 4.191594271789614e-6 2.0957971358948143e-6 1.0083129807262293e-6 4.667857454627196e-7 2.079294111394597e-7 8.912311827874256e-8 3.67570741562095e-8; 5.001865894743796e-8 1.212778481380712e-7 2.8294826342072373e-7 6.351973746496115e-7 1.3721022212182733e-6 2.8519397849197576e-6 5.703879569839505e-6 1.0976817769746173e-5 2.0326315988787573e-5 3.6217377717852536e-5 6.209425824669222e-5 0.00010243821352435266 0.00016261052791030013 0.0002483770329867688 0.0003650482924697282 0.0005162562461318299 0.0007025163372637559 0.0009198640558215816 0.0011589560869712824 0.0014050326745275112 0.0016390114163896403 0.0018397281116431623 0.0019870162638941495 0.0020650249845273186 0.0020650249845273186 0.0019870162638941495 0.0018397281116431623 0.0016390114163896401 0.0014050326745275118 0.0011589560869712824 0.0009198640558215815 0.0007025163372637562 0.0005162562461318299 0.0003650482924697282 0.00024837703298676836 0.00016261052791030042 0.00010243821352435266 6.209425824669217e-5 3.621737771785266e-5 2.0326315988787573e-5 1.0976817769746173e-5 5.703879569839505e-6 2.851939784919763e-6 1.3721022212182733e-6 6.351973746496115e-7 2.8294826342072373e-7 1.212778481380712e-7 5.001865894743796e-8; 6.549366049671296e-8 1.587993436624031e-7 3.704880916956411e-7 8.317176445578389e-7 1.7966094840265248e-6 3.734285963701764e-6 7.46857192740351e-6 1.43728758722122e-5 2.6614964625850903e-5 4.742247573704192e-5 8.130526395515017e-5 0.00013413101669726773 0.00021291971700680657 0.00032522105582060074 0.00047798860335382485 0.000675977965522753 0.0009198640558215816 0.00120445580595067 0.0015175192235853427 0.001839728111643163 0.0021460962671562816 0.002408911611901339 0.0026017684465648046 0.002703911862091011 0.002703911862091011 0.0026017684465648046 0.002408911611901339 0.002146096267156281 0.0018397281116431638 0.0015175192235853427 0.00120445580595067 0.000919864055821582 0.000675977965522753 0.00047798860335382485 0.0003252210558206002 0.00021291971700680693 0.00013413101669726773 8.130526395515009e-5 4.742247573704209e-5 2.6614964625850903e-5 1.43728758722122e-5 7.46857192740351e-6 3.734285963701767e-6 1.7966094840265248e-6 8.317176445578389e-7 3.704880916956411e-7 1.587993436624031e-7 6.549366049671296e-8; 8.25168414944769e-8 2.0007463578975186e-7 4.667857454627196e-7 1.0478985679474042e-6 2.263586107365787e-6 4.704905491994815e-6 9.409810983989606e-6 1.810868885892628e-5 3.353275417431698e-5 5.9748575419228004e-5 0.00010243821352435257 0.00016899449138068818 0.0002682620333945348 0.00040975285409741024 0.0006022279029753348 0.0008516788680272268 0.0011589560869712824 0.0015175192235853427 0.001911954413415297 0.002317912173942565 0.002703911862091009 0.0030350384471706845 0.0032780228327792802 0.003406715472108905 0.003406715472108905 0.0032780228327792802 0.0030350384471706845 0.002703911862091008 0.0023179121739425657 0.001911954413415297 0.0015175192235853427 0.0011589560869712833 0.0008516788680272268 0.0006022279029753348 0.00040975285409740943 0.0002682620333945353 0.00016899449138068818 0.00010243821352435248 5.97485754192282e-5 3.353275417431698e-5 1.810868885892628e-5 9.409810983989606e-6 4.704905491994819e-6 2.263586107365787e-6 1.0478985679474042e-6 4.667857454627196e-7 2.0007463578975186e-7 8.25168414944769e-8; 1.0003731789487592e-7 2.4255569627614246e-7 5.658965268414476e-7 1.2703947492992233e-6 2.7442044424365453e-6 5.703879569839515e-6 1.1407759139679001e-5 2.1953635539492345e-5 4.065263197757515e-5 7.243475543570508e-5 0.0001241885164933844 0.00020487642704870515 0.0003252210558206002 0.0004967540659735375 0.0007300965849394566 0.0010325124922636595 0.0014050326745275112 0.001839728111643163 0.002317912173942565 0.0028100653490550223 0.0032780228327792802 0.0036794562232863232 0.003974032527788298 0.004130049969054636 0.004130049969054636 0.003974032527788298 0.0036794562232863232 0.0032780228327792794 0.0028100653490550223 0.002317912173942565 0.0018397281116431623 0.0014050326745275118 0.0010325124922636595 0.0007300965849394566 0.0004967540659735367 0.00032522105582060074 0.00020487642704870515 0.00012418851649338426 7.243475543570533e-5 4.065263197757515e-5 2.1953635539492345e-5 1.1407759139679001e-5 5.703879569839521e-6 2.7442044424365453e-6 1.2703947492992233e-6 5.658965268414476e-7 2.4255569627614246e-7 1.0003731789487592e-7; 1.166964363656801e-7 2.8294826342072426e-7 6.60134731955813e-7 1.481952366782562e-6 3.2011941726360163e-6 6.653741156462725e-6 1.330748231292542e-5 2.5609553381088117e-5 4.742247573704209e-5 8.449724569034398e-5 0.0001448695108714104 0.00023899430167691242 0.0003793798058963355 0.0005794780434856418 0.0008516788680272268 0.0012044558059506696 0.0016390114163896403 0.0021460962671562816 0.002703911862091009 0.0032780228327792802 0.003823908826830595 0.004292192534312561 0.0046358243478851315 0.004817823223802676 0.004817823223802676 0.0046358243478851315 0.004292192534312561 0.003823908826830595 0.0032780228327792807 0.002703911862091009 0.002146096267156281 0.0016390114163896412 0.0012044558059506696 0.0008516788680272268 0.0005794780434856407 0.0003793798058963361 0.00023899430167691242 0.0001448695108714103 8.449724569034428e-5 4.742247573704209e-5 2.5609553381088117e-5 1.330748231292542e-5 6.653741156462732e-6 3.2011941726360163e-6 1.481952366782562e-6 6.60134731955813e-7 2.8294826342072426e-7 1.166964363656801e-7; 1.3098732099342573e-7 3.175986873248062e-7 7.40976183391281e-7 1.6634352891156783e-6 3.5932189680530496e-6 7.468571927403523e-6 1.4937143854807011e-5 2.8745751744424377e-5 5.322992925170177e-5 9.484495147408377e-5 0.0001626105279103003 0.0002682620333945352 0.0004258394340136128 0.000650442111641201 0.0009559772067076492 0.0013519559310455054 0.0018397281116431623 0.002408911611901339 0.0030350384471706845 0.0036794562232863232 0.004292192534312561 0.004817823223802676 0.005203536893129606 0.005407823724182016 0.005407823724182018 0.005203536893129606 0.004817823223802676 0.00429219253431256 0.0036794562232863245 0.0030350384471706845 0.0024089116119013388 0.0018397281116431631 0.0013519559310455054 0.0009559772067076492 0.0006504421116412 0.00042583943401361353 0.0002682620333945352 0.00016261052791030005 9.484495147408411e-5 5.322992925170177e-5 2.8745751744424377e-5 1.4937143854807011e-5 7.46857192740353e-6 3.5932189680530496e-6 1.6634352891156783e-6 7.40976183391281e-7 3.175986873248062e-7 1.3098732099342573e-7; 1.4147413171036213e-7 3.430255553045695e-7 8.002985431590045e-7 1.7966094840265274e-6 3.88089114041826e-6 8.066503845809857e-6 1.6133007691619673e-5 3.104712912334606e-5 5.749150348884896e-5 0.0001024382135243524 0.00017562908431593906 0.0002897390217428208 0.00045993202791079024 0.0007025163372637559 0.0010325124922636595 0.0014601931698789128 0.0019870162638941495 0.0026017684465648046 0.0032780228327792802 0.003974032527788298 0.0046358243478851315 0.005203536893129606 0.005620130698110045 0.005840772679515649 0.005840772679515649 0.005620130698110045 0.005203536893129606 0.00463582434788513 0.003974032527788299 0.0032780228327792802 0.002601768446564804 0.0019870162638941503 0.0014601931698789128 0.0010325124922636595 0.0007025163372637547 0.000459932027910791 0.0002897390217428208 0.0001756290843159388 0.00010243821352435276 5.749150348884896e-5 3.104712912334606e-5 1.6133007691619673e-5 8.066503845809865e-6 3.88089114041826e-6 1.7966094840265274e-6 8.002985431590045e-7 3.430255553045695e-7 1.4147413171036213e-7; 1.4702829662483802e-7 3.564924731149702e-7 8.317176445578389e-7 1.8671429818508787e-6 4.033251922904918e-6 8.383188543579245e-6 1.6766377087158443e-5 3.2266015383239325e-5 5.97485754192282e-5 0.00010645985850340317 0.00018252414623486416 0.0003011139514876675 0.00047798860335382403 0.0007300965849394566 0.0010730481335781406 0.0015175192235853429 0.0020650249845273186 0.002703911862091011 0.003406715472108905 0.004130049969054636 0.004817823223802676 0.005407823724182016 0.005840772679515649 0.006070076894341369 0.006070076894341369 0.005840772679515649 0.005407823724182016 0.004817823223802676 0.004130049969054637 0.003406715472108905 0.00270391186209101 0.00206502498452732 0.0015175192235853429 0.0010730481335781406 0.0007300965849394552 0.00047798860335382485 0.0003011139514876675 0.00018252414623486392 0.00010645985850340356 5.97485754192282e-5 3.2266015383239325e-5 1.6766377087158443e-5 8.38318854357925e-6 4.033251922904918e-6 1.8671429818508787e-6 8.317176445578389e-7 3.564924731149702e-7 1.4702829662483802e-7; 1.4702829662483802e-7 3.564924731149702e-7 8.317176445578389e-7 1.8671429818508787e-6 4.033251922904918e-6 8.383188543579245e-6 1.6766377087158443e-5 3.2266015383239325e-5 5.97485754192282e-5 0.00010645985850340317 0.00018252414623486416 0.0003011139514876675 0.00047798860335382403 0.0007300965849394566 0.0010730481335781406 0.0015175192235853429 0.0020650249845273186 0.002703911862091011 0.003406715472108905 0.004130049969054636 0.004817823223802676 0.005407823724182018 0.005840772679515649 0.006070076894341369 0.006070076894341371 0.005840772679515649 0.005407823724182018 0.004817823223802676 0.004130049969054637 0.003406715472108905 0.00270391186209101 0.00206502498452732 0.0015175192235853429 0.0010730481335781406 0.0007300965849394552 0.00047798860335382485 0.0003011139514876675 0.00018252414623486392 0.00010645985850340356 5.97485754192282e-5 3.2266015383239325e-5 1.6766377087158443e-5 8.38318854357925e-6 4.033251922904918e-6 1.8671429818508787e-6 8.317176445578389e-7 3.564924731149702e-7 1.4702829662483802e-7; 1.4147413171036213e-7 3.430255553045695e-7 8.002985431590045e-7 1.7966094840265274e-6 3.88089114041826e-6 8.066503845809857e-6 1.6133007691619673e-5 3.104712912334606e-5 5.749150348884896e-5 0.0001024382135243524 0.00017562908431593906 0.0002897390217428208 0.00045993202791079024 0.0007025163372637559 0.0010325124922636595 0.0014601931698789128 0.0019870162638941495 0.0026017684465648046 0.0032780228327792802 0.003974032527788298 0.0046358243478851315 0.005203536893129606 0.005620130698110045 0.005840772679515649 0.005840772679515649 0.005620130698110045 0.005203536893129606 0.00463582434788513 0.003974032527788299 0.0032780228327792802 0.002601768446564804 0.0019870162638941503 0.0014601931698789128 0.0010325124922636595 0.0007025163372637547 0.000459932027910791 0.0002897390217428208 0.0001756290843159388 0.00010243821352435276 5.749150348884896e-5 3.104712912334606e-5 1.6133007691619673e-5 8.066503845809865e-6 3.88089114041826e-6 1.7966094840265274e-6 8.002985431590045e-7 3.430255553045695e-7 1.4147413171036213e-7; 1.3098732099342573e-7 3.175986873248062e-7 7.40976183391281e-7 1.6634352891156783e-6 3.5932189680530496e-6 7.468571927403523e-6 1.4937143854807011e-5 2.8745751744424377e-5 5.322992925170177e-5 9.484495147408377e-5 0.0001626105279103003 0.0002682620333945352 0.0004258394340136128 0.000650442111641201 0.0009559772067076492 0.0013519559310455054 0.0018397281116431623 0.002408911611901339 0.0030350384471706845 0.0036794562232863232 0.004292192534312561 0.004817823223802676 0.005203536893129606 0.005407823724182016 0.005407823724182018 0.005203536893129606 0.004817823223802676 0.00429219253431256 0.0036794562232863245 0.0030350384471706845 0.0024089116119013388 0.0018397281116431631 0.0013519559310455054 0.0009559772067076492 0.0006504421116412 0.00042583943401361353 0.0002682620333945352 0.00016261052791030005 9.484495147408411e-5 5.322992925170177e-5 2.8745751744424377e-5 1.4937143854807011e-5 7.46857192740353e-6 3.5932189680530496e-6 1.6634352891156783e-6 7.40976183391281e-7 3.175986873248062e-7 1.3098732099342573e-7; 1.166964363656801e-7 2.8294826342072426e-7 6.60134731955813e-7 1.481952366782562e-6 3.2011941726360163e-6 6.653741156462725e-6 1.330748231292542e-5 2.5609553381088117e-5 4.742247573704209e-5 8.449724569034398e-5 0.0001448695108714104 0.00023899430167691231 0.0003793798058963353 0.0005794780434856414 0.0008516788680272265 0.0012044558059506696 0.0016390114163896401 0.002146096267156281 0.002703911862091008 0.0032780228327792794 0.003823908826830595 0.00429219253431256 0.00463582434788513 0.004817823223802676 0.004817823223802676 0.00463582434788513 0.00429219253431256 0.003823908826830594 0.0032780228327792807 0.002703911862091008 0.0021460962671562807 0.0016390114163896412 0.0012044558059506696 0.0008516788680272265 0.0005794780434856405 0.00037937980589633605 0.00023899430167691231 0.00014486951087141022 8.449724569034428e-5 4.742247573704209e-5 2.5609553381088117e-5 1.330748231292542e-5 6.653741156462732e-6 3.2011941726360163e-6 1.481952366782562e-6 6.60134731955813e-7 2.8294826342072426e-7 1.166964363656801e-7; 1.0003731789487592e-7 2.4255569627614246e-7 5.658965268414476e-7 1.2703947492992233e-6 2.7442044424365453e-6 5.703879569839515e-6 1.1407759139679001e-5 2.1953635539492345e-5 4.065263197757515e-5 7.243475543570508e-5 0.00012418851649338445 0.00020487642704870529 0.00032522105582060036 0.0004967540659735377 0.0007300965849394566 0.0010325124922636597 0.0014050326745275118 0.0018397281116431638 0.0023179121739425657 0.0028100653490550223 0.0032780228327792807 0.0036794562232863245 0.003974032527788299 0.004130049969054637 0.004130049969054637 0.003974032527788299 0.0036794562232863245 0.0032780228327792807 0.0028100653490550236 0.0023179121739425657 0.001839728111643163 0.0014050326745275127 0.0010325124922636597 0.0007300965849394566 0.0004967540659735367 0.00032522105582060085 0.00020487642704870529 0.00012418851649338431 7.243475543570533e-5 4.065263197757515e-5 2.1953635539492345e-5 1.1407759139679001e-5 5.703879569839521e-6 2.7442044424365453e-6 1.2703947492992233e-6 5.658965268414476e-7 2.4255569627614246e-7 1.0003731789487592e-7; 8.25168414944769e-8 2.0007463578975186e-7 4.667857454627196e-7 1.0478985679474042e-6 2.263586107365787e-6 4.704905491994815e-6 9.409810983989606e-6 1.810868885892628e-5 3.353275417431698e-5 5.9748575419228004e-5 0.00010243821352435257 0.00016899449138068818 0.0002682620333945348 0.00040975285409741024 0.0006022279029753348 0.0008516788680272268 0.0011589560869712824 0.0015175192235853427 0.001911954413415297 0.002317912173942565 0.002703911862091009 0.0030350384471706845 0.0032780228327792802 0.003406715472108905 0.003406715472108905 0.0032780228327792802 0.0030350384471706845 0.002703911862091008 0.0023179121739425657 0.001911954413415297 0.0015175192235853427 0.0011589560869712833 0.0008516788680272268 0.0006022279029753348 0.00040975285409740943 0.0002682620333945353 0.00016899449138068818 0.00010243821352435248 5.97485754192282e-5 3.353275417431698e-5 1.810868885892628e-5 9.409810983989606e-6 4.704905491994819e-6 2.263586107365787e-6 1.0478985679474042e-6 4.667857454627196e-7 2.0007463578975186e-7 8.25168414944769e-8; 6.549366049671285e-8 1.587993436624031e-7 3.704880916956404e-7 8.317176445578389e-7 1.7966094840265248e-6 3.734285963701764e-6 7.46857192740351e-6 1.43728758722122e-5 2.6614964625850903e-5 4.742247573704192e-5 8.130526395515017e-5 0.00013413101669726773 0.00021291971700680657 0.00032522105582060074 0.00047798860335382463 0.0006759779655227528 0.0009198640558215815 0.00120445580595067 0.0015175192235853427 0.0018397281116431623 0.002146096267156281 0.0024089116119013388 0.002601768446564804 0.00270391186209101 0.00270391186209101 0.002601768446564804 0.0024089116119013388 0.0021460962671562807 0.001839728111643163 0.0015175192235853427 0.0012044558059506696 0.000919864055821582 0.0006759779655227528 0.00047798860335382463 0.0003252210558206002 0.00021291971700680693 0.00013413101669726773 8.130526395515009e-5 4.742247573704209e-5 2.6614964625850903e-5 1.43728758722122e-5 7.46857192740351e-6 3.734285963701767e-6 1.7966094840265248e-6 8.317176445578389e-7 3.704880916956404e-7 1.587993436624031e-7 6.549366049671285e-8; 5.001865894743796e-8 1.212778481380712e-7 2.8294826342072373e-7 6.351973746496115e-7 1.3721022212182733e-6 2.85193978491976e-6 5.703879569839505e-6 1.0976817769746181e-5 2.032631598878759e-5 3.621737771785257e-5 6.209425824669228e-5 0.00010243821352435266 0.00016261052791030018 0.00024837703298676895 0.00036504829246972854 0.0005162562461318303 0.0007025163372637562 0.000919864055821582 0.0011589560869712833 0.0014050326745275118 0.0016390114163896412 0.0018397281116431631 0.0019870162638941503 0.00206502498452732 0.00206502498452732 0.0019870162638941503 0.0018397281116431631 0.0016390114163896412 0.0014050326745275127 0.0011589560869712833 0.000919864055821582 0.0007025163372637566 0.0005162562461318303 0.00036504829246972854 0.0002483770329867685 0.0001626105279103005 0.00010243821352435266 6.209425824669217e-5 3.621737771785269e-5 2.032631598878759e-5 1.0976817769746181e-5 5.703879569839505e-6 2.851939784919763e-6 1.3721022212182733e-6 6.351973746496115e-7 2.8294826342072373e-7 1.212778481380712e-7 5.001865894743796e-8; 3.67570741562095e-8 8.912311827874256e-8 2.079294111394597e-7 4.667857454627196e-7 1.0083129807262293e-6 2.095797135894812e-6 4.191594271789614e-6 8.066503845809836e-6 1.493714385480706e-5 2.661496462585082e-5 4.563103655871605e-5 7.527848787191693e-5 0.00011949715083845609 0.00018252414623486427 0.0002682620333945353 0.00037937980589633605 0.0005162562461318299 0.000675977965522753 0.0008516788680272268 0.0010325124922636595 0.0012044558059506696 0.0013519559310455054 0.0014601931698789128 0.0015175192235853429 0.0015175192235853429 0.0014601931698789128 0.0013519559310455054 0.0012044558059506696 0.0010325124922636597 0.0008516788680272268 0.0006759779655227528 0.0005162562461318303 0.00037937980589633605 0.0002682620333945353 0.00018252414623486392 0.0001194971508384563 7.527848787191693e-5 4.563103655871602e-5 2.6614964625850903e-5 1.493714385480706e-5 8.066503845809836e-6 4.191594271789614e-6 2.0957971358948143e-6 1.0083129807262293e-6 4.667857454627196e-7 2.079294111394597e-7 8.912311827874256e-8 3.67570741562095e-8; 2.599117639243255e-8 6.301956129538964e-8 1.4702829662483776e-7 3.3006736597790707e-7 7.129849462299381e-7 1.4819523667825646e-6 2.963904733565124e-6 5.7038795698395e-6 1.056215571129304e-5 1.8819621967979195e-5 3.2266015383239386e-5 5.322992925170173e-5 8.449724569034405e-5 0.00012906406153295754 0.00018968990294816794 0.0002682620333945353 0.0003650482924697282 0.00047798860335382485 0.0006022279029753348 0.0007300965849394566 0.0008516788680272268 0.0009559772067076492 0.0010325124922636595 0.0010730481335781406 0.0010730481335781406 0.0010325124922636595 0.0009559772067076492 0.0008516788680272265 0.0007300965849394566 0.0006022279029753348 0.00047798860335382463 0.00036504829246972854 0.0002682620333945353 0.00018968990294816794 0.0001290640615329573 8.449724569034421e-5 5.322992925170173e-5 3.226601538323935e-5 1.8819621967979263e-5 1.056215571129304e-5 5.7038795698395e-6 2.963904733565124e-6 1.4819523667825673e-6 7.129849462299381e-7 3.3006736597790707e-7 1.4702829662483776e-7 6.301956129538964e-8 2.599117639243255e-8; 1.768426646379523e-8 4.287819441307118e-8 1.0003731789487539e-7 2.2457618550331592e-7 4.851113925522823e-7 1.0083129807262312e-6 2.016625961452459e-6 3.880891140418253e-6 7.186437936106112e-6 1.2804776690544035e-5 2.1953635539492345e-5 3.621737771785257e-5 5.749150348884871e-5 8.78145421579694e-5 0.0001290640615329573 0.00018252414623486392 0.00024837703298676836 0.0003252210558206002 0.00040975285409740943 0.0004967540659735367 0.0005794780434856407 0.0006504421116412 0.0007025163372637547 0.0007300965849394552 0.0007300965849394552 0.0007025163372637547 0.0006504421116412 0.0005794780434856405 0.0004967540659735367 0.00040975285409740943 0.0003252210558206002 0.0002483770329867685 0.00018252414623486392 0.0001290640615329573 8.781454215796923e-5 5.749150348884881e-5 3.621737771785257e-5 2.195363553949232e-5 1.280477669054408e-5 7.186437936106112e-6 3.880891140418253e-6 2.016625961452459e-6 1.0083129807262312e-6 4.851113925522823e-7 2.2457618550331592e-7 1.0003731789487539e-7 4.287819441307118e-8 1.768426646379523e-8; 1.15777528654888e-8 2.807202318791463e-8 6.549366049671273e-8 1.4702829662483802e-7 3.1759868732480573e-7 6.601347319558153e-7 1.3202694639116285e-6 2.5407894985984416e-6 4.704905491994819e-6 8.383188543579222e-6 1.4372875872212226e-5 2.371123786852102e-5 3.7639243935958425e-5 5.7491503488848916e-5 8.449724569034421e-5 0.0001194971508384563 0.00016261052791030042 0.00021291971700680693 0.0002682620333945353 0.00032522105582060074 0.0003793798058963361 0.00042583943401361353 0.000459932027910791 0.00047798860335382485 0.00047798860335382485 0.000459932027910791 0.00042583943401361353 0.00037937980589633605 0.00032522105582060085 0.0002682620333945353 0.00021291971700680693 0.0001626105279103005 0.0001194971508384563 8.449724569034421e-5 5.749150348884881e-5 3.763924393595849e-5 2.371123786852102e-5 1.4372875872212212e-5 8.38318854357925e-6 4.704905491994819e-6 2.5407894985984416e-6 1.3202694639116285e-6 6.601347319558153e-7 3.1759868732480573e-7 1.4702829662483802e-7 6.549366049671273e-8 2.807202318791463e-8 1.15777528654888e-8; 7.293527272855018e-9 1.7684266463795293e-8 4.125842074723837e-8 9.262202292391026e-8 2.0007463578975112e-7 4.158588222789203e-7 8.317176445578389e-7 1.6005970863180064e-6 2.963904733565132e-6 5.2810778556465e-6 9.054344429463164e-6 1.4937143854807035e-5 2.371123786852098e-5 3.621737771785263e-5 5.322992925170173e-5 7.527848787191693e-5 0.00010243821352435266 0.00013413101669726773 0.00016899449138068818 0.00020487642704870515 0.00023899430167691242 0.0002682620333945352 0.0002897390217428208 0.0003011139514876675 0.0003011139514876675 0.0002897390217428208 0.0002682620333945352 0.00023899430167691231 0.00020487642704870529 0.00016899449138068818 0.00013413101669726773 0.00010243821352435266 7.527848787191693e-5 5.322992925170173e-5 3.621737771785257e-5 2.371123786852102e-5 1.4937143854807035e-5 9.054344429463147e-6 5.28107785564652e-6 2.963904733565132e-6 1.6005970863180064e-6 8.317176445578389e-7 4.1585882227892094e-7 2.0007463578975112e-7 9.262202292391026e-8 4.125842074723837e-8 1.7684266463795293e-8 7.293527272855018e-9; 4.421066615948815e-9 1.0719548603267791e-8 2.5009329473718887e-8 5.614404637582897e-8 1.2127784813807057e-7 2.5207824518155776e-7 5.041564903631147e-7 9.702227851045628e-7 1.796609484026531e-6 3.2011941726360103e-6 5.488408884873096e-6 9.054344429463147e-6 1.4372875872212187e-5 2.1953635539492362e-5 3.226601538323935e-5 4.563103655871602e-5 6.209425824669217e-5 8.130526395515009e-5 0.00010243821352435248 0.00012418851649338426 0.0001448695108714103 0.00016261052791030005 0.0001756290843159388 0.00018252414623486392 0.00018252414623486392 0.0001756290843159388 0.00016261052791030005 0.00014486951087141022 0.00012418851649338431 0.00010243821352435248 8.130526395515009e-5 6.209425824669217e-5 4.563103655871602e-5 3.226601538323935e-5 2.195363553949232e-5 1.4372875872212212e-5 9.054344429463147e-6 5.488408884873085e-6 3.2011941726360218e-6 1.796609484026531e-6 9.702227851045628e-7 5.041564903631147e-7 2.5207824518155823e-7 1.2127784813807057e-7 5.614404637582897e-8 2.5009329473718887e-8 1.0719548603267791e-8 4.421066615948815e-9; 2.5786512967024106e-9 6.252332368429754e-9 1.4587054545710036e-8 3.274683024835643e-8 7.073706585518093e-8 1.4702829662483828e-7 2.9405659324967556e-7 5.658965268414476e-7 1.047898567947408e-6 1.8671429818508787e-6 3.201194172636025e-6 5.28107785564652e-6 8.383188543579237e-6 1.2804776690544102e-5 1.8819621967979263e-5 2.6614964625850903e-5 3.621737771785266e-5 4.742247573704209e-5 5.97485754192282e-5 7.243475543570533e-5 8.449724569034428e-5 9.484495147408411e-5 0.00010243821352435276 0.00010645985850340356 0.00010645985850340356 0.00010243821352435276 9.484495147408411e-5 8.449724569034428e-5 7.243475543570533e-5 5.97485754192282e-5 4.742247573704209e-5 3.621737771785269e-5 2.6614964625850903e-5 1.8819621967979263e-5 1.280477669054408e-5 8.38318854357925e-6 5.28107785564652e-6 3.2011941726360218e-6 1.8671429818508857e-6 1.047898567947408e-6 5.658965268414476e-7 2.9405659324967556e-7 1.4702829662483855e-7 7.073706585518093e-8 3.274683024835643e-8 1.4587054545710036e-8 6.252332368429754e-9 2.5786512967024106e-9; 1.447219108186102e-9 3.509002898489335e-9 8.186707562089107e-9 1.837853707810478e-8 3.9699835915600776e-8 8.251684149447703e-8 1.650336829889538e-7 3.1759868732480573e-7 5.881131864993532e-7 1.0478985679474042e-6 1.796609484026531e-6 2.963904733565132e-6 4.70490549199481e-6 7.186437936106125e-6 1.056215571129304e-5 1.493714385480706e-5 2.0326315988787573e-5 2.6614964625850903e-5 3.353275417431698e-5 4.065263197757515e-5 4.742247573704209e-5 5.322992925170177e-5 5.749150348884896e-5 5.97485754192282e-5 5.97485754192282e-5 5.749150348884896e-5 5.322992925170177e-5 4.742247573704209e-5 4.065263197757515e-5 3.353275417431698e-5 2.6614964625850903e-5 2.032631598878759e-5 1.493714385480706e-5 1.056215571129304e-5 7.186437936106112e-6 4.704905491994819e-6 2.963904733565132e-6 1.796609484026531e-6 1.047898567947408e-6 5.881131864993532e-7 3.1759868732480573e-7 1.650336829889538e-7 8.251684149447703e-8 3.9699835915600776e-8 1.837853707810478e-8 8.186707562089107e-9 3.509002898489335e-9 1.447219108186102e-9; 7.815415460537175e-10 1.8949663771573584e-9 4.421066615948807e-9 9.924958978900158e-9 2.143909720653547e-8 4.456155913937119e-8 8.912311827874205e-8 1.7151277765228381e-7 3.1759868732480573e-7 5.658965268414456e-7 9.702227851045645e-7 1.6005970863180064e-6 2.5407894985984374e-6 3.88089114041826e-6 5.7038795698395e-6 8.066503845809836e-6 1.0976817769746173e-5 1.43728758722122e-5 1.810868885892628e-5 2.1953635539492345e-5 2.5609553381088117e-5 2.8745751744424377e-5 3.104712912334606e-5 3.2266015383239325e-5 3.2266015383239325e-5 3.104712912334606e-5 2.8745751744424377e-5 2.5609553381088117e-5 2.1953635539492345e-5 1.810868885892628e-5 1.43728758722122e-5 1.0976817769746181e-5 8.066503845809836e-6 5.7038795698395e-6 3.880891140418253e-6 2.5407894985984416e-6 1.6005970863180064e-6 9.702227851045628e-7 5.658965268414476e-7 3.1759868732480573e-7 1.7151277765228381e-7 8.912311827874205e-8 4.4561559139371265e-8 2.143909720653547e-8 9.924958978900158e-9 4.421066615948807e-9 1.8949663771573584e-9 7.815415460537175e-10; 4.061121311317591e-10 9.846806452404627e-10 2.297317134763089e-9 5.1573025934047955e-9 1.1140389784842775e-8 2.315550573097756e-8 4.631101146195504e-8 8.912311827874205e-8 1.650336829889538e-7 2.9405659324967445e-7 5.041564903631155e-7 8.317176445578389e-7 1.320269463911624e-6 2.016625961452459e-6 2.963904733565124e-6 4.191594271789614e-6 5.703879569839505e-6 7.46857192740351e-6 9.409810983989606e-6 1.1407759139679001e-5 1.330748231292542e-5 1.4937143854807011e-5 1.6133007691619673e-5 1.6766377087158443e-5 1.6766377087158443e-5 1.6133007691619673e-5 1.4937143854807011e-5 1.330748231292542e-5 1.1407759139679001e-5 9.409810983989606e-6 7.46857192740351e-6 5.703879569839505e-6 4.191594271789614e-6 2.963904733565124e-6 2.016625961452459e-6 1.3202694639116285e-6 8.317176445578389e-7 5.041564903631147e-7 2.9405659324967556e-7 1.650336829889538e-7 8.912311827874205e-8 4.631101146195504e-8 2.31555057309776e-8 1.1140389784842775e-8 5.1573025934047955e-9 2.297317134763089e-9 9.846806452404627e-10 4.061121311317591e-10; 2.0305606556587994e-10 4.923403226202322e-10 1.1486585673815487e-9 2.578651296702407e-9 5.570194892421407e-9 1.1577752865488842e-8 2.31555057309776e-8 4.4561559139371265e-8 8.251684149447703e-8 1.4702829662483802e-7 2.5207824518155866e-7 4.1585882227892094e-7 6.601347319558153e-7 1.0083129807262346e-6 1.4819523667825673e-6 2.0957971358948143e-6 2.851939784919763e-6 3.734285963701767e-6 4.704905491994819e-6 5.703879569839521e-6 6.653741156462732e-6 7.46857192740353e-6 8.066503845809865e-6 8.38318854357925e-6 8.38318854357925e-6 8.066503845809865e-6 7.46857192740353e-6 6.653741156462732e-6 5.703879569839521e-6 4.704905491994819e-6 3.734285963701767e-6 2.851939784919763e-6 2.0957971358948143e-6 1.4819523667825673e-6 1.0083129807262312e-6 6.601347319558153e-7 4.1585882227892094e-7 2.5207824518155823e-7 1.4702829662483855e-7 8.251684149447703e-8 4.4561559139371265e-8 2.31555057309776e-8 1.1577752865488842e-8 5.570194892421407e-9 2.578651296702407e-9 1.1486585673815487e-9 4.923403226202322e-10 2.0305606556587994e-10; 9.769269325671449e-11 2.368707971446701e-10 5.526333269935996e-10 1.2406198723625213e-9 2.6798871508169387e-9 5.5701948924213965e-9 1.1140389784842775e-8 2.143909720653547e-8 3.9699835915600776e-8 7.073706585518068e-8 1.2127784813807075e-7 2.0007463578975112e-7 3.1759868732480456e-7 4.851113925522823e-7 7.129849462299381e-7 1.0083129807262293e-6 1.3721022212182733e-6 1.7966094840265248e-6 2.263586107365787e-6 2.7442044424365453e-6 3.2011941726360163e-6 3.5932189680530496e-6 3.88089114041826e-6 4.033251922904918e-6 4.033251922904918e-6 3.88089114041826e-6 3.5932189680530496e-6 3.2011941726360163e-6 2.7442044424365453e-6 2.263586107365787e-6 1.7966094840265248e-6 1.3721022212182733e-6 1.0083129807262293e-6 7.129849462299381e-7 4.851113925522823e-7 3.1759868732480573e-7 2.0007463578975112e-7 1.2127784813807057e-7 7.073706585518093e-8 3.9699835915600776e-8 2.143909720653547e-8 1.1140389784842775e-8 5.570194892421407e-9 2.6798871508169387e-9 1.2406198723625213e-9 5.526333269935996e-10 2.368707971446701e-10 9.769269325671449e-11; 4.522559713081575e-11 1.0965634057779147e-10 2.5583461131528484e-10 5.743292836907731e-10 1.2406198723625213e-9 2.5786512967024023e-9 5.1573025934047955e-9 9.924958978900158e-9 1.837853707810478e-8 3.274683024835631e-8 5.6144046375829066e-8 9.262202292391026e-8 1.4702829662483752e-7 2.2457618550331592e-7 3.3006736597790707e-7 4.667857454627196e-7 6.351973746496115e-7 8.317176445578389e-7 1.0478985679474042e-6 1.2703947492992233e-6 1.481952366782562e-6 1.6634352891156783e-6 1.7966094840265274e-6 1.8671429818508787e-6 1.8671429818508787e-6 1.7966094840265274e-6 1.6634352891156783e-6 1.481952366782562e-6 1.2703947492992233e-6 1.0478985679474042e-6 8.317176445578389e-7 6.351973746496115e-7 4.667857454627196e-7 3.3006736597790707e-7 2.2457618550331592e-7 1.4702829662483802e-7 9.262202292391026e-8 5.614404637582897e-8 3.274683024835643e-8 1.837853707810478e-8 9.924958978900158e-9 5.1573025934047955e-9 2.578651296702407e-9 1.2406198723625213e-9 5.743292836907731e-10 2.5583461131528484e-10 1.0965634057779147e-10 4.522559713081575e-11; 2.0145713255487537e-11 4.884634662835752e-11 1.1396136363835936e-10 2.5583461131528484e-10 5.526333269935996e-10 1.1486585673815487e-9 2.297317134763089e-9 4.421066615948807e-9 8.186707562089107e-9 1.4587054545709983e-8 2.5009329473718933e-8 4.125842074723837e-8 6.549366049671273e-8 1.0003731789487575e-7 1.4702829662483776e-7 2.079294111394597e-7 2.8294826342072373e-7 3.704880916956411e-7 4.667857454627196e-7 5.658965268414476e-7 6.60134731955813e-7 7.40976183391281e-7 8.002985431590045e-7 8.317176445578389e-7 8.317176445578389e-7 8.002985431590045e-7 7.40976183391281e-7 6.60134731955813e-7 5.658965268414476e-7 4.667857454627196e-7 3.704880916956404e-7 2.8294826342072373e-7 2.079294111394597e-7 1.4702829662483776e-7 1.0003731789487539e-7 6.549366049671273e-8 4.125842074723837e-8 2.5009329473718887e-8 1.4587054545710036e-8 8.186707562089107e-9 4.421066615948807e-9 2.297317134763089e-9 1.1486585673815487e-9 5.526333269935996e-10 2.5583461131528484e-10 1.1396136363835936e-10 4.884634662835752e-11 2.0145713255487537e-11; 8.634895734275082e-12 2.0936618365757505e-11 4.884634662835752e-11 1.0965634057779147e-10 2.368707971446701e-10 4.923403226202322e-10 9.846806452404627e-10 1.8949663771573584e-9 3.509002898489335e-9 6.252332368429731e-9 1.0719548603267813e-8 1.7684266463795293e-8 2.8072023187914533e-8 4.287819441307118e-8 6.301956129538964e-8 8.912311827874256e-8 1.212778481380712e-7 1.587993436624031e-7 2.0007463578975186e-7 2.4255569627614246e-7 2.8294826342072426e-7 3.175986873248062e-7 3.430255553045695e-7 3.564924731149702e-7 3.564924731149702e-7 3.430255553045695e-7 3.175986873248062e-7 2.8294826342072426e-7 2.4255569627614246e-7 2.0007463578975186e-7 1.587993436624031e-7 1.212778481380712e-7 8.912311827874256e-8 6.301956129538964e-8 4.287819441307118e-8 2.807202318791463e-8 1.7684266463795293e-8 1.0719548603267791e-8 6.252332368429754e-9 3.509002898489335e-9 1.8949663771573584e-9 9.846806452404627e-10 4.923403226202322e-10 2.368707971446701e-10 1.0965634057779147e-10 4.884634662835752e-11 2.0936618365757505e-11 8.634895734275082e-12; 3.5612926136987482e-12 8.634895734275082e-12 2.0145713255487537e-11 4.522559713081575e-11 9.769269325671449e-11 2.0305606556587994e-10 4.061121311317591e-10 7.815415460537175e-10 1.447219108186102e-9 2.5786512967024023e-9 4.4210666159488215e-9 7.293527272855018e-9 1.15777528654888e-8 1.7684266463795293e-8 2.599117639243255e-8 3.67570741562095e-8 5.001865894743796e-8 6.549366049671296e-8 8.25168414944769e-8 1.0003731789487592e-7 1.166964363656801e-7 1.3098732099342573e-7 1.4147413171036213e-7 1.4702829662483802e-7 1.4702829662483802e-7 1.4147413171036213e-7 1.3098732099342573e-7 1.166964363656801e-7 1.0003731789487592e-7 8.25168414944769e-8 6.549366049671285e-8 5.001865894743796e-8 3.67570741562095e-8 2.599117639243255e-8 1.768426646379523e-8 1.15777528654888e-8 7.293527272855018e-9 4.421066615948815e-9 2.5786512967024106e-9 1.447219108186102e-9 7.815415460537175e-10 4.061121311317591e-10 2.0305606556587994e-10 9.769269325671449e-11 4.522559713081575e-11 2.0145713255487537e-11 8.634895734275082e-12 3.5612926136987482e-12])

To make the Gaussian Markov random field efficient we first precompute a bunch of quantities that allow us to scale things linearly with the number of image pixels. The returns a functional that accepts a single argument related to the correlation length of the field. The second argument defines the underlying random field of the Markov process. Here we are using a zero mean and unit variance Gaussian Markov random field. For this tutorial we will use the first order random field

julia

cprior = corr_image_prior(grid, dvis)

HierarchicalPrior(
	map: 
	ConditionalMarkov(
Random Field: VLBIImagePriors.GaussMarkovRandomField
Graph: MarkovRandomFieldGraph{1}(
dims: (48, 48)
)
)	hyper prior: 
	Truncated(Distributions.InverseGamma{Float64}(
invd: Distributions.Gamma{Float64}(α=1.0, θ=0.03623874772901475)
θ: 27.594772520225487
)
; lower=1.0, upper=96.0)

)

Putting everything together the total prior is then our image prior, a prior on the standard deviation of the MRF, and a prior on the fractional flux of the Gaussian component.

julia

prior = (
    c = cprior,
    σimg = truncated(Normal(0.0, 0.5); lower = 0.0),
    fg = Uniform(0.0, 1.0),
)

(c = HierarchicalPrior(
	map: 
	ConditionalMarkov(
Random Field: VLBIImagePriors.GaussMarkovRandomField
Graph: MarkovRandomFieldGraph{1}(
dims: (48, 48)
)
)	hyper prior: 
	Truncated(Distributions.InverseGamma{Float64}(
invd: Distributions.Gamma{Float64}(α=1.0, θ=0.03623874772901475)
θ: 27.594772520225487
)
; lower=1.0, upper=96.0)

)
, σimg = Truncated(Distributions.Normal{Float64}(μ=0.0, σ=0.5); lower=0.0), fg = Distributions.Uniform{Float64}(a=0.0, b=1.0))

Now we can construct our sky model.

julia

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

SkyModel
  with map: sky
   on grid: 
RectiGrid(
executor: ComradeBase.Serial()
Dimensions: 
(↓ X Sampled{Float64} LinRange{Float64}(-4.747133960864206e-10, 4.747133960864206e-10, 48) ForwardOrdered Regular Points,
→ Y Sampled{Float64} LinRange{Float64}(-4.747133960864206e-10, 4.747133960864206e-10, 48) ForwardOrdered Regular Points)
)
   )

Unlike other imaging examples (e.g., Imaging a Black Hole using only Closure Quantities) we also need to include a model for the instrument, i.e., gains. The gains will be broken into two components

Gain amplitudes which are typically known to 10-20%, except for LMT, which has amplitudes closer to 50-100%.
Gain phases which are more difficult to constrain and can shift rapidly.

julia

G = SingleStokesGain() do x
    lg = x.lg
    gp = x.gp
    return exp(lg + 1im * gp)
end

intpr = (
    lg = ArrayPrior(IIDSitePrior(ScanSeg(), Normal(0.0, 0.2)); LM = IIDSitePrior(ScanSeg(), Normal(0.0, 1.0))),
    gp = ArrayPrior(IIDSitePrior(ScanSeg(), DiagonalVonMises(0.0, inv(π^2))); refant = SEFDReference(0.0), phase = true),
)
intmodel = InstrumentModel(G, intpr)

InstrumentModel
  with Jones: SingleStokesGain
  with reference basis: PolarizedTypes.CirBasis()

To form the posterior we just combine the skymodel, instrument model and the data. To utilize gradients of the posterior we also need to load Enzyme.jl. Under the hood, Comrade will use Enzyme to compute the gradients of the posterior.

julia

using Enzyme
post = VLBIPosterior(skym, intmodel, dvis)

VLBIPosterior
ObservedSkyModel
  with map: sky
   on grid: 
FourierDualDomain(
Algorithm: VLBISkyModels.NFFTAlg{Float64, AbstractNFFTs.PrecomputeFlags, UInt32}(1, 1.0e-9, AbstractNFFTs.TENSOR, 0x00000000)
Image Domain: RectiGrid(
executor: ComradeBase.Serial()
Dimensions: 
(↓ X Sampled{Float64} LinRange{Float64}(-4.747133960864206e-10, 4.747133960864206e-10, 48) ForwardOrdered Regular Points,
→ Y Sampled{Float64} LinRange{Float64}(-4.747133960864206e-10, 4.747133960864206e-10, 48) ForwardOrdered Regular Points)
)
Visibility Domain: UnstructredDomain(
executor: ComradeBase.Serial()
Dimensions: 
274-element StructArray(::Vector{Float64}, ::Vector{Float64}, ::Vector{Float64}, ::Vector{Float64}) with eltype @NamedTuple{U::Float64, V::Float64, Ti::Float64, Fr::Float64}:
 (U = -4.405690154666661e9, V = -4.523017159111106e9, Ti = 0.9166666567325592, Fr = 2.27070703125e11)
 (U = 787577.6145833326, V = -1.6838098888888871e6, Ti = 1.2166666388511658, Fr = 2.27070703125e11)
 (U = -4.444299918222218e9, V = -4.597825294222218e9, Ti = 1.2166666388511658, Fr = 2.27070703125e11)
 (U = 1.337045162666665e9, V = -3.765300401777774e9, Ti = 1.2166666388511658, Fr = 2.27070703125e11)
 (U = -1.336260540444443e9, V = 3.763616127999996e9, Ti = 1.2166666388511658, Fr = 2.27070703125e11)
 (U = 4.445088654222218e9, V = 4.596145080888884e9, Ti = 1.2166666388511658, Fr = 2.27070703125e11)
 (U = 5.781345607111105e9, V = 8.325259893333325e8, Ti = 1.2166666388511658, Fr = 2.27070703125e11)
 (U = 757554.6649305547, V = -1.6707483020833314e6, Ti = 1.516666665673256, Fr = 2.27070703125e11)
 (U = 1.4806382151111097e9, V = -3.741479615999996e9, Ti = 1.516666665673256, Fr = 2.27070703125e11)
 (U = -4.455366328888884e9, V = -4.673060451555551e9, Ti = 1.516666665673256, Fr = 2.27070703125e11)
 (U = -1.4798758791111097e9, V = 3.739809735111107e9, Ti = 1.516666665673256, Fr = 2.27070703125e11)
 (U = 4.456123861333328e9, V = 4.671391715555551e9, Ti = 1.516666665673256, Fr = 2.27070703125e11)
 (U = 5.936013027555549e9, V = 9.315912497777768e8, Ti = 1.516666665673256, Fr = 2.27070703125e11)
 (U = 722830.7065972214, V = -1.6582321493055536e6, Ti = 1.816666603088379, Fr = 2.27070703125e11)
 (U = -4.438811278222218e9, V = -4.748261176888884e9, Ti = 1.816666603088379, Fr = 2.27070703125e11)
 (U = 1.615060401777776e9, V = -3.715306943999996e9, Ti = 1.816666603088379, Fr = 2.27070703125e11)
 (U = -1.6143345528888872e9, V = 3.713649343999996e9, Ti = 1.816666603088379, Fr = 2.27070703125e11)
 (U = 4.439536184888884e9, V = 4.746598001777773e9, Ti = 1.816666603088379, Fr = 2.27070703125e11)
 (U = 6.053865329777771e9, V = 1.0329465084444433e9, Ti = 1.816666603088379, Fr = 2.27070703125e11)
 (U = 683620.5937499993, V = -1.6463396631944429e6, Ti = 2.1166666746139526, Fr = 2.27070703125e11)
 (U = -4.394744988444439e9, V = -4.822934485333328e9, Ti = 2.1166666746139526, Fr = 2.27070703125e11)
 (U = 1.7394626631111093e9, V = -3.686947740444441e9, Ti = 2.1166666746139526, Fr = 2.27070703125e11)
 (U = -1.7387777884444425e9, V = 3.6853017173333297e9, Ti = 2.1166666746139526, Fr = 2.27070703125e11)
 (U = 4.395427669333329e9, V = 4.821289443555551e9, Ti = 2.1166666746139526, Fr = 2.27070703125e11)
 (U = 6.134203007999993e9, V = 1.1359907128888876e9, Ti = 2.1166666746139526, Fr = 2.27070703125e11)
 (U = -1.8643712995555537e9, V = 3.651452543999996e9, Ti = 2.449999988079071, Fr = 2.27070703125e11)
 (U = 6.178850375111104e9, V = 1.2516714417777765e9, Ti = 2.449999988079071, Fr = 2.27070703125e11)
 (U = 4.314477112888885e9, V = 4.903116785777773e9, Ti = 2.449999988079071, Fr = 2.27070703125e11)
 (U = 587276.6180555549, V = -1.6236182569444429e6, Ti = 2.7500000596046448, Fr = 2.27070703125e11)
 (U = 1.9658264106666646e9, V = -3.6206990791111073e9, Ti = 2.7500000596046448, Fr = 2.27070703125e11)
 (U = -4.212778346666662e9, V = -4.976839239111106e9, Ti = 2.7500000596046448, Fr = 2.27070703125e11)
 (U = -1.965238759111109e9, V = 3.6190756053333297e9, Ti = 2.7500000596046448, Fr = 2.27070703125e11)
 (U = 6.178607487999993e9, V = 1.356139594666665e9, Ti = 2.7500000596046448, Fr = 2.27070703125e11)
 (U = 4.213364039111107e9, V = 4.975216568888884e9, Ti = 2.7500000596046448, Fr = 2.27070703125e11)
 (U = 535806.2881944438, V = -1.6141240694444429e6, Ti = 3.0500001311302185, Fr = 2.27070703125e11)
 (U = 2.0544588266666646e9, V = -3.586710200888885e9, Ti = 3.0500001311302185, Fr = 2.27070703125e11)
 (U = -4.085601863111107e9, V = -5.046995640888884e9, Ti = 3.0500001311302185, Fr = 2.27070703125e11)
 (U = -2.0539194026666646e9, V = 3.5850976213333297e9, Ti = 3.0500001311302185, Fr = 2.27070703125e11)
 (U = 4.0861423288888845e9, V = 5.045379299555551e9, Ti = 3.0500001311302185, Fr = 2.27070703125e11)
 (U = 6.140063729777771e9, V = 1.4602817422222207e9, Ti = 3.0500001311302185, Fr = 2.27070703125e11)
 (U = 481010.2230902773, V = -1.6055275520833316e6, Ti = 3.3499998450279236, Fr = 2.27070703125e11)
 (U = 2.130349845333331e9, V = -3.5513328995555515e9, Ti = 3.3499998450279236, Fr = 2.27070703125e11)
 (U = -3.933104519111107e9, V = -5.114784767999995e9, Ti = 3.3499998450279236, Fr = 2.27070703125e11)
 (U = -2.1298682808888867e9, V = 3.5497276302222185e9, Ti = 3.3499998450279236, Fr = 2.27070703125e11)
 (U = 6.063457521777771e9, V = 1.5634511359999983e9, Ti = 3.3499998450279236, Fr = 2.27070703125e11)
 (U = 3.9335861404444404e9, V = 5.113178993777773e9, Ti = 3.3499998450279236, Fr = 2.27070703125e11)
 (U = 416648.93055555515, V = -1.5970948472222206e6, Ti = 3.6833333373069763, Fr = 2.27070703125e11)
 (U = 2.578938929777775e9, V = -4.734788238222218e9, Ti = 3.6833333373069763, Fr = 2.27070703125e11)
 (U = -3.735105663999996e9, V = -5.186825713777772e9, Ti = 3.6833333373069763, Fr = 2.27070703125e11)
 (U = 2.1991734328888865e9, V = -3.5106569031111073e9, Ti = 3.6833333373069763, Fr = 2.27070703125e11)
 (U = -2.5785191964444413e9, V = 4.733192191999995e9, Ti = 3.6833333373069763, Fr = 2.27070703125e11)
 (U = -2.198756124444442e9, V = 3.509060273777774e9, Ti = 3.6833333373069763, Fr = 2.27070703125e11)
 (U = 3.797641537777774e8, V = -1.2241308906666653e9, Ti = 3.6833333373069763, Fr = 2.27070703125e11)
 (U = 3.735526620444441e9, V = 5.185227192888883e9, Ti = 3.6833333373069763, Fr = 2.27070703125e11)
 (U = 6.314048611555549e9, V = 4.520333351111106e8, Ti = 3.6833333373069763, Fr = 2.27070703125e11)
 (U = 5.934280931555549e9, V = 1.6761668835555537e9, Ti = 3.6833333373069763, Fr = 2.27070703125e11)
 (U = 364284.04464285675, V = -1.591376535714284e6, Ti = 3.98333340883255, Fr = 2.27070703125e11)
 (U = -3.5324296035555515e9, V = -5.248264732444439e9, Ti = 3.98333340883255, Fr = 2.27070703125e11)
 (U = 2.246733795555553e9, V = -3.4730721351111073e9, Ti = 3.98333340883255, Fr = 2.27070703125e11)
 (U = 2.7040352142222195e9, V = -4.690125752888884e9, Ti = 3.98333340883255, Fr = 2.27070703125e11)
 (U = -2.687901513142854e9, V = 4.694663167999995e9, Ti = 3.98333340883255, Fr = 2.27070703125e11)
 (U = -2.2408562468571405e9, V = 3.476569380571425e9, Ti = 3.98333340883255, Fr = 2.27070703125e11)
 (U = 4.573011991111106e8, V = -1.2170540408888876e9, Ti = 3.98333340883255, Fr = 2.27070703125e11)
 (U = 3.561408036571425e9, V = 5.238641298285708e9, Ti = 3.98333340883255, Fr = 2.27070703125e11)
 (U = 5.779150734222217e9, V = 1.7751998862222204e9, Ti = 3.98333340883255, Fr = 2.27070703125e11)
 (U = 6.236464440888882e9, V = 5.581377528888882e8, Ti = 3.98333340883255, Fr = 2.27070703125e11)
 (U = 292810.0842803027, V = -1.5850529772727257e6, Ti = 4.283333480358124, Fr = 2.27070703125e11)
 (U = 2.81236712533333e9, V = -4.643490716444439e9, Ti = 4.283333480358124, Fr = 2.27070703125e11)
 (U = 2.280368668444442e9, V = -3.43479938133333e9, Ti = 4.283333480358124, Fr = 2.27070703125e11)
 (U = -3.307854179555552e9, V = -5.306092316444439e9, Ti = 4.283333480358124, Fr = 2.27070703125e11)
 (U = -2.812399111529409e9, V = 4.641748931764702e9, Ti = 4.283333480358124, Fr = 2.27070703125e11)
 (U = -2.280159088941174e9, V = 3.43308967905882e9, Ti = 4.283333480358124, Fr = 2.27070703125e11)
 (U = 5.3200019644444394e8, V = -1.2086904924444432e9, Ti = 4.283333480358124, Fr = 2.27070703125e11)
 (U = 3.307383679999996e9, V = 5.304687194352936e9, Ti = 4.283333480358124, Fr = 2.27070703125e11)
 (U = 6.120220771555549e9, V = 6.626024835555549e8, Ti = 4.283333480358124, Fr = 2.27070703125e11)
 (U = 5.588211185777772e9, V = 1.8712975182222202e9, Ti = 4.283333480358124, Fr = 2.27070703125e11)
 (U = 238138.85156249977, V = -1.5812476874999984e6, Ti = 4.583333194255829, Fr = 2.27070703125e11)
 (U = 5.092881095111105e9, V = -4.199598421333329e9, Ti = 4.583333194255829, Fr = 2.27070703125e11)
 (U = 2.903268209777775e9, V = -4.595170062222218e9, Ti = 4.583333194255829, Fr = 2.27070703125e11)
 (U = 2.299867178666664e9, V = -3.396077226666663e9, Ti = 4.583333194255829, Fr = 2.27070703125e11)
 (U = -3.062767004444441e9, V = -5.359950862222217e9, Ti = 4.583333194255829, Fr = 2.27070703125e11)
 (U = 5.050181503999994e9, V = -4.210810559999995e9, Ti = 4.583333194255829, Fr = 2.27070703125e11)
 (U = -2.8907034879999967e9, V = 4.600886271999995e9, Ti = 4.583333194255829, Fr = 2.27070703125e11)
 (U = 2.189606492444442e9, V = 3.955690115555551e8, Ti = 4.583333194255829, Fr = 2.27070703125e11)
 (U = -2.2976930559999976e9, V = 3.400280575999996e9, Ti = 4.583333194255829, Fr = 2.27070703125e11)
 (U = 2.79301371733333e9, V = -8.035218719999992e8, Ti = 4.583333194255829, Fr = 2.27070703125e11)
 (U = 6.034046399999994e8, V = -1.199091349333332e9, Ti = 4.583333194255829, Fr = 2.27070703125e11)
 (U = 3.1008046079999967e9, V = 5.350639359999994e9, Ti = 4.583333194255829, Fr = 2.27070703125e11)
 (U = 8.155639423999991e9, V = 1.1603480675555544e9, Ti = 4.583333194255829, Fr = 2.27070703125e11)
 (U = 5.966032455111105e9, V = 7.64783736888888e8, Ti = 4.583333194255829, Fr = 2.27070703125e11)
 (U = 5.36262916266666e9, V = 1.963875015111109e9, Ti = 4.583333194255829, Fr = 2.27070703125e11)
 (U = 163232.87499999983, V = -1.5773965972222206e6, Ti = 4.916666686534882, Fr = 2.27070703125e11)
 (U = 5.388079544888883e9, V = -4.101135303111107e9, Ti = 4.916666686534882, Fr = 2.27070703125e11)
 (U = 2.3048084479999976e9, V = -3.3528182257777743e9, Ti = 4.916666686534882, Fr = 2.27070703125e11)
 (U = -2.768253631999997e9, V = -5.414730879999994e9, Ti = 4.916666686534882, Fr = 2.27070703125e11)
 (U = 5.388055679999994e9, V = -4.101249678222218e9, Ti = 4.916666686534882, Fr = 2.27070703125e11)
 (U = 2.9831417315555525e9, V = -4.539868501333328e9, Ti = 4.916666686534882, Fr = 2.27070703125e11)
 (U = 5.356877311999994e9, V = -4.110922695111107e9, Ti = 4.916666686534882, Fr = 2.27070703125e11)
 (U = -2.975202047999997e9, V = 4.544589937777773e9, Ti = 4.916666686534882, Fr = 2.27070703125e11)
 (U = 2.4049356302222195e9, V = 4.387324826666662e8, Ti = 4.916666686534882, Fr = 2.27070703125e11)
 (U = 2.4049011413333306e9, V = 4.3861474933333284e8, Ti = 4.916666686534882, Fr = 2.27070703125e11)
 (U = -2.305032220444442e9, V = 3.356104334222219e9, Ti = 4.916666686534882, Fr = 2.27070703125e11)
 (U = 3.083277411555552e9, V = -7.48315953777777e8, Ti = 4.916666686534882, Fr = 2.27070703125e11)
 (U = 3.0832264319999967e9, V = -7.484370453333325e8, Ti = 4.916666686534882, Fr = 2.27070703125e11)
 (U = 6.783317208888881e8, V = -1.187050453333332e9, Ti = 4.916666686534882, Fr = 2.27070703125e11)
 (U = 2.8027259164444413e9, V = 5.407289457777772e9, Ti = 4.916666686534882, Fr = 2.27070703125e11)
 (U = 8.156328149333324e9, V = 1.3135985066666653e9, Ti = 4.916666686534882, Fr = 2.27070703125e11)
 (U = 8.156297400888881e9, V = 1.3134812622222207e9, Ti = 4.916666686534882, Fr = 2.27070703125e11)
 (U = 5.751382712888883e9, V = 8.74867279999999e8, Ti = 4.916666686534882, Fr = 2.27070703125e11)
 (U = 5.073054606222218e9, V = 2.0619144959999979e9, Ti = 4.916666686534882, Fr = 2.27070703125e11)
 (U = -5.356795278222217e9, V = 4.1110577777777734e9, Ti = 4.916666686534882, Fr = 2.27070703125e11)
 (U = 30750.899088541635, V = 116905.95138888876, Ti = 4.916666686534882, Fr = 2.27070703125e11)
 (U = 99881.10751488084, V = -1.575297547619046e6, Ti = 5.183333337306976, Fr = 2.27070703125e11)
 (U = 5.59465437866666e9, V = -4.018611114666662e9, Ti = 5.183333337306976, Fr = 2.27070703125e11)
 (U = 3.030624810666663e9, V = -4.494683207111106e9, Ti = 5.183333337306976, Fr = 2.27070703125e11)
 (U = 2.2960602239999976e9, V = -3.3182473599999967e9, Ti = 5.183333337306976, Fr = 2.27070703125e11)
 (U = 5.594626830222217e9, V = -4.0187286115555515e9, Ti = 5.183333337306976, Fr = 2.27070703125e11)
 (U = -2.5172720142222195e9, V = -5.45444669155555e9, Ti = 5.183333337306976, Fr = 2.27070703125e11)
 (U = 5.579493814857137e9, V = -4.023605065142853e9, Ti = 5.183333337306976, Fr = 2.27070703125e11)
 (U = -3.0274068479999967e9, V = 4.496662674285709e9, Ti = 5.183333337306976, Fr = 2.27070703125e11)
 (U = 2.5640353279999976e9, V = 4.7607165511111057e8, Ti = 5.183333337306976, Fr = 2.27070703125e11)
 (U = 2.564000689777775e9, V = 4.7595243111111057e8, Ti = 5.183333337306976, Fr = 2.27070703125e11)
 (U = -2.2970926933333306e9, V = 3.3193651078095202e9, Ti = 5.183333337306976, Fr = 2.27070703125e11)
 (U = 3.2985927679999967e9, V = -7.003643235555549e8, Ti = 5.183333337306976, Fr = 2.27070703125e11)
 (U = 7.345665155555549e8, V = -1.1764336924444432e9, Ti = 5.183333337306976, Fr = 2.27070703125e11)
 (U = 3.2985603128888855e9, V = -7.004829137777771e8, Ti = 5.183333337306976, Fr = 2.27070703125e11)
 (U = 2.5374800944761877e9, V = 5.44992148723809e9, Ti = 5.183333337306976, Fr = 2.27070703125e11)
 (U = 8.111924807111103e9, V = 1.4358354488888874e9, Ti = 5.183333337306976, Fr = 2.27070703125e11)
 (U = 4.813340984888884e9, V = 2.136196949333331e9, Ti = 5.183333337306976, Fr = 2.27070703125e11)
 (U = 5.547893447111105e9, V = 9.597677279999989e8, Ti = 5.183333337306976, Fr = 2.27070703125e11)
 (U = 8.111898709333324e9, V = 1.4357162488888874e9, Ti = 5.183333337306976, Fr = 2.27070703125e11)
 (U = -5.579455049142851e9, V = 4.0237272990476146e9, Ti = 5.183333337306976, Fr = 2.27070703125e11)
 (U = 27098.63357204858, V = 117340.62348090265, Ti = 5.183333337306976, Fr = 2.27070703125e11)
 (U = 45690.53906249995, V = -1.5743099999999984e6, Ti = 5.449999988079071, Fr = 2.27070703125e11)
 (U = 5.773825464888883e9, V = -3.933187214222218e9, Ti = 5.449999988079071, Fr = 2.27070703125e11)
 (U = 2.2760630257777753e9, V = -3.283890659555552e9, Ti = 5.449999988079071, Fr = 2.27070703125e11)
 (U = 3.0632664319999967e9, V = -4.448890581333328e9, Ti = 5.449999988079071, Fr = 2.27070703125e11)
 (U = 5.773791857777772e9, V = -3.9333102648888845e9, Ti = 5.449999988079071, Fr = 2.27070703125e11)
 (U = -2.253957283555553e9, V = -5.490298367999994e9, Ti = 5.449999988079071, Fr = 2.27070703125e11)
 (U = 5.744240127999994e9, V = -3.946954495999996e9, Ti = 5.449999988079071, Fr = 2.27070703125e11)
 (U = -3.0586252799999967e9, V = 4.455449599999995e9, Ti = 5.449999988079071, Fr = 2.27070703125e11)
 (U = 2.710558015999997e9, V = 5.1570311111111057e8, Ti = 5.449999988079071, Fr = 2.27070703125e11)
 (U = 2.7105363982222195e9, V = 5.1558591999999946e8, Ti = 5.449999988079071, Fr = 2.27070703125e11)
 (U = -2.2804378879999976e9, V = 3.2883723519999967e9, Ti = 5.449999988079071, Fr = 2.27070703125e11)
 (U = 3.4977656675555515e9, V = -6.492943324444438e8, Ti = 5.449999988079071, Fr = 2.27070703125e11)
 (U = 7.872002471111102e8, V = -1.164999438222221e9, Ti = 5.449999988079071, Fr = 2.27070703125e11)
 (U = 3.497740145777774e9, V = -6.494127164444437e8, Ti = 5.449999988079071, Fr = 2.27070703125e11)
 (U = 2.3015257599999976e9, V = 5.482690303999994e9, Ti = 5.449999988079071, Fr = 2.27070703125e11)
 (U = 8.027784391111103e9, V = 1.5571086684444427e9, Ti = 5.449999988079071, Fr = 2.27070703125e11)
 (U = 5.317234005333327e9, V = 1.0414041991111101e9, Ti = 5.449999988079071, Fr = 2.27070703125e11)
 (U = 8.027758634666658e9, V = 1.556993443555554e9, Ti = 5.449999988079071, Fr = 2.27070703125e11)
 (U = 4.530021788444439e9, V = 2.206406378666664e9, Ti = 5.449999988079071, Fr = 2.27070703125e11)
 (U = -5.744238079999993e9, V = 3.9470612479999957e9, Ti = 5.449999988079071, Fr = 2.27070703125e11)
 (U = 23313.64084201386, V = 117719.41384548598, Ti = 5.449999988079071, Fr = 2.27070703125e11)
 (U = 5.924709319111105e9, V = -3.8452848639999957e9, Ti = 5.716666638851166, Fr = 2.27070703125e11)
 (U = 5.924688497777772e9, V = -3.845403854222218e9, Ti = 5.716666638851166, Fr = 2.27070703125e11)
 (U = 3.08089890133333e9, V = -4.402725802666661e9, Ti = 5.716666638851166, Fr = 2.27070703125e11)
 (U = 2.244918698666664e9, V = -3.2499217208888855e9, Ti = 5.716666638851166, Fr = 2.27070703125e11)
 (U = 2.8438056746666636e9, V = 5.57437258666666e8, Ti = 5.716666638851166, Fr = 2.27070703125e11)
 (U = 2.8437909475555525e9, V = 5.573208195555549e8, Ti = 5.716666638851166, Fr = 2.27070703125e11)
 (U = 3.6797841706666627e9, V = -5.953657653333327e8, Ti = 5.716666638851166, Fr = 2.27070703125e11)
 (U = 3.679770524444441e9, V = -5.954820355555549e8, Ti = 5.716666638851166, Fr = 2.27070703125e11)
 (U = 8.35980113777777e8, V = -1.1528026026666653e9, Ti = 5.716666638851166, Fr = 2.27070703125e11)
 (U = 19414.168077256923, V = 118040.4856770832, Ti = 5.716666638851166, Fr = 2.27070703125e11)
 (U = -92093.64322916657, V = -1.5751088611111094e6, Ti = 6.049999952316284, Fr = 2.27070703125e11)
 (U = 3.081715143111108e9, V = -4.344833365333329e9, Ti = 6.049999952316284, Fr = 2.27070703125e11)
 (U = 6.07242564266666e9, V = -3.732578396444441e9, Ti = 6.049999952316284, Fr = 2.27070703125e11)
 (U = 2.190544647111109e9, V = -3.208251676444441e9, Ti = 6.049999952316284, Fr = 2.27070703125e11)
 (U = 6.05836305066666e9, V = -3.7438256071111073e9, Ti = 6.049999952316284, Fr = 2.27070703125e11)
 (U = -3.0829770239999967e9, V = 4.34976796444444e9, Ti = 6.049999952316284, Fr = 2.27070703125e11)
 (U = 2.99071138133333e9, V = 6.122509119999994e8, Ti = 6.049999952316284, Fr = 2.27070703125e11)
 (U = -2.1976609848888865e9, V = 3.211308942222219e9, Ti = 6.049999952316284, Fr = 2.27070703125e11)
 (U = 3.881882190222218e9, V = -5.2432489066666615e8, Ti = 6.049999952316284, Fr = 2.27070703125e11)
 (U = 8.911678026666657e8, V = -1.1365769386666653e9, Ti = 6.049999952316284, Fr = 2.27070703125e11)
 (U = -150875.45052083317, V = -1.5769051249999984e6, Ti = 6.316666603088379, Fr = 2.27070703125e11)
 (U = 6.157184056888882e9, V = -3.6406887111111073e9, Ti = 6.316666603088379, Fr = 2.27070703125e11)
 (U = 6.157173617777771e9, V = -3.6408075306666627e9, Ti = 6.316666603088379, Fr = 2.27070703125e11)
 (U = 2.1349524657777758e9, V = -3.175750087111108e9, Ti = 6.316666603088379, Fr = 2.27070703125e11)
 (U = 3.065380067555552e9, V = -4.298637105777774e9, Ti = 6.316666603088379, Fr = 2.27070703125e11)
 (U = 6.14621098666666e9, V = -3.6535959466666627e9, Ti = 6.316666603088379, Fr = 2.27070703125e11)
 (U = -3.0691708159999967e9, V = 4.304277162666662e9, Ti = 6.316666603088379, Fr = 2.27070703125e11)
 (U = 3.091809799111108e9, V = 6.579570684444437e8, Ti = 6.316666603088379, Fr = 2.27070703125e11)
 (U = 3.091795150222219e9, V = 6.578363697777771e8, Ti = 6.316666603088379, Fr = 2.27070703125e11)
 (U = -2.1445705386666644e9, V = 3.1792032426666636e9, Ti = 6.316666603088379, Fr = 2.27070703125e11)
 (U = 4.0222341617777734e9, V = -4.6493393599999946e8, Ti = 6.316666603088379, Fr = 2.27070703125e11)
 (U = 4.0222294044444404e9, V = -4.650499119999995e8, Ti = 6.316666603088379, Fr = 2.27070703125e11)
 (U = 9.304277457777768e8, V = -1.1228893119999988e9, Ti = 6.316666603088379, Fr = 2.27070703125e11)
 (U = -6.146197247999993e9, V = 3.653717759999996e9, Ti = 6.316666603088379, Fr = 2.27070703125e11)
 (U = 10322.971652560753, V = 118544.06944444432, Ti = 6.316666603088379, Fr = 2.27070703125e11)
 (U = 6.211784632888882e9, V = -3.5477483875555515e9, Ti = 6.583333253860474, Fr = 2.27070703125e11)
 (U = 6.211780195555549e9, V = -3.547863480888885e9, Ti = 6.583333253860474, Fr = 2.27070703125e11)
 (U = 2.068902503111109e9, V = -3.1441626808888855e9, Ti = 6.583333253860474, Fr = 2.27070703125e11)
 (U = 3.034031331555552e9, V = -4.2528075377777734e9, Ti = 6.583333253860474, Fr = 2.27070703125e11)
 (U = 3.177752227555552e9, V = 7.050634364444437e8, Ti = 6.583333253860474, Fr = 2.27070703125e11)
 (U = 3.177752810666663e9, V = 7.049488924444437e8, Ti = 6.583333253860474, Fr = 2.27070703125e11)
 (U = 4.1428866915555515e9, V = -4.035797448888885e8, Ti = 6.583333253860474, Fr = 2.27070703125e11)
 (U = 9.651305439999989e8, V = -1.108645550222221e9, Ti = 6.583333253860474, Fr = 2.27070703125e11)
 (U = 4.1428811519999957e9, V = -4.036979662222218e8, Ti = 6.583333253860474, Fr = 2.27070703125e11)
 (U = 6186.6973334418335, V = 118668.12304687487, Ti = 6.583333253860474, Fr = 2.27070703125e11)
 (U = -266936.3645833331, V = -1.583141166666665e6, Ti = 6.849999904632568, Fr = 2.27070703125e11)
 (U = 6.235954645333326e9, V = -3.4543274026666627e9, Ti = 6.849999904632568, Fr = 2.27070703125e11)
 (U = 2.987818652444441e9, V = -4.207558826666662e9, Ti = 6.849999904632568, Fr = 2.27070703125e11)
 (U = 6.235956423111104e9, V = -3.454211199999996e9, Ti = 6.849999904632568, Fr = 2.27070703125e11)
 (U = 1.9927185599999979e9, V = -3.1136442239999967e9, Ti = 6.849999904632568, Fr = 2.27070703125e11)
 (U = 6.23440708266666e9, V = -3.468812543999996e9, Ti = 6.849999904632568, Fr = 2.27070703125e11)
 (U = -2.9972223146666636e9, V = 4.213749674666662e9, Ti = 6.849999904632568, Fr = 2.27070703125e11)
 (U = 3.2481352177777743e9, V = 7.533482702222215e8, Ti = 6.849999904632568, Fr = 2.27070703125e11)
 (U = 3.2481338097777743e9, V = 7.532299875555547e8, Ti = 6.849999904632568, Fr = 2.27070703125e11)
 (U = -2.006915669333331e9, V = 3.1172445866666636e9, Ti = 6.849999904632568, Fr = 2.27070703125e11)
 (U = 4.243231630222218e9, V = -3.4069073333333296e8, Ti = 6.849999904632568, Fr = 2.27070703125e11)
 (U = 4.243239466666662e9, V = -3.4056780622222185e8, Ti = 6.849999904632568, Fr = 2.27070703125e11)
 (U = 9.951025973333323e8, V = -1.0939171839999988e9, Ti = 6.849999904632568, Fr = 2.27070703125e11)
 (U = -6.23440452266666e9, V = 3.4689294506666627e9, Ti = 6.849999904632568, Fr = 2.27070703125e11)
 (U = 2020.1427137586786, V = 118729.78884548598, Ti = 6.849999904632568, Fr = 2.27070703125e11)
 (U = -343950.03624999966, V = -1.589413219999998e6, Ti = 7.183333396911621, Fr = 2.27070703125e11)
 (U = 6.223212003555549e9, V = -3.337163527111108e9, Ti = 7.183333396911621, Fr = 2.27070703125e11)
 (U = 6.223215416888882e9, V = -3.337282766222219e9, Ti = 7.183333396911621, Fr = 2.27070703125e11)
 (U = 2.909506687999997e9, V = -4.152156949333329e9, Ti = 7.183333396911621, Fr = 2.27070703125e11)
 (U = 1.8838133048888867e9, V = -3.07722722133333e9, Ti = 7.183333396911621, Fr = 2.27070703125e11)
 (U = 6.225276682239994e9, V = -3.340931624959996e9, Ti = 7.183333396911621, Fr = 2.27070703125e11)
 (U = -2.913978071039997e9, V = 4.1530716876799955e9, Ti = 7.183333396911621, Fr = 2.27070703125e11)
 (U = 3.313704135111108e9, V = 8.14994561777777e8, Ti = 7.183333396911621, Fr = 2.27070703125e11)
 (U = 3.3137085368888855e9, V = 8.148773351111102e8, Ti = 7.183333396911621, Fr = 2.27070703125e11)
 (U = -1.889491799039998e9, V = 3.077254901759997e9, Ti = 7.183333396911621, Fr = 2.27070703125e11)
 (U = 4.339403889777773e9, V = -2.5993809599999973e8, Ti = 7.183333396911621, Fr = 2.27070703125e11)
 (U = 4.339408696888885e9, V = -2.6005522088888863e8, Ti = 7.183333396911621, Fr = 2.27070703125e11)
 (U = 1.0256997831111101e9, V = -1.0749325617777767e9, Ti = 7.183333396911621, Fr = 2.27070703125e11)
 (U = -6.225278402559994e9, V = 3.3410462617599964e9, Ti = 7.183333396911621, Fr = 2.27070703125e11)
 (U = -3199.665127224389, V = 118718.70724826376, Ti = 7.183333396911621, Fr = 2.27070703125e11)
 (U = -401459.5972222218, V = -1.5953378194444429e6, Ti = 7.450000047683716, Fr = 2.27070703125e11)
 (U = 6.178709461333326e9, V = -3.243972359111108e9, Ti = 7.450000047683716, Fr = 2.27070703125e11)
 (U = 2.8307861475555525e9, V = -4.1090220657777734e9, Ti = 7.450000047683716, Fr = 2.27070703125e11)
 (U = 6.178720483555549e9, V = -3.244096682666663e9, Ti = 7.450000047683716, Fr = 2.27070703125e11)
 (U = 1.786251743999998e9, V = -3.049647466666663e9, Ti = 7.450000047683716, Fr = 2.27070703125e11)
 (U = 6.179113016888882e9, V = -3.242380138666663e9, Ti = 7.450000047683716, Fr = 2.27070703125e11)
 (U = 3.347918890666663e9, V = 8.650494186666657e8, Ti = 7.450000047683716, Fr = 2.27070703125e11)
 (U = -2.8311915662222195e9, V = 4.107428693333329e9, Ti = 7.450000047683716, Fr = 2.27070703125e11)
 (U = 3.347926620444441e9, V = 8.649315555555546e8, Ti = 7.450000047683716, Fr = 2.27070703125e11)
 (U = -1.7866596017777758e9, V = 3.048053795555552e9, Ti = 7.450000047683716, Fr = 2.27070703125e11)
 (U = 4.392453788444439e9, V = -1.9432878977777755e8, Ti = 7.450000047683716, Fr = 2.27070703125e11)
 (U = 4.392462890666661e9, V = -1.944448759999998e8, Ti = 7.450000047683716, Fr = 2.27070703125e11)
 (U = 1.0445360586666656e9, V = -1.0593764586666656e9, Ti = 7.450000047683716, Fr = 2.27070703125e11)
 (U = -6.179119544888882e9, V = 3.2424975217777743e9, Ti = 7.450000047683716, Fr = 2.27070703125e11)
 (U = -7359.556450737839, V = 118639.36892361098, Ti = 7.450000047683716, Fr = 2.27070703125e11)
 (U = -453840.1640624995, V = -1.6017645034722206e6, Ti = 7.7166666984558105, Fr = 2.27070703125e11)
 (U = 2.7382116835555525e9, V = -4.0671804444444404e9, Ti = 7.7166666984558105, Fr = 2.27070703125e11)
 (U = 6.103944163555549e9, V = -3.151683946666663e9, Ti = 7.7166666984558105, Fr = 2.27070703125e11)
 (U = 1.6799544319999983e9, V = -3.023603527111108e9, Ti = 7.7166666984558105, Fr = 2.27070703125e11)
 (U = 6.103956750222216e9, V = -3.1518036479999967e9, Ti = 7.7166666984558105, Fr = 2.27070703125e11)
 (U = 6.10439825066666e9, V = -3.150082410666663e9, Ti = 7.7166666984558105, Fr = 2.27070703125e11)
 (U = 3.3657344497777743e9, V = 9.154982524444435e8, Ti = 7.7166666984558105, Fr = 2.27070703125e11)
 (U = -2.7386646755555525e9, V = 4.065578346666662e9, Ti = 7.7166666984558105, Fr = 2.27070703125e11)
 (U = 3.365744974222219e9, V = 9.153749564444435e8, Ti = 7.7166666984558105, Fr = 2.27070703125e11)
 (U = -1.6804040142222204e9, V = 3.02200078933333e9, Ti = 7.7166666984558105, Fr = 2.27070703125e11)
 (U = 4.423988807111106e9, V = -1.2808080399999987e8, Ti = 7.7166666984558105, Fr = 2.27070703125e11)
 (U = 4.424001663999995e9, V = -1.2819521311111099e8, Ti = 7.7166666984558105, Fr = 2.27070703125e11)
 (U = 1.0582553617777767e9, V = -1.0435755999999989e9, Ti = 7.7166666984558105, Fr = 2.27070703125e11)
 (U = -6.104410581333326e9, V = 3.150201891555552e9, Ti = 7.7166666984558105, Fr = 2.27070703125e11)
 (U = -11483.574490017349, V = 118497.78168402765, Ti = 7.7166666984558105, Fr = 2.27070703125e11)
 (U = -504000.3645833328, V = -1.6089620833333316e6, Ti = 7.983333349227905, Fr = 2.27070703125e11)
 (U = 5.999277468444439e9, V = -3.060741980444441e9, Ti = 7.983333349227905, Fr = 2.27070703125e11)
 (U = 2.6322108302222195e9, V = -4.0268253795555515e9, Ti = 7.983333349227905, Fr = 2.27070703125e11)
 (U = 5.999288348444438e9, V = -3.060856782222219e9, Ti = 7.983333349227905, Fr = 2.27070703125e11)
 (U = 1.565423818666665e9, V = -2.999217870222219e9, Ti = 7.983333349227905, Fr = 2.27070703125e11)
 (U = 5.999780252444438e9, V = -3.0591319537777743e9, Ti = 7.983333349227905, Fr = 2.27070703125e11)
 (U = -2.6327225742222195e9, V = 4.025219150222218e9, Ti = 7.983333349227905, Fr = 2.27070703125e11)
 (U = 3.3670615039999967e9, V = 9.660856764444433e8, Ti = 7.983333349227905, Fr = 2.27070703125e11)
 (U = 3.367077162666663e9, V = 9.659668408888879e8, Ti = 7.983333349227905, Fr = 2.27070703125e11)
 (U = 4.433867007999995e9, V = -6.164012344444438e7, Ti = 7.983333349227905, Fr = 2.27070703125e11)
 (U = -1.5659301759999983e9, V = 2.99760936533333e9, Ti = 7.983333349227905, Fr = 2.27070703125e11)
 (U = 4.433851406222218e9, V = -6.1523151111111045e7, Ti = 7.983333349227905, Fr = 2.27070703125e11)
 (U = 1.0667900266666656e9, V = -1.0276084337777767e9, Ti = 7.983333349227905, Fr = 2.27070703125e11)
 (U = -5.999801315555549e9, V = 3.059254300444441e9, Ti = 7.983333349227905, Fr = 2.27070703125e11)
 (U = -15551.297851562484, V = 118294.64453124987, Ti = 7.983333349227905, Fr = 2.27070703125e11)
)
   )
ObservedInstrumentModel
  with Jones: SingleStokesGain
  with reference basis: PolarizedTypes.CirBasis()Data Products: Comrade.EHTVisibilityDatum

Optimization and Sampling

Now we need to actually compute our image. For this we will first follow standard approaches in VLBI and find the maximum a posteriori (MAP) estimate of the image and instrument model. This is done using the comrade_opt function which accepts the posterior, an optimization algorithm, and some keyword arguments. For this tutorial we will use the L-BFGS algorithm. For more information about the optimization algorithms available see the Optimization.jl docs.

julia

using Optimization, OptimizationLBFGSB
xopt, sol = comrade_opt(
    post, LBFGSB(); initial_params = prior_sample(rng, post),
    maxiters = 2000, g_tol = 1.0e-1
);

Warning

Fitting gains tends to be very difficult, meaning that optimization can take a lot longer. The upside is that we usually get nicer images.

First we will evaluate our fit by plotting the residuals

julia

using CairoMakie
using DisplayAs
res = residuals(post, xopt)
plotfields(res[1], :uvdist, :res) |> DisplayAs.PNG |> DisplayAs.Text

These look reasonable, although there may be some minor overfitting. This could be improved in a few ways, but that is beyond the goal of this quick tutorial. Plotting the image, we see that we have a much cleaner version of the closure-only image from Imaging a Black Hole using only Closure Quantities.

julia

g = imagepixels(fovx, fovy, 128, 128)
img = intensitymap(skymodel(post, xopt), g)
imageviz(img, size = (500, 400)) |> DisplayAs.PNG |> DisplayAs.Text

Because we also fit the instrument model, we can inspect their parameters. First, let's query the posterior object with the optimal parameters to get the instrument model.

julia

intopt = instrumentmodel(post, xopt)

126-element Comrade.SiteArray{ComplexF64, 1, Vector{ComplexF64}, Vector{Comrade.IntegrationTime{Int64, Float64}}, Vector{Comrade.FrequencyChannel{Float64, Int64}}, Vector{Symbol}}:
   1.0128324263169777 + 0.0im
 -0.15770103207980632 - 1.0034335563546242im
     1.00715221306648 + 0.0im
  -0.5984400223187923 - 0.8363093078848629im
   0.6310817337383087 + 0.4654408791611354im
 -0.36435194162008283 - 1.051459150051846im
   1.0038265716525712 + 0.0im
  -0.6400975626890201 - 0.8070752103148986im
   0.6934459052920989 + 0.6094490627935497im
  -0.5576251235986284 - 0.9043093671495606im
   1.0192903932639523 + 0.0im
  -0.6600798764074701 - 0.7766074888109628im
   0.6094003798814168 + 0.5914629477487924im
  -0.7005196423199624 - 0.7577909687140906im
   0.9970453918373701 + 0.0im
    -0.71275498469719 - 0.7536021351108658im
   0.5689415804823951 + 0.6258665975190575im
  -0.7920541122536097 - 0.5699355526022637im
   1.0361701470240006 + 0.0im
  0.13844214040422087 + 0.5808647399898206im
   0.8528672540524315 - 0.2938205635668648im
   1.0219094613861677 + 0.0im
  -0.7532366411444555 - 0.6864473567733708im
   0.3564774581533805 + 0.394089328170687im
  -0.9796358667276772 - 0.22989620248190443im
   1.0106353524508074 + 0.0im
  -0.7746156088114069 - 0.6615490139087048im
   0.3224407477047667 + 0.3546796328969029im
  -1.0355512340477382 - 0.033990622017430805im
   1.0224863433935587 + 0.0im
  -0.7939589342640039 - 0.6371427328165785im
  0.32362562206324813 + 0.3356799705438963im
  -1.0211584560161377 + 0.1804531353759631im
    1.007751260991969 + 0.0im
  -0.8026910681959919 - 0.6393032988393835im
   0.8986625534921335 + 0.33364217504209054im
  0.30003143122040515 + 0.27825962377609764im
  -0.9763489192539099 + 0.4203463587413953im
   0.9842574545074682 + 0.0im
  -0.8020403581745975 - 0.6805727964891728im
   0.9020897850101944 + 0.18886338574478367im
  0.35197204484969613 + 0.2724632406212819im
  -0.8261168800827081 + 0.6040906103300759im
   1.0419607736892011 + 0.0im
  -0.7767952185448218 - 0.6347797487061446im
   0.9484697364013773 + 0.029467919565401988im
    0.460464194273073 + 0.2656714598996466im
  -0.7161223550432674 + 0.8072844740898627im
   0.9741711199339589 + 0.0im
  -0.7632069596179887 - 0.7390959649201451im
   0.9496369620888456 - 0.13946492260936147im
  0.39561926025826505 - 0.9650504065401782im
  0.42820000307429273 + 0.19682135431485823im
  -0.5186478359993152 + 0.9795113508975535im
    1.030912400618893 + 0.0im
  -0.7561156900131101 - 0.658975894971847im
   0.9107699820121835 - 0.3336125082393986im
   0.5803820984365956 - 0.8405957122945993im
   0.4376373176938829 + 0.12070521817206992im
 -0.25850406358956074 + 1.0706399118338565im
  -0.6089313232544407 - 0.8117921732311348im
   0.9594585970281224 + 0.0im
  -0.7752661259582635 - 0.7323286232971853im
    0.842367360587754 - 0.48887798095079366im
   0.6930560012121669 - 0.7280423654830743im
   0.4211806491143241 + 0.031773971920243796im
 -0.07415580280961577 + 1.2216172743529135im
  -0.7456581394101057 - 0.6877660617135617im
   0.8763622196208114 + 0.0im
  -0.7874819318048909 - 0.8750603544725678im
   0.7535289899309744 - 0.6048989313864533im
   0.8140927381267441 - 0.6132938136242346im
   0.4210937556470165 - 0.06185600601972403im
  0.10436857688594268 + 0.900358624916992im
  -0.8829306300644443 - 0.5262640737049606im
   1.0137260492188291 + 0.0im
   0.6129466449414646 - 0.7633802522536864im
   0.8836670677592738 - 0.49918508660833305im
   0.3729282031711632 - 0.2018858181068971im
  -0.9610297375794468 - 0.3378430094498872im
   1.0653736495467894 + 0.0im
  -0.5958541304525377 - 0.7771900430509742im
   0.3753289941004512 - 0.909399471016624im
   0.9555692987426572 - 0.3769467886626964im
  0.13783754663105105 - 0.4226369974448228im
   1.0149127918833019 + 0.0im
  -0.5699006738884422 - 0.8474742840869786im
  0.18565507424149522 - 0.954585044511965im
   0.9679243531447982 - 0.305043975888851im
 -0.13502581550172182 - 0.35412399483957396im
  -1.0022268645391217 - 0.12112420520997405im
   1.0859294662547727 + 0.0im
 -0.07759227979582019 - 0.9852391673412179im
   0.9889997453782464 - 0.23810127367358505im
 -0.29277472332827953 - 0.282287869790243im
   -1.015577382618973 - 0.010422417040714678im
   1.1140735487832354 + 0.0im
  -0.4768627206409162 - 0.7994851606018073im
 -0.35117340238280415 - 0.8941476678387069im
   0.9967648424176495 - 0.19314410729733442im
  -0.3158969007651894 - 0.14631444442344058im
  -1.0071865318534763 + 0.10604689932435692im
   1.0358294000542552 + 0.0im
  -0.4256500658099036 - 0.907839203595079im
  -0.5648961306955103 - 0.7743587136884738im
   1.0038051381966104 - 0.16132164256475837im
 -0.38314911571082433 - 0.040973493164552434im
  -1.0126739801804154 + 0.1360479562321885im
   1.0098432239195876 + 0.0im
  -0.4004722469751987 - 0.9521415761715858im
   -0.726346141143774 - 0.5959199503815598im
   1.0047034756129742 - 0.14922738958391008im
 -0.49345845020759216 + 0.05264219394923949im
    -1.00827931538545 + 0.16827215176807536im
   0.9830349817354289 + 0.0im
  -0.3286000684762058 - 0.9915119341513305im
  -0.8223614433517974 - 0.40160911538006644im
    1.001918456068092 - 0.15924099843307535im
  -0.4712997437670932 + 0.14850434650371158im
  -1.0029330605780067 + 0.1389719308211633im
   1.0061395756565013 + 0.0im
  -0.2644638051855477 - 0.9966635598345622im
  -0.7997223035171793 - 0.27953499923495917im
   1.0111343552371506 - 0.14433505261767443im
 -0.42865817063832334 + 0.21589465536850588im
  -1.0010739626596934 + 0.1259013550292596im

This returns a SiteArray object which contains the gains as a flat vector with metadata about the sites, time, and frequency. To visualize the gains we can use the plotcaltable function which automatically plots the gains. Since the gains are complex we will first plot the phases.

julia

plotcaltable(angle.(intopt)) |> DisplayAs.PNG |> DisplayAs.Text

Due to the a priori calibration of the data, the gain phases are quite stable and just drift over time. Note that the one outlier around 2.5 UT is when ALMA leaves the array. As a result we shift the reference antenna to APEX on that scan causing the gain phases to jump.

We can also plot the gain amplitudes

julia

plotcaltable(abs.(intopt)) |> DisplayAs.PNG |> DisplayAs.Text

Here we find relatively stable gains for most stations. The exception is LMT which has a large offset after 2.5 UT. This is a known issue with the LMT in 2017 and is due to pointing issues. However, we can see the power of simultaneous imaging and instrument modeling as we are able to solve for the gain amplitudes and get a reasonable image.

One problem with the MAP estimate is that it does not provide uncertainties on the image. That is we are unable to statistically assess which components of the image are certain. Comrade is really a Bayesian imaging and calibration package for VLBI. Therefore, our goal is to sample from the posterior distribution of the image and instrument model. This is a very high-dimensional distribution with typically 1,000 - 100,000 parameters. To sample from this very high dimensional distribution, Comrade has an array of samplers that can be used. However, note that `Comrade` also satisfies the `LogDensityProblems.jl` interface. Therefore, you can use any package that supports `LogDensityProblems.jl` if you have your own fancy sampler.

For this example, we will use HMC, specifically the NUTS algorithm. For However, due to the need to sample a large number of gain parameters, constructing the posterior can take a few minutes. Therefore, for this tutorial, we will only do a quick preliminary run.

julia

using AdvancedHMC
chain = sample(rng, post, NUTS(0.8), 1_000; n_adapts = 500, initial_params = xopt)

PosteriorSamples
  Samples size: (1000,)
  sampler used: AHMC
Mean
┌─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────┬────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────┐
│ sky                                                                                                                                                                                                                                                             │ instrument                                                                                                                                                                                                                                                                                                                                                                                                                                 │
│ @NamedTuple{c::@NamedTuple{params::Matrix{Float64}, hyperparams::Float64}, σimg::Float64, fg::Float64}                                                                                                                                                          │ @NamedTuple{lg::Comrade.SiteArray{Float64, 1, Vector{Float64}, Vector{Comrade.IntegrationTime{Int64, Float64}}, Vector{Comrade.FrequencyChannel{Float64, Int64}}, Vector{Symbol}}, gp::Comrade.SiteArray{Float64, 1, Vector{Float64}, Vector{Comrade.IntegrationTime{Int64, Float64}}, Vector{Comrade.FrequencyChannel{Float64, Int64}}, Vector{Symbol}}}                                                                                  │
├─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────┼────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────┤
│ (c = (params = [-0.0291094 -0.00225421 … 0.0216955 0.0236794; -0.0181094 -0.00155772 … 0.0302226 0.0108692; … ; 0.0209937 0.0423676 … 0.0626501 0.0611539; 0.0181415 0.0146608 … 0.0265551 -0.00889806], hyperparams = 34.3558), σimg = 1.99357, fg = 0.421652) │ (lg = [0.00594672, 0.0276283, 0.00480356, 0.0284994, -0.350063, 0.126732, 0.00192549, 0.0304128, -0.168388, 0.0877074  …  -0.0434778, 0.0143653, -0.837161, 0.0126092, 0.00503483, 0.0301617, -0.109414, 0.0217365, -0.886427, 0.0103596], gp = [0.0, -0.766224, 0.0, -2.1922, 0.947843, -0.931819, 0.0, -2.24157, 1.0096, -1.14359  …  -2.35496, -0.111776, -0.302359, 0.906863, 0.0, -1.82973, -2.48048, -0.0844765, 2.20203, 0.624686]) │
└─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────┴────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────┘
Std. Dev.
┌────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────┬────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────┐
│ sky                                                                                                                                                                                                                                    │ instrument                                                                                                                                                                                                                                                                                                                                                                                                                     │
│ @NamedTuple{c::@NamedTuple{params::Matrix{Float64}, hyperparams::Float64}, σimg::Float64, fg::Float64}                                                                                                                                 │ @NamedTuple{lg::Comrade.SiteArray{Float64, 1, Vector{Float64}, Vector{Comrade.IntegrationTime{Int64, Float64}}, Vector{Comrade.FrequencyChannel{Float64, Int64}}, Vector{Symbol}}, gp::Comrade.SiteArray{Float64, 1, Vector{Float64}, Vector{Comrade.IntegrationTime{Int64, Float64}}, Vector{Comrade.FrequencyChannel{Float64, Int64}}, Vector{Symbol}}}                                                                      │
├────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────┼────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────┤
│ (c = (params = [0.510558 0.560139 … 0.555882 0.540653; 0.594193 0.634351 … 0.65004 0.58107; … ; 0.604785 0.622589 … 0.618572 0.583014; 0.56876 0.581371 … 0.569766 0.557317], hyperparams = 21.6653), σimg = 0.202466, fg = 0.0501662) │ (lg = [0.148407, 0.148837, 0.0198994, 0.0212641, 0.0734756, 0.0649633, 0.0185407, 0.020678, 0.0663128, 0.0629422  …  0.0475308, 0.0178749, 0.0616498, 0.0175936, 0.0168237, 0.0199535, 0.0478732, 0.0190105, 0.0610418, 0.0175704], gp = [0.0, 0.361329, 0.0, 0.0171294, 0.165646, 0.364756, 0.0, 0.0173166, 0.16153, 0.367835  …  0.16562, 0.253086, 3.01074, 2.79097, 0.0, 0.0180674, 0.163787, 0.251966, 1.99781, 2.87411]) │
└────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────┴────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────┘

Note

The above sampler will store the samples in memory, i.e. RAM. For large models this can lead to out-of-memory issues. To avoid this we recommend using the saveto = DiskStore() kwargs which periodically saves the samples to disk limiting memory useage. You can load the chain using load_samples(diskout) where diskout is the object returned from sample.

Now we prune the adaptation phase

julia

chain = chain[501:end]

PosteriorSamples
  Samples size: (500,)
  sampler used: AHMC
Mean
┌───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────┬─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────┐
│ sky                                                                                                                                                                                                                                                           │ instrument                                                                                                                                                                                                                                                                                                                                                                                                                                  │
│ @NamedTuple{c::@NamedTuple{params::Matrix{Float64}, hyperparams::Float64}, σimg::Float64, fg::Float64}                                                                                                                                                        │ @NamedTuple{lg::Comrade.SiteArray{Float64, 1, Vector{Float64}, Vector{Comrade.IntegrationTime{Int64, Float64}}, Vector{Comrade.FrequencyChannel{Float64, Int64}}, Vector{Symbol}}, gp::Comrade.SiteArray{Float64, 1, Vector{Float64}, Vector{Comrade.IntegrationTime{Int64, Float64}}, Vector{Comrade.FrequencyChannel{Float64, Int64}}, Vector{Symbol}}}                                                                                   │
├───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────┼─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────┤
│ (c = (params = [-0.0243213 0.0245655 … -0.0216799 -0.011827; -0.010915 0.00367128 … -0.0345316 -0.0290166; … ; 0.0460715 0.0419414 … 0.0266523 0.0202475; 0.0348695 0.0147173 … -0.0114064 -0.0148543], hyperparams = 34.8259), σimg = 2.0071, fg = 0.435175) │ (lg = [0.0172179, 0.0235857, 0.00494417, 0.0281882, -0.345205, 0.132622, 0.00297411, 0.0294023, -0.162546, 0.0936986  …  -0.0357132, 0.0155571, -0.831933, 0.0124704, 0.00507508, 0.0296937, -0.101636, 0.0232079, -0.88144, 0.00960929], gp = [0.0, -0.746245, 0.0, -2.19201, 1.00095, -0.91044, 0.0, -2.24065, 1.05997, -1.12134  …  -2.30558, -0.0262927, -0.932174, -0.208057, 0.0, -1.83069, -2.43101, 0.00106522, 1.7873, -0.370299]) │
└───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────┴─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────┘
Std. Dev.
┌────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────┬─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────┐
│ sky                                                                                                                                                                                                                                    │ instrument                                                                                                                                                                                                                                                                                                                                                                                                                  │
│ @NamedTuple{c::@NamedTuple{params::Matrix{Float64}, hyperparams::Float64}, σimg::Float64, fg::Float64}                                                                                                                                 │ @NamedTuple{lg::Comrade.SiteArray{Float64, 1, Vector{Float64}, Vector{Comrade.IntegrationTime{Int64, Float64}}, Vector{Comrade.FrequencyChannel{Float64, Int64}}, Vector{Symbol}}, gp::Comrade.SiteArray{Float64, 1, Vector{Float64}, Vector{Comrade.IntegrationTime{Int64, Float64}}, Vector{Comrade.FrequencyChannel{Float64, Int64}}, Vector{Symbol}}}                                                                   │
├────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────┼─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────┤
│ (c = (params = [0.541964 0.570475 … 0.547423 0.564993; 0.592064 0.641471 … 0.648574 0.562291; … ; 0.612713 0.654454 … 0.599796 0.573622; 0.581939 0.589508 … 0.574155 0.576854], hyperparams = 21.9001), σimg = 0.205603, fg = 0.0457) │ (lg = [0.151529, 0.154518, 0.0201872, 0.0227134, 0.0754855, 0.0621, 0.0176351, 0.0205681, 0.0687485, 0.0576778  …  0.0452482, 0.0181772, 0.0628307, 0.0170137, 0.0167872, 0.019709, 0.0469676, 0.019031, 0.0623365, 0.0164684], gp = [0.0, 0.310844, 0.0, 0.0167329, 0.155484, 0.314228, 0.0, 0.0175852, 0.154411, 0.31623  …  0.161631, 0.281152, 2.87415, 2.91047, 0.0, 0.0187262, 0.160577, 0.277016, 2.41482, 2.89561]) │
└────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────┴─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────┘

Warning

This should be run for likely an order of magnitude more steps to properly estimate expectations of the posterior

Now that we have our posterior, we can put error bars on all of our plots above. Let's start by finding the mean and standard deviation of the gain phases

julia

mchain = Comrade.rmap(mean, chain);
schain = Comrade.rmap(std, chain);

Now we can use the measurements package to automatically plot everything with error bars. First we create a caltable the same way but making sure all of our variables have errors attached to them.

julia

using Measurements
gmeas = instrumentmodel(post, (; instrument = map((x, y) -> Measurements.measurement.(x, y), mchain.instrument, schain.instrument)))
ctable_am = caltable(abs.(gmeas))
ctable_ph = caltable(angle.(gmeas))

─────────┬───────────────────┬────────────────────────────────────────────────────────────────────────────────────
      Ti │                Fr │      AA            AP          AZ           JC          LM          PV          SM
─────────┼───────────────────┼────────────────────────────────────────────────────────────────────────────────────
 0.92 hr │ 227.070703125 GHz │ 0.0±0.0       missing     missing      missing     missing  -0.75±0.31     missing
 1.22 hr │ 227.070703125 GHz │ 0.0±0.0  -2.192±0.017     missing      missing    1.0±0.16  -0.91±0.31     missing
 1.52 hr │ 227.070703125 GHz │ 0.0±0.0  -2.241±0.018     missing      missing   1.06±0.15  -1.12±0.32     missing
 1.82 hr │ 227.070703125 GHz │ 0.0±0.0  -2.276±0.017     missing      missing   1.09±0.15   -1.3±0.32     missing
 2.12 hr │ 227.070703125 GHz │ 0.0±0.0   -2.33±0.017     missing      missing   1.15±0.15   -1.5±0.32     missing
 2.45 hr │ 227.070703125 GHz │ missing       0.0±0.0     missing      missing   1.66±0.16    0.7±0.33     missing
 2.75 hr │ 227.070703125 GHz │ 0.0±0.0  -2.402±0.018     missing      missing   1.15±0.16  -1.87±0.32     missing
 3.05 hr │ 227.070703125 GHz │ 0.0±0.0  -2.433±0.019     missing      missing   1.16±0.15  -2.05±0.32     missing
 3.35 hr │ 227.070703125 GHz │ 0.0±0.0  -2.463±0.019     missing      missing   1.14±0.16  -2.25±0.32     missing
 3.68 hr │ 227.070703125 GHz │ 0.0±0.0   -2.47±0.017   0.78±0.17      missing    1.1±0.16  -2.44±0.47     missing
 3.98 hr │ 227.070703125 GHz │ 0.0±0.0   -2.439±0.02   0.65±0.17      missing   1.03±0.16    -2.2±1.6     missing
 4.28 hr │ 227.070703125 GHz │ 0.0±0.0  -2.455±0.021   0.48±0.17      missing    0.9±0.16    -1.6±2.3     missing
 4.58 hr │ 227.070703125 GHz │ 0.0±0.0   -2.373±0.02   0.31±0.17   -1.04±0.23   0.82±0.16   -0.16±2.9     missing
 4.92 hr │ 227.070703125 GHz │ 0.0±0.0   -2.425±0.02    0.1±0.17   -0.88±0.24   0.67±0.16     1.6±2.3  -2.12±0.24
 5.18 hr │ 227.070703125 GHz │ 0.0±0.0  -2.385±0.019  -0.08±0.17   -0.76±0.26   0.47±0.16     2.2±1.6  -2.34±0.26
 5.45 hr │ 227.070703125 GHz │ 0.0±0.0   -2.304±0.02  -0.24±0.17   -0.61±0.26   0.25±0.16   2.34±0.96  -2.56±0.36
 5.72 hr │ 227.070703125 GHz │ 0.0±0.0       missing  -0.46±0.18   -0.49±0.27   -0.1±0.16     missing    -2.0±1.9
 6.05 hr │ 227.070703125 GHz │ 0.0±0.0  -2.224±0.019  -0.76±0.17   -0.35±0.28  -0.87±0.16     missing     missing
 6.32 hr │ 227.070703125 GHz │ 0.0±0.0  -2.163±0.019  -0.96±0.17   -0.26±0.28  -1.55±0.16     missing    -1.2±2.6
 6.58 hr │ 227.070703125 GHz │ 0.0±0.0       missing  -1.24±0.17   -0.18±0.29   -2.0±0.16     missing    -0.7±2.8
 6.85 hr │ 227.070703125 GHz │ 0.0±0.0   -2.108±0.02  -1.54±0.17   -0.11±0.29  -2.34±0.15     missing  -0.015±2.9
 7.18 hr │ 227.070703125 GHz │ 0.0±0.0   -2.01±0.019   -1.8±0.17  -0.059±0.29  -2.68±0.15     missing   0.004±2.9
 7.45 hr │ 227.070703125 GHz │ 0.0±0.0  -1.968±0.019  -2.06±0.16   -0.03±0.28    -2.5±1.4     missing    0.12±2.9
 7.72 hr │ 227.070703125 GHz │ 0.0±0.0   -1.89±0.018  -2.31±0.16  -0.026±0.28   -0.93±2.9     missing   -0.21±2.9
 7.98 hr │ 227.070703125 GHz │ 0.0±0.0  -1.831±0.019  -2.43±0.16   0.001±0.28     1.8±2.4     missing   -0.37±2.9
─────────┴───────────────────┴────────────────────────────────────────────────────────────────────────────────────

Now let's plot the phase curves

julia

plotcaltable(ctable_ph) |> DisplayAs.PNG |> DisplayAs.Text

and now the amplitude curves

julia

plotcaltable(ctable_am) |> DisplayAs.PNG |> DisplayAs.Text

Finally let's construct some representative image reconstructions.

julia

samples = skymodel.(Ref(post), chain[begin:5:end])
imgs = intensitymap.(samples, Ref(g))

mimg = mean(imgs)
simg = std(imgs)
fig = Figure(; resolution = (700, 700));
axs = [Axis(fig[i, j], xreversed = true, aspect = 1) for i in 1:2, j in 1:2]
image!(axs[1, 1], mimg, colormap = :afmhot); axs[1, 1].title = "Mean"
image!(axs[1, 2], simg ./ (max.(mimg, 1.0e-8)), colorrange = (0.0, 2.0), colormap = :afmhot);axs[1, 2].title = "Frac. Uncer."
image!(axs[2, 1], imgs[1], colormap = :afmhot);
image!(axs[2, 2], imgs[end], colormap = :afmhot);
hidedecorations!.(axs)
fig |> DisplayAs.PNG |> DisplayAs.Text

And viola, you have just finished making a preliminary image and instrument model reconstruction. In reality, you should run the sample step for many more MCMC steps to get a reliable estimate for the reconstructed image and instrument model parameters.

This page was generated using Literate.jl.

Stokes I Simultaneous Image and Instrument Modeling ​

Load the Data ​

Optimization and Sampling ​

Stokes I Simultaneous Image and Instrument Modeling

Load the Data

Optimization and Sampling