Commit 190f1129 authored by GZhao's avatar GZhao
Browse files

pep 8 check.

parent b6e8bc32
Pipeline #7115 failed with stage
in 0 seconds
......@@ -15,4 +15,4 @@ __all__ = [
"quick_run_v2",
"vis_observation",
"__version__"
]
\ No newline at end of file
]
......@@ -2,7 +2,6 @@ import math
import numpy as np
import scipy.ndimage as nd
from astropy.io import fits
import matplotlib.pyplot as plt
from .config import config, S
from .utils import region_replace, random_seed_select
......@@ -13,6 +12,8 @@ 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.
......@@ -22,7 +23,7 @@ def sky_frame_maker(band, skybg, platescale, shape):
band : str
The band of the sky background.
skybg : str
The sky background file name.
The sky background file name.
platescale : float
The platescale of the camera in arcsec/pixel.
shape : tuple
......@@ -288,7 +289,7 @@ class CpicVisEmccd(object):
"""
def __init__(self, config_dict=None):
"""Initialize the camera.
Parameters
----------
config_dict : dict, optional
......@@ -301,7 +302,7 @@ class CpicVisEmccd(object):
self._defaut_config()
if config_dict is not None:
old_switch = self.switch
self.__dict__.update(config_dict) #not safe, be careful
self.__dict__.update(config_dict) # not safe, be careful
old_switch.update(self.switch)
self.switch = old_switch
self.config_init()
......@@ -345,17 +346,17 @@ class CpicVisEmccd(object):
self.cic = None
self.dark = None
self.flat = None
self.fullwell = 80_000
self.em_fullwell = 500_000 #780_000
self.em_fullwell = 500_000 # 780_000
self.em_cte = 0.9996
self.emreg_cal_num = 10 # 用来加速计算
self.emreg_cal_num = 10 # 用来加速计算
self.emreg_num = 604
self.readout_speed = 6.25e6 # Hz
self.readout_time = 0.365 # s
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.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
......@@ -372,14 +373,12 @@ class CpicVisEmccd(object):
self.nonlinear_coefficient = -0.1
self.detector_name = 'EMCCD'
self.ccd_label= 'CCD201-20'
self.detector_name = 'CCD201-20-EM'
self.ccd_label = 'EMCCD'
self.pitch_size = 13
def config_init(self):
"""initialize the camera.
"""initialize the camera.
If the config is set, call this function to update the config.
"""
......@@ -418,7 +417,7 @@ class CpicVisEmccd(object):
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
......@@ -431,7 +430,6 @@ class CpicVisEmccd(object):
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):
"""Calculate the emgain fix coeficient to fix the gamma distribution.
The coeficient is from fixing of ideal emgain distribution.
......@@ -440,7 +438,7 @@ class CpicVisEmccd(object):
----------
emgain : float
The emgain.
Returns
-------
float
......@@ -448,6 +446,7 @@ class CpicVisEmccd(object):
"""
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)
......@@ -465,10 +464,9 @@ class CpicVisEmccd(object):
output.append(kernel(em))
return np.array(output)
def bias_frame(self):
"""Generate bias frame
The bias frame contains vertical, horizontal, peper-salt noise, bias drift effect.
The bias frame contains vertical, horizontal, peper-salt noise, bias drift effect.
Can be configurable using self.switch.
Returns
......@@ -479,7 +477,7 @@ class CpicVisEmccd(object):
shape = self.bias_shape
TPI = np.pi * 2
# vertical pattern
# 使用一维的曲线描述竖条纹的截面
vp_1d = np.zeros(shape[1])
......@@ -505,11 +503,11 @@ class CpicVisEmccd(object):
# horizontal pattern
# 图像上的横条纹是梳状,分为两个部分,左边大约77个像素是周期小一点,其余的会大一点
boundary = 77 # boundary between left and width
boundary_width = 5 # 左右需要平滑过度一下
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_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
......@@ -518,10 +516,10 @@ class CpicVisEmccd(object):
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
......@@ -565,22 +563,22 @@ class CpicVisEmccd(object):
# plt.title('vertical pattern')
# 接上制冷机后,会有亮暗点
#cooler interfence effect
# 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_interval = 250 # 6.25MHz readout / 2.5KHz interfence
ci_dn = self.cooler_interfence
npix = shape[0] * shape[1]
n_ci_event = npix // ci_interval
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)
......@@ -589,7 +587,7 @@ class CpicVisEmccd(object):
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)
......@@ -607,7 +605,7 @@ class CpicVisEmccd(object):
bias_frame += rn_adu * np.random.randn(shape[0], shape[1]) + bias_shift
return bias_frame
def nonlinear_effect(self, image):
"""
nonlinear effect
......@@ -619,11 +617,11 @@ class CpicVisEmccd(object):
image += (image / fullwell)**2 * nonlinear_coefficient * fullwell
return image
def time_syn(self, t, readout=True, initial=False):
"""
time synchronization and update the system time and ccd temperature
Parameters
----------
t : float
......@@ -633,26 +631,26 @@ class CpicVisEmccd(object):
if True, the ccd temperature will increase, otherwise, it will decrease
initial : bool, optional
If inital is True, the ccd will be intialized to the cooler temperature
"""
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: #
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)
......@@ -672,14 +670,15 @@ class CpicVisEmccd(object):
if True, update the emgain and emset. Default is True.
if False, only return the emgain.
"""
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
def volt_func(es):
return volt_coe_a * es + volt_coe_b
self.volt = volt_func(em_set)
......@@ -692,16 +691,16 @@ class CpicVisEmccd(object):
# 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
self.emgain = emgain
self.emset = em_set
return emgain
def vertical_blooming(self, image):
"""
vertical blooming effect
......@@ -764,7 +763,7 @@ class CpicVisEmccd(object):
def readout(self, image_focal, em_set, expt_set, image_cosmic_ray=False, emgain=None):
"""From focal planet image to ccd output.
Interface function for emccd. Simulate the readout process.
Parameters
----------
image_focal : np.ndarray
......@@ -787,7 +786,7 @@ class CpicVisEmccd(object):
expt = expt_set
if expt_set == 0:
expt = 0.001
dt = self.readout_time + expt
self.time_syn(dt, readout=True)
......@@ -832,7 +831,7 @@ class CpicVisEmccd(object):
image[deadpix_y:, deadpix_x] = 0
img_bias = np.zeros(self.bias_shape, dtype=int)
img_bias[
self.pscan2: -self.oscan2,
self.pscan1: -self.oscan1
......@@ -845,8 +844,7 @@ class CpicVisEmccd(object):
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['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)
......@@ -857,11 +855,10 @@ class CpicVisEmccd(object):
# >>> for _ in range(self.emreg_num):
# >>> image += np.random.binomial(image, pEM)
# This code is too slow, so we used a modified gamma
em_fix = self.em_fix_fuc_fit(emgain) * emgain
image = np.random.gamma(image, em_fix) + image * (emgain - em_fix)
if self.switch['em_blooming']:
image = self.emregester_blooming(image)
......
import os, yaml
import os
import yaml
import warnings
from datetime import datetime
import numpy as np
......@@ -16,7 +17,7 @@ config_aim = os.path.join(config_aim, 'data/refdata_path.yaml')
# with open(config_aim, 'w') as f:
# yaml.dump(refdata_path, f)
# return refdata_path
# try:
# with open(config_aim, 'r') as f:
......@@ -37,8 +38,6 @@ def load_refdata_path(config_aim):
----------
config_aim : str
config_aim file path
"""
with open(config_aim, 'r') as f:
refdata_list = yaml.load(f, Loader=yaml.FullLoader)
......@@ -46,7 +45,7 @@ def load_refdata_path(config_aim):
for refdata in refdata_list:
if os.path.isdir(refdata):
return refdata
print("csst_cpic_sim refdata folder not found, please input cpism refencence data folder")
refdata = input()
refdata = os.path.abspath(refdata)
......@@ -55,7 +54,7 @@ def load_refdata_path(config_aim):
with open(config_aim, 'w') as f:
yaml.dump(refdata_list, f)
exit()
cpism_refdata = load_refdata_path(config_aim)
......@@ -80,13 +79,13 @@ config['bands'] = {
'f850': f'{cpism_refdata}/throughtput/f850.fits',
'f720': f'{cpism_refdata}/throughtput/f720.fits',
}
config['diameter'] = 2 # in meters
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['log_level'] = 'info'
config['output'] = './'
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'
......@@ -95,22 +94,25 @@ 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'
'apm_file', 'actuator_file', 'aberration',
'log_dir', 'log_level', 'catalog_folder',
'nsample', 'csst_format', 'output', 'check_fits_header'
]
def replace_cpism_refdata(
config: dict,
config: dict,
output: str = '$') -> None:
"""Replace the cpism_refdata in the config.
"""Replace the cpism_refdata in the config.
In the config file, we use ${cpism_refdata} to indicate the cpism_refdata.
This function is used to replace the cpism_refdata in the config, or replace back.
Parameters
----------
config: dict
config dict.
config dict.
output: str
'$' or 'other'. If output is '$', then replace the cpism_refdata in the config with ${cpism_refdata}.
'$' or 'other'. If output is '$', then replace the cpism_refdata in the config with ${cpism_refdata}.
If output is 'other', then replace the ${cpism_refdata} in the config file with the real path.
'$' is used meanly to generate a demo config file.
"""
......@@ -138,17 +140,18 @@ __version__ = '2.0.0'
with warnings.catch_warnings(): # pragma: no cover
warnings.filterwarnings("ignore")
import pysynphot as S
_ = S # S will be used in other modules, but I need to use it once here to pass the lint.
def setup_config(new_config):
"""Set up config from a dict.
"""Set up config from a dict.
Some of the configs need to calcuate.
Parameters
----------
new_config: dict
new config dict.
new config dict.
Returns
-------
None
......@@ -156,12 +159,14 @@ 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['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']]
setup_config({})
def which_focalplane(band):
"""
Return the name of the focalplane which the band belongs to.
......@@ -171,7 +176,6 @@ def which_focalplane(band):
-----------
band: str
The name of the band.
Returns
--------
......@@ -194,30 +198,32 @@ def which_focalplane(band):
return 'vis'
# raise ValueError(f"未知的波段{band}")
def iso_time(time):
"""Transfer relative time to iso time format
Parameters
----------
time: str or float
The relative time in seconds.
Returns
-------
str
The iso time format.
"""
if isinstance(time, str):
_ = datetime.fromisoformat(time)
_ = datetime.fromisoformat(time) # check if it is a iso time
return time
utc0 = config['utc0']
time0 = datetime.timestamp(datetime.fromisoformat(utc0))
time = datetime.fromtimestamp(time0 + time)
return time.isoformat()
def relative_time(time):
"""Transfer iso time format to relative time in seconds
......@@ -225,7 +231,7 @@ def relative_time(time):
----------
time: str or float
The iso time format.
Returns
-------
float
......@@ -236,8 +242,7 @@ def relative_time(time):
return time
if isinstance(time, int):
return float(time)
utc0 = config['utc0']
time0 = datetime.timestamp(datetime.fromisoformat(utc0))
return datetime.timestamp(datetime.fromisoformat(time)) - time0
\ No newline at end of file
import yaml, os, re
import yaml
import os
from datetime import datetime
import numpy as np
import pandas as pd
from astropy.io import fits
from astropy.coordinates import SkyCoord
import astropy.units as u
from .config import __version__, which_focalplane
from .config import __version__, which_focalplane
from .utils import Logger
from .config import config, iso_time
......@@ -14,13 +16,16 @@ default_output_dir = config['output']
log_level = config['log_level']
header_check = config['check_fits_header']
def set_up_logger(log_dir):
if not os.path.exists(log_dir):
os.makedirs(log_dir)
return Logger(log_dir+'/cpism_pack.log', log_level).logger
log = set_up_logger(config['log_dir'])
def check_and_update_fits_header(header):
"""
Check the header keywords and update the description according to the data model.
......@@ -43,7 +48,7 @@ def check_and_update_fits_header(header):
data_model = yaml.load(fid, Loader=yaml.FullLoader)
if 'FILETYPE' in header.keys():
header_model =data_model['HDU0']
header_model = data_model['HDU0']
hdu = 'hdu0'
else:
header_model = data_model['HDU1']
......@@ -88,9 +93,14 @@ def check_and_update_fits_header(header):
else:
key_type = 'ukn'
# print(f"keyword: {keyword} type: {key_type}, datamodel: {dtype}, value: {value}")
if key_type != dtype[0:3]:
print_warning(
f"Keyword {keyword} has wrong type in [{hdu}]. {dtype} expected, {key_type} found.")
if key_type == 'int' and dtype[0:3] == 'flo':
header[keyword] = float(header[keyword])
# print('transfer int to float with header keyword:', keyword)
else:
print_warning(
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():
......@@ -98,7 +108,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.")
elif keyword != 'COMMENT': # comment keyword is not allowed to be updated
elif keyword != 'COMMENT': # comment keyword is not allowed to be updated
header[keyword] = (header[keyword], dm_comment[keyword])
return header
......@@ -146,6 +156,27 @@ def obsid_parser(
return obstype
def datetime_obj_to_iso(time_obj):
"""
transfer datetime object to iso format used in csst fits file
example '2014-02-10T12:32:12.4'
Parameters
----------
time_obj: datetime.datetime
The datetime object.
Returns
-------
str
The iso format of the datetime object.
"""
isotime = time_obj.isoformat(sep='T', timespec='milliseconds')
subsec = str(round(float(isotime[19:22]), 1))
return isotime[:18] + subsec
def datetime_obj_to_mjd(time_obj):
"""
transfer datetime object to mean julian date (MJD).
......@@ -225,7 +256,7 @@ def primary_hdu(
header['NEXTEND'] = 1 # + parameters['nframe']
# header['GROUPS'] = False
header['DATE'] = datetime.now().isoformat(timespec='seconds')
header['DATE'] = datetime_obj_to_iso(datetime.now())
heaer_filename = filename[:-4]
if len(heaer_filename) > 68:
heaer_filename = heaer_filename[:68]
......@@ -243,8 +274,10 @@ def primary_hdu(
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)
ra_str = radec.ra.to_string(unit=u.hour, sep='', precision=1, pad=True)
dec_str = radec.dec.to_string(unit=u.deg, sep='', alwayssign=True, pad=True)
target_name = ra_str + dec_str
header['OBJECT'] = cstar.get('name', target_name)
header['TARGET'] = target_name
header['OBSID'] = str(obsid)
......@@ -253,7 +286,7 @@ def primary_hdu(
# telescope information
header['REFFRAME'] = 'CSSTGSC-1.0'
header['DATE-OBS'] = exp_start.isoformat(timespec='seconds')
header['DATE-OBS'] = datetime_obj_to_iso(exp_start)
header['SATESWV'] = '1'
header['EXPSTART'] = datetime_obj_to_mjd(exp_start)
......@@ -299,10 +332,10 @@ def primary_hdu(
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['EPOCH'] = float(exp_start.year)
header['EXPTIME'] = (exp_end - exp_start).total_seconds()
header['CHECKSUM'] = '0000000000000000'
header['DATASUM'] = '0000000000'
......@@ -368,12 +401,12 @@ def frame_header(obs_info, index, primary_header, camera_dict):
header['BUNIT'] = 'ADU'
header['FILTER'] = obs_info['band']
header['DETSN'] = '0'
header['DETSN'] = '00000000000'
header['DETNAME'] = camera_config['detector_name']
header['CHIPLAB'] = camera_config['ccd_label']
header['DEWTEMP'] = float(camera_config['cooler_temp'])
header['DEWTEMP'] = float(camera_config['cooler_temp']) + 273.15
frame_info = obs_info['frame_info'][index]
header['CHIPTEMP'] = float(frame_info['chiptemp'])
header['CHIPTEMP'] = float(frame_info['chiptemp']) + 273.15
header['DETSIZE'] = f"{imgszx} * {imgszy}"
header['IMGINDEX'] = index + 1
......@@ -513,7 +546,6 @@ def save_fits_simple(images, obs_info, output_folder='./'):
shift = obs_info['shift']
header['shift'] = f"x:{shift[0]},y:{shift[1]}"
fullname = os.path.join(output_folder, filename)
print(fullname)
if not os.path.exists(output_folder):
......
import argparse, sys, tqdm, time, os, yaml
import argparse
import tqdm
import time
import os
import yaml
from glob import glob
from datetime import datetime
import traceback
......@@ -25,8 +30,8 @@ def vis_observation(
shift: list = [0, 0],
gnc_info: dict = {},
csst_format: bool = True,
camera = CpicVisEmccd(),
crmaker = CosmicRayFrameMaker(),
camera=CpicVisEmccd(),
crmaker=CosmicRayFrameMaker(),
nsample: int = 1,
emgain=None,
prograss_bar=None,
......@@ -49,7 +54,8 @@ def vis_observation(
emset: int
EM gain setting value. 1023(0x3FF) for ~1.0× EM gain.
obsid: int
observation ID. Start from 4 for CPIC, 01 for science observation. See the input of io.obsid_parser for more details.
observation ID. Start from 4 for CPIC, 01 for science observation.
See the input of io.obsid_parser for more details.
rotation: float
rotation of the telescope. in unit of degree. 0 means North is up.
shift: list
......@@ -174,6 +180,7 @@ def vis_observation(
print(f'\r Done [{time.time() - start_time:.1f}s] ')
return image_cube
def quick_run_v2(
target_str: str,
band: str,
......@@ -183,12 +190,12 @@ def quick_run_v2(
skybg: float = None,
rotation: float = 0,
shift: list = [0, 0],
emset_input: bool=False,
cr_frame: bool=True,
camera_effect: bool=True,
emset_input: bool = False,
cr_frame: bool = True,
camera_effect: bool = True,
prograss_bar=False,
output='./') -> np.ndarray:
"""A quick run function for CPIC.
Parameters
......@@ -226,7 +233,6 @@ def quick_run_v2(
"""
print(f'Quick Run: {target_str}')
log.debug(
f"""input parameters:
target_str: {target_str}
......@@ -257,7 +263,7 @@ def quick_run_v2(
camera.switch['bias_vp'] = False
camera.switch['bias_hp'] = False
camera.switch['bias_ci'] = False
camera.switch['bias_shift'] = False
camera.switch['bias_shift'] = False
camera.switch['badcolumn'] = False
camera.switch['flat'] = False
camera.switch['cte'] = False
......@@ -283,10 +289,9 @@ def quick_run_v2(
)
def deduplicate_names_add_count(names: list):
"""remove duplicate names and add count"""
for i in range(len(names)-1,-1,-1):
for i in range(len(names)-1, -1, -1):
if names.count(names[i]) > 1:
names[i] = names[i] + '_' + str(names.count(names[i]))
......@@ -302,7 +307,7 @@ def observation_simulation_from_config(obs_file, config_file):
config file.
Examples:
see examples in `example` folder.
see examples in `example` folder.
"""
config = {}
if config_file:
......@@ -315,17 +320,17 @@ def observation_simulation_from_config(obs_file, config_file):
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]
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 = []
......@@ -333,11 +338,10 @@ def observation_simulation_from_config(obs_file, config_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):
try:
with open(file, 'r') as fid:
......@@ -405,15 +409,16 @@ def observation_simulation_from_config(obs_file, config_file):
except Exception as e:
log.error(f"{info_text} failed with {type(e).__name__}{e}.\n\n {traceback.format_exc()}")
def main(argv=None):
"""
Command line interface of csst_cpic_sim
Parameters
-----------
argv: list
input arguments. Default is None. If None, sys.argv is used.
"""
parser = argparse.ArgumentParser(description='Cpic obsevation image simulation')
parser.set_defaults(func=lambda _: parser.print_usage())
......@@ -426,10 +431,18 @@ def main(argv=None):
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(
'-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):
......@@ -464,7 +477,7 @@ def main(argv=None):
args = parser.parse_args(argv)
args.func(args)
# if __name__ == '__main__': # pragma: no cover
# sys.exit(main())
......
......@@ -79,9 +79,9 @@ def ideal_focus_image(
np.ndarray
The ideal focus image.
"""
focal_image = np.zeros(platesize)
focal_shape = np.array(platesize)[::-1] # x, y
focal_shape = np.array(platesize)[::-1] # x, y
if not targets:
return focal_image
......@@ -109,7 +109,7 @@ def ideal_focus_image(
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
......@@ -133,9 +133,9 @@ def focal_convolve(
rotation: float = 0,
nsample: int = 5,
error: float = 0,
platesize: list = [1024, 1024]) -> np.ndarray :
platesize: list = [1024, 1024]) -> np.ndarray:
"""PSF convolution of the ideal focus image.
Parameters
----------
band: str
......@@ -163,10 +163,10 @@ def focal_convolve(
area = config['aperature_area']
filter = filter_throughput(band)
throughput = filter.throughput
wave = filter.wave
throughput_criterion = throughput.max() * 0.1
wave_criterion = wave[throughput > throughput_criterion]
min_wave = wave_criterion[0]
......@@ -179,33 +179,33 @@ def focal_convolve(
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]
all_fp_image = []
if not targets:
return np.zeros((platesize[1], platesize[0]))
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[1:], platescale, platesize, init_shifts, rotation)
psf = single_band_psf(center_wavelength, error=error)
_, _, cstar_sp, _ = targets[0]
cstar_flux = (cstar_sp * i_band).integrate()
cstar_psf = single_band_masked_psf(center_wavelength, error=error, shift=init_shifts)
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
all_fp_image.append(c_fp_image * area) # trans to photon/second
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[1:], platescale, platesize, init_shifts, rotation)
psf = single_band_psf(center_wavelength, error=error)
_, _, cstar_sp, _ = targets[0]
cstar_flux = (cstar_sp * i_band).integrate()
cstar_psf = single_band_masked_psf(center_wavelength, error=error, shift=init_shifts)
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
all_fp_image.append(c_fp_image * area) # trans to photon/second
return np.array(all_fp_image).sum(axis=0)
......
import os, pickle
import os
import pickle
import numpy as np
from astropy.io import fits
......@@ -23,12 +24,12 @@ pupil_diameter = 2 # m
F_number = 83
focal_length = pupil_diameter * F_number
pupil_shape = apm.shape
pupil_rate = apm_header['PHRATE'] # meter/pixel
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 * 2 # just gauss a number
second_pupil_size = pupil_diameter * 2 # just gauss a number
second_aperture = make_circular_aperture(second_pupil_size)(pupil_grid)
lyot_stop = ComplexSurfaceApodizer(second_aperture, second_aperture*0, lambda _: 2)
......@@ -73,12 +74,13 @@ aberration = SurfaceApodizer(
aberration_distance = 80 * focal_length
aberration = SurfaceAberrationAtDistance(aberration, aberration_distance)
def single_band_masked_psf(
wavelength: float,
error: float = 0,
wavelength: float,
error: float = 0,
shift: list = [0, 0]) -> np.ndarray:
"""CPIC PSF considering the focal plane mask.
Parameters
-----------
wavelength : float
......@@ -91,7 +93,7 @@ def single_band_masked_psf(
Returns
----------
psf : np.ndarray
psf in the focal plane. Normalized as the input flux is 1.
psf in the focal plane. Normalized as the input flux is 1.
(Note that total flux of the psf is not 1, because it is masked)
"""
error = np.random.normal(0, error*1e-9, actuator.shape)
......@@ -99,8 +101,10 @@ def single_band_masked_psf(
wf = Wavefront(aperture, wavelength)
shift = np.array(shift) * ARCSEC2RAD / 2
tiptilt_mirror.actuators = shift
wf = tiptilt_mirror(wf)
wf = aberration(wf)
wf = tiptilt_mirror(wf)
first_focal = prop_full_frame(deformable_mirror(wf))
strength = first_focal.intensity.shaped.sum()
......@@ -109,8 +113,43 @@ def single_band_masked_psf(
psf = second_focal.intensity.shaped
return psf / strength
def single_band_shift_psf(
wavelength: float,
error: float = 0,
shift: list = [0, 0]) -> np.ndarray:
"""CPIC PSF considering the focal plane mask.
Parameters
-----------
wavelength : float
observation wavelength in meter
error : float
deformable mirror control error in nm
shift : list
angular shift of the target in arcsec.
Returns
----------
psf : np.ndarray
psf in the focal plane. Normalized as the input flux is 1.
(Note that total flux of the psf is not 1, because it is masked)
"""
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 = aberration(wf)
wf = tiptilt_mirror(wf)
first_focal = prop_full_frame(deformable_mirror(wf))
image = np.array(first_focal.intensity.shaped)
return image / image.sum()
def single_band_psf(
wavelength: float,
wavelength: float,
error: float = 0) -> np.ndarray:
"""CPIC PSF without considering the focal plane mask.
Used to simulate the off-axis PSF.
......@@ -126,7 +165,6 @@ def single_band_psf(
----------
psf : np.ndarray
psf in the focal plane. Normalized. The total flux is 1.
"""
error = np.random.normal(0, error*1e-9, actuator.shape)
deformable_mirror.actuators = actuator + error
......@@ -136,7 +174,6 @@ def single_band_psf(
image = np.array(first_focal.intensity.shaped)
return image / image.sum()
# def single_band_psf(wavelength, error=0, aber_phase=None):
# error = np.random.normal(0, error*1e-9, actuator.shape)
# deformable_mirror.actuators = actuator + error
......@@ -148,4 +185,4 @@ def single_band_psf(
# wf = other_error(wf)
# 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
# return image / image.sum()
import os, re, json, yaml
import os
import re
import json
import yaml
from typing import Union
import numpy as np
from scipy import constants
......@@ -24,9 +27,10 @@ BCC_MODEL_FOLDER = config['bcc_model']
MAG_SYSTEM = config['mag_system']
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
......@@ -45,16 +49,17 @@ class AlbedoCat(S.spectrum.TabularSpectralElement, S.catalog.Icat):
---------
Color Classification of Extrasolar Giant Planets: Prospects and Cautions
Natasha E. Batalha et al 2018 AJ 156 158
"""
def __init__(self,
phase: float,
metallicity: float,
f_sed: float,
):
def __init__(
self,
phase: float,
metallicity: float,
f_sed: float,
):
catname = 'BCCalbedo'
self.isAnalytic=False
self.isAnalytic = False
self.name = f"{catname}_{phase}_{metallicity}_{f_sed}"
self.parameter_names = ['phase', 'metallicity', 'f_sed']
self.catalog_folder = BCC_MODEL_FOLDER
......@@ -66,19 +71,19 @@ class AlbedoCat(S.spectrum.TabularSpectralElement, S.catalog.Icat):
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)
list0, list1 = self._breakList(indices, 0, phase)
list2,list3 = self._breakList(list0, 1, metallicity)
list4,list5 = self._breakList(list1, 1, metallicity)
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)
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)
......@@ -113,7 +118,7 @@ class AlbedoCat(S.spectrum.TabularSpectralElement, S.catalog.Icat):
column = name.split('[')[1][:-1]
filename = f"{self.catalog_folder}/{filename}"
# sp = S.spectrum.TabularSpectralElement(filename, thrucol=column)
with fits.open(filename) as td:
......@@ -255,21 +260,22 @@ def star_photlam(
star_sp.convert('photlam')
return star_sp
def standard_spectrum(
magnitude: float) -> S.ArraySpectrum:
"""Standard spectrum of magnitude system.
For example, the standard_spectrum of vega megnitude is vega spectrum.
"""Standard spectrum of magnitude system.
For example, the standard_spectrum of vega megnitude is vega spectrum.
The standard spectrum of ab magnitude is a flat spectrum.
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(
......@@ -282,7 +288,8 @@ def standard_spectrum(
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) -> S.spectrum.ArraySpectralElement:
......@@ -407,7 +414,7 @@ def detect_template_path(template: str) -> str:
----------
template: str
template file name
Returns
-------
str
......@@ -417,14 +424,14 @@ def detect_template_path(template: str) -> str:
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
......@@ -459,7 +466,7 @@ class TargetOjbect(object):
center star object is used to calculate the projection x, y of each target according to its ra and dec
also center star's distance is used to calculate seperation of planet
Examples:
Examples:
--------
cstar = {
......@@ -470,7 +477,7 @@ class TargetOjbect(object):
'sptype': 'G0III'
'sp_model': 'star'
}
stars = {
'magnitude': 15,
'ra': '120.001d',
......@@ -478,27 +485,27 @@ class TargetOjbect(object):
'sptype': 'F0III',
'sp_model': 'blackbody'
}
planets = {
'radius': 2,
'pangle': 60,
'coe_cloud': 0.3,
'coe_metal': 0.7,
'separation': 0.5,
'phase_angle': 90,
'phase_angle': 90,
'sp_model': 'hybrid_planet',
'image': 'extend_planet.fits'
}
# planet using input costum albedo spectrum!
# Note that the albedo spectrum is not normalized!
# 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.
# See the documentation of pysynphot for details.
planets = {
'pangle': 60,
'separation': 0.5,
'separation': 0.5,
'sp_spectrum': 'template_planet',
'template': 'planet_albedo.fits'
}
......@@ -556,26 +563,27 @@ class TargetOjbect(object):
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'
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,
phase_angle: float,
radius: float) -> float:
"""
calculate the contrast of a planet
Calculate the contrast of a planet.
Parameters
----------
......@@ -610,16 +618,17 @@ def planet_contrast(
contrast = (radius / sep_3d * R_JUPITER_METER / AU_METER)**2 * phase
return contrast
def spectrum_generator(
targets: dict) -> list:
"""
generate the spectrum due to the input target list
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:
The values are:
- cstar: dict
- center star information. must contain keywords ra, dec, distance, magnitude, sptype
- objects: list of dict, optional, recommended! V2.0 new!
......@@ -654,17 +663,17 @@ def target_file_load(
----------
target: Union[dict, str]
target file name, formated string or target dict.
Outputs
--------
target: dict
dictionary of target. Format of input
dictionary of target. Format of input
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.
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.
......@@ -679,20 +688,20 @@ def target_file_load(
if isinstance(target, dict):
return target
if not target: # None or empty string
if not target: # None or empty string
return {}
if isinstance(target, str): #filename or formatted string
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_name = os.path.basename(target_file)[:-5]
file_search = [target_file, os.path.join(catalog_folder, target_file)]
for file in file_search:
if os.path.isfile(file):
with open(file) as fid:
......
......@@ -2,7 +2,7 @@ import numpy as np
import scipy.ndimage as nd
import logging
import random
import matplotlib.pyplot as plt
import matplotlib.pyplot as plt
# DO NOT IMPORT CPICIMGSIM MODULES HERE
......@@ -66,7 +66,7 @@ def region_replace(
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.
[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.
......@@ -133,7 +133,7 @@ def psf_imshow(psf, vmin=1e-8, vmax=0.1, log=True, region=1):
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
......
band: f661
band: f662
emgain: 1
emset: -1
expt: 10
......
band: f661
band: f662
emgain: 200
emset: -1
expt: 100
......
band: f565
band: f520
emgain: 1
emset: -1
expt: 10
......
band: f565
band: f520
emgain: 50
emset: -1
expt: 50
......
band: f661
band: f662
emgain: 50
emset: -1
expt: 20
......
band: f565
band: f520
emgain: 1
emset: -1
expt: 10
......
band: f565
band: f520
emgain: 2
emset: -1
expt: 10
......
band: f565
band: f720
emgain: 2
emset: -1
expt: 10
......
band: f565
band: f520
emgain: 50
emset: -1
expt: 10
......
......@@ -12,6 +12,10 @@ bands:
f743: ${cpism_refdata}/throughtput/f743_total.fits
f883: ${cpism_refdata}/throughtput/f883_total.fits
f565: ${cpism_refdata}/throughtput/f565_total.fits
f520: ${cpism_refdata}/throughtput/f520.fits
f662: ${cpism_refdata}/throughtput/f662.fits
f850: ${cpism_refdata}/throughtput/f850.fits
f720: ${cpism_refdata}/throughtput/f720.fits
diameter: 2
platescale: 0.016153
datamodel: ${cpism_refdata}/io/csst-cpic-l0.yaml
......
import sys
from csst_cpic_sim.main import main
sys.exit(main())
\ No newline at end of file
sys.exit(main())
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