Vortex Coronagraph

[zhangderby]

I've been after a GPU-compatible vortex coronagraph for a while and know that dygdug's is currently in progress/prysm is GPU compatible. What is the status of dygdug's vortex coronagraph?

[brandondube]

(there is a tl;dr at the bottom)

As a whole, dygdug development is a little stalled out at the moment. It will be resumed eventually, but my brain lacks the spare energy to work on it (or prysm, or...) for a while now.

All of the coronagraph models in dygdug have (maybe?) a scaling issue, in scaling the radiometry properly to units of normalized intensity (let alone contrast). Or I misunderstand something somewhere, which is also possible.

The Vortex model is usable, but has an error somewhere, here is some 7 month old code that works, and produces a DH of ~1e-8 residual. The intensity is asymmetric, so there is an asymmetry error somewhere (incomplete destructive interference, which is how a vortex coro works).

%load_ext autoreload
%autoreload 2

from pathlib import Path

import numpy as np

from dygdug.models import (
    Pupil,
    FPM,
    ImgSamplingSpec,
    SingleDMCoronagraph,
    TwoDMCoronagraph,
    SingleDMVortexCoronagraph,
    curry_bb,
    curry_nb,
    incoh_sum,
    one_sided_annulus,
    plottable
)
from dygdug.efc import (
    EFC,
    calculate_jacobian_fwddiff,
)

from prysm.propagation import Wavefront as WF
from prysm import fttools, coordinates, geometry
from prysm.x.dm import DM
from prysm._richdata import RichData
from prysm.conf import config

from astropy.io import fits

from matplotlib import pyplot as plt

efl = 500
D = 30
fno = efl/D
wvl0 = 0.550
ls_size = 0.95
lamD = fno*wvl0
Nmodel = 512
Npup = 512

iss = ImgSamplingSpec.from_N_lamD_px_per_lamD(N=256, lamD=lamD, px_per_lamD=4)


pu = Pupil.circle(Nmodel, D, Npup)
pu.data = pu.data.astype(config.precision)
aperture = pu.data
aperture2 = aperture * (iss.lamD/iss.dx *np.sqrt(aperture.size)/aperture.sum())
pu.data = aperture2


fpm = FPM.vortex(Nmodel, pu.dx, efl, 4)
ls = Pupil.circle(Nmodel, D*ls_size, Npup*ls_size)

fpm = FPM.vortex(Nmodel, pu.dx, efl, 2, padsamples=0)
c = SingleDMVortexCoronagraph(pu, dm1, fpm, ls, iss, efl)


pak=c.fwd(.550,debug=True)
pak['lyot'][1].intensity.plot2d(xlim=D/2, log=True, clim=(1e-12,1e-6), interpolation='nearest')
pak['img'].intensity.slices().plot('x', yscale='log')

My guess is that I tried to do things with "FFT style" martix DFTs, where the "center" pixel is DC and for even sized arrays there is an asymmetric support in the Fourier domain, and by shifting the matrix DFTs by 1/2 "pixel" this can be remedied. But the outermost scale (highest frequencies) are done by an FFT still, which will have this issue for even sized arrays.

tl;dr there is a functional VC model, which produces ~1e-8 peak residual contrast, which means its errors are at the ~1e-8 level. You could use dygdug's EFC to "conjugate" those residuals with the DM, but if you tried to implement the same DM settings in the lab, they would be wrong.

It is possible that the model is really much better than I think it is, and it is just a matter of broken scale factors. I deeply regret the choices I made in the most recent changes to matrix DFT normalization. They make getting properly parseval's theorem satisfying integrals such a royal pain in the ass. I am probably going to revert that choice in the next version of prysm, so we can just divide ap = ap / ap.size like civilized people to have properly scaled radiometry.