Skip to content
GitLab
Commit 95f81f91 authored by GZhao's avatar GZhao
Browse files

v2.0beta update

parent b7f8c4fa
Pipeline #4357 failed with stage
in 0 seconds
Show whitespace changes
Inline Side-by-side
.gitignore
View file @ 95f81f91
......@@ -8,6 +8,10 @@ dist/
starmodel/
cpism_refdata/
*.egg-info
example/example_output
refdata/starmodel
refdata/target_model
# Other files and folders
.settings/
......@@ -26,6 +30,7 @@ docs/notebooks/image_files/_*
# unitest output
tests/.coverage
tests/htmlcov/
tests/*.xml
# Executables
*.swf
......
.gitignore copy 0 → 100644
View file @ 95f81f91
# Build and Release Folders
bin-debug/
bin-release/
[Oo]bj/
[Bb]in/
build/
dist/
starmodel/
cpism_refdata/
*.egg-info
# Other files and folders
.settings/
_*/
~*
.VSCodeCounter
.vscode
*.log
output/
*.log.*
/*.fits
testrun/
example/example_output/
docs/notebooks/image_files/_*
# unitest output
tests/.coverage
tests/htmlcov/
# Executables
*.swf
*.air
*.ipa
*.apk
# Project files, i.e. `.project`, `.actionScriptProperties` and `.flexProperties`
# should NOT be excluded as they contain compiler settings and other important
# information for Eclipse / Flash Builder.
CpicImgSim/__init__.py deleted 100644 → 0
View file @ b7f8c4fa
# from .main import quick_run, observation_simulation
from .optics import make_focus_image, focal_mask
from .target import star_photlam, planet_contrast, extract_target_x_y, spectrum_generator
from .camera import EMCCD, CosmicRayFrameMaker, sky_frame_maker
from .config import __version__
__all__ = [
"EMCCD",
"CosmicRayFrameMaker",
"sky_frame_maker",
"star_photlam",
"planet_contrast",
"extract_target_x_y",
"spectrum_generator",
"make_focus_image",
"focal_mask",
# "quick_run",
# "observation_simulation"
]
\ No newline at end of file
CpicImgSim/camera.py deleted 100644 → 0
View file @ b7f8c4fa
import os
import yaml
import numpy as np
import scipy.ndimage as nd
from astropy.io import fits
import matplotlib.pyplot as plt
import math
from CpicImgSim.config import cpism_refdata, solar_spectrum, MAG_SYSTEM
from CpicImgSim.utils import region_replace, random_seed_select
from CpicImgSim.io import log
from CpicImgSim.optics import filter_throughput
def sky_frame_maker(band, skybg, platescale, shape):
"""
generate a sky background frame.
Parameters
----------
band : str
The band of the sky background.
skybg : str
The sky background file name.
platescale : float
The platescale of the camera in arcsec/pixel.
shape : tuple
The shape of the output frame. (y, x)
Returns
-------
sky_bkg_frame : numpy.ndarray
The sky background frame.
"""
filter = filter_throughput(band)
sk_spec = solar_spectrum.renorm(skybg, MAG_SYSTEM, filter)
sky_bkg_frame = np.zeros(shape)
sky_bkg_frame += (sk_spec * filter).integrate() * platescale**2
return sky_bkg_frame
class CRobj(object):
"""
Cosmic ray object.
Attributes
----------
flux : float
The flux of the cosmic ray.
angle : float
The angle of the cosmic ray.
sigma : float
The width of the cosmic ray.
length : int
The length of the cosmic ray.
"""
def __init__(self, flux, angle, sigma, length) -> None:
self.flux = flux
self.angle = angle
self.sigma = sigma
self.length = length
class CosmicRayFrameMaker(object):
"""
Cosmic ray frame maker.
Parameters
----------
depth : float
The depth of the camera pixel in um.
pitch : float
The pitch of the camera pixel in um.
cr_rate : float
The cosmic ray rate per second per cm2.
"""
def __init__(self) -> None:
self.tmp_size = [7, 101]
self.freq_power = -0.9
self.trail_std = 0.1
self.depth = 10 # um
self.pitch = 13 # um
self.cr_rate = 1 # particle per s per cm2 from Miles et al. 2021
def make_CR(self, length, sigma, seed=-1):
"""
make a image of cosmic ray with given length and sigma.
Parameters
----------
length : int
The length of the cosmic ray in pixel.
sigma : float
The width of the cosmic ray in pixel.
Returns
-------
output : numpy.ndarray
The image of cosmic ray.
"""
h = self.tmp_size[0]
w = self.tmp_size[1]
freq = ((w-1)/2-np.abs(np.arange(w)-(w-1)/2)+1)**(self.freq_power)
x = np.arange(w) - (w-1)/2
hl = (length-1)/2
x_wing = np.exp(-(np.abs(x)-hl)**2/sigma/sigma/2)
x_wing[np.abs(x) < hl] = 1
cr = np.zeros([h, w])
center = (h-1)/2
for i in range(h):
phase = np.random.rand(w)*2*np.pi
trail0 = abs(np.fft.fft(freq*np.sin(phase) + 1j*x*np.cos(phase)))
# TODO maybe somthing wrong
trail_norm = (trail0 - trail0.mean())/trail0.std()
cr[i, :] = np.exp(-(i - center)**2/sigma/sigma/2) \
* (trail_norm * self.trail_std + 1) * x_wing
output = np.zeros([w, w])
d = (w-h)//2
output[d:d+h, :] = cr
return output
def _length_rand(self, N, seed=-1):
"""
randomly generate N cosmic ray length.
"""
len_out = []
seed = random_seed_select(seed=seed)
log.debug(f"cr length seed: {seed}")
for i in range(N):
x2y2 = 2
while x2y2 > 1:
lx = 1 - 2 * np.random.rand()
ly = 1 - 2 * np.random.rand()
x2y2 = lx * lx + ly * ly
z = 1 - 2 * x2y2
r = 2 * np.sqrt(x2y2 * (1 - x2y2))
length = abs(r / z * self.depth)
pitch = self.pitch
len_out.append(int(length/pitch))
return np.array(len_out)
def _number_rand(self, expt, pixsize, random=False, seed=-1):
"""
randomly generate the number of cosmic rays.
"""
area = (self.pitch / 1e4)**2 * pixsize[0] * pixsize[1]
ncr = area * expt * self.cr_rate
if random:
seed = random_seed_select(seed=seed)
log.debug(f"cr count: {seed}")
ncr = np.random.poisson(ncr)
else:
ncr = int(ncr)
self.ncr = ncr
return ncr
def _sigma_rand(self, N, seed=-1):
"""
randomly generate N cosmic ray sigma.
"""
sig_sig = 0.5 # asuming the sigma of size of cosmic ray is 0.5
seed = random_seed_select(seed=seed)
log.debug(f"cr width seed: {seed}")
sig = abs(np.random.randn(N))*sig_sig + 1/np.sqrt(12) * 1.2
# assume sigma is 1.2 times of pictch sig
return sig
def _flux_rand(self, N, seed=-1):
"""
randomly generate N cosmic ray flux.
"""
seed = random_seed_select(seed=seed)
log.debug(f"cr flux seed: {seed}")
u = np.random.rand(N)
S0 = 800
lam = 0.57
S = (-np.log(1-u)/lam + S0**0.25)**4
return S
def random_CR_parameter(self, expt, pixsize):
"""
randomly generate cosmic ray parameters, including number, length, flux, sigma and angle.
Parameters
----------
expt : float
The exposure time in second.
pixsize : list
The size of the image in pixel.
Returns
-------
CRs : list
A list of cosmic ray objects.
"""
N = self._number_rand(expt, pixsize)
log.debug(f"cr count: {N}")
length = self._length_rand(N)
if N > 0:
log.debug(f"cr length, max: {length.max()}, min: {length.min()}")
flux = self._flux_rand(N)
log.debug(f"cr flux, max: {flux.max()}, min: {flux.min()}")
sig = self._sigma_rand(N)
log.debug(f"cr width, max: {sig.max()}, min: {sig.min()}")
seed = random_seed_select(seed=-1)
log.debug(f"cr angle seed: {seed}")
angle = np.random.rand(N) * 180
CRs = []
for i in range(N):
CRs.append(CRobj(flux[i], angle[i], sig[i], length[i]))
return CRs
def make_cr_frame(self, shape, expt, seed=-1):
"""
make a cosmic ray frame.
Parameters
----------
shape : list
The size of the image in pixel.
expt : float
The exposure time in second.
seed : int, optional
The random seed. The default is -1. If seed is -1, the seed will be randomly selected.
Returns
-------
image : numpy.ndarray
The cosmic ray frame.
"""
image = np.zeros(shape)
sz = shape
cr_array = self.random_CR_parameter(expt, shape)
cr_center = (self.tmp_size[1] - 1)/2
seed = random_seed_select(seed=seed)
log.debug(f"cr position seed: {seed}")
for i in range(len(cr_array)):
cr = cr_array[i]
x = np.random.rand() * sz[1]
y = np.random.rand() * sz[0]
cr_img = self.make_CR(cr.length, cr.sigma)
cr_img *= cr.flux
cr_img = abs(nd.rotate(cr_img, cr.angle, reshape=False))
if i == 0:
pdin = False
else:
pdin = True
if i == len(cr_array) - 1:
pdout = False
else:
pdout = True
image = region_replace(
image, cr_img,
[y-cr_center, x-cr_center],
padded_in=pdin,
padded_out=pdout, subpix=True
)
image = np.maximum(image, 0)
log.debug(f"cr image max: {image.max()}, min: {image.min()}")
return image
class CpicVisEmccd(object):
def __init__(self):
self._defaut_config()
self.config_init()
@classmethod
def from_config(cls, config):
obj = cls.__new__()
obj._defaut_config()
obj.__dict__.update(config)
obj._refresh_config()
return obj
def _defaut_config(self):
self.plszx = 1024
self.plszy = 1024
self.pscan1 = 8
self.pscan2 = 0
self.oscan1 = 16
self.oscan2 = 18
self.udark = 6
self.bdark = 2
self.ldark = 16
self.rdark = 16
self.max_adu = 16_383
self.switch = {
'flat': True,
'dark': True,
'stripe': True,
'cic': True,
'cte': True,
'badcolumn': True,
'nonlinear': True,
'cosmicray': True,
'blooming': False,
'bias_vp': True,
'bias_hp': True,
'bias_ci': True,
'bias_shift': True,
}
self.dark_file = cpism_refdata + '/camera/emccd_dark_current.fits'
self.flat_file = cpism_refdata + '/camera/emccd_flat_field.fits'
self.cic_file = cpism_refdata + '/camera/emccd_cic.fits'
self.bad_col_file = cpism_refdata + '/camera/emccd_bad_columns.fits'
self.cic = 0.2
self.fullwell = 80_000
self.em_fullwell = 500_000 #780_000
self.em_cte = 0.9996
self.emreg_cal_num = 10 # 用来加速计算
self.emreg_num = 604
self.readout_speed = 6.25e6 # Hz
self.readout_time = 0.365 # s
self.heat_speed = 1000 # voltage / 1000 degree per frame
self.temper_speed = 0.05 # degree per second
self.cooler_temp = -80
self.readout_noise = 160
self.ph_per_adu = 59
self.bias_level = 200
self.vertical_param = [0, 14]
self.horizontal1_param = [5.3/2, 12, 5/2, 50]
self.horizontal2_param = [5.3/4, 16.17, 5/4, 76.65]
self.bias_hp_resd = 2
self.cooler_interfence = 20
self.bias_level_std = 3
self.bias_shift_per_volt = -1.34
self.shift_time = 1 / 1000
def shutter_effect(self, image, expt):
image_mean = image.mean(axis=0)
image_shutter = np.broadcast_to(image_mean, image.shape)
image_shutter /= image_shutter.mean()
image_shutter *= image.mean() * self.shift_time / expt
return image_shutter.astype(int)
def config_init(self):
self._img_index = 0
self._ccd_temp = self.cooler_temp
self._vertical_i0 = np.random.randint(0, 2)
self._frame_read_time = 2080 * 1055 / 6.25e6
self.flat_shape = [self.plszy, self.plszx]
darksz_x = self.plszx + self.rdark + self.ldark
darksz_y = self.plszy + self.udark + self.bdark
self.dark_shape = [darksz_y, darksz_x]
biassz_x = darksz_x + self.pscan1 + self.oscan1
biassz_y = darksz_y + self.pscan2 + self.oscan2
self.image_shape = [biassz_y, biassz_x]
self.flat = fits.getdata(self.flat_file)
# self.cic = fits.getdata(self.cic_file)
self.dark = fits.getdata(self.dark_file)
self.bad_col = fits.getdata(self.bad_col_file)
self.ccd_temp = self.cooler_temp
self.system_time = 0
self.time_syn(0, initial=True)
self.emgain_set(1024, -30)
def CTE_cal(n, N_p, CTE, S_0):
'''
CTE_cal(order of trail pixels, number of pixel transfers, CTE, initial intensity of target pixel)
'''
CTI = 1 - CTE
S_0_n = S_0 * ((N_p * CTI) ** n) / math.factorial(n) * math.exp(-N_p * CTI)
return S_0_n
def cte_201(cte, start=0, length=10):
N_p = 604
S_0 = 1
res = [0]*start
for n in range(length):
s = CTE_cal(n, N_p, cte, S_0)
res.append(s)
return np.array(res)
cti_trail = cte_201(self.em_cte, start=0, length=10)
self.cti_trail = cti_trail / cti_trail.sum()
def bias_frame(self, show=False):
shape = self.image_shape
TPI = np.pi * 2
# vertical pattern
# 使用一维的曲线描述竖条纹的截面
vp_1d = np.zeros(shape[1])
# 以下代码用于模拟垂直间隔的竖条纹在奇偶幅图像时的不同表现
# 后续相机更新过程中，已经将改特性进行修改，因此不再使用此代码
# 这个代码写的还是有点意思，所以保留下来了
# vp_1d[0::2] = self.vertical_param[self._vertical_i0]
# self._vertical_i0 = 1 - self._vertical_i0
# vp_1d[1::2] = self.vertical_param[self._vertical_i0]
vp_1d[0::2] = self.vertical_param[self._vertical_i0]
vp_1d[1::2] = self.vertical_param[1 - self._vertical_i0]
vp_frame = np.zeros(shape)
if self.switch['bias_vp']:
vp_frame += vp_1d
if show: # pragma: no cover
plt.figure(figsize=(10, 3))
plt.plot(vp_1d)
plt.xlim([0, 100])
plt.xlabel('x-axis')
plt.ylabel('ADU')
plt.title('vertical pattern')
# horizontal pattern
# 图像上的横条纹是梳状，分为两个部分，左边大约77个像素是周期小一点，其余的会大一点
boundary = 77 # boundary between left and width
boundary_width = 5 # 左右需要平滑过度一下
y = np.arange(self.image_shape[0])
hp_left_param = self.horizontal1_param # 实测数据拟合得到的
hp_left_1d = hp_left_param[0] * np.sin(TPI * (y / hp_left_param[1] + np.random.rand()))
hp_left_1d += hp_left_param[2] * np.sin(TPI * (y / hp_left_param[3] + np.random.rand()))
hp_left_frame = np.broadcast_to(hp_left_1d, [boundary+boundary_width, len(hp_left_1d),]).T
hp_right_param = self.horizontal2_param
hp_right_1d = hp_right_param[0] * np.sin(TPI * (y / hp_right_param[1] + np.random.rand()))
hp_right_1d += hp_right_param[2] * np.sin(TPI * (y / hp_right_param[3] + np.random.rand()))
hp_right_frame = np.broadcast_to(hp_right_1d, [shape[1] - boundary, len(hp_right_1d)]).T
combine_profile_left = np.ones(boundary + boundary_width)
combine_profile_left[-boundary_width:] = (boundary_width - np.arange(boundary_width) - 1) / boundary_width
combine_profile_right = np.ones(shape[1] - boundary)
combine_profile_right[:boundary_width] = np.arange(boundary_width)/boundary_width
hp_frame = np.zeros(shape)
hp_frame[:, :boundary+boundary_width] = hp_left_frame * combine_profile_left
hp_frame[:, boundary:] = hp_right_frame * combine_profile_right
# residual frame 横条纹外，还有一个垂直方向的梯度，根据测试数据，使用直线加指数函数的方式来拟合
exp_a, exp_b, lin_a, lin_b = (-1.92377463e+00, 1.32698365e-01, 8.39509583e-04, 4.25384480e-01)
y_trend = exp_a * np.exp(-(y + 1) * exp_b) + y * lin_a - lin_b
# random horizontal pattern generated in frequency domain
# 除了规则横条纹外，还有随机的横条纹，随机的横条纹在频域空间生成，具有相同的频率谱和随机的相位
rsd_freq_len = len(y_trend) * 4
red_freq = np.arange(rsd_freq_len)
red_freq = red_freq - (len(red_freq) - 1)/2
red_freq = np.exp(-red_freq**2 / 230**2) * 240 + np.random.randn(rsd_freq_len)*30
red_freq = np.fft.fftshift(red_freq)
phase = np.random.rand(rsd_freq_len) * TPI
hp_rsd_1d = np.fft.ifft(np.exp(1j * phase) * red_freq)
hp_rsd_1d = np.abs(hp_rsd_1d) * self.bias_hp_resd
hp_rsd_1d = hp_rsd_1d[:rsd_freq_len//4]
hp_rsd_1d = hp_rsd_1d - np.mean(hp_rsd_1d)
hp_rsd_1d = y_trend + hp_rsd_1d
hp_frame = (hp_frame.T + hp_rsd_1d).T
if not self.switch['bias_hp']:
hp_frame *= 0
if show:
plt.figure(figsize=(10, 3))
# plt.plot(hp_right_1d)
plt.plot(hp_rsd_1d)
# plt.xlim([0, 200])
plt.xlabel('y-axis')
plt.ylabel('ADU')
plt.title('vertical pattern')
# 接上制冷机后，会有亮暗点
#cooler interfence effect
ci_position = 10
ci_sub_struct = 80
ci_sub_exp = 2.5
ci_x_shft = 3
ci_interval = 250 # 6.25MHz readout / 2.5KHz interfence
ci_dn = self.cooler_interfence
npix = shape[0] * shape[1]
n_ci_event = npix // ci_interval
ci_align = np.zeros((n_ci_event, ci_interval))
ci_align[:, ci_position] = np.random.randn(n_ci_event) * ci_dn
ci_align[:, ci_position+1] = np.random.randn(n_ci_event) * ci_dn
yi0 = np.random.randint(0, ci_sub_struct)
xs0 = (ci_interval - ci_position) / (ci_sub_struct / 2)**ci_sub_exp
for yi in range(n_ci_event):
sub_yi = (yi - yi0) % ci_sub_struct
sub_yi = abs(sub_yi - ci_sub_struct / 2)
shiftx = int(sub_yi**ci_sub_exp * xs0)
ci_align[yi, :] = np.roll(ci_align[yi, :], shiftx)
ci_align = np.pad(ci_align.flatten(), (0, npix-n_ci_event*ci_interval))
ci_frame = ci_align.reshape(shape[0], shape[1])
for yi in range(shape[0]):
ci_frame[yi, :] = np.roll(ci_frame[yi, :], yi * ci_x_shft)
if not self.switch['bias_ci']:
ci_frame *= 0
bias_shift = 0
if self.switch['bias_shift']:
bias_shift = (self.volt - 25) * self.bias_shift_per_volt
# 混合在一起
rn_adu = self.readout_noise / self.ph_per_adu
bias_level = self.bias_level + np.random.randn() * self.bias_level_std
bias_frame = vp_frame + ci_frame + hp_frame + bias_level
bias_frame += rn_adu * np.random.randn(shape[0], shape[1]) + bias_shift
return bias_frame
def nonlinear_effect(self, image):
"""
nonlinear effect
"""
fullwell = self.fullwell
nonlinear_coefficient = 0.1
log.debug(
f"nonlinear effect added with coefficient {nonlinear_coefficient}")
image += (image / fullwell)**2 * nonlinear_coefficient * fullwell
return image
def time_syn(self, t, readout=True, initial=False):
if initial:
self.ccd_temp = self.cooler_temp
self.system_time = t
return
dt = t
heat = 0
if readout:
heat = self.volt / self.heat_speed
self.ccd_temp = heat + self.cooler_temp + (self.ccd_temp - self.cooler_temp) * np.exp(-dt * self.temper_speed)
if self.ccd_temp < self.cooler_temp: #
self.ccd_temp = self.cooler_temp
self.system_time += dt
# def em_cte(self, img):
# i_shift = 1
# cte_coe = 0.1
# img_shift_i = np.zeros_like(img)
# img_shift_i = np.random.poisson(img * cte_coe)
# pass
def emgain_set(self, em_set, ccd_temp=None, self_update=True):
if ccd_temp is None:
ccd_temp = self.ccd_temp
volt_coe_a = -0.01828
volt_coe_b = 43.61
volt_func = lambda es: volt_coe_a * es + volt_coe_b
self.volt = volt_func(em_set)
volt_3ff = volt_func(int('3ff', 16))
volt_190 = volt_func(int('190', 16))
em_coe_c = 0.24
# using the expression of em_b = ln(g199) / constant to make fitting easier
constant = (np.exp(em_coe_c * volt_190) - np.exp(em_coe_c * volt_3ff))
# fitting from the ccd test result
ln_g190 = (-ccd_temp - 7) * 0.0325
em_coe_b = ln_g190 / constant
emgain = np.exp(em_coe_b * np.exp(em_coe_c * self.volt))
emgain = np.maximum(1, emgain)
# print(emgain, em_coe_b, em_coe_c * self.volt, self.volt, np.exp(em_coe_c * self.volt))
if self_update:
self.emgain = emgain
self.emset = em_set
return emgain
def add_em_cte_conv(self, image):
shape = image.shape
img_line = np.convolve(image.flatten(), self.cti_trail, mode='full')
return img_line[:shape[0]*shape[1]].reshape(shape)
def readout(self, image_focal, em_set, expt_set, image_cosmic_ray=False):
expt = expt_set
if expt_set == 0:
expt = 0.001
dt = self.readout_time + expt
self.time_syn(dt, readout=True)
emgain = self.emgain_set(em_set, self.ccd_temp)
image = image_focal
if self.switch['flat']:
image = image * self.flat
if self.switch['nonlinear']:
image = self.nonlinear_effect(image)
darksz_x = self.plszx + self.rdark + self.ldark
darksz_y = self.plszy + self.udark + self.bdark
img_dark = np.zeros((darksz_y, darksz_x))
img_dark[
self.bdark:self.plszy+self.bdark,
self.ldark:self.ldark+self.plszx
] = image
image = img_dark
if self.switch['dark']:
image += self.dark
image *= expt
if self.switch['nonlinear']:
image = self.nonlinear_effect(image)
if self.switch['cic']:
image += self.cic
if image_cosmic_ray is not None and self.switch['cosmicray']:
image += image_cosmic_ray
if self.switch['blooming']:
image = self.vertical_blooming(image)
if self.switch['badcolumn']:
for i in range(self.bad_col.shape[1]):
deadpix_x = self.bad_col[0, i]
deadpix_y = self.bad_col[1, i]
image[deadpix_y:, deadpix_x] = 0
biassz_x = darksz_x + self.pscan1 + self.oscan1
biassz_y = darksz_y + self.pscan2 + self.oscan2
img_bias = np.zeros((biassz_y, biassz_x), dtype=int)
img_bias[
self.pscan2:self.pscan2+darksz_y,
self.pscan1:self.pscan1+darksz_x
] = np.random.poisson(image)
sub_img = image[100:-100, 100:-100]
print(np.random.poisson(sub_img).std(), sub_img.mean())
image = img_bias
lower_limit_n = int(np.log(emgain)/np.log(1 + 0.2)) # to make pEM < 0.2
emregister_num = np.maximum(lower_limit_n, self.emreg_cal_num)
emregister_num = np.minimum(emregister_num, self.emreg_num)
pEM = np.exp(np.log(emgain)/emregister_num) - 1
big_cte = self.em_cte ** (self.emreg_num / emregister_num)
image += np.random.poisson(image )
# for _ in range(emregister_num):
# # image += np.random.binomial(image, pEM)
# residual = np.random.binomial(image, 1-big_cte)
# image[:, 1:] += residual[:, :-1] - residual[:, 1:]
sub_img = image[100:-100, 100:-100] / emgain
print(sub_img.std(), sub_img.mean(), emgain)
bias = self.bias_frame()
image = image / self.ph_per_adu + bias
image = np.minimum(image, self.max_adu)
image = np.maximum(image, 0)
return image#.astype(dtype=np.uint16)
class EMCCD(object):
"""
EMCCD camera class
Parameters
----------
config_file : str
config file name
Attributes
----------
switch : dict
switch for each camera effects, including:
- 'flat': bool,
- 'dark': bool,
- 'stripe': bool,
- 'cic': bool,
- 'cte': bool,
- 'badcolumn': bool,
- 'nonlinear': bool,
- 'cosmicray': bool,
- 'blooming': bool,
"""
def __init__(self, config_file="ccd201_config.yaml"):
self.plszx = 1024
self.plszy = 1024
self.pscan1 = 8
self.pscan2 = 0
self.oscan1 = 16
self.oscan2 = 18
self.udark = 6
self.bdark = 2
self.ldark = 16
self.rdark = 16
self.fullwell = 80_000
self.em_fullwell = 780_000
# if config file exists, load it, otherwise use default values
config_file = cpism_refdata + '/camera/' + config_file
log.debug(f"emccd config file: {config_file}")
if os.path.exists(config_file):
self.load_config(config_file)
else:
raise(ValueError('config_file error'))
self.flat_shape = [self.plszy, self.plszx]
darksz_x = self.plszx + self.rdark + self.ldark
darksz_y = self.plszy + self.udark + self.bdark
self.dark_shape = [darksz_y, darksz_x]
biassz_x = darksz_x + self.pscan1 + self.oscan1
biassz_y = darksz_y + self.pscan2 + self.oscan2
self.image_shape = [biassz_y, biassz_x]
self.flat = fits.getdata(self.flat_file)
self.cic = fits.getdata(self.cic_file)
self.dark = fits.getdata(self.dark_file)
self.bad_col = fits.getdata(self.bad_col_file)
def load_config(self, config_file):
"""
load config file. Only for internal use.
"""
with open(config_file, 'r') as f:
config = yaml.load(f, Loader=yaml.FullLoader)
self.switch = config['switch']
self.readout_noise = config['readout_noise']
self.ph_per_adu = config['ph_per_adu']
self.bias_level = config['bias_level']
self.max_adu = config['max_adu']
self.dark_file = cpism_refdata + "/camera/" + config['dark_file']
self.flat_file = cpism_refdata + "/camera/" + config['flat_file']
self.cic_file = cpism_refdata + "/camera/" + config['cic_file']
self.bad_col_file = cpism_refdata + \
"/camera/" + config['bad_col_file']
def vertical_blooming(self, image):
"""
vertical blooming effect
"""
fullwell = self.fullwell
line = np.arange(image.shape[0])
yp, xp = np.where(image > fullwell)
n_saturated = len(xp)
log.debug(f"{len(xp)} pixels are saturated!")
if n_saturated > 5000:
log.warning(f"More than 5000({len(xp)}) pixels are saturated!")
img0 = image.copy()
for x, y in zip(xp, yp):
image[:, x] += np.exp(-(line-y)**2/20**2) * img0[y, x] * 0.2
return np.minimum(image, fullwell)
def nonlinear_effect(self, image):
"""
nonlinear effect
"""
fullwell = self.fullwell
nonlinear_coefficient = 0.1
log.debug(
f"nonlinear effect added with coefficient {nonlinear_coefficient}")
image += (image / fullwell)**2 * nonlinear_coefficient * fullwell
return image
def emregester_blooming(self, image, max_iteration=5):
"""
emregester blooming effect
"""
line = image.flatten().copy()
curve_x = np.arange(1300)+2
curve_y = np.exp(11*curve_x**(-0.19)-11)
curve_y[0] = 0
curve_y /= curve_y.sum()
over_limit_coe = 0.999
saturated = image > self.em_fullwell
n_saturated = saturated.sum()
if n_saturated > 0:
log.debug(f"{n_saturated} pixels are saturated during EM process.")
if n_saturated > 2000:
log.warning(
f"More than 2000 ({n_saturated}) pixels are saturated during EM process!")
for index in range(max_iteration):
over_limit = np.maximum(
line - self.em_fullwell * over_limit_coe, 0)
line = np.minimum(line, self.em_fullwell * over_limit_coe)
blooming = np.convolve(over_limit, curve_y, mode='full')[
:len(line)]
line = line + blooming
n_over = (line > self.em_fullwell).sum()
if n_over <= 0:
break
log.debug(
f'{index}/{max_iteration} loop: saturated pixel number: {n_over}')
return line.reshape(image.shape)
def cte(self, image):
"""
cte effect
"""
image = self.emregester_blooming(image, max_iteration=5)
return image
def readout(self, image_focal, emgain, expt, image_cosmic_ray=None):
"""
emccd readout. Generate a image with emccd readout effect.
Parameters
----------
image_focal : numpy.ndarray
image at focal plane. Unit: electron / second
emgain : float
emgain of emccd
expt : float
exposure time. Unit: second
image_cosmic_ray : numpy.ndarray, optional
cosmic ray image. Unit: electron, by default None
Returns
-------
numpy.ndarray
image with emccd readout effect. Unit: ADU
Notes
-----
1. effects include: dark, flat, cte, blooming, nonlinear, etc. Can be turned on/off by switch.
2. size of input image_focal must be 1024x1024
3. size of output image is 1080x1056 (including overscan and dark reference region)
4. Q.E is not included in this function. It should be included in image_focal. See optics.py for details.
"""
log.debug(
fr"EMCCD readout: {emgain=:}, {expt=:}, image_comic_ray:{'None' if image_cosmic_ray is None else 'Not None'}")
log.debug(
f"camera effects switch={self.switch}"
)
image = image_focal
if self.switch['flat']:
image = image * self.flat
if self.switch['nonlinear']:
image = self.nonlinear_effect(image)
darksz_x = self.plszx + self.rdark + self.ldark
darksz_y = self.plszy + self.udark + self.bdark
img_dark = np.zeros((darksz_y, darksz_x))
img_dark[
self.bdark:self.plszy+self.bdark,
self.ldark:self.ldark+self.plszx
] = image
image = img_dark
if self.switch['dark']:
image += self.dark
image *= expt
if self.switch['cic']:
image += self.cic
if image_cosmic_ray is not None and self.switch['cosmicray']:
image += image_cosmic_ray
if self.switch['blooming']:
image = self.vertical_blooming(image)
if self.switch['badcolumn']:
for i in range(self.bad_col.shape[1]):
deadpix_x = self.bad_col[0, i]
deadpix_y = self.bad_col[1, i]
image[deadpix_y:, deadpix_x] = 0
biassz_x = darksz_x + self.pscan1 + self.oscan1
biassz_y = darksz_y + self.pscan2 + self.oscan2
img_bias = np.zeros((biassz_y, biassz_x), dtype=int)
seed = random_seed_select()
log.debug(f"photon noise seed: {seed}")
print(image.mean())
img_bias[
self.pscan2:self.pscan2+darksz_y,
self.pscan1:self.pscan1+darksz_x
] = np.random.poisson(image)
image = img_bias
if self.switch['cte']:
image = self.cte(image * emgain) / emgain
seed = random_seed_select()
log.debug(f"gamma noise seed: {seed}")
if emgain != 1:
image = np.random.gamma(image, emgain)
image = np.minimum(image, self.em_fullwell)
seed = random_seed_select()
log.debug(f"readout noise seed: {seed}")
image += np.random.randn(biassz_y, biassz_x) * self.readout_noise
image = image / self.ph_per_adu + self.bias_level
if self.switch['stripe']:
image += self.add_stripe_effect(image)
image = np.minimum(image, self.max_adu)
image = np.maximum(image, 0)
return image.astype(np.uint16)
def add_stripe_effect(self, image):
"""
add stripe effect
"""
shape = image.shape
v_width = 1
v_amplitude = 30
v_limit = 0.01
v_base = 10
h_width = 20
h_amplitude = 10
h_limit = 0.9
h_base = 20
index = np.linspace(0, np.pi, shape[0] * shape[1])
def stripe(width, limit, amplitude, base, axis=0):
seed = random_seed_select()
log.debug(f"stripe noise seed: {seed}")
dim_axis = shape[axis]
dim_other = shape[0] * shape[1] // shape[axis]
value = np.sin(index / width * dim_axis + np.pi *
dim_axis / width * np.random.randint(1024))
value = np.maximum(value, -limit)
value = np.minimum(value, limit)
value = (value / limit + limit) / 2 * amplitude + base
value = value.reshape(dim_axis, dim_other)
if axis == 1:
value = value.T
return value
output = stripe(v_width, v_limit, v_amplitude, v_base, axis=1)
output += stripe(h_width, h_limit, h_amplitude, h_base, axis=0)
return output
# # plt.plot(horizontal_index, horizontal_value)
# # # plt.xlim([0, 6.28])
# # plt.show()
# fits.writeto('horizontal_value.fits', output, overwrite=True)
# if __name__ == '__main__':
# import matplotlib.pyplot as plt
# emccd = EMCCD()
# image_focal = np.zeros((emccd.plszy, emccd.plszx)) + 1000
# image_focal[100:105, 100:105] = 10_000_000
# after_cte = emccd.emregester_blooming(image_focal, max_iteration=100)
# print(after_cte.sum(), image_focal.sum())
# fits.writeto('after_cte.fits', after_cte, overwrite=True)
# # darksz_x = emccd.plszx + emccd.rdark + emccd.ldark
# # darksz_y = emccd.plszy + emccd.udark + emccd.bdark
# # iamge_cosmic_ray = np.zeros((darksz_y, darksz_x))
# # emgain = 10
# # expt = 10
# # image = emccd.readout(image_focal, emgain, expt, iamge_cosmic_ray)
# # fits.writeto('test.fits', image, overwrite=True)
# image = np.zeros((1000, 1000))
# make_cosmic_ray_frame = CosmicRayFrameMaker()
# crimage = make_cosmic_ray_frame(image.shape, 3000)
# fits.writeto('crimage.fits', crimage, overwrite=True)
# # emccd.add_stripe_effect(image)
CpicImgSim/config.py deleted 100644 → 0
View file @ b7f8c4fa
import os, yaml
import warnings
cpism_refdata = os.environ.get('cpism_refdata', './cpism_refdata')
if not os.path.exists(cpism_refdata): # pragma: no cover
raise Exception(
"Can not find CPISM reference data.")
if not os.path.exists(os.environ.get('PYSYN_CDBS', './trd')): # pragma: no cover
raise Exception(
"Can not find PYSYN Stellar reference data.")
# we need to ignore the warning from pysynphot, because we only use the star models.
with warnings.catch_warnings(): # pragma: no cover
warnings.filterwarnings("ignore")
import pysynphot as S
solar_spectrum = S.FileSpectrum(
f"{os.environ['PYSYN_CDBS']}/grid/solsys/solar_spec.fits")
solar_spectrum.convert('photlam')
config_file = cpism_refdata + '/optics/optics_config.yaml'
if not os.path.exists(config_file): # pragma: no cover
raise FileNotFoundError(f"光学配置文件不存在({config_file})")
with open(config_file, 'r') as f:
optics_config = yaml.load(f, Loader=yaml.FullLoader)
MAG_SYSTEM = 'abmag'
__version__ = '2.0.0'
def which_focalplane(band):
"""
Return the name of the focalplane which the band belongs to.
Parameters
-----------
band: str
The name of the band.
from ['f565', 'f661', 'f743', 'f883', 'f940', 'f1265', 'f1425', 'f1542', 'wfs']
Returns
--------
str
The name of the focalplane.
'vis' or 'nir' or 'wfs'
Raises
-------
ValueError
If the band is not in ['f565', 'f661', 'f743', 'f883', 'f940', 'f1265', 'f1425', 'f1542', 'wfs']
"""
band = band.lower()
if band in ['f565', 'f661', 'f743', 'f883']:
return 'vis'
if band in ['f940', 'f1265', 'f1425', 'f1542']:
return 'nir'
if band in ['wfs']:
return 'wfs'
raise ValueError(f"未知的波段{band}")
\ No newline at end of file
CpicImgSim/io.py deleted 100644 → 0
View file @ b7f8c4fa
import numpy as np
from astropy.io import fits
import yaml
from CpicImgSim.config import cpism_refdata, __version__, which_focalplane
from CpicImgSim.utils import Logger
import os
from datetime import datetime, timedelta
from astropy.coordinates import SkyCoord
import re
import json
import pandas as pd
config_file = f'{cpism_refdata}/cpism_config.yaml'
with open(config_file, 'r') as fid:
config = yaml.load(fid, Loader=yaml.FullLoader)
output_dir = config['output_dir']
if config['relative_path']:
ref_dir_base = os.path.dirname(cpism_refdata)
output_dir = f'{ref_dir_base}/{output_dir}'
log_dir = output_dir + '/LOG'
tmp_dir = config['tmp_dir']
log_level = config['log_level']
header_check = config['check_fits_header']
for dir in ['', 'TMP', 'CAL', 'SCI', 'LOG']:
sub_dir = f"{output_dir}/{dir}"
if not os.path.exists(sub_dir):
os.makedirs(sub_dir)
tmp_folder_path = '.'
if tmp_dir == 'TMP':
tmp_folder_path = output_dir + '/TMP'
log = Logger(log_dir+'/cpism_pack.log', log_level).logger
def check_and_update_fits_header(header):
"""
Check the header keywords and update the description according to the data model.
Parameters
-----------
header: astropy.io.fits.header.Header
The header to be checked.
Returns
--------
None
"""
hdu = 'image'
if 'FILETYPE' in header.keys():
hdu = 'primary'
model_file = f"{cpism_refdata}/io/image_header.json"
if hdu == 'primary':
model_file = f"{cpism_refdata}/io/primary_header.json"
with open(model_file, 'r', encoding='utf-8') as fid:
data_model = json.load(fid)
dm_comment = {}
def print_warning(info):
if header_check:
log.warning(info)
# check existance and format of keyword in fits header
for keyword, comment, format, _, _ in data_model:
if pd.isnull(comment):
comment = ''
if len(comment) > 46:
comment = comment[:46]
print_warning(
f"Keyword {keyword} has a comment longer than 46 characters. It will be truncated to 46 characters.")
dm_comment[keyword] = comment
if keyword not in header.keys():
print_warning(f"Keyword {keyword} not found in [{hdu}] header.")
elif not pd.isnull(format):
value = header[keyword]
# check the type of the value, I for int, R for float, C for str
if isinstance(value, str):
type = 'C'
elif isinstance(value, float):
type = 'R'
elif isinstance(value, bool):
type = 'L'
elif isinstance(value, int):
type = 'I'
else:
type = 'U'
if type != format[0]:
print_warning(
f"Keyword {keyword} has wrong type in [{hdu}]. {format[0]} expected, {type} found.")
# check if there are extral keyword in header, and update the comment
for keyword in header.keys():
# print(keyword)
if keyword not in dm_comment.keys():
print_warning(
f"Keyword {keyword} not found in the [{hdu}] data model.")
else:
header[keyword] = (header[keyword], dm_comment[keyword])
return header
def obsid_parser(
obsid: int):
"""
Parse the obsid to get the obstype.
Parameters
----------
obsid: str
The obsid of the observation.
Obsid must be 11 digits and start with 5 for CPIC.
Returns
-------
str
The obstype of the observation.
Raises
------
ValueError
If the obsid is not 11 digits or does not start with 5.
"""
obsid = str(obsid)
if len(obsid) != 11:
raise ValueError('Obsid must be 11 digits.')
if obsid[0] != '5':
raise ValueError('Obsid must start with 5 for CPIC')
obstype_dict = {
'00': 'BIAS',
'01': 'DARK',
'02': 'FLAT',
'03': 'BKGD',
'04': 'LASR',
'10': 'SCIE',
'11': 'DENF',
'12': 'CALS',
'15': 'TEMP'
}
obstype = obstype_dict.get(obsid[1:3], 'DEFT')
return obstype
def datetime_obj_to_mjd(time_obj):
"""
transfer datetime object to mean julian date (MJD).
Parameters
----------
time_obj: datetime.datetime
The datetime object.
Returns
-------
float
The mean julian date (MJD).
"""
return (time_obj - datetime(1858, 11, 17)).total_seconds() / 86400
def primary_hdu(
obs_info: dict,
gnc_info: dict,
filename_output=False):
"""
Generate the primary hdu of the fits file.
Parameters
----------
obs_info: dict
The parameters of the observation. See `save_fits` function.
gnc_info: dict
The gnc information of the observation.
filename_output: bool (optional)
If True, the folder and the filename will be returned.
Returns
-------
fits.PrimaryHDU
The primary hdu of the fits file.
str, str
The folder and filename of the fits file. Only returned if filename_output is True.
Notes
-----
The gnc_info dict should contain the information of orbit and observation.
these informations are used to genrated the hdu header. Refer to the data model for more information.
"""
camera_config, _ = load_camera_and_optics_config(obs_info['band'])
obsid = obs_info['obsid']
exp_start = gnc_info.get(
'EXPSTART', datetime.now().isoformat(timespec='seconds'))
exp_start = datetime.fromisoformat(exp_start)
duartion = (obs_info['expt'] +
camera_config['readout_time']) * obs_info['nframe']
default_end = exp_start + timedelta(seconds=duartion)
exp_end = gnc_info.get('EXPEND', default_end.isoformat(timespec='seconds'))
exp_end = datetime.fromisoformat(exp_end)
filename = "CSST_CPIC"
filename += "_" + which_focalplane(obs_info['band']).upper()
filename += "_" + obsid_parser(obsid)
filename += "_" + exp_start.strftime("%Y%m%d%H%M%S")
filename += "_" + exp_end.strftime("%Y%m%d%H%M%S")
filename += f"_{obsid}_X_L0_V01.fits"
type_dir = 'CAL'
if str(obsid)[1] == '1':
type_dir = 'SCI'
mjd_dir = f"{datetime_obj_to_mjd(exp_start):.0f}"
folder = f"{type_dir}/{mjd_dir}"
header = fits.Header()
# General keywords
header['SIMPLE'] = True
header['BITPIX'] = 8
header['NAXIS'] = 0
header['EXTEND'] = True
header['NEXTEND'] = 1 # + parameters['nframe']
header['GROUPS'] = False
header['DATE'] = datetime.now().isoformat(timespec='seconds')
heaer_filename = filename[:-4]
if len(heaer_filename) > 68:
heaer_filename = heaer_filename[:68]
header['FILENAME'] = heaer_filename
header['FILETYPE'] = obsid_parser(obsid)
header['TELESCOP'] = 'CSST'
header['INSTRUME'] = 'CPIC'
header['RADECSYS'] = 'ICRS'
header['EQUINOX'] = 2000.0
header['FITSCREA'] = f'CPISM V{__version__}'
cstar = {'ra': '0d', 'dec': '0d'}
if obs_info['target'] != {}:
cstar = obs_info['target']['cstar']
radec = SkyCoord(cstar['ra'], cstar['dec'])
target_name = radec.to_string('hmsdms')
target_name = re.sub(R'[hdms\s]', '', target_name)
header['OBJECT'] = cstar.get('name', target_name)
header['TARGET'] = target_name
header['OBSID'] = str(obsid)
header['OBJ_RA'] = radec.ra.degree
header['OBJ_DEC'] = radec.dec.degree
# telescope information
header['REFFRAME'] = 'CSSTGSC-1.0'
header['DATE-OBS'] = exp_start.isoformat(timespec='seconds')
header['SATESWV'] = 'SIMULATION'
header['EXPSTART'] = datetime_obj_to_mjd(exp_start)
header['CABSTART'] = gnc_info.get('CABSTART', header['EXPSTART'])
header['SUNANGL0'] = gnc_info.get('SUNANGL0', -1.0)
header['MOONANG0'] = gnc_info.get('MOONANG0', -1.0)
header['TEL_ALT0'] = gnc_info.get('TEL_ALT0', -1.0)
header['POS_ANG0'] = gnc_info.get(
'POS_ANG0', float(obs_info['rotation']))
header['POSI0_X'] = gnc_info.get('POSI0_X', -1.0)
header['POSI0_Y'] = gnc_info.get('POSI0_Y', -1.0)
header['POSI0_Z'] = gnc_info.get('POSI0_Z', -1.0)
header['VELO0_X'] = gnc_info.get('VELO0_X', -1.0)
header['VELO0_Y'] = gnc_info.get('VELO0_Y', -1.0)
header['VELO0_Z'] = gnc_info.get('VELO0_Z', -1.0)
header['EULER0_1'] = gnc_info.get('EULER0_1', -1.0)
header['EULER0_2'] = gnc_info.get('EULER0_2', -1.0)
header['EULER0_3'] = gnc_info.get('EULER0_3', -1.0)
header['RA_PNT0'] = gnc_info.get('RA_PNT0', header['OBJ_RA'])
header['DEC_PNT0'] = gnc_info.get('DEC_PNT0', header['OBJ_DEC'])
header['EXPEND'] = datetime_obj_to_mjd(exp_end)
header['CABEND'] = gnc_info.get('CABEDN', header['EXPEND'])
header['SUNANGL1'] = gnc_info.get('SUNANGL1', header['SUNANGL0'])
header['MOONANG1'] = gnc_info.get('MOONANG1', header['MOONANG0'])
header['TEL_ALT1'] = gnc_info.get('TEL_ALT1', header['TEL_ALT0'])
header['POS_ANG1'] = gnc_info.get('POS_ANG1', header['POS_ANG0'])
header['POSI1_X'] = gnc_info.get('POSI1_X', header['POSI0_x'])
header['POSI1_Y'] = gnc_info.get('POSI1_Y', header['POSI0_y'])
header['POSI1_Z'] = gnc_info.get('POSI1_Z', header['POSI0_z'])
header['VELO1_X'] = gnc_info.get('VELO1_X', header['VELO0_x'])
header['VELO1_Y'] = gnc_info.get('VELO1_Y', header['VELO0_y'])
header['VELO1_Z'] = gnc_info.get('VELO1_Z', header['VELO0_z'])
header['EULER1_1'] = gnc_info.get('EULER1_1', header['EULER0_1'])
header['EULER1_2'] = gnc_info.get('EULER1_2', header['EULER0_2'])
header['EULER1_3'] = gnc_info.get('EULER1_3', header['EULER0_3'])
header['RA_PNT1'] = gnc_info.get('RA_PNT1', header['RA_PNT0'])
header['DEC_PNT1'] = gnc_info.get('DEC_PNT1', header['DEC_PNT0'])
header['EXPTIME'] = (exp_end - exp_start).total_seconds()
header['EPOCH'] = float(exp_start.year)
header['CHECKSUM'] = '0000000000000000'
header['DATASUM'] = '0000000000'
check_and_update_fits_header(header)
# other information
hdu = fits.PrimaryHDU(header=header)
hdu.add_checksum()
if filename_output:
return hdu, folder, filename
else:
return hdu
def load_camera_and_optics_config(band):
"""
Load camera and optics configuration from reference data.
Parameters
----------
band : str
Band name.
Returns camera_config, optics_config : dict, dict
"""
camera = which_focalplane(band)
if camera == 'vis':
config_file = 'emccd_config.yaml'
elif camera == 'nir':
raise ValueError('NIR camera is not supported yet')
config_file = 'nir_config.yaml'
with open(f"{cpism_refdata}/camera/{config_file}", 'r') as fid:
camera_config = yaml.load(fid, Loader=yaml.FullLoader)
with open(f"{cpism_refdata}/optics/optics_config.yaml", 'r') as fid:
optics_config = yaml.load(fid, Loader=yaml.FullLoader)[camera]
return camera_config, optics_config
def frame_header(obs_info, index, bunch_start, primary_header={}):
"""
Generate the header for a single frame.
Parameters
----------
obs_info : dict
Dictionary of parameters. See `save_fits` function.
index : int
Frame index.
bunch_start : str
Start time of the bunch.
primary_header : dict (optional)
Primary header. default: {}
Returns
-------
astropy.io.fits.Header
"""
header = fits.Header()
camera_config, optics_config = load_camera_and_optics_config(
obs_info['band'])
plszx = camera_config['plszx']
plszy = camera_config['plszy']
pscan1 = camera_config['pscan1']
pscan2 = camera_config['pscan2']
oscan1 = camera_config['oscan1']
oscan2 = camera_config['oscan2']
udark = camera_config['udark']
bdark = camera_config['bdark']
ldark = camera_config['ldark']
rdark = camera_config['rdark']
imgszx = plszx + pscan1 + oscan1 + ldark + rdark
imgszy = plszy + pscan2 + oscan2 + udark + bdark
header['XTENSION'] = 'IMAGE'
header['BITPIX'] = 16
header['NAXIS'] = 2
header['NAXIS1'] = 1080
header['NAXIS2'] = 1050
header['EXTNAME'] = 'IMAGE'
header['EXTVER'] = 1
header['BSCALE'] = 1.0
header['BZERO'] = 32768.0
header['BUNIT'] = 'ADU'
header['FILTER'] = obs_info['band']
header['DETSN'] = '00000000000'
header['DETNAME'] = camera_config['detector_name']
header['CHIPLAB'] = camera_config['ccd_label']
header['CHIPTEMP'] = float(camera_config['chip_temp'])
header['DEWTEMP'] = float(camera_config['dewar_temp'])
header['DETSIZE'] = f"{imgszx} * {imgszy}"
header['IMGINDEX'] = index
frame_time = obs_info['expt'] + camera_config['readout_time']
bunch_start = datetime.fromisoformat(bunch_start)
expstart = bunch_start + timedelta(seconds=frame_time * (index - 1))
bunch_start_mjd = datetime_obj_to_mjd(bunch_start)
ra0 = primary_header.get('RA_PNT0', -1.0)
dec0 = primary_header.get('DEC_PNT0', -1.0)
pa0 = primary_header.get('POS_ANG0', -1.0)
cab0 = primary_header.get('CABSTART', bunch_start_mjd)
delta_t = frame_time * index
bunch_end = bunch_start + timedelta(seconds=delta_t)
bunch_end_mjd = datetime_obj_to_mjd(bunch_end)
ra1 = primary_header.get('RA_PNT1', ra0)
dec1 = primary_header.get('DEC_PNT1', dec0)
pa1 = primary_header.get('POS_ANG1', pa0)
cab1 = primary_header.get('CABEND', bunch_end_mjd)
img_cab = datetime_obj_to_mjd(expstart)
ratio = (img_cab - cab0)/(cab1 - cab0)
ra = ra0 + (ra1 - ra0) * ratio
dec = dec0 + (dec1 - dec0) * ratio
pa = pa0 + (pa1 - pa0) * ratio
header['IMG_EXPT'] = expstart.isoformat(timespec='seconds')
header['IMG_CABT'] = header['IMG_EXPT']
header['IMG_DUR'] = float(obs_info['expt'])
header['IMG_PA'] = ra
header['IMG_RA'] = dec
header['IMG_DEC'] = pa
header['DATASECT'] = f"{plszx} * {plszy}"
header['PIXSCAL'] = optics_config['platescale']
header['PIXSIZE'] = camera_config['pitch_size']
header['NCHAN'] = 1
header['PSCAN1'] = pscan1
header['PSCAN2'] = pscan2
header['OSCAN1'] = oscan1
header['OSCAN2'] = oscan2
header['UDARK'] = udark
header['BDARK'] = bdark
header['LDARK'] = ldark
header['RDARK'] = rdark
# WCS
cstar = {'ra': '0d', 'dec': '0d'}
if obs_info['target'] != {}:
cstar = obs_info['target']['cstar']
radec = SkyCoord(cstar['ra'], cstar['dec'])
shift = obs_info['shift']
platescale = optics_config['platescale']
rotation = np.radians(obs_info['rotation'])
header['WCSAXES'] = 2
header['CRPIX1'] = (plszx + 1)/2 + pscan1 + ldark + shift[0] / platescale
header['CRPIX2'] = (plszy + 1)/2 + pscan2 + udark + shift[0] / platescale
header['CRVAL1'] = radec.ra.degree
header['CRVAL2'] = radec.dec.degree
header['CTYPE1'] = 'RA---TAN'
header['CTYPE2'] = 'DEC--TAN'
header['CD1_1'] = np.cos(rotation)
header['CD1_2'] = -np.sin(rotation)
header['CD2_1'] = np.sin(rotation)
header['CD2_2'] = np.cos(rotation)
header['others'] = 'other'
# Readout information
header['EMGAIN'] = float(obs_info['emgain'])
header['GAIN'] = float(camera_config['ph_per_adu'])
header['DET_BIAS'] = float(camera_config['bias_level'])
header['RON'] = float(camera_config['readout_noise'])
header['READTIME'] = float(camera_config['readout_time'])
header['ROSPEED'] = float(camera_config['readout_speed'])
# CPIC information
header['LS_STAT'] = 'OFF'
header['IWA'] = optics_config['mask_width'] / 2
header['WFSINFO1'] = -1.0
header['WFSINFO2'] = -1.0
header['CHECKSUM'] = '0000000000000000'
header['DATASUM'] = '0000000000'
header = check_and_update_fits_header(header)
return header
def save_fits_simple(images, obs_info):
"""
Save the image to a fits file with a simple header to TMP directory.
Parameters
----------
images : numpy.ndarray
Image array to be written.
obs_info : dict
Dictionary of observation informations. See `save_fits` function.
Returns
----------
Filename of the saved fits file.
"""
target = obs_info['target']
target_info = 'NO_TARGET'
if 'cstar' in target.keys():
target_info = ''
target_info = f"S{target['cstar']['magnitude']:.1f}"
target_info += f"_P{len(target.get('planets', '[]'))}"
name = target_info
if 'name' in target.keys():
name = target['name']
name = name.replace('/', '_')
name = name.replace(',', '_')
now = datetime.now()
time = now.strftime("%Y%m%d%H%M%S")
filename = f"{name}_{time}.fits"
header = fits.Header()
header['skybg'] = obs_info['skybg']
header['name'] = name
header['exptime'] = obs_info['expt']
header['nframe'] = obs_info['nframe']
header['band'] = obs_info['band']
header['emgain'] = obs_info['emgain']
header['obsid'] = obs_info['obsid']
header['rotation'] = obs_info['rotation']
shift = obs_info['shift']
header['shift'] = f"x:{shift[0]},y:{shift[1]}"
fullname = f"{tmp_folder_path}/{filename}"
fits.writeto(fullname, images, overwrite=True, header=header)
return fullname
def save_fits(images, obs_info, gnc_info, csst_format=True):
"""
Save the image to a fits file.
Parameters
----------
images : numpy.ndarray
Image array to be saved.
obs_info : dict
Dictionary of observation informations.
Must contain the following keys
- band: str
- Band of the image.
- expt: float
- Exposure time of the each image.
- nframe: int
- Number of frames in the image.
- emgain: int
- EM gain of the camera.
- obsid: str
- Observation ID. Obsid must be 11 digits and start with 5 for CPIC. See pharse_obsid() for details.
- rotation: float
- Rotation angle of the image.
- shift: list
- Shift of the image.
gnc_info : dict
Dictionary of GNC information.
Contains the keywords in the primary header. See primary_hdu() for details.
csst_format : bool, optional
Whether to save the fits file in CSST format, by default True.
"""
if not csst_format:
save_fits_simple(images, obs_info)
return
hdu_header, folder, filename = primary_hdu(obs_info, gnc_info, True)
hdu_list = fits.HDUList([hdu_header])
if len(images.shape) == 2:
images = np.array([images])
for index in range(images.shape[0]):
header = frame_header(
obs_info,
index + 1,
hdu_header.header['DATE-OBS'],
primary_header=hdu_header.header
)
frame_hdu = fits.ImageHDU(images[index, :, :], header=header)
frame_hdu.add_checksum()
hdu_list.append(frame_hdu)
folder = f"{output_dir}/{folder}"
if not os.path.exists(folder):
os.makedirs(folder)
hdu_list.writeto(f"{folder}/{filename}", overwrite=True)
CpicImgSim/main.py deleted 100644 → 0
View file @ b7f8c4fa
import numpy as np
import re
from .target import spectrum_generator
from .optics import make_focus_image, focal_mask, optics_config
from .psf_simulation import simulate_psf
from .camera import EMCCD, CosmicRayFrameMaker, sky_frame_maker
from .io import save_fits, log
from .config import which_focalplane
def psf_function(band, cstar_spectrum, shape, error=0.1):
cstar = True
if shape < 300:
cstar = False
return simulate_psf(error, band, cstar_spectrum, nsample=1, cstar=cstar)
def observation_simulation(
target: dict,
skybg: float,
expt: float,
nframe: int,
band: str,
emgain: float,
obsid: int = 51900000000,
rotation: float = 0,
shift: list = [0, 0],
gnc_info: dict = {},
csst_format: bool = True,
psf_function: callable = psf_function):
"""
Simulate the observation. All-In-One function of the package.
Parameters
-----------
target: dict
The target information. See target.py for details.
skybg: float
magnitude of the skybackground at the input b and. (abmag system)
expt: float
exposure time in second.
nframe: int
number of frames to be simulated.
band: str
the band of the observation. (e.g. 'f661')
emgain: float
the EM gain of the camera.
obsid: int
the observation id. Default is 51900000000.
rotation: float
the rotation angle of the target in degree. Default is 0.
shift: list
the shift of the target in arcsec. Default is [0, 0].
gnc_info: dict
the gnc information. Default is {}. See io.py for details.
csst_format: bool
whether to save the fits file in CSST format. Default is True.
psf_function: callable
the function to generate the psf. See optics.py for details.
Returns
-----------
np.ndarray of the simulated images with shape (nframe, 1088, 1050).
"""
target_list = []
if 'cstar' in target.keys():
target_list = spectrum_generator(target)
focal_name = which_focalplane(band)
this_focal_config = optics_config[focal_name]
telescope_config = optics_config['telescope']
area = telescope_config['aperature_area']
if focal_name == 'vis':
camera = EMCCD()
else:
raise ValueError('Only VIS focal plane is supported.')
platescale = this_focal_config['platescale']
iwa = this_focal_config['mask_width'] / 2
crmaker = CosmicRayFrameMaker()
images = []
params = {
'target': target,
'skybg': skybg,
'expt': expt,
'nframe': nframe,
'band': band,
'emgain': emgain,
'obsid': obsid,
'rotation': rotation,
'shift': shift,
}
paramstr = ', '.join([f'{k}={v}' for k, v in params.items()])
log.debug(f"parameters: {paramstr}")
for i in range(nframe):
log.info(f'Simulation Running: Frame {i+1}/{nframe}')
focal_frame = make_focus_image(
band,
target_list,
psf_function,
platesize=camera.flat_shape,
rotation=rotation,
init_shifts=shift,
)
if skybg is None or skybg > 100:
sky_bkg_frame = 0
else:
sky_bkg_frame = sky_frame_maker(
band,
skybg,
platescale,
camera.flat_shape
)
focal_frame = (focal_frame + sky_bkg_frame) * area
focal_frame = focal_mask(focal_frame, iwa, platescale)
cr_frame = crmaker.make_cr_frame(camera.dark_shape, expt)
img = camera.readout(
focal_frame,
emgain,
expt,
image_cosmic_ray=cr_frame
)
images.append(img)
images = np.array(images)
save_fits(images, params, gnc_info, csst_format=csst_format)
return images
def quick_run(
target_str: str,
skymag: float,
band: str,
expt: float,
nframe: int,
emgain: float,
rotation: float = 0,
shift: list = [0, 0]) -> np.ndarray:
"""
A quick run function to simulate the observation.
Parameters
-----------
target_str: str
The target information in string format.
In the format of "\*5.1/25.3(1.3,1.5)/22.1(2.3,-4.5)" which means a central star
with magnitude 5.1, and two substellar with magnitude 25.3 and 22.1, respectively.
The first number in the parenthesis is the x position in arcsec, and the second is the y position.
skybg: float
magnitude of the skybackground at the input band. (abmag system)
band: str
the band of the observation. (e.g. 'f661')
expt: float
exposure time in second.
nframe: int
number of frames to be simulated.
emgain: float
the EM gain of the camera.
rotation: float (optional)
the rotation angle of the target in degree. Default is 0.
shift: list (optional)
the shift of the target in arcsec. Default is [0, 0].
Returns
-----------
np.ndarray of the simulated images, with shape (nframe, 1088, 1050)
Notes
-----------
1. stars are assumed to be G0III with distance 10pc.
2. magnitude of the star and substellar are assumed to be in the same band.
3. Csst format is not supported.
4. The psf is assumed to be the default one.
5. fits file will be saved in the current directory.
"""
log.info(f'Quick Run: {target_str}')
target_dict = {
'name': 'cal',
}
if (target_str != '') and (target_str[0] == '*'):
objects = target_str[1:].split('/')
cstar_mag = float(objects[0])
cstar = {
'magnitude': cstar_mag,
'ra': '0d',
'dec': '0d',
'sptype': 'G0III',
'distance': 10,
'mag_input_band': band
}
stars = []
for sub_stellar in objects[1:]:
float_regex = R"[+-]?\d+(?:\.\d+)?"
match = re.match(
rf"({float_regex})\(({float_regex}),({float_regex})\)", sub_stellar)
if not match:
raise ValueError('Wrong format for sub stellar.')
mag = float(match.group(1))
x = float(match.group(2))
y = float(match.group(3))
pangle = np.arctan2(x, y) * 180 / np.pi
separation = np.sqrt(x**2 + y**2)
stars.append({
'magnitude': mag,
'pangle': pangle,
'separation': separation,
'sptype': 'G0III',
'mag_input_band': band
})
target_dict = {
'name': target_str[1:],
'cstar': cstar,
'stars': stars,
}
return observation_simulation(
target=target_dict,
skybg=skymag,
expt=expt,
nframe=nframe,
band=band,
emgain=emgain,
csst_format=False,
shift=shift,
rotation=rotation,
)
# observation_simulation(
# target={},
# skybg=15,
# expt=10,
# nframe=2,
# band='f661',
# emgain=30,
# obsid=50112345678,
# )
# quick_run('*5.1/25.3(0.8,0.8)', None, 'f661', 10, 1, 10)
# quick_run('*5/20(0.8,0.8)', None, 'f883', 10, 1, 10)
# # quick *5.1/25.3(1.3,1.5) expt nframe emgain band rotation shift
# # quick target_name expt nframe emgain band rotation shift
# # plan plan_file_or_folder
if __name__ == '__main__': # pragma: no cover
target_example = {
'cstar': {
'magnitude': 1,
'ra': '120d',
'dec': '40d',
'distance': 10,
'sptype': 'F0III',
},
'stars': [
{
'magnitude': 20,
'pangle': 60,
'separation': 1,
'sptype': 'F0III'
}
]
}
# quick_run('', 10, 'f661', 1, 1, 30)
# quick_run('*2.4/10(3,5)/15(-4,2)', 21, 'f661', 1, 1, 30)
# # normal target
observation_simulation(
target=target_example,
skybg=21,
expt=1,
nframe=2,
band='f661',
emgain=30,
obsid=51012345678,
)
# # bias
# observation_simulation(
# target=target_example,
# skybg=999,
# expt=1,
# nframe=2,
# band='f661',
# emgain=1,
# obsid=51012345678,
# shift=[3, 3],
# rotation=60
# )
# # bias-gain
# observation_simulation(
# target={},
# skybg=999,
# expt=0.01,
# nframe=2,
# band='f661',
# emgain=1000,
# obsid=50012345678,
# )
# # dark
# observation_simulation(
# target={},
# skybg=999,
# expt=100,
# nframe=2,
# band='f661',
# emgain=30,
# obsid=50112345678,
# )
# # flat
# observation_simulation(
# target={},
# skybg=15,
# expt=10,
# nframe=2,
# band='f661',
# emgain=30,
# obsid=50112345678,
# )
CpicImgSim/optics.py deleted 100644 → 0
View file @ b7f8c4fa
import os
import yaml
import time
import scipy as sp
import numpy as np
from CpicImgSim.config import cpism_refdata, which_focalplane, S # S is synphot
from CpicImgSim.config import optics_config
from CpicImgSim.utils import region_replace
from CpicImgSim.io import log
from astropy.convolution import convolve_fft
from scipy.signal import fftconvolve
FILTERS = {
'f565': S.FileBandpass(f'{cpism_refdata}/throughtput/f565_total.fits'),
'f661': S.FileBandpass(f'{cpism_refdata}/throughtput/f661_total.fits'),
'f743': S.FileBandpass(f'{cpism_refdata}/throughtput/f743_total.fits'),
'f883': S.FileBandpass(f'{cpism_refdata}/throughtput/f883_total.fits'),
'f940': S.FileBandpass(f'{cpism_refdata}/throughtput/f940_total.fits'),
'f1265': S.FileBandpass(f'{cpism_refdata}/throughtput/f1265_total.fits'),
'f1425': S.FileBandpass(f'{cpism_refdata}/throughtput/f1425_total.fits'),
'f1542': S.FileBandpass(f'{cpism_refdata}/throughtput/f1542_total.fits'),
}
def filter_throughput(filter_name):
"""
Totally throughput of the each CPIC band.
Including the throughput of the filter, telescope, cpic, and camera QE.
If the filter_name is not supported, return the throughput of the default filter(f661).
Parameters
-----------
filter_name: str
The name of the filter.
One of ['f565', 'f661'(default), 'f743', 'f883', 'f940', 'f1265', 'f1425', 'f1542']
Returns
--------
synphot.Bandpass
The throughput of the filter.
"""
filter_name = filter_name.lower()
filter_name = 'f661' if filter_name == 'default' else filter_name
if filter_name not in FILTERS.keys():
log.warning(f"滤光片名称错误({filter_name})，返回默认滤光片(f661)透过率")
filter_name = 'f661'
return FILTERS[filter_name]
def example_psf_func(band, spectrum, frame_size, error=0.1):
"""
Example psf generating function.
Parameters
-------------
band: str
The name of the band.
spectrum: synphot.Spectrum or synphot.SourceSpectrum
The spectrum of the target.
frame_size: int
The size of the frame.
error: float
Phase RMS error.
Returns
---------------
2D array
psf image with shape of `frame_size`
"""
pass
def example_monochromatic_psf(wavelength, error=0.1):
pass
def rotate_and_shift(shift, rotation, init_shifts):
rotation_rad = rotation / 180 * np.pi
return np.array([
shift[0] * np.cos(rotation_rad) + shift[1] * np.sin(rotation_rad),
-shift[0] * np.sin(rotation_rad) + shift[1] * np.cos(rotation_rad)
]) + np.array(init_shifts)
from scipy.ndimage import rotate
def ideal_focus_image(
bandpass: S.spectrum.SpectralElement,
targets: list,
platescale,
platesize: list = [1024, 1024],
init_shifts: list = [0, 0],
rotation: float = 0,):
focal_image = np.zeros(platesize)
focal_shape = np.array(platesize)[::-1] # x, y
if not targets:
return focal_image
for target in targets:
sub_x, sub_y, sub_spectrum, sub_image = target
sub_shift = rotate_and_shift([sub_x, sub_y], rotation, init_shifts) / platescale
sed = (sub_spectrum * bandpass).integrate()
if sub_image is None:
x = (focal_shape[0] - 1)/2 + sub_shift[0]
y = (focal_shape[1] - 1)/2 + sub_shift[1]
int_x = int(x)
int_y = int(y)
if int_x < 0 or int_x >= focal_shape[0] - 1 or int_y < 0 or int_y >= focal_shape[1] - 1:
continue
dx1 = x - int_x
dx0 = 1 - dx1
dy1 = y - int_y
dy0 = 1 - dy1
sub = np.array([
[dx0*dy0, dx1*dy0],
[dx0*dy1, dx1*dy1]]) * sed
focal_image[int_y: int_y+2, int_x: int_x+2] += sub
else:
# sub_image = sub_image
sub_image = np.abs(rotate(sub_image, rotation, reshape=False))
sub_image = sub_image / sub_image.sum()
sub_img_shape = np.array(sub_image.shape)[::-1]
sub_shift += (focal_shape-1)/2 - (sub_img_shape-1)/2
focal_image = region_replace(
focal_image,
sub_image * sed,
sub_shift,
subpix=True
)
return focal_image
from scipy.signal import fftconvolve
def sp_convole_fft(image, kernal):
kernal = kernal / kernal.sum()
# y0 = kernal.shape[0] // 2
# x0 = kernal.shape[1] // 2
outimg = fftconvolve(image, kernal, mode='same')
# return outimg[y0:y0+image.shape[0], x0:x0+image.shape[1]]
return outimg
def convolve_psf(
band: str,
targets: list,
psf_function: callable,
init_shifts: list = [0, 0],
rotation: float = 0,
nsample: int = 5,
error: float = 1,
platesize: list = [1024, 1024]) -> np.ndarray :
config = optics_config[which_focalplane(band)]
platescale = config['platescale']
filter = filter_throughput(band)
wave = filter.wave
throughput = filter.throughput
min_wave = wave[0]
max_wave = wave[-1]
all_fp_image = []
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
center_wavelength = (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_fp_image = ideal_focus_image(i_band, targets, platescale, platesize, init_shifts, rotation)
psf = psf_function(center_wavelength, error=error)
t0 = time.time()
# c_fp_image = convolve_fft(i_fp_image, psf, allow_huge=True)
c_fp_image = sp_convole_fft(i_fp_image, psf)
print(f"Convolution time: {time.time()-t0}")
all_fp_image.append(c_fp_image)
return np.array(all_fp_image).mean(axis=0)
def make_focus_image(
band: str,
targets: list,
psf_function: callable,
init_shifts: list = [0, 0],
rotation: float = 0,
platesize: list = [1024, 1024]) -> np.ndarray:
"""
Make the focus image of the targets.
Parameters
-----------
band: str
The name of the band.
targets: list
The list of the targets.
Each element of the list is a tuple of (x, y, spectrum).
- x, y: float
- The position of the target in the focal plane.
- spectrum: synphot.Spectrum or synphot.SourceSpectrum
- The spectrum of the target.
psf_function: callable
The function to generate the PSF, with same parameters and return as `example_psf_func`.
init_shifts: list
The initial shifts of the center targets. Unit: arcsec.
The default is [0, 0].
rotation: float
The rotation of the focal plane. Unit: degree.
The default is 0 degree.
platesize: list
The size of the focal plane. Unit: pixel.
The default is [1024, 1024].
Returns
--------
np.ndarray
The focus image of the targets.
2D array with the shape of platesize.
"""
config = optics_config[which_focalplane(band)]
platescale = config['platescale']
focal_image = np.zeros(platesize)
if not targets:
return focal_image
cstar_x, cstar_y, cstar_spectrum = targets[0]
cstar_shift = rotate_and_shift([cstar_x, cstar_y]) / platescale
error_value = 0 # nm
cstar_psf = psf_function(band, cstar_spectrum, config['cstar_frame_size'],
error=error_value)
platesize = np.array(platesize)[::-1]
psf_shape = np.array(cstar_psf.shape)[::-1]
cstar_shift += (platesize-1)/2 - (psf_shape-1)/2
focal_image = region_replace(
focal_image,
cstar_psf,
cstar_shift,
padded_in=False,
padded_out=False,
subpix=True)
for i_target in range(1, len(targets)):
sub_x, sub_y, sub_spectrum = targets[i_target]
pdout = False if i_target == len(targets)-1 else True
pdin = False if i_target == 1 else True
log.debug(f"input target {sub_x=:}, {sub_y=:}")
sub_shift = rotate_and_shift([sub_x, sub_y], rotation, init_shifts) / platescale
log.debug(f"after rotate and shift {sub_shift=:}")
sub_psf = psf_function(
band,
sub_spectrum,
config['substellar_frame_size'],
error=error_value
)
psf_shape = np.array(sub_psf.shape)[::-1]
sub_shift += (platesize-1)/2 - (psf_shape-1)/2
log.debug(f"input shift of region_replace: {sub_shift=:}")
focal_image = region_replace(
focal_image,
sub_psf,
sub_shift,
padded_in=pdin,
padded_out=pdout,
subpix=True
)
return focal_image
def focal_mask(image, iwa, platescale, throughtput=1e-6):
"""
Mask the image outside the inner working angle.
Parameters
-----------
image: np.ndarray
The image to be masked.
iwa: float
The inner working angle. Unit: arcsec.
platescale: float
The platescale of the image. Unit: arcsec/pixel.
throughtput: float
The throughtput of the mask. The default is 1e-6.
Returns
--------
np.ndarray
The masked image.
"""
xx, yy = np.mgrid[0:image.shape[0], 0:image.shape[1]]
center = np.array([(image.shape[0]-1)/2, (image.shape[1]-1)/2])
mask = (abs(xx - center[0]) < iwa /
platescale) | (abs(yy - center[1]) < iwa / platescale)
image_out = image.copy()
image_out[mask] *= throughtput
return image_out
CpicImgSim/psf_simulation.py deleted 100644 → 0
View file @ b7f8c4fa
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)
CpicImgSim/target.py deleted 100644 → 0
View file @ b7f8c4fa
import os
import re
import json
import numpy as np
from scipy import constants
from astropy.io import fits
from astropy.coordinates import SkyCoord
from CpicImgSim.config import cpism_refdata, MAG_SYSTEM
from CpicImgSim.config import S # pysynphot
from CpicImgSim.optics import filter_throughput
from CpicImgSim.io import log
PLANK_H = constants.h # ~6.62607015e-34
LIGHT_SPEED = constants.c # ~3e8
DEG2RAD = np.pi / 180
R_JUPITER_METER = 6.99e4
AU_METER = 1.49e8
DEFAULT_FILTER_FILE = f'{cpism_refdata}\\throughtput\\f661_total.fits'
HYBRID_MODEL_FILE = f'{cpism_refdata}\\target_model\\hybrid_model.fits'
BCC_MODEL_FOLDER = f'{cpism_refdata}\\target_model\\bccmodels'
CATALOG_CACHE = {}
class AlbedoCat(S.spectrum.TabularSpectralElement, S.catalog.Icat):
"""
This class is used to generate a spectrum from the planet reflection model.
Parameters
----------
phase : float
[degree], The phase angle of the planet, range from 0 to 180.
metallicity : float
The metallicity of the planet. log(Z/Zsun), range from 0 to 2.
f_sed : float
The sedimentation efficiency of the cloud. log(f_sed), range from -2 to 2. 2 is for cloud free case.
Notes
-----
f_sed is only reliable between -2 and log10(6) and the cloud free case (2).
Values between log10(6) and 2 are interpolated from the cloud free case (2) and log10(6).
"""
def __init__(self,
phase: float,
metallicity: float,
f_sed: float,
):
catname = 'BCCalbedo'
self.isAnalytic=False
self.name = f"{catname}_{phase}_{metallicity}_{f_sed}"
self.parameter_names = ['phase', 'metallicity', 'f_sed']
self.catalog_folder = BCC_MODEL_FOLDER
filename = os.path.join(self.catalog_folder, 'catalog.fits')
if filename in CATALOG_CACHE:
indices = CATALOG_CACHE[filename]
else:
with fits.open(filename) as table:
indexList = table[1].data.field('INDEX')
filenameList = table[1].data.field('FILENAME')
indices = self._getArgs(indexList, filenameList)
CATALOG_CACHE[filename] = indices
list0,list1 = self._breakList(indices, 0, phase)
list2,list3 = self._breakList(list0, 1, metallicity)
list4,list5 = self._breakList(list1, 1, metallicity)
list6,list7 = self._breakList(list2, 2, f_sed)
list8,list9 = self._breakList(list3, 2, f_sed)
list10,list11 = self._breakList(list4, 2, f_sed)
list12,list13 = self._breakList(list5, 2, f_sed)
sp1, wave, waveunits = self._getSpectrum(list6[0], catname, wave_output=True)
sp2 = self._getSpectrum(list7[0], catname)
sp3 = self._getSpectrum(list8[0], catname)
sp4 = self._getSpectrum(list9[0], catname)
sp5 = self._getSpectrum(list10[0], catname)
sp6 = self._getSpectrum(list11[0], catname)
sp7 = self._getSpectrum(list12[0], catname)
sp8 = self._getSpectrum(list13[0], catname)
spa1 = self._interpolateSpectrum(sp1, sp2, f_sed)
spa2 = self._interpolateSpectrum(sp3, sp4, f_sed)
spa3 = self._interpolateSpectrum(sp5, sp6, f_sed)
spa4 = self._interpolateSpectrum(sp7, sp8, f_sed)
spa5 = self._interpolateSpectrum(spa1, spa2, metallicity)
spa6 = self._interpolateSpectrum(spa3, spa4, metallicity)
spa7 = self._interpolateSpectrum(spa5, spa6, phase)
sp = spa7[0]
self._wavetable = wave * 1e4
self._throughputtable = sp
self.waveunits = S.units.Units(waveunits.lower())
self.warnings = {}
def _getSpectrum(self, parlist, catdir, wave_output=False):
name = parlist[3]
filename = name.split('[')[0]
column = name.split('[')[1][:-1]
filename = f"{self.catalog_folder}/{filename}"
# sp = S.spectrum.TabularSpectralElement(filename, thrucol=column)
with fits.open(filename) as td:
sp = td[1].data.field(column)
wave = td[1].data.field('wavelength')
# waveunits = td[1].header['tunit1']
waveunits = 'micron'
result = []
for member in parlist:
result.append(member)
result.pop()
result.append(sp)
if wave_output:
return result, wave, waveunits
return result
def _sptype2num(spectral_type):
"""
convert spectral type string to number, for interpretation
case0: normal case
- M1V: 6.1, 5
- O5IV: 0.5, 4
- F3V: 3.3, 5
- K4.5II: 5.45, 2
case 1: star type or subtype missing
zero or V will return
- M1: 6.1, 5
- M: 6.0, 5
case 2: spectral type + subtype
subtype will be ignored
- K3Vvar: 5.3, 5
- F6Vbwvar: 3.6, 5
- K0IV SB: 5.0, 4
- F5V+: 3.5, 5
case 3: multi spectral type
only the first sptype is used
- G5IV/V +K1IV: 4.5, 4
- F7IV-V: 3.7, 4
- O4/O5IV: 0.4, 0
- G3/G5V: 4.3, 0
case 4: illegal input
ValueError will be raised
"""
obafgkm = 'OBAFGKML'
spectral_type = spectral_type.upper()
# match spectral type such as M1V, O5IV, F3V, K4.5II
matched = re.match(
RF'^([{obafgkm}])([0-9]\d*\.?\d*)*([IV]*)', spectral_type)
if not matched:
raise ValueError(f"illegal spectral type input: {spectral_type}")
shorttype = obafgkm.find(matched.group(1))
subtype = 0.0
if matched.group(2):
subtype = float(matched.group(2))
stlist = ['O', 'I', 'II', 'III', 'IV', 'V']
startype = dict(zip(stlist, range(len(stlist)))).get(matched.group(3), 5)
return shorttype + subtype / 10, startype
def _spnum2teff(spn, stn):
"""
interpret of spectral number (by __sptype2num) to get t_eff and log_g
look up table from the document of ck04model
"""
with open(f'{cpism_refdata}\\target_model\\sptype2teff_lut.json', 'r') as fid:
sptype_teff_lut = json.load(fid)
def _interp(spn, stn):
lut = sptype_teff_lut[f'{stn}']
teff = np.interp(spn, lut[0], lut[1])
logg = np.interp(spn, lut[0], lut[2])
return teff, logg
stn = 5 if stn not in [1, 2, 3, 4, 5] else stn
if stn in [1, 3, 5]:
return _interp(spn, stn)
else:
teff_low, logg_low = _interp(spn, stn-1)
teff_high, logg_high = _interp(spn, stn+1)
return (teff_high + teff_low)/2, (logg_high + logg_low)/2
def star_photlam(
magnitude: float,
sptype: str,
is_blackbody: bool = False,
mag_input_band: str = 'f661'):
"""
genrate flux spectrum of a star by its observal magnitude and spectral type
Parameters
----------
magnitude: float
magnitude of the star
sptype: str
spectral type of the star
is_blackbody: bool
if True, use blackbody spectrum instead of ck04model
mag_input_band: str
bandpass of the input magnitude, default is f661
Returns
-------
pysynphot.spectrum.ArraySpectrum
spectrum of the star in photlam unit
"""
spn, stn = _sptype2num(sptype)
t_eff, log_g = _spnum2teff(spn, stn)
log.debug(f"{sptype} star => [{t_eff=:}, {log_g=:}]; {is_blackbody=:}")
filter = filter_throughput(mag_input_band)
if not is_blackbody:
METALLICITY = 0
spectrum = S.Icat('ck04models', t_eff, METALLICITY, log_g)
else:
spectrum = S.BlackBody(t_eff)
star_sp = spectrum.renorm(magnitude, MAG_SYSTEM, filter)
star_sp.convert('photlam')
return star_sp
def bcc_spectrum(
coe_cloud: float,
coe_metalicity: float
):
"""Albedo spectrum of planet using BCC model (Batalha et al. 2018),
Parameters
----------
coe_cloud: float
The sedimentation efficiency of the cloud. log(f_sed), range from -2 to 2. 2 is for cloud free case.
coe_metalicity: float
The metallicity of the planet. log(Z/Zsun), range from 0 to 2.
Returns
-------
pysynphot.spectrum.ArrayBandpass
albedo spectrum of the planet
"""
spec = AlbedoCat(0, coe_metalicity, coe_cloud)
spec.convert('angstrom')
return spec
def hybrid_albedo_spectrum(
coe_b: float,
coe_r: float):
"""Albedo spectrum of planet using hybrid-jupiter-neptune model (Lacy et al. 2018)
jupiter and neptune spectrum is from Karkoschka’s 1994
Parameters
----------
coe_b: float
coefficient of blue spectrum, 1 for jupiter, 0 for neptune
coe_r: float
coefficient of red spectrum, 1 for jupiter, 0 for neptune
Returns
-------
pysynphot.spectrum.ArrayBandpass
albedo spectrum of the planet
"""
log.debug(f"planet hybrid spectrum with {coe_b=:}, {coe_r=:}")
model = fits.getdata(HYBRID_MODEL_FILE)
spec = model[1, :] * coe_r
spec += model[2, :] * coe_b
spec += model[3, :] * (1 - coe_r)
spec += model[4, :] * (1 - coe_b)
albedo = S.ArrayBandpass(
wave=model[0, :],
throughput=spec,
waveunits='nm'
)
albedo.convert('angstrom')
return albedo
def extract_target_x_y(
target: dict,
ra0: str = None,
dec0: str = None):
"""
extract x, y of target from target dict
Parameters
----------
target: dict
target dict. must contain either (ra, dec) or (pangle, spearation)
ra0: str
ra of center star. must be provided if (ra, dec) of target is used
dec0: str
dec of center star. must be provided if (ra, dec) of target is used
Returns
-------
x, y: float
x, y of target in arcsec
Raises
------
ValueError
if (ra, dec) of target is used but (ra, dec) of center star is not provided.
ValueError
one of (ra, dec) or (pangle, spearation) is not provided.
"""
def _pa2xy(p_angle, separation):
p_angle_rad = p_angle * DEG2RAD
x = separation * np.sin(p_angle_rad)
y = separation * np.cos(p_angle_rad)
log.debug(f"({p_angle=:}, {separation=:}) => ({x=:}, {y=:})")
return x, y
if 'pangle' in target.keys() and 'separation' in target.keys():
return _pa2xy(target['pangle'], target['separation'])
if 'ra' not in target.keys() or 'dec' not in target.keys():
raise ValueError(
'either (ra, dec) or (pangle, separation) needed in target dict')
if ra0 is None or dec0 is None:
raise ValueError(
'(ra, dec) of center star must be provided if (ra, dec) of bkstar is used'
)
ra, dec = target['ra'], target['dec']
log.debug(f"target: {ra=:}, {dec=:}, center star: {ra0=:}, {dec0=:}")
cstar = SkyCoord(ra0, dec0)
bkstar = SkyCoord(ra, dec)
separation = cstar.separation(bkstar).arcsec
p_angle = cstar.position_angle(bkstar).degree
x, y = _pa2xy(p_angle, separation)
return x, y
class TargetOjbect(object):
"""A helper class to generate target spectrum and albedo spectrum
Attributes
----------
x: float
x of target in arcsec
y: fload
y of target in arcsec
ra: str
ra string for center star, such as '15d23m05s'
dec: str
dec string for center star
distance: float
distance of center star in pc
image: 2d np.array
image of the target
spectrum: pysynphot.spectrum.Spectrum
spectrum of the target
"""
def __init__(self, info, cstar=None, sp_model=None):
"""Initialize a target object
Parameters
----------
info: dict
target info, see Example for more details
cstar: TargetOjbect or None
center star object bounded, if None, means the target is the center star itself
center star object is used to calculate the x, y of each target according to its ra and dec
also center star's distance is used to calculate seperation of planet
sp_model: str
keyword to compatible with V1.0 input. See Notes
spectral type of the target, can be 'star' or 'planet'
Notes
-----
keyargv 中的sp_model 完全用来用来兼容V1.0版的输入,执行逻辑如下:
1. 如果没有设置sp_model，对应了V2.0的输入，则sp_model从info中获取
2. 如果设置了sp_model，则区分为star和planet
3. 如果sp_model为star, 则判断是否设置了isblackbody，如果是的话，则sp_model改为blackbody
4. 如果sp_model为planet，则进一步区分为hybrid_planet 和 bcc_planet 和 template_planet
注意：2.0版本的输入尽量不要使用sp_model关键字,逻辑还是有点复杂的，容易乱，而应该直接在info中设置
Examples:
--------
cstar = {
'magnitude': 0,
'ra': '120d',
'dec': '40d',
'distance': 10,
'sptype': 'G0III'
}
stars = {
'magnitude': 15,
'ra': '120.001d',
'dec': '40.001d',
'sptype': 'F0III',
'sp_spectrum': 'blackbody',
}
planets = {
'radius': 2,
'pangle': 60,
'coe_cloud': 0.3,
'coe_metal': 0.7,
'separation': 0.5,
'phase_angle': 90,
'sp_spectrum': 'hybrid_planet',
'image': 'extend_planet.fits'
}
# planet using input costum albedo spectrum!
# Note that the albedo spectrum is not normalized!
# so the contrast of the planet sould be considered in the input file by user!
# The input file is in pysynphot.FileSpectralElement format.
# See the documentation of pysynphot for details.
planets = {
'pangle': 60,
'separation': 0.5,
'sp_spectrum': 'template_planet',
'template': 'planet_albedo.fits'
}
"""
if sp_model is None:
sp_model = info.get('sp_type', 'star') # 如果不设置，默认为star
elif sp_model == 'star' and info.get('isblackbody', False):
#原先的star 现在分为star和blackbody两个类型
sp_model = 'blackbody'
elif sp_model == 'planet':
# 原先的planet 现在分为hybrid_planet 和 bcc_planet两类
sp_model = info.get('model', 'hybrid_planet')
if sp_model not in ['bcc_planet', 'hybrid_planet', 'template_planet']:
Warning(f'Unknown planet model {sp_model}, using default[hybrid] planet model' )
self.sp_model = sp_model
if cstar is None:
self.x, self.y = 0, 0
self.ra = info['ra']
self.dec = info['dec']
self.distance = info.get('distance', None)
else:
self.x, self.y = extract_target_x_y(info, cstar.ra, cstar.dec)
self.image = None
if 'image' in info.keys():
file = info['image']
if os.path.exists(file):
self.image = fits.getdata(file)
else:
self.image = None
Warning(f"extend Ojbect File {file} not found, return None!")
if self.sp_model == 'blackbody':
self.spectrum = star_photlam(
info['magnitude'],
info['sptype'],
is_blackbody=True,
mag_input_band=info.get('mag_input_band', 'f661')
)
if self.sp_model == 'template_star':
self.spectrum = S.FileSpectrum(info['template'])
if self.sp_model == 'star':
self.spectrum = star_photlam(
info['magnitude'],
info['sptype'],
is_blackbody=False,
mag_input_band=info.get('mag_input_band', 'f661')
)
if self.sp_model in ['hybrid_planet', 'bcc_planet']:
planet = info
phase_angle = planet.get('phase_angle', 90)
radius = planet.get('radius', 1)
if cstar.distance is None:
raise ValueError('distance of center star must be provided if planet is added')
if planet.get('contrast', None) is not None:
contrast = planet['contrast']
else:
contrast = planet_contrast(
self.x * cstar.distance,
self.y * cstar.distance,
phase_angle,
radius,
)
if sp_model == 'hybrid_planet':
coe_blue, coe_red = planet.get('coe_b', 1), planet.get('coe_r', 1)
albedo_spect = hybrid_albedo_spectrum(coe_blue, coe_red) * contrast
elif sp_model == 'bcc_planet':
coe_c, coe_m = planet.get('coe_cloud', 2), planet.get('coe_metal', 0)
albedo_spect = bcc_spectrum(coe_c, coe_m) * contrast
if sp_model == 'template_planet':
albedo_spect = S.FileBandpass(planet['template'])
self.spectrum = cstar.spectrum * albedo_spect
def planet_contrast(
planet_x_au: float,
planet_y_au: float,
phase_angle: float,
radius: float):
"""
calculate the contrast of a planet
Parameters
----------
planet_x_au: float
x position of the planet in au
planet_y_au: float
y position of the planet in au
phase_angle: float
phase angle of the planet in degree
radius: float
radius of the planet in jupiter radius
Returns
-------
contrast: float
contrast of the planet
"""
separation = np.sqrt(planet_x_au**2 + planet_y_au**2)
phase_angle = phase_angle * DEG2RAD
if np.sin(phase_angle) < 1e-9:
raise ValueError('sin(phase_angle) can not be 0')
sep_3d = separation / np.sin(phase_angle)
# Lambert Scattering phase function
# from Madhusudhan and Burrows 2012 equation 33.
phase = (np.sin(phase_angle) + (np.pi - phase_angle)
* np.cos(phase_angle)) / np.pi
log.debug(f'alpha: {phase_angle/np.pi*180} {phase=}')
contrast = (radius / sep_3d * R_JUPITER_METER / AU_METER)**2 * phase
return contrast
def spectrum_generator(
targets: dict):
"""
generate the spectrum due to the input target list
Parameters
----------
targets: dict
target dictionary which contains keyword 'cstar' (necessary), 'stars'(optional), 'planets' (optional).
The values are:
- cstar: dict
- center star information. must contain keywords ra, dec, distance, magnitude, sptype
- stars: list of dict, optional, not preferred, compatible with old version
- list of background stars. each dict must contain ra, dec, magnitude, sptype
- planets: list of dict, optional, not preferred, compatible with old version
- list of planets. each dict must contain pangle, separation, magnitude, radius
- objects: list of dict, optional, recommended! V2.0 new!
- list of targets. each dict must contain ra, dec, magnitude, sptype
Returns
-------
obj_sp_list: list
list of [x, y, spectrum, image] of each target
"""
cstar = targets['cstar']
stars = targets.get('stars', [])
planets = targets.get('planets', [])
objects = targets.get('objects', [])
obj_sp_list = []
cstar_obj = TargetOjbect(cstar)
obj_sp_list.append([cstar_obj.x, cstar_obj.y, cstar_obj.spectrum, cstar_obj.image])
for star in stars:
star_obj = TargetOjbect(star, sp_model='star', cstar=cstar_obj)
obj_sp_list.append([star_obj.x, star_obj.y, star_obj.spectrum, star_obj.image])
for planet in planets:
planet_obj = TargetOjbect(planet, sp_model='planet', cstar=cstar_obj)
obj_sp_list.append([planet_obj.x, planet_obj.y, planet_obj.spectrum, planet_obj.image])
for target in objects:
target_obj = TargetOjbect(target, cstar=cstar_obj)
obj_sp_list.append([target_obj.x, target_obj.y, target_obj.spectrum, target_obj.image])
return obj_sp_list
\ No newline at end of file
CpicImgSim/utils.py deleted 100644 → 0
View file @ b7f8c4fa
import numpy as np
import scipy.ndimage as nd
import logging
import random
# DO NOT IMPORT CPICIMGSIM MODULES HERE
class Logger(object):
def __init__(self, filename, level='INFO'):
self.logger = logging.getLogger('cpism_log')
self.logger.setLevel(logging.DEBUG)
shinfo = logging.StreamHandler()
onlyinfo = logging.Filter()
onlyinfo.filter = lambda record: (record.levelno < logging.WARNING)
fmtstr = '%(message)s'
shinfo.setFormatter(logging.Formatter(fmtstr)) # 设置屏幕上显示的格式
shinfo.setLevel(logging.INFO)
shinfo.addFilter(onlyinfo)
sh = logging.StreamHandler()
fmtstr = '!%(levelname)s!: %(message)s [%(filename)s - %(funcName)s (line: %(lineno)d)]: '
sh.setFormatter(logging.Formatter(fmtstr)) # 设置屏幕上显示的格式
sh.setLevel(logging.WARNING)
th = logging.FileHandler(filename) # 往文件里写入#指定间隔时间自动生成文件的处理器
fmtstr = '%(asctime)s %(filename)s [%(funcName)s] - %(levelname)s: %(message)s'
th.setFormatter(logging.Formatter(fmtstr)) # 设置文件里写入的格式
th.setLevel(logging.__dict__.get(level.upper()))
self.logger.addHandler(shinfo)
self.logger.addHandler(sh)
self.logger.addHandler(th)
def random_seed_select(seed=-1):
"""
Select a random seed for numpy.random and return it.
"""
if seed == -1:
seed = random.randint(0, 2**32-1)
np.random.seed(seed)
return seed
def region_replace(
background: np.ndarray,
front: np.ndarray,
shift: list,
bmask: float = 1.0,
fmask: float = 1.0,
padded_in: bool = False,
padded_out: bool = False,
subpix: bool = False
):
"""
replace a region of the background with the front image.
Parameters
----------
background: np.ndarray
The background image.
front: np.ndarray
The front image.
shift: list
The [x, y] shift of the front image. Unit: pixel.
Relative to the lower-left corner of the background image.
[0, 0] means the lower-left corner of the front image is at the lower-left corner of the background image.
bmask: float
The mask of the background image. Default: 1.0
0.0 means the background image is masked.
1.0 means the background image is fully added.
fmask: float
The mask of the front image. Default: 1.0
0.0 means the front image is masked (not added).
1.0 means the front image is fully added.
padded_in: bool
Whether the input background image is padded. Default: False
In the function, the background image is padded by the size of the front image.
If True, means the background image is padded.
padded_out: bool
Whether the output image is padded. Default: False
In the function, the background image is padded by the size of the front image.
If True, means the output image is padded.
padded_in and padded_out are designed for the case that replace_region fuction is called multiple times.
subpix: bool
Whether the shift is subpixel. Default: False
If True, the shift is subpixel, using scipy.ndimage.shift to shift the front image.
If False, the shift is integer, using numpy slicing to shift the front image.
Returns
-------
np.ndarray
The output image.
shape = background.shape if padded_out = False
shape = background.shape + 2 * front.shape if padded_out = True
"""
int_shift = np.array(shift).astype(int)
b_sz = np.array(background.shape)
f_sz = np.array(front.shape)
if padded_in:
padded = background
b_sz = b_sz - f_sz * 2
else:
padded = np.pad(background, ((f_sz[0], f_sz[0]), (f_sz[1], f_sz[1])))
if np.any((int_shift < -b_sz) | (int_shift > b_sz)):
if padded_out:
return padded
return background
if subpix:
subs = np.array(shift) - int_shift
front = nd.shift(front, (subs[0], subs[1]))
int_shift += f_sz
roi_y = int_shift[1]
roi_x = int_shift[0]
padded[roi_y: roi_y+f_sz[0], roi_x:roi_x+f_sz[1]] *= bmask
padded[roi_y: roi_y+f_sz[0], roi_x:roi_x+f_sz[1]] += fmask * front
if padded_out:
return padded
return padded[f_sz[0]:b_sz[0]+f_sz[0], f_sz[1]:b_sz[1]+f_sz[1]]
csst_cpic_sim/__init__.py
View file @ 95f81f91
from .main import quick_run, observation_simulation
from .optics import make_focus_image, focal_mask
# from .main import quick_run, observation_simulation
from .optics import focal_mask
from .target import star_photlam, planet_contrast, extract_target_x_y, spectrum_generator
from .camera import EMCCD, CosmicRayFrameMaker, sky_frame_maker
from .camera import CosmicRayFrameMaker, sky_frame_maker
from .config import __version__
__all__ = [
"EMCCD",
"CosmicRayFrameMaker",
"sky_frame_maker",
"star_photlam",
"planet_contrast",
"extract_target_x_y",
"spectrum_generator",
"make_focus_image",
"focal_mask",
"quick_run",
"observation_simulation"
]
\ No newline at end of file
csst_cpic_sim/camera.py
View file @ 95f81f91
import os
import yaml
import math
import numpy as np
import scipy.ndimage as nd
from astropy.io import fits
import matplotlib.pyplot as plt
from .config import cpism_refdata, solar_spectrum, MAG_SYSTEM
from .config import config, S
from .utils import region_replace, random_seed_select
from .io import log
from .optics import filter_throughput
cpism_refdata = config['cpism_refdata']
MAG_SYSTEM = config['mag_system']
solar_spectrum = S.FileSpectrum(config['solar_spectrum'])
solar_spectrum.convert('photlam')
def sky_frame_maker(band, skybg, platescale, shape):
"""
generate a sky background frame.
......@@ -274,32 +277,17 @@ class CosmicRayFrameMaker(object):
return image
class EMCCD(object):
"""
EMCCD camera class
Parameters
----------
config_file : str
config file name
Attributes
----------
switch : dict
switch for each camera effects, including:
- 'flat': bool,
- 'dark': bool,
- 'stripe': bool,
- 'cic': bool,
- 'cte': bool,
- 'badcolumn': bool,
- 'nonlinear': bool,
- 'cosmicray': bool,
- 'blooming': bool,
"""
class CpicVisEmccd(object):
def __init__(self, config_dict=None):
self._defaut_config()
if config_dict is not None:
old_switch = self.switch
self.__dict__.update(config_dict) #not safe, be careful
old_switch.update(self.switch)
self.switch = old_switch
self.config_init()
def __init__(self, config_file="emccd_config.yaml"):
def _defaut_config(self):
self.plszx = 1024
self.plszy = 1024
self.pscan1 = 8
......@@ -311,43 +299,70 @@ class EMCCD(object):
self.ldark = 16
self.rdark = 16
self.fullwell = 80_000
self.em_fullwell = 780_000
# if config file exists, load it, otherwise use default values
config_file = cpism_refdata + '/camera/' + config_file
log.debug(f"emccd config file: {config_file}")
if os.path.exists(config_file):
self.load_config(config_file)
else: # pragma: no cover
# set default values for EMCCD
# note: these values are default values, you can change them by load_config()
# ↓↓↓↓↓↓↓start default values setting↓↓↓↓↓↓
self.readout_noise = 40
self.ph_per_adu = 8
self.bias_level = 30
self.max_adu = 16_383
self.switch = {
'flat': True,
'bias_vp': True,
'bias_hp': True,
'bias_ci': True,
'bias_shift': True,
'cic': True,
'dark': True,
'stripe': True,
'cic': False,
'cte': False,
'flat': True,
'badcolumn': True,
'nonlinear': False,
'shutter': True,
'cte': True,
'nonlinear': True,
'cosmicray': True,
'blooming': False,
'blooming': True,
'em_blooming': True,
}
self.dark_file = cpism_refdata + '/camera/emccd_dark_current.fits'
self.flat_file = cpism_refdata + '/camera/emccd_flat_field.fits'
self.cic_file = cpism_refdata + '/camera/emcid_cic.fits'
self.cic_file = cpism_refdata + '/camera/emccd_cic2.fits'
self.bad_col_file = cpism_refdata + '/camera/emccd_bad_columns.fits'
# ↑↑↑↑↑↑↑end default values setting↑↑↑↑↑↑
# note: these values are default values, you can change them by load_config()
self.cic = None
self.dark = None
self.flat = None
self.fullwell = 80_000
self.em_fullwell = 500_000 #780_000
self.em_cte = 0.9996
self.emreg_cal_num = 10 # 用来加速计算
self.emreg_num = 604
self.readout_speed = 6.25e6 # Hz
self.readout_time = 0.365 # s
self.heat_speed = 1 / 1000 # voltage / 1000 degree per frame
self.temper_speed = 0.05 # degree per second
self.cooler_temp = -80
self.readout_noise = 160
self.ph_per_adu = 59
self.bias_level = 200
self.vertical_param = [0, 14]
self.horizontal1_param = [5.3/2, 12, 5/2, 50]
self.horizontal2_param = [5.3/4, 16.17, 5/4, 76.65]
self.bias_hp_resd = 2
self.cooler_interfence = 20
self.bias_level_std = 3
self.bias_shift_per_volt = -1.34
self.shift_time = 1 / 1000
self.nonlinear_coefficient = -0.1
self.detector_name = 'EMCCD'
self.ccd_label= 'CCD201-20'
self.pitch_size = 13
def config_init(self):
self._img_index = 0
self._ccd_temp = self.cooler_temp
self._vertical_i0 = np.random.randint(0, 2)
self._frame_read_time = 2080 * 1055 / 6.25e6
self.flat_shape = [self.plszy, self.plszx]
......@@ -357,33 +372,262 @@ class EMCCD(object):
biassz_x = darksz_x + self.pscan1 + self.oscan1
biassz_y = darksz_y + self.pscan2 + self.oscan2
self.image_shape = [biassz_y, biassz_x]
self.bias_shape = [biassz_y, biassz_x]
if self.flat is None:
self.flat = fits.getdata(self.flat_file)
if self.cic is None:
self.cic = fits.getdata(self.cic_file)
if self.dark is None:
self.dark = fits.getdata(self.dark_file)
self.bad_col = fits.getdata(self.bad_col_file)
self.ccd_temp = self.cooler_temp
self.system_time = 0
self.time_syn(0, initial=True)
self.emgain_set(1024, -80)
def CTE_cal(n, N_p, CTE, S_0):
'''
CTE_cal(order of trail pixels, number of pixel transfers, CTE, initial intensity of target pixel)
'''
CTI = 1 - CTE
S_0_n = S_0 * ((N_p * CTI) ** n) / math.factorial(n) * math.exp(-N_p * CTI)
return S_0_n
def cte_201(cte, start=0, length=10):
N_p = 604
S_0 = 1
res = [0]*start
for n in range(length):
s = CTE_cal(n, N_p, cte, S_0)
res.append(s)
return np.array(res)
cti_trail = cte_201(self.em_cte, start=0, length=10)
self.cti_trail = cti_trail / cti_trail.sum()
def em_fix_fuc_fit(self, emgain):
emgain = np.array([emgain]).flatten()
p = [0.01014486, -0.00712984, -0.17163414, 0.09523666, -0.53926089]
def kernel(em):
log_em = np.log10(em)
loglog_g = np.log10(log_em)
loglog_t = np.polyval(p, loglog_g)
log_t = 10**loglog_t
t = 10**log_t - 1
return t
output = []
for em in emgain:
if em <= 1:
output.append(0)
elif em > 80:
output.append(kernel(80))
else:
output.append(kernel(em))
return np.array(output)
def bias_frame(self):
shape = self.bias_shape
TPI = np.pi * 2
# vertical pattern
# 使用一维的曲线描述竖条纹的截面
vp_1d = np.zeros(shape[1])
# 以下代码用于模拟垂直间隔的竖条纹在奇偶幅图像时的不同表现
# 后续相机更新过程中，已经将改特性进行修改，因此不再使用此代码
# vp_1d[0::2] = self.vertical_param[self._vertical_i0]
# self._vertical_i0 = 1 - self._vertical_i0
# vp_1d[1::2] = self.vertical_param[self._vertical_i0]
vp_1d[0::2] = self.vertical_param[self._vertical_i0]
vp_1d[1::2] = self.vertical_param[1 - self._vertical_i0]
vp_frame = np.zeros(shape)
if self.switch['bias_vp']:
vp_frame += vp_1d
# if show: # pragma: no cover
# plt.figure(figsize=(10, 3))
# plt.plot(vp_1d)
# plt.xlim([0, 100])
# plt.xlabel('x-axis')
# plt.ylabel('ADU')
# plt.title('vertical pattern')
# horizontal pattern
# 图像上的横条纹是梳状，分为两个部分，左边大约77个像素是周期小一点，其余的会大一点
boundary = 77 # boundary between left and width
boundary_width = 5 # 左右需要平滑过度一下
y = np.arange(self.bias_shape[0])
hp_left_param = self.horizontal1_param # 实测数据拟合得到的
hp_left_1d = hp_left_param[0] * np.sin(TPI * (y / hp_left_param[1] + np.random.rand()))
hp_left_1d += hp_left_param[2] * np.sin(TPI * (y / hp_left_param[3] + np.random.rand()))
hp_left_frame = np.broadcast_to(hp_left_1d, [boundary+boundary_width, len(hp_left_1d),]).T
hp_right_param = self.horizontal2_param
hp_right_1d = hp_right_param[0] * np.sin(TPI * (y / hp_right_param[1] + np.random.rand()))
hp_right_1d += hp_right_param[2] * np.sin(TPI * (y / hp_right_param[3] + np.random.rand()))
hp_right_frame = np.broadcast_to(hp_right_1d, [shape[1] - boundary, len(hp_right_1d)]).T
combine_profile_left = np.ones(boundary + boundary_width)
combine_profile_left[-boundary_width:] = (boundary_width - np.arange(boundary_width) - 1) / boundary_width
combine_profile_right = np.ones(shape[1] - boundary)
combine_profile_right[:boundary_width] = np.arange(boundary_width)/boundary_width
hp_frame = np.zeros(shape)
hp_frame[:, :boundary+boundary_width] = hp_left_frame * combine_profile_left
hp_frame[:, boundary:] = hp_right_frame * combine_profile_right
# residual frame 横条纹外，还有一个垂直方向的梯度，根据测试数据，使用直线加指数函数的方式来拟合
exp_a, exp_b, lin_a, lin_b = (-1.92377463e+00, 1.32698365e-01, 8.39509583e-04, 4.25384480e-01)
y_trend = exp_a * np.exp(-(y + 1) * exp_b) + y * lin_a - lin_b
# random horizontal pattern generated in frequency domain
# 除了规则横条纹外，还有随机的横条纹，随机的横条纹在频域空间生成，具有相同的频率谱和随机的相位
rsd_freq_len = len(y_trend) * 4
red_freq = np.arange(rsd_freq_len)
red_freq = red_freq - (len(red_freq) - 1)/2
red_freq = np.exp(-red_freq**2 / 230**2) * 240 + np.random.randn(rsd_freq_len)*30
red_freq = np.fft.fftshift(red_freq)
phase = np.random.rand(rsd_freq_len) * TPI
hp_rsd_1d = np.fft.ifft(np.exp(1j * phase) * red_freq)
hp_rsd_1d = np.abs(hp_rsd_1d) * self.bias_hp_resd
hp_rsd_1d = hp_rsd_1d[:rsd_freq_len//4]
hp_rsd_1d = hp_rsd_1d - np.mean(hp_rsd_1d)
hp_rsd_1d = y_trend + hp_rsd_1d
hp_frame = (hp_frame.T + hp_rsd_1d).T
if not self.switch['bias_hp']:
hp_frame *= 0
# if show: # pragma: no cover
# plt.figure(figsize=(10, 3))
# # plt.plot(hp_right_1d)
# plt.plot(hp_rsd_1d)
# # plt.xlim([0, 200])
# plt.xlabel('y-axis')
# plt.ylabel('ADU')
# plt.title('vertical pattern')
# 接上制冷机后，会有亮暗点
#cooler interfence effect
ci_position = 10
ci_sub_struct = 80
ci_sub_exp = 2.5
ci_x_shft = 3
ci_interval = 250 # 6.25MHz readout / 2.5KHz interfence
ci_dn = self.cooler_interfence
npix = shape[0] * shape[1]
n_ci_event = npix // ci_interval
ci_align = np.zeros((n_ci_event, ci_interval))
ci_align[:, ci_position] = np.random.randn(n_ci_event) * ci_dn
ci_align[:, ci_position+1] = np.random.randn(n_ci_event) * ci_dn
yi0 = np.random.randint(0, ci_sub_struct)
xs0 = (ci_interval - ci_position) / (ci_sub_struct / 2)**ci_sub_exp
for yi in range(n_ci_event):
sub_yi = (yi - yi0) % ci_sub_struct
sub_yi = abs(sub_yi - ci_sub_struct / 2)
shiftx = int(sub_yi**ci_sub_exp * xs0)
ci_align[yi, :] = np.roll(ci_align[yi, :], shiftx)
ci_align = np.pad(ci_align.flatten(), (0, npix-n_ci_event*ci_interval))
ci_frame = ci_align.reshape(shape[0], shape[1])
for yi in range(shape[0]):
ci_frame[yi, :] = np.roll(ci_frame[yi, :], yi * ci_x_shft)
if not self.switch['bias_ci']:
ci_frame *= 0
bias_shift = 0
if self.switch['bias_shift']:
bias_shift = (self.volt - 25) * self.bias_shift_per_volt
# 混合在一起
rn_adu = self.readout_noise / self.ph_per_adu
bias_level = self.bias_level + np.random.randn() * self.bias_level_std
bias_frame = vp_frame + ci_frame + hp_frame + bias_level
bias_frame += rn_adu * np.random.randn(shape[0], shape[1]) + bias_shift
return bias_frame
def load_config(self, config_file):
def nonlinear_effect(self, image):
"""
load config file. Only for internal use.
nonlinear effect
"""
fullwell = self.fullwell
nonlinear_coefficient = self.nonlinear_coefficient
log.debug(
f"nonlinear effect added with coefficient {nonlinear_coefficient}")
image += (image / fullwell)**2 * nonlinear_coefficient * fullwell
with open(config_file, 'r') as f:
config = yaml.load(f, Loader=yaml.FullLoader)
return image
def time_syn(self, t, readout=True, initial=False):
if initial:
self.ccd_temp = self.cooler_temp
self.system_time = t
return
dt = np.maximum(t, 0)
heat = 0
if readout:
heat = self.volt * self.heat_speed
self.ccd_temp = heat + self.cooler_temp + (self.ccd_temp - self.cooler_temp) * np.exp(-dt * self.temper_speed)
if self.ccd_temp < self.cooler_temp: #
self.ccd_temp = self.cooler_temp
self.system_time += dt
# def em_cte(self, img):
# i_shift = 1
# cte_coe = 0.1
# img_shift_i = np.zeros_like(img)
# img_shift_i = np.random.poisson(img * cte_coe)
# pass
def emgain_set(self, em_set, ccd_temp=None, self_update=True):
if ccd_temp is None:
ccd_temp = self.ccd_temp
volt_coe_a = -0.01828
volt_coe_b = 43.61
volt_func = lambda es: volt_coe_a * es + volt_coe_b
self.volt = volt_func(em_set)
self.switch = config['switch']
volt_3ff = volt_func(int('3ff', 16))
volt_190 = volt_func(int('190', 16))
self.readout_noise = config['readout_noise']
self.ph_per_adu = config['ph_per_adu']
self.bias_level = config['bias_level']
self.max_adu = config['max_adu']
em_coe_c = 0.24
# using the expression of em_b = ln(g199) / constant to make fitting easier
constant = (np.exp(em_coe_c * volt_190) - np.exp(em_coe_c * volt_3ff))
# fitting from the ccd test result
ln_g190 = (-ccd_temp - 7) * 0.0325
em_coe_b = ln_g190 / constant
self.dark_file = cpism_refdata + "/camera/" + config['dark_file']
self.flat_file = cpism_refdata + "/camera/" + config['flat_file']
self.cic_file = cpism_refdata + "/camera/" + config['cic_file']
self.bad_col_file = cpism_refdata + \
"/camera/" + config['bad_col_file']
emgain = np.exp(em_coe_b * np.exp(em_coe_c * self.volt))
emgain = np.maximum(1, emgain)
# print(emgain, em_coe_b, em_coe_c * self.volt, self.volt, np.exp(em_coe_c * self.volt))
if self_update:
self.emgain = emgain
self.emset = em_set
return emgain
def vertical_blooming(self, image):
"""
......@@ -401,18 +645,6 @@ class EMCCD(object):
image[:, x] += np.exp(-(line-y)**2/20**2) * img0[y, x] * 0.2
return np.minimum(image, fullwell)
def nonlinear_effect(self, image):
"""
nonlinear effect
"""
fullwell = self.fullwell
nonlinear_coefficient = 0.1
log.debug(
f"nonlinear effect added with coefficient {nonlinear_coefficient}")
image += (image / fullwell)**2 * nonlinear_coefficient * fullwell
return image
def emregester_blooming(self, image, max_iteration=5):
"""
emregester blooming effect
......@@ -448,59 +680,36 @@ class EMCCD(object):
log.debug(
f'{index}/{max_iteration} loop: saturated pixel number: {n_over}')
line = np.minimum(line, self.em_fullwell)
return line.reshape(image.shape)
def cte(self, image):
"""
cte effect
"""
image = self.emregester_blooming(image, max_iteration=5)
return image
# def add_em_cte_conv(self, image):
# shape = image.shape
# img_line = np.convolve(image.flatten(), self.cti_trail, mode='full')
# return img_line[:shape[0]*shape[1]].reshape(shape)
def readout(self, image_focal, emgain, expt, image_cosmic_ray=None):
"""
emccd readout. Generate a image with emccd readout effect.
def readout(self, image_focal, em_set, expt_set, image_cosmic_ray=False, emgain=None):
expt = expt_set
if expt_set == 0:
expt = 0.001
Parameters
----------
image_focal : numpy.ndarray
image at focal plane. Unit: electron/second
emgain : float
emgain of emccd
expt : float
exposure time. Unit: second
image_cosmic_ray : numpy.ndarray, optional
cosmic ray image. Unit: electron/second, by default None
dt = self.readout_time + expt
self.time_syn(dt, readout=True)
Returns
-------
numpy.ndarray
image with emccd readout effect. Unit: ADU
Notes
-----
1. effects include: dark, flat, cte, blooming, nonlinear, etc. Can be turned on/off by switch.
2. size of input image_focal must be 1024x1024
3. size of output image is 1080x1056 (including overscan and dark reference region)
4. Q.E is not included in this function. It should be included in image_focal. See optics.py for details.
"""
log.debug(
fr"EMCCD readout: {emgain=:}, {expt=:}, image_comic_ray:{'None' if image_cosmic_ray is None else 'Not None'}")
if emgain is None:
self.emgain_set(em_set, self.ccd_temp)
else:
self.emgain = emgain
self.emset = 0
log.debug(
f"camera effects switch={self.switch}"
)
image = image_focal * expt
emgain = self.emgain
image = image_focal.astype(float)
log.debug(f"image total photon: {image.sum()}")
if self.switch['flat']:
image = image * self.flat
if self.switch['nonlinear']:
image = self.nonlinear_effect(image)
darksz_x = self.plszx + self.rdark + self.ldark
darksz_y = self.plszy + self.udark + self.bdark
img_dark = np.zeros((darksz_y, darksz_x))
img_dark = np.zeros(self.dark_shape)
img_dark[
self.bdark:self.plszy+self.bdark,
self.ldark:self.ldark+self.plszx
......@@ -508,14 +717,16 @@ class EMCCD(object):
image = img_dark
if self.switch['dark']:
image += self.dark * expt
image = image + self.dark
if self.switch['cic']:
image += self.cic
image *= expt
if image_cosmic_ray is not None and self.switch['cosmicray']:
image += image_cosmic_ray
if self.switch['nonlinear']:
image = self.nonlinear_effect(image)
if self.switch['blooming']:
image = self.vertical_blooming(image)
......@@ -525,86 +736,64 @@ class EMCCD(object):
deadpix_y = self.bad_col[1, i]
image[deadpix_y:, deadpix_x] = 0
biassz_x = darksz_x + self.pscan1 + self.oscan1
biassz_y = darksz_y + self.pscan2 + self.oscan2
img_bias = np.zeros((biassz_y, biassz_x), dtype=int)
img_bias = np.zeros(self.bias_shape, dtype=int)
seed = random_seed_select()
log.debug(f"photon noise seed: {seed}")
img_bias[
self.pscan2:self.pscan2+darksz_y,
self.pscan1:self.pscan1+darksz_x
self.pscan2: -self.oscan2,
self.pscan1: -self.oscan1
] = np.random.poisson(image)
image = img_bias
if self.switch['cte']:
image = self.cte(image * emgain) / emgain
seed = random_seed_select()
log.debug(f"gamma noise seed: {seed}")
if emgain != 1:
image = np.random.gamma(image, emgain)
image = np.minimum(image, self.em_fullwell)
seed = random_seed_select()
log.debug(f"readout noise seed: {seed}")
image += np.random.randn(biassz_y, biassz_x) * self.readout_noise
image = image / self.ph_per_adu + self.bias_level
if self.switch['shutter']:
img_line = image_focal.sum(axis=0)
image_shutter = np.broadcast_to(img_line, [self.bias_shape[0], self.flat_shape[1]])
image_shutter = image_shutter / self.flat_shape[1] * self.shift_time
image_shutter = np.random.poisson(image_shutter)
image[:, self.pscan1+self.ldark:-self.oscan1-self.rdark] += image_shutter
if self.switch['stripe']:
image += self.add_stripe_effect(image)
image = np.minimum(image, self.max_adu)
image = np.maximum(image, 0)
return image.astype(np.uint16)
def add_stripe_effect(self, image):
"""
add stripe effect
"""
shape = image.shape
if self.switch['cic']:
cic_frame = np.zeros((self.dark_shape[0], self.bias_shape[1])) + self.cic
image[self.pscan2:-self.oscan2, :] += np.random.poisson(cic_frame)
v_width = 1
v_amplitude = 30
v_limit = 0.01
v_base = 10
# em_fix found to fitting the gamma distribution to theoritical one.
# here is the theoritical distribution. See Robbins and Hadwen 2003 for more details.
# >>> pEM = np.exp(np.log(emgain)/self.emreg_num) - 1
# >>> for _ in range(self.emreg_num):
# >>> image += np.random.binomial(image, pEM)
# This code is too slow, so we used a modified gamma
h_width = 20
h_amplitude = 10
h_limit = 0.9
h_base = 20
em_fix = self.em_fix_fuc_fit(emgain) * emgain
image = np.random.gamma(image, em_fix) + image * (emgain - em_fix)
index = np.linspace(0, np.pi, shape[0] * shape[1])
def stripe(width, limit, amplitude, base, axis=0):
seed = random_seed_select()
log.debug(f"stripe noise seed: {seed}")
dim_axis = shape[axis]
dim_other = shape[0] * shape[1] // shape[axis]
value = np.sin(index / width * dim_axis + np.pi *
dim_axis / width * np.random.randint(1024))
if self.switch['em_blooming']:
image = self.emregester_blooming(image)
value = np.maximum(value, -limit)
value = np.minimum(value, limit)
value = (value / limit + limit) / 2 * amplitude + base
value = value.reshape(dim_axis, dim_other)
image = np.maximum(image, 0)
image = image.astype(int)
if axis == 1:
value = value.T
if self.switch['cte']:
big_cte = self.em_cte ** (self.emreg_num / self.emreg_cal_num)
for _ in range(self.emreg_cal_num):
residual = np.random.binomial(image, 1-big_cte)
image[:, 1:] += residual[:, :-1] - residual[:, 1:]
return value
bias = self.bias_frame()
output = stripe(v_width, v_limit, v_amplitude, v_base, axis=1)
output += stripe(h_width, h_limit, h_amplitude, h_base, axis=0)
return output
image = image / self.ph_per_adu + bias
image = np.minimum(image, self.max_adu)
image = np.maximum(image, 0)
# # plt.plot(horizontal_index, horizontal_value)
# # # plt.xlim([0, 6.28])
# # plt.show()
log.debug(
f"emccd paramters: \
emset: {em_set} \
emgain: {emgain} \
expt: {expt} \
readout time: {dt}"
)
# fits.writeto('horizontal_value.fits', output, overwrite=True)
return image.astype(dtype=np.uint16)
# if __name__ == '__main__':
......
csst_cpic_sim/config.py
View file @ 95f81f91
import os, yaml
import warnings
from datetime import datetime
import numpy as np
cpism_refdata = '/nfsdata/share/simulation-unittest/cpic_sim/cpism_refdata/'
print(cpism_refdata)
if not os.path.exists(cpism_refdata): # pragma: no cover
raise Exception(
"Can not find CPISM reference data.")
os.environ['PYSYN_CDBS'] = '/nfsdata/share/simulation-unittest/cpic_sim/starmodel/grp/redcat/trds/'
if not os.path.exists(os.environ.get('PYSYN_CDBS', '/home/ubuntu/Downloads/cpic-img-sim-master/cpic-img-sim-master/trd/grp/redcat/trds/')): # pragma: no cover
raise Exception(
"Can not find PYSYN Stellar reference data.")
config_aim = os.path.dirname(os.path.dirname(__file__))
config_aim = os.path.join(config_aim, 'data/refdata_path.yaml')
config_set = False
def set_config(refdata_path=None):
if refdata_path is None:
print("input cpism refencence data folder")
refdata_path = input()
refdata_path = os.path.abspath(refdata_path)
with open(config_aim, 'w') as f:
yaml.dump(refdata_path, f)
return refdata_path
try:
with open(config_aim, 'r') as f:
cpism_refdata = yaml.load(f, Loader=yaml.FullLoader)
if not os.path.isdir(cpism_refdata):
raise FileNotFoundError('cpism refdata path not found')
config_set = True
except FileNotFoundError:
warnings.warn(f'refdata not setup yet, set it before use')
cpism_refdata = set_config()
config = {}
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
config['bands'] = {
'f661': f'{cpism_refdata}/throughtput/f661_total.fits',
'f743': f'{cpism_refdata}/throughtput/f743_total.fits',
'f883': f'{cpism_refdata}/throughtput/f883_total.fits',
'f565': f'{cpism_refdata}/throughtput/f565_total.fits',
}
config['diameter'] = 2 # in meters
config['platescale'] = 0.016153
config['datamodel'] = f'{cpism_refdata}/io/csst-cpic-l0.yaml'
config['log_dir'] = f'{cpism_refdata}/log'
config['log_level'] = f'info'
config['output'] = f'./'
config['sp2teff_model'] = f'{cpism_refdata}/target_model/sptype2teff_lut.json'
config['dm_pickle'] = f'{cpism_refdata}/optics/dm_model.pkl'
config['pysyn_refdata'] = f'{cpism_refdata}/starmodel/grp/redcat/trds'
config['catalog_folder'] = f'{cpism_refdata}/demo_catalog'
config['csst_format'] = True
config['nsample'] = 5
update_able_keys = [
'apm_file', 'actuator_file', 'aberration', 'log_dir', 'log_level', 'catalog_folder', 'nsample', 'csst_format', 'output', 'check_fits_header'
]
def replace_cpism_refdata(config, output='$'):
aim = cpism_refdata
target = '${cpism_refdata}'
if output != '$':
aim, target = target, aim
for key, value in config.items():
if isinstance(value, str):
config[key] = value.replace(aim, target)
if isinstance(value, dict):
replace_cpism_refdata(value, output)
with open(cpism_refdata + '/cpism_config.yaml', 'r') as f:
new_config = yaml.load(f, Loader=yaml.FullLoader)
replace_cpism_refdata(new_config, None)
config.update(new_config)
os.environ['PYSYN_CDBS'] = config['pysyn_refdata']
__version__ = '2.0.0'
# we need to ignore the warning from pysynphot, because we only use the star models.
with warnings.catch_warnings(): # pragma: no cover
warnings.filterwarnings("ignore")
import pysynphot as S
solar_spectrum = S.FileSpectrum(
f"{os.environ['PYSYN_CDBS']}/grid/solsys/solar_spec.fits")
solar_spectrum.convert('photlam')
config_file = cpism_refdata + '/optics/optics_config.yaml'
if not os.path.exists(config_file): # pragma: no cover
raise FileNotFoundError(f"光学配置文件不存在({config_file})")
with open(config_file, 'r') as f:
optics_config = yaml.load(f, Loader=yaml.FullLoader)
def setup_config(new_config):
config.update(new_config)
config['utc0_float'] = datetime.timestamp(datetime.fromisoformat(config['utc0']))
config['solar_spectrum'] = f"{os.environ['PYSYN_CDBS']}/grid/solsys/solar_spec.fits"
config['aperature_area'] = (config['diameter'] * 50)**2 * np.pi # cm^2
config['default_band'] = list(config['bands'].keys())[0]
config['default_filter'] = config['bands'][config['default_band']]
MAG_SYSTEM = 'abmag'
__version__ = '1.0.0'
setup_config({})
def which_focalplane(band):
"""
......@@ -61,3 +130,24 @@ def which_focalplane(band):
return 'wfs'
raise ValueError(f"未知的波段{band}")
def iso_time(time):
if isinstance(time, str):
_ = datetime.fromisoformat(time)
return time
utc0 = config['utc0']
time0 = datetime.timestamp(datetime.fromisoformat(utc0))
time = datetime.fromtimestamp(time0 + time)
return time.isoformat()
def relative_time(time):
if isinstance(time, float):
return time
if isinstance(time, int):
return float(int)
utc0 = config['utc0']
time0 = datetime.timestamp(datetime.fromisoformat(utc0))
return datetime.timestamp(datetime.fromisoformat(time)) - time0
\ No newline at end of file
csst_cpic_sim/io.py
View file @ 95f81f91
import yaml, os, re
from datetime import datetime
import numpy as np
import pandas as pd
from astropy.io import fits
import yaml
from .config import cpism_refdata, __version__, which_focalplane
from .utils import Logger
import os
from datetime import datetime, timedelta
from astropy.coordinates import SkyCoord
import re
import json
import pandas as pd
config_file = f'{cpism_refdata}/cpism_config.yaml'
with open(config_file, 'r') as fid:
config = yaml.load(fid, Loader=yaml.FullLoader)
output_dir = config['output_dir']
if config['relative_path']:
ref_dir_base = os.path.dirname(cpism_refdata)
output_dir = f'{ref_dir_base}/{output_dir}'
from .config import __version__, which_focalplane
from .utils import Logger
from .config import config, iso_time
log_dir = output_dir + '/LOG'
tmp_dir = config['tmp_dir']
default_output_dir = config['output']
log_level = config['log_level']
header_check = config['check_fits_header']
for dir in ['', 'TMP', 'CAL', 'SCI', 'LOG']:
sub_dir = f"{output_dir}/{dir}"
if not os.path.exists(sub_dir):
os.makedirs(sub_dir)
tmp_folder_path = '.'
if tmp_dir == 'TMP':
tmp_folder_path = output_dir + '/TMP'
log_dir = config['log_dir']
if not os.path.exists(log_dir):
os.makedirs(log_dir)
log = Logger(log_dir+'/cpism_pack.log', log_level).logger
......@@ -53,17 +34,19 @@ def check_and_update_fits_header(header):
None
"""
hdu = 'image'
if 'FILETYPE' in header.keys():
hdu = 'primary'
model_file = f"{cpism_refdata}/io/image_header.json"
if hdu == 'primary':
model_file = f"{cpism_refdata}/io/primary_header.json"
# model_file = f"{cpism_refdata}/io/image_header.json"
# if hdu == 'primary':
# model_file = f"{cpism_refdata}/io/primary_header.json"
model_file = config['datamodel']
with open(model_file, 'r', encoding='utf-8') as fid:
data_model = json.load(fid)
data_model = yaml.load(fid, Loader=yaml.FullLoader)
if 'FILETYPE' in header.keys():
header_model =data_model['HDU0']
hdu = 'hdu0'
else:
header_model = data_model['HDU1']
hdu = 'hdu1'
dm_comment = {}
def print_warning(info):
......@@ -71,37 +54,42 @@ def check_and_update_fits_header(header):
log.warning(info)
# check existance and format of keyword in fits header
for keyword, comment, format, _, _ in data_model:
for _, content in header_model.items():
comment = content['comment']
keyword = content['key']
dtype = content['dtype']
if pd.isnull(comment):
comment = ''
if len(comment) > 46:
comment = comment[:46]
if len(comment) > 47:
# comment = comment[:46]
print_warning(
f"Keyword {keyword} has a comment longer than 46 characters. It will be truncated to 46 characters.")
f"Keyword {keyword} has a comment longer than 47 characters.")
dm_comment[keyword] = comment
if keyword not in header.keys():
print_warning(f"Keyword {keyword} not found in [{hdu}] header.")
print_warning(f"Keyword {keyword} not found in header.")
elif not pd.isnull(format):
elif not pd.isnull(dtype):
value = header[keyword]
# check the type of the value, I for int, R for float, C for str
if isinstance(value, str):
type = 'C'
key_type = 'str'
elif isinstance(value, float):
type = 'R'
key_type = 'flo'
elif isinstance(value, bool):
type = 'L'
key_type = 'boo'
elif isinstance(value, int):
type = 'I'
key_type = 'int'
elif isinstance(value, type(header['COMMENT'])):
key_type = 'str'
else:
type = 'U'
key_type = 'ukn'
if type != format[0]:
if key_type != dtype[0:3]:
print_warning(
f"Keyword {keyword} has wrong type in [{hdu}]. {format[0]} expected, {type} found.")
f"Keyword {keyword} has wrong type in [{hdu}]. {dtype} expected, {key_type} found.")
# check if there are extral keyword in header, and update the comment
for keyword in header.keys():
......@@ -109,7 +97,7 @@ def check_and_update_fits_header(header):
if keyword not in dm_comment.keys():
print_warning(
f"Keyword {keyword} not found in the [{hdu}] data model.")
else:
elif keyword != 'COMMENT': # comment keyword is not allowed to be updated
header[keyword] = (header[keyword], dm_comment[keyword])
return header
......@@ -140,19 +128,19 @@ def obsid_parser(
if len(obsid) != 11:
raise ValueError('Obsid must be 11 digits.')
if obsid[0] != '5':
raise ValueError('Obsid must start with 5 for CPIC')
if obsid[0] != '4':
raise ValueError('Obsid must start with 4 for CPIC')
obstype_dict = {
'00': 'BIAS',
'01': 'DARK',
'02': 'FLAT',
'03': 'BKGD',
'04': 'LASR',
'10': 'SCIE',
'11': 'DENF',
'12': 'CALS',
'15': 'TEMP'
'20': 'BIAS',
'21': 'DARK',
'22': 'FLAT',
'23': 'BKG',
'24': 'LASER',
'01': 'SCI',
'02': 'DSF',
'10': 'CALS',
'00': 'TEMP'
}
obstype = obstype_dict.get(obsid[1:3], 'DEFT')
return obstype
......@@ -204,19 +192,14 @@ def primary_hdu(
these informations are used to genrated the hdu header. Refer to the data model for more information.
"""
camera_config, _ = load_camera_and_optics_config(obs_info['band'])
# camera_config, _ = load_camera_and_optics_config(obs_info['band'])
obsid = obs_info['obsid']
exp_start = gnc_info.get(
'EXPSTART', datetime.now().isoformat(timespec='seconds'))
exp_start = obs_info.get('EXPSTART')
exp_start = datetime.fromisoformat(exp_start)
duartion = (obs_info['expt'] +
camera_config['readout_time']) * obs_info['nframe']
default_end = exp_start + timedelta(seconds=duartion)
exp_end = gnc_info.get('EXPEND', default_end.isoformat(timespec='seconds'))
exp_end = obs_info['EXPEND']
exp_end = datetime.fromisoformat(exp_end)
filename = "CSST_CPIC"
......@@ -227,7 +210,7 @@ def primary_hdu(
filename += f"_{obsid}_X_L0_V01.fits"
type_dir = 'CAL'
if str(obsid)[1] == '1':
if int(f'{obsid}'[1:3]) <= 10:
type_dir = 'SCI'
mjd_dir = f"{datetime_obj_to_mjd(exp_start):.0f}"
......@@ -240,7 +223,7 @@ def primary_hdu(
header['NAXIS'] = 0
header['EXTEND'] = True
header['NEXTEND'] = 1 # + parameters['nframe']
header['GROUPS'] = False
# header['GROUPS'] = False
header['DATE'] = datetime.now().isoformat(timespec='seconds')
heaer_filename = filename[:-4]
......@@ -252,7 +235,8 @@ def primary_hdu(
header['INSTRUME'] = 'CPIC'
header['RADECSYS'] = 'ICRS'
header['EQUINOX'] = 2000.0
header['FITSCREA'] = f'CPISM V{__version__}'
header['FITSSWV'] = f'CPISM V{__version__}'
header['COMMENT'] = ''
cstar = {'ra': '0d', 'dec': '0d'}
if obs_info['target'] != {}:
......@@ -264,16 +248,19 @@ def primary_hdu(
header['OBJECT'] = cstar.get('name', target_name)
header['TARGET'] = target_name
header['OBSID'] = str(obsid)
header['OBJ_RA'] = radec.ra.degree
header['OBJ_DEC'] = radec.dec.degree
header['RA_OBJ'] = radec.ra.degree
header['DEC_OBJ'] = radec.dec.degree
# telescope information
header['REFFRAME'] = 'CSSTGSC-1.0'
header['DATE-OBS'] = exp_start.isoformat(timespec='seconds')
header['SATESWV'] = 'SIMULATION'
header['SATESWV'] = '1'
header['EXPSTART'] = datetime_obj_to_mjd(exp_start)
header['CABSTART'] = gnc_info.get('CABSTART', header['EXPSTART'])
cabstart = gnc_info.get('CABSTART', exp_start.isoformat(timespec='seconds'))
cabstart = iso_time(cabstart)
cabstart_mjd = datetime_obj_to_mjd(datetime.fromisoformat(cabstart))
header['CABSTART'] = cabstart_mjd
header['SUNANGL0'] = gnc_info.get('SUNANGL0', -1.0)
header['MOONANG0'] = gnc_info.get('MOONANG0', -1.0)
header['TEL_ALT0'] = gnc_info.get('TEL_ALT0', -1.0)
......@@ -288,10 +275,14 @@ def primary_hdu(
header['EULER0_1'] = gnc_info.get('EULER0_1', -1.0)
header['EULER0_2'] = gnc_info.get('EULER0_2', -1.0)
header['EULER0_3'] = gnc_info.get('EULER0_3', -1.0)
header['RA_PNT0'] = gnc_info.get('RA_PNT0', header['OBJ_RA'])
header['DEC_PNT0'] = gnc_info.get('DEC_PNT0', header['OBJ_DEC'])
header['RA_PNT0'] = gnc_info.get('RA_PNT0', header['RA_OBJ'])
header['DEC_PNT0'] = gnc_info.get('DEC_PNT0', header['DEC_OBJ'])
header['EXPEND'] = datetime_obj_to_mjd(exp_end)
cabend = gnc_info.get('CABEND', exp_end.isoformat(timespec='seconds'))
cabend = iso_time(cabend)
cabend_mjd = datetime_obj_to_mjd(datetime.fromisoformat(cabend))
header['CABEND'] = gnc_info.get('CABEDN', header['EXPEND'])
header['SUNANGL1'] = gnc_info.get('SUNANGL1', header['SUNANGL0'])
header['MOONANG1'] = gnc_info.get('MOONANG1', header['MOONANG0'])
......@@ -309,8 +300,8 @@ def primary_hdu(
header['RA_PNT1'] = gnc_info.get('RA_PNT1', header['RA_PNT0'])
header['DEC_PNT1'] = gnc_info.get('DEC_PNT1', header['DEC_PNT0'])
header['EXPTIME'] = (exp_end - exp_start).total_seconds()
header['EPOCH'] = float(exp_start.year)
header['EXPTIME'] = (exp_end - exp_start).total_seconds()
header['CHECKSUM'] = '0000000000000000'
header['DATASUM'] = '0000000000'
......@@ -327,35 +318,7 @@ def primary_hdu(
return hdu
def load_camera_and_optics_config(band):
"""
Load camera and optics configuration from reference data.
Parameters
----------
band : str
Band name.
Returns camera_config, optics_config : dict, dict
"""
camera = which_focalplane(band)
if camera == 'vis':
config_file = 'emccd_config.yaml'
elif camera == 'nir':
raise ValueError('NIR camera is not supported yet')
config_file = 'nir_config.yaml'
with open(f"{cpism_refdata}/camera/{config_file}", 'r') as fid:
camera_config = yaml.load(fid, Loader=yaml.FullLoader)
with open(f"{cpism_refdata}/optics/optics_config.yaml", 'r') as fid:
optics_config = yaml.load(fid, Loader=yaml.FullLoader)[camera]
return camera_config, optics_config
def frame_header(obs_info, index, bunch_start, primary_header={}):
def frame_header(obs_info, index, primary_header, camera_dict={}):
"""
Generate the header for a single frame.
......@@ -365,10 +328,8 @@ def frame_header(obs_info, index, bunch_start, primary_header={}):
Dictionary of parameters. See `save_fits` function.
index : int
Frame index.
bunch_start : str
Start time of the bunch.
primary_header : dict (optional)
Primary header. default: {}
primary_header : dict
Primary header
Returns
-------
......@@ -377,9 +338,9 @@ def frame_header(obs_info, index, bunch_start, primary_header={}):
"""
header = fits.Header()
camera_config, optics_config = load_camera_and_optics_config(
obs_info['band'])
# camera_config, optics_config = load_camera_and_optics_config(
# obs_info['band'])
camera_config = camera_dict
plszx = camera_config['plszx']
plszy = camera_config['plszy']
pscan1 = camera_config['pscan1']
......@@ -398,48 +359,49 @@ def frame_header(obs_info, index, bunch_start, primary_header={}):
header['NAXIS'] = 2
header['NAXIS1'] = 1080
header['NAXIS2'] = 1050
header['PCOUNT'] = 0
header['GCOUNT'] = 1
header['BSCALE'] = 1
header['BZERO'] = 32768
header['EXTNAME'] = 'IMAGE'
header['EXTVER'] = 1
header['BSCALE'] = 1.0
header['BZERO'] = 32768.0
header['BUNIT'] = 'ADU'
header['FILTER'] = obs_info['band']
header['DETSN'] = '00000000000'
header['DETSN'] = '0'
header['DETNAME'] = camera_config['detector_name']
header['CHIPLAB'] = camera_config['ccd_label']
header['CHIPTEMP'] = float(camera_config['chip_temp'])
header['DEWTEMP'] = float(camera_config['dewar_temp'])
header['DEWTEMP'] = float(camera_config['cooler_temp'])
frame_info = obs_info['frame_info'][index]
header['CHIPTEMP'] = float(frame_info['chiptemp'])
header['DETSIZE'] = f"{imgszx} * {imgszy}"
header['IMGINDEX'] = index
frame_time = obs_info['expt'] + camera_config['readout_time']
bunch_start = datetime.fromisoformat(bunch_start)
expstart = bunch_start + timedelta(seconds=frame_time * (index - 1))
bunch_start_mjd = datetime_obj_to_mjd(bunch_start)
header['IMGINDEX'] = index + 1
utc0_float = config['utc0_float']
ra0 = primary_header.get('RA_PNT0', -1.0)
dec0 = primary_header.get('DEC_PNT0', -1.0)
pa0 = primary_header.get('POS_ANG0', -1.0)
cab0 = primary_header.get('CABSTART', bunch_start_mjd)
delta_t = frame_time * index
bunch_end = bunch_start + timedelta(seconds=delta_t)
bunch_end_mjd = datetime_obj_to_mjd(bunch_end)
cab0 = primary_header.get('CABSTART', obs_info['frame_info'][0]['expt_start'])
ra1 = primary_header.get('RA_PNT1', ra0)
dec1 = primary_header.get('DEC_PNT1', dec0)
pa1 = primary_header.get('POS_ANG1', pa0)
cab1 = primary_header.get('CABEND', bunch_end_mjd)
cab1 = primary_header.get('CABEND', cab0)
frame_stamp = frame_info['expt_start'] + utc0_float
frame_mjd = datetime_obj_to_mjd(datetime.fromtimestamp(frame_stamp))
img_cab = datetime_obj_to_mjd(expstart)
cab0 = iso_time(cab0)
cab1 = iso_time(cab1)
cab0 = datetime_obj_to_mjd(datetime.fromisoformat(cab0))
cab1 = datetime_obj_to_mjd(datetime.fromisoformat(cab1))
ratio = (img_cab - cab0)/(cab1 - cab0)
ratio = (frame_mjd - cab0)/(cab1 - cab0)
ra = ra0 + (ra1 - ra0) * ratio
dec = dec0 + (dec1 - dec0) * ratio
pa = pa0 + (pa1 - pa0) * ratio
header['IMG_EXPT'] = expstart.isoformat(timespec='seconds')
header['IMG_EXPT'] = datetime.fromtimestamp(frame_stamp).isoformat(timespec='seconds')
header['IMG_CABT'] = header['IMG_EXPT']
header['IMG_DUR'] = float(obs_info['expt'])
......@@ -448,8 +410,8 @@ def frame_header(obs_info, index, bunch_start, primary_header={}):
header['IMG_DEC'] = pa
header['DATASECT'] = f"{plszx} * {plszy}"
header['PIXSCAL'] = optics_config['platescale']
header['PIXSIZE'] = camera_config['pitch_size']
header['PIXSCAL'] = frame_info['platescale']
header['PIXSIZE'] = float(camera_config['pitch_size'])
header['NCHAN'] = 1
header['PSCAN1'] = pscan1
header['PSCAN2'] = pscan2
......@@ -467,7 +429,7 @@ def frame_header(obs_info, index, bunch_start, primary_header={}):
radec = SkyCoord(cstar['ra'], cstar['dec'])
shift = obs_info['shift']
platescale = optics_config['platescale']
platescale = frame_info['platescale']
rotation = np.radians(obs_info['rotation'])
header['WCSAXES'] = 2
......@@ -481,7 +443,6 @@ def frame_header(obs_info, index, bunch_start, primary_header={}):
header['CD1_2'] = -np.sin(rotation)
header['CD2_1'] = np.sin(rotation)
header['CD2_2'] = np.cos(rotation)
header['others'] = 'other'
# Readout information
header['EMGAIN'] = float(obs_info['emgain'])
......@@ -493,7 +454,7 @@ def frame_header(obs_info, index, bunch_start, primary_header={}):
# CPIC information
header['LS_STAT'] = 'OFF'
header['IWA'] = optics_config['mask_width'] / 2
header['IWA'] = frame_info['iwa']
header['WFSINFO1'] = -1.0
header['WFSINFO2'] = -1.0
......@@ -505,7 +466,7 @@ def frame_header(obs_info, index, bunch_start, primary_header={}):
return header
def save_fits_simple(images, obs_info):
def save_fits_simple(images, obs_info, output_folder=None):
"""
Save the image to a fits file with a simple header to TMP directory.
......@@ -551,12 +512,20 @@ def save_fits_simple(images, obs_info):
shift = obs_info['shift']
header['shift'] = f"x:{shift[0]},y:{shift[1]}"
if output_folder is None:
fullname = f"{tmp_folder_path}/{filename}"
else:
fullname = f"{output_folder}/{filename}"
if os.path.exists(output_folder) is False:
os.makedirs(output_folder)
log.debug(f"Output folder {output_folder} is created.")
log.debug(f"save fits file to {fullname}")
fits.writeto(fullname, images, overwrite=True, header=header)
return fullname
def save_fits(images, obs_info, gnc_info, csst_format=True):
def save_fits(images, obs_info, gnc_info, camera_dict={}, csst_format=True, output_folder=None):
"""
Save the image to a fits file.
......@@ -592,7 +561,7 @@ def save_fits(images, obs_info, gnc_info, csst_format=True):
"""
if not csst_format:
save_fits_simple(images, obs_info)
save_fits_simple(images, obs_info, output_folder=output_folder)
return
hdu_header, folder, filename = primary_hdu(obs_info, gnc_info, True)
......@@ -604,16 +573,23 @@ def save_fits(images, obs_info, gnc_info, csst_format=True):
for index in range(images.shape[0]):
header = frame_header(
obs_info,
index + 1,
hdu_header.header['DATE-OBS'],
primary_header=hdu_header.header
index,
hdu_header.header,
camera_dict=camera_dict,
)
frame_hdu = fits.ImageHDU(images[index, :, :], header=header)
frame_hdu.add_checksum()
hdu_list.append(frame_hdu)
folder = f"{output_dir}/{folder}"
if output_folder is None:
folder = f"{default_output_dir}/{folder}"
else:
folder = f"{output_folder}/{folder}"
if not os.path.exists(folder):
os.makedirs(folder)
log.debug(f'make new folder {folder}')
hdu_list.writeto(f"{folder}/{filename}", overwrite=True)
full_path = f"{folder}/{filename}"
log.debug(f'save fits file: {full_path}')
hdu_list.writeto(full_path, overwrite=True)
csst_cpic_sim/main.py
View file @ 95f81f91
import numpy as np
import re
from .target import spectrum_generator
from .optics import make_focus_image, focal_mask, optics_config
from .psf_simulation import simulate_psf
from .camera import EMCCD, CosmicRayFrameMaker, sky_frame_maker
from .io import save_fits, log
from .config import which_focalplane
import argparse, sys, tqdm, time, os, yaml
from glob import glob
from datetime import datetime
def psf_function(band, cstar_spectrum, shape, error=0.1):
cstar = True
if shape < 300:
cstar = False
return simulate_psf(error, band, cstar_spectrum, nsample=1, cstar=cstar)
import numpy as np
from .target import spectrum_generator, target_file_load
from .optics import focal_mask, focal_convolve
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
def observation_simulation(
def vis_observation(
target: dict,
skybg: float,
expt: float,
nframe: int,
band: str,
emgain: float,
obsid: int = 51900000000,
emset: int,
obsid: int = 41000000000,
rotation: float = 0,
shift: list = [0, 0],
gnc_info: dict = {},
csst_format: bool = True,
psf_function: callable = psf_function):
"""
Simulate the observation. All-In-One function of the package.
Parameters
-----------
target: dict
The target information. See target.py for details.
skybg: float
magnitude of the skybackground at the input b and. (abmag system)
expt: float
exposure time in second.
nframe: int
number of frames to be simulated.
band: str
the band of the observation. (e.g. 'f661')
emgain: float
the EM gain of the camera.
obsid: int
the observation id. Default is 51900000000.
rotation: float
the rotation angle of the target in degree. Default is 0.
shift: list
the shift of the target in arcsec. Default is [0, 0].
gnc_info: dict
the gnc information. Default is {}. See io.py for details.
csst_format: bool
whether to save the fits file in CSST format. Default is True.
psf_function: callable
the function to generate the psf. See optics.py for details.
Returns
-----------
np.ndarray of the simulated images with shape (nframe, 1088, 1050).
"""
camera = CpicVisEmccd(),
crmaker = CosmicRayFrameMaker(),
nsample: int = 1,
emgain=None,
prograss_bar=None,
output='./'
):
start_time = time.time()
platescale = default_config['platescale']
iwa = default_config['mask_width'] / 2
area = default_config['aperature_area']
expt_start = camera.system_time
target_list = []
if 'cstar' in target.keys():
target_list = spectrum_generator(target)
focal_name = which_focalplane(band)
this_focal_config = optics_config[focal_name]
telescope_config = optics_config['telescope']
area = telescope_config['aperature_area']
image_cube = []
if focal_name == 'vis':
camera = EMCCD()
if emgain is None:
emgain_value = camera.emgain_set(emset, self_update=False)
else:
raise ValueError('Only VIS focal plane is supported.')
platescale = this_focal_config['platescale']
iwa = this_focal_config['mask_width'] / 2
crmaker = CosmicRayFrameMaker()
images = []
emgain_value, emset = emgain, -emgain
params = {
'target': target,
......@@ -90,25 +54,33 @@ def observation_simulation(
'expt': expt,
'nframe': nframe,
'band': band,
'emgain': emgain,
'emset': emset,
'emgain': emgain_value,
'obsid': obsid,
'rotation': rotation,
'shift': shift,
}
paramstr = ', '.join([f'{k}={v}' for k, v in params.items()])
log.debug(f"parameters: {paramstr}")
for i in range(nframe):
all_frame_info = []
log.info(f'Simulation Running: Frame {i+1}/{nframe}')
if prograss_bar:
pg_bar = tqdm.tqdm(total=nframe, leave=False)
focal_frame = make_focus_image(
for i in range(nframe):
if prograss_bar:
pg_bar.update(1)
else:
print(f'\r Simulation Running: Frame {i+1}/{nframe}', end='')
frame_info = {}
frame_info['expt_start'] = camera.system_time
focal_frame = focal_convolve(
band,
target_list,
psf_function,
platesize=camera.flat_shape,
rotation=rotation,
init_shifts=shift,
rotation=rotation,
nsample=nsample,
platesize=camera.flat_shape
)
if skybg is None or skybg > 100:
......@@ -120,181 +92,310 @@ def observation_simulation(
platescale,
camera.flat_shape
)
focal_frame = (focal_frame + sky_bkg_frame) * expt * area
focal_frame = focal_mask(focal_frame, iwa, platescale)
sky_bkg_frame = sky_bkg_frame * area * expt
sky_bkg_frame = focal_mask(sky_bkg_frame, iwa, platescale)
cr_frame = crmaker.make_cr_frame(camera.dark_shape, expt)
img = camera.readout(
focal_frame,
emgain,
camera.time_syn(expt, readout=True)
image = camera.readout(
focal_frame + sky_bkg_frame,
emset,
expt,
image_cosmic_ray=cr_frame
image_cosmic_ray=cr_frame,
emgain=emgain
)
image_cube.append(image)
frame_info['expt_end'] = camera.system_time
frame_info['chiptemp'] = camera.ccd_temp
frame_info['platescale'] = platescale
frame_info['iwa'] = iwa
all_frame_info.append(frame_info)
image_cube = np.array(image_cube)
expt_end = camera.system_time
utc0 = default_config['utc0']
utc0_stamp = datetime.timestamp(datetime.fromisoformat(utc0))
expt_start_iso = datetime.fromtimestamp(utc0_stamp + expt_start)
expt_end_iso = datetime.fromtimestamp(utc0_stamp + expt_end)
params['EXPSTART'] = expt_start_iso.isoformat()
params['EXPEND'] = expt_end_iso.isoformat()
params['frame_info'] = all_frame_info
save_fits(image_cube, params, gnc_info, camera_dict=camera.__dict__.copy(), csst_format=csst_format, output_folder=output)
if prograss_bar:
pg_bar.close()
print(f' Done [{time.time() - start_time:.1f}s] ')
else:
print(f'\r Done [{time.time() - start_time:.1f}s] ')
return image_cube
images.append(img)
images = np.array(images)
save_fits(images, params, gnc_info, csst_format=csst_format)
return images
def quick_run(
def quick_run_v2(
target_str: str,
skymag: float,
band: str,
expt: float,
nframe: int,
emgain: float,
nframe: int,
skybg: float = None,
rotation: float = 0,
shift: list = [0, 0]) -> np.ndarray:
"""
A quick run function to simulate the observation.
Parameters
-----------
target_str: str
The target information in string format.
In the format of "\*5.1/25.3(1.3,1.5)/22.1(2.3,-4.5)" which means a central star
with magnitude 5.1, and two substellar with magnitude 25.3 and 22.1, respectively.
The first number in the parenthesis is the x position in arcsec, and the second is the y position.
skybg: float
magnitude of the skybackground at the input band. (abmag system)
band: str
the band of the observation. (e.g. 'f661')
expt: float
exposure time in second.
nframe: int
number of frames to be simulated.
emgain: float
the EM gain of the camera.
rotation: float (optional)
the rotation angle of the target in degree. Default is 0.
shift: list (optional)
the shift of the target in arcsec. Default is [0, 0].
Returns
-----------
np.ndarray of the simulated images, with shape (nframe, 1088, 1050)
Notes
-----------
1. stars are assumed to be G0III with distance 10pc.
2. magnitude of the star and substellar are assumed to be in the same band.
3. Csst format is not supported.
4. The psf is assumed to be the default one.
5. fits file will be saved in the current directory.
"""
log.info(f'Quick Run: {target_str}')
target_dict = {
'name': 'cal',
}
if (target_str != '') and (target_str[0] == '*'):
objects = target_str[1:].split('/')
cstar_mag = float(objects[0])
cstar = {
'magnitude': cstar_mag,
'ra': '0d',
'dec': '0d',
'sptype': 'G0III',
'distance': 10,
'mag_input_band': band
}
stars = []
for sub_stellar in objects[1:]:
float_regex = R"[+-]?\d+(?:\.\d+)?"
match = re.match(
rf"({float_regex})\(({float_regex}),({float_regex})\)", sub_stellar)
if not match:
raise ValueError('Wrong format for sub stellar.')
mag = float(match.group(1))
x = float(match.group(2))
y = float(match.group(3))
pangle = np.arctan2(x, y) * 180 / np.pi
separation = np.sqrt(x**2 + y**2)
stars.append({
'magnitude': mag,
'pangle': pangle,
'separation': separation,
'sptype': 'G0III',
'mag_input_band': band
})
target_dict = {
'name': target_str[1:],
'cstar': cstar,
'stars': stars,
}
return observation_simulation(
shift: list = [0, 0],
emset_input: bool=False,
cr_frame: bool=True,
camera_effect: bool=True,
prograss_bar=False,
output='./') -> np.ndarray:
print(f'Quick Run: {target_str}')
log.debug(
f"""input parameters:
target_str: {target_str}
skybg: {skybg}
band: {band}
expt: {expt}
nframe: {nframe}
emgain: {emgain}
rotation: {rotation}
shift: {shift}
emset_input: {emset_input}
cr_frame: {cr_frame}
camera_effect: {camera_effect}
prograss_bar: {prograss_bar}
output: {output}
""")
target_dict = target_file_load(target_str)
if target_dict == {}:
target_dict['name'] = 'blank'
emgain_value = emgain
if emset_input:
emgain_value = None
camera = CpicVisEmccd()
if not camera_effect:
camera.switch['bias_vp'] = False
camera.switch['bias_hp'] = False
camera.switch['bias_ci'] = False
camera.switch['bias_shift'] = False
camera.switch['badcolumn'] = False
camera.switch['flat'] = False
camera.switch['cte'] = False
camera.switch['shutter'] = False
camera.switch['nonlinear'] = False
if not cr_frame:
camera.switch['cosmicray'] = False
return vis_observation(
target=target_dict,
skybg=skymag,
skybg=skybg,
expt=expt,
nframe=nframe,
band=band,
emgain=emgain,
emset=emgain,
emgain=emgain_value,
csst_format=False,
shift=shift,
rotation=rotation,
camera=camera,
prograss_bar=prograss_bar,
output=output
)
def deduplicate_names_add_count(names: list):
for i in range(len(names)-1,-1,-1):
if names.count(names[i]) > 1:
names[i] = names[i] + '_' + str(names.count(names[i]))
def observation_simulation_from_config(obs_file, config_file):
config = {}
if config_file:
with open(config_file) as fid:
config = yaml.load(fid, Loader=yaml.FullLoader)
camera_config = config.get('camera', [{}])
all_camera = []
all_camera_name = []
for c_config in camera_config:
all_camera.append(CpicVisEmccd(c_config))
all_camera_name.append(c_config.get('name', 'camera'))
deduplicate_names_add_count(all_camera_name)
for key in config.keys():
if key in update_able_keys:
default_config[key] = config[key]
output_folder = default_config['output']
csst_format = default_config['csst_format']
nsample = default_config['nsample']
obs_file = os.path.abspath(obs_file)
file_list = []
if os.path.isdir(obs_file):
file_list = glob(f"{obs_file}/*.yaml")
elif '.yaml' == obs_file[-5:]:
file_list = [obs_file]
if not file_list:
log.warning(f"No observation file found in {obs_file}")
for ind_target, file in enumerate(file_list):
with open(file, 'r') as fid:
obs_info = yaml.load(fid, Loader=yaml.FullLoader)
target = target_file_load(obs_info.get('target', {}))
skybg = obs_info.get('skybg', None)
expt = obs_info['expt']
band = obs_info['band']
emset = obs_info['emset']
nframe = obs_info['nframe']
obsid = obs_info['obsid']
rotation = obs_info.get('rotation', 0)
shift = obs_info.get('shift', [0, 0])
gnc_info = obs_info.get('gnc_info', {})
time = obs_info.get('time', 0)
emgain = obs_info.get('emgain', None)
time = relative_time(time)
ind_camera = 0
for camera_name, camera in zip(all_camera_name, all_camera):
try:
ind_camera += 1
ind_run = ind_target * len(all_camera) + ind_camera
all_run = len(all_camera) * len(file_list)
info_text = f"({ind_run}/{all_run}) obsid[{obsid}]/{os.path.basename(file)[:-5]} with {camera_name}"
log.info(info_text)
if time == 0:
camera.time_syn(time, initial=True)
else:
dt = time - camera.system_time
if dt < 0:
log.warning(f'Time is not synced. {dt} seconds are added.')
dt = 0
camera.time_syn(dt, readout=False)
if len(all_camera) > 1:
output = os.path.join(output_folder, camera_name)
else:
output = output_folder
vis_observation(
target,
skybg,
expt,
nframe,
band,
emset,
obsid,
emgain=emgain,
rotation=rotation,
shift=shift,
gnc_info=gnc_info,
camera=camera,
output=output,
nsample=nsample,
csst_format=csst_format,
prograss_bar=True)
except Exception as e:
raise(e)
def main(argv=None):
parser = argparse.ArgumentParser(description='Cpic obsevation image simulation')
parser.set_defaults(func=lambda _: parser.print_usage())
subparsers = parser.add_subparsers(help='type of runs')
parser_quickrun = subparsers.add_parser('quickrun', help='a quick observation with no configration file')
parser_quickrun.add_argument('target_string', type=str, help='example: \*5.1/25.3(1.3,1.5)/22.1(2.3,-4.5)')
parser_quickrun.add_argument('expt', type=float, help='exposure time [ms]')
parser_quickrun.add_argument('emgain', type=float, help='emgain or emgain set value if emgain_input is False')
parser_quickrun.add_argument('nframe', type=int, help='number of frames')
parser_quickrun.add_argument('-b', '--band', type=str, default='f661', help='band, one of f565/f661/f743/f883')
parser_quickrun.add_argument('-r', '--rotation', type=float, default=0, help='rotation angle [degree]')
parser_quickrun.add_argument('-s', '--skybk', type=float, default=21, help='magnitude of sky background [mag/arcsec^2]')
parser_quickrun.add_argument('-f', '--cr_frame', action='store_true', help='if True, cosmic ray frame will be added')
parser_quickrun.add_argument('-e', '--emset', action='store_true', help='if True, emgain set value will be used as input')
parser_quickrun.add_argument('-c', '--camera_effect', action='store_true', help='if True, camera effect will be added')
parser_quickrun.add_argument('-o', '--output', type=str, default='./', help='output folder')
def quick_run_call(args):
quick_run_v2(
target_str=args.target_string,
expt=args.expt,
emgain=args.emgain,
nframe=args.nframe,
band=args.band,
rotation=args.rotation,
skybg=args.skybk,
cr_frame=args.cr_frame,
emset_input=args.emset,
camera_effect=args.camera_effect,
output=os.path.abspath(args.output),
prograss_bar=True
)
parser_quickrun.set_defaults(func=quick_run_call)
parser_run = subparsers.add_parser('run', help='observation with configration files or folder')
parser_run.add_argument('target_config', type=str, help='configration file or folder')
parser_run.add_argument('-c', '--config', type=str, help='instrument configration file')
def run_all(args):
observation_simulation_from_config(
args.target_config,
args.config,
)
parser_run.set_defaults(func=run_all)
args = parser.parse_args(argv)
args.func(args)
if __name__ == '__main__': # pragma: no cover
sys.exit(main())
# target_example = {
# 'cstar': {
# 'magnitude': 1,
# 'ra': '120d',
# 'dec': '40d',
# 'distance': 10,
# 'sptype': 'F0III',
# },
# 'stars': [
# {
# 'magnitude': 20,
# 'pangle': 60,
# 'separation': 1,
# 'sptype': 'F0III'
# }
# ]
# }
# # quick_run('', 10, 'f661', 1, 1, 30)
# # quick_run('*2.4/10(3,5)/15(-4,2)', 21, 'f661', 1, 1, 30)
# # # normal target
# observation_simulation(
# target={},
# skybg=15,
# expt=10,
# target=target_example,
# skybg=21,
# expt=1,
# nframe=2,
# band='f661',
# emgain=30,
# obsid=50112345678,
# obsid=51012345678,
# )
# quick_run('*5.1/25.3(0.8,0.8)', None, 'f661', 10, 1, 10)
# quick_run('*5/20(0.8,0.8)', None, 'f883', 10, 1, 10)
# # quick *5.1/25.3(1.3,1.5) expt nframe emgain band rotation shift
# # quick target_name expt nframe emgain band rotation shift
# # plan plan_file_or_folder
if __name__ == '__main__': # pragma: no cover
target_example = {
'cstar': {
'magnitude': 1,
'ra': '120d',
'dec': '40d',
'distance': 10,
'sptype': 'F0III',
},
'stars': [
{
'magnitude': 20,
'pangle': 60,
'separation': 1,
'sptype': 'F0III'
}
]
}
# quick_run('', 10, 'f661', 1, 1, 30)
# quick_run('*2.4/10(3,5)/15(-4,2)', 21, 'f661', 1, 1, 30)
# # normal target
observation_simulation(
target=target_example,
skybg=21,
expt=1,
nframe=2,
band='f661',
emgain=30,
obsid=51012345678,
)
# # bias
# observation_simulation(
# target=target_example,
......
csst_cpic_sim/optics.py
View file @ 95f81f91
import os
import yaml
import numpy as np
from .config import cpism_refdata, which_focalplane, S # S is synphot
from .config import optics_config
from scipy.signal import fftconvolve
from scipy.ndimage import rotate
from .config import config, S # S is synphot
from .utils import region_replace
from .io import log
from .psf_simulation import single_band_masked_psf, single_band_psf
FILTERS = {
"f565": S.FileBandpass(f"{cpism_refdata}/throughtput/f565_total.fits"),
"f661": S.FileBandpass(f"{cpism_refdata}/throughtput/f661_total.fits"),
"f743": S.FileBandpass(f"{cpism_refdata}/throughtput/f743_total.fits"),
"f883": S.FileBandpass(f"{cpism_refdata}/throughtput/f883_total.fits"),
"f940": S.FileBandpass(f"{cpism_refdata}/throughtput/f940_total.fits"),
"f1265": S.FileBandpass(f"{cpism_refdata}/throughtput/f1265_total.fits"),
"f1425": S.FileBandpass(f"{cpism_refdata}/throughtput/f1425_total.fits"),
"f1542": S.FileBandpass(f"{cpism_refdata}/throughtput/f1542_total.fits")
}
FILTERS = {}
for key, value in config['bands'].items():
FILTERS[key] = S.FileBandpass(value)
default_band = config['default_band']
def filter_throughput(filter_name):
"""
Totally throughput of each CPIC band.
Totally throughput of the each CPIC band.
Including the throughput of the filter, telescope, cpic, and camera QE.
If the filter_name is not supported, return the throughput of the default filter(f661).
......@@ -39,142 +33,132 @@ def filter_throughput(filter_name):
"""
filter_name = filter_name.lower()
filter_name = "f661" if filter_name == "default" else filter_name
filter_name = default_band if filter_name == 'default' else filter_name
if filter_name not in FILTERS.keys():
log.warning(f"滤光片名称错误({filter_name})，返回默认滤光片(f661)透过率")
filter_name = "f661"
log.warning(f"滤光片名称错误({filter_name})，返回默认滤光片({default_band})透过率")
filter_name = default_band
return FILTERS[filter_name]
def example_psf_func(band, spectrum, frame_size, error=0.1):
"""
Example psf generating function.
def rotate_and_shift(shift, rotation, init_shifts):
rotation_rad = rotation / 180 * np.pi
return np.array([
shift[0] * np.cos(rotation_rad) + shift[1] * np.sin(rotation_rad),
-shift[0] * np.sin(rotation_rad) + shift[1] * np.cos(rotation_rad)
]) + np.array(init_shifts)
Parameters
-------------
band: str
The name of the band.
spectrum: synphot.Spectrum or synphot.SourceSpectrum
The spectrum of the target.
frame_size: int
The size of the frame.
error: float
Phase RMS error.
Returns
---------------
2D array
psf image with shape of `frame_size`
def ideal_focus_image(
bandpass: S.spectrum.SpectralElement,
targets: list,
platescale,
platesize: list = [1024, 1024],
init_shifts: list = [0, 0],
rotation: float = 0,):
"""
pass
focal_image = np.zeros(platesize)
focal_shape = np.array(platesize)[::-1] # x, y
if not targets:
return focal_image
for target in targets:
sub_x, sub_y, sub_spectrum, sub_image = target
sub_shift = rotate_and_shift([sub_x, sub_y], rotation, init_shifts) / platescale
sed = (sub_spectrum * bandpass).integrate()
if sub_image is None:
x = (focal_shape[0] - 1)/2 + sub_shift[0]
y = (focal_shape[1] - 1)/2 + sub_shift[1]
int_x = int(x)
int_y = int(y)
if int_x < 0 or int_x >= focal_shape[0] - 1 or int_y < 0 or int_y >= focal_shape[1] - 1:
continue
dx1 = x - int_x
dx0 = 1 - dx1
dy1 = y - int_y
dy0 = 1 - dy1
sub = np.array([
[dx0*dy0, dx1*dy0],
[dx0*dy1, dx1*dy1]]) * sed
focal_image[int_y: int_y+2, int_x: int_x+2] += sub
else:
# sub_image = sub_image
sub_image = np.abs(rotate(sub_image, rotation, reshape=False))
sub_image = sub_image / sub_image.sum()
sub_img_shape = np.array(sub_image.shape)[::-1]
sub_shift += (focal_shape-1)/2 - (sub_img_shape-1)/2
focal_image = region_replace(
focal_image,
sub_image * sed,
sub_shift,
subpix=True
)
return focal_image
def make_focus_image(
def focal_convolve(
band: str,
targets: list,
psf_function: callable,
init_shifts: list = [0, 0],
rotation: float = 0,
platesize: list = [1024, 1024],
) -> np.ndarray:
"""
Make the focus image of the targets.
nsample: int = 5,
error: float = 0,
platesize: list = [1024, 1024]) -> np.ndarray :
Parameters
-----------
band: str
The name of the band.
targets: list
The list of the targets.
Each element of the list is a tuple of (x, y, spectrum).
- x, y: float
- The position of the target in the focal plane.
- spectrum: synphot.Spectrum or synphot.SourceSpectrum
- The spectrum of the target.
psf_function: callable
The function to generate the PSF, with same parameters and return as `example_psf_func`.
init_shifts: list
The initial shifts of the center targets. Unit: arcsec.
The default is [0, 0].
rotation: float
The rotation of the focal plane. Unit: degree.
The default is 0 degree.
platesize: list
The size of the focal plane. Unit: pixel.
The default is [1024, 1024].
# config = optics_config[which_focalplane(band)]
platescale = config['platescale']
Returns
--------
np.ndarray
The focus image of the targets.
2D array with the shape of platesize.
"""
# telescope_config = optics_config['telescope']
area = config['aperature_area']
config = optics_config[which_focalplane(band)]
platescale = config["platescale"]
filter = filter_throughput(band)
wave = filter.wave
throughput = filter.throughput
min_wave = wave[0]
max_wave = wave[-1]
focal_image = np.zeros(platesize)
if not targets:
return focal_image
platescale = config['platescale']
iwa = config['mask_width'] / 2
def rotate_and_shift(shift):
rotation_rad = rotation / 180 * np.pi
return np.array(
[
shift[0] * np.cos(rotation_rad) + shift[1] * np.sin(rotation_rad),
-shift[0] * np.sin(rotation_rad) + shift[1] * np.cos(rotation_rad),
]
) + np.array(init_shifts)
if abs(init_shifts[0]) > 4 or abs(init_shifts[1]) > 4:
print('Input shifts are too large, and are set to zero')
init_shifts = [0, 0]
cstar_x, cstar_y, cstar_spectrum = targets[0]
cstar_shift = rotate_and_shift([cstar_x, cstar_y]) / platescale
all_fp_image = []
if not targets:
return np.zeros((platesize[1], platesize[0]))
error_value = 0 # nm
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
center_wavelength = (wave0 + wave1) / 2 * 1e-10
cstar_psf = psf_function(
band, cstar_spectrum, config["cstar_frame_size"], error=error_value
)
i_throughput = throughput.copy()
i_throughput[(wave > wave1) | (wave < wave0)] = 0
i_band = S.ArrayBandpass(wave, i_throughput, waveunits='angstrom')
platesize = np.array(platesize)[::-1]
psf_shape = np.array(cstar_psf.shape)[::-1]
cstar_shift += (platesize - 1) / 2 - (psf_shape - 1) / 2
i_fp_image = ideal_focus_image(i_band, targets[1:], platescale, platesize, init_shifts, rotation)
psf = single_band_psf(center_wavelength, error=error)
focal_image = region_replace(
focal_image,
cstar_psf,
cstar_shift,
padded_in=False,
padded_out=False,
subpix=True,
)
_, _, cstar_sp, _ = targets[0]
cstar_flux = (cstar_sp * i_band).integrate()
cstar_psf = single_band_masked_psf(center_wavelength, error=error, shift=init_shifts)
for i_target in range(1, len(targets)):
sub_x, sub_y, sub_spectrum = targets[i_target]
pdout = False if i_target == len(targets) - 1 else True
pdin = False if i_target == 1 else True
log.debug(f"input target {sub_x=:}, {sub_y=:}")
sub_shift = rotate_and_shift([sub_x, sub_y]) / platescale
log.debug(f"after rotate and shift {sub_shift=:}")
sub_psf = psf_function(
band, sub_spectrum, config["substellar_frame_size"], error=error_value
)
psf_shape = np.array(sub_psf.shape)[::-1]
sub_shift += (platesize - 1) / 2 - (psf_shape - 1) / 2
log.debug(f"input shift of region_replace: {sub_shift=:}")
focal_image = region_replace(
focal_image,
sub_psf,
sub_shift,
padded_in=pdin,
padded_out=pdout,
subpix=True,
)
c_fp_image = fftconvolve(i_fp_image, psf, mode='same')
c_fp_image = focal_mask(c_fp_image, iwa, platescale)
c_fp_image = c_fp_image + cstar_flux * cstar_psf
return focal_image
all_fp_image.append(c_fp_image * area) # trans to photon/second
return np.array(all_fp_image).sum(axis=0)
def focal_mask(image, iwa, platescale, throughtput=1e-6):
......@@ -197,11 +181,10 @@ def focal_mask(image, iwa, platescale, throughtput=1e-6):
np.ndarray
The masked image.
"""
xx, yy = np.mgrid[0 : image.shape[0], 0 : image.shape[1]]
center = np.array([(image.shape[0] - 1) / 2, (image.shape[1] - 1) / 2])
mask = (abs(xx - center[0]) < iwa / platescale) | (
abs(yy - center[1]) < iwa / platescale
)
xx, yy = np.mgrid[0:image.shape[0], 0:image.shape[1]]
center = np.array([(image.shape[0]-1)/2, (image.shape[1]-1)/2])
mask = (abs(xx - center[0]) < iwa /
platescale) | (abs(yy - center[1]) < iwa / platescale)
image_out = image.copy()
image_out[mask] *= throughtput
return image_out
csst_cpic_sim/psf_simulation.py
View file @ 95f81f91
import os, pickle
import numpy as np
from astropy.io import fits
from hcipy import Field, Wavefront, DeformableMirror, FraunhoferPropagator
from hcipy import Field, Wavefront, 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 hcipy import TipTiltMirror, DeformableMirror
from hcipy import ComplexSurfaceApodizer
from hcipy import make_pupil_grid, make_circular_aperture, make_obstructed_circular_aperture
from hcipy import make_xinetics_influence_functions, make_uniform_grid
from .config import cpism_refdata, S
from .optics import filter_throughput
from .config import config
DMCACHE = True
ARCSEC2RAD = np.radians(1 / 3600)
# 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')
apm, apm_header = fits.getdata(config['apm_file'], header=True)
actuator = fits.getdata(config['actuator_file'])
surface_aberration = fits.getdata(config['aberration'])
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'])
F_number = 83
focal_length = pupil_diameter * F_number
pupil_shape = apm.shape
pupil_rate = apm_header['PHRATE'] # meter/pixel
pupil_grid = make_pupil_grid(pupil_shape[0], pupil_shape[0] * pupil_rate)
aperture = make_circular_aperture(pupil_diameter)(pupil_grid)
aperture = aperture * Field(apm.flatten(), pupil_grid)
second_pupil_size = pupil_diameter * 0.8 # just gauss a number
second_aperture = make_circular_aperture(second_pupil_size)(pupil_grid)
lyot_stop = ComplexSurfaceApodizer(second_aperture, second_aperture*0, lambda _: 2)
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)
platescale = emccd_pixel_size / (ARCSEC2RAD * focal_length) # arcsec / pixel
focal_shape = np.array([1024, 1024])
focal_size = focal_shape * emccd_pixel_size
# platescale = vis_focal_config['platescale']
mask_width = config['mask_width']
focal_full_frame = make_uniform_grid(focal_shape, focal_size, has_center=True)
prop_full_frame = FraunhoferPropagator(pupil_grid, focal_full_frame, focal_length)
spider_width = mask_width/platescale*emccd_pixel_size
spider = make_obstructed_circular_aperture(focal_size[0]*2, 0, 4, spider_width)(focal_full_frame)
pupil_cross_mask = ComplexSurfaceApodizer(spider, spider*0, lambda _: 2)
num_actuators_across = 32
# dm spacing is little smaller than pupil
actuator_spacing = 0.95 / 32 * pupil_diameter
influence_functions = make_xinetics_influence_functions(
actuator_spacing = 0.95 * pupil_diameter / num_actuators_across
pickle_file = config['dm_pickle']
if os.path.exists(pickle_file) and DMCACHE:
with open(pickle_file, 'rb') as fid:
influence_functions = pickle.load(fid)
else:
influence_functions = make_xinetics_influence_functions(
pupil_grid, num_actuators_across, actuator_spacing)
with open(pickle_file, 'wb') as fid:
pickle.dump(influence_functions, fid)
deformable_mirror = DeformableMirror(influence_functions)
tiptilt_mirror = TipTiltMirror(pupil_grid)
aberration = SurfaceApodizer(
surface_sag=surface_aberration.flatten(), refractive_index=-1)
......@@ -59,48 +73,31 @@ aberration = SurfaceApodizer(
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 = []
def single_band_masked_psf(wavelength, error=0, shift=[0, 0]):
error = np.random.normal(0, error*1e-9, actuator.shape)
deformable_mirror.actuators = actuator + error
wf = Wavefront(aperture, wavelength)
shift = np.array(shift) * ARCSEC2RAD / 2
tiptilt_mirror.actuators = shift
wf = tiptilt_mirror(wf)
wf = aberration(wf)
first_focal = prop_full_frame(deformable_mirror(wf))
prop = prop_full_frame
if not cstar:
prop = prop_sub_frame
strength = first_focal.intensity.shaped.sum()
second_pupil = prop_full_frame.backward(pupil_cross_mask(first_focal))
second_focal = prop_full_frame(lyot_stop(second_pupil))
psf = second_focal.intensity.shaped
return psf / strength
for i_sample in range(nsample):
wf = Wavefront(aperture, sed_center_wavelength[i_sample])
def single_band_psf(wavelength, error=0, aber_phase=None, intensity=True):
error = np.random.normal(0, error*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(deformable_mirror(wf)).intensity.shaped
imgs.append(img / img.sum() * sed[i_sample])
return np.array(imgs).sum(axis=0)
first_focal = prop_full_frame(deformable_mirror(wf))
image = np.array(first_focal.intensity.shaped)
return image / image.sum()
\ No newline at end of file
csst_cpic_sim/target.py
View file @ 95f81f91
import re
import json
import os, re, json, yaml
from typing import Union
import numpy as np
from scipy import constants
from astropy.io import fits
from astropy.coordinates import SkyCoord
from pysynphot.renorm import StdRenorm
from .config import cpism_refdata, MAG_SYSTEM
from .config import config
from .config import S # pysynphot
from .optics import filter_throughput
from .io import log
PLANK_H = constants.h # ~6.62607015e-34
LIGHT_SPEED = constants.c # ~3e8
DEG2RAD = np.pi / 180
R_JUPITER_METER = 6.99e4
AU_METER = 1.49e8
DEFAULT_FILTER_FILE = f'{cpism_refdata}/throughtput/f661_total.fits'
HYBRID_MODEL_FILE = f'{cpism_refdata}/target_model/hybrid_model.fits'
DEFAULT_FILTER_FILE = config['default_band']
HYBRID_MODEL_FILE = config['hybrid_model']
BCC_MODEL_FOLDER = config['bcc_model']
MAG_SYSTEM = config['mag_system']
DEFAULT_FILTER = S.ObsBandpass('johnson,v')
CATALOG_CACHE = {}
class AlbedoCat(S.spectrum.TabularSpectralElement, S.catalog.Icat):
"""Generate albedo spectrum from the planet reflection model in Batalha et al. 2018
Parameters
----------
phase : float
[degree], The phase angle of the planet, range from 0 to 180.
metallicity : float
The metallicity of the planet. log(Z/Zsun), range from 0 to 2.
f_sed : float
The sedimentation efficiency of the cloud. log(f_sed), range from -2 to 2. 2 is for cloud free case.
Notes
-----
f_sed is only reliable between -2 and log10(6) and the cloud free case (2).
Values between log10(6) and 2 are interpolated from the cloud free case (2) and log10(6).
Reference
---------
Color Classification of Extrasolar Giant Planets: Prospects and Cautions
Natasha E. Batalha et al 2018 AJ 156 158
def _sptype2num(spectral_type):
"""
def __init__(self,
phase: float,
metallicity: float,
f_sed: float,
):
catname = 'BCCalbedo'
self.isAnalytic=False
self.name = f"{catname}_{phase}_{metallicity}_{f_sed}"
self.parameter_names = ['phase', 'metallicity', 'f_sed']
self.catalog_folder = BCC_MODEL_FOLDER
filename = os.path.join(self.catalog_folder, 'catalog.fits')
if filename in CATALOG_CACHE:
indices = CATALOG_CACHE[filename]
else:
with fits.open(filename) as table:
indexList = table[1].data.field('INDEX')
filenameList = table[1].data.field('FILENAME')
indices = self._getArgs(indexList, filenameList)
CATALOG_CACHE[filename] = indices
list0,list1 = self._breakList(indices, 0, phase)
list2,list3 = self._breakList(list0, 1, metallicity)
list4,list5 = self._breakList(list1, 1, metallicity)
list6,list7 = self._breakList(list2, 2, f_sed)
list8,list9 = self._breakList(list3, 2, f_sed)
list10,list11 = self._breakList(list4, 2, f_sed)
list12,list13 = self._breakList(list5, 2, f_sed)
sp1, wave, waveunits = self._getSpectrum(list6[0], catname, wave_output=True)
sp2 = self._getSpectrum(list7[0], catname)
sp3 = self._getSpectrum(list8[0], catname)
sp4 = self._getSpectrum(list9[0], catname)
sp5 = self._getSpectrum(list10[0], catname)
sp6 = self._getSpectrum(list11[0], catname)
sp7 = self._getSpectrum(list12[0], catname)
sp8 = self._getSpectrum(list13[0], catname)
spa1 = self._interpolateSpectrum(sp1, sp2, f_sed)
spa2 = self._interpolateSpectrum(sp3, sp4, f_sed)
spa3 = self._interpolateSpectrum(sp5, sp6, f_sed)
spa4 = self._interpolateSpectrum(sp7, sp8, f_sed)
spa5 = self._interpolateSpectrum(spa1, spa2, metallicity)
spa6 = self._interpolateSpectrum(spa3, spa4, metallicity)
spa7 = self._interpolateSpectrum(spa5, spa6, phase)
sp = spa7[0]
self._wavetable = wave * 1e4
self._throughputtable = sp
self.waveunits = S.units.Units(waveunits.lower())
self.warnings = {}
def _getSpectrum(self, parlist, catdir, wave_output=False):
name = parlist[3]
filename = name.split('[')[0]
column = name.split('[')[1][:-1]
filename = f"{self.catalog_folder}/{filename}"
# sp = S.spectrum.TabularSpectralElement(filename, thrucol=column)
with fits.open(filename) as td:
sp = td[1].data.field(column)
wave = td[1].data.field('wavelength')
# waveunits = td[1].header['tunit1']
waveunits = 'micron'
result = []
for member in parlist:
result.append(member)
result.pop()
result.append(sp)
if wave_output:
return result, wave, waveunits
return result
def _sptype2num(
spectral_type: str) -> tuple:
"""
convert spectral type string to number, for interpretation
......@@ -73,12 +190,14 @@ def _sptype2num(spectral_type):
return shorttype + subtype / 10, startype
def _spnum2teff(spn, stn):
def _spnum2teff(
spn: int,
stn: int) -> tuple:
"""
interpret of spectral number (by __sptype2num) to get t_eff and log_g
look up table from the document of ck04model
"""
with open(f'{cpism_refdata}/target_model/sptype2teff_lut.json', 'r') as fid:
with open(config['sp2teff_model'], 'r') as fid:
sptype_teff_lut = json.load(fid)
def _interp(spn, stn):
......@@ -101,7 +220,7 @@ def star_photlam(
magnitude: float,
sptype: str,
is_blackbody: bool = False,
mag_input_band: str = 'f661'):
mag_input_band: str = 'f661') -> S.ArraySpectrum:
"""
genrate flux spectrum of a star by its observal magnitude and spectral type
......@@ -137,12 +256,59 @@ def star_photlam(
star_sp.convert('photlam')
return star_sp
def standard_spectrum(magnitude: float):
"""Standard spectrum of magnitude system.
Parameters
-----------
magnitude : float
magnitude of the standard spectrum
Returns
----------
star_sp : S.spectrum.Spectrum
"""
inner_std = S.units.Units(MAG_SYSTEM).StdSpectrum
std_spectrum = S.ArraySpectrum(
inner_std.wave,
inner_std.flux,
inner_std.waveunits,
inner_std.fluxunits
)
filter = filter_throughput(DEFAULT_FILTER_FILE)
std_spectrum = StdRenorm(std_spectrum, filter, magnitude, MAG_SYSTEM)
std_spectrum.convert('photlam')
return std_spectrum
def bcc_spectrum(
coe_cloud: float,
coe_metalicity: float
):
"""Albedo spectrum of planet using BCC model (Batalha et al. 2018),
Parameters
----------
coe_cloud: float
The sedimentation efficiency of the cloud. log(f_sed), range from -2 to 2. 2 is for cloud free case.
coe_metalicity: float
The metallicity of the planet. log(Z/Zsun), range from 0 to 2.
Returns
-------
pysynphot.spectrum.ArrayBandpass
albedo spectrum of the planet
"""
spec = AlbedoCat(0, coe_metalicity, coe_cloud)
spec.convert('angstrom')
return spec
def hybrid_albedo_spectrum(
coe_b: float,
coe_r: float):
"""
albedo spectrum of planet using hybrid-jupiter-neptune model (Lacy et al. 2018)
"""Albedo spectrum of planet using hybrid-jupiter-neptune model (Lacy et al. 2018)
jupiter and neptune spectrum is from Karkoschka’s 1994
Parameters
......@@ -233,6 +399,185 @@ def extract_target_x_y(
return x, y
def detect_template_path(template: str) -> str:
"""Find template file in catalog folder or current folder.
Parameters
----------
template: str
template file name
Returns
-------
str
absolute path of template file
"""
if os.path.isfile(template):
return os.path.abspath(template)
catalog = config['catalog_folder']
cat_temp = os.path.join(catalog, template)
if os.path.isfile(cat_temp):
return cat_temp
raise FileExistsError(f'cant find {template} in ./ or catalog folder')
class TargetOjbect(object):
"""A helper class to generate target spectrum and albedo spectrum
Attributes
----------
x: float
x of target in arcsec
y: fload
y of target in arcsec
ra: str
ra string for center star, such as '15d23m05s'
dec: str
dec string for center star
distance: float
distance of center star in pc
image: 2d np.array
image of the target
spectrum: pysynphot.spectrum.Spectrum
spectrum of the target
"""
def __init__(self, info, cstar=None):
"""Initialize a target object
Parameters
----------
info: dict
target info, see Example for more details
cstar: TargetOjbect or None
center star object bounded, if None, means the target is the center star itself
center star object is used to calculate the x, y of each target according to its ra and dec
also center star's distance is used to calculate seperation of planet
sp_model: str
keyword to compatible with V1.0 input. See Notes
spectral type of the target, can be 'star' or 'planet'
Notes
-----
keyargv 中的sp_model 完全用来用来兼容V1.0版的输入,执行逻辑如下:
1. 如果没有设置sp_model, 对应了V2.0的输入, 则sp_model从info中获取
2. 如果设置了sp_model, 则区分为star和planet
3. 如果sp_model为star, 则判断是否设置了isblackbody, 如果是的话, 则sp_model改为blackbody
4. 如果sp_model为planet, 则进一步区分为hybrid_planet 和 bcc_planet 和 template_planet
Examples:
--------
cstar = {
'magnitude': 0,
'ra': '120d',
'dec': '40d',
'distance': 10,
'sptype': 'G0III'
}
stars = {
'magnitude': 15,
'ra': '120.001d',
'dec': '40.001d',
'sptype': 'F0III',
'sp_spectrum': 'blackbody',
}
planets = {
'radius': 2,
'pangle': 60,
'coe_cloud': 0.3,
'coe_metal': 0.7,
'separation': 0.5,
'phase_angle': 90,
'sp_spectrum': 'hybrid_planet',
'image': 'extend_planet.fits'
}
# planet using input costum albedo spectrum!
# Note that the albedo spectrum is not normalized!
# so the contrast of the planet sould be considered in the input file by user!
# The input file is in pysynphot.FileSpectralElement format.
# See the documentation of pysynphot for details.
planets = {
'pangle': 60,
'separation': 0.5,
'sp_spectrum': 'template_planet',
'template': 'planet_albedo.fits'
}
"""
self.sp_model = info.get('sp_model', 'star')
if cstar is None:
self.x, self.y = 0, 0
self.ra = info['ra']
self.dec = info['dec']
self.distance = info.get('distance', None)
else:
self.x, self.y = extract_target_x_y(info, cstar.ra, cstar.dec)
self.image = None
if 'image' in info.keys():
self.image = fits.getdata(detect_template_path(info['image']))
if self.sp_model == 'blackbody':
self.spectrum = star_photlam(
info['magnitude'],
info['sptype'],
is_blackbody=True,
mag_input_band=info.get('mag_input_band', 'f661')
)
if self.sp_model == 'reference':
self.spectrum = standard_spectrum(info['magnitude'])
if self.sp_model == 'template_star':
self.spectrum = S.FileSpectrum(detect_template_path(info['template']))
if self.sp_model == 'star':
self.spectrum = star_photlam(
info['magnitude'],
info['sptype'],
is_blackbody=False,
mag_input_band=info.get('mag_input_band', 'f661')
)
if self.sp_model in ['hybrid_planet', 'bcc_planet', 'template_planet']:
planet = info
phase_angle = planet.get('phase_angle', 90)
sp_model = self.sp_model
radius = planet.get('radius', 1)
if cstar.distance is None:
raise ValueError('distance of center star must be provided if planet is added')
if planet.get('contrast', None) is not None:
contrast = planet['contrast']
else:
contrast = planet_contrast(
self.x * cstar.distance,
self.y * cstar.distance,
phase_angle,
radius,
)
if sp_model == 'hybrid_planet':
coe_blue, coe_red = planet.get('coe_b', 1), planet.get('coe_r', 1)
albedo_spect = hybrid_albedo_spectrum(coe_blue, coe_red)
elif sp_model == 'bcc_planet':
coe_c, coe_m = planet.get('coe_cloud', 2), planet.get('coe_metal', 0)
albedo_spect = bcc_spectrum(coe_c, coe_m)
else: #sp_model == 'template_planet'
albedo_spect = S.FileBandpass(detect_template_path(planet['template']))
contrast = 1
self.spectrum = cstar.spectrum * albedo_spect * contrast
def planet_contrast(
planet_x_au: float,
planet_y_au: float,
......@@ -274,7 +619,6 @@ def planet_contrast(
contrast = (radius / sep_3d * R_JUPITER_METER / AU_METER)**2 * phase
return contrast
def spectrum_generator(
targets: dict):
"""
......@@ -285,111 +629,122 @@ def spectrum_generator(
targets: dict
target dictionary which contains keyword 'cstar' (necessary), 'stars'(optional), 'planets' (optional).
The values are:
- cstar: dict
- center star information. must contain keywords ra, dec, distance, magnitude, sptype
- stars: list of dict, optional
- list of background stars. each dict must contain ra, dec, magnitude, sptype
- planets: list of dict, optional
- list of planets. each dict must contain pangle, separation, magnitude, radius
- objects: list of dict, optional, recommended! V2.0 new!
- list of targets. each dict must contain ra, dec, magnitude, sptype
Returns
-------
obj_sp_list: list
list of [x, y, spectrum] of each target
Examples
--------
>>> targets = {
... target = {
... 'name': 'TEST1',
... 'cstar': {
... 'magnitude': 0,
... 'ra': '120d',
... 'dec': '40d',
... 'distance': 10,
... 'sptype': 'G0III',
... },
... 'stars': [
... {
... 'magnitude': 15,
... 'ra': '120.001d',
... 'dec': '40.001d',
... 'sptype': 'F0III',
... 'isblackbody': True
... },
... {
... 'magnitude': 10,
... 'pangle': -60,
... 'separation': 2,
... 'sptype': 'A2.5II',
... 'isblackbody': False
... },
...
... ],
... 'planets': [
... {
... 'radius': 2,
... 'pangle': 60,
... 'coe_b': 0.3,
... 'coe_r': 0.7,
... 'separation': 0.5,
... 'phase_angle': 90,
... }
... ]
... }
>>> spectrum_generator(targets)
list of [x, y, spectrum, image] of each target
"""
cstar = targets['cstar']
stars = targets.get('stars', [])
planets = targets.get('planets', [])
objects = targets.get('objects', [])
obj_sp_list = []
cstar_spectrum = star_photlam(
cstar['magnitude'],
cstar['sptype'],
is_blackbody=cstar.get('isblackbody', False),
mag_input_band=cstar.get('mag_input_band', 'f661')
)
cstar_obj = TargetOjbect(cstar)
obj_sp_list.append([cstar_obj.x, cstar_obj.y, cstar_obj.spectrum, cstar_obj.image])
cstar_ra = cstar['ra']
cstar_dec = cstar['dec']
cstar_dist = cstar['distance']
obj_sp_list.append([0, 0, cstar_spectrum])
for target in objects:
target_obj = TargetOjbect(target, cstar=cstar_obj)
obj_sp_list.append([target_obj.x, target_obj.y, target_obj.spectrum, target_obj.image])
for star in stars:
star_x, star_y = extract_target_x_y(star, cstar_ra, cstar_dec)
return obj_sp_list
spectrum = star_photlam(
star['magnitude'],
star['sptype'],
is_blackbody=star.get('isblackbody', False),
mag_input_band=cstar.get('mag_input_band', 'f661')
)
obj_sp_list.append([star_x, star_y, spectrum])
def target_file_load(
target: Union[dict, str]) -> dict:
"""Generate target dict from file, string or dict.
for planet in planets:
planet_x, planet_y = extract_target_x_y(planet)
Parameters
----------
target: Union[dict, str]
target file name, formated string or target dict.
# angle between observation direction and star-planet direction
phase_angle = planet.get('phase_angle', 90)
Outputs
--------
target: dict
dictionary of target. Format of input
radius = planet.get('radius', 1)
Note
-------
If target is a string start with *, it will be treated as a formatted string.
e.g. "*5.1/25.3(1.3,1.5)/22.1(2.3,-4.5)" which means a central object
with magnitude 5.1, and two substellar with magnitude 25.3 and 22.1, respectively.
The spectrum of each object is standard refernece spectrum of ABmag.
The first number in the parenthesis is the x position in arcsec, and the second is the y position.
contrast = planet_contrast(
planet_x * cstar_dist,
planet_y * cstar_dist,
phase_angle,
radius,
)
If target is a string without *, it will be treated as file name. And find the file in catalog folder.
The file need to be in yaml format.
And end with .yaml (note .yml not work). If not .yaml will be added.
coe_blue, coe_red = planet.get('coe_b', 1), planet.get('coe_r', 1)
albedo_spect = hybrid_albedo_spectrum(coe_blue, coe_red)
If target is a dict, it will be returned directly.
spectrum = albedo_spect * cstar_spectrum * contrast
obj_sp_list.append([planet_x, planet_y, spectrum])
If all the above conditions are not met, an empty dict will be returned.
"""
return obj_sp_list
if not target: # None or empty string or {}
return {}
if isinstance(target, dict):
return target
if isinstance(target, str): #filename or formatted string
target = target.strip()
if not target:
return {}
catalog_folder = config['catalog_folder']
target_file = target
target_file += '.yaml' if target_file[-5:].lower() != '.yaml' else ""
target_full_path = os.path.join(catalog_folder, target_file)
if os.path.isfile(target_full_path):
with open(target_full_path) as fid:
target = yaml.load(fid, Loader=yaml.FullLoader)
target['name'] = target_file[:-5]
return target
target_str = target
if (target_str[0] == '*'):
objects = target_str[1:].split('/')
cstar_mag = float(objects[0])
cstar = {
'magnitude': cstar_mag,
'ra': '0d',
'dec': '0d',
'sp_model': 'reference',
'distance': 10,
}
stars = []
for sub_stellar in objects[1:]:
float_regex = R"[+-]?\d+(?:\.\d+)?"
match = re.match(
rf"({float_regex})\(({float_regex}),({float_regex})\)", sub_stellar)
if not match:
log.error(f'Wrong format for sub stellar: {sub_stellar}, Skip it')
continue
mag = float(match.group(1))
x = float(match.group(2))
y = float(match.group(3))
pangle = np.arctan2(x, y) * 180 / np.pi
separation = np.sqrt(x**2 + y**2)
stars.append({
'magnitude': mag,
'pangle': pangle,
'separation': separation,
'sp_model': 'reference',
})
target_dict = {
'name': target_str[1:],
'cstar': cstar,
'objects': stars,
}
return target_dict
log.error(f'Wrong format for target: {target}, using blank target instead')
return {}
csst_cpic_sim/utils.py
View file @ 95f81f91
......@@ -2,6 +2,7 @@ import numpy as np
import scipy.ndimage as nd
import logging
import random
import matplotlib.pyplot as plt
# DO NOT IMPORT CPICIMGSIM MODULES HERE
......@@ -125,3 +126,16 @@ def region_replace(
return padded
return padded[f_sz[0]:b_sz[0]+f_sz[0], f_sz[1]:b_sz[1]+f_sz[1]]
def psf_imshow(psf, vmin=1e-8, vmax=0.1, log=True, region=1):
focal_img = psf.copy()
focal_img = (focal_img - focal_img.min()) / (focal_img.max() - focal_img.min())
if log:
focal_img = np.log10(focal_img * 9 + 1)
plt.imshow(focal_img, origin='lower', cmap='gray', vmin=vmin, vmax=vmax)
shape = psf.shape
plt.xlim(shape[1] * (1 - region) / 2, shape[1] * (1 + region) / 2)
plt.ylim(shape[0] * (1 - region) / 2, shape[0] * (1 + region) / 2)
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment