MICADO limited magnitude comparison with Davies+2010

[hugobuddel]

Based on slack discussion with Francesco Valentino

I’ve started looking into how PSFs are currently implemented in MICADO simulations in ScopeSim.
The conversation stemmed from my attempt to reproduce something similar to _get_limiting_mags in SIMCADO and compute limiting magnitude for SCAO, MCAO simulations in ScopeSim. Since I get some different results from Davies+2010, we were wondering if this is due to the background, PSF, or other factors.

When I implement some simple Micado simulations with micado = sim.OpticalTrain(micado_cmds), I can see that the code calls for MORFEO.yaml (MCAO) and MICADO_Standalone_RO.yaml (SCAO). There I can see two non-varying PSFs that are used:

-  name : relay_psf
  description : SCAO PSF
  class : FieldConstantPSF
  kwargs:
    filename : PSF_SCAO_ConstPSF_0_5off.fits
    warning : “Default PSF is NOT field varying. See documentation.”

and

-  name : morfeo_generic_psf
  description : MORFEO field varying MCAO PSF
  class : FieldConstantPSF
  kwargs:
    filename : PSF_MCAO_ConstPSF_40_18_6.fits
    warning : “Default PSF is not Field Varying. See Documentation”So here’s some questions:

• How were these PSFs generated? With Anisocado? In the header of the psf fits files, I can read some of the parameters (e.g., for PSF_MCAO_ConstPSF_40_18_6.fits: wavelength=1.2 (um), EsoMedian (atmospheric conditions from Anisocado?), Seeing=0.67, Strehl = 0.047, etc.)
• How’s the wavelength dependence accounted for? It’s an analytical rescaling?
• Looking at Anisocado, one can create a different PSF and pass it on in fits format, I think. Does one have to add the information to a yaml file or can the psf be passed on directly within python, after creating the PSF?
• Is .fits format from Anisocado the standard that will be kept? I ask because other groups are working on realistic PSFs that depend on the availability of nearby stars, distance from the target, etc. We were wondering if / how these PSFs can be implemented in future simulations.

For reference, here’s the plot I was trying to recreate from Davies+10

And here’s my version

Davies+10 numbers (AB):

[hugobuddel]

I don't actually know the answers to these questions directly, @teutoburg or @astronomyk can you chime in?

Nevertheless, my attempt:

• I do believe the SCAO PSF is generated with AnisoCADO, but not sure. The MCAO PSF is not I think.
• As for the wavelength dependence, there are multiple FITS extensions for different wavelenghts, but not sure how it is interpolated. Could be just nearest neighbour
• I think it should be possible to just use the AnisoCADO PSFs directly, as in, the filename of the FieldConstantPSF class can be overwritten
• We don't have any plans for AnisoCADO or the file format at the moment. You are suggesting to stabilize it so other tools can use it too?

As for the limiting magnitude, perhaps we should first try to verify the results that Francesco got.

[hugobuddel]

Assigning to @teutoburg for now, because I have to prepare the METIS FDR. @teutoburg can you check whether we get the same numbers as Francesco? And perhaps get some answers out of @astronomyk.

[teutoburg]

Answers from @astronomyk :

• I do believe the SCAO PSF is generated with AnisoCADO, but not sure. The MCAO PSF is not I think.

By default, the SCAO PSF is from a static file, which was generated using AnisoCADO. The MCAO PSF is also generated fom AnisoCADO using static strehl numbers (K,H,J=40,15,7%). By default, both SCAO & MCAO PSFs are read from static FITS files.

• As for the wavelength dependence, there are multiple FITS extensions for different wavelenghts, but not sure how it is interpolated. Could be just nearest neighbour

"Could be just nearest neighbour" - correct

The PSF FITS extensions are chosen as the midpoint for each broadband filter.

• I think it should be possible to just use the AnisoCADO PSFs directly, as in, the filename of the FieldConstantPSF class can be overwritten

In principle, this is possible and in theory straightforward (quote by @astronomyk ). See the docs.

• We don't have any plans for AnisoCADO or the file format at the moment. You are suggesting to stabilize it so other tools can use it too?

For the moment we're using FITS files, but we're open to creating an effect which uses numpy arrays as an input, for example...

[carmeloarci]

Hi, for the "limiting magnitude comparison" it is essential to cross-check that Davies+2010 and ScopeSim use the same input, the essential ingredient being:

• Sky Background
• Transmissivity/quantum efficiency of the system telescope-optical relay (MICADO stand-alone or Morfeo)-MICADO-Detector
• Detector Dark/RON
• PSF

As reminder. The latest MCAO PSF are here: https://drive.google.com/drive/folders/1-WxTTyrQ5oZVxO420UuNkYRsjeuQe5D0?usp=drive_link

We need to define/check the Signal to Noise ratio definition

• area of the core of the PSF used for the computation of the signal
• the way the noise is estimated

For SIMCADO the function _get_limiting_mag is described here:
https://simcado.readthedocs.io/en/latest/reference/simcado.simulation.html
and this computes the signal from the measurements and the noise using basic principles:

The signal is computed on the simulated image of a grid of equal magnitude stars, by computing the total of the counts found within an aperture (i.e. 4x4 or 2x2 pixel centered on the star peak).

From an observational point of view, it is more practical to compute both signal and noise for the detection of a point source from the image (background subtracted image). In the case (as in the grid of stars) you have multiple identical sources, you may fix an aperture size and on this aperture compute the average of the total of the counts (this is the signal), while the noise is the sdev of the total of the counts you measure on multiple aperture realizations placed on an area of detector free of light sources.

[hugobuddel]

Attached a notebook by Francesco to produce the figures. Based on https://nbviewer.org/github/astronomyk/SimCADO/blob/master/docs/source/_static/python_notebooks/Limiting_Magnitudes.ipynb

test_performance.zip

Produced final figure:

In code:

# # Test on MICADO performances
# (F. Valentino, ESO)
#
# See the example in SimCado: https://nbviewer.org/github/astronomyk/SimCADO/blob/master/docs/source/_static/python_notebooks/Limiting_Magnitudes.ipynb

# +
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.colors import LogNorm
import matplotlib.ticker as mtick

# Update plot parameters
plt.rcParams.update({'font.size': 14, 'axes.labelsize': 14, 'axes.linewidth': 2, 
                     'xtick.minor.width': 2, 'xtick.major.width': 2, 'ytick.major.width': 2, 
                     'ytick.minor.width': 2})


from astropy import units as u
from sklearn.linear_model import LinearRegression
from photutils.aperture import CircularAperture, CircularAnnulus, aperture_photometry
from scipy import interpolate
from scipy.optimize import curve_fit


# ScopeSim
import scopesim as sim
import scopesim_templates as sim_tp
import os

from copy import deepcopy


# This is needed to download the ELT packages - only once in the selected folder
# (this line can be commented out after the first usage):
if os.path.isdir("inst_pkgs") == False:
    sim.download_package(["instruments/MICADO", "instruments/MORFEO", "telescopes/ELT", "locations/Armazones"])
# -

import numpy as np
def get_snr_aperture_photometry(image, x_pos, y_pos, circle_radius, annulus_inner_radius, annulus_outer_radius):
    
    """
    Perform aperture photometry to calculate signal, noise, background, and SNR values for given positions.
    
    Parameters:
        image (numpy.ndarray): The input image.
        x_pos (array-like): X-coordinates of aperture positions.
        y_pos (array-like): Y-coordinates of aperture positions.
        circle_radius (array-like): Radii of circular apertures.
        annulus_inner_radius (array-like): Inner radii of annuli for background estimation.
        annulus_outer_radius (array-like): Outer radii of annuli for background estimation.
        
    Returns:
        tuple: A tuple containing lists of calculated values (flux, err, std, snr_values, sky).
    """
    
    # Initialize lists to store calculated values
    flux = []
    err = []
    std = []
    sky = []
    snr_values = []
    
    i = 0  # Initialize a counter for iterating through positions
    
    # Loop over x_pos and y_pos, which are arrays of x and y coordinates
    for x_cen, y_cen in zip(x_pos.astype(int), y_pos.astype(int)):
        
        y, x = np.ogrid[0:image.shape[0], 0:image.shape[1]]
        
        mask_circle = ((x - x_cen) ** 2 + (y - y_cen) ** 2 <= circle_radius[i] ** 2)
        mask_annulus = ((x - x_cen) ** 2 + (y - y_cen) ** 2 <= annulus_outer_radius ** 2) & \
                       ((x - x_cen) ** 2 + (y - y_cen) ** 2 > annulus_inner_radius ** 2)
        
        background = np.median(image[mask_annulus])
        noise = np.std(image[mask_annulus]) * np.sqrt(len(image[mask_circle]))
        signal = np.sum(image[mask_circle] - background)
        
        flux.append(signal)
        err.append(noise)
        std.append(np.std(image[mask_annulus]))
        sky.append(background)
        snr_values.append([signal / noise])
        
        i = i + 1
    
    # Return the lists of calculated values
    return flux, err, std, snr_values, sky


# Define the function to fit (square root function)
def sqrt_function(x, a):
    return a / np.sqrt(x)


# # Grid of stars
#
# We simulate a grid of point sources with decreasing magnitudes. \
# Pay attention to the system in use to start with (these will be A0V stars, thus Vega!)

# +
n = 100  # number of stars to simulate
mmin = 25  # minimum magnitude [mag]
mmax = 33  # maximum magnitude [mag]

# This is in VEGA!
mags = np.linspace(mmin, mmax, n) * u.mag  # array for mags / snr comparison


circle_radius = np.ones(n) * 5 # aperture for signal computation [pix]
annulus_inner_radius = 35 # inner aperture for noise computation [pix]
annulus_outer_radius = 40 # outer aperture for noise computation [pix]


# Future test for A. Zanella: diameter=28 mas = 7 pixel in 4mas, MCAO
# -

# # Simulation
#
# Ideally, we'll simulate 3 filters over a range of integration times. The resulting S/N should increase as the sqrt(time).
#
# For simplicity, we switch off the linearity of the response and deal with single long integrations only - but this can be fixed. \
# We rely on the code that M. Ginolfi wrote with K. Leschinski to estimate the aperture photometry.

# +
# Let's set a grid of observation times, expressed as ``NDIT`` x ``DIT``
nsteps = 50

# E.g., Log-space steps from 10 mins to 5 hours
dit = np.logspace(np.log10(600.), np.log10(3600.*5.), num=nsteps, endpoint=True)

# Or linear steps:
dit = np.linspace(600., 3600.*5., num=nsteps, endpoint=True)
ndit = np.ones(nsteps, dtype=np.int8)

# Filters:
target_filter = ["J", "H", "Ks"]
filter_wheel_1 = ["J", "open", "open"] 
filter_wheel_2 = ["open", "H", "Ks"] 

# Vega to AB conversion (Vega+conv = AB)
vega_to_AB = np.array([0.9, 1.4, 1.85]) 
refmag = 'ABmag'

# # You want to see the results in Vega?
# vega_to_AB = np.zeros(len(target_filter)) 
# refmag = 'Vega'

# Mode name
mode_name = ["SCAO", "IMG_1.5mas"]

# How many sigmas?
sigma_lim = 5.

# Image size for MICADO
npixx, npixy = 1024, 1024

# +
# We initialize the array of N-sigma limiting magnitudes here
lim_mag = np.ndarray((len(target_filter), nsteps))

for filt_ind, filt in enumerate(target_filter):
    
    # Set the basic commands for the simulation (image size, filter, etc.)
    micado_cmds = sim.UserCommands(use_instrument="MICADO",
                               set_modes=mode_name,
                               properties={"!OBS.filter_name": target_filter[filt_ind],
                                           "!OBS.filter_name_fw1": filter_wheel_1[filt_ind],
                                           "!OBS.filter_name_fw2": filter_wheel_2[filt_ind],
                                           "!DET.width": npixx,
                                           "!DET.height": npixy})

    micado = sim.OpticalTrain(micado_cmds)

    # This is the pixel scale
    pixel_scale = micado.cmds["!INST.pixel_scale"]
    
    # Shut down some of the effects, particularly the linearity of the response
    micado["micado_adc_3D_shift"].include = False
    micado["detector_linearity"].include = False
    micado["full_detector_array"].include = False
    micado["detector_window"].include = True

    
    # Image the stellar grid (in the right filter):
    target_micado = sim_tp.star_grid(n=n,
                                     mmin=mmin,
                                     mmax=mmax,
                                     filter_name=target_filter[filt_ind],
                                     separation=100.*pixel_scale)
    
    
    # Position of the stars in the field
    offset = 0 # Well known bug in the creation of the image. This could be 1
    x_pos = target_micado.fields[0]["x"] / pixel_scale + npixx / 2  + offset
    y_pos = target_micado.fields[0]["y"] / pixel_scale + npixy / 2  + offset
    
    
    # Now we set the expoure time and then measure our limiting magnitudes:
    for idx, (dd, nn) in enumerate(zip(dit, ndit)):
        
        # Update NDIT, DIT
        micado.cmds["!OBS.dit"] = dd
        micado.cmds["!OBS.ndit"] = nn
        micado.update()
        
        # Observe the target:
        micado.observe(target_micado, update=True)
        # Read out the image
        hdus_micado = micado.readout()
        micado_im = hdus_micado[0][1].data


        # Aperture photometry
        flux, err, std, snr_values, sky = get_snr_aperture_photometry(image=micado_im, x_pos=x_pos, y_pos=y_pos,
                                             circle_radius=circle_radius,
                                             annulus_inner_radius=annulus_inner_radius,
                                             annulus_outer_radius=annulus_outer_radius)


        # Fit of mag vs SNR to determine where SNR=5.
        mask_nan = np.copy(np.isnan(np.log10(snr_values)))
        x_fit = np.log10(snr_values)[~np.transpose(mask_nan)[0]]
        y_fit = mags.value[~np.transpose(mask_nan)[0]] + vega_to_AB[filt_ind]

        # We can mask the low-SNR measurements in some fashion: a constant limit (not my favorite)
        # model = LinearRegression().fit(x_fit[y_fit<28], y_fit[y_fit<28])
        # or above a SNR threshold (e.g., SNR>5 to fit only robust measurements).
        model = LinearRegression().fit(x_fit[10**(x_fit.flatten())>5], y_fit[10**(x_fit.flatten())>5])
        # This would be useful only to plot the model for mag vs SNR
        y_pred = model.predict(x_fit)

        # This is the limiting magnitude based on the best-fit model:
        lim_mag[filt_ind, idx] = model.predict([[np.log10(sigma_lim)]])

        
# -

# # Print results and plot a figure

# +
# Print some results:
print('')
print('Limiting magnitude at '+str(round(sigma_lim,1))+'-sigma ('+refmag+' system):')

# Let's extrapolate the results at 5 hr as in Davies+2010 (see below):
nhr_lim = np.zeros(len(target_filter))
nhours = 5.

for ind, i in enumerate(target_filter):
    lim = np.interp(nhours*3600., ndit*dit, lim_mag[ind,:])
    nhr_lim[ind] = lim
    
    print(i+': '+str(round(lim,1))+' ('+str(round(nhours,1))+' hours)')
print('')


# Prepare the figure
fig = plt.figure(figsize=(8,5))
ax = fig.add_subplot(111)

# Colors:
clr = ['coral', 'red', 'firebrick']

# Generate x values for plotting the best-fit 1/sqrt(t) trend:
x_fit = np.linspace(min(dit*ndit / 3600.), max(dit*ndit / 3600.), 100)

for filt_ind, filt in enumerate(target_filter):
    
    # This scales as 1/sqrt(t)
    rms_sigma_lim =    (10**(-0.4*(lim_mag[filt_ind,:]-23.9))) / sigma_lim
    # Fit the function to the data using curve_fit
    params, covariance = curve_fit(sqrt_function,  dit*ndit / 3600., rms_sigma_lim)
    # Extract the fitted parameter
    fitted_a = params[0]
    # Calculate the fitted y values using the fitted parameter
    y_fit = sqrt_function(x_fit, fitted_a)

    
    #Scatter point
    ax.scatter(dit*ndit / 3600., lim_mag[filt_ind,:], color=clr[filt_ind])
    
    # Lines:
    # Connect points:
    #ax.plot(dit*ndit / 3600., lim_mag[filt_ind,:], color=clr[filt_ind], label=filt+' ('+str(round(nhr_lim[filt_ind],1))+', '+str(round(nhours,1))+'hr)')
    
    # Best-fit:
    ax.plot(x_fit, -2.5*np.log10(y_fit*sigma_lim)+23.9, color=clr[filt_ind], label=filt+' ('+str(round(nhr_lim[filt_ind],1))+', '+str(round(nhours,1))+'hr)')

# Set labels, legend
ax.set_xlabel('Exposure time (DIT x NDIT) [hours]')
ax.set_ylabel(r'Limiting magnitude at $5\sigma$')
ax.grid(alpha=0.3)
ax.legend(loc='lower right', frameon=False)
ax.set_title(mode_name[0]+' '+mode_name[1]+' - '+refmag)

# Manage ticks:
ax.yaxis.MinorTick='on'
ax.xaxis.MinorTick='on'
ax.yaxis.set_minor_locator(mtick.AutoMinorLocator(10))
ax.xaxis.set_minor_locator(mtick.AutoMinorLocator(10))
ax.tick_params(which='minor', length=3, direction='in')
ax.tick_params(which='major', length=6, direction='in')

# Save the figure:
fig.tight_layout();
fig.savefig('plots/performance_micado_'+mode_name[0]+'_'+mode_name[1]+'_'+refmag+'.png', format='png')
# -

# This is the plot in Davies+2010:

# ![Screenshot 2023-08-23 at 16.46.34.png](attachment:08459127-1c22-49f3-a020-8f6e95b65180.png)![Screenshot 2023-08-23 at 16.46.41.png](attachment:94933064-259a-4671-80a0-e0bcb5561ee7.png)

# +
#for i in micado.effects: print(i)
# -

micado.cmds