different dm volt can be choose for each exposure

README.md

@@ -8,3 +8,8 @@ CpicImgSim is a Python library for numerical simulation of the observational ima
 数据文件：https://pan.cstcloud.cn/s/gjMlzXRuRsM	提取密码：i3Zs
 程序文档：https://csst-tb.bao.ac.cn/simulation/cpic/index.html
 联系方式：南京天光所 赵刚 gzhao@niaot.ac.cn 025-85482316
+
+## 软件更新日志
+2025年6月5日：添加了新的宇宙线生成模块，使用主巡天内的模块代码。新的模块可以更好的兼容通用的宇宙线识别模块。设置方法：在run_config 里添加msc_cr_model: true条目可以启用新的宇宙线生成模块。
+
+2025年6月10日：增加了波前电压设置功能。现在对于每一次曝光可以设置不同的波前电压。在观测yaml中增加dm_volt_index: xxx选项。xxx为整形，例如565, 661等。另外支持1, 2, 3, 4选项，为四个不同方向的暗区。实验功能：dm_volt_index可设置-1，自动选择和波段匹配的暗区（该功能谨慎使用）。

csst_cpic_sim/camera.py

@@ -4,13 +4,13 @@ import numpy as np
 import scipy.ndimage as nd
 from astropy.io import fits
 
-from .config import config, S
+from .config import config, S, cpism_refdata
 from .utils import region_replace, random_seed_select
 from .io import log
 from .optics import filter_throughput
 from .msc_cr import produceCR_Map
 
-cpism_refdata = config['cpism_refdata']
+# cpism_refdata = config['cpism_refdata']
 MAG_SYSTEM = config['mag_system']
 solar_spectrum = S.FileSpectrum(config['solar_spectrum'])
 solar_spectrum.convert('photlam')
@@ -260,7 +260,7 @@ class CosmicRayFrameMaker(object):
                 exTime=expt,
                 cr_pixelRatio=0.003*expt/600.,  # 5e-6/pixel , 1.3e-6 for cpic package
                 gain=1,  # gain = 1, means the unit of output image is electron
-                attachedSizes=self.attachedSizes)
+                attachedSizes=self.attachedSize)
             return cr_map
 
         for i in range(len(cr_array)):
@@ -881,8 +881,6 @@ class CpicVisEmccd(object):
         em_fix = self.em_fix_fuc_fit(emgain) * emgain
         image = np.random.gamma(image, em_fix) + image * (emgain - em_fix)
 
-        fits.writeto("image.fits", image, overwrite=True)
-
         if self.switch['em_blooming']:
             image = self.emregester_blooming(image)
 

csst_cpic_sim/config.py

@@ -51,21 +51,24 @@ def load_refdata_path(config_aim):
     refdata = os.path.abspath(refdata)
     if os.path.isdir(refdata):
         refdata_list.append(refdata)
-    with open(config_aim, 'w') as f:
-        yaml.dump(refdata_list, f)
+        with open(config_aim, 'w') as f:
+            yaml.dump(refdata_list, f)
+        print(f"refdata path successfully set to be {refdata}")
+    else:
+        print(f"{refdata} is not a valid path")
     exit()
 
 
 cpism_refdata = load_refdata_path(config_aim)
 
 config = {}
-config['cpism_refdata'] = cpism_refdata
+# config['cpism_refdata'] = cpism_refdata
 config['utc0'] = '2024-05-01T00:00:00'
 config['hybrid_model'] = f'{cpism_refdata}/target_model/hybrid_model.fits'
 config['bcc_model'] = f'{cpism_refdata}/target_model/bccmodels'
 config['mag_system'] = 'abmag'
 config['apm_file'] = f'{cpism_refdata}/optics/apm.fits'
-config['actuator_file'] = f'{cpism_refdata}/optics/actuator.fits'
+
 config['aberration'] = f'{cpism_refdata}/optics/initial_phase_aberration.fits'
 config['mask_width'] = 0.4
 config['check_fits_header'] = False
@@ -93,10 +96,25 @@ config['catalog_folder'] = f'{cpism_refdata}/demo_catalog'
 config['csst_format'] = True
 config['nsample'] = 1
 config['msc_cr_model'] = False
+config['dm_volt_group'] = {
+    1: f'{cpism_refdata}/optics/dark_zone_0011_volt.fits',
+    2: f'{cpism_refdata}/optics/dark_zone_0110_volt.fits',
+    3: f'{cpism_refdata}/optics/dark_zone_1100_volt.fits',
+    4: f'{cpism_refdata}/optics/dark_zone_1001_volt.fits',
+    565: f'{cpism_refdata}/optics/dark_zone_1100_F565_volt.fits',
+    661: f'{cpism_refdata}/optics/dark_zone_1100_F661_volt.fits',
+    743: f'{cpism_refdata}/optics/dark_zone_1100_F743_volt.fits',
+    883: f'{cpism_refdata}/optics/dark_zone_1100_F883_volt.fits',
+    520: f'{cpism_refdata}/optics/dark_zone_1100_F520_volt.fits',
+    662: f'{cpism_refdata}/optics/dark_zone_1100_F662_volt.fits',
+    850: f'{cpism_refdata}/optics/dark_zone_1100_F850_volt.fits',
+    720: f'{cpism_refdata}/optics/dark_zone_1100_F720_volt.fits',
+}
+config['wavefront_error'] = 0.075  # nm
 
 update_able_keys = [
-    'apm_file', 'actuator_file', 'aberration',
-    'log_dir', 'log_level', 'catalog_folder', 'msc_cr_model'
+    'apm_file', 'aberration', 'wavefront_error',
+    'log_dir', 'log_level', 'catalog_folder', 'msc_cr_model',
     'nsample', 'csst_format', 'output', 'check_fits_header'
 ]
 

csst_cpic_sim/main.py

@@ -10,7 +10,6 @@ from datetime import datetime
 import traceback
 
 import numpy as np
-from astropy.io import fits
 
 from .target import spectrum_generator, target_file_load
 from .optics import focal_mask, focal_convolve
@@ -18,7 +17,7 @@ from .camera import CosmicRayFrameMaker, sky_frame_maker, CpicVisEmccd
 from .io import save_fits, log
 from .config import update_able_keys, relative_time
 from .config import config as default_config, replace_cpism_refdata
-from .psf_simulation import actuator
+# from .psf_simulation import actuators
 
 
 def vis_observation(
@@ -39,7 +38,8 @@ def vis_observation(
         emgain=None,
         prograss_bar=None,
         output='./',
-        dataset='default'
+        dataset='default',
+        dm_volt_index=1,
 ) -> np.ndarray:
     """Observation simulation main process on visable band using EMCCD.
 
@@ -87,6 +87,7 @@ def vis_observation(
     platescale = default_config['platescale']
     iwa = default_config['mask_width'] / 2
     area = default_config['aperature_area']
+    error = default_config['wavefront_error']
     expt_start = camera.system_time
 
     target_list = []
@@ -133,7 +134,9 @@ def vis_observation(
             init_shifts=shift,
             rotation=rotation,
             nsample=nsample,
-            platesize=camera.flat_shape
+            error=error,
+            platesize=camera.flat_shape,
+            dm_volt_index=dm_volt_index
         )
 
         if skybg is None or skybg > 100:
@@ -330,9 +333,9 @@ def observation_simulation_from_config(obs_file, config_file):
 
     deduplicate_names_add_count(all_camera_name)
 
-    if 'actuator_file' in config.keys():
-        actuator_data = fits.getdata(config['actuator_file'])
-        actuator[:] = actuator_data[:]
+    # if 'actuator_file' in config.keys():
+    #     actuator_data = fits.getdata(config['actuator_file'])
+    #     actuator[:] = actuator_data[:]
 
     for key in config.keys():
         if key in update_able_keys:
@@ -377,6 +380,7 @@ def observation_simulation_from_config(obs_file, config_file):
             emgain = obs_info.get('emgain', None)
             time = relative_time(time)
             dataset = obs_info.get('dataset', 'default')
+            dm_volt_index = obs_info.get('dm_volt_index', 1)
 
         except Exception as e:
             log.error(f"{file} is not a valid yaml file.")
@@ -424,7 +428,9 @@ def observation_simulation_from_config(obs_file, config_file):
                     csst_format=csst_format,
                     crmaker=crmaker,
                     dataset=dataset,
-                    prograss_bar=True)
+                    prograss_bar=True,
+                    dm_volt_index=dm_volt_index,
+                    )
             except Exception as e:
                 log.error(f"{info_text} failed with {type(e).__name__}{e}.\n\n {traceback.format_exc()}")
 

csst_cpic_sim/optics.py

@@ -133,7 +133,8 @@ def focal_convolve(
         rotation: float = 0,
         nsample: int = 5,
         error: float = 0,
-        platesize: list = [1024, 1024]) -> np.ndarray:
+        platesize: list = [1024, 1024],
+        dm_volt_index: int = 1) -> np.ndarray:
     """PSF convolution of the ideal focus image.
 
     Parameters
@@ -184,6 +185,9 @@ def focal_convolve(
     if not targets:
         return np.zeros((platesize[1], platesize[0]))
 
+    if dm_volt_index == -1:  # experimental function, be careful to use it.
+        dm_volt_index = int(band[1:])
+
     for i_wave in range(nsample):
         d_wave = (max_wave - min_wave) / nsample
         wave0 = min_wave + i_wave * d_wave
@@ -195,11 +199,15 @@ def focal_convolve(
         i_band = S.ArrayBandpass(wave, i_throughput, waveunits='angstrom')
 
         i_fp_image = ideal_focus_image(i_band, targets[1:], platescale, platesize, init_shifts, rotation)
-        psf = single_band_psf(center_wavelength, error=error)
+        psf = single_band_psf(center_wavelength, error=error, actuator_index=dm_volt_index)
 
         _, _, cstar_sp, _ = targets[0]
         cstar_flux = (cstar_sp * i_band).integrate()
-        cstar_psf = single_band_masked_psf(center_wavelength, error=error, shift=init_shifts)
+        cstar_psf = single_band_masked_psf(
+            center_wavelength,
+            error=error,
+            shift=init_shifts,
+            actuator_index=dm_volt_index)
 
         c_fp_image = fftconvolve(i_fp_image, psf, mode='same')
         c_fp_image = focal_mask(c_fp_image, iwa, platescale)

csst_cpic_sim/psf_simulation.py

@@ -17,7 +17,13 @@ ARCSEC2RAD = np.radians(1 / 3600)
 
 # initial psf simulation
 apm, apm_header = fits.getdata(config['apm_file'], header=True)
-actuator = fits.getdata(config['actuator_file'])
+
+actuator_files = config['dm_volt_group']
+actuators = {}
+for key, actuator_file in actuator_files.items():
+    actuator_data = fits.getdata(actuator_file)
+    actuators[key] = actuator_data
+
 surface_aberration = fits.getdata(config['aberration'])
 
 pupil_diameter = 2  # m
@@ -78,7 +84,8 @@ aberration = SurfaceAberrationAtDistance(aberration, aberration_distance)
 def single_band_masked_psf(
         wavelength: float,
         error: float = 0,
-        shift: list = [0, 0]) -> np.ndarray:
+        shift: list = [0, 0],
+        actuator_index: int = 1) -> np.ndarray:
     """CPIC PSF considering the focal plane mask.
 
     Parameters
@@ -89,13 +96,15 @@ def single_band_masked_psf(
         deformable mirror control error in nm
     shift : list
         angular shift of the target in arcsec.
-
+    actuator_index : int
+        actuator index for the observation.
     Returns
     ----------
     psf : np.ndarray
         psf in the focal plane. Normalized as the input flux is 1.
         (Note that total flux of the psf is not 1, because it is masked)
     """
+    actuator = actuators[actuator_index]
     error = np.random.normal(0, error*1e-9, actuator.shape)
     deformable_mirror.actuators = actuator + error
     wf = Wavefront(aperture, wavelength)
@@ -117,7 +126,8 @@ def single_band_masked_psf(
 def single_band_shift_psf(
         wavelength: float,
         error: float = 0,
-        shift: list = [0, 0]) -> np.ndarray:
+        shift: list = [0, 0],
+        actuator_index: int = 1) -> np.ndarray:
     """CPIC PSF considering the focal plane mask.
 
     Parameters
@@ -135,6 +145,7 @@ def single_band_shift_psf(
         psf in the focal plane. Normalized as the input flux is 1.
         (Note that total flux of the psf is not 1, because it is masked)
     """
+    actuator = actuators[actuator_index]
     error = np.random.normal(0, error*1e-9, actuator.shape)
     deformable_mirror.actuators = actuator + error
     wf = Wavefront(aperture, wavelength)
@@ -150,7 +161,8 @@ def single_band_shift_psf(
 
 def single_band_psf(
         wavelength: float,
-        error: float = 0) -> np.ndarray:
+        error: float = 0,
+        actuator_index: int = 1) -> np.ndarray:
     """CPIC PSF without considering the focal plane mask.
     Used to simulate the off-axis PSF.
 
@@ -166,6 +178,7 @@ def single_band_psf(
     psf : np.ndarray
         psf in the focal plane. Normalized. The total flux is 1.
     """
+    actuator = actuators[actuator_index]
     error = np.random.normal(0, error*1e-9, actuator.shape)
     deformable_mirror.actuators = actuator + error
     wf = Wavefront(aperture, wavelength)

example/observations/05_sci.yaml

-band: f661
+band: f520
 # emgain: 100
 emset: 180
 expt: 30
@@ -10,6 +10,7 @@ shift:
 - 0
 - 0
 skybg: 20.31
+dm_volt_index: -1
 target:
   cstar:
     dec: 84.2875d