From a8e88b2eefd239490f4aebd898b0b7653e50c579 Mon Sep 17 00:00:00 2001
From: Chen Yili <ylchen@niaot.ac.cn>
Date: Mon, 15 Apr 2024 03:34:04 +0000
Subject: [PATCH] Upload New File

---
 CpicImgSim/psf_simulation.py | 106 +++++++++++++++++++++++++++++++++++
 1 file changed, 106 insertions(+)
 create mode 100644 CpicImgSim/psf_simulation.py

diff --git a/CpicImgSim/psf_simulation.py b/CpicImgSim/psf_simulation.py
new file mode 100644
index 0000000..0807813
--- /dev/null
+++ b/CpicImgSim/psf_simulation.py
@@ -0,0 +1,106 @@
+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 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)
-- 
GitLab