psf_simulation.py

import numpy as np
from astropy.io import fits

from hcipy import Field, Wavefront, DeformableMirror, FraunhoferPropagator
from hcipy import SurfaceApodizer, SurfaceAberrationAtDistance
from hcipy import make_pupil_grid, make_circular_aperture, make_focal_grid
from hcipy import make_xinetics_influence_functions
from hcipy import read_fits

from .config import cpism_refdata, S
from .optics import filter_throughput


# initial psf simulation
apm, apm_header = fits.getdata(
    f'{cpism_refdata}/optics/apm.fits', header=True)
actuator = read_fits(f'{cpism_refdata}/optics/actuator.fits')
surface_aberration = read_fits(
    f'{cpism_refdata}/optics/initial_phase_aberration.fits')

wavelength = 625e-9  # m
pupil_diameter = 2  # m
focal_length = pupil_diameter * 83
pupil_grid = make_pupil_grid(apm.shape[0], apm.shape[0] * apm_header['PHRATE'])
aperture = make_circular_aperture(pupil_diameter)(pupil_grid)
aperture = aperture * Field(apm.flatten(), pupil_grid)


emccd_pixel_size = 13e-6  # m
arcsec2rad = np.radians(1 / 3600)
emccd_pixel_scale = emccd_pixel_size / \
    (arcsec2rad * focal_length)  # arcsec / pixel
spatial_resolution = focal_length * wavelength / \
    pupil_diameter  # meter per airy disk

q = spatial_resolution / emccd_pixel_size  # pixels per airy disk
num_airy = 512 / q  # airy disk per frame (2 * 512 = 1024 pix)
focal_full_frame = make_focal_grid(
    q, num_airy, spatial_resolution=spatial_resolution)
prop_full_frame = FraunhoferPropagator(
    pupil_grid, focal_full_frame, focal_length)

num_airy = 128 / q  # make a small frame for the planets
focal_sub_frame = make_focal_grid(
    q, num_airy, spatial_resolution=spatial_resolution)
prop_sub_frame = FraunhoferPropagator(
    pupil_grid, focal_sub_frame, focal_length)

num_actuators_across = 32
# dm spacing is little smaller than pupil
actuator_spacing = 0.95 / 32 * pupil_diameter
influence_functions = make_xinetics_influence_functions(
    pupil_grid, num_actuators_across, actuator_spacing)
deformable_mirror = DeformableMirror(influence_functions)

aberration = SurfaceApodizer(
    surface_sag=surface_aberration.flatten(), refractive_index=-1)
# arbitrary distance for the aberration to propagate
aberration_distance = 80 * focal_length
aberration = SurfaceAberrationAtDistance(aberration, aberration_distance)


def single_band_psf(wavelength, waveerror=0, aber_phase=None):
    error = np.random.normal(0, waveerror*1e-9, actuator.shape)
    deformable_mirror.actuators = actuator + error
    wf = Wavefront(aperture, wavelength)
    wf = aberration(wf)
    if aber_phase is not None:
        other_error = SurfaceApodizer(
            surface_sag=aber_phase.flatten(), refractive_index=-1)
        wf = other_error(wf)
    img = prop_full_frame(deformable_mirror(wf)).intensity.shaped 
    return img

def simulate_psf(error_value, band, spectrum, nsample=5, cstar=True, aber_phase=None):
    filter = filter_throughput(band)
    wave = filter.wave
    throughput = filter.throughput
    min_wave = wave[0]
    max_wave = wave[-1]

    sed = []
    sed_center_wavelength = []
    for i_wave in range(nsample):
        d_wave = (max_wave - min_wave) / nsample
        wave0 = min_wave + i_wave * d_wave
        wave1 = min_wave + (i_wave + 1) * d_wave
        sed_center_wavelength.append((wave0 + wave1) / 2 * 1e-10)

        i_throughput = throughput.copy()
        i_throughput[(wave > wave1) | (wave < wave0)] = 0
        i_band = S.ArrayBandpass(wave, i_throughput, waveunits='angstrom')
        i_sed = (spectrum * i_band).integrate()
        sed.append(i_sed)

    error = np.random.normal(0, error_value*1e-9, actuator.shape)


    imgs = []
    deformable_mirror.actuators = actuator + error

    prop = prop_full_frame
    if not cstar:
        prop = prop_sub_frame

    for i_sample in range(nsample):
        wf = Wavefront(aperture, sed_center_wavelength[i_sample])
        wf = aberration(wf)
        if aber_phase is not None:
            other_error = SurfaceApodizer(
                surface_sag=aber_phase.flatten(), refractive_index=-1)
            wf = other_error(wf)

        img = prop(deformable_mirror(wf)).intensity.shaped
        imgs.append(img / img.sum() * sed[i_sample])

    return np.array(imgs).sum(axis=0)