PHARO data kernel-phase analysis example

Table of Contents

1 Introduction

This page features all of the code that was used for the analysis of the the P3K/PHARO dataset featured in the Martinache et al (2020) paper. This page should therefore be a good introduction into how to process your own datasets using XARA and the KERNEL project software tools in general.

2 PHARO medium-cross pupil mask model

Kernel-phase data analysis requires a decent model of the aperture. This dataset was acquired using a pupil mask called the medium-cross Lyot-stop so as to limit some of the flux for this fairly bright star α-Ophiucus. Pupil specifications for this mask are outlined in the Hayward et al (2001) paper describing the PHARO camera.

The script uses XAOSIM (another software product of the KERNEL project) to produce an image of the pupil of PHARO. To avoid rounding errors as much as possible in the creation of this pupil image that should be rotated 2 degrees to actually match the instrument pupil, the discrete representation of the aperture was computed with the spider veins aligned with the aligned with the main image axes.

After all (x,y) coordinates are computed, it is easy to get the right model by multiplying them by a rotation matrix. These coordinates, along with the unchanged pupil transmission model are what will serve as the basis to produce the kernel-phase model.

import numpy as np
import matplotlib.pyplot as plt
import as cm
import as pf
import xaosim as xs
from xaosim.pupil import PHARO
from scipy.ndimage import rotate
import xara

binary_model = False # use a finer transmission model for the aperture
PSZ     = 4978*2   # size of the array for the model
pdiam   = 4.978    # telescope diameter in meters
mstep   = 0.160    # step size in meters

pmask   = PHARO(PSZ, PSZ/2, mask="med")
pmask2  = PHARO(PSZ, PSZ/2, mask="med", ang=-2) # rotated!
ppscale = pdiam / PSZ

if binary_model:
        p3k_model = xara.core.create_discrete_model(
                pmask, ppscale, mstep, binary=True, tmin=0.4)
        p3k_model = xara.core.create_discrete_model(
                pmask, ppscale, mstep, binary=False, tmin=0.05)
        p3k_model[:,2] = np.round(p3k_model[:,2],2)

# rotate the model by two degrees
# --------------------------------

th0 = -2.0 * np.pi / 180.0 # rotation angle
rmat = np.array([[np.cos(th0), -np.sin(th0)],
                                 [np.sin(th0), np.cos(th0)]])

p3k_model[:,:2] = p3k_model[:,:2].dot(rmat)

# -------------------------
#      simple plot
# -------------------------
f0 = plt.figure(0)
ax = f0.add_subplot(111)
                 PSZ/2+p3k_model[:,1]/ppscale, 'b.')
f0.set_size_inches(5,5, forward=True)

# -------------------------


This image shows a high-fidelity (10,000 x 10,000 pixels!) image of the pupil rotated and the points overlaid show the location of the center of the virtual sub-apertures that will be used to produce the kernel-phase model of the PHARO data.

These (x,y) coordinates and the transmission function that is associated to each location are then used as input to compute the intial structure that will be used to do the kernel-phase analysis. The result is then saved as a multi-extension fits file such as described in the dedicated documentation page.

# compute the kernel-phase data structure
kpo_0 = xara.KPO(array=p3k_model, bmax=4.646)

# show side by side, the pupil model and its associated uv-coverage
kpo_0.kpi.plot_pupil_and_uv(xymax=2.5, cmap=cm.plasma_r, ssize=9, figsize=(10,5), marker='o')

# and save to a multi-extension kernel-phase fits file for later use
fname = "p3k_med_%s_model.fits" % (mtype,)
print("saving %s" % (fname,))

Attempting to build KPI from array None
1362 distinct baselines were identified
1312 baselines were preserved after filtering
10 first singular values for this array:
[ 185.70474  185.69751  130.22834  130.22309   86.37561   86.3731
   67.30037   67.29565   44.92332   44.92221]

        Summary of properties for no_name
        - 528 sub-apertures
        - 1312 distinct baselines
        - 1048 Ker-phases (79.9 % target phase information recovery)
        - 264 Eig-phases (50.0 % wavefront information recovery)
KPI data successfully created
File provided is not a fits file
saving p3k_med_grey_model.fits


Back to the documentation index

3 PHARO data-set kernel-phase analysis

3.1 Kernel-phase extraction

In order to be able to do the extraction of the Fourier-phase information, the software needs two additional pieces of information: what is the central wavelength of the filter used and what is the plate-scale of the camera. These two pieces of information make it possible to translate the model baselines, thus far computed in meters into radians, and to determine the parameters of the Fourier transform that is required to read the phase in the right place!

The script pre-loads two data-cubes containing 100 images of α-Oph and ε-Her (follow the links to download the data). As explained in the paper, these images were trimmed to 64x64 pixels.

import numpy as np
import matplotlib.pyplot as plt
import as pf
import xara

tgt_cube = pf.getdata('tgt_cube.fits') # alpha Ophiuchi
ca2_cube = pf.getdata('ca2_cube.fits') # epsilon Herculis

pscale = 25.0              # plate scale of the image in mas/pixels
wl     = 2.145e-6          # central wavelength in meters (Hayward paper)
ISZ    = tgt_cube.shape[1] # image size
kpo1   = xara.KPO(fname="p3k_med_grey_model.fits")
kpo2   = kpo1.copy()

        tgt_cube, pscale, wl,target="alpha Ophiuchi", recenter=True)
        ca2_cube, pscale, wl, target="epsilon Herculis", recenter=True)

3.2 Calibration of the kernel-phase

At this point, we have two complete kpo data structures, containing multiple kernel-phase measurements for both sources. Here, we'll simply average these multiple realizations of the kernel-phase vector and store the standard deviation of these realizations. The calibrated kernel-phase vector is the subtraction of the median kernel-phase obtained on the target of interest and that obtained on the calibration star. Some amount of error is added in quadrature to the experimental standard deviation so as to lead to a reduced χ2 = 1 at the later stages of this analysis.

data1 = np.array(kpo1.KPDT)[0]
data2 = np.array(kpo2.KPDT)[0]

mydata = np.median(data1, axis=0) - np.median(data2, axis=0)
myerr  = np.sqrt(np.var(data1, axis=0) / (kpo1.KPDT[0].shape[0] - 1) + np.var(data2, axis=0) / (kpo2.KPDT[0].shape[0] - 1))
myerr = np.sqrt(myerr**2 + 1.2**2)

3.3 Colinearity maps

In the high-contrast approximation, the kernel-phase signature of a binary companion is directly proportional to the contrast (using the secondary/primary convention). It is then somewhat covenient to look in the right ascension - declination parameter space surrounding the central object where the calibrated kernel-phase signal best mimics the signature of a companion. This is what the computation of a colinearity map achieves: for a finite number of laid out on a regular grid, one computes the theoretical signal that a contrasted companion would have, and then look at the result of the dot product between this prediction and the experimental kernel-phase. If a companion is present, the colinearity map will feature a maximum.

Note that due to the properties of kernel-phase that is sensitive to the presence of asymmetric features, the kernel-phase colinearity maps are always anti-symmetric.

print("\ncomputing colinearity map...")
gsize = 100 # gsize x gsize grid
gstep = 10 # grid step in mas
xx, yy = np.meshgrid(
        np.arange(gsize) - gsize/2, np.arange(gsize) - gsize/2)
azim = -np.arctan2(xx, yy) * 180.0 / np.pi
dist = np.hypot(xx, yy) * gstep

#mmap = kpo1.kpd_binary_match_map(100, 10, mydata/myerr, norm=True)
mmap = kpo1.kpd_binary_match_map(100, 10, mydata, norm=True)
x0, y0 = np.argmax(mmap) % gsize, np.argmax(mmap) // gsize
print("max colinearity found for sep = %.2f mas and ang = %.2f deg" % (
        dist[y0, x0], azim[y0, x0]))

f1 = plt.figure(figsize=(5,5))
ax1 = f1.add_subplot(111)
ax1.imshow(mmap, extent=(
        gsize/2*gstep, -gsize/2*gstep, -gsize/2*gstep, gsize/2*gstep))
ax1.set_xlabel("right ascension (mas)")
ax1.set_ylabel("declination (mas)")
ax1.plot([0,0], [0,0], "w*", ms=16)
ax1.set_title("Calibrated signal colinearity map")


Colinearity map for the calibrated kernel-phase. A clear and unambiguous maximum is visible to the left of the primary, slightly less than 100 mas away from the center.

3.4 χ2 minimization

The solution identified at the previous stage makes for a very good first guess to be used by a traditional χ2 minimization procedure.

print("\nbinary model fitting...")
p0 = [dist[y0, x0], azim[y0, x0], mmap.max()] # good starting point

mfit = kpo1.binary_model_fit(p0, calib=kpo2)
p1 = mfit[0] # the best fit parameter vector (sep, P.A., contrast)

cvis_b = xara.core.cvis_binary(
        kpo1.kpi.UVC[:,0], kpo1.kpi.UVC[:,1], wl, p1) # binary
ker_theo =

fig = plt.figure(figsize=(6,6))
ax = fig.add_subplot(111)

ax.errorbar(ker_theo, mydata, yerr=myerr, fmt="none", ecolor='c')
ax.plot(ker_theo, mydata, 'b.')
mmax = np.round(np.abs(mydata).max())
ax.plot([-mmax,mmax],[-mmax,mmax], 'r')
ax.set_ylabel("data kernel-phase")
ax.set_xlabel("model kernel-phase")
ax.set_title('kernel-phase correlation diagram')
ax.axis([-11, 11, -11, 11])

if myerr is not None:
        chi2 = np.sum(((mydata - ker_theo)/myerr)**2) / kpo1.kpi.nbkp
        chi2 = np.sum(((mydata - ker_theo))**2) / kpo1.kpi.nbkp

print("sep = %3f, ang=%3f, con=%3f => chi2 = %.3f" % (p1[0], p1[1], p1[2], chi2))
print("correlation matrix of parameters")
print(np.round(mfit[1], 2))


4 Conclusion

This concludes the basic kernel-phase data analysis tutorial. More refined analysis would require the computation of covariance matrices which this tutorial did not bother with.

Back to the documentation index

Author: Frantz

Created: 2020-05-04 Mon 20:32