initial disc update

[IainHammond]

Updated 23/07:

Adding disc.py for the purposes of disc forward modelling. I've kept it tidy, verbose and to the point. It should work with most RT codes. The process is pretty self explanatory:

• user provides a disc model and their conf dictionary, and everything is automatic
• PSFs are pulled from the grid, both for the science and reference (obviously the parameters they'd want to test need to be in the PSF grid)
• user can then run their own post-processing

Caveat regarding RDI: due to the 0.5 mag increment of the PSF grid, it's possible that the reference star will have the same PSF sequence as the science cube. In this case, the code takes a segment from the end of the science sequence to use as the reference, so they are still independent.

Possibly the next thing to add:

• making sure the pixel scale of the RT model matches METIS, but this has been tricky since not all RT codes include the pixel scale in the header (maybe we just support the main ones?).
• let users bring their own on- and off-axis sequences instead of relying on the grid, but I don't trust this flexibility yet
• make a notebook and post-processing example, but given that Exeter is next week I wont have time for this

Bonus: adding noise is now much faster, and psf_template had a clean up

[GillesOrban]

Thanks for this PR!
I had a first look and have some remarks and questions:

• In create_disc_sequence, I think there is one mistake on the default for transmission (line 137-138): this is supposed to be the off-axis transmission with the vortex coronagraph (1st column with radial separation, 2nd column with off-axis transmission at the given radial separation) and not f_vc_trans (which is the transmission of the bulk material and is taken into account by background). You can check cc_adi.py and search for OAT.

• It might be appropriate to document (in the docstring of the method) how the relative photometry between the disc and the on-axis star is handled. Actually, if I understand that correctly, this is handled upstream, with your radiative transfer tool which generate a disc image where the central pixel is the star (-> starval). But then the magnitude of the star should match what is used in your RT tool. Adding band, dit, mag, add_bckg, as explicit parameters (rather than just relying on conf) would make that more explicit ?

• I think having two small notebooks illustrating 1/ the forward modelling, 2/ the post-processing, would be strong assets for people who get started. The forward modelling notebook could go to 'notebooks/sources/'. For post-processing, one could put it there: notebooks/post_processing.
  Related to that, I am wondering if it is actually useful to have postproc_disc, which is a pure VIP helper function ? I tend to think that it could live entirely in the post_processing notebook and be removed from disc.py

[IainHammond]

Thanks Gilles, I will take a look into your first two points tomorrow. I re-named a lot of my parameters when converting this so it's possible I got the transmission wrong. I agree about the the post-processing just going into a demo which I will get started on soon

[IainHammond]

I might change the RDI section to be its own function

[IainHammond]

In create_disc_sequence, I think there is one mistake on the default for transmission (line 137-138): this is supposed to be the off-axis transmission with the vortex coronagraph (1st column with radial separation, 2nd column with off-axis transmission at the given radial separation) and not f_vc_trans (which is the transmission of the bulk material and is taken into account by background). You can check cc_adi.py and search for OAT.

@GillesOrban
I'm just trying to address this now, f_oat refers to 'optics/vc/oat_L_RAVC.fits' which is the RAVC and not the CVC, or are they identical here?

[GillesOrban]

f_oat refers to the RAVC because the default config simulate the RAVC. In the datapack (downloaded from the gdrive), there is a oat_L_CVC.fits file.

[oabsil]

Hi @IainHammond, I've finally found the time to take a look at your pull request. I've dropped four comments on the disc.py function. In addition to these detailed comments, can I also please ask you to provide the example notebook requested by @GillesOrban to illustrate the use of disc.py? This would be most useful to wrap up this pull request and start using your contribution.

[oabsil]

Not sure why the off-axis PSF is used for cropping here, as the disc will be injected into the on-axis PSF cube. More generally, I'm not sure what's the use of the off-axis PSF in RDI except for flux normalisation, which probably doesn't need any cropping to be done?

[oabsil]

I must admit that this option is a bit mysterious to me. Can you clarify what the use case is? Would that pertain to the location of an injected protoplanet, for which you would want to explore the effect of extinction without having to produce a new disc image each time you start a simulation?

[oabsil]

This assumes that psf_ON is normalised in a way where the stellar flux is equal to one in the off-axis PSF (total intensity in the off-axis PSF image = 1). This is however not the case: in the PSF data base, the on- and off-axis PSFs are normalised so that the total intensity of the non-coronagraphic PSF is equal to one. The difference between the non-coronagraphic PSF and the off-axis PSF in the coronagraphic mode comes from the throughtput of the Lyot stop and from the intrinsic transmission of the vortex phase mask. The disc cube needs to be multiplied by the throughput of the coronagraphic mode relative to the non-coronagraphic mode to have the same normalisation as psf_ON. This can be done through something like: psf_ON += (cube * np.sum(psf_OFF))

[oabsil]

@GillesOrban just warned me that the scaling with the off-axis PSF flux may already be taken care of through the convolution with pst_OFF in the call to cube_inject_fakedisk with the option normalize_psf=False. It seems indeed to be the case, as far as I can tell from a quick look at that function. It would still be useful to add a note in the code to keep that in mind.

[oabsil]

I think this requires a bit more explanation. From what I understand this is the off-axis stellar photon flux per sec at the detector in an aperture of diameter equal to the FWHM of the (off-axis) PSF.