From 1aa0db4820aa9a3ee3a184b2d79f5c557e0ac319 Mon Sep 17 00:00:00 2001 From: yan Date: Wed, 8 Jan 2025 15:49:33 +0800 Subject: [PATCH] update --- csst_ifs_sim/csst_ifs_sim_bug_250104.py | 6240 ----------------------- 1 file changed, 6240 deletions(-) delete mode 100644 csst_ifs_sim/csst_ifs_sim_bug_250104.py diff --git a/csst_ifs_sim/csst_ifs_sim_bug_250104.py b/csst_ifs_sim/csst_ifs_sim_bug_250104.py deleted file mode 100644 index 8803553..0000000 --- a/csst_ifs_sim/csst_ifs_sim_bug_250104.py +++ /dev/null @@ -1,6240 +0,0 @@ -#!/usr/bin/env python3 -# -*- coding: utf-8 -*- -""" -Created on Thu Apr 11 15:18:57 2024 - -@author: yan -""" -from astropy.utils.iers import conf -import galsim -from scipy.integrate import simps -import pandas as pd -from datetime import datetime, timedelta -from astropy.coordinates import get_sun -from astropy.time import Time -from astropy import units as u -from astropy.coordinates import SkyCoord, EarthLocation -from scipy import interpolate -import numpy as np -import scipy.io as sio -from astropy.io import fits -from scipy.signal import fftconvolve -from scipy import ndimage -import cmath -import configparser as ConfigParser -import os -from scipy.interpolate import interp1d -import astropy.coordinates as coord -import ctypes -import sys - -##sys.path.append('./csst_ifs_sim') -conf.auto_max_age = None - -""" -""" -""" -Contact Information: zhaojunyan@shao.ac.cn - -The CSST IFS Image Simulator -============================================= - -This file contains an image simulator for the CSST IFS. - -The approximate sequence of events in the simulator is as follows: - -#. Read in a configuration file, which defines for example, - detector characteristics (bias, dark and readout noise, gain, - pixel scale,exposure time etc.). -#. Read in another file containing charge trap definitions (for CTI modelling). -#. Read in a file defining the cosmic rays (trail lengths and - cumulative distributions). -#. Load the wavefront aberration data used to calculate PSF . -#. Apply a multiplicative flat-field map to emulate pixel-to-pixel -non-uniformity. -#. Add cosmic ray tracks onto the CCD with random positions but known -distribution. -#. Apply detector charge bleeding in column direction. -#. Add constant dark current and background light from Zodiacal -light. -#. Add photon (Poisson) noise -#. Add cosmetic defects from an input file. -#. Add pre- and overscan regions in the serial direction. -#. Apply the CDM03 radiation damage model. -#. Apply non-linearity model to the pixel data. -#. Add readout noise selected from a Gaussian distribution. -#. Convert from electrons to ADUs using a given gain factor. -#. Add a given bias level and discretise the counts -#. Finally the simulated image is converted to a FITS file. - -Warning:: The code is still work in progress and new features are being added. -The code has been tested, but nevertheless bugs may be lurking in corners, so -please report any weird or inconsistent simulations to the author. - -Note:: This class is Python 3 compatible. -2023.06.13 , add new lamp hole simulation -2023.12.26 add date frame transfer effect and software version -2024.1.19 add multiple exposure mode function; -2024.3.20 --- update the fits header as defined order 0 data format -2024.5.10 ---updata and correct the bug of frame transfer effect simulation - -""" -######################## functions definition ################################ - - -class CDM03bidir(): - """ - Charge Transfer Inefficiency - ============================ - - This file contains a simple class to run a CDM03 CTI model developed by Alex Short (ESA). - - This now contains both the official CDM03 and a new version that allows different trap - parameters in parallel and serial direction. - - :requires: NumPy - :requires: CDM03 (FORTRAN code, f2py -c -m cdm03bidir cdm03bidir.f90) - - Class to run CDM03 CTI model, class Fortran routine to perform the actual CDM03 calculations. - - :param settings: input parameters - :type settings: dict - :param data: input data to be radiated - :type data: ndarray - :param log: instance to Python logging - :type log: logging instance - """ - def __init__(self, settings, data, log=None): - """ - Class constructor. - - :param settings: input parameters - :type settings: dict - :param data: input data to be radiated - :type data: ndarray - :param log: instance to Python logging - :type log: logging instance - """ - self.data = data - self.values = dict(quads=(0,1,2,3), xsize=2048, ysize=2066, dob=0.0, rdose=8.0e9) - self.values.update(settings) - self.log = log - self._setupLogger() - - # #default CDM03 settings - # self.params = dict(beta_p=0.6, beta_s=0.6, fwc=200000., vth=1.168e7, vg=6.e-11, t=20.48e-3, - # sfwc=730000., svg=1.0e-10, st=5.0e-6, parallel=1., serial=1.) - - self.params = dict(beta_p=0.6, beta_s=0.6, fwc=100000., vth=1.168e7, vg=6.e-11, t=1.0e-3, - sfwc=700000., svg=1.0e-10, st=1.0e-6, parallel=1., serial=0.) - - #update with inputs - self.params.update(self.values) - - #read in trap information - trapdata = np.loadtxt(self.values['dir_path']+self.values['paralleltrapfile']) - if trapdata.ndim > 1: - self.nt_p = trapdata[:, 0] - self.sigma_p = trapdata[:, 1] - self.taur_p = trapdata[:, 2] - else: - #only one trap species - self.nt_p = [trapdata[0],] - self.sigma_p = [trapdata[1],] - self.taur_p = [trapdata[2],] - - trapdata = np.loadtxt(self.values['dir_path']+self.values['serialtrapfile']) - if trapdata.ndim > 1: - self.nt_s = trapdata[:, 0] - self.sigma_s = trapdata[:, 1] - self.taur_s = trapdata[:, 2] - else: - #only one trap species - self.nt_s = [trapdata[0],] - self.sigma_s = [trapdata[1],] - self.taur_s = [trapdata[2],] - - #scale thibaut's values - if 'thibaut' in self.values['parallelTrapfile']: - self.nt_p /= 0.576 #thibaut's values traps / pixel - self.sigma_p *= 1.e4 #thibaut's values in m**2 - if 'thibaut' in self.values['serialTrapfile']: - self.nt_s *= 0.576 #thibaut's values traps / pixel #should be division? - self.sigma_s *= 1.e4 #thibaut's values in m**2 - - - def _setupLogger(self): - """ - Set up the logger. - """ - self.logger = True - - - def applyRadiationDamage(self, data, iquadrant=0): - """ - Apply radian damage based on FORTRAN CDM03 model. The method assumes that - input data covers only a single quadrant defined by the iquadrant integer. - - :param data: imaging data to which the CDM03 model will be applied to. - :type data: ndarray - - :param iquandrant: number of the quadrant to process - :type iquandrant: int - - cdm03 - Function signature:: - - sout = cdm03(sinp,iflip,jflip,dob,rdose,in_nt,in_sigma,in_tr,[xdim,ydim,zdim]) - Required arguments: - sinp : input rank-2 array('d') with bounds (xdim,ydim) - iflip : input int - jflip : input int - dob : input float - rdose : input float - in_nt : input rank-1 array('d') with bounds (zdim) - in_sigma : input rank-1 array('d') with bounds (zdim) - in_tr : input rank-1 array('d') with bounds (zdim) - Optional arguments: - xdim := shape(sinp,0) input int - ydim := shape(sinp,1) input int - zdim := len(in_nt) input int - Return objects: - sout : rank-2 array('d') with bounds (xdim,ydim) - - .. Note:: Because Python/NumPy arrays are different row/column based, one needs - to be extra careful here. NumPy.asfortranarray will be called to get - an array laid out in Fortran order in memory. Before returning the - array will be laid out in memory in C-style (row-major order). - - :return: image that has been run through the CDM03 model - :rtype: ndarray - """"" - #return data - - iflip = iquadrant / 2 - jflip = iquadrant % 2 - - params = [self.params['beta_p'], self.params['beta_s'], self.params['fwc'], self.params['vth'], - self.params['vg'], self.params['t'], self.params['sfwc'], self.params['svg'], self.params['st'], - self.params['parallel'], self.params['serial']] - - if self.logger: - self.log.info('nt_p=' + str(self.nt_p)) - self.log.info('nt_s=' + str(self.nt_s)) - self.log.info('sigma_p= ' + str(self.sigma_p)) - self.log.info('sigma_s= ' + str(self.sigma_s)) - self.log.info('taur_p= ' + str(self.taur_p)) - self.log.info('taur_s= ' + str(self.taur_s)) - self.log.info('dob=%f' % self.values['dob']) - self.log.info('rdose=%e' % self.values['rdose']) - self.log.info('xsize=%i' % data.shape[1]) - self.log.info('ysize=%i' % data.shape[0]) - self.log.info('quadrant=%i' % iquadrant) - self.log.info('iflip=%i' % iflip) - self.log.info('jflip=%i' % jflip) - ################ - ############################################ - -################################################################################# -###modify - #sys.path.append('../so') - #from .ifs_so import cdm03bidir - try: - from ifs_so import cdm03bidir - except: - from .ifs_so import cdm03bidir - - # from ifs_so.cdm03.cpython-38-x86_64-linux-gnu import cdm03bidir - # import cdm03bidir - CTIed = cdm03bidir.cdm03(np.asfortranarray(data), - jflip, iflip, - self.values['dob'], self.values['rdose'], - self.nt_p, self.sigma_p, self.taur_p, - self.nt_s, self.sigma_s, self.taur_s, - params, - [data.shape[0], data.shape[1], len(self.nt_p), len(self.nt_s), len(self.params)]) - - - return np.asanyarray(CTIed) - -################################################################################################################# - -################################################################################################################# - -""" -These functions can be used for logging information. - -.. Warning:: logger is not multiprocessing safe. - -:version: 0.3 -""" -import logging -import logging.handlers - - -def lg(log_filename, loggername='logger'): - """ - Sets up a logger. - - :param: log_filename: name of the file to save the log. - :param: loggername: name of the logger - - :return: logger instance - """ - # create logger - logger = logging.getLogger(loggername) - logger.setLevel(logging.DEBUG) - # Add the log message handler to the logger - handler = logging.handlers.RotatingFileHandler(log_filename) - #maxBytes=20, backupCount=5) - # create formatter - formatter = logging.Formatter('%(asctime)s - %(module)s - %(funcName)s - %(levelname)s - %(message)s') - # add formatter to ch - handler.setFormatter(formatter) - # add handler to logger - - if (logger.hasHandlers()): - logger.handlers.clear() - - logger.addHandler(handler) - - return logger - - -############################################################################## - -def IFSinformation(): - """ - IFS Instrument Model - ==================== - - The file provides a function that returns IFS related information such as pixel - size, dark current, gain... - - Returns a dictionary describing VIS. The following information is provided - (id: value - reference):: - - :return: instrument model parameters - :rtype: dict - """ - - ######################################################################################################### - #out = dict(readnoise=4, pixel_size=0.1, dark=0.0008333, fullwellcapacity=90000, bluesize=4000, redsize=6000, readtime=300.) - out=dict() - out.update({'dob' : 0, 'rdose' : 8.0e9, - 'parallelTrapfile' : 'cdm_euclid_parallel.dat', 'serialTrapfile' : 'cdm_euclid_serial.dat', - 'beta_s' : 0.6, 'beta_p': 0.6, 'fwc' : 90000, 'vth' : 1.168e7, 't' : 20.48e-3, 'vg' : 6.e-11, - 'st' : 5.0e-6, 'sfwc' : 730000., 'svg' : 1.0e-10}) - return out - - -def CCDnonLinearityModel(data, beta=6e-7): - """ - - The non-linearity is modelled based on the results presented. - :param data: data to which the non-linearity model is being applied to - :type data: ndarray - - :return: input data after conversion with the non-linearity model - :rtype: float or ndarray - """ - out = data-beta*data**2 - - return out - -############################################################################# - -from scipy.interpolate import InterpolatedUnivariateSpline - -class cosmicrays(): - """ - - Cosmic Rays - =========== - - This simple class can be used to include cosmic ray events to an image. - By default the cosmic ray events are drawn from distributions describing - the length and energy of the events. Such distributions can be generated - for example using Stardust code (http://ieeexplore.ieee.org/stamp/stamp.jsp?arnumber=04636917). - The energy of the cosmic ray events can also be set to constant for - testing purposes. The class can be used to draw a single cosmic ray - event or up to a covering fraction. - - - Cosmic ray generation class. Can either draw events from distributions or - set the energy of the events to a constant. - - :param log: logger instance - :param image: image to which cosmic rays are added to (a copy is made not to change the original numpy array) - :param crInfo: column information (cosmic ray file) - :param information: cosmic ray track information (file containing track length and energy information) and - exposure time. - """ - def __init__(self, log, image, crInfo=None, information=None): - """ - Cosmic ray generation class. Can either draw events from distributions or - set the energy of the events to a constant. - - :param log: logger instance - :param image: image to which cosmic rays are added to (a copy is made not to change the original numpy array) - :param crInfo: column information (cosmic ray file) - :param information: cosmic ray track information (file containing track length and energy information) and - exposure time. - """ - #setup logger - self.log = log - - #image and size - self.image = image.copy() - self.ysize, self.xsize = self.image.shape - - #set up the information dictionary, first with defaults and then overwrite with inputs if given - self.information = (dict(cosmicraylengths='/home/yan/csst-master/data/cdf_cr_length.dat', - cosmicraydistance='/home/yan/csst-master/data/cdf_cr_total.dat', - exptime=565)) - if information is not None: - self.information.update(information) - - if crInfo is not None: - self.cr = crInfo - else: - self._readCosmicrayInformation() -############################################################################## - - def _cosmicRayIntercepts(self, lum, x0, y0, l, phi): - """ - Derive cosmic ray streak intercept points. - - :param lum: luminosities of the cosmic ray tracks - :param x0: central positions of the cosmic ray tracks in x-direction - :param y0: central positions of the cosmic ray tracks in y-direction - :param l: lengths of the cosmic ray tracks - :param phi: orientation angles of the cosmic ray tracks - - :return: cosmic ray map (image) - :rtype: nd-array - """ - #create empty array - crImage = np.zeros((self.ysize, self.xsize), dtype=np.float64) - - #x and y shifts - dx = l * np.cos(phi) / 2. - dy = l * np.sin(phi) / 2. - mskdx = np.abs(dx) < 1e-8 - mskdy = np.abs(dy) < 1e-8 - dx[mskdx] = 0. - dy[mskdy] = 0. - - #pixels in x-direction - ilo = np.round(x0.copy() - dx) - msk = ilo < 0. - ilo[msk] = 0 - ilo = ilo.astype(int) - - ihi = 1 + np.round(x0.copy() + dx) - msk = ihi > self.xsize - ihi[msk] = self.xsize - ihi = ihi.astype(int) - - #pixels in y-directions - jlo = np.round(y0.copy() - dy) - msk = jlo < 0. - jlo[msk] = 0 - jlo = jlo.astype(int) - - jhi = 1 + np.round(y0.copy() + dy) - msk = jhi > self.ysize - jhi[msk] = self.ysize - jhi = jhi.astype(int) - - #loop over the individual events - for i, luminosity in enumerate(lum): - n = 0 # count the intercepts - - u = [] - x = [] - y = [] - - #Compute X intercepts on the pixel grid - if ilo[i] < ihi[i]: - for xcoord in range(ilo[i], ihi[i]): - ok = (xcoord - x0[i]) / dx[i] - if np.abs(ok) <= 0.5: - n += 1 - u.append(ok) - x.append(xcoord) - y.append(y0[i] + ok * dy[i]) - else: - for xcoord in range(ihi[i], ilo[i]): - ok = (xcoord - x0[i]) / dx[i] - if np.abs(ok) <= 0.5: - n += 1 - u.append(ok) - x.append(xcoord) - y.append(y0[i] + ok * dy[i]) - - #Compute Y intercepts on the pixel grid - if jlo[i] < jhi[i]: - for ycoord in range(jlo[i], jhi[i]): - ok = (ycoord - y0[i]) / dy[i] - if np.abs(ok) <= 0.5: - n += 1 - u.append(ok) - x.append(x0[i] + ok * dx[i]) - y.append(ycoord) - else: - for ycoord in range(jhi[i], jlo[i]): - ok = (ycoord - y0[i]) / dy[i] - if np.abs(ok) <= 0.5: - n += 1 - u.append(ok) - x.append(x0[i] + ok * dx[i]) - y.append(ycoord) - - #check if no intercepts were found - if n < 1: - xc = int(np.floor(x0[i])) - yc = int(np.floor(y0[i])) - crImage[yc, xc] += luminosity - - #Find the arguments that sort the intersections along the track - u = np.asarray(u) - x = np.asarray(x) - y = np.asarray(y) - - args = np.argsort(u) - - u = u[args] - x = x[args] - y = y[args] - - #Decide which cell each interval traverses, and the path length - for i in range(1, n - 1): - w = u[i + 1] - u[i] - cx = int(1 + np.floor((x[i + 1] + x[i]) / 2.)) - cy = int(1 + np.floor((y[i + 1] + y[i]) / 2.)) - - if 0 <= cx < self.xsize and 0 <= cy < self.ysize: - crImage[cy, cx] += (w * luminosity) - - return crImage - - - - def _drawEventsToCoveringFactor(self, coveringFraction=3.0, limit=1000, verbose=False): - """ - Generate cosmic ray events up to a covering fraction and include it to a cosmic ray map (self.cosmicrayMap). - - :param coveringFraction: covering fraction of cosmic rya events in per cent of total number of pixels - :type coveringFraction: float - :param limit: limiting energy for the cosmic ray event [None = draw from distribution] - :type limit: None or float - :param verbose: print out information to stdout - :type verbose: bool - - - :return: None - """ - - if coveringFraction>1 or coveringFraction<0.001: - self.log.error( - 'coveringFraction error, it shoub be in [0.1, 1]!') - raise ValueError( - 'coveringFraction error, it shoub be in [0.1, 1]!') - - - self.cosmicrayMap = np.zeros((self.ysize, self.xsize)) - - #how many events to draw at once, too large number leads to exceeding the covering fraction - cr_n = int(500 * self.information['exptime'] / 565. * coveringFraction / 1.4) - - if cr_n<1: - cr_n=1 - - - covering = 0.0 - - while covering < coveringFraction: - #pseudo-random numbers taken from a uniform distribution between 0 and 1 - np.random.seed() - luck = np.random.rand(cr_n) - - #draw the length of the tracks - ius = InterpolatedUnivariateSpline(self.cr['cr_cdf'], self.cr['cr_u']) - self.cr['cr_l'] = ius(luck) - - if limit is None: - ius = InterpolatedUnivariateSpline(self.cr['cr_cde'], self.cr['cr_v']) - self.cr['cr_e'] = ius(luck) - else: - #set the energy directly to the limit - self.cr['cr_e'] = np.asarray([limit,]) - - #Choose the properties such as positions and an angle from a random Uniform dist - np.random.seed() - cr_x = self.xsize * np.random.rand(int(np.floor(cr_n))) - - np.random.seed() - cr_y = self.ysize * np.random.rand(int(np.floor(cr_n))) - - np.random.seed() - cr_phi = np.pi * np.random.rand(int(np.floor(cr_n))) - - #find the intercepts - self.cosmicrayMap += self._cosmicRayIntercepts(self.cr['cr_e'], cr_x, cr_y, self.cr['cr_l'], cr_phi) - - #count the covering factor - area_cr = np.count_nonzero(self.cosmicrayMap) - covering = 100.*area_cr / (self.xsize*self.ysize) - - - def addUpToFraction(self, coveringFraction, limit=None, verbose=False): - """ - Add cosmic ray events up to the covering Fraction. - - :param coveringFraction: covering fraction of cosmic rya events in per cent of total number of pixels - :type coveringFraction: float - :param limit: limiting energy for the cosmic ray event [None = draw from distribution] - :type limit: None or float - :param verbose: print out information to stdout - :type verbose: bool - - :return: image with cosmic rays - :rtype: ndarray - """ - self._drawEventsToCoveringFactor(coveringFraction, limit=limit, verbose=verbose) - - #paste cosmic rays - self.image += self.cosmicrayMap - - return self.image - - -############################################################################### -def transRaDec2D(ra, dec): - """ - - ra,dec转为竞天程序里的ob, 赤道坐标系下的笛卡尔三维坐标xyz. - Parameters - ---------- - ra : float - ra in deg. - dec : float - dec in deg. - - Returns - ------- - TYPE - DESCRIPTION. - - """ - x1 = np.cos(dec / 57.2957795) * np.cos(ra / 57.2957795) - y1 = np.cos(dec / 57.2957795) * np.sin(ra / 57.2957795) - z1 = np.sin(dec / 57.2957795) - return np.array([x1, y1, z1]) -############################################################################### - -################################################################ - - -def ill2flux(path, E): - """ - - - Parameters - ---------- - path : TYPE - DESCRIPTION. - E : TYPE - DESCRIPTION. - - Returns - ------- - wave0 : TYPE - DESCRIPTION. - spec_scaled : TYPE - DESCRIPTION. - - """ - - # use template from sky_bkg (background_spec_hst.dat) - filename = path+'IFS_inputdata/refs/background_spec_hst.dat' - cat_spec = pd.read_csv(filename, sep='\s+', header=None, comment='#') - wave0 = cat_spec[0].values # A - spec0 = cat_spec[2].values # erg/s/cm^2/A/arcsec^2 - - # convert erg/s/cm^2/A/arcsec^2 to erg/s/cm^2/A/sr - flux1 = spec0 / (1/4.25452e10) - # convert erg/s/cm^2/A/sr to W/m^2/sr/um - flux2 = flux1 * 10 - - # 对接收面积积分，输出单位 W/m^2/nm - D = 2 # meter - f = 28 # meter, 焦距，转换关系来源于王维notes. - flux3 = flux2 * np.pi * D**2 / 4 / f**2 / 10**3 - - f = interp1d(wave0, flux3) - wave_range = np.arange(3800, 7800) - flux3_mean = f(wave_range) - delta_lamba = 0.1 # nm - E0 = np.sum(flux3_mean * delta_lamba) - - factor = E / E0 - spec_scaled = factor * spec0 - - return wave0, spec_scaled - -############################################################## - -############################################################################## - - -def earth_angle(time_jd, x_sat, y_sat, z_sat, ra_obj, dec_obj): - """ - # - - Parameters - ---------- - time_jd : TYPE - DESCRIPTION. - x_sat : TYPE - DESCRIPTION. - y_sat : TYPE - DESCRIPTION. - z_sat : TYPE - DESCRIPTION. - ra_obj : TYPE - DESCRIPTION. - dec_obj : TYPE - DESCRIPTION. - - Returns - ------- - angle : TYPE - DESCRIPTION. - - """ - - ra_sat = np.arctan2(y_sat, x_sat) / np.pi * 180 - dec_sat = np.arctan2(z_sat, np.sqrt(x_sat**2+y_sat**2)) / np.pi * 180 - radec_sat = SkyCoord(ra=ra_sat*u.degree, - dec=dec_sat*u.degree, frame='gcrs') - lb_sat = radec_sat.transform_to('geocentrictrueecliptic') - - # get the obj location - radec_obj = SkyCoord(ra=ra_obj*u.degree, - dec=dec_obj*u.degree, frame='gcrs') - lb_obj = radec_obj.transform_to('geocentrictrueecliptic') - - # calculate the angle between sub-satellite point and the earth side - earth_radius = 6371 # km - sat_height = np.sqrt(x_sat**2 + y_sat**2 + z_sat**2) - angle_a = np.arcsin(earth_radius/sat_height) / np.pi * 180 - - # calculate the angle between satellite position and the target position - angle_b = lb_sat.separation(lb_obj) - - # calculat the earth angle - angle = 180 - angle_a - angle_b.degree - - return angle - -##################################################################### - - -class StrayLight(object): - - def __init__(self, path, jtime=2460843., sat=np.array([0, 0, 0]), radec=np.array([0, 0])): - self.path = path - self.jtime = jtime - self.sat = sat - self.equator = coord.SkyCoord( - radec[0]*u.degree, radec[1]*u.degree, frame='icrs') - self.ecliptic = self.equator.transform_to('barycentrictrueecliptic') - self.pointing = transRaDec2D(radec[0], radec[1]) - self.slcdll = ctypes.CDLL( - self.path+'IFS_inputdata/refs/libstraylight.so') # dylib - - self.slcdll.Calculate.argtypes = [ctypes.c_double, ctypes.POINTER(ctypes.c_double), - ctypes.POINTER(ctypes.c_double), ctypes.POINTER( - ctypes.c_double), - ctypes.POINTER(ctypes.c_double), ctypes.c_char_p] - - self.slcdll.PointSource.argtypes = [ctypes.c_double, ctypes.POINTER(ctypes.c_double), - ctypes.POINTER(ctypes.c_double), ctypes.POINTER( - ctypes.c_double), - ctypes.POINTER(ctypes.c_double), ctypes.c_char_p] - - self.slcdll.EarthShine.argtypes = [ctypes.c_double, ctypes.POINTER(ctypes.c_double), - ctypes.POINTER(ctypes.c_double), ctypes.POINTER( - ctypes.c_double), - ctypes.POINTER(ctypes.c_double)] - - self.slcdll.Zodiacal.argtypes = [ctypes.c_double, ctypes.POINTER(ctypes.c_double), - ctypes.POINTER(ctypes.c_double)] - self.slcdll.ComposeY.argtypes = [ctypes.POINTER(ctypes.c_double), ctypes.POINTER(ctypes.c_double), - ctypes.POINTER(ctypes.c_double)] - self.slcdll.Init.argtypes = [ - ctypes.c_char_p, ctypes.c_char_p, ctypes.c_char_p, ctypes.c_char_p] - self.deFn = self.path+"IFS_inputdata/refs/DE405" - self.PSTFn = self.path+"IFS_inputdata/refs/PST" - self.RFn = self.path + "IFS_inputdata/refs/R" - self.ZolFn = self.path+"IFS_inputdata/refs/Zodiacal" - self.brightStarTabFn = self.path+"IFS_inputdata/refs/BrightGaia_with_csst_mag" - - self.slcdll.Init(str.encode(self.deFn), str.encode( - self.PSTFn), str.encode(self.RFn), str.encode(self.ZolFn)) -############################################################################ - - def caculateStarLightFilter(self, filter='i'): - """ - - - Parameters - ---------- - filter : TYPE, optional - DESCRIPTION. The default is 'i'. - - Returns - ------- - TYPE - DESCRIPTION. - - """ - filterIndex = {'nuv': 0, 'u': 1, 'g': 2, - 'r': 3, 'i': 4, 'z': 5, 'y': 6} - - sat = (ctypes.c_double*3)() - sat[:] = self.sat - ob = (ctypes.c_double*3)() - ob[:] = self.pointing - - py1 = (ctypes.c_double*3)() - py2 = (ctypes.c_double*3)() - self.slcdll.ComposeY(ob, py1, py2) - - star_e1 = (ctypes.c_double*7)() - self.slcdll.PointSource(self.jtime, sat, ob, py1, - star_e1, str.encode(self.brightStarTabFn)) - star_e2 = (ctypes.c_double*7)() - self.slcdll.PointSource(self.jtime, sat, ob, py2, - star_e2, str.encode(self.brightStarTabFn)) - - band_star_e1 = star_e1[:][filterIndex[filter]] - band_star_e2 = star_e2[:][filterIndex[filter]] - - return max(band_star_e1, band_star_e2) -############################################################################### - - def caculateEarthShineFilter(self, filter='i'): - """ - # - - Parameters - ---------- - filter : TYPE, optional - DESCRIPTION. The default is 'i'. - - Returns - ------- - TYPE - DESCRIPTION. - - """ - filterIndex = {'nuv': 0, 'u': 1, 'g': 2, - 'r': 3, 'i': 4, 'z': 5, 'y': 6} - - sat = (ctypes.c_double*3)() - sat[:] = self.sat - ob = (ctypes.c_double*3)() - ob[:] = self.pointing - - py1 = (ctypes.c_double*3)() - py2 = (ctypes.c_double*3)() - self.slcdll.ComposeY(ob, py1, py2) - - earth_e1 = (ctypes.c_double*7)() - self.slcdll.EarthShine(self.jtime, sat, ob, py1, earth_e1) - # e[7]代表7个波段的照度 - - earth_e2 = (ctypes.c_double*7)() - self.slcdll.EarthShine(self.jtime, sat, ob, py2, earth_e2) - - band_earth_e1 = earth_e1[:][filterIndex[filter]] - band_earth_e2 = earth_e2[:][filterIndex[filter]] - - return max(band_earth_e1, band_earth_e2) - -############################################################################### -############################################################################### - - -def time2mjd(dateT): - """ - - - Parameters - ---------- - dateT : TYPE - DESCRIPTION. - - Returns - ------- - TYPE - DESCRIPTION. - - """ - t0 = datetime(1858, 11, 17, 0, 0, 0, 0) # 简化儒略日起始日 - mjd = (dateT-t0).days - mjd_s = dateT.hour*3600.0+dateT.minute*60.0 + \ - dateT.second+dateT.microsecond/1000000.0 - return mjd+mjd_s/86400.0 - -########################################################################### - - -def time2jd(dateT): - """ - - - Parameters - ---------- - dateT : TYPE - DESCRIPTION. - - Returns - ------- - TYPE - DESCRIPTION. - - """ - t0 = datetime(1858, 11, 17, 0, 0, 0, 0) # 简化儒略日起始日 - mjd = (dateT-t0).days - mjd_s = dateT.hour*3600.0+dateT.minute*60.0 + \ - dateT.second+dateT.microsecond/1000000.0 - return mjd+mjd_s/86400.0++2400000.5 - -# mjd转datetime类 -############################################################################# - - -def mjd2time(mjd): - """ - - - Parameters - ---------- - mjd : TYPE - DESCRIPTION. - - Returns - ------- - TYPE - DESCRIPTION. - - """ - t0 = datetime(1858, 11, 17, 0, 0, 0, 0) # 简化儒略日起始日 - return t0+timedelta(days=mjd) - -############################################################################# - - -def jd2time(jd): - """ - - - Parameters - ---------- - jd : TYPE - DESCRIPTION. - - Returns - ------- - TYPE - DESCRIPTION. - - """ - mjd = jd2mjd(jd) - return mjd2time(mjd) - -# mjd和jd互转 - -############################################################################ - - -def mjd2jd(mjd): - """ - - - Parameters - ---------- - mjd : TYPE - DESCRIPTION. - - Returns - ------- - TYPE - DESCRIPTION. - - """ - return mjd+2400000.5 - -########################################################################### - - -def jd2mjd(jd): - """ - - - Parameters - ---------- - jd : TYPE - DESCRIPTION. - - Returns - ------- - TYPE - DESCRIPTION. - - """ - return jd-2400000.5 - -########################################################################## - - -def dt2hmd(dt): - """ - - - Parameters - ---------- - dt : TYPE - DESCRIPTION. - - Returns - ------- - TYPE - DESCRIPTION. - - """ - ## dt is datetime - hour = dt.hour - minute = dt.minute - second = dt.second - if hour < 10: - str_h = '0'+str(hour) - else: - str_h = str(hour) - - if minute < 10: - str_m = '0'+str(minute) - else: - str_m = str(minute) - - if second < 10: - str_d = '0'+str(second+dt.microsecond*1e-6) - else: - str_d = str(second+dt.microsecond*1e-6) - - return str_h+':'+str_m+':'+str_d - -############################################################################### -############################################################################### - - -def deg2HMS(ra0, dec0): - """ - - - Parameters - ---------- - ra0 : TYPE - DESCRIPTION. - dec0 : TYPE - DESCRIPTION. - - Returns - ------- - TYPE - DESCRIPTION. - - """ - '''convert deg to ra's HMS and dec's DMS''' - c = SkyCoord(ra=ra0*u.degree, dec=dec0*u.degree) - ss = c.to_string('hmsdms') - print(ss) - for k in range(5, len(ss)-1): - - if ss[k] == 's': - s_idx = k - - if ss[k] == 'd': - d_idx = k - - if ss[k] == 'm': - m_idx = k - - if c.ra.hms.s < 10 and c.ra.hms.s > 0.1: - temp = str(c.ra.hms.s) - ra_ss = '0'+temp[:3] - if c.ra.hms.s > 10: - temp = str(c.ra.hms.s) - ra_ss = temp[:4] - - if c.ra.hms.s < 0.1: - temp = str(c.ra.hms.s) - ra_ss = '00.0' - - dms_d = ss[s_idx+2:d_idx] - dms_m = ss[d_idx+1:m_idx] - dms_s = ss[m_idx+1:m_idx+3] - - hhmmss = ss[0:2]+ss[3:5]+ra_ss - return hhmmss+dms_d+dms_m+dms_s - - -############################################################################ - -def beta_angle(x_sat, y_sat, z_sat, vx_sat, vy_sat, vz_sat, ra_obj, dec_obj): - """ - - - Parameters - ---------- - x_sat : TYPE - DESCRIPTION. - y_sat : TYPE - DESCRIPTION. - z_sat : TYPE - DESCRIPTION. - vx_sat : TYPE - DESCRIPTION. - vy_sat : TYPE - DESCRIPTION. - vz_sat : TYPE - DESCRIPTION. - ra_obj : TYPE - DESCRIPTION. - dec_obj : TYPE - DESCRIPTION. - - Returns - ------- - angle : TYPE - DESCRIPTION. - - """ - - # get the vector for next motion - x_sat2 = x_sat + vx_sat - y_sat2 = y_sat + vy_sat - z_sat2 = z_sat + vz_sat - - vector1 = np.stack((x_sat, y_sat, z_sat), axis=0) - vector2 = np.stack((x_sat2, y_sat2, z_sat2), axis=0) - vector_normal = np.cross(vector1, vector2) - - location_normal = EarthLocation.from_geocentric(vector_normal[0], - vector_normal[1], - vector_normal[2], - unit=u.km) - radec_normal = SkyCoord(ra=location_normal.lon, - dec=location_normal.lat, - frame='gcrs') - lb_normal = radec_normal.transform_to('geocentrictrueecliptic') - - radec_sun = SkyCoord(ra=ra_obj*u.degree, - dec=dec_obj*u.degree, frame='gcrs') - lb_sun = radec_sun.transform_to('geocentrictrueecliptic') - - # get the angle between normal and solar - angle = 90 - lb_normal.separation(lb_sun).degree - - return angle - -############################################################################### -def find_min(arr): - min_val = arr[0] - min_index = 0 - for i in range(1, len(arr)): - if arr[i] < min_val: - min_val = arr[i] - min_index = i - return min_val, min_index - -################################# -def find_max(arr): - max_val = arr[0] - max_index = 0 - for i in range(1, len(arr)): - if arr[i] > max_val: - max_val = arr[i] - max_index = i - return max_val, max_index - -########################################################## - - -def LSR_velocity(ra, dec, velocity, Obstime): - """ - - - Parameters - ---------- - ra : TYPE - DESCRIPTION. - dec : TYPE - DESCRIPTION. - velocity : TYPE - DESCRIPTION. - Obstime : TYPE - DESCRIPTION. - - Returns - ------- - TYPE - DESCRIPTION. - - """ - # ''' - # 输入参数为ra，dec 原始视向速度in km/s,观测时间 - # local为观测地位置，c为利用输入坐标建立的坐标系对象。 - # barycorr为程序输出的地球坐标系与BSR间的修正，单位为m/s - # 修正后再将BSR中的速度转化为LSR中的速度。 - # ''' - local = EarthLocation.from_geodetic( - lat=22.349368*u.deg, lon=113.584068*u.deg, height=10*u.m) - # convert ra and dec to - source = SkyCoord(ra*u.deg, dec*u.deg, frame='icrs', - unit=(u.hourangle, u.deg)) - l = source.galactic.l.deg - b = source.galactic.b.deg - c = SkyCoord(l=l*u.degree, b=b*u.degree, frame='galactic') - c_icrs = c.transform_to('icrs') - barycorr = c_icrs.radial_velocity_correction( - obstime=Time(Obstime), location=local) - velocity = velocity + barycorr.value/1000 - # print(barycorr.value/1000) - l = l * np.pi / 180 - b = b * np.pi / 180 - return velocity + 9 * np.cos(l) * np.cos(b) + 12 * np.sin(l) * np.cos(b) + 7 * np.sin(b) - -############################################################################### - - -def rotation_yan(x0, y0, x1, y1, angle): - """ - - - Parameters - ---------- - x0 : TYPE - DESCRIPTION. - y0 : TYPE - DESCRIPTION. - x1 : TYPE - DESCRIPTION. - y1 : TYPE - DESCRIPTION. - angle : TYPE - DESCRIPTION. - - Returns - ------- - x2 : TYPE - DESCRIPTION. - y2 : TYPE - DESCRIPTION. - - """ - # % point A (x0,y0) - # % point B(x1,y1) - # % roate angle ,point B rotate with point A - alpha = angle/180*3.1415926 # % in radians - x2 = (x1-x0)*np.cos(alpha)-(y1-y0)*np.sin(alpha)+x0 - y2 = (x1-x0)*np.sin(alpha)+(y1-y0)*np.cos(alpha)+y0 - return x2, y2 -###################################################################### - - -def centroid(data): - """ - - - Parameters - ---------- - data : TYPE - DESCRIPTION. - - Returns - ------- - TYPE - DESCRIPTION. - - """ - h, w = np.shape(data) - x = np.arange(0, w) - y = np.arange(0, h) - x1 = np.ones((1, h)) - y1 = np.ones((w, 1)) - cx = (np.dot(np.dot(x1, data), y)) / (np.dot(np.dot(x1, data), y1)) - cy = (np.dot(np.dot(x, data), y1)) / (np.dot(np.dot(x1, data), y1)) - return float(cx), float(cy) - -############################################################################### - - -def SpecCube2photon(data, wave): - """ - - - Parameters - ---------- - data : TYPE - DESCRIPTION. - wave : TYPE - DESCRIPTION. - - Returns - ------- - Nphoton : TYPE - DESCRIPTION. - - """ - # calcutle photons from original spec cube data, - # data: input data, two-dimentional, unit : 1e-17 erg/s/A/cm^2 - # wave: the relative wavefront, unit in nm; - planckh = 6.62620*10**-27 # % erg s; - cc = 2.99792458*10**17 # in nm/s - # fov2=0.01 # 单个像素的面积 arcsec^2 - # telarea=3.1415926*100*100 # in cm^2，望远镜聚光面积； - fluxlam = 1e-17*data # convert original unit to unit of erg/s/A/cm^2 - lam = wave # wave in nm ;; - ephoton = planckh*cc/lam # in erg/photon, cc与lambda单位需要一致； - Nphoton = fluxlam/ephoton # in unit of photons/cm2/s/A - return Nphoton - -###################################################################### -###################################################################### - - -def mag2flux(starmag, lam): - """ - - - Parameters - ---------- - starmag : TYPE - DESCRIPTION. - lam : TYPE - DESCRIPTION. - - Returns - ------- - fluxlam : TYPE - DESCRIPTION. - photonu : TYPE - DESCRIPTION. - - """ - # starmag: input AB mag - # lambda: wavelength in nm - - planckh = 6.62620*10**-27 # % erg s; - cc = 2.99792458*10**17 # % in nm/s - - iskyv0 = starmag # mag/arcsec^2 unit - fluxnu = 10.0**((iskyv0+48.6)/-2.5) # erg/s/cm2/Hz/arcsec^2 - fluxlam = fluxnu*cc/lam**2/10 # erg/s/cm2/A/arcsec^2 - - Ee = planckh*cc/lam # single photon energy in erg/photon; - photonu = fluxlam/Ee # photons/cm2/s/A/arcsec^2 - - return fluxlam, photonu - -############################################################################### - - -def flux2photons(flux, lam): - """ - - - Parameters - ---------- - flux : TYPE - DESCRIPTION. - lam : TYPE - DESCRIPTION. - - Returns - ------- - photonu : TYPE - DESCRIPTION. - - """ - # flux: # erg/s/cm2/A/ - # lam: # wavelength in nm - - planckh = 6.62620*10**-27 # % erg s; - cc = 2.99792458*10**17 # % in nm/s - Ee = planckh*cc/lam # singe photon energy in erg/photon; - photonu = flux/Ee # photons/cm2/s/A - return photonu - -############################################################################### - - -def fftrange(n, dtype=float): - """ - - - Parameters - ---------- - n : TYPE - DESCRIPTION. - dtype : TYPE, optional - DESCRIPTION. The default is float. - - Returns - ------- - TYPE - DESCRIPTION. - - """ - """FFT-aligned coordinate grid for n samples.""" - return np.arange(-n//2, -n//2+n, dtype=dtype) - -############################################################################### - - -def dft2(ary, Q, samples, shift=None): - """ - - - Parameters - ---------- - ary : TYPE - DESCRIPTION. - Q : TYPE - DESCRIPTION. - samples : TYPE - DESCRIPTION. - shift : TYPE, optional - DESCRIPTION. The default is None. - - Returns - ------- - out : TYPE - DESCRIPTION. - - """ - """Compute the two dimensional Discrete Fourier Transform of a matrix. - - Parameters - ---------- - ary : `numpy.ndarray` - an array, 2D, real or complex. Not fftshifted. - Q : `float` - oversampling / padding factor to mimic an FFT. If Q=2, Nyquist sampled - samples : `int` or `Iterable` - number of samples in the output plane. - If an int, used for both dimensions. If an iterable, used for each dim - shift : `float`, optional - shift of the output domain, as a frequency. Same broadcast - rules apply as with samples. - - Returns - ------- - `numpy.ndarray` - 2D array containing the shifted transform. - Equivalent to ifftshift(fft2(fftshift(ary))) modulo output - sampling/grid differences - - """ - - # this is for dtype stabilization - Q = float(Q) # Q=lambda*f/D/pixelsize - - n, m = ary.shape # ary maybe is the pupil function - N, M = samples, samples - - X, Y, U, V = (fftrange(n) for n in (m, n, M, N)) - - a = 1 / Q - - ################################################################### - Eout_fwd = np.exp(-1j * 2 * np.pi * a / n * np.outer(Y, V).T) - Ein_fwd = np.exp(-1j * 2 * np.pi * a / m * np.outer(X, U)) - - ################################################################### - - out = Eout_fwd @ ary @ Ein_fwd # ary is the input pupil function - - return out -############################################################################## -############################################################################### - - -def anySampledPSFnew(wavefront, pupil, Q, sizeout): - """ - - - Parameters - ---------- - wavefront : TYPE - DESCRIPTION. - pupil : TYPE - DESCRIPTION. - Q : TYPE - DESCRIPTION. - sizeout : TYPE - DESCRIPTION. - - Returns - ------- - psf : TYPE - DESCRIPTION. - - """ - ''' - written on 2020.12.04 by Yan @Shao - % wavefront sampled in the united circle; - % pupil function samped as wavefront; - % sample ratio Q=lambda/D/pixel; - % lambda is wavelength; - % D is diameter; - % pixel is in arcsec; - % or sample ratio Q=lambda*f/D/pixelsize - % f is focal length; - % pixelsize is the actural size of the detector; - % make sure all the varia have the same unit; - % the returned PSF has sum value of 1 - ''' - - m, n = np.shape(wavefront) - - phase = 2.0*np.pi*wavefront # % wavefront phase in radians - - # %generalized pupil function of channel1; - pk = pupil*np.exp(cmath.sqrt(-1)*(phase)) - - pf = dft2(pk, Q, sizeout) - - pf2 = abs(pf)**2 - psfout = pf2/pf2.sum() - - cx, cy = centroid(psfout) - Nt = sizeout - - # for convolve method - psf = ndimage.shift( - psfout, [Nt/2-cy-1, Nt/2-cx-1], order=1, mode='nearest') - psf = psf/psf.sum() - - return psf - -############################################################################## - -############################################################################## - - -def get_dx_dy_blue(wave): - """ - - - Parameters - ---------- - wave : TYPE - DESCRIPTION. - - Returns - ------- - dx : TYPE - DESCRIPTION. - dy : TYPE - DESCRIPTION. - - """ - # wave is the wavelength in nm; - # dx is in dispersion direction, dy is in vertical direction; - # dydl = np.array([-423.256, 0.001, 0.00075, - # 0.0000078, -0.0000000000007, 0.0, 0.0]) - # # 色散方向 - # dxdl = 0.2*np.array([-9.1519, -1.00000000e-06, 3.50000000e-08, -5.00000000e-09, - # -1.70000000e-11, 4.00949787e-12, -6.16873452e-15]) - - ##### update @2024.10.16 - dydl=np.array([ 2447.9, -1/0.141, 0.0000075, 0.00000078, -0.000000000007] ) ; #色散方向 - dxdl=np.array([ 5.46, -1.5e-02, 3.5e-08, -5.0e-09]) #垂直方向 - - dx = 0.0 - dy = 0.0 - for i in range(len(dxdl)): - dx = dx+dxdl[i]*wave**(i) - dy = dy+dydl[i]*wave**(i) - return dx, dy -############################################################################## - - -def get_dx_dy_red(wave): - """ - - - Parameters - ---------- - wave : TYPE - DESCRIPTION. - - Returns - ------- - dx : TYPE - DESCRIPTION. - dy : TYPE - DESCRIPTION. - - """ - # wave is the wavelength in nm; - # dydl = np.array([-640.0239901372472, 0.0018, 0.00048, - # 0.0000028, -0.0000000000007, 0.0, 0.0]) # 色散方向 - # dxdl = 0.00325*np.array([-1638.8, 4.0e-2, 5.500e-3, - - # 5.2e-10, 1.7000e-10, 7.1e-13, -5.16e-15]) # 垂直方向 - ## update @2014.10.17 - dydl=np.array([3519.78622, -1/0.1555, 0.0000048, 0.00000028, -0.0000000000007] ) #色散方向 - dxdl=np.array([-5.6305, 1.0e-2, 5.500e-7, -5.2e-10]) #垂直方向 - - dx = 0.0 - dy = 0.0 - for i in range(len(dxdl)): - dx = dx+dxdl[i]*wave**(i) - dy = dy+dydl[i]*wave**(i) - return dx, dy - -############################################################################################## - - -def getSpectrum(Spectrum0, lam, sigma): - """ - - - Parameters - ---------- - Spectrum0 : TYPE - DESCRIPTION. - lam : TYPE - DESCRIPTION. - sigma : TYPE - DESCRIPTION. - - Returns - ------- - SpmOut : TYPE - DESCRIPTION. - - """ - - # %获得输入波长lambda的光谱强度 - wave = Spectrum0[:, 0] # %原始光谱数据给定的波长； - Qt = Spectrum0[:, 1] # %原始光谱数据给定的光谱强度； - d = abs(lam-wave) # %计算波长距离； - if min(d) > 3*sigma: - - SpmOut = 1.0/20000 - else: - - column = np.where(d == min(d)) - SpmOut = Qt[column[0]]*np.exp(-(lam-wave[column[0]])**2/2/sigma**2) - - return SpmOut - - -############################################################################## -############################################################################## - -class IFSsimulator(): - """ - CSST IFS Simulator - - The image that is being build is in:: - - self.image - self.image_b blue channel - self.image_r red channel - - :param opts: OptionParser instance - :type opts: OptionParser instance - """ - - def __init__(self, configfile): - """ - Class Constructor. - - :param opts: OptionParser instance - :type opts: OptionParser instance - - """ - - #################################### - - self.configfile = configfile - self.section = 'TEST' # simulation section; - - # load instrument model - self.information = IFSinformation() - - # update settings with defaults - self.information.update(dict(quadrant=int(0), - ccdx=int(0), - ccdy=int(0), - - ccdxgap=1.643, - ccdygap=8.116, - bluesize=4000, - redsize=6000, - # prescanx=50, - # ovrscanx=20, - fullwellcapacity=90000, - dark=0.001, - readout=4, - bias=500.0, - e_adu=1.5, - readouttime=0.13824, - rdose=8.0e9, - ra=0.0, - dec=0.0, - coveringfraction=0.1, - # CR: 3.0% is for 565s exposure - ######################################### - - mode='same')) - - ##################################### - - def readConfigs(self, simnumber): - """ - - - Parameters - ---------- - simnumber : TYPE - DESCRIPTION. - - Returns - ------- - None. - - """ - """ - Reads the config file information using configParser and sets up a logger. - """ - self.config = ConfigParser.RawConfigParser() - - self.config.read_file(open(self.configfile)) - -############################################################################### - - def processConfigs(self): - """ - Processes configuration information and save the information to a dictionary self.information. - - The configuration file may look as follows:: - - [TEST] - - - For explanation of each field, see /data/test.config. Note that if an input field does not exist, - then the values are taken from the default instrument model as described in - support.IFSinstrumentModel.VISinformation(). Any of the defaults can be overwritten by providing - a config file with a correct field name. - """ - # parse options and update the information dictionary - options = self.config.options(self.section) - settings = {} - for option in options: - try: - settings[option] = self.config.getint(self.section, option) - except ValueError: - try: - settings[option] = self.config.getfloat( - self.section, option) - except ValueError: - settings[option] = self.config.get(self.section, option) - - self.information.update(settings) - - # force gain to be float - self.information['e_adu'] = float(self.information['e_adu']) - - self.cosmicRays = self.config.getboolean(self.section, 'cosmicRays') - self.darknoise = self.config.getboolean(self.section, 'darknoise') - self.cosmetics = self.config.getboolean(self.section, 'cosmetics') - self.radiationDamage = self.config.getboolean( - self.section, 'radiationDamage') - - self.bleeding = self.config.getboolean(self.section, 'bleeding') - self.save_cosmicrays = self.config.getboolean( - self.section, 'save_cosmicrays') - self.sky_noise = self.config.getboolean(self.section, 'sky_noise') - - self.nonlinearity = self.config.getboolean( - self.section, 'nonlinearity') - - self.flatfieldM = self.config.getboolean( - self.section, 'flatfieldM') - - self.readoutNoise = self.config.getboolean( - self.section, 'readoutnoise') - - self.appbianpai = self.config.getboolean( - self.section, 'appbianpai') - - self.intscale = True - ###################################################################### - - self.information['variablePSF'] = False - - self.booleans = dict(nonlinearity=self.nonlinearity, - flatfieldM=self.flatfieldM, - - cosmicRays=self.cosmicRays, - darknoise=self.darknoise, - cosmetics=self.cosmetics, - radiationDamage=self.radiationDamage, - bleeding=self.bleeding, - sky_noise=self.sky_noise, - intscale=self.intscale, - save_cosmicrays=self.save_cosmicrays, - appbianpai=self.appbianpai) - - ############################################################################ - - def _createEmpty(self): - """ - Creates and empty array of a given x and y size full of zeros. - add g r i channel images; - - Creates lensing parameters; - - - """ - - self.imgtemp = dict() - - self.pixel = 0.1 # arcsec, pixel scale size; - - -############################################################################## -############################################################################## - - def zodiacal(self, ra, dec, time): - """ - - - Parameters - ---------- - ra : TYPE - DESCRIPTION. - dec : TYPE - DESCRIPTION. - time : TYPE - DESCRIPTION. - - Returns - ------- - wave_A : TYPE - DESCRIPTION. - spec_erg2 : TYPE - DESCRIPTION. - - """ - """ - For given RA, DEC and TIME, return the interpolated zodical - spectrum in Leinert-1998. - - :param ra: RA in unit of degree, ICRS frame - :param dec: DEC in unit of degree, ICRS frame - :param time: the specified string that in ISO format i.e., yyyy-mm-dd. - :return: - wave_A: wavelength of the zodical spectrum - spec_mjy: flux of the zodical spectrum, in unit of MJy/sr - spec_erg: flux of the zodical spectrum, in unit of erg/s/cm^2/A/sr - - """ - - # get solar position - dt = datetime.fromisoformat(time) - jd = time2jd(dt) - ## jd = julian.to_jd(dt, fmt='jd') - - t = Time(jd, format='jd', scale='utc') - - astro_sun = get_sun(t) - ra_sun, dec_sun = astro_sun.gcrs.ra.deg, astro_sun.gcrs.dec.deg - - radec_sun = SkyCoord(ra=ra_sun*u.degree, - dec=dec_sun*u.degree, frame='gcrs') - lb_sun = radec_sun.transform_to('geocentrictrueecliptic') - - # get offsets between the target and sun. - radec_obj = SkyCoord(ra=ra*u.degree, dec=dec*u.degree, frame='icrs') - lb_obj = radec_obj.transform_to('geocentrictrueecliptic') - - beta = abs(lb_obj.lat.degree) - lamda = abs(lb_obj.lon.degree - lb_sun.lon.degree) - - # interpolated zodical surface brightness at 0.5 um - zodi = pd.read_csv( - self.information['dir_path']+'IFS_inputdata/refs/zodi_map.dat', sep='\s+', header=None, comment='#') - beta_angle = np.array([0, 5, 10, 15, 20, 25, 30, 45, 60, 75]) - lamda_angle = np.array([0, 5, 10, 15, 20, 25, 30, 35, 40, 45, - 60, 75, 90, 105, 120, 135, 150, 165, 180]) - xx, yy = np.meshgrid(beta_angle, lamda_angle) - #xx, yy = np.meshgrid(beta_angle, lamda_angle,indexing='ij', sparse=True) - - f = interpolate.interp2d(xx, yy, zodi, kind='linear') - #f = interpolate.RegularGridInterpolator((xx, yy), zodi, method='linear') - - zodi_obj = f(beta, lamda) # - - # read the zodical spectrum in the ecliptic - cat_spec = pd.read_csv( - self.information['dir_path']+'IFS_inputdata/refs/solar_spec.dat', sep='\s+', header=None, comment='#') - wave = cat_spec[0].values # A - spec0 = cat_spec[1].values # - zodi_norm = 252 # - - spec = spec0 * (zodi_obj / zodi_norm) * 1e-8 # - - # convert to the commonly used unit of MJy/sr, erg/s/cm^2/A/sr - wave_A = wave # A - # spec_mjy = spec * 0.1 * wave_A**2 / 3e18 * 1e23 * 1e-6 # MJy/sr - spec_erg = spec * 0.1 # erg/s/cm^2/A/sr - # erg/s/cm^2/A/arcsec^2 - spec_erg2 = spec_erg / 4.25452e10 - - # self.zodiacal_wave=wave_A # in A - - # self.zodiacal_flux=spec_erg2 - - return wave_A, spec_erg2 - ########################################################################## - - # def smoothingWithChargeDiffusion(self, image, sigma=(0.32, 0.32)): - # """ - - # Parameters - # ---------- - # image : TYPE - # DESCRIPTION. - # sigma : TYPE, optional - # DESCRIPTION. The default is (0.32, 0.32). - - # Returns - # ------- - # TYPE - # DESCRIPTION. - - # """ - # """ - # Smooths a given image with a gaussian kernel with widths given as sigmas. - # This smoothing can be used to mimic charge diffusion within the CCD. - - # The default values are from Table 8-2 of CCD_273_Euclid_secification_1.0.130812.pdf converted - # to sigmas (FWHM / (2sqrt(2ln2)) and rounded up to the second decimal. - - # .. Note:: This method should not be called for the full image if the charge spreading - # has already been taken into account in the system PSF to avoid double counting. - - # :param image: image array which is smoothed with the kernel - # :type image: ndarray - # :param sigma: widths of the gaussian kernel that approximates the charge diffusion [0.32, 0.32]. - # :param sigma: tuple - - # :return: smoothed image array - # :rtype: ndarray - # """ - # return ndimage.filters.gaussian_filter(image, sigma) -############################################################################### - - def readCosmicRayInformation(self): - """ - - - Returns - ------- - None. - - """ - """ - Reads in the cosmic ray track information from two input files. - - Stores the information to a dictionary called cr. - """ - self.log.info('Reading in cosmic ray information from %s and %s' % (self.information['cosmicraylengths'], - self.information['cosmicraydistance'])) - # print(self.information) - - crLengths = np.loadtxt( - self.information['dir_path']+self.information['cosmicraylengths']) - crDists = np.loadtxt( - self.information['dir_path']+self.information['cosmicraydistance']) - - self.cr = dict(cr_u=crLengths[:, 0], cr_cdf=crLengths[:, 1], cr_cdfn=np.shape(crLengths)[0], - cr_v=crDists[:, 0], cr_cde=crDists[:, 1], cr_cden=np.shape(crDists)[0]) - -############################################################################### -############################################################################### - def configure(self, sourcein, dir_path, result_path, simnumber, debug, applyhole): - """ - - - Parameters - ---------- - simnumber : TYPE - DESCRIPTION. - - Returns - ------- - None. - - """ - """ - Configures the simulator with input information and creates and - empty array to which the final image will - be build on. - """ - self.readConfigs(simnumber) - - self.processConfigs() - - self._createEmpty() - - self.information['dir_path']=dir_path - - self.information['result_path']=result_path - - self.information['debug'] = debug - - self.source=sourcein - - self.simnumber=simnumber - - ############################################################ - - now = datetime.utcnow() - - result_day = now.strftime("%Y-%m-%d") - - - if self.source == 'LAMP': - if applyhole == 'yes': - ss = '_with_hole_' - else: - ss = '_no_hole_' - - else: - ss = '_' - - # if self.information['dir_path'] == '/nfsdata/share/simulation-unittest/ifs_sim/': - # self.result_path = self.information['dir_path'] + \ - # self.information['result_path']+'/'+self.source+ss+result_day - # else: - - # home_path = os.environ['HOME'] - - # if home_path == '/home/yan': - - # self.result_path = '../IFS_simData/'+self.source+ss+result_day - # else: - # self.result_path = '/data/ifspip/CCD_ima/'+self.source+ss+result_day - - - - self.result_path= self.information['result_path']+'/'+self.source+ss+result_day - print(self.information['result_path']) - - if os.path.isdir(self.result_path) == False: - os.mkdir(self.result_path) - os.mkdir(self.result_path+'/calibration_Data') - os.mkdir(self.result_path+'/log_file') - os.mkdir(self.result_path+'/original_Calibration') - os.mkdir(self.result_path+'/original_sky') - os.mkdir(self.result_path+'/shift_rot_sky') - os.mkdir(self.result_path+'/sky_Data') - - ############################################################################## - now = datetime.utcnow() - - data_time = now.strftime("%Y-%m-%d-%H-%M-%S") - ############################################################ - - self.log = lg(self.result_path+'/log_file/IFS_' + - self.source+'_'+data_time+'_Num_'+str(simnumber)+'.log') - - self.log.info('STARTING A NEW SIMULATION') - - self.log.info('The exposure order is %i ' % self.simnumber) - - ################################################################################ - # public system information - self.pixelscale = 0.1 # the pixel size is 0.1 arcsec - self.pixelsize = 15 # the pixel physical size is 15 micron - self.Fnum = 15.469875 # the f number= focal_length/D; - self.telD = 2 # tht telescope size is 2 meter, diamter size; - self.telarea = 3.1415926*(self.telD/2)**2 * \ - 100*100 # telescope square in cm^2 - self.fov2 = 0.01 # the pixel square - self.planckh = 6.62620*1e-27 # % erg s; - self.cc = 2.99792458*1e17 # in nm/s - self.light_FWHM = 0.175 - self.HgArsigma = self.light_FWHM/2.35 # sigma value of the light source - - # load system optical and CCD efficiency data - matfn0 = self.information['dir_path']+'IFS_inputdata/TotalQ200923.mat' - self.log.info('Optical optical efficiency file path is: %s' % (matfn0)) - da0 = sio.loadmat(matfn0) - # optical efficiency of blue channel - self.optical_blue_Q = da0['opticalQ_1'] - # optical efficiency of red channel - self.optical_red_Q = da0['opticalQ_2'] - self.CCD_Qe = da0['ccdQE'] # red channel - - # load all useful data; - # load wavefront data; - matfn2 = self.information['dir_path']+'IFS_inputdata/opd/opd_638nm.mat' - self.log.info('OPD0 file path is: %s' % (matfn2)) - da2 = sio.loadmat(matfn2) - opd0 = da2['opd'] # opd unit is in meter - self.opd0 = opd0*1e9 # convert opd to nm - self.pupil = abs(opd0) > 0.0 - - #################### - matfn2 = self.information['dir_path']+'IFS_inputdata/opd/opd1.fits' - da = fits.open(matfn2) - self.opd1 = da[0].data - self.log.info('OPD1 file path is: %s' % (matfn2)) - - matfn2 = self.information['dir_path']+'IFS_inputdata/opd/opd2.fits' - self.log.info('OPD2 file path is: %s' % (matfn2)) - da = fits.open(matfn2) - self.opd2 = da[0].data - - matfn2 = self.information['dir_path']+'IFS_inputdata/opd/opd3.fits' - self.log.info('OPD3 file path is: %s' % (matfn2)) - da = fits.open(matfn2) - self.opd3 = da[0].data - - ################################################################################ - # slice information; the slice is put in vertial direction; - # the silce long size is 64 pixels - slice_blue = dict() - slice_red = dict() - - randRedpos = np.array([0.52835362, 1.1105926, -0.81667794, 0.88860884, -2.78092636, - -0.15810022, -1.56726852, -0.71601855, -1.31768647, 1.73107581, - 0.4933876, 2.83673377, 0.22226286, -0.02634998, 0.35539383, - -0.91989574, -2.44856212, 0.91020484, -3.03097852, -1.11638071, - 1.28360669, -0.12521128, -0.3907698, 0.70183575, 1.00578099, - 1.67339662, 0.18067182, -0.56303075, 0.40959616, 1.45518379, - -0.93194744, 0.41492972]) - - randBluepos = np.array([0.97449008, -0.21371406, -1.62513338, -3.06938604, 1.72615283, - 0.73333374, 0.80923919, -0.9418576, -0.16806578, -1.04416631, - 2.20155068, -0.0900156, 0.07597706, 0.76208373, 0.29426592, - -0.89434703, 0.34017826, 1.16936499, 0.10977829, -1.31179539, - -0.50859787, -1.01891651, -0.95791581, -1.53018041, 0.88358827, - 0.07837641, -0.86157157, -1.18070438, 0.53970599, 1.4913733, - 2.10938775, 1.23213412]) - - ##################### - for i in range(32): - if i == 0: - slice_blue['py'] = np.zeros(32) - slice_blue['px'] = np.zeros(32) - - slice_red['py'] = np.zeros(32) - slice_red['px'] = np.zeros(32) - - if i < 16: - slice_blue['py'][i] = 50+randBluepos[i]*4 - slice_blue['px'][i] = 3.55/0.015*i+166.0 - - slice_red['py'][i] = 50+randRedpos[i]*4 - slice_red['px'][i] = 3.55/0.015*i+1190.0 - - else: - slice_blue['py'][i] = 50+250+randBluepos[i]*4 - slice_blue['px'][i] = 3.55/0.015*(i-16)+166.0+118 - - slice_red['py'][i] = 50+250+randRedpos[i]*4 - slice_red['px'][i] = 3.55/0.015*(i-16)+1190.0+118 - ####### - ###### flip the fringe up to down,down to up@2024.10.16 - self.slice_blue = dict() - self.slice_red = dict() - self.slice_blue['py'] =2000-slice_blue['py'] - self.slice_blue['px'] = slice_blue['px'] - - self.slice_red['py'] =3000-slice_red['py'] - self.slice_red['px'] =slice_red['px'] - - ####################################################################### - maskSlice = dict() - maskSlit = dict() - sizeout = 100 - for k in range(32): - maskSlice[str(k)] = np.zeros((sizeout, sizeout)) - maskSlice[str(k)][2*k+18:2*k+20, int(sizeout/2) - - 32:int(sizeout/2)+32] = 1 - - maskSlit[str(k)] = np.zeros((sizeout, sizeout)) - maskSlit[str(k)][2*k+18:2*k+20, int(sizeout/2) - - 37:int(sizeout/2)+37] = 1 - - self.maskSlice = maskSlice - self.maskSlit = maskSlit -###################################################################### - -############################################################################# - - def MakeFlatMatrix(self, img, seed): - """ - - - Parameters - ---------- - img : TYPE - DESCRIPTION. - seed : TYPE - DESCRIPTION. - - Returns - ------- - FlatMat : TYPE - DESCRIPTION. - - """ - #### - ysize, xsize = img.shape - np.random.seed(seed) - r1, r2, r3, r4 = np.random.random(4) - a1 = -0.5 + 0.2*r1 - a2 = -0.5 + 0.2*r2 - a3 = r3+5 - a4 = r4+5 - xmin, xmax, ymin, ymax = 0, xsize, 0, ysize - Flty, Fltx = np.mgrid[ymin:ymax, xmin:xmax] - np.random.seed(seed) - p1, p2, bg = np.random.poisson(1000, 3) - Fltz = 1e-6*(a1 * (Fltx-p1) ** 2 + a2 * (Flty-p2) - ** 2 - a3*Fltx - a4*Flty) + bg*20 - FlatMat = Fltz/np.mean(Fltz) - - return FlatMat - #################################################################### - - def getFrameTransferImg(self): - """ - - - Returns - ------- - None. - - """ - """ - get frame transfer image from original image of self.image_b and self.image_r - """ - overscan=320 - self.frame_b_4 = np.zeros((1024+overscan, 2048), dtype=float) - self.frame_b_3 = np.zeros((1024+overscan, 2048), dtype=float) - self.frame_b_2 = np.zeros((1024+overscan, 2048), dtype=float) - self.frame_b_1 = np.zeros((1024+overscan, 2048), dtype=float) - - self.frame_r_4 = np.zeros((1536+overscan, 3072), dtype=float) - self.frame_r_3 = np.zeros((1536+overscan, 3072), dtype=float) - self.frame_r_2 = np.zeros((1536+overscan, 3072), dtype=float) - self.frame_r_1 = np.zeros((1536+overscan, 3072), dtype=float) - - # OSH, part2 ,no change - b_2 = np.sum(self.image_b[0:1024, 0:2048], axis=0) - - # flip right to left of part1 - temp_b_1 = np.fliplr(self.image_b[0:1024, 2048:4096]) # left to right - b_1 = np.sum(temp_b_1, axis=0) - - # part 4, left to right and up to down - temp_b_4 = np.fliplr( - self.image_b[1024:1024+1024, 2048:2048+2048]) # left to right - b_4 = np.sum(temp_b_4, axis=0) - - # part 3,OSE, up to dwon - b_3 = np.sum(self.image_b[1024:1024+1024, 0:2048], axis=0) - #### - - for k in range(1024+overscan): - # part - self.frame_b_4[k, :] += b_4*0.09216/self.information['exptime'] - self.frame_b_3[k, :] += b_3*0.09216/self.information['exptime'] - self.frame_b_2[k, :] += b_2*0.09216/self.information['exptime'] - self.frame_b_1[k, :] += b_1*0.09216/self.information['exptime'] - - ################### - # red channle - # OSH, part2 ,no change - r_2 = np.sum(self.image_r[0:1536, 0:3072], axis=0) - - # # flip right to left of part1 - temp_r_1 = np.fliplr(self.image_r[0:1536, 3072:6144]) # left to right - r_1 = np.sum(temp_r_1, axis=0) - - # part 4, left to right and up to down - temp_r_4 = np.fliplr( - self.image_r[1536:1536+1536, 3072:6144]) # left to right - r_4 = np.sum(temp_r_4, axis=0) - - # part 3,OSE, up to dwon - r_3 = np.sum(self.image_r[1536:1536+1536, 0:3072], axis=0) - - for k in range(1536+overscan): - # part 4, OSH zone; - self.frame_r_4[k, :] += r_4*0.13824/self.information['exptime'] - self.frame_r_3[k, :] += r_3*0.13824/self.information['exptime'] - self.frame_r_2[k, :] += r_2*0.13824/self.information['exptime'] - self.frame_r_1[k, :] += r_1*0.13824/self.information['exptime'] - # end - -############################################################################### - def applyflatfield(self, simnumber): - """ - - - Parameters - ---------- - simnumber : TYPE - DESCRIPTION. - - Returns - ------- - None. - - """ - """ - Applies multiplicative flat field to emulate pixel-to-pixel non-uniformity. - - Because the pixel-to-pixel non-uniformity effect (i.e. multiplicative) flat fielding takes place - before CTI and other effects, the flat field file must be the same size as the pixels that see - the sky. - """ - # apply Flat field to image - ### - flat_b = self.MakeFlatMatrix(self.image_b, 100) - - flat_r = self.MakeFlatMatrix(self.image_r, 200) - - self.image_b *= flat_b - - self.image_r *= flat_r - - self.log.info('Applied flatfield to images.') - - return - -############################################################################### - -############################################################################### - def addCosmicRays(self, idk): - """ - - - Parameters - ---------- - idk : TYPE - DESCRIPTION. - - Returns - ------- - None. - - """ - """ - Add cosmic rays to the arrays based on a power-law intensity - distribution for tracks. - Cosmic ray properties (such as location and angle) are chosen from - random Uniform distribution. - For details, see the documentation for the cosmicrays class in the - support package. - """ - self.readCosmicRayInformation() - # to scale the number of cosmics with exposure time - self.cr['exptime'] = self.information['exptime'] - - # cosmic ray image - crImage_b = np.zeros((2048, 4096), dtype=float) - - crImage_r = np.zeros((3072, 6144), dtype=float) - - # cosmic ray instance - cosmics_b = cosmicrays(self.log, crImage_b, crInfo=self.cr) - cosmics_r = cosmicrays(self.log, crImage_r, crInfo=self.cr) - - CCD_cr_b = cosmics_b.addUpToFraction( - self.information['coveringfraction']*self.information['exptime']/300.0, limit=None) - CCD_cr_r = cosmics_r.addUpToFraction( - self.information['coveringfraction']*self.information['exptime']/300.0, limit=None) - - # paste the information - #self.image += CCD_cr - self.image_b += CCD_cr_b - self.image_r += CCD_cr_r - # count the covering factor - area_cr_b = np.count_nonzero(CCD_cr_b) - area_cr_r = np.count_nonzero(CCD_cr_r) - - #self.log.info('The cosmic ray covering factor is %i pixels ' % area_cr) - self.log.info( - 'The cosmic ray in blue channel covering factor is %i pixels ' % area_cr_b) - self.log.info( - 'The cosmic ray in red channel covering factor is %i pixels ' % area_cr_r) - - # if self.save_cosmicrays: - - # self.log.info('Saved the cosmicRays fits...') - # fits.writeto(self.result_path+'/calibration_Data/cosmicMap_B_SN_'+str( - # self.simnumber)+'_exp_'+str(idk)+'.fits', np.int32(CCD_cr_b), overwrite=True) - # fits.writeto(self.result_path+'/calibration_Data/cosmicMap_R_SN_'+str( - # self.simnumber)+'_exp_'+str(idk)+'.fits', np.int32(CCD_cr_r), overwrite=True) - - ########################################################## -######################################################################### - - -############################################################################## - - def applyDarkCurrent(self): - """ - - - Returns - ------- - None. - - """ - """ - Apply dark current. Scales the dark with the exposure time. - - Additionally saves the image without noise to a FITS file. - """ - - self.log.info('Added dark current to bule and red channel') - - - if self.information['dark1_b']>0.001 or self.information['dark1_b']>0.001: - self.log.error( - 'dark1_b value error, it shoub be in [0.0001, 0.001]!') - raise ValueError( - 'dark1_b value error, it shoub be in [0.0001, 0.001]!') - - if self.information['dark2_b']>0.001 or self.information['dark2_b']>0.001: - self.log.error( - 'dark2_b value error, it shoub be in [0.0001, 0.001]!') - raise ValueError( - 'dark2_b value error, it shoub be in [0.0001, 0.001]!') - - if self.information['dark3_b']>0.001 or self.information['dark3_b']>0.001: - self.log.error( - 'dark3_b value error, it shoub be in [0.0001, 0.001]!') - raise ValueError( - 'dark3_b value error, it shoub be in [0.0001, 0.001]!') - - if self.information['dark4_b']>0.001 or self.information['dark4_b']>0.001: - self.log.error( - 'dark4_b value error, it shoub be in [0.0001, 0.001]!') - raise ValueError( - 'dark4_b value error, it shoub be in [0.0001, 0.001]!') - - #################################################################### - if self.information['dark1_r']>0.001 or self.information['dark1_r']>0.001: - self.log.error( - 'dark1_r value error, it shoub be in [0.0001, 0.001]!') - raise ValueError( - 'dark1_r value error, it shoub be in [0.0001, 0.001]!') - - if self.information['dark2_r']>0.001 or self.information['dark2_r']>0.001: - self.log.error( - 'dark2_r value error, it shoub be in [0.0001, 0.001]!') - raise ValueError( - 'dark2_r value error, it shoub be in [0.0001, 0.001]!') - - if self.information['dark3_r']>0.001 or self.information['dark3_r']>0.001: - self.log.error( - 'dark3_r value error, it shoub be in [0.0001, 0.001]!') - raise ValueError( - 'dark3_r value error, it shoub be in [0.0001, 0.001]!') - - if self.information['dark4_r']>0.001 or self.information['dark4_r']>0.001: - self.log.error( - 'dark4_r value error, it shoub be in [0.0001, 0.001]!') - raise ValueError( - 'dark4_r value error, it shoub be in [0.0001, 0.001]!') - ########################################################################### - - # blue zone 4 - self.image_b[0:1024, 0:2048] += self.information['exptime'] * \ - self.information['dark2_b'] - - ########## zone 1 ################# - self.image_b[1024:2048, 0:2048] += self.information['exptime'] * \ - self.information['dark3_b'] - - ########## zone 3 ################### - self.image_b[0:1024, 2048:4096] += self.information['exptime'] * \ - self.information['dark1_b'] - - # zone 2 - self.image_b[1024:2048, 2048:4096] += self.information['exptime'] * \ - self.information['dark4_b'] - - # red zone 4 - self.image_r[0:1536, 0:3072] += self.information['exptime'] * \ - self.information['dark2_r'] - ########## zone 1 ################# - - self.image_r[1536:3712, 0:3072] += self.information['exptime'] * \ - self.information['dark3_r'] - - ########## zone 3 ################### - self.image_r[0:1536, 3072:6144] += self.information['exptime'] * \ - self.information['dark1_r'] - - # zone 2 - self.image_r[1536:3072, 3072:6144] += self.information['exptime'] * \ - self.information['dark4_r'] - - -############################################################################## -############################################################################## - - def applyPoissonNoise(self): - """ - - - Returns - ------- - None. - - """ - """ - Add Poisson noise to the image. - """ - - rounded = np.rint(self.image_b) # round to - # ugly workaround for multiple rounding operations... - residual = self.image_b.copy() - rounded - rounded[rounded < 0.0] = 0.0 - - np.random.seed() - self.image_b = np.random.poisson(rounded).astype(float) - self.log.info('Added Poisson noise on channel blue') - self.image_b += residual - - rounded = np.rint(self.image_r) # round to - # ugly workaround for multiple rounding operations... - residual = self.image_r.copy() - rounded - rounded[rounded < 0.0] = 0.0 - np.random.seed() - self.image_r = np.random.poisson(rounded).astype(float) - self.log.info('Added Poisson noise on channel red') - self.image_r += residual - - -############################################################################## - - def applyCosmetics(self): - """ - - - Returns - ------- - None. - - """ - """ - Apply cosmetic defects described in the input file. - - #Number of hot and dead pixels from MSSL/Euclid/TR/12003 Issue 2 Draft b - - .. Warning:: This method does not work if the input file has exactly one line. - """ - cosmetics = np.loadtxt( - self.information['dir_path']+self.information['cosmeticsfile_b']) - self.log.info('cosmeticsfile_b path is: %s' % - (self.information['cosmeticsfile_b'])) - - x = np.round(cosmetics[:, 0]).astype(int) - y = np.round(cosmetics[:, 1]).astype(int) - value = cosmetics[:, 2] - - # cosmetics_b=np.zeros((1344,9642)) - - # cosmetics_r=np.zeros((1856,13768)) - - self.log.info('Adding cosmetic defects to blue channel:') - for xc, yc, val in zip(x, y, value): - if 0 <= xc <= 1024 and 50 < yc % 2418 <= 2048+50: - #self.image[yc, xc] = val - self.image_b[xc, yc] = val - # cosmetics_b[xc,yc]=val - self.log.info('x=%i, y=%i, value=%f' % (xc, yc, val)) -############################################################################# - - cosmetics = np.loadtxt( - self.information['dir_path']+self.information['cosmeticsfile_r']) - - self.log.info('cosmeticsfile_r path is: %s' % - (self.information['cosmeticsfile_r'])) - - x = np.round(cosmetics[:, 0]).astype(int) - y = np.round(cosmetics[:, 1]).astype(int) - value = cosmetics[:, 2] - - self.log.info('Adding cosmetic defects to red channel:') - - for xc, yc, val in zip(x, y, value): - if 0 <= xc <= 1536 and 50 < yc % 3442 <= 3072+50: - #self.image[yc, xc] = val - self.image_r[xc, yc] = val - # cosmetics_r[yc,xc]=val - self.log.info('x=%i, y=%i, value=%f' % (xc, yc, val)) - - ####### save cosmetics image code #### - # if self.simnumber<2: - # #save cosmic ray image map - - # obsid=self.simnumber - # bluefile=self.result_path+'/calibration_Data/IFS_cosmetics'+'_'+str(obsid)+'_B.fits' - # redfile =self.result_path+'/calibration_Data/IFS_cosmetics'+'_'+str(obsid)+'_R.fits' - - # fits.writeto(bluefile, cosmetics_b, overwrite=True) - # fits.writeto(redfile, cosmetics_r, overwrite=True) - - # #output information to a FITS file - - # self.log.info('Save cosmetics images fits file to calibration_Data file ') -################################################################################ - - def applyRadiationDamage(self): - """ - - - Returns - ------- - None. - - """ - """ - Applies CDM03 radiation model to the image being constructed. - - .. seealso:: Class :`CDM03` - """ - - self.log.debug('Starting to apply radiation damage model...') - # at this point we can give fake data... - cti = CDM03bidir(self.information, [], log=self.log) - - # here we need the right input data - self.image_b = cti.applyRadiationDamage(self.image_b.copy( - ).transpose(), iquadrant=self.information['quadrant']).transpose() - self.log.info('Radiation damage added.') - - self.log.debug('Starting to apply radiation damage model...') - # at this point we can give fake data... - cti = CDM03bidir(self.information, [], log=self.log) - # here we need the right input data - self.image_r = cti.applyRadiationDamage(self.image_r.copy( - ).transpose(), iquadrant=self.information['quadrant']).transpose() - self.log.info('Radiation damage added.') -############################################################################## - - def applyNonlinearity(self): - """ - - - Returns - ------- - None. - - """ - """ - Applies a CCD273 non-linearity model to the image being constructed. - """ - - self.log.debug('Starting to apply non-linearity model...') - self.image_b = CCDnonLinearityModel( - self.image_b.copy()) - - self.log.info('Non-linearity effects included.') - - self.log.debug('Starting to apply non-linearity model...') - self.image_r = CCDnonLinearityModel( - self.image_r.copy()) - - self.log.info('Non-linearity effects included.') - ###################################################################### - - def applyReadoutNoise(self): - """ - - - Returns - ------- - None. - - """ - """ - Applies readout noise to the image being constructed. - - The noise is drawn from a Normal (Gaussian) distribution with - average=0.0 and std=readout noise. - """ - self.log.info('readnoise added in blue channel') - - if self.information['rn1_b']>10 or self.information['rn1_b']<3: - self.log.error( - 'rn1_b value error, it shoub be in [3, 10]!') - raise ValueError( - 'rn1_b value error, it shoub be in [3, 10]!') - - if self.information['rn2_b']>10 or self.information['rn2_b']<3: - self.log.error( - 'rn2_b value error, it shoub be in [3, 10]!') - raise ValueError( - 'rn2_b value error, it shoub be in [3, 10]!') - - if self.information['rn3_b']>10 or self.information['rn3_b']<3: - self.log.error( - 'rn3_b value error, it shoub be in [3, 10]!') - raise ValueError( - 'rn3_b value error, it shoub be in [3, 10]!') - - if self.information['rn4_b']>10 or self.information['rn4_b']<3: - self.log.error( - 'rn4_b value error, it shoub be in [3, 10]!') - raise ValueError( - 'rn4_b value error, it shoub be in [3, 10]!') - - #################################################################### - if self.information['rn1_r']>10 or self.information['rn1_r']<3: - self.log.error( - 'rn1_r value error, it shoub be in [3, 10]!') - raise ValueError( - 'rn1_r value error, it shoub be in [3, 10]!') - - if self.information['rn2_r']>10 or self.information['rn2_r']<3: - self.log.error( - 'rn2_r value error, it shoub be in [3, 10]!') - raise ValueError( - 'rn2_r value error, it shoub be in [3, 10]!') - - if self.information['rn3_r']>10 or self.information['rn3_r']<3: - self.log.error( - 'rn3_r value error, it shoub be in [3, 10]!') - raise ValueError( - 'rn3_r value error, it shoub be in [3, 10]!') - - if self.information['rn4_r']>10 or self.information['rn4_r']<3: - self.log.error( - 'rn4_r value error, it shoub be in [3, 10]!') - raise ValueError( - 'rn4_r value error, it shoub be in [3, 10]!') - ########################################################33 - - - - # blue zone 1 - np.random.seed() - prescan = 50 - overscan = 320 - # zone 4 , OSH #### 1024, 2048 - self.image_b[0:1024+overscan, 0:2048+prescan+overscan] += np.random.normal( - loc=0.0, scale=self.information['rn2_b'], size=(1344, 2418)) - - ########## zone 3, OSG ################# - np.random.seed() - self.image_b[0:1024+overscan, 2418:2418+2048+prescan+overscan] += np.random.normal( - loc=0.0, scale=self.information['rn1_b'], size=(1344, 2418)) - - ########## zone 1, OSE ################### - np.random.seed() - self.image_b[0:1024+overscan, 2418*2:2418*2+2048+prescan+overscan] += np.random.normal( - loc=0.0, scale=self.information['rn3_b'], size=(1344, 2418)) - - ########## zone 2, OSF ############### - np.random.seed() - self.image_b[0:1024+overscan, 2418*3:2418*3+2048+prescan+overscan] += np.random.normal( - loc=0.0, scale=self.information['rn4_b'], size=(1344, 2418)) - - self.log.info('readnoise added in blue channel') - - # red zone 4 , OSH### 1536* 3072 - np.random.seed() - self.image_r[0:1536+overscan, 0:3072+prescan+overscan] += np.random.normal( - loc=0.0, scale=self.information['rn2_r'], size=(1856, 3442)) - - ########## zone 3 ,OSG ################# - np.random.seed() - self.image_r[0:1536+overscan, 3442:3442+3072+prescan+overscan] += np.random.normal( - loc=0.0, scale=self.information['rn1_r'], size=(1856, 3442)) - - ########## zone 1 ,OSE ################### - np.random.seed() - self.image_r[0:1536+overscan, 3442*2:3442*2+3072+prescan+overscan] += np.random.normal( - loc=0.0, scale=self.information['rn3_r'], size=(1856, 3442)) - - ########## zone 2,OSF ########### - np.random.seed() - self.image_r[0:1536+overscan, 3442*3:3442*3+3072+prescan+overscan] += np.random.normal( - loc=0.0, scale=self.information['rn4_r'], size=(1856, 3442)) - -########################################################################################## - - def appFrameTransferEffect(self): - """ - - - Returns - ------- - None. - - """ - """ - apply frame transfer effect to the four part images - """ - prescan = 50 - overscan = 0 - - # zone 1 , OSG #### 1024, 2048 - self.image_b[0:1024+overscan, prescan:2048 + - prescan] += self.frame_b_1[0:1024+overscan, 0:2048] - - ########## zone 2, OSH ################# - - self.image_b[0:1024+overscan, prescan+2418:2418+2048 + - prescan] += self.frame_b_2[0:1024+overscan, 0:2048] - - ########## zone 3, OSE ################### - - self.image_b[0:1024+overscan, prescan+2418*2:2418*2 + - 2048+prescan] += self.frame_b_3[0:1024+overscan, 0:2048] - - ########## zone 4, OSF ############### - - self.image_b[0:1024+overscan, prescan+2418*3:2418*3 + - 2048+prescan] += self.frame_b_4[0:1024+overscan, 0:2048] - - # 3 - - self.log.info('readnoise added in blue channel') - - # red zone 1 , OSG### 1536* 3072 - - self.image_r[0:1536+overscan, prescan+3442*0:3442*0 + - 3072+prescan] += self.frame_r_1[0:1536+overscan, 0:3072] - - ########## zone 2 ,OSh ################# - - self.image_r[0:1536+overscan, prescan+3442*1:3442*1 + - 3072+prescan] += self.frame_r_2[0:1536+overscan, 0:3072] - - ########## zone 1 ,OSE ################### - self.image_r[0:1536+overscan, prescan+3442*2:3442*2 + - 3072+prescan] += self.frame_r_3[0:1536+overscan, 0:3072] - - ########## zone 4,OSF ########### - - self.image_r[0:1536+overscan, prescan+3442*3:3442*3 + - 3072+prescan] += self.frame_r_4[0:1536+overscan, 0:3072] - - ########################################################################## - - def electrons2ADU(self): - """ - - - Returns - ------- - None. - - """ - """ - Convert from electrons to ADUs using the value read from the configuration file. - """ - if self.information['gain1_b']>2 or self.information['gain1_b']<1: - self.log.error( - 'gain1_b value error, it shoub be in [1, 2]!') - raise ValueError( - 'gain1_b value error, it shoub be in [1, 2]!') - - if self.information['gain2_b']>2 or self.information['gain2_b']<1: - self.log.error( - 'gain2_b value error, it shoub be in [1, 2]!') - raise ValueError( - 'gain2_b value error, it shoub be in [1, 2]!') - - if self.information['gain3_b']>2 or self.information['gain3_b']<1: - self.log.error( - 'gain3_b value error, it shoub be in [1, 2]!') - raise ValueError( - 'gain3_b value error, it shoub be in [1, 2]!') - - if self.information['gain4_b']>2 or self.information['gain4_b']<1: - self.log.error( - 'gain4_b value error, it shoub be in [1, 2]!') - raise ValueError( - 'gain4_b value error, it shoub be in [1, 2]!') - - #################################################################### - if self.information['gain1_r']>2 or self.information['gain1_r']<1: - self.log.error( - 'gain1_r value error, it shoub be in [1, 2]!') - raise ValueError( - 'gain1_r value error, it shoub be in [1, 2]!') - - if self.information['gain2_r']>2 or self.information['gain2_r']<1: - self.log.error( - 'gain2_r value error, it shoub be in [1, 2]!') - raise ValueError( - 'gain2_r value error, it shoub be in [1, 2]!') - - if self.information['gain3_r']>2 or self.information['gain3_r']<1: - self.log.error( - 'gain3_r value error, it shoub be in [1, 2]!') - raise ValueError( - 'gain3_r value error, it shoub be in [1, 2]!') - - if self.information['gain4_r']>2 or self.information['gain4_r']<1: - self.log.error( - 'gain4_r value error, it shoub be in [1, 2]!') - raise ValueError('gain4_r value error, it shoub be in [1, 2]!') - #################################################################### - - - self.log.info( - 'Converting from electrons to ADUs using a factor of gain') - - # # blue zone 4 - # self.image_b[0:1344, 0:2418] /= self.information['gain4_b'] - - # ########## zone 3 ################# - # self.image_b[0:1344, 2418:2418*2] /= self.information['gain3_b'] - - # ########## zone 1 ################### - # self.image_b[0:1344, 2418*2:2418*3] /= self.information['gain1_b'] - - # # zone 2 - # self.image_b[0:1344, 2418*3:2418*4] /= self.information['gain2_b'] - -################# update @2024.10.18 ### - # first part, menas old blue zone 4 - self.image_b[0:1344, 0:2418] /= self.information['gain1_b'] - - ##########second part, means old zone 3 ################# - self.image_b[0:1344, 2418:2418*2] /= self.information['gain2_b'] - - ##########third part, means old zone 1 ################### - self.image_b[0:1344, 2418*2:2418*3] /= self.information['gain3_b'] - - #### fourth part, means old zone 2 - self.image_b[0:1344, 2418*3:2418*4] /= self.information['gain4_b'] - ############################################################################ - -########################## - # red zone 4 - self.image_r[0:1856, 0:3442] /= self.information['gain1_r'] - - ########## zone 3 ################# - self.image_r[0:1856, 3442:3442*2] /= self.information['gain2_r'] - - ########## zone 1 ################### - self.image_r[0:1856, 3442*2:3442*3] /= self.information['gain3_r'] - - # zone 2 - self.image_r[0:1856, 3442*3:3442*4] /= self.information['gain4_r'] - - # 3 - - #################################################################################### - - def applyBias(self): - """ - - - Returns - ------- - None. - - """ - """ - Adds a bias level to the image being constructed. - - The value of bias is read from the configure file and stored - in the information dictionary (key bias). - - """ - - if self.information['bias1_b']>2000 or self.information['bias1_b']<100: - self.log.error( - 'bias1_b value error, it shoub be in [100, 2000]!') - raise ValueError( - 'bias1_b value error, it shoub be in [100, 2000]!') - - if self.information['bias2_b']>2000 or self.information['bias2_b']<100: - self.log.error( - 'bias2_b value error, it shoub be in [100, 2000]!') - raise ValueError( - 'bias2_b value error, it shoub be in [100, 2000]!') - - if self.information['bias3_b']>2000 or self.information['bias3_b']<100: - self.log.error( - 'bias3_b value error, it shoub be in [100, 2000]!') - raise ValueError( - 'bias3_b value error, it shoub be in [100, 2000]!') - - if self.information['bias4_b']>2000 or self.information['bias4_b']<100: - self.log.error( - 'bias4_b value error, it shoub be in [100, 2000]!') - raise ValueError( - 'bias4_b value error, it shoub be in [100, 2000]!') - - #################################################################### - if self.information['bias1_r']>2000 or self.information['bias1_r']<100: - self.log.error( - 'bias1_r value error, it shoub be in [100, 2000]!') - raise ValueError( - 'bias1_r value error, it shoub be in [100, 2000]!') - - if self.information['bias2_r']>2000 or self.information['bias2_r']<100: - self.log.error( - 'bias2_r value error, it shoub be in [100, 2000]!') - raise ValueError( - 'bias2_r value error, it shoub be in [100, 2000]!') - - if self.information['bias3_r']>2000 or self.information['bias3_r']<100: - self.log.error( - 'bias3_r value error, it shoub be in [100, 2000]!') - raise ValueError( - 'bias3_r value error, it shoub be in [100, 2000]!') - - if self.information['bias4_r']>2000 or self.information['bias4_r']<100: - self.log.error( - 'bias4_r value error, it shoub be in [100, 2000]!') - raise ValueError( - 'bias4_r value error, it shoub be in [100, 2000]!') - - ######################################################################## - - # blue zone 4 - self.image_b[0:1344, 0:2418] += self.information['bias2_b'] - - ########## zone 3 ################# - self.image_b[0:1344, 2418:2418*2] += self.information['bias1_b'] - - ########## zone 1 ################### - self.image_b[0:1344, 2418*2:2418*3] += self.information['bias3_b'] - - # zone 2 - self.image_b[0:1344, 2418*3:2418*4] += self.information['bias4_b'] - - ############################################################################ - - self.log.info('Bias counts were added to the blue image') - - ####################################################################### - - # red zone 4 - self.image_r[0:1856, 0:3442] += self.information['bias2_r'] - - ########## zone 3 ################# - self.image_r[0:1856, 3442:3442*2] += self.information['bias1_r'] - - ########## zone 1 ################### - self.image_r[0:1856, 3442*2:3442*3] += self.information['bias3_r'] - - # zone 2 - self.image_r[0:1856, 3442*3:3442*4] += self.information['bias4_r'] - - ####################################################################### - - self.log.info('Bias counts were added to the red image') - -############################################################################### - def addPreOverScans(self): - """ - - - Returns - ------- - None. - - """ - """ - Add pre- and overscan regions to the self.image. These areas are added only in the serial direction. - Because the 1st and 3rd quadrant are read out in to a different serial direction than the nominal - orientation, in these images the regions are mirrored. - - The size of prescan and overscan regions are defined by the prescanx and overscanx keywords, respectively. - """ - self.log.info('Adding pre- and overscan regions') - - canvas = np.zeros((self.information['bluesize'], - (self.information['bluesize'] + self.information['prescanx'] + self.information['ovrscanx']))) - - # because the pre- and overscans are in x-direction this needs to be taken into account for the - # 1st and 3rd quadrant - if self.information['quadrant'] in (0, 2): - canvas[:, self.information['prescanx']: self.information['prescanx'] + - self.information['bluesize']] = self.image_b - self.image_b = canvas - elif self.information['quadrant'] in (1, 3): - canvas[:, self.information['ovrscanx']: self.information['ovrscanx'] + - self.information['bluesize']] = self.image_b - self.image_b = canvas - else: - self.log.error( - 'Cannot include pre- and overscan because of an unknown quadrant!') - - ##################################################################### - - canvas = np.zeros((self.information['redsize'], - (self.information['redsize'] + self.information['prescanx'] + self.information['ovrscanx']))) - - # because the pre- and overscans are in x-direction this needs to be taken into account for the - # 1st and 3rd quadrant - if self.information['quadrant'] in (0, 2): - canvas[:, self.information['prescanx']: self.information['prescanx'] + - self.information['redsize']] = self.image_r - self.image_r = canvas - elif self.information['quadrant'] in (1, 3): - canvas[:, self.information['ovrscanx']: self.information['ovrscanx'] + - self.information['redsize']] = self.image_r - self.image_r = canvas - else: - self.log.error( - 'Cannot include pre- and overscan because of an unknown quadrant!') - -############################################################################### - - def applyBleeding_yan(self): - """ - - - Returns - ------- - None. - - """ - """ - Apply bleeding along the CCD columns if the number of electrons in - a pixel exceeds the full-well capacity. - - Bleeding is modelled in the parallel direction only, because the - CCD273s are assumed not to bleed in - serial direction. - - :return: None - """ - self.log.info('Applying column bleeding...') - - # loop over each column, as bleeding is modelled column-wise - for i, column in enumerate(self.image_b.T): - sum = 0. - for j, value in enumerate(column): - # first round - from bottom to top (need to half the bleeding) - overload = value - self.information['fullwellcapacity'] - if overload > 0.: - overload /= 2. - #self.image[j, i] -= overload - self.image_b[j, i] -= overload - - sum += overload - - elif sum > 0.: - if -overload > sum: - overload = -sum - - self.image_b[j, i] -= overload - sum += overload - - for i, column in enumerate(self.image_b.T): - sum = 0. - for j, value in enumerate(column[::-1]): - # second round - from top to bottom (bleeding was half'd already, so now full) - overload = value - self.information['fullwellcapacity'] - if overload > 0.: - #self.image[-j-1, i] -= overload - self.image_b[-j-1, i] -= overload - - sum += overload - elif sum > 0.: - if -overload > sum: - overload = -sum - #self.image[-j-1, i] -= overload - self.image_b[-j-1, i] -= overload - - sum += overload - print('Applying column bleeding to blue image finished.......') - # 3333 - for i, column in enumerate(self.image_r.T): - sum = 0. - for j, value in enumerate(column): - # first round - from bottom to top (need to half the bleeding) - overload = value - self.information['fullwellcapacity'] - if overload > 0.: - overload /= 2. - #self.image[j, i] -= overload - self.image_r[j, i] -= overload - - sum += overload - elif sum > 0.: - if -overload > sum: - overload = -sum - #self.image[j, i] -= overload - self.image_r[j, i] -= overload - - sum += overload - - for i, column in enumerate(self.image_r.T): - sum = 0. - for j, value in enumerate(column[::-1]): - # second round - from top to bottom (bleeding was half'd already, so now full) - overload = value - self.information['fullwellcapacity'] - if overload > 0.: - #self.image[-j-1, i] -= overload - self.image_r[-j-1, i] -= overload - - sum += overload - elif sum > 0.: - if -overload > sum: - overload = -sum - #self.image[-j-1, i] -= overload - self.image_r[-j-1, i] -= overload - - sum += overload - print('Applying column bleeding to red image finished.......') - - ########################################################################### - - ############################################################################ - - def discretise(self, max=2**16-1): - """ - - - Parameters - ---------- - max : TYPE, optional - DESCRIPTION. The default is 2**16-1. - - Returns - ------- - None. - - """ - """ - Converts a floating point image array (self.image) to an integer array with max values - defined by the argument max. - - :param max: maximum value the the integer array may contain [default 65k] - :type max: float - - :return: None - """ - # avoid negative numbers in case bias level was not added - self.image_b[self.image_b < 0.0] = 0. - # cut of the values larger than max - self.image_b[self.image_b > max] = max - - self.image_b = np.rint(self.image_b).astype(int) - self.log.info('Maximum and total values of the image are %i and %i, respectively' % (np.max(self.image_b), - np.sum(self.image_b))) - - # avoid negative numbers in case bias level was not added - self.image_r[self.image_r < 0.0] = 0. - # cut of the values larger than max - self.image_r[self.image_r > max] = max - - self.image_r = np.rint(self.image_r).astype(int) - self.log.info('Maximum and total values of the image are %i and %i, respectively' % (np.max(self.image_r), - np.sum(self.image_r))) - - ################################################################################################## - - # def applyImageShift(self): - # """ - - - # Returns - # ------- - # None. - - # """ - - # np.random.seed(9*self.simnumber) - # ud = np.random.random() # Choose a random rotation - # dx = 2 * (ud-0.5) * self.information['shiftmax'] - - # np.random.seed(99*self.simnumber) - # ud = np.random.random() # Choose a random rotation - # dy = 2 * (ud-0.5) * self.information['shiftmax'] - - # self.image_b = ndimage.shift(self.image_b.copy(), [ - # dy+self.information['shift_b_y'], dx+self.information['shift_b_x']], order=0, mode='nearest') - # self.image_r = ndimage.shift(self.image_r.copy(), [ - # dy+self.information['shift_r_y'], dx+self.information['shift_r_x']], order=0, mode='nearest') - - # self.log.info('Applied image shifting to g r i channels.') - # self.information['ra'] = dx*self.information['pixel_size'] - # self.information['dec'] = dy*self.information['pixel_size'] - - -############################################################################## - - # def applyImageRotate(self): - # """ - - # Returns - # ------- - # None. - - # """ - # np.random.seed(10*self.simnumber) - # ud = np.random.random() # Choose a random rotation - # angle = 2 * (ud-0.5) * self.information['tel_rotmax'] - - # inputimg = self.image_b.copy() - # # here we choose reshape=False, the rotated image will - # rotimg = ndimage.rotate( - # inputimg, angle+self.information['rotate_b'], order=1, reshape=False) - # self.image_b = rotimg - - # inputimg = self.image_r.copy() - # # here we choose reshape=False, the rotated image will - # rotimg = ndimage.rotate( - # inputimg, angle+self.information['rotate_r'], order=1, reshape=False) - # self.image_r = rotimg - - # self.information['Tel_rot'] = angle - - # self.log.info( - # 'Applied telescope rotation with angle (in degree)= %f.', angle) -############################################################################### - - - def CCDreadout(self): - """ - - - Returns - ------- - None. - - """ - - imgb = self.image_b.copy() - prescan = int(self.information['prescan']) - overscan = int(self.information['overscan']) - - temp = np.zeros((1344, 9672)) - # # zone 4, OSH ######0:1024, 50:2048+50 - # x1 = 0 - # x2 = x1+1024 - - # y1 = 0+prescan - # y2 = y1+2048 - # temp[x1:x2, y1:y2] = imgb[0:1024, 0:2048] - # # zone 3, OSG , left to right ################# - # # np.fliplr(b2) ## left to right - # # np.flipud(b3) ## down to up - # x1 = 0 - # x2 = x1+1024 - - # y1 = 2418+prescan - # y2 = y1+2048 - # temp[x1:x2, y1:y2] = np.fliplr(imgb[0:1024, 2048:4096]) - # # zone 1, OSE,down to up ################### - # x1 = 0 - # x2 = x1+1024 - - # y1 = 2418*2+prescan - # y2 = y1+2048 - # temp[x1:x2, y1:y2] = np.flipud(imgb[1024:2048, 0:2048]) - # ########## zone 2, OSF down to yp ,left to right ####### - # x1 = 0 - # x2 = x1+1024 - - # y1 = 2418*3+prescan - # y2 = y1+2048 - # temp[x1:x2, y1:y2] = np.flipud(np.fliplr(imgb[1024:2048, 2048:4096])) - - ### update 2024.10.18 - # first part, old OSG part ,## shift: left to right - x1 = 0 - x2 = x1+1024 - - y1 = 0+prescan - y2 = y1+2048 - temp[x1:x2, y1:y2] = np.fliplr(imgb[0:1024, 2048:4096]) - # second part, old OSH, no shift ################# - # np.fliplr(b2) ## left to right - # np.flipud(b3) ## down to up - x1 = 0 - x2 = x1+1024 - - y1 = 2418+prescan - y2 = y1+2048 - temp[x1:x2, y1:y2] = imgb[0:1024, 0:2048] - ### third part, old OSE, down to up ################### - x1 = 0 - x2 = x1+1024 - - y1 = 2418*2+prescan - y2 = y1+2048 - temp[x1:x2, y1:y2] = np.flipud(imgb[1024:2048, 0:2048]) - ########## fourth part, old OSF part; down to yp ,left to right ####### - x1 = 0 - x2 = x1+1024 - - y1 = 2418*3+prescan - y2 = y1+2048 - temp[x1:x2, y1:y2] = np.flipud(np.fliplr(imgb[1024:2048, 2048:4096])) - - - self.image_b = temp - - ####################################################################### - imgr = self.image_r.copy() - temp = np.zeros((1856, 13768)) - - # # zone 4, OSH ######0:1024, 50:2048+50 - # x1 = 0 - # x2 = x1+1536 - - # y1 = 0+prescan - # y2 = y1+3072 - # temp[x1:x2, y1:y2] = imgr[0:1536, 0:3072] - # ########## zone 3, OSG , left to right ################# - # # np.fliplr(b2) ## left to right - # # np.flipud(b3) ## down to up - # x1 = 0 - # x2 = x1+1536 - - # y1 = 3442+prescan - # y2 = y1+3072 - # temp[x1:x2, y1:y2] = np.fliplr(imgr[0:1536, 3072:6144]) - # ########## zone 1, OSE,down to up ############################ - # x1 = 0 - # x2 = x1+1536 - - # y1 = 3442*2+prescan - # y2 = y1+3072 - # temp[x1:x2, y1:y2] = np.flipud(imgr[1536:3072, 0:3072]) - # ########## zone 2, OSF down to yp ,left to right ################ - # x1 = 0 - # x2 = x1+1536 - - # y1 = 3442*3+prescan - # y2 = y1+3072 - # temp[x1:x2, y1:y2] = np.flipud(np.fliplr(imgr[1536:3072, 3072:6144])) - - ###### update @2024.10.18 - # readout image ,first part, old OSG, shift: left to right - x1 = 0 - x2 = x1+1536 - - y1 = 0+prescan - y2 = y1+3072 - temp[x1:x2, y1:y2] = np.fliplr(imgr[0:1536, 3072:6144]) - ########## second part, old OSH ,no change; ################# - # np.fliplr(b2) ## left to right - # np.flipud(b3) ## down to up - x1 = 0 - x2 = x1+1536 - - y1 = 3442+prescan - y2 = y1+3072 - temp[x1:x2, y1:y2] = imgr[0:1536, 0:3072] - ########## third part , old OSE , down to up ############################ - x1 = 0 - x2 = x1+1536 - - y1 = 3442*2+prescan - y2 = y1+3072 - temp[x1:x2, y1:y2] = np.flipud(imgr[1536:3072, 0:3072]) - ########## fourth part, old OSF, down to up ,left to right ################ - x1 = 0 - x2 = x1+1536 - - y1 = 3442*3+prescan - y2 = y1+3072 - temp[x1:x2, y1:y2] = np.flipud(np.fliplr(imgr[1536:3072, 3072:6144])) - - - self.image_r = temp - - return -############################################################################## - - def writeOutputs(self, obnum): - """ - - - Parameters - ---------- - obnum : TYPE - DESCRIPTION. - - Returns - ------- - None. - - """ - """ - Writes out a FITS file using PyFITS and converts the image array to 16bit unsigned integer as - appropriate for VIS. - - Updates header with the input values and flags used during simulation. - """ - HeaderTest = 'no' - sim_ver = str(self.information['sim_ver']) - - # Readout information - now = self.dt.utcnow() - - # exposre start time - exp_start_utc = now - exp_start_str = now.strftime("%Y%m%d%H%M%S") - - # exsposure end time,including mulipe exposure times and readout times; - # but minus one readout_time - tt = (self.information['exptime']+self.information['readouttime']) * \ - self.information['exposuretimes']-self.information['readouttime'] - end = self.dt.utcnow()+timedelta(seconds=tt) - - exp_end_str = end.strftime("%Y%m%d%H%M%S") - exp_end_utc = end - - # write the actual file - tt = (self.information['exptime']+self.information['readouttime'] - )*self.information['exposuretimes'] - write_end = self.dt.utcnow()+timedelta(seconds=tt) - - # write_time_utc = write_end.strftime("%Y-%m-%dT%H:%M:%S") - # write_time_str = end.strftime("%Y-%m-%dT%H:%M:%S") - - if self.source == 'SCI' or self.source == 'COMP': - if self.source == 'SCI': - obsid = 30100000000+obnum # 3+OBSTYPE+8位曝光编号； - else: - obsid = 30200000000+obnum # 3+OBSTYPE+8位曝光编号； - - filename_b = 'CSST_IFS_B_'+self.source+'_'+exp_start_str + \ - '_'+exp_end_str+'_'+str(obsid)+'_B_L0_V'+sim_ver - self.file_b = self.result_path+'/sky_Data/'+filename_b+'.fits' - - filename_r = 'CSST_IFS_R_'+self.source+'_'+exp_start_str + \ - '_'+exp_end_str+'_'+str(obsid)+'_R_L0_V'+sim_ver - self.file_r = self.result_path+'/sky_Data/'+filename_r + '.fits' - - else: - - if self.source == 'BIAS': - obsid = 30300000000+obnum # 3+OBSTYPE+8位曝光编号； - if self.source == 'DARK': - obsid = 30400000000+obnum # 3+OBSTYPE+8位曝光编号； - if self.source == 'FLAT': - obsid = 30500000000+obnum # 3+OBSTYPE+8位曝光编号； - if self.source == 'LAMP': - obsid = 30600000000+obnum # 3+OBSTYPE+8位曝光编号； - - filename_b = 'CSST_IFS_B_'+self.source+'_'+exp_start_str + \ - '_'+exp_end_str+'_'+str(obsid)+'_B_L0_V'+sim_ver - self.file_b = self.result_path+'/calibration_Data/'+filename_b+'.fits' - - filename_r = 'CSST_IFS_R_'+self.source+'_'+exp_start_str + \ - '_'+exp_end_str+'_'+str(obsid)+'_R_L0_V'+sim_ver - self.file_r = self.result_path+'/calibration_Data/'+filename_r+'.fits' - - # create a new FITS file, using HDUList instance - ##### layer 0 #### - - ofd_b = fits.PrimaryHDU() - - ofd_b.header['SIMPLE'] = (True, 'conforms to FITS standard') - ofd_b.header['BITPIX'] = (0, 'array data type') - ofd_b.header['NAXIS'] = (0, 'number of array dimensions') - ofd_b.header['NEXTEND'] = (np.int16(1), 'number of array dimensions') - - ############################################################################################ - temp = write_end.utcnow() - data_time = temp.strftime("%Y-%m-%dT%H:%M:%S.%f") - ofd_b.header['DATE'] = ( - data_time[:21], 'written date (yyyy-mm-ddThh:mm:ss.s)') - ofd_b.header['FILENAME'] = (filename_b[:68], '') - ofd_b.header['FILETYPE'] = (self.source[:12], 'observation type') - ofd_b.header['TELESCOP'] = ('CSST', 'telescope name') - ofd_b.header['INSTRUME'] = ('IFS', 'instrument name') - ofd_b.header['RADECSYS'] = ( - 'ICRS', 'coordinate system of the object') - ofd_b.header['EQUINOX'] = (float(2000.0), '') - ofd_b.header['FITSSWV'] = ( - 'ifs_sim_0.8.03', 'FITS creating software version') - ######### Object information ############# - if self.source == 'SCI' or self.source == 'COMP': - ofd_b.header['OBJECT'] = ( - self.information['name_obj'][:30], 'object name') - ofd_b.header['TARGET'] = ( - (self.information['target']), 'target name (hhmmss.s+ddmmss)') - - ofd_b.header['OBSID'] = (str(obsid), 'observation ID') - - ofd_b.header['RA_OBJ'] = (np.float64( - self.information['ra_obj']), 'object RA (deg)') - ofd_b.header['DEC_OBJ'] = (np.float64( - self.information['dec_obj']), 'object Dec (deg)') - - else: - ofd_b.header['OBJECT'] = (self.source, ' object name') - ofd_b.header['TARGET'] = (0, 'target name (hhmmss.s+ddmmss)') - - ofd_b.header['OBSID'] = (str(obsid), 'observation ID') - - ofd_b.header['RA_OBJ'] = (np.float64(0), 'object RA (deg)') - ofd_b.header['DEC_OBJ'] = (np.float64(0), 'object Dec (deg)') - - # , 'the orbital position of CSST in X direction at EXPSTART') - self.information['POSI0_X'] = 0.0 - # , 'the orbital position of CSST in Y direction at EXPSTART') - self.information['POSI0_Y'] = 0.0 - # , 'the orbital position of CSST in Z direction at EXPSTART') - self.information['POSI0_Z'] = 0.0 - # , 'the orbital velocity of CSST in X direction at EXPSTART') - self.information['VELO0_X'] = 0.0 - # , 'the orbital velocity of CSST in Y direction at EXPSTART') - self.information['VELO0_Y'] = 0.0 - # , 'the orbital velocity of CSST in Z direction at EXPSTART') - self.information['VELO0_Z'] = 0.0 - - # , 'the orbital position of CSST in X direction at EXPSTART') - self.information['POSI1_X'] = 0.0 - # , 'the orbital position of CSST in Y direction at EXPSTART') - self.information['POSI1_Y'] = 0.0 - # , 'the orbital position of CSST in Z direction at EXPSTART') - self.information['POSI1_Z'] = 0.0 - # , 'the orbital velocity of CSST in X direction at EXPSTART') - self.information['VELO1_X'] = 0.0 - # , 'the orbital velocity of CSST in Y direction at EXPSTART') - self.information['VELO1_Y'] = 0.0 - # , 'the orbital velocity of CSST in Z direction at EXPSTART') - self.information['VELO1_Z'] = 0.0 - - ####################################################################### - - ######## Telescope information ############### - - ofd_b.header['REFFRAME'] = ('CSSTGSC-1.0', 'guiding catalog version') - ofd_b.header['DATE-OBS'] = (data_time[:21], - 'observation date (yyyy-mm-ddThh:mm:ss.s)') - ofd_b.header['SATESWV'] = ( - 'softwave-1.0', 'satellite software version') - ofd_b.header['EXPSTART'] = (np.float64( - time2mjd(exp_start_utc)), 'exposure start time (MJD)') - ofd_b.header['CABSTART'] = (np.float64( - time2jd(exp_start_utc)), 'first cabin time after exposure start (MJD)') - ofd_b.header['SUNANGL0'] = (np.float32( - 0.0), 'angle between the Sun and opt axis at CABSTART') - ofd_b.header['MOONANG0'] = (np.float32( - 0.0), 'angle between the Moon and opt axis at CABSTART') - ofd_b.header['TEL_ALT0'] = (np.float64( - 0.0), 'angle between opt axis and Elimb at CABSTART') - ofd_b.header['POS_ANG0'] = (np.float64( - 0.0), 'angle between y axis and North Pole at CABSTART') - ofd_b.header['POSI0_X'] = (np.float64( - self.information['POSI0_X']), 'orbital position in X at CABSTART (km)') - ofd_b.header['POSI0_Y'] = (np.float64( - self.information['POSI0_Y']), 'orbital position in Y at CABSTART (km)') - ofd_b.header['POSI0_Z'] = (np.float64( - self.information['POSI0_Z']), 'orbital position in Z at CABSTART (km)') - ofd_b.header['VELO0_X'] = (np.float64( - self.information['VELO0_X']), 'orbital velocity in X at CABSTART (km/s)') - ofd_b.header['VELO0_Y'] = (np.float64( - self.information['VELO0_Y']), 'orbital velocity in Y at CABSTART (km/s)') - ofd_b.header['VELO0_Z'] = (np.float64( - self.information['VELO0_Z']), 'orbital velocity in Z at CABSTART (km/s)') - ofd_b.header['Euler0_1'] = (np.float64( - 0.0), 'Euler angle 1 at CABSTART (deg)') - ofd_b.header['Euler0_2'] = (np.float64( - 0.0), 'Euler angle 2 at CABSTART (deg)') - ofd_b.header['Euler0_3'] = (np.float64( - 0.0), 'Euler angle 3 at CABSTART (deg)') - - ##### - if self.source == 'SCI' or self.source == 'COMP': - - ofd_b.header['RA_PNT0'] = (np.float64( - self.information['ra_pnt0']), 'pointing RA at CABSTART (deg)') - ofd_b.header['DEC_PNT0'] = (np.float64( - self.information['dec_pnt0']), 'pointing Dec at CABSTART (deg)') - - else: - - ofd_b.header['RA_PNT0'] = (np.float64( - 0.0), 'pointing RA at CABSTART (deg)') - ofd_b.header['DEC_PNT0'] = (np.float64( - 0.0), 'pointing Dec at CABSTART (deg)') - - ofd_b.header['EXPEND'] = (np.float64( - time2mjd(exp_end_utc)), 'exposure end time (MJD)') - - ofd_b.header['CABEND'] = (np.float64( - time2jd(exp_end_utc)), 'first cabin time after exposure end (MJD)') - - ofd_b.header['SUNANGL1'] = (np.float32( - 0.0), 'angle between the Sun and opt axis at CABEND') - ofd_b.header['MOONANG1'] = (np.float32( - 0.0), 'angle between the Moon and opt axis at CABEND') - - ofd_b.header['TEL_ALT1'] = (np.float64( - 0.0), 'angle between opt axis and Elimb at CABEND') - ofd_b.header['POS_ANG1'] = (np.float64( - self.information['sky_rot']), 'angle between y axis and North Pole at CABEND') - - ofd_b.header['POSI1_X'] = (np.float64( - self.information['POSI1_X']), 'orbital position in X at CABEND (km)') - ofd_b.header['POSI1_Y'] = (np.float64( - self.information['POSI1_Y']), 'orbital position in Y at CABEND (km)') - ofd_b.header['POSI1_Z'] = (np.float64( - self.information['POSI1_Z']), 'orbital position in Z at CABEND (km)') - ofd_b.header['VELO1_X'] = (np.float64( - self.information['VELO1_X']), 'orbital velocity in X at CABEND (km/s)') - ofd_b.header['VELO1_Y'] = (np.float64( - self.information['VELO1_Y']), 'orbital velocity in Y at CABEND (km/s)') - ofd_b.header['VELO1_Z'] = (np.float64( - self.information['VELO1_Z']), 'orbital velocity in Z at CABEND (km/s)') - ofd_b.header['Euler1_1'] = (np.float64( - 0.0), 'Euler angle 1 at CABEND (deg)') - ofd_b.header['Euler1_2'] = (np.float64( - 0.0), 'Euler angle 2 at CABEND (deg)') - ofd_b.header['Euler1_3'] = (np.float64( - 0.0), 'Euler angle 3 at CABEND (deg)') - ofd_b.header['RA_PNT1'] = (np.float64( - ofd_b.header['RA_PNT0']), 'pointing RA at CABEND (deg)') - ofd_b.header['DEC_PNT1'] = (np.float64( - ofd_b.header['DEC_PNT0']), 'pointing Dec at CABEND (deg)') - ofd_b.header['EPOCH'] = (np.float32( - time2jd(end.utcnow())), 'equinox of pointing RA and Dec') - ofd_b.header['EXPTIME'] = (np.float32( - self.information['exptime']*self.information['exposuretimes']), 'exposure time (s)') - ss1 = 'HDU checksum' - ss2 = 'data unit checksum' - ofd_b.header['CHECKSUM'] = (data_time[:19], ss1) - ofd_b.header['DATASUM'] = (data_time[:19], ss2) - - ########## - ########## finish header for 0 layer ###################### - ############################################################# - # header - # new image HDU, blue channel, layer 1 - if HeaderTest == 'yes': - hdu_b = fits.ImageHDU(data=np.zeros((10, 10, 3))) - else: - hdu_b = fits.ImageHDU(data=np.uint16(self.image_b_N)) - - ##### - - hdu_b.header['XTENSION'] = ('IMAGE', 'extension type') - hdu_b.header['BITPIX'] = (np.int16(16), 'bits per data value') - hdu_b.header['NAXIS'] = (np.int16(2), 'number of data axies') - - hdu_b.header['NAXIS1'] = (np.int32(1344), 'length of first data axis') - - hdu_b.header['NAXIS2'] = (np.int32(9672), 'length of second data axis') - - if self.information['exposuretimes'] > 1: - hdu_b.header['NAXIS3'] = ( - np.int32(self.information['exposuretimes']), 'length of third data axis') - - if self.source == 'SCI' or self.source == 'COMP': - hdu_b.header['EXTNAME'] = ('SCI', '') - - else: - - hdu_b.header['EXTNAME'] = ('CAL', '') - - hdu_b.header['EXTVER'] = (np.int16(1), '') - hdu_b.header['BSCALE'] = (np.float64(1), '') - hdu_b.header['BZERO'] = (np.float64(32768), '') - hdu_b.header['BUNIT'] = ('ADU', 'physical unit of array values') - - ####################################################################### - ######### instrument information ################################### - if self.source == 'SCI' or self.source == 'COMP': - - hdu_b.header['CMIRRPOS'] = ( - bool(False), 'position of calibration switch mirror') - else: - - hdu_b.header['CMIRRPOS'] = ( - bool(True), 'position of calibration switch mirror') - - if self.source == 'FLAT': - hdu_b.header['FLAMP'] = (np.int16(1), 'status of flat lamp') - else: - hdu_b.header['FLAMP'] = (np.int16(0), 'status of flat lamp') - - if self.source == 'LAMP': - hdu_b.header['ALAMP'] = ( - np.int16(1), 'status of flat atomic emission line lamp') - else: - hdu_b.header['ALAMP'] = ( - np.int16(0), 'status of flat atomic emission line lamp') - - ############# - hdu_b.header['IFSMODE'] = (np.int32(0), 'IFS working mode') - hdu_b.header['IFSTEMP'] = (np.float32( - 0), 'IFS components temperature in degC') - hdu_b.header['IFSSTAT'] = ( - np.int32(0), 'IFS components status parameter') - - ###################################################################### - - #################### detector information########################## - - hdu_b.header['CAMERA'] = ('Blue', 'camera of IFS') - hdu_b.header['DETSN'] = ('CCD231-c4-00', 'detector serial number') - hdu_b.header['DETNAME'] = ('blue', 'detector name') - hdu_b.header['DETSIZE'] = ('2048*4096', 'detector size') - hdu_b.header['DATASECT'] = ('data-section', 'data section') - hdu_b.header['PIXSCAL1'] = (np.float32(1344), 'pixel scale for axis 1') - hdu_b.header['PIXSCAL2'] = (np.float32(9672), 'pixel scale for axis 2') - hdu_b.header['PIXSIZE1'] = (np.float32( - 15), 'pixel size for axis 1 (micron)') - hdu_b.header['PIXSIZE2'] = (np.float32( - 15), 'pixel size for axis 2 (micron)') - - hdu_b.header['NCHAN'] = (np.int16(4), 'number of readout channels') - hdu_b.header['PSCAN1'] = ( - np.int32(50), 'horizontal prescan width, per readout channel') - hdu_b.header['PSCAN2'] = ( - np.int32(0), 'vertical prescan height, per readout channel') - hdu_b.header['OSCAN1'] = ( - np.int32(320), 'horizontal overscan width, per readout channel') - hdu_b.header['OSCAN2'] = ( - np.int32(320), 'vertical overscan height, per readout channel') - - ####################################################################### - # Readout information, revised on 2024.2.27 - - frame_time_b = 0.09216 # data frame transfer time in Bule camera - exptimes = self.information['exposuretimes'] - for k in range(exptimes): - - tt = (self.information['exptime']+frame_time_b)*(k) - temp = self.dt+timedelta(seconds=tt) - read_start = temp.strftime("%Y-%m-%dT%H:%M:%S.%f") - - temp2 = temp+timedelta(seconds=self.information['exptime']) - read_end = temp2.strftime("%Y-%m-%dT%H:%M:%S.%f") - - str1 = 'EXPT0_'+str(k+1) - str2 = 'EXPT1_'+str(k+1) - - if k == 0: - com1 = 'exposure start time of the 1st frame (UTC)' - com2 = 'exposure end time of the 1st frame (UTC)' - if k == 1: - com1 = 'exposure start time of the 2nd frame (UTC)' - com2 = 'exposure end time of the 2nd frame (UTC)' - if k == 2: - com1 = 'exposure start time of the 3rd frame (UTC)' - com2 = 'exposure end time of the 3rd frame (UTC)' - if k > 2: - com1 = 'exposure start time of the '+str(k+1)+'th frame' - com2 = 'exposure end time of the '+str(k+1)+'th frame' - - hdu_b.header[str1] = (read_start[:21], com1) - hdu_b.header[str2] = (read_end[:21], com2) - - hdu_b.header['DETTEMP0'] = (np.float32( - 0.0), 'detector temperature at EXPT0_1 (K)') - hdu_b.header['DETTEMP1'] = (np.float32( - 0.0), 'detector temperature at EXPT1_'+str(k+1)+' (K)') - - #######################end revised on 2024.2.27 ### - - hdu_b.header['BIN_X'] = ( - np.int16(1), 'bin number in X (wavelength)') - hdu_b.header['BIN_Y'] = (np.int16(1), 'bin number in Y (spatial)') - - hdu_b.header['GAIN1'] = (np.float32( - self.information['gain1_b']), 'CCD gain (e-/ADU, channel 1)') - hdu_b.header['GAIN2'] = (np.float32( - self.information['gain2_b']), 'CCD gain (e-/ADU, channel 2)') - hdu_b.header['GAIN3'] = (np.float32( - self.information['gain3_b']), 'CCD gain (e-/ADU, channel 3)') - hdu_b.header['GAIN4'] = (np.float32( - self.information['gain4_b']), 'CCD gain (e-/ADU, channel 4)') - - hdu_b.header['DARK1'] = (np.float32( - self.information['dark1_b']), 'CCD dark (e-/s, channel 1)') - hdu_b.header['DARK2'] = (np.float32( - self.information['dark2_b']), 'CCD dark (e-/s, channel 2)') - hdu_b.header['DARK3'] = (np.float32( - self.information['dark3_b']), 'CCD dark (e-/s, channel 3)') - hdu_b.header['DARK4'] = (np.float32( - self.information['dark4_b']), 'CCD dark (e-/s, channel 4)') - - hdu_b.header['RON1'] = (np.float32( - self.information['rn1_b']), 'read noise (e-, channel 1)') - hdu_b.header['RON2'] = (np.float32( - self.information['rn2_b']), 'read noise (e-, channel 2)') - hdu_b.header['RON3'] = (np.float32( - self.information['rn3_b']), 'read noise (e-, channel 3)') - hdu_b.header['RON4'] = (np.float32( - self.information['rn4_b']), 'read noise (e-, channel 4)') - - hdu_b.header['DETBIA1'] = (np.float32( - self.information['bias1_b']), 'amplifier bias voltage (ADU, channel1)') - hdu_b.header['DETBIA2'] = (np.float32( - self.information['bias2_b']), 'amplifier bias voltage (ADU, channel2)') - hdu_b.header['DETBIA3'] = (np.float32( - self.information['bias3_b']), 'amplifier bias voltage (ADU, channel3)') - hdu_b.header['DETBIA4'] = (np.float32( - self.information['bias4_b']), 'amplifier bias voltage (ADU, channel4)') - - hdu_b.header['ROSPEED'] = (np.float32(100), 'readout speed (MHz)') - - hdu_b.header['TOML_TAG'] = ('csst-ifs-l0-c09', 'the tag of toml file') - - tt = write_end.utcnow() - date_time = tt.strftime("%Y-%m-%dT%H:%M:%S") - - hdu_b.header['CHECKSUM'] = (date_time, 'HDU checksum') - - hdu_b.header['DATASUM'] = (data_time, 'data unit checksum') - - if self.source == 'LAMP' and self.information['holemask'] == 'yes': - hdu_b.header['Hole'] = ('yes', 'apply hole to LAMP') - - ##################################################################### - #################### red camera ###################### - - # create a new FITS file, using HDUList instance - ofd_r = fits.PrimaryHDU() - - ofd_r.header['SIMPLE'] = (True, 'conforms to FITS standard') - ofd_r.header['BITPIX'] = (0, 'array data type') - ofd_r.header['NAXIS'] = (0, 'number of array dimensions') - ofd_r.header['NEXTEND'] = (np.int16(1), 'number of array dimensions') - - temp = write_end.utcnow() - data_time = temp.strftime("%Y-%m-%dT%H:%M:%S.%f") - - ofd_r.header['DATE'] = ( - data_time[:21], 'written date (yyyy-mm-ddThh:mm:ss.s)') - ofd_r.header['FILENAME'] = (filename_r[:68], '') - - ofd_r.header['FILETYPE'] = (self.source[:12], 'observation type') - - ofd_r.header['TELESCOP'] = ('CSST', 'telescope name') - - ofd_r.header['INSTRUME'] = ('IFS', 'instrument name') - - ofd_r.header['RADECSYS'] = ( - 'ICRS', 'coordinate system of the object') - - ofd_r.header['EQUINOX'] = (float(2000.0), '') - ofd_r.header['FITSSWV'] = ( - 'ifs_sim_0.8.03', 'FITS creating software version') - - ######### Object information ############# - if self.source == 'SCI' or self.source == 'COMP': - ofd_r.header['OBJECT'] = ( - self.information['name_obj'][:30], 'object name') - ofd_r.header['TARGET'] = ( - (self.information['target']), 'target name (hhmmss.s+ddmmss)') - - ofd_r.header['OBSID'] = (str(obsid), 'observation ID') - - ofd_r.header['RA_OBJ'] = (np.float64( - self.information['ra_obj']), 'object RA (deg)') - ofd_r.header['DEC_OBJ'] = (np.float64( - self.information['dec_obj']), 'object Dec (deg)') - - else: - ofd_r.header['OBJECT'] = (self.source, ' object name') - ofd_r.header['TARGET'] = (0, 'target name (hhmmss.s+ddmmss)') - - ofd_r.header['OBSID'] = (str(obsid), 'observation ID') - - ofd_r.header['RA_OBJ'] = (np.float64(0), 'object RA (deg)') - ofd_r.header['DEC_OBJ'] = (np.float64(0), 'object Dec (deg)') - - # , 'the orbital position of CSST in X direction at EXPSTART') - self.information['POSI0_X'] = 0.0 - # , 'the orbital position of CSST in Y direction at EXPSTART') - self.information['POSI0_Y'] = 0.0 - # , 'the orbital position of CSST in Z direction at EXPSTART') - self.information['POSI0_Z'] = 0.0 - # , 'the orbital velocity of CSST in X direction at EXPSTART') - self.information['VELO0_X'] = 0.0 - # , 'the orbital velocity of CSST in Y direction at EXPSTART') - self.information['VELO0_Y'] = 0.0 - # , 'the orbital velocity of CSST in Z direction at EXPSTART') - self.information['VELO0_Z'] = 0.0 - - # , 'the orbital position of CSST in X direction at EXPSTART') - self.information['POSI1_X'] = 0.0 - # , 'the orbital position of CSST in Y direction at EXPSTART') - self.information['POSI1_Y'] = 0.0 - # , 'the orbital position of CSST in Z direction at EXPSTART') - self.information['POSI1_Z'] = 0.0 - # , 'the orbital velocity of CSST in X direction at EXPSTART') - self.information['VELO1_X'] = 0.0 - # , 'the orbital velocity of CSST in Y direction at EXPSTART') - self.information['VELO1_Y'] = 0.0 - # , 'the orbital velocity of CSST in Z direction at EXPSTART') - self.information['VELO1_Z'] = 0.0 - - ####################################################################### - - ######## Telescope information ############### - - ofd_r.header['REFFRAME'] = ('CSSTGSC-1.0', 'guiding catalog version') - ofd_r.header['DATE-OBS'] = (data_time[:21], - 'observation date (yyyy-mm-ddThh:mm:ss.s) ') - - ofd_r.header['SATESWV'] = ( - 'softwave-1.0', 'satellite software version') - - ofd_r.header['EXPSTART'] = (np.float64( - time2mjd(exp_start_utc)), 'exposure start time (MJD)') - - ofd_r.header['CABSTART'] = (np.float64( - time2jd(exp_start_utc)), 'first cabin time after exposure start (MJD)') - - ofd_r.header['SUNANGL0'] = (np.float32( - 0.0), 'angle between the Sun and opt axis at CABSTART') - ofd_r.header['MOONANG0'] = (np.float32( - 0.0), 'angle between the Moon and opt axis at CABSTART') - - ofd_r.header['TEL_ALT0'] = (np.float64( - 0.0), 'angle between opt axis and Elimb at CABSTART') - ofd_r.header['POS_ANG0'] = (np.float64( - 0.0), 'angle between y axis and North Pole at CABSTART') - - ofd_r.header['POSI0_X'] = (np.float64( - self.information['POSI0_X']), 'orbital position in X at CABSTART (km)') - ofd_r.header['POSI0_Y'] = (np.float64( - self.information['POSI0_Y']), 'orbital position in Y at CABSTART (km)') - ofd_r.header['POSI0_Z'] = (np.float64( - self.information['POSI0_Z']), 'orbital position in Z at CABSTART (km)') - ofd_r.header['VELO0_X'] = (np.float64( - self.information['VELO0_X']), 'orbital velocity in X at CABSTART (km/s)') - ofd_r.header['VELO0_Y'] = (np.float64( - self.information['VELO0_Y']), 'orbital velocity in Y at CABSTART (km/s)') - ofd_r.header['VELO0_Z'] = (np.float64( - self.information['VELO0_Z']), 'orbital velocity in Z at CABSTART (km/s)') - - ofd_r.header['Euler0_1'] = (np.float64( - 0.0), 'Euler angle 1 at CABSTART (deg)') - ofd_r.header['Euler0_2'] = (np.float64( - 0.0), 'Euler angle 2 at CABSTART (deg)') - ofd_r.header['Euler0_3'] = (np.float64( - 0.0), 'Euler angle 3 at CABSTART (deg)') - - ##### - if self.source == 'SCI' or self.source == 'COMP': - - ofd_r.header['RA_PNT0'] = (np.float64( - self.information['ra_pnt0']), 'pointing RA at CABSTART (deg)') - ofd_r.header['DEC_PNT0'] = (np.float64( - self.information['dec_pnt0']), 'pointing Dec at CABSTART (deg)') - - else: - - ofd_r.header['RA_PNT0'] = (np.float64( - 0.0), 'pointing RA at CABSTART (deg)') - ofd_r.header['DEC_PNT0'] = (np.float64( - 0.0), 'pointing Dec at CABSTART (deg)') - - ofd_r.header['EXPEND'] = (np.float64( - time2mjd(exp_end_utc)), 'exposure end time (MJD)') - - ofd_r.header['CABEND'] = (np.float64( - time2jd(exp_end_utc)), 'first cabin time after exposure end (MJD)') - - ofd_r.header['SUNANGL1'] = (np.float32( - 0.0), 'angle between the Sun and opt axis at CABEND') - ofd_r.header['MOONANG1'] = (np.float32( - 0.0), 'angle between the Moon and opt axis at CABEND') - - ofd_r.header['TEL_ALT1'] = (np.float64( - 0.0), 'angle between opt axis and Elimb at CABEND') - ofd_r.header['POS_ANG1'] = (np.float64( - self.information['sky_rot']), 'angle between y axis and North Pole at CABEND') - - ofd_r.header['POSI1_X'] = (np.float64( - self.information['POSI1_X']), 'orbital position in X at CABEND (km)') - ofd_r.header['POSI1_Y'] = (np.float64( - self.information['POSI1_Y']), 'orbital position in Y at CABEND (km)') - ofd_r.header['POSI1_Z'] = (np.float64( - self.information['POSI1_Z']), 'orbital position in Z at CABEND (km)') - - ofd_r.header['VELO1_X'] = (np.float64( - self.information['VELO1_X']), 'orbital velocity in X at CABEND (km/s)') - ofd_r.header['VELO1_Y'] = (np.float64( - self.information['VELO1_Y']), 'orbital velocity in Y at CABEND (km/s)') - ofd_r.header['VELO1_Z'] = (np.float64( - self.information['VELO1_Z']), 'orbital velocity in Z at CABEND (km/s)') - - ofd_r.header['Euler1_1'] = (np.float64( - 0.0), 'Euler angle 1 at CABEND (deg)') - ofd_r.header['Euler1_2'] = (np.float64( - 0.0), 'Euler angle 2 at CABEND (deg)') - ofd_r.header['Euler1_3'] = (np.float64( - 0.0), 'Euler angle 3 at CABEND (deg)') - - ofd_r.header['RA_PNT1'] = (np.float64( - ofd_r.header['RA_PNT0']), 'pointing RA at CABEND (deg)') - ofd_r.header['DEC_PNT1'] = (np.float64( - ofd_r.header['DEC_PNT0']), 'pointing Dec at CABEND (deg)') - - ofd_r.header['EPOCH'] = (np.float32( - time2jd(end.utcnow())), 'equinox of pointing RA and Dec') - ofd_r.header['EXPTIME'] = (np.float32( - self.information['exptime']*self.information['exposuretimes']), 'exposure time (s)') - - ss1 = 'HDU checksum' - ss2 = 'data unit checksum' - ofd_r.header['CHECKSUM'] = (data_time[:19], ss1) - ofd_r.header['DATASUM'] = (data_time[:19], ss2) - - ########## - ########## finish header for 0 layer ###################### - ############################################################# - - # header - # new image HDU, red channel, layer 1 - if HeaderTest == 'yes': - hdu_r = fits.ImageHDU(data=np.zeros((10, 10, 3))) - else: - hdu_r = fits.ImageHDU(data=np.uint16(self.image_r_N)) - # - hdu_r.header['XTENSION'] = ('IMAGE', 'extension type') - hdu_r.header['BITPIX'] = (np.int16(16), 'bits per data value') - hdu_r.header['NAXIS'] = (np.int16(2), 'number of data axies') - - hdu_r.header['NAXIS1'] = (np.int32(1856), 'length of first data axis') - - hdu_r.header['NAXIS2'] = ( - np.int32(13768), 'length of second data axis') - - if self.information['exposuretimes'] > 1: - hdu_r.header['NAXIS3'] = ( - np.int32(self.information['exposuretimes']), 'length of third data axis') - - if self.source == 'SCI' or self.source == 'COMP': - hdu_r.header['EXTNAME'] = ('SCI', '') - - else: - - hdu_r.header['EXTNAME'] = ('CAL', '') - - hdu_r.header['EXTVER'] = (np.int16(1), '') - hdu_r.header['BSCALE'] = (np.float64(1), '') - hdu_r.header['BZERO'] = (np.float64(32768), '') - hdu_r.header['BUNIT'] = ('ADU', 'physical unit of array values') - - ######### instrument information ###### - - if self.source == 'SCI' or self.source == 'COMP': - - hdu_r.header['CMIRRPOS'] = ( - bool(False), 'position of calibration switch mirror') - else: - - hdu_r.header['CMIRRPOS'] = ( - bool(True), 'position of calibration switch mirror') - - if self.source == 'FLAT': - hdu_r.header['FLAMP'] = (np.int16(1), 'status of flat lamp') - else: - hdu_r.header['FLAMP'] = (np.int16(0), 'status of flat lamp') - - if self.source == 'LAMP': - hdu_r.header['ALAMP'] = ( - np.int16(1), 'status of flat atomic emission line lamp') - else: - hdu_r.header['ALAMP'] = ( - np.int16(0), 'status of flat atomic emission line lamp') - - ############# - hdu_r.header['IFSMODE'] = (np.int32(0), 'IFS working mode') - hdu_r.header['IFSTEMP'] = (np.float32( - 0), 'IFS components temperature in degC') - hdu_r.header['IFSSTAT'] = ( - np.int32(0), 'IFS components status parameter') - - ################### detector information############################ - - hdu_r.header['CAMERA'] = ('Red', 'camera of IFS') - hdu_r.header['DETSN'] = ('CCD231-c6-00', 'detector serial number') - hdu_r.header['DETNAME'] = ('red', 'detector name') - - hdu_r.header['DETSIZE'] = ('3072*6144', 'detector size') - hdu_r.header['DATASECT'] = ('data-section', 'data section') - hdu_r.header['PIXSCAL1'] = (np.float32(1856), 'pixel scale for axis 1') - hdu_r.header['PIXSCAL2'] = ( - np.float32(13768), 'pixel scale for axis 2') - hdu_r.header['PIXSIZE1'] = (np.float32( - 15), 'pixel size for axis 1 (micron)') - hdu_r.header['PIXSIZE2'] = (np.float32( - 15), 'pixel size for axis 2 (micron)') - - hdu_r.header['NCHAN'] = (np.int16(4), 'number of readout channels') - - hdu_r.header['PSCAN1'] = ( - np.int32(50), 'horizontal prescan width, per readout channel') - hdu_r.header['PSCAN2'] = ( - np.int32(0), 'vertical prescan height, per readout channel') - hdu_r.header['OSCAN1'] = ( - np.int32(320), 'horizontal overscan width, per readout channel') - hdu_r.header['OSCAN2'] = ( - np.int32(320), 'vertical overscan height, per readout channel') - -############################################################################## - # Readout information - frame_time_r = 0.13824 # data frame transfer time in red camera - - for k in range(exptimes): - tt = (self.information['exptime']+frame_time_r)*(k) - - temp = self.dt+timedelta(seconds=tt) - read_start = temp.strftime("%Y-%m-%dT%H:%M:%S.%f") - - temp2 = temp+timedelta(seconds=self.information['exptime']) - read_end = temp2.strftime("%Y-%m-%dT%H:%M:%S.%f") - - str1 = 'EXPT0_'+str(k+1) - str2 = 'EXPT1_'+str(k+1) - - if k == 0: - com1 = 'exposure start time of the 1st frame (UTC)' - com2 = 'exposure end time of the 1st frame (UTC)' - if k == 1: - com1 = 'exposure start time of the 2nd frame (UTC)' - com2 = 'exposure end time of the 2nd frame (UTC)' - if k == 2: - com1 = 'exposure start time of the 3rd frame (UTC)' - com2 = 'exposure end time of the 3rd frame (UTC)' - if k > 2: - com1 = 'exposure start time of the '+str(k+1)+'th frame' - com2 = 'exposure end time of the '+str(k+1)+'th frame' - - hdu_r.header[str1] = (read_start[:21], com1) - hdu_r.header[str2] = (read_end[:21], com2) - #################################################################################### - hdu_r.header['DETTEMP0'] = (np.float32( - 0.0), 'detector temperature at EXPT0_1 (K)') - hdu_r.header['DETTEMP1'] = (np.float32( - 0.0), 'detector temperature at EXPT1_'+str(k+1)+' (K)') - - hdu_r.header['BIN_X'] = ( - np.int16(1), 'bin number in X (wavelength)') - hdu_r.header['BIN_Y'] = (np.int16(1), 'bin number in Y (spatial)') - - hdu_r.header['GAIN1'] = (np.float32( - self.information['gain1_r']), 'CCD gain (e-/ADU, channel 1)') - hdu_r.header['GAIN2'] = (np.float32( - self.information['gain2_r']), 'CCD gain (e-/ADU, channel 2)') - hdu_r.header['GAIN3'] = (np.float32( - self.information['gain3_r']), 'CCD gain (e-/ADU, channel 3)') - hdu_r.header['GAIN4'] = (np.float32( - self.information['gain4_r']), 'CCD gain (e-/ADU, channel 4)') - - hdu_r.header['DARK1'] = (np.float32( - self.information['dark1_r']), 'CCD dark (e-/s, channel 1)') - hdu_r.header['DARK2'] = (np.float32( - self.information['dark2_r']), 'CCD dark (e-/s, channel 2)') - hdu_r.header['DARK3'] = (np.float32( - self.information['dark3_r']), 'CCD dark (e-/s, channel 3)') - hdu_r.header['DARK4'] = (np.float32( - self.information['dark4_r']), 'CCD dark (e-/s, channel 4)') - - hdu_r.header['RON1'] = (np.float32( - self.information['rn1_r']), 'read noise (e-, channel 1)') - hdu_r.header['RON2'] = (np.float32( - self.information['rn2_r']), 'read noise (e-, channel 2)') - hdu_r.header['RON3'] = (np.float32( - self.information['rn3_r']), 'read noise (e-, channel 3)') - hdu_r.header['RON4'] = (np.float32( - self.information['rn4_r']), 'read noise (e-, channel 4)') - - hdu_r.header['DETBIA1'] = (np.float32( - self.information['bias1_r']), 'amplifier bias voltage (ADU, channel1)') - hdu_r.header['DETBIA2'] = (np.float32( - self.information['bias2_r']), 'amplifier bias voltage (ADU, channel2)') - hdu_r.header['DETBIA3'] = (np.float32( - self.information['bias3_r']), 'amplifier bias voltage (ADU, channel3)') - hdu_r.header['DETBIA4'] = (np.float32( - self.information['bias4_r']), 'amplifier bias voltage (ADU, channel4)') - - hdu_r.header['ROSPEED'] = (np.float32(100), 'readout speed (MHz)') - hdu_r.header['TOML_TAG'] = ('csst-ifs-l0-c09', 'the tag of toml file') - - date_timer = write_end.strftime("%Y-%m-%dT%H:%M:%S") - - hdu_r.header['CHECKSUM'] = (date_timer, 'HDU checksum') - - hdu_r.header['DATASUM'] = (date_timer, 'data unit checksum') - - if self.source == 'LAMP' and self.information['holemask'] == 'yes': - hdu_r.header['Hole'] = ('yes', 'apply hole to LAMP') - - ######################### - - hdu1 = fits.PrimaryHDU(header=ofd_b.header) - hdu2 = fits.PrimaryHDU(header=ofd_r.header) - ##################################################################### - hdu1.header.add_comment( - '========================================================================', after='FITSSWV') - hdu1.header.add_comment('OBJECT INFORMATION', after='FITSSWV') - hdu1.header.add_comment( - '========================================================================', after='FITSSWV') - - hdu1.header.add_comment( - '========================================================================', after='DEC_OBJ') - hdu1.header.add_comment('TELESCOPE INFORMATION', after='DEC_OBJ') - hdu1.header.add_comment( - '========================================================================', after='DEC_OBJ') - - hdu1.header.add_comment( - '========================================================================', after='EPOCH') - hdu1.header.add_comment('VERIFICATION INFORMATION', after='EPOCH') - hdu1.header.add_comment( - '========================================================================', after='EPOCH') - - ###################################################################### - hdu2.header.add_comment( - '========================================================================', after='FITSSWV') - - hdu2.header.add_comment( - '========================================================================', after='DEC_OBJ') - hdu2.header.add_comment('TELESCOPE INFORMATION', after='DEC_OBJ') - hdu2.header.add_comment( - '========================================================================', after='DEC_OBJ') - - hdu2.header.add_comment( - '========================================================================', after='EPOCH') - hdu2.header.add_comment('VERIFICATION INFORMATION', after='EPOCH') - hdu2.header.add_comment( - '========================================================================', after='EPOCH') - - hdulist_b = fits.HDUList([hdu1, hdu_b]) - - hdu_b.header.add_comment( - '========================================================================', before='CMIRRPOS') - hdu_b.header.add_comment('INSTRUMENT INFORMATION', before='CMIRRPOS') - hdu_b.header.add_comment( - '========================================================================', before='CMIRRPOS') - - hdu_b.header.add_comment( - '========================================================================', before='CAMERA') - hdu_b.header.add_comment('DETECTOR INFORMATION', before='CAMERA') - hdu_b.header.add_comment( - '========================================================================', before='CAMERA') - - hdu_b.header.add_comment( - '========================================================================', before='EXPT0_1') - hdu_b.header.add_comment('READOUT INFORMATION', before='EXPT0_1') - hdu_b.header.add_comment( - '========================================================================', before='EXPT0_1') - - hdulist_b.writeto(self.file_b, output_verify='ignore', checksum=True) - ######################################################################### - hdulist_r = fits.HDUList([hdu2, hdu_r]) - - hdu_r.header.add_comment( - '========================================================================', before='CMIRRPOS') - hdu_r.header.add_comment('INSTRUMENT INFORMATION', before='CMIRRPOS') - hdu_r.header.add_comment( - '========================================================================', before='CMIRRPOS') - - hdu_r.header.add_comment( - '========================================================================', before='CAMERA') - hdu_r.header.add_comment('DETECTOR INFORMATION', before='CAMERA') - hdu_r.header.add_comment( - '========================================================================', before='CAMERA') - - hdu_r.header.add_comment( - '========================================================================', before='EXPT0_1') - hdu_r.header.add_comment('READOUT INFORMATION', before='EXPT0_1') - hdu_r.header.add_comment( - '========================================================================', before='EXPT0_1') - - hdulist_r.writeto(self.file_r, output_verify='ignore', checksum=True) - - ###################################################################### - - def earthshine(self, theta): - """ - For given theta angle, return the earth-shine spectrum. - - :param theta: angle (in degree) from the target to earth limb. - :return: the scaled solar spectrum - template_wave: unit in A - template_flux: unit in erg/s/cm^2/A/arcsec^2 - - """ - - # read solar template - solar_template = pd.read_csv(self.information['dir_path']+'IFS_inputdata/refs/solar_spec.dat', sep='\s+', - header=None, comment='#') - template_wave = solar_template[0].values - template_flux = solar_template[1].values - - # read earth shine surface brightness - earthshine_curve = pd.read_csv(self.information['dir_path']+'IFS_inputdata/refs/earthshine.dat', - header=None, comment='#') - angle = earthshine_curve[0].values - surface_brightness = earthshine_curve[1].values - - # read V-band throughtput - cat_filter_V = pd.read_csv(self.information['dir_path']+'IFS_inputdata/refs/filter_Bessell_V.dat', sep='\s+', - header=None, comment='#') - filter_wave = cat_filter_V[0].values - filter_response = cat_filter_V[1].values - - # interplate to the target wavelength in V-band - ind_filter = (template_wave >= np.min(filter_wave)) & ( - template_wave <= np.max(filter_wave)) - filter_wave_interp = template_wave[ind_filter] - filter_response_interp = np.interp( - filter_wave_interp, filter_wave, filter_response) - - filter_constant = simps( - filter_response_interp * filter_wave_interp, filter_wave_interp) - template_constant = simps(filter_response_interp * template_wave[ind_filter] * template_flux[ind_filter], - template_wave[ind_filter]) - dwave = filter_wave_interp[1:] - filter_wave_interp[:-1] - wave_eff = np.nansum(dwave * filter_wave_interp[1:] * filter_response_interp[1:]) / \ - np.nansum(dwave * filter_response_interp[1:]) - - # get the normalized value at theta. - u0 = np.interp(theta, angle, surface_brightness) # mag/arcsec^2 - # target flux in erg/s/cm^2/Hz unit - u0 = 10**((u0 + 48.6)/(-2.5)) - u0 = u0 * 3e18 / wave_eff**2 # erg/s/cm^2/A/arcsec^2 - - factor = u0 * filter_constant / template_constant - norm_flux = template_flux * factor # erg/s/cm^2/A/arcsec^2 - - self.earthshine_wave = template_wave # A - self.earthshine_flux = norm_flux - - return - - -######################################################################################################################################################################################################################################################## - - ################################################################################### - ################################################################################## - - def CalskyNoise(self, lam): - """ - - - Parameters - ---------- - lam : TYPE - DESCRIPTION. - - Returns - ------- - Ns_skynoise : TYPE - DESCRIPTION. - - """ - - # calculate the reference flux - # calculate sky noise; - planckh = 6.62620*10**-27 # % erg s; - cc = 2.99792458*10**17 # in nm/s - fov2 = 0.01 # 单个像素的面积 arcsec^2 - # telarea=3.1415926*100*100 # in cm^2，望远镜聚光面积； - # lam is input wavelength in nm - - ########################################## - # self.earthshine_wave # A - # self.earthshine_flux # erg/s/cm^2/A/arcsec^2 - earthshine_flux = np.interp( - lam*10.0, self.earthshine_wave, self.earthshine_flux) # flux from zodiacal - - # self.zodiacal_wave # in A - # self.zodiacal_flux # erg/s/cm^2/A/arcsec^2 - zodiacal_flux = np.interp( - lam*10.0, self.zodiacal_wave, self.zodiacal_flux) # flux from zodiacal - fluxlam_sky = (earthshine_flux+zodiacal_flux) * \ - fov2 # erg/s/cm2/A - ############### - - ephoton = planckh*cc/lam # in erg/photon, cc与lambda单位需要一致； - Ns_skynoise = fluxlam_sky/ephoton # in unit of photons/cm2/s/A - - return Ns_skynoise - - ################################################################################# - ############## - - def sim_lamp_hole_img(self, exposuretime): - """ - - - Parameters - ---------- - exposuretime : TYPE - DESCRIPTION. - - Returns - ------- - None. - - """ - - sizeout = 1000 # here image has the pixelscale =0.01 arcsec - ################# create the slicer and slit file for high resolution 0.01arcsec simulation ################### - - maskSlice = dict() - maskSlit = dict() - - for k in range(32): - maskSlice[str(k)] = np.zeros((sizeout, sizeout)) - maskSlice[str(k)][20*k+180:20*k+200, int(sizeout/2) - - 320:int(sizeout/2)+320] = 1 - - maskSlit[str(k)] = np.zeros((sizeout, sizeout)) - maskSlit[str(k)][20*k+180:20*k+200, int(sizeout/2) - - 370:int(sizeout/2)+370] = 1 - - self.maskSliceHole = maskSlice - self.maskSlitHole = maskSlit - ######################################################################### - ##### load slicer_Qe1K.fits #### - slicerfile = self.information['dir_path'] + \ - 'IFS_inputdata/Flatfield/slicer_QE_1K_230625.fits' - da = fits.open(slicerfile) - self.log.info('hole Sim and slicer 1K QE file:%s' % (slicerfile)) - slicer_Qe = da[0].data - - ##### load hole mask ###### - # load hole mask matrix - file = self.information['dir_path'] + \ - 'IFS_inputdata/Hole/holemask_230612.fits' - self.log.info('hole mask file is %s' % (file)) - da = fits.open(file) - HoleMask = da[self.simnumber].data - - # load HgAr data; - - matfn3 = self.information['dir_path']+'IFS_inputdata/HgAr.mat' - self.log.info('lamp flux file is : %s' % (matfn3)) - da3 = sio.loadmat(matfn3) - HgArst = da3['HgArst'] # opd unit is in meter - # load flat data - - exptime = self.information['exptime'] # exposure time - - blue_img = galsim.Image(np.zeros((2048, 4096)), copy=True) - blue_img.scale = self.pixelscale - blue_img.setOrigin = (0, 0) - - red_img = galsim.Image(np.zeros((3072, 6144)), copy=True) - red_img.scale = self.pixelscale - red_img.setOrigin = (0, 0) - - slice_image = np.zeros((sizeout, sizeout)) - - width_blue = 0 - - LightSt = np.zeros(6000) - flux_1 = np.zeros(6000) - Width_lambda = (1000-350)/6000 - Wave = np.zeros(6000) - - if self.debug: - alfa = 1000 - else: - alfa = 1 - - ######################################### - - for klam in range(6000): - - ilam = klam*alfa - - if self.debug: - print(ilam) - - if ilam % 100 == 0: - self.log.info('ilam = %i' % ilam) - self.log.info('blue image max value = %i' % - blue_img.array.max()) - self.log.info('red image max value = %i' % red_img.array.max()) - - if ilam >= 6000: - break - - # the wavelength of the i-th frame data - lam = 350.0+(1000-350)/6000.0 * ilam - - Wave[ilam] = lam - - #print('i = ',ilam ) - ############################################### - - sumt = 0.0 - - for k in range(11): - t = getSpectrum(HgArst, lam+(k-5) * - Width_lambda/11, self.HgArsigma) - sumt = sumt+t[0]/11 - #print('t=', t) - - # the normalized HgAr light source in the range of 0 and 1.0 - LightSt[ilam] = sumt - - skymag = 6 # mag/arcsec^2 unit - - if ilam == 0: - # fluxnu1,fluxlam1,fluxphoton1=mag2flux(skymag,633) - fluxlam1, fluxphoton1 = mag2flux(skymag, 633) - - # save flux information - flux_1[ilam] = LightSt[ilam]*fluxlam1*self.fov2 # erg/s/cm2/A - - nphoton = flux2photons(flux_1[ilam], lam) - - ns = nphoton*np.ones((sizeout, sizeout)) # photons/cm2/s/A/arcsec2 - - slice_image = ns*self.telarea*Width_lambda*10*exptime - - ############################################################# - - if slice_image.max() < 0.1: - continue - - if slice_image.min() < 0: - self.log.error('Error in datacube, negative values !!!!!!!') - continue - - ######################## - ############## - # here Q is 10 times - Q = 10*lam*1e-3*self.Fnum/self.pixelsize - - # for calibration , there is no primay CSST system, but diffraction effect should consider - # opd0 rms = 0.0 @632.8nm - wavefront0 = self.opd0/lam/100000.0 - - psf0 = anySampledPSFnew(wavefront0, self.pupil, Q, 256) - - # do convolve with psf0 - conv = fftconvolve(slice_image, psf0, mode='same') - - conv[conv < 0.0] = 0.0 - - slice_image2 = conv - - # 3 - - # applyHoleMask only for lamp - - img0 = slice_image2*HoleMask - - # opd1 rms =0.075 @632.8nm - wavefront = self.opd1/lam - - psf1 = anySampledPSFnew(wavefront, self.pupil, Q, 256) - - # do convolve with psf1 - conv = fftconvolve(img0, psf1, mode='same') - conv[conv < 0.0] = 0.0 - img0 = conv - - ####################################################### - - - # CCD quantum efficiency - CCD_Qe_lam = np.interp(lam, self.CCD_Qe[:, 0], self.CCD_Qe[:, 1]) - - # optical efficiency , convert the wavelength to A - Qe_blue = np.interp( - lam, self.optical_blue_Q[:, 0], self.optical_blue_Q[:, 1]) - - # optical efficiency, convert the wavelength to A - Qe_red = np.interp( - lam, self.optical_red_Q[:, 0], self.optical_red_Q[:, 1]) - - Qe_blue = Qe_blue*CCD_Qe_lam - Qe_red = Qe_red * CCD_Qe_lam - - # consider the slice optical efficiency in different slicer channel - - img0 = img0*slicer_Qe - ########## do the slice effect ################### - for k in range(32): - - # do slice effect to get slice image - img1 = img0*self.maskSliceHole[str(k)] - - ############ get opd2 and PSF2 ###################### - wavefront2 = self.opd2/lam - psf2 = anySampledPSFnew(wavefront2, self.pupil, Q, 256) - # do convolve - ################## - conv = fftconvolve(img1, psf2, mode='same') - # suppress negative numbers - conv[conv < 0.0] = 0.0 - img2 = conv - - ##############do Slit Mask ########################### - - img2 = img2*self.maskSlitHole[str(k)] - - ######### get opd3 and PSF3 ########################## - - wavefront3 = self.opd3/lam - psf3 = anySampledPSFnew(wavefront3, self.pupil, Q, 256) - - # do convolve - ######################################## - conv = fftconvolve(img2, psf3, mode='same') - # suppress negative numbers - conv[conv < 0.0] = 0.0 - - img3 = conv - - ######################## get subimag ##################### - subimage = galsim.Image(800, 800) - subimage.array[:, :] = img3[int( - sizeout/2)-400:int(sizeout/2)+400, int(sizeout/2)-400:int(sizeout/2)+400] - - # fits.writeto('subimage.fits',subimage.array,overwrite=True) - - ####### get photons from sub-image ##################### - subimage.scale = 0.01 # here pixelscale=0.01arcsec - - subimage.setOrigin(0, 0) - - photons = galsim.PhotonArray.makeFromImage( - subimage, max_flux=max(img3.max()/10000.0, 1)) - - ############################################################## - # do something for each photons; - ############################################################# - idx0 = np.where(photons.flux > 1e-3) - # energy=energy+sum(photons.flux[idx0]) ### totla energy for slice image - p_num = len(idx0[0]) - - ############################################################## - - # find photons for blue channel, and make the flux multiple the optical and CCD efficiency - - if (lam >= 350.0 and lam <= 650.0): - # bulue channel - photons_blue = galsim.PhotonArray( - N=p_num, x=photons.x[idx0], y=photons.y[idx0], flux=Qe_blue*photons.flux[idx0], wavelength=lam) - - dx_blue, dy_blue = get_dx_dy_blue(lam) - - photons_blue.x = photons_blue.x / \ - (self.pixelscale)+dx_blue+self.slice_blue['px'][k] - - photons_blue.y = photons_blue.y / \ - (self.pixelscale)+dy_blue+self.slice_blue['py'][k] - - photons_blue.addTo(blue_img) - - - # fits.writeto('blueImg.fits',blue_img.array,overwrite=True) - - if (lam >= 560.0) & (lam <= 1000.0): - # red channel - photons_red = galsim.PhotonArray( - N=p_num, x=photons.x[idx0], y=photons.y[idx0], flux=Qe_red*photons.flux[idx0], wavelength=lam) - - dx_red, dy_red = get_dx_dy_red(lam) - - photons_red.x = photons_red.x / \ - (self.pixelscale)+dx_red+self.slice_red['px'][k] - - photons_red.y = photons_red.y / \ - (self.pixelscale)+dy_red+self.slice_red['py'][k] - - #red_sensor.accumulate(photons_red, red_img) - - photons_red.addTo(red_img) - - # energy_red=energy_red+sum(photons_red.flux) - - #################################################################################### - - self.image_b = blue_img.array - self.image_r = red_img.array - - now = datetime.utcnow() - - data_time = now.strftime("%Y-%m-%d") - # write the actual file - - bluefile = self.result_path+'/original_Calibration/IFS_'+self.source + \ - '_'+data_time+'_'+str(self.simnumber)+'_B_original.fits' - redfile = self.result_path+'/original_Calibration/IFS_'+self.source + \ - '_'+data_time+'_'+str(self.simnumber)+'_R_original.fits' - - fits.writeto(bluefile, blue_img.array, overwrite=True) - fits.writeto(redfile, red_img.array, overwrite=True) - - #### save flux #### - Flux = flux_1 - - hdu1 = fits.PrimaryHDU(Flux) - - dtime = datetime.utcnow().strftime('%Y -%m -%d %H: %M: %S') - hdu1.header.add_history('Flux data is generated on :'+dtime) - hdu1.header['FluxUnit'] = ('erg/s/cm2/A', 'erg/s/cm2/A') - - hdu2 = fits.ImageHDU(Wave) - - hdu2.header.append(('WaveUnit', 'nm')) - - newd = fits.HDUList([hdu1, hdu2]) - - filename = self.result_path+'/calibration_Data/'+self.source+'_flux.fits' - - if self.simnumber == 1: - newd.writeto(filename, overwrite=True) - - return - ################################################################################## - - def sim_sky_img(self, skyfitsfilename, skyRa, skyDec, sky_rot, telRa, telDec, exposuretime, simnumber): - """ - - - Parameters - ---------- - skyfitsfilename : TYPE - DESCRIPTION. - skyRa : TYPE - DESCRIPTION. - skyDec : TYPE - DESCRIPTION. - sky_rot : TYPE - DESCRIPTION. - telRa : TYPE - DESCRIPTION. - telDec : TYPE - DESCRIPTION. - exposuretime : TYPE - DESCRIPTION. - simnumber : TYPE - DESCRIPTION. - - Returns - ------- - None. - - """ - - # load the input original data cube simulated by Fengshuai - - ############################################################################ - # consider the slice optical efficiency in different slicer channel - slicer_QE_file = self.information['dir_path'] + \ - 'IFS_inputdata/Flatfield/slicer_QE230625.fits' - da = fits.open(slicer_QE_file) - self.log.info('slicer_QE_file path is: %s' % (slicer_QE_file)) - slicer_Qe = da[0].data - - ####################################################################### - - indatafile = skyfitsfilename # 'CSSTMock-M49_HST-FOV10R01-Z00.v202012.fits' - - self.log.info('skyfitsin file is:%s' % (skyfitsfilename)) - a = fits.open(indatafile) - # hdr = a[0].header - # print(a.info()) - # hdr = a[1].header - ####################################################################### - self.information['name_obj'] = a[0].header['OBJECT'] - self.information['ra_obj'] = a[0].header['RA'] # in degree - self.information['dec_obj'] = a[0].header['DEC'] # in degree - - disRa = (telRa-skyRa)/3600.0 # convert unit of arcsec to degree - disDec = (telDec-skyDec)/3600.0 # convert unit of arcsec to degree - # RA_PNT0 =OBJ_RA + DIS_RA/cos(OBJ_DEC), DEC_PNT0= OBJ_DEC DIS_DEC - self.information['ra_pnt0'] = a[0].header['RA'] + \ - disRa/np.cos(a[0].header['DEC']/180.0*np.pi) - self.information['dec_pnt0'] = a[0].header['DEC'] + disDec - - ############### calculate the earthshine and zodiacal noise ,new code 2023.11.1 ############ - ############### - - ra = self.information['ra_pnt0'] - dec = self.information['dec_pnt0'] - - time_jd = time2jd(self.dt) - - if self.appbianpai: - - - sn=self.simnumber-1; - x_sat = float(self.bianpai_data['x_sat'][sn*5+self.exptime_start_index]) - y_sat = float(self.bianpai_data['y_sat'][sn*5+self.exptime_start_index]) - z_sat = float(self.bianpai_data['z_sat'][sn*5+self.exptime_start_index]) - - ### - else: - - x_sat = float(self.orbit_pars[self.orbit_exp_num, 1]) - y_sat = float(self.orbit_pars[self.orbit_exp_num, 2]) - z_sat = float(self.orbit_pars[self.orbit_exp_num, 3]) - - wave0, zodi0 = self.zodiacal( - ra, dec, self.TianCe_day) # erg/s/cm^2/A/arcsec^2 - - # EarthShine from straylight - sl = StrayLight(self.information['dir_path'], jtime=time_jd, sat=np.array([x_sat, y_sat, z_sat]), - radec=np.array([(ra*u.degree).value, (dec*u.degree).value])) - - earth_e = sl.caculateEarthShineFilter(filter='r') - star_e = sl.caculateStarLightFilter(filter='r') - angle_earth = earth_angle(time_jd, x_sat, y_sat, z_sat, ra, dec) - - if angle_earth < 0: - earth_e = 0 - - earthshine_wave0, earthshine_flux0 = ill2flux( - self.information['dir_path'], earth_e+star_e) - - # sample as ifs wavelength - wave_ifs = np.arange(3500, 10000, 1.0) - f1 = interp1d(wave0, zodi0) - zodi_ifs = f1(wave_ifs) - - f2 = interp1d(earthshine_wave0, earthshine_flux0) - earthshine_ifs = f2(wave_ifs) - - self.zodiacal_wave = wave_ifs # in A - self.zodiacal_flux = zodi_ifs - - self.earthshine_wave = wave_ifs # A - self.earthshine_flux = earthshine_ifs - - ############################################################################# - - self.information['target'] = deg2HMS( - self.information['ra_obj'], self.information['dec_obj']) - - # main input data - SpecCube = a[1].data # spectrum data cube; - # Cubesize=len(SpecCube[:,0,0]) - # #print('SpecCube data header', hdr) - # hdr = a[2].header - # the main input data - # the relatived wavelength which is converted from Unit A to nm - Wave = 0.1*a[2].data - #print('Wave data header', hdr) - - ###################################################################################### - - exptime = self.information['exptime'] # exposure time - - dis_dx = ((telRa-skyRa)) # image shift Ra in arcsec - dis_dy = ((telDec-skyDec)) # image shift Dec in arcsec - - sizeout = len(SpecCube[:, 0, 0]) - - blue_img = galsim.Image(np.zeros((2048, 4096)), copy=True) - blue_img.scale = self.pixelscale - blue_img.setOrigin = (0, 0) - - red_img = galsim.Image(np.zeros((3072, 6144)), copy=True) - red_img.scale = self.pixelscale - red_img.setOrigin = (0, 0) - - blue_sensor = galsim.Sensor() - - red_sensor = galsim.Sensor() - - deltalam = np.mean(np.diff(Wave)) - - energy = 0.0 - - energy_blue = 0.0 - energy_red = 0.0 - - width_blue = 0 - ################################ - ############## doppler effect to photons.wavelength ############# - if self.appbianpai: - - sn=self.simnumber-1; - x_sat = float(self.bianpai_data['x_sat'][sn*5+self.exptime_start_index]) - y_sat = float(self.bianpai_data['y_sat'][sn*5+self.exptime_start_index]) - z_sat = float(self.bianpai_data['z_sat'][sn*5+self.exptime_start_index]) - - vx_sat = float(self.bianpai_data['vx_sat'][sn*5+self.exptime_start_index]) - vy_sat = float(self.bianpai_data['vy_sat'][sn*5+self.exptime_start_index]) - vz_sat = float(self.bianpai_data['vz_sat'][sn*5+self.exptime_start_index]) - - else: - - # self.orbit_pars - x_sat = float(self.orbit_pars[self.orbit_exp_num, 1]) - y_sat = float(self.orbit_pars[self.orbit_exp_num, 2]) - z_sat = float(self.orbit_pars[self.orbit_exp_num, 3]) - vx_sat = float(self.orbit_pars[self.orbit_exp_num, 4]) - vy_sat = float(self.orbit_pars[self.orbit_exp_num, 5]) - vz_sat = float(self.orbit_pars[self.orbit_exp_num, 6]) - - self.information['POSI0_X'] = x_sat - self.information['POSI0_Y'] = y_sat - self.information['POSI0_Z'] = z_sat - self.information['VELO0_X'] = vx_sat - self.information['VELO0_Y'] = vy_sat - self.information['VELO0_Z'] = vz_sat - - theta1 = beta_angle(x_sat, y_sat, z_sat, vx_sat, vy_sat, vz_sat, - self.information['ra_obj'], self.information['dec_obj']) - - v1 = np.sqrt(vx_sat**2+vy_sat**2+vz_sat**2)*np.cos(theta1 / - 180.0*np.pi) # velocity at stat exposure time - vv1 = LSR_velocity( - self.information['ra_obj'], self.information['dec_obj'], v1, self.TianCe_day) - - ################################################# - if self.appbianpai: - sn=self.simnumber-1; - - p1x = float(self.bianpai_data['x_sat'][sn*5+self.exptime_end_index]) - p1y = float(self.bianpai_data['y_sat'][sn*5+self.exptime_end_index]) - p1z = float(self.bianpai_data['z_sat'][sn*5+self.exptime_end_index]) - - p1vx = float(self.bianpai_data['vx_sat'][sn*5+self.exptime_end_index]) - p1vy = float(self.bianpai_data['vy_sat'][sn*5+self.exptime_end_index]) - p1vz = float(self.bianpai_data['vz_sat'][sn*5+self.exptime_end_index]) - - else: - - # exposure end time is t2 ; - t2 = self.dt+timedelta(seconds=self.information['exptime']) - ### data read time is the exposure end time plus readouttime ### - # data read end time is t3 - - t2jd = time2jd(t2) - - if self.orbit_pars[-1, 0] < t2jd: # orbit parameters are not in currenct txt file - self.orbit_file_num = self.orbit_file_num+1 - fn = self.information['dir_path'] + \ - 'IFS_inputdata/TianCe/orbit20160925/' + \ - str(self.orbit_file_num)+'.txt' - self.orbit_pars = np.loadtxt(fn) - self.orbit_exp_num = 0 - - for k in range(self.orbit_exp_num, len(self.orbit_pars), 1): - - if t2jd-self.orbit_pars[k, 0] < 0: - break - if k == 0: - deltaT = jd2time(self.orbit_pars[k, 0])-t2 - p1x = self.orbit_pars[k, 1]-(self.orbit_pars[k+1, 1] - - self.orbit_pars[k, 1])*deltaT.seconds/120 - p1y = self.orbit_pars[k, 2]-(self.orbit_pars[k+1, 2] - - self.orbit_pars[k, 2])*deltaT.seconds/120 - p1z = self.orbit_pars[k, 3]-(self.orbit_pars[k+1, 3] - - self.orbit_pars[k, 3])*deltaT.seconds/120 - - p1vx = self.orbit_pars[k, 4]-(self.orbit_pars[k+1, 4] - - self.orbit_pars[k, 4])*deltaT.seconds/120 - p1vy = self.orbit_pars[k, 5]-(self.orbit_pars[k+1, 5] - - self.orbit_pars[k, 5])*deltaT.seconds/120 - p1vz = self.orbit_pars[k, 6]-(self.orbit_pars[k+1, 6] - - self.orbit_pars[k, 6])*deltaT.seconds/120 - - else: - deltaT = jd2time(self.orbit_pars[k, 0])-t2 - p1x = self.orbit_pars[k-1, 1]+(self.orbit_pars[k, 1] - - self.orbit_pars[k-1, 1])*deltaT.seconds/120 - p1y = self.orbit_pars[k-1, 2]+(self.orbit_pars[k, 2] - - self.orbit_pars[k-1, 2])*deltaT.seconds/120 - p1z = self.orbit_pars[k-1, 3]+(self.orbit_pars[k, 3] - - self.orbit_pars[k-1, 3])*deltaT.seconds/120 - - p1vx = self.orbit_pars[k-1, 4]+(self.orbit_pars[k, 4] - - self.orbit_pars[k-1, 4])*deltaT.seconds/120 - p1vy = self.orbit_pars[k-1, 5]+(self.orbit_pars[k, 5] - - self.orbit_pars[k-1, 5])*deltaT.seconds/120 - p1vz = self.orbit_pars[k-1, 6]+(self.orbit_pars[k, 6] - - self.orbit_pars[k-1, 6])*deltaT.seconds/120 - - ####### - self.information['POSI1_X'] = p1x - self.information['POSI1_Y'] = p1y - self.information['POSI1_Z'] = p1z - self.information['VELO1_X'] = p1vx - self.information['VELO1_Y'] = p1vy - self.information['VELO1_Z'] = p1vz - - theta2 = beta_angle(p1x, p1y, p1z, p1vx, p1vy, p1vz, - self.information['ra_obj'], self.information['dec_obj']) - - v2 = np.sqrt(p1vx**2+p1vy**2+p1vz**2)*np.cos(theta2 / - 180.0*np.pi) # velocity at end exposure time - vv2 = LSR_velocity( - self.information['ra_obj'], self.information['dec_obj'], v2, self.TianCe_day) - - # get slice and slit mask - - ######################################################################################## - - if self.debug: - alfa = 1000 - else: - alfa = 1 - - for klam in range(6000): - - ilam = klam*alfa - - if ilam >= 6000: - break - - # print(ilam) - - if self.debug: - print(ilam) - - if ilam % 100 == 0: - - self.log.info('ilam = %i' % ilam) - self.log.info('blue image max value = %i' % - blue_img.array.max()) - self.log.info('red image max value = %i' % red_img.array.max()) - - # print('ilam=',ilam) - - # the wavelength of the i-th frame data - lam = Wave[ilam] - - if lam < 350: - continue - - ############################################### - # get the i-th frame of the input SpecCube data; - Specimg = SpecCube[:, :, ilam] - # convert to photons/cm2/s/A - Nspecimg = SpecCube2photon(Specimg, lam) - Nskynoise = self.CalskyNoise(lam) - - if self.sky_noise: - # add sky noise - Nimg = Nspecimg+Nskynoise - - else: - Nimg = Nspecimg+0.0 - - # multipe the tel area, exposure time, and bandwidh; - # photons/cm2/s/A to photons, here Nimg size is 100*100 - Nimg = Nimg*self.telarea*exptime*deltalam*10.0 - - if Nimg.max() < 0.1: - continue - - if Nimg.min() < 0: - self.log.error('Error in datacube, negative values !!!!!!!') - continue - - ############################################################################## - # if ilam==100: - # fits.writeto('/media/yan//IFSsim/IFSsim Data/original_Img.fits', Nimg, overwrite=True) - - ########################################################################## - # shift image with photons position and rotate them round the image true center - img = galsim.Image(Nimg, copy=True) - img.scale = self.pixelscale - img.setOrigin(0, 0) - - photons = galsim.PhotonArray.makeFromImage( - img, max_flux=max(Nimg.max()/10000.0, 1)) - - # now shift and rotated photons - self.information['shift_dx'] = dis_dx - self.information['shift_dy'] = dis_dy - - photons.x = photons.x-dis_dx # apply shift to photons position in ra direction - photons.y = photons.y-dis_dy # apply shift to photons position in dec direction - # print(dis_dx) - # print(dis_dy) - - x2, y2 = rotation_yan(img.true_center.x, img.true_center.y, - photons.x/img.scale, photons.y/img.scale, sky_rot) - - #print('rotation time=', t2-t1) - - photons2 = galsim.PhotonArray( - N=len(x2), x=x2*img.scale, y=y2*img.scale, flux=photons.flux) - photons = photons2 - # convert the photon image to galsim photons; - ################################################################### - - image0 = galsim.Image(sizeout, sizeout) - image0.scale = self.pixelscale - image0.setOrigin(0, 0) - #image.setCenter(int(np.mean(photons.x)), int(np.mean(photons.y))) - - rotphotons = galsim.PhotonArray( - N=len(photons.x), x=x2, y=y2, flux=photons2.flux) - - sensor = galsim.Sensor() - sensor.accumulate(rotphotons, image0) - - ##################################### - - # ####### - if ilam == 1000: - hdu1 = fits.PrimaryHDU(image0.array[50-32:50+32, 50-32:50+32]) - - hdu1.header['dis_Ra'] = (dis_dx, 'image shift Ra in arcsec') - - hdu1.header['dis_Dec'] = (dis_dy, 'image shift Ra in arcsec') - - hdu1.header['rotAngle'] = ( - sky_rot, 'sky rotation relative to telescope anticlockwise in degree') - - #PSFfilename='/media/yan//IFSsim/IFSsim Data/rot_shift_Img.fits' - - #fits.writeto('/media/yan//IFSsim/IFSsim Data/rot_shift_subImg.fits',image0.array[50-32:50+32,50-32:50+32],overwrite=True ) - hdu1.writeto(self.result_path+'/shift_rot_sky/rot_shift_subImg_' + - str(simnumber)+'.fits', overwrite=True) - fits.writeto(self.result_path+'/shift_rot_sky/original_Img_' + - str(simnumber)+'.fits', Nimg, overwrite=True) - - ##################################################################### - ### do convolve image0 with PSF0 from primay CSST ### - # calculate the PSF0 at this wavelength - - Q = lam*1e-3*self.Fnum/self.pixelsize - wavefront0 = self.opd0/lam - # here we choose reshape=False, rotate the wavefront - wavefront = ndimage.rotate( - wavefront0, -sky_rot, order=1, reshape=False) - - psf0 = anySampledPSFnew(wavefront, self.pupil, Q, 64) - - #psf01=anySampledPSFnew(wavefront, self.pupil , Q , 64) - - conv = fftconvolve(image0.array, psf0, mode='same') - conv[conv < 0.0] = 0.0 - img0 = conv - - # opd1 rms =0.075 @632.8nm - wavefront = self.opd1/lam - # here we choose reshape=False, rotate the wavefront - wavefront = ndimage.rotate( - wavefront, -sky_rot, order=1, reshape=False) - psf1 = anySampledPSFnew(wavefront, self.pupil, Q, 32) - # psf1=ndimage.rotate(psf1, -sky_rot, order=1, reshape=False) # here we choose reshape=False, the rotated image will - - # do convolve with psf1 - conv = fftconvolve(img0, psf1, mode='same') - conv[conv < 0.0] = 0.0 - img0 = conv - - ####################################################### - - # energy=energy+img0[50-32:50+32,50-32:50+32].sum() ## calculate the slice image energy; - - # CCD quantum efficiency - CCD_Qe_lam = np.interp(lam, self.CCD_Qe[:, 0], self.CCD_Qe[:, 1]) - - # optical efficiency , convert the wavelength to A - Qe_blue = np.interp( - lam, self.optical_blue_Q[:, 0], self.optical_blue_Q[:, 1]) - - # optical efficiency, convert the wavelength to A - Qe_red = np.interp( - lam, self.optical_red_Q[:, 0], self.optical_red_Q[:, 1]) - - Qe_blue = Qe_blue*CCD_Qe_lam - - Qe_red = Qe_red * CCD_Qe_lam - - lam1 = lam*(1+vv1/(3.0*1e5)) - - lam2 = lam*(1+vv2/(3.0*1e5)) - - img0 = img0*slicer_Qe - ########## do the slice effect ################### - for k in range(32): - - # do slice effect to get slice image - img1 = img0*self.maskSlice[str(k)] - - ############ get opd2 and PSF2 ###################### - wavefront2 = self.opd2/lam - psf2 = anySampledPSFnew(wavefront2, self.pupil, Q, 64) - # here we choose reshape=False, the rotated image will - psf2 = ndimage.rotate(psf2, -sky_rot, order=1, reshape=False) - # do convolve - ################## - conv = fftconvolve(img1, psf2, mode='same') - # suppress negative numbers - conv[conv < 0.0] = 0.0 - img2 = conv - - ##############do Slit Mask ########################### - - img2 = img2*self.maskSlit[str(k)] - - ######### get opd3 and PSF3 ########################## - - wavefront3 = self.opd3/lam - psf3 = anySampledPSFnew(wavefront3, self.pupil, Q, 64) - # here we choose reshape=False, the rotated image will - psf3 = ndimage.rotate(psf3, -sky_rot, order=1, reshape=False) - - # do convolve - ######################################## - conv = fftconvolve(img2, psf3, mode='same') - # suppress negative numbers - conv[conv < 0.0] = 0.0 - - img3 = conv - - #print('covtime =', t2-t1) - - ######################## get subimage ##################### - subimage = galsim.Image(80, 80) - subimage.array[:, :] = img3[int( - sizeout/2)-40:int(sizeout/2)+40, int(sizeout/2)-40:int(sizeout/2)+40] - - # fits.writeto('subimage.fits',subimage.array,overwrite=True) - - ######################## get photons from sub-image ##################### - subimage.scale = self.pixelscale - - subimage.setOrigin(0, 0) - - photons = galsim.PhotonArray.makeFromImage( - subimage, max_flux=max(img3.max()/10000.0, 1)) - - ############################################################################# - # do something for each photons; - ############################################################################## - idx0 = np.where(photons.flux > 1e-3) - # totla energy for slice image - energy = energy+sum(photons.flux[idx0]) - p_num = len(idx0[0]) - - ############################################################################################### - ####################### code for test slice is right or not ################################### - ############################################################################################### - # photons_2=galsim.PhotonArray(N=p_num, x=photons.x[idx0], y=photons.y[idx0], flux = photons.flux[idx0], wavelength=lam) - # if k==0: - # sliceimage = galsim.Image(100,100) - # sliceimage.scale = pixelscale - # sliceimage.setOrigin(0,0) - # #image.setCenter(int(np.mean(photons.x)), int(np.mean(photons.y))) - - # rotphotons=galsim.PhotonArray(N=p_num, x=photons_2.x/pixelscale, y=photons_2.y/pixelscale, flux = photons_2.flux) - - # sensor = galsim.Sensor() - # sensor.accumulate(rotphotons, sliceimage) - ######################### finish testing ####################################################### - ####################################################################################################### - ##################################### - # find photons for blue channel, and make the flux multiple the optical and CCD efficiency - - np.random.seed(100*self.simnumber) - wavesample = lam1+(lam2-lam1)*np.random.rand(p_num) - - if (lam >= 350.0 and lam <= 650.0): - # bulue channel - photons_blue = galsim.PhotonArray( - N=p_num, x=photons.x[idx0], y=photons.y[idx0], flux=Qe_blue*photons.flux[idx0], wavelength=wavesample) - - dx_blue, dy_blue = get_dx_dy_blue(wavesample) - - photons_blue.x = photons_blue.x/self.pixelscale + \ - dx_blue+self.slice_blue['px'][k] - - photons_blue.y = photons_blue.y/self.pixelscale + \ - dy_blue+self.slice_blue['py'][k] - - blue_sensor.accumulate(photons_blue, blue_img) - - energy_blue = energy_blue+sum(photons_blue.flux) - - - ################# - - # fits.writeto('blueImg.fits',blue_img.array,overwrite=True) - - if (lam >= 560.0) & (lam <= 1000.0): - # red channel - photons_red = galsim.PhotonArray( - N=p_num, x=photons.x[idx0], y=photons.y[idx0], flux=Qe_red*photons.flux[idx0], wavelength=wavesample) - - dx_red, dy_red = get_dx_dy_red(wavesample) - - photons_red.x = photons_red.x/self.pixelscale + \ - dx_red+self.slice_red['px'][k] - - photons_red.y = photons_red.y/self.pixelscale + \ - dy_red+self.slice_red['py'][k] - - red_sensor.accumulate(photons_red, red_img) - - energy_red = energy_red+sum(photons_red.flux) - - ######################################################## - - ########################### test code ################################ - ###################################################################### - # fits.writeto('originalImg.fits',Nimg,overwrite=True) - # fits.writeto('PSFedImg.fits',temp,overwrite=True) - # fits.writeto('rotedphotonsImg.fits',image.array,overwrite=True) - # fits.writeto('newsub_rotedphotonsImg.fits',newimage.array,overwrite=True) - # fits.writeto('slicephotonsImg.fits',sliceimage.array,overwrite=True) - - # pause() - ###################################################################### - ###################### finish test code ########################### - - #################################################################################### - # stray light will cover 2% of input total light; - if self.sky_noise: - - blue_img.array[:, :] = 0.99 * \ - blue_img.array[:, :]+0.01*energy/2048/4096 - - red_img.array[:, :] = 0.99 * \ - red_img.array[:, :] + 0.01*energy/3072/6144 - - self.image_b = blue_img.array - self.image_r = red_img.array - - obsid = simnumber - bluefile = self.result_path+'/original_sky/IFS_' + \ - self.source+'_'+str(obsid)+'_B.fits' - redfile = self.result_path+'/original_sky/IFS_' + \ - self.source+'_'+str(obsid)+'_R.fits' - - fits.writeto(bluefile, blue_img.array, overwrite=True) - fits.writeto(redfile, red_img.array, overwrite=True) - - return - -################################################################################################# - def sim_calibration_img(self, exposuretime, source): - """ - - - Parameters - ---------- - exposuretime : TYPE - DESCRIPTION. - source : TYPE - DESCRIPTION. - - Returns - ------- - None. - - """ - - sizeout = 100 - ################# load the hole mask file ################### - - # load HgAr data; - if self.source == 'LAMP': - matfn3 = self.information['dir_path']+'IFS_inputdata/HgAr.mat' - self.log.info('lamp flux file is : %s' % (matfn3)) - da3 = sio.loadmat(matfn3) - HgArst = da3['HgArst'] # opd unit is in meter - # load flat data - if self.source == 'FLAT': - matfn3 = self.information['dir_path'] + \ - 'IFS_inputdata/flat_light.mat' - self.log.info('flat flux file is : %s' % (matfn3)) - da3 = sio.loadmat(matfn3) - flat_light = da3['flat_light'] # opd unit is in meter - - exptime = self.information['exptime'] # exposure time - - blue_img = galsim.Image(np.zeros((2048, 4096)), copy=True) - blue_img.scale = self.pixelscale - # blue_img.setOrigin=(0,0) - - red_img = galsim.Image(np.zeros((3072, 6144)), copy=True) - red_img.scale = self.pixelscale - # red_img.setOrigin=(0,0) - - blue_sensor = galsim.Sensor() - - red_sensor = galsim.Sensor() - - slice_image = np.zeros((sizeout, sizeout)) - - LightSt = np.zeros(6000) - flux_1 = np.zeros(6000) - Width_lambda = (1000-350)/6000 - Wave = np.zeros(6000) - - # code for slicer position #########33 - #### - - if self.debug: - alfa = 1000 - else: - alfa = 1 - - for klam in range(6000): - - ilam = klam*alfa - - if self.debug: - print(ilam) - - # ilam=ilam*200 - if ilam % 100 == 0: - self.log.info('ilam = %i' % ilam) - self.log.info('blue image max value = %i' % - blue_img.array.max()) - self.log.info('red image max value = %i' % red_img.array.max()) - - if ilam >= 6000: - break - - # the wavelength of the i-th frame data - lam = 350.0+(1000-350)/6000.0 * ilam - - Wave[ilam] = lam - - # tab['wavelength'][ilam]=lam - #print('i = ',ilam ) - ############################################### - if self.source == 'LAMP': - - sumt = 0.0 - - for k in range(11): - t = getSpectrum(HgArst, lam+(k-5) * - Width_lambda/11, self.HgArsigma) - sumt = sumt+t[0]/11 - #print('t=', t) - - # the normalized HgAr light source in the range of 0 and 1.0 - LightSt[ilam] = sumt - - skymag = 9 # mag/arcsec^2 unit - - if ilam == 0: - # fluxnu1,fluxlam1,fluxphoton1=mag2flux(skymag,633) - fluxlam1, fluxphoton1 = mag2flux(skymag, 633) - - # save flux information - flux_1[ilam] = LightSt[ilam] * \ - fluxlam1*self.fov2 # erg/s/cm2/A - - nphoton = flux2photons(flux_1[ilam], lam) - - # photons/cm2/s/A/arcsec2 - ns = nphoton*np.ones((sizeout, sizeout)) - - slice_image = ns*self.telarea*Width_lambda*10*exptime - -# #####################3 transfer: photons/cm2/s/A/arcsec2 to photons - if self.source == 'FLAT': - ########## - flatS = np.interp(lam, flat_light[:, 0], flat_light[:, 1]) - - # the normalized HgAr light source in the range of 0 and 1.0 - LightSt[ilam] = flatS - - skymag = 9 # mag/arcsec^2 unit - - if ilam == 0: - fluxlam1, fluxphoton1 = mag2flux(skymag, 633) - - # save flux information - flux_1[ilam] = LightSt[ilam] * \ - fluxlam1*self.fov2 # erg/s/cm2/A - - nphoton = flux2photons(flux_1[ilam], lam) - - # photons/cm2/s/A/arcsec2 - ns = nphoton*np.ones((sizeout, sizeout)) - - slice_image = ns*self.telarea*Width_lambda*10*exptime - -############################################################################ - - if slice_image.max() < 0.1: - continue - - if slice_image.min() < 0: - self.log.error('Error in datacube, negative values !!!!!!!') - continue -############################################################################# - - # for calibration , there is no primay CSST system, but diffraction effect should consider ; - - Q = lam*1e-3*self.Fnum/self.pixelsize - - wavefront0 = self.opd0/lam/100000.0 - - psf0 = anySampledPSFnew(wavefront0, self.pupil, Q, 64) - - #psf01=anySampledPSFnew(wavefront, self.pupil , Q , 64) - - conv = fftconvolve(slice_image, psf0, mode='same') - - conv[conv < 0.0] = 0.0 - - slice_image2 = conv - - img0 = slice_image2 - - ########################################## - # opd1 rms =0.075 @632.8nm - wavefront = self.opd1/lam - - psf1 = anySampledPSFnew(wavefront, self.pupil, Q, 64) - - # do convolve with psf1 - conv = fftconvolve(img0, psf1, mode='same') - conv[conv < 0.0] = 0.0 - img0 = conv - - # CCD quantum efficiency - CCD_Qe_lam = np.interp(lam, self.CCD_Qe[:, 0], self.CCD_Qe[:, 1]) - - # optical efficiency , convert the wavelength to A - Qe_blue = np.interp( - lam, self.optical_blue_Q[:, 0], self.optical_blue_Q[:, 1]) - - # optical efficiency, convert the wavelength to A - Qe_red = np.interp( - lam, self.optical_red_Q[:, 0], self.optical_red_Q[:, 1]) - - Qe_blue = Qe_blue*CCD_Qe_lam - Qe_red = Qe_red * CCD_Qe_lam - - # print('time=',(end-start)) - # consider the slice optical efficiency in different slicer channel - da = fits.open( - self.information['dir_path']+'IFS_inputdata/Flatfield/slicer_QE230625.fits') - slicer_Qe = da[0].data - img0 = img0*slicer_Qe - ########## do the slice effect ################### - for k in range(32): - - # do slice effect to get slice image - img1 = img0*self.maskSlice[str(k)] - - ############ get opd2 and PSF2 ###################### - wavefront2 = self.opd2/lam - psf2 = anySampledPSFnew(wavefront2, self.pupil, Q, 64) - # do convolve - ################## - conv = fftconvolve(img1, psf2, mode='same') - # suppress negative numbers - conv[conv < 0.0] = 0.0 - img2 = conv - - ##############do Slit Mask ########################### - - img2 = img2*self.maskSlit[str(k)] - - ######### get opd3 and PSF3 ########################## - - wavefront3 = self.opd3/lam - psf3 = anySampledPSFnew(wavefront3, self.pupil, Q, 64) - - # do convolve - ######################################## - conv = fftconvolve(img2, psf3, mode='same') - # suppress negative numbers - conv[conv < 0.0] = 0.0 - - img3 = conv - - ######################## get subimage ##################### - subimage = galsim.Image(80, 80) - subimage.array[:, :] = img3[int( - sizeout/2)-40:int(sizeout/2)+40, int(sizeout/2)-40:int(sizeout/2)+40] - - # fits.writeto('subimage.fits',subimage.array,overwrite=True) - - ######################## get photons from sub-image ##################### - subimage.scale = self.pixelscale - - subimage.setOrigin(0, 0) - - photons = galsim.PhotonArray.makeFromImage( - subimage, max_flux=max(img3.max()/10000.0, 1)) - - ############################################################################# - # do something for each photons; - ############################################################################## - idx0 = np.where(photons.flux > 1e-2) - # energy=energy+sum(photons.flux[idx0]) ### totla energy for slice image - p_num = len(idx0[0]) - - ############################################################################################### - ####################### code for test slice is right or not ################################### - ############################################################################################### - # photons_2=galsim.PhotonArray(N=p_num, x=photons.x[idx0], y=photons.y[idx0], flux = photons.flux[idx0], wavelength=lam) - # if k==0: - # sliceimage = galsim.Image(100,100) - # sliceimage.scale = pixelscale - # sliceimage.setOrigin(0,0) - # #image.setCenter(int(np.mean(photons.x)), int(np.mean(photons.y))) - - # rotphotons=galsim.PhotonArray(N=p_num, x=photons_2.x/pixelscale, y=photons_2.y/pixelscale, flux = photons_2.flux) - - # sensor = galsim.Sensor() - # sensor.accumulate(rotphotons, sliceimage) - # ######################## finish testing ####################################################### - # ###################################################################################################### - - ##################################### - # find photons for blue channel, and make the flux multiple the optical and CCD efficiency - - if (lam >= 350.0 and lam <= 650.0): - # bulue channel - photons_blue = galsim.PhotonArray( - N=p_num, x=photons.x[idx0], y=photons.y[idx0], flux=Qe_blue*photons.flux[idx0], wavelength=lam) - - dx_blue, dy_blue = get_dx_dy_blue(lam) - - photons_blue.x = photons_blue.x/self.pixelscale + \ - dx_blue+self.slice_blue['px'][k] - - photons_blue.y = photons_blue.y/self.pixelscale + \ - dy_blue+self.slice_blue['py'][k] - - #photons_blue.y =2048-photons_blue.y - - blue_sensor.accumulate(photons_blue, blue_img) - - # energy_blue=energy_blue+sum(photons_blue.flux) - - # width_blue=width_blue+Width_lambda/32.0 - - # fits.writeto('blueImg.fits',blue_img.array,overwrite=True) - - if (lam >= 560.0) & (lam <= 1000.0): - # red channel - photons_red = galsim.PhotonArray( - N=p_num, x=photons.x[idx0], y=photons.y[idx0], flux=Qe_red*photons.flux[idx0], wavelength=lam) - - dx_red, dy_red = get_dx_dy_red(lam) - - photons_red.x = photons_red.x/self.pixelscale + \ - dx_red+self.slice_red['px'][k] - - photons_red.y = photons_red.y/self.pixelscale + \ - dy_red+self.slice_red['py'][k] - - #photons_red.y =3072-photons_red.y - - red_sensor.accumulate(photons_red, red_img) - - # energy_red=energy_red+sum(photons_red.flux) - - #################################################################################### - - self.image_b = blue_img.array - self.image_r = red_img.array - - now = datetime.utcnow() - - data_time = now.strftime("%Y-%m-%d") - # write the actual file - - bluefile = self.result_path+'/original_Calibration/IFS_'+self.source + \ - '_'+data_time+'_'+str(self.simnumber)+'_B_original.fits' - redfile = self.result_path+'/original_Calibration/IFS_'+self.source + \ - '_'+data_time+'_'+str(self.simnumber)+'_R_original.fits' - - fits.writeto(bluefile, blue_img.array, overwrite=True) - fits.writeto(redfile, red_img.array, overwrite=True) - - #### save flux #### - Flux = flux_1 - - hdu1 = fits.PrimaryHDU(Flux) - - dtime = datetime.utcnow().strftime('%Y -%m -%d %H: %M: %S') - hdu1.header.add_history('Flux data is generated on :'+dtime) - hdu1.header['FluxUnit'] = ('erg/s/cm2/A', 'erg/s/cm2/A') - - hdu2 = fits.ImageHDU(Wave) - - hdu2.header.append(('WaveUnit', 'nm')) - - newd = fits.HDUList([hdu1, hdu2]) - - filename = self.result_path+'/calibration_Data/'+self.source+'_flux.fits' - - if self.simnumber == 1: - newd.writeto(filename, overwrite=True) - - return - - ########################################################################## - # def saveCCDfullimage(self): - - # # - # bluefile=self.result_path+'/original_sky/IFS_CCD_'+self.source+'_'+str(self.simnumber)+'_B.fits' - # redfile =self.result_path+'/original_sky/IFS_CCD_'+self.source+'_'+str(self.simnumber)+'_R.fits' - - # fits.writeto(bluefile, self.image_b, overwrite=True) - # fits.writeto(redfile, self.image_r, overwrite=True) -################################################################################################# - - def simulate(self, sourcein, simnumber): - """ - - - Parameters - ---------- - sourcein : TYPE - DESCRIPTION. - simnumber : TYPE - DESCRIPTION. - - Returns - ------- - None. - - """ - # def simulate(self, skyfitsfilename, skyRa, skyDec, sky_rot, telRa, telDec, exposuretime,simnumber): - """ - Create a single simulated image of a quadrant defined by the configuration file. - Will do all steps defined in the config file sequentially. - - :return: None - """ - - self.source = sourcein - - self.simnumber = simnumber - - #self.configure(simnumber,dir_path,result_path) # print the configfile name and path; - - self.debug = self.information['debug'] - - - - - if self.information['exptime']>2000 or self.information['exptime']<0: - self.log.error( - 'exptime value error, it shoub be in [0, 2000]!') - exit(2) - - if self.source == 'SCI' or self.source == 'COMP': - - if simnumber <= 50 and simnumber > 20: - self.information['exptime'] = 300 - - if simnumber <= 100 and simnumber > 50: - self.information['exptime'] = 600 - - if simnumber <= 150 and simnumber > 100: - self.information['exptime'] = 900 - - if simnumber <= 200 and simnumber > 150: - self.information['exptime'] = 1200 - - self.skyfilepath = self.information['dir_path'] + \ - self.information['sky_fitsin'] - - print('self.skyfilepath = ', self.skyfilepath) - - ################################################################### - if self.appbianpai: - ### load yunxingbianpai csv file - ############ load star data catlog ##################### - starcat=self.information['bianpai_file'] - - ##starcat='selection_20230517_concat.fits' - ################################################### - - self.log.info('Stat catlog file name is %s' % (starcat)) - - ########################################## - df=pd.read_csv(self.information['dir_path']+'IFS_inputdata/TianCe/'+starcat) - - ################################################################### - sn=self.simnumber-1; - arr=np.array(df['time'][sn*5:sn*5+5]); - self.exptime_start_jd,self.exptime_start_index=find_min(arr); - self.exptime_end_jd, self.exptime_end_index =find_max(arr); - - ###self.earthshine_theta=df['earth_angle'][sn*5+index] # in degree - self.dt = jd2time(self.exptime_start_jd); - self.bianpai_data=df; - ################################################################## - else: - - #### load orbit parameters ##### - flag = 0 - self.dt = datetime.utcnow() - self.dt_num = int( - self.simnumber*(self.information['exptime']+self.information['readouttime']+125)/120) - now_dt = datetime.utcnow() - now_jd = time2jd(now_dt) - - - - for k in range(1, 50, 1): - - # fn=father_path+'/IFS_inputdata/TianCe/orbit20160925/'+str(k)+'.txt'; - - fn = self.information['dir_path'] + \ - 'IFS_inputdata/TianCe/orbit20160925/'+str(k)+'.txt' - d = np.loadtxt(fn) - - for kk in range(len(d[:, 0])): - if now_jd-d[kk, 0] <= 0: - flag = 1 - break - if flag == 1: - break - - ################################################################## - - if kk + self.dt_num < len(d[:, 0]): - - exptime_start_jd = d[kk+self.dt_num, 0] - - self.orbit_pars = d - self.orbit_file_num = k - self.orbit_exp_num = kk - - else: - - fn = self.information['dir_path'] + \ - 'IFS_inputdata/TianCe/orbit20160925/'+str(k+1)+'.txt' - d = np.loadtxt(fn) - - self.orbit_pars = d - self.orbit_file_num = k+1 - self.orbit_exp_num = self.dt_num - self.dt = jd2time(exptime_start_jd) - - ################################################################## - ################################################################## - - - self.TianCe_day = self.dt.strftime("%Y-%m-%d") - - self.TianCe_exp_start = dt2hmd(self.dt) - - self.zodiacal_time = self.TianCe_day - ################################################################### - - skyRa = 0 - skyDec = 0 - - if simnumber > 1: - np.random.seed(self.simnumber) - ud = np.random.random() # Choose a random rotation in degree - sky_rot = 2 * (ud-0.5) * 5 - - dsmax = np.floor(50*np.cos(sky_rot/180*np.pi))-34 - np.random.seed(11*self.simnumber) - ud = np.random.random() # Choose a random shift in arcsec - telRa = 2 * (ud-0.5) * dsmax * 0.1 # in arcsec - - np.random.seed(123*self.simnumber) - ud = np.random.random() # Choose a random shift in arcsec - telDec = 2 * (ud-0.5) * dsmax * 0.1 # in arcsec - - sky_rot = 0.0 - telRa = 0 - telDec = 0 - - self.sim_sky_img(self.skyfilepath, skyRa, skyDec, 0*sky_rot, - 0*telRa, 0*telDec, self.information['exptime'], simnumber) - - self.skyRa = skyRa - self.skyDec = skyDec - self.sky_rot = sky_rot - - self.telRa = telRa - self.telDec = telDec - - self.information['sky_rot'] = sky_rot - - ####################################################################### - elif self.source == 'FLAT': - self.dt = datetime.utcnow() - self.information['exptime'] = 200+simnumber*100 - self.sim_calibration_img(self.information['exptime'], 'FLAT') - self.information['sky_rot'] = 0 - - ######## - - elif self.source == 'LAMP': - self.dt = datetime.utcnow() - self.information['exptime'] = 200+simnumber*100 - self.information['sky_rot'] = 0 - - if self.information['holemask'] == 'yes': - self.information['exptime'] = 300 - self.sim_lamp_hole_img(self.information['exptime']) - else: - self.sim_calibration_img(self.information['exptime'], 'LAMP') - - ######### - elif self.source == 'DARK': - self.dt = datetime.utcnow() - self.information['sky_rot'] = 0 - self.information['exptime'] = simnumber*3600*24 - self.image_b = np.zeros((2048, 4096)) - self.image_r = np.zeros((3072, 6144)) - - elif self.source == 'BIAS': - self.dt = datetime.utcnow() - self.information['sky_rot'] = 0 - self.information['exptime'] = 0 - self.image_b = np.zeros((2048, 4096)) - self.image_r = np.zeros((3072, 6144)) - - else: - print('Souce is not correct and programe will return') - raise ValueError('Souce is not correct and programe will return') -############################################################################### - # save the original image firstly; - self.image_b_ori = self.image_b+10.0 - self.image_r_ori = self.image_r+10.0 - - exptimes = self.information['exposuretimes'] - - if self.source == 'SCI' or self.source == 'COMP': - exptimes = self.information['exposuretimes'] - else: - exptimes = 1 - self.information['exposuretimes'] = 1 - - if self.debug: - exptimes = 1 - - if exptimes > 1: - self.image_b_N = np.zeros((1344, 9672, exptimes)) - - self.image_r_N = np.zeros((1856, 13768, exptimes)) - else: - self.image_b_N = np.zeros((1344, 9672)) - - self.image_r_N = np.zeros((1856, 13768)) - - for idk in range(exptimes): - - ############ add some effect to images ####################################### - - self.image_b = self.image_b_ori-10.0 - self.image_r = self.image_r_ori-10.0 - - self.image_b=np.zeros((2048, 4096)) - self.image_b[500,500:550]=150000.0 - self.image_b[500,555]=250000.0 - - if self.information['exptime'] > 0.09216: - - self.getFrameTransferImg() - print('Applying getFrameTransferImg finished.......') - - if self.flatfieldM: - self.applyflatfield(simnumber) - print('Applying flatfieldM finished.......') - - fits.writeto('after_flatfieldM.fits', self.image_b, overwrite=True) - - # if self.source == 'SCI' or self.source == 'COMP': - - # if self.cosmicRays: - # self.addCosmicRays(idk) - # print('Applying cosmicRays finished.......') - # fits.writeto('after_cosmicRays.fits', self.image_b, overwrite=True) - - if self.bleeding: - self.applyBleeding_yan() - print('Applying bleeding finished.......') - - fits.writeto('after_bleeding.fits', self.image_b, overwrite=True) - - - if self.darknoise: - self.applyDarkCurrent() - print('Applying dark noise finished.......') - - fits.writeto('after_darknoise.fits', self.image_b, overwrite=True) - - - ################################## - - if self.source != 'SCI' and self.source != 'COMP': - self.applyPoissonNoise() - print('Applying Possion noise finished.......') - - fits.writeto('after_PossionNoise.fits', self.image_b, overwrite=True) - - - ################################# - - if self.nonlinearity: - self.applyNonlinearity() - print('Applying nonlinearity finished.......') - - fits.writeto('after_Nonlinearity.fits', self.image_b, overwrite=True) - - - - # cut original CCD image to four parts by four read out channels - self.CCDreadout() - - fits.writeto('after_CCDreadout.fits', self.image_b, overwrite=True) - - - - if self.radiationDamage: - self.applyRadiationDamage() - print('Applying radiationDamage finished.......') - - fits.writeto('after_RadiationDamage.fits', self.image_b, overwrite=True) - - - if self.readoutNoise: - self.applyReadoutNoise() - print('Applying readoutNoise finished.......') - - fits.writeto('after_readoutNoise.fits', self.image_b, overwrite=True) - - - if self.information['exptime'] > 0.09216: - self.appFrameTransferEffect() - print('Applying Frame Transfer Effect finished.......') - - fits.writeto('after_FrameTransfer.fits', self.image_b, overwrite=True) - - - self.electrons2ADU() - - fits.writeto('after_ADU.fits', self.image_b, overwrite=True) - - - self.applyBias() - - fits.writeto('after_Bias.fits', self.image_b, overwrite=True) - - - if self.cosmetics: - self.applyCosmetics() - print('Applying cosmetics finished.......') - - fits.writeto('after_Cosmetics.fits', self.image_b, overwrite=True) - - self.discretise() - fits.writeto('after_discretise.fits', self.image_b, overwrite=True) - - - if exptimes > 1: - - self.image_b_N[:, :, idk] = self.image_b[:, :].copy() - self.image_r_N[:, :, idk] = self.image_r[:, :].copy() - else: - self.image_b_N[:, :] = self.image_b[:, :].copy() - self.image_r_N[:, :] = self.image_r[:, :].copy() - - self.writeOutputs(simnumber) - - fits.writeto('after_writeOutputs.fits', self.image_b, overwrite=True) - - ################################################# - - self.log.info('Using the following input values:') - - for key, value in self.information.items(): - self.log.info('%s = %s' % (key, value)) - self.log.info('Using the following booleans:') - - for key, value in self.booleans.items(): - self.log.info('%s = %s' % (key, value)) - - self.log.info('Finished the ith_Exposure = %i' % (simnumber)) -# print('The iLoop= % d simlaiton finished. ' %simnumber) - -############################################################################ -############################################################################ - - -def runIFSsim(sourcein, configfile, dir_path, result_path, iLoop, debug, applyhole): - """ - - - Parameters - ---------- - sourcein : TYPE - DESCRIPTION. - configfile : TYPE - DESCRIPTION. - iLoop : TYPE - DESCRIPTION. - applyhole : TYPE, optional - DESCRIPTION. The default is 'no'. - - Returns - ------- - int - DESCRIPTION. - - """ - - simulate = dict() - simulate[iLoop] = IFSsimulator(configfile) - - simulate[iLoop].configure(sourcein, dir_path, result_path, iLoop, debug, applyhole) # load the configfile; - - - if applyhole == 'yes' and sourcein == 'LAMP': - simulate[iLoop].information['holemask'] = 'yes' - else: - simulate[iLoop].information['holemask'] = 'no' - - ### dir_path = os.path.join(os.environ['UNIT_TEST_DATA_ROOT'], 'ifs_sim/') - # simulate[iLoop].information['dir_path'] = dir_path - - # simulate[iLoop].information['debug'] = debug - - # simulate[iLoop].information['result_path'] = result_path - - simulate[iLoop].simulate(sourcein, iLoop) - - return 1 - -############################## end ########################################## -- GitLab