GitHub

Extracting Traditional Photometry and Astrometry

Though not its primary purpose, you can use Octofitter to extract traditional astrometry and photometry from one or more images. This uses the functionality in the Fit Orbits to Images tutorial, but with a much simpler model.

Instead of fitting an entire orbit, we will simply fit an X / Y position and brightness.

Start by loading your images:

using Octofitter
using OctofitterImages
using Distributions
using Pigeons
using AstroImages
using CairoMakie

# Load individual iamges
# image1 = load("image1.fits")
# image2 = load("image2.fits")

# Or slices from a cube:
# cube = load("cube1.fits")
# image1 = cube[:,:,1]

# Download sample images from GitHub
download(
    "https://zenodo.org/records/6823071/files/HR8799.2021.fits?download=1",
    "HR8799-2021.fits"
)

# Or multi-extension FITS (this example)
image = AstroImages.load("HR8799-2021.fits")
Example block output

You can preview the image using imview from AstroImages:

imview(image)
Example block output

Note that to accurately extract astrometry and photometry, the input image should have already been convolved with the star or planet point spread function. If this isn't available, a convolution by a Gaussian or Airy disk might be an acceptable approximation.

Build the model

First, we create a table of our image data that will be attached to the Planet:

imglike = ImageLikelihood(
    (
        band=:L,
        image=AstroImages.recenter(image), platescale=9.971,
        epoch=mjd("2021")
    ),
)
OctofitterImages.ImageLikelihood Table with 6 columns and 1 row:
     band  image                 platescale  epoch    contrast              ⋯
   ┌─────────────────────────────────────────────────────────────────────────
 1 │ L     [NaN NaN NaN NaN Na…  9.971       59215.0  218-element extrapo…  ⋯

Note that you can also supply a contrast curve or map directly. If not provided, a simple contrast curve will be calculated directly from the data.

Next create the simplest possible model of 2D position, plus a contrast variable matching the band name used in the ImageLikelihood above:

@planet b Visual{Octofitter.FixedPosition} begin
    sep ~ Uniform(0, 2000)
    pa ~ Uniform(0,2pi)
    # Contrast ratio
    L ~ Uniform(0, 1)
end imglike

@system sys begin
    plx = 24.4620
end b

model = Octofitter.LogDensityModel(sys, verbosity=4)
LogDensityModel for System sys of dimension 3 and 1 epochs with fields .ℓπcallback and .∇ℓπcallback

Sample from the model (locally)

If you already know where the planet is and you only want to extract astrometry from that known location, you can specify a starting point and use hamiltonian monte carlo as follows. This will be very very fast.

initialize!(model, (;
    sep=1704,
    pa=deg2rad(70.63),
    L=1e-4,
))
chain = octofit(model, iterations=10000)

Sample from the model (globally)

You could also try sampling across the entire image, without necessarily specifying a starting position. Note that if there are multiple candidates, taking the naive mean and standard deviation will average across all planets.

using Pigeons
initialize!(model)
chain, pt = octofit_pigeons(model, n_rounds=11)
(chain = MCMC chain (2048×8×1 Array{Float64, 3}), pt = PT(checkpoint = false, ...))

Access results

samples_sep = chain[:b_sep]
samples_pa = chain[:b_pa]
println("The median separation is ", median(samples_sep))
The median separation is 1703.7058646958212

Visualize

using CairoMakie, PairPlots
octocorner(model,chain)
Example block output