diff --git a/.DS_Store b/.DS_Store new file mode 100644 index 0000000000000000000000000000000000000000..f02c8b3595aebcacb892df267db6274e320c77ba Binary files /dev/null and b/.DS_Store differ diff --git a/csst_ifs_sim/csst_ifs_sim.py b/csst_ifs_sim/csst_ifs_sim.py index dfd300fda1716d014f74745ccc92f4aabec30fcb..c06caa7dcad92693b5f9e29e0b5baa306254bc8d 100644 --- a/csst_ifs_sim/csst_ifs_sim.py +++ b/csst_ifs_sim/csst_ifs_sim.py @@ -1,12 +1,49 @@ #!/usr/bin/env python3 # -*- coding: utf-8 -*- """ -Created on Mon Mar 7 14:53:41 2022 +Created on Thu Apr 11 15:18:57 2024 @author: yan """ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +# import sys +import os +from scipy.interpolate import interp1d +import astropy.coordinates as coord +import ctypes +from csst_ifs_sim.support import IFSinstrumentModel +from csst_ifs_sim.support import cosmicrays +from csst_ifs_sim.support import logger as lg +from csst_ifs_sim.CTI import CTI +# from optparse import OptionParser +import configparser as ConfigParser +import cmath +from scipy import ndimage +from scipy.signal import fftconvolve +from astropy.io import fits +import scipy.io as sio +import numpy as np +from scipy import interpolate +from astropy.coordinates import SkyCoord, EarthLocation +from astropy import units as u +from astropy.time import Time +from astropy.coordinates import get_sun +from datetime import datetime, timedelta +import pandas as pd +from scipy.integrate import simps +import galsim +# sys.path.append('../') +from joblib import Parallel, delayed +from astropy.utils.iers import conf +conf.auto_max_age = None + +""" +""" """ +Contact Information: zhaojunyan@shao.ac.cn + The CSST IFS Image Simulator ============================================= @@ -40,11 +77,14 @@ a star or an extended source (i.e. a galaxy). (first step), and finally overlay onto the detector according to its position. #. Apply calibration unit flux to mimic flat field exposures [optional]. -#. Apply a multiplicative flat-field map to emulate pixel-to-pixel non-uniformity [optional]. +#. Apply a multiplicative flat-field map to emulate pixel-to-pixel +non-uniformity [optional]. #. Add a charge injection line (horizontal and/or vertical) [optional]. -#. Add cosmic ray tracks onto the CCD with random positions but known distribution [optional]. +#. Add cosmic ray tracks onto the CCD with random positions but known +distribution [optional]. #. Apply detector charge bleeding in column direction [optional]. -#. Add constant dark current and background light from Zodiacal light [optional]. +#. Add constant dark current and background light from Zodiacal +light [optional]. #. Include spatially uniform scattered light to the pixel grid [optional]. #. Add photon (Poisson) noise [optional] #. Add cosmetic defects from an input file [optional]. @@ -53,135 +93,52 @@ a star or an extended source (i.e. a galaxy). #. Apply CCD273 non-linearity model to the pixel data [optional]. #. Add readout noise selected from a Gaussian distribution [optional]. #. Convert from electrons to ADUs using a given gain factor. -#. Add a given bias level and discretise the counts (the output is going to be in 16bit unsigned integers). +#. Add a given bias level and discretise the counts (the output is going to be +in 16bit unsigned integers). #. Finally the simulated image is converted to a FITS file, a WCS is assigned and the output is saved to the current working directory. -.. 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. - - -Dependencies ------------- - -This script depends on the following packages: - -:requires: PyFITS (tested with 3.0.6 and 3.3) -:requires: NumPy (tested with 1.6.1, 1.7.1, 1.8.0, 1.9.1, and 1.9.2) -:requires: numexpr (tested with 2.0.1 and 2.3.1) -:requires: SciPy (tested with 0.10.1, 0.12, 0.13, 0.14, and 0.15.1) -:requires: vissim-python package - -.. Note:: This class is Python 3 compatible. - -.. Note:: CUDA acceleration requires an NVIDIA GPU that supports CUDA and PyFFT and PyCUDA packages. - Note that the CUDA acceleration does not speed up the computations unless oversampled PSF - is being used. If > 2GB of memory is available on the GPU, speed up up to a factor of 50 is - possible. - -Testing -------- - -the input confiefile as follows: - -config_file=['IFS_test_210412.config' ] - -you can modify the configfile to change the simulaiton parameters. - -please run the as follows:: - - python IFS_sim_Pro.py - -This will run the unittest and compare the result to a previously calculated image. -Please inspect the standard output for results. - -Running the test will produce an image representing IFS lower left (0th) quadrant of the CCD (x, y) = (0, 0). Because -noise and cosmic rays are randomised one cannot directly compare the science -outputs but we must rely on the outputs that are free from random effects. The data subdirectory -contains a file named "NoNoiseNoCr_IFS_b.fits, NoNoiseNoCr_IFS_r.fits, NoNoiseNoCr_IFS_i.fits", which is the comparison image without -any noise or cosmic rays. +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. -:version: 2022.12.15 +Note:: This class is Python 3 compatible. - -@230613 ,add new lamp hole simulation - -Contact Information: zhaojunyan@shao.ac.cn - -------- - ---23.4.10 --- update the fits header as defined order 0 data format - -2023.12.26 add date frame transfer effect and software version from the config file - -2024.1.19... add multiple exposure mode function; +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 """ -import sys -import os -from scipy.interpolate import interp1d -import astropy.coordinates as coord -import ctypes - -from csst_ifs_sim.support import IFSinstrumentModel -from csst_ifs_sim.support import cosmicrays -from csst_ifs_sim.support import logger as lg - - -from csst_ifs_sim.CTI import CTI -from optparse import OptionParser -import configparser as ConfigParser -import cmath -from scipy import ndimage -from scipy.signal import fftconvolve -from astropy.io import fits -import scipy.io as sio -import numpy as np -from scipy import interpolate -from astropy.coordinates import SkyCoord, EarthLocation -from astropy import units as u -from astropy.time import Time -from astropy.coordinates import get_sun -from datetime import datetime, timedelta -import pandas as pd -from scipy.integrate import simps -from astropy.utils.iers import conf -conf.auto_max_age = None - -import galsim -###import julian - - -sys.path.append('../') - - -#from joblib import Parallel, delayed # from astropy.table import Table ############################################################################### - -#from sky_bkg import sky_bkg - - filterPivotWave = {'nuv': 2875.5, 'u': 3629.6, 'g': 4808.4, 'r': 6178.2, 'i': 7609.0, 'z': 9012.9, 'y': 9627.9} filterIndex = {'nuv': 0, 'u': 1, 'g': 2, 'r': 3, 'i': 4, 'z': 5, 'y': 6} -# filterCCD = {'nuv':'UV0','u':'UV0','g':'Astro_MB','r':'Astro_MB', 'i':'Basic_NIR', 'z':'Basic_NIR','y':'Basic_NIR'} -# bandRange = {'nuv':[2504.0,3230.0],'u':[3190.0,4039.0],'g':[3989.0,5498.0],'r':[5438.0,6956.0], 'i':[6886.0,8469.0], -# 'z':[8379.0,10855.0],'y':[9217.0, 10900.0], 'GU':[2550, 4000],'GV':[4000, 6200],'GI':[6200,10000]} -# Instrument_dir = '/Users/linlin/Data/csst/straylightsim-master/Instrument/' -# SpecOrder = ['-2','-1','0','1','2'] -# -# filterMirrorEff = {'nuv':0.54,'u':0.68,'g':0.8,'r':0.8, 'i':0.8, 'z':0.8,'y':0.8} def transRaDec2D(ra, dec): - # radec转为竞天程序里的ob, 赤道坐标系下的笛卡尔三维坐标xyz. + """ + + 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) @@ -190,6 +147,22 @@ def transRaDec2D(ra, dec): def flux2ill(wave, flux): + """ + + + Parameters + ---------- + wave : TYPE + DESCRIPTION. + flux : TYPE + DESCRIPTION. + + Returns + ------- + E : TYPE + DESCRIPTION. + + """ # erg/s/cm^2/A/arcsec^2 to W/m^2 @@ -208,22 +181,6 @@ def flux2ill(wave, flux): f = 28 # meter flux3 = flux2 * np.pi * D**2 / 4 / f**2 / 10**3 - # # # u_lambda: W/m^2/nm - # # obj0: - # V = 73 * 1e-8 # W/m^2/nm - # Es = V * np.pi * D**2 / 4 / f**2 / 10**3 - # c1 = 3.7418e-16 - # c2 = 1.44e-2 - # t = 5700 - # # wave需要代入meter. - # wave0 = np.arange(380, 780) # nm - # delta_lamba0 = 1 # nm - # u_lambda = c1 / (wave0*1e-9)**5 / (np.exp(c2 / (wave0*1e-9) / t) - 1) - # f_lambda = (u_lambda / u_lambda[wave0 == 500]) * Es - # E0 = np.sum(f_lambda * 1) - # # plt.plot(wave, u_lambda) - # # plt.show() - # 对波长积分 f = interp1d(wave, flux3) wave_interp = np.arange(3800, 7800) @@ -232,14 +189,30 @@ def flux2ill(wave, flux): delta_lamba = 0.1 # nm E = np.sum(flux3_interp * delta_lamba) - # pdb.set_trace() - return E ################################################################ 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' @@ -270,10 +243,33 @@ def ill2flux(path, E): ############################################################## - ########################################################################################## 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 @@ -299,8 +295,7 @@ def earth_angle(time_jd, x_sat, y_sat, z_sat, ra_obj, dec_obj): return angle - -# 3 +##################################################################### class StrayLight(object): @@ -344,8 +339,22 @@ class StrayLight(object): 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. + + """ sat = (ctypes.c_double*3)() sat[:] = self.sat @@ -367,8 +376,22 @@ class StrayLight(object): 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. + + """ sat = (ctypes.c_double*3)() sat[:] = self.sat ob = (ctypes.c_double*3)() @@ -379,8 +402,8 @@ class StrayLight(object): 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个波段的照度 + 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) @@ -395,23 +418,65 @@ class StrayLight(object): def str2time(strTime): + """ + + + Parameters + ---------- + strTime : TYPE + DESCRIPTION. + + Returns + ------- + TYPE + DESCRIPTION. + + """ if len(strTime) > 20: # 暂时未用到 msec = int(float('0.'+strTime[20:])*1000000) # 微秒 str2 = strTime[0:19]+' '+str(msec) return datetime.strptime(str2, '%Y %m %d %H %M %S %f') # datetime类转mjd - +########################################################################## 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 + \ @@ -419,29 +484,99 @@ def time2jd(dateT): 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): - return mjd+2400000.5 + """ + Parameters + ---------- + mjd : TYPE + DESCRIPTION. + + Returns + ------- + TYPE + DESCRIPTION. + + """ + return mjd+2400000.5 + +########################################################################### def jd2mjd(jd): - return jd-2400000.5 + """ + + + 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 @@ -464,9 +599,23 @@ def dt2hmd(dt): return str_h+':'+str_m+':'+str_d ############################################################################### +def float2char(a, z): + """ -def float2char(a, z): + Parameters + ---------- + a : TYPE + DESCRIPTION. + z : TYPE + DESCRIPTION. + + Returns + ------- + TYPE + DESCRIPTION. + + """ # a is input float value # transfer float a to chars and save z bis b = str(a) @@ -481,6 +630,22 @@ def float2char(a, z): ############################################################################### 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') @@ -517,9 +682,35 @@ def deg2HMS(ra0, dec0): ############################################################################ - -########################################################### 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 @@ -549,10 +740,28 @@ def beta_angle(x_sat, y_sat, z_sat, vx_sat, vy_sat, vz_sat, ra_obj, dec_obj): return angle ############################################################################### - - ########################################################## 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为利用输入坐标建立的坐标系对象。 @@ -580,6 +789,30 @@ def LSR_velocity(ra, dec, velocity, Obstime): 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 @@ -591,6 +824,20 @@ def rotation_yan(x0, y0, x1, y1, angle): 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) @@ -600,10 +847,24 @@ def centroid(data): cy = (np.dot(np.dot(x, data), y1)) / (np.dot(np.dot(x1, data), y1)) return float(cx), float(cy) -# 33 +############################################################################### +def centroidN(data): + """ -def centroidN(data): + Parameters + ---------- + data : TYPE + DESCRIPTION. + + Returns + ------- + cx : TYPE + DESCRIPTION. + cy : TYPE + DESCRIPTION. + + """ ''' calculate the centroid of the input two-dimentional image data @@ -624,8 +885,23 @@ def centroidN(data): #################################################################### - 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; @@ -640,11 +916,26 @@ def SpecCube2photon(data, wave): return Nphoton ###################################################################### +###################################################################### +def mag2flux(starmag, lam): + """ -#################################################################################### + Parameters + ---------- + starmag : TYPE + DESCRIPTION. + lam : TYPE + DESCRIPTION. + + Returns + ------- + fluxlam : TYPE + DESCRIPTION. + photonu : TYPE + DESCRIPTION. -def mag2flux(starmag, lam): + """ # starmag: input AB mag # lambda: wavelength in nm @@ -662,8 +953,23 @@ def mag2flux(starmag, lam): ############################################################################### - def flux2photons(flux, lam): + """ + + + Parameters + ---------- + flux : TYPE + DESCRIPTION. + lam : TYPE + DESCRIPTION. + + Returns + ------- + photonu : TYPE + DESCRIPTION. + + """ # flux: # erg/s/cm2/A/ # lam: # wavelength in nm @@ -674,17 +980,48 @@ def flux2photons(flux, lam): return photonu ############################################################################### -################################################## +def fftrange(n, dtype=float): + """ -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): + """ -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 @@ -723,59 +1060,35 @@ def dft2(ary, Q, samples, shift=None): 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): + """ -def idft2(ary, Q, samples, shift=None): - """Compute the two dimensional inverse 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. + wavefront : TYPE + DESCRIPTION. + pupil : TYPE + DESCRIPTION. + Q : TYPE + DESCRIPTION. + sizeout : TYPE + DESCRIPTION. Returns ------- - `numpy.ndarray` - 2D array containing the shifted transform. - Equivalent to ifftshift(ifft2(fftshift(ary))) modulo output - sampling/grid differences + psf : TYPE + DESCRIPTION. """ - # 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 - X, Y, U, V = (fftrange(n) for n in (m, n, M, N)) - ############################################################################### - nm = n*m - NM = N*M - r = NM/nm - a = 1 / Q - ################################################################### - Eout_rev = np.exp(1j * 2 * np.pi * a / n * np.outer(Y, V).T) * (1/r) - Ein_rev = np.exp(1j * 2 * np.pi * a / m * np.outer(X, U)) * (1/nm) - ########################################################################### - out = Eout_rev @ ary @ Ein_rev - return out -############################################################################### - - -def anySampledPSFnew(wavefront, pupil, Q, sizeout): - ''' + ''' written on 2020.12.04 by Yan @Shao % wavefront sampled in the united circle; % pupil function samped as wavefront; @@ -804,7 +1117,7 @@ def anySampledPSFnew(wavefront, pupil, Q, sizeout): cx, cy = centroid(psfout) Nt = sizeout - # psf=ndimage.shift(psfout,[Nt/2-cy, Nt/2-cx],order=1, mode='nearest' ) ## for my fft method + # for convolve method psf = ndimage.shift( psfout, [Nt/2-cy-1, Nt/2-cx-1], order=1, mode='nearest') @@ -812,79 +1125,60 @@ def anySampledPSFnew(wavefront, pupil, Q, sizeout): return psf -########################################################################################### +############################################################################## +############################################################################## -def zero_pad(image, shape, position='corner'): +def get_dx_dy_blue(wave): """ - Extends image to a certain size with zeros + + Parameters ---------- - image: real 2d `np.ndarray` - Input image - shape: tuple of int - Desired output shape of the image - position : str, optional - The position of the input image in the output one: - * 'corner' - top-left corner (default) - * 'center' - centered + wave : TYPE + DESCRIPTION. + Returns ------- - padded_img: real `np.ndarray` - The zero-padded image - """ - shape = np.asarray(shape, dtype=int) - imshape = np.asarray(image.shape, dtype=int) - - if np.alltrue(imshape == shape): - return image - - if np.any(shape <= 0): - raise ValueError("ZERO_PAD: null or negative shape given") - - dshape = shape - imshape - if np.any(dshape < 0): - raise ValueError("ZERO_PAD: target size smaller than source one") - - pad_img = np.zeros(shape, dtype=image.dtype) - - idx, idy = np.indices(imshape) - - if position == 'center': - if np.any(dshape % 2 != 0): - raise ValueError("ZERO_PAD: source and target shapes " - "have different parity.") - offx, offy = dshape // 2 - else: - offx, offy = (0, 0) - - pad_img[idx + offx, idy + offy] = image - - return pad_img -################################################################################ - -############################################################################## + dx : TYPE + DESCRIPTION. + dy : TYPE + DESCRIPTION. - -def get_dx_dy_blue(wave): + """ # 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]) # 垂直方向 - + 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]) + # 垂直方向 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]) # 色散方向 @@ -900,8 +1194,25 @@ def get_dx_dy_red(wave): ############################################################################################## - def getSpectrum(Spectrum0, lam, sigma): + """ + + + Parameters + ---------- + Spectrum0 : TYPE + DESCRIPTION. + lam : TYPE + DESCRIPTION. + sigma : TYPE + DESCRIPTION. + + Returns + ------- + SpmOut : TYPE + DESCRIPTION. + + """ # %获得输入波长lambda的光谱强度 wave = Spectrum0[:, 0] # %原始光谱数据给定的波长； @@ -918,39 +1229,8 @@ def getSpectrum(Spectrum0, lam, sigma): return SpmOut -##################################################################################### - - -#################################################################################################################### - -def processArgs(printHelp=False): - """ - Processes command line arguments. - """ - parser = OptionParser() - - parser.add_option('-c', '--configfile', dest='configfile', - help="Name of the configuration file", metavar="string") - parser.add_option('-s', '--section', dest='section', - help="Name of the section of the config file [SCIENCE]", metavar="string") - parser.add_option('-q', '--quadrant', dest='quadrant', help='CCD quadrant to simulate [0, 1, 2, 3]', - metavar='int') - parser.add_option('-x', '--xCCD', dest='xCCD', help='CCD number in X-direction within the FPA matrix', - metavar='int') - parser.add_option('-y', '--yCCD', dest='yCCD', help='CCD number in Y-direction within the FPA matrix', - metavar='int') - parser.add_option('-d', '--debug', dest='debug', action='store_true', - help='Debugging mode on') - parser.add_option('-t', '--test', dest='test', action='store_true', - help='Run unittest') - parser.add_option('-f', '--fixed', dest='fixed', help='Use a fixed seed for the random number generators', - metavar='int') - if printHelp: - parser.print_help() - else: - return parser.parse_args() -############################################################################################## - +############################################################################## +############################################################################## class IFSsimulator(): """ @@ -978,27 +1258,9 @@ class IFSsimulator(): #################################### self.configfile = configfile + self.section = 'TEST' # simulation section; - # if opts.section is None: - # ####self.section = 'DEFAULT' - # self.section = 'TEST' # simulation section; - - # else: - # self.section = opts.section - - self.section = 'TEST' # simulation section; - - # if opts.debug is None: - # self.debug = False - # else: - # self.debug = opts.debug - - # try: - # self.random = opts.testing - # except: - # self.random = False - - # load instrument model, these values are also stored in the FITS header + # load instrument model self.information = IFSinstrumentModel.IFSinformation() # update settings with defaults @@ -1026,20 +1288,29 @@ class IFSsimulator(): ra=0.0, dec=0.0, injection=80000.0, - coveringfraction=0.5, # CR: 3.0% is for 565s exposure - ################################################### - # cosmicraylengths=self.information['dir_path']+'IFS_inputdata/cdf_cr_length.dat', - # cosmicraydistance=self.information['dir_path']+'IFS_inputdata/cdf_cr_total.dat', - # parallelTrapfile=self.information['dir_path']+'IFS_inputdata/cdm_euclid_parallel.dat', - # serialTrapfile=self.information['dir_path']+'IFS_inputdata/cdm_euclid_serial.dat', - # cosmeticsFile='../IFS_inputdata/cosmetics.dat', + coveringfraction=0.5, + # CR: 3.0% is for 565s exposure + ######################################### + mode='same')) ##################################### -# 3 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() @@ -1048,9 +1319,9 @@ class IFSsimulator(): # self.log.info(self.information) - ############################################################################### + def processConfigs(self): """ Processes configuration information and save the information to a dictionary self.information. @@ -1083,15 +1354,6 @@ class IFSsimulator(): # force gain to be float self.information['e_adu'] = float(self.information['e_adu']) - # #name of the output file, include quadrants and CCDs - ################################################################################################### - #self.information['output'] = self.config.get(self.section, 'output') - ################################################################################################## - - # #booleans to control the flow - # self.chargeInjectionx = self.config.getboolean(self.section, 'chargeInjectionx') - # self.chargeInjectiony = self.config.getboolean(self.section, 'chargeInjectiony') - self.cosmicRays = self.config.getboolean(self.section, 'cosmicRays') self.darknoise = self.config.getboolean(self.section, 'darknoise') self.cosmetics = self.config.getboolean(self.section, 'cosmetics') @@ -1126,7 +1388,7 @@ class IFSsimulator(): except: self.intscale = True - ############################################################################## + ###################################################################### # 3 @@ -1146,7 +1408,7 @@ class IFSsimulator(): ##################################################################### - ######################################################################################### + ############################################################################ def _createEmpty(self): """ @@ -1163,12 +1425,32 @@ class IFSsimulator(): self.pixel = 0.1 # arcsec, pixel scale size; -######################################################################################################### -########################################################## - +############################################################################## +############################################################################## def zodiacal(self, ra, dec, time): """ - For given RA, DEC and TIME, return the interpolated zodical spectrum in Leinert-1998. + + + 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 @@ -1183,7 +1465,7 @@ class IFSsimulator(): # get solar position dt = datetime.fromisoformat(time) jd = time2jd(dt) - ##jd = julian.to_jd(dt, fmt='jd') + ## jd = julian.to_jd(dt, fmt='jd') t = Time(jd, format='jd', scale='utc') @@ -1236,10 +1518,26 @@ class IFSsimulator(): # 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. @@ -1258,21 +1556,28 @@ class IFSsimulator(): :rtype: ndarray """ return ndimage.filters.gaussian_filter(image, sigma) -################################################################################# +############################################################################### + + def readCosmicRayInformation(self): + """ + + + Returns + ------- + None. - def readCosmicRayInformation(self): + """ """ 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'])) + 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']) + self.information['dir_path']+self.information['cosmicraylengths']) crDists = np.loadtxt( self.information['dir_path']+self.information['cosmicraydistance']) @@ -1280,92 +1585,106 @@ class IFSsimulator(): cr_v=crDists[:, 0], cr_cde=crDists[:, 1], cr_cden=np.shape(crDists)[0]) ############################################################################### - def _writeFITSfile(self, image, filename): - """ - :param image: image array to save - :type image: ndarray - :param filename: name of the output file, e.g. file.fits - :type filename: str + # def _writeFITSfile(self, image, filename): + # """ + # :param image: image array to save + # :type image: ndarray + # :param filename: name of the output file, e.g. file.fits + # :type filename: str - :return: None - """ - if os.path.isfile(filename): - os.remove(filename) + # :return: None + # """ + # if os.path.isfile(filename): + # os.remove(filename) - # create a new FITS file, using HDUList instance - ofd = fits.HDUList(fits.PrimaryHDU()) + # # create a new FITS file, using HDUList instance + # ofd = fits.HDUList(fits.PrimaryHDU()) - # new image HDU - hdu = fits.ImageHDU(data=image) + # # new image HDU + # hdu = fits.ImageHDU(data=image) - # update and verify the header - hdu.header.add_history('Created by IFSsim at %s' % - datetime.isoformat(datetime.utcnow())) - hdu.verify('fix') + # # update and verify the header + # hdu.header.add_history('Created by IFSsim at %s' % + # datetime.isoformat(datetime.utcnow())) + # hdu.verify('fix') - ofd.append(hdu) + # ofd.append(hdu) - # write the actual file - ofd.writeto(filename) + # # write the actual file + # ofd.writeto(filename) ############################################################################### - def writeFITSfile(self, data, filename, unsigned16bit=False): - """ - Writes out a simple FITS file. - - :param data: data to be written - :type data: ndarray - :param filename: name of the output file - :type filename: str - :param unsigned16bit: whether to scale the data using bzero=32768 - :type unsigned16bit: bool + # def writeFITSfile(self, data, filename, unsigned16bit=False): + # """ + # Writes out a simple FITS file. + + # :param data: data to be written + # :type data: ndarray + # :param filename: name of the output file + # :type filename: str + # :param unsigned16bit: whether to scale the data using bzero=32768 + # :type unsigned16bit: bool + + # :return: None + # """ + # if os.path.isfile(filename): + # os.remove(filename) + + # # create a new FITS file, using HDUList instance + # ofd = fits.HDUList(fits.PrimaryHDU()) + + # # new image HDU + # hdu = fits.ImageHDU(data=data) + + # # add input keywords to the header + # for key, value in self.information.items(): + # # truncate long keys + # if len(key) > 8: + # key = key[:7] + # try: + # hdu.header.set(key.upper(), value) + # except: + # try: + # hdu.header.set(key.upper(), str(value)) + # except: + # pass + + # # write booleans + # for key, value in self.booleans.items(): + # # truncate long keys + # if len(key) > 8: + # key = key[:7] + # hdu.header.set(key.upper(), str(value), 'Boolean Flags') + + # # update and verify the header + # hdu.header.add_history( + # 'This is an itermediate data product no the final output!') + + # hdu.verify('fix') + + # ofd.append(hdu) + + # # write the actual file + # ofd.writeto(filename) - :return: None +############################################################################### + def configure(self, simnumber): """ - if os.path.isfile(filename): - os.remove(filename) - - # create a new FITS file, using HDUList instance - ofd = fits.HDUList(fits.PrimaryHDU()) - - # new image HDU - hdu = fits.ImageHDU(data=data) - - # add input keywords to the header - for key, value in self.information.items(): - # truncate long keys - if len(key) > 8: - key = key[:7] - try: - hdu.header.set(key.upper(), value) - except: - try: - hdu.header.set(key.upper(), str(value)) - except: - pass - # write booleans - for key, value in self.booleans.items(): - # truncate long keys - if len(key) > 8: - key = key[:7] - hdu.header.set(key.upper(), str(value), 'Boolean Flags') - - # update and verify the header - hdu.header.add_history( - 'This is an itermediate data product no the final output!') - hdu.verify('fix') + Parameters + ---------- + simnumber : TYPE + DESCRIPTION. - ofd.append(hdu) - - # write the actual file - ofd.writeto(filename) + Returns + ------- + None. -############################################################################### - def configure(self, simnumber): """ - Configures the simulator with input information and creates and empty array to which the final image will + """ + Configures the simulator with input information and creates and + empty array to which the final image will be build on. """ self.readConfigs(simnumber) @@ -1391,8 +1710,9 @@ class IFSsimulator(): else: ss = '_' - if self.information['dir_path']=='/nfsdata/share/simulation-unittest/ifs_sim/': - self.result_path = self.information['dir_path']+'ifs_sim_result/'+self.source+ss+result_day + if self.information['dir_path'] == '/nfsdata/share/simulation-unittest/ifs_sim/': + self.result_path = self.information['dir_path'] + \ + 'ifs_sim_result/'+self.source+ss+result_day else: home_path = os.environ['HOME'] @@ -1402,9 +1722,6 @@ class IFSsimulator(): self.result_path = '../IFS_simData_'+self.source+ss+result_day else: self.result_path = '/data/ifspip/CCD_ima/IFS_simData_'+self.source+ss+result_day - - - if os.path.isdir(self.result_path) == False: os.mkdir(self.result_path) @@ -1526,7 +1843,7 @@ class IFSsimulator(): self.slice_blue = slice_blue self.slice_red = slice_red - ############################################################################### + ####################################################################### maskSlice = dict() maskSlit = dict() sizeout = 100 @@ -1541,19 +1858,37 @@ class IFSsimulator(): self.maskSlice = maskSlice self.maskSlit = maskSlit - ################################################################################################ +###################################################################### -################################################################################################################# +############################################################################## def generateflat(self, ave=1.0, sigma=0.01): """ + + + Parameters + ---------- + ave : TYPE, optional + DESCRIPTION. The default is 1.0. + sigma : TYPE, optional + DESCRIPTION. The default is 0.01. + + Returns + ------- + TYPE + DESCRIPTION. + TYPE + DESCRIPTION. + + """ + """ Creates a flat field image with given properties. :return: flat field image :rtype: ndarray """ self.log.info('Generating a flat field...') - self.log.info( - 'The flat field has mean value of 1 and a given fluctuations, usually either 1 or 2 percent defined by sigma= %d...' % sigma) + self.log.info('The flat field has mean value of 1 and a given fluctuations, + usually either 1 or 2 percent defined by sigma= %d...' % sigma) np.random.seed(5*self.simnumber) self.flat_b = np.random.normal(loc=ave, scale=sigma, size=(2048, 4096)) @@ -1589,10 +1924,18 @@ class IFSsimulator(): return self.flat_b, self.flat_r - ######################################################################################################### + ########################################################################## def addLampFlux(self): """ + + + Returns + ------- + None. + + """ + """ Include flux from the calibration source. """ @@ -1604,6 +1947,22 @@ class IFSsimulator(): ############################################################################# def MakeFlatMatrix(self, img, seed): + """ + + + Parameters + ---------- + img : TYPE + DESCRIPTION. + seed : TYPE + DESCRIPTION. + + Returns + ------- + FlatMat : TYPE + DESCRIPTION. + + """ #### ysize, xsize = img.shape np.random.seed(seed) @@ -1625,6 +1984,14 @@ class IFSsimulator(): def getFrameTransferImg(self): """ + + + Returns + ------- + None. + + """ + """ get frame transfer image from original image of self.image_b and self.image_r """ self.frame_b_4 = np.zeros((1024+320, 2048), dtype=float) @@ -1670,11 +2037,24 @@ class IFSsimulator(): ############################################################################### 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. + the sky. """ # apply Flat field to image ### @@ -1690,31 +2070,30 @@ class IFSsimulator(): return -################################################################################# +############################################################################### - def addChargeInjection(self): - """ - Add either horizontal or vertical charge injection line to the image. +############################################################################### + def addCosmicRays(self, idk): """ - if self.chargeInjectionx: - #self.image[1500:1511, :] = self.information['injection'] - self.image_b[1500:1511, :] = self.information['injection'] - self.image_r[1500:1511, :] = self.information['injection'] + - self.log.info('Adding vertical charge injection line') + Parameters + ---------- + idk : TYPE + DESCRIPTION. - if self.chargeInjectiony: - #self.image[:, 1500:1511] = self.information['injection'] - self.image_b[:, 1500:1511] = self.information['injection'] - self.image_r[:, 1500:1511] = self.information['injection'] - self.log.info('Adding horizontal charge injection line') + Returns + ------- + None. -############################################################################### - def addCosmicRays(self, idk): """ - 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. + """ + 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 @@ -1761,6 +2140,14 @@ class IFSsimulator(): def addReadoutTrails(self): """ + + + Returns + ------- + None. + + """ + """ Add readout trails resulting from reading out the shutter open. Quadrants assumed to be numbered: @@ -1791,7 +2178,7 @@ class IFSsimulator(): else: self.image_b = data - # 3 + flux_ratio = self.information['readouttime'] / float( self.information['redsize']) / self.information['exptime'] # make a copy, this will be updated @@ -1815,10 +2202,18 @@ class IFSsimulator(): else: self.image_r = data -############################################################################################################## +############################################################################## 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. @@ -1859,11 +2254,18 @@ class IFSsimulator(): self.information['dark2_r'] -# 3 - +############################################################################## def applyCosmicBackground(self): """ + + + Returns + ------- + None. + + """ + """ Apply dark the cosmic background. Scales the background with the exposure time. Additionally saves the image without noise to a FITS file. @@ -1882,10 +2284,18 @@ class IFSsimulator(): self.imagenoCR_b += bcgr self.imagenoCR_r += bcgr - ########################################################################################## + ########################################################################## def applyScatteredLight(self): """ + + + Returns + ------- + None. + + """ + """ Adds spatially uniform scattered light to the image. """ sl = self.information['exptime'] * self.information['scattered_light'] @@ -1898,6 +2308,14 @@ class IFSsimulator(): ############################################################################## def applyPoissonNoise(self): """ + + + Returns + ------- + None. + + """ + """ Add Poisson noise to the image. """ @@ -1923,9 +2341,16 @@ class IFSsimulator(): ################################################################################################################### - 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 @@ -1952,7 +2377,7 @@ class IFSsimulator(): 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']) @@ -1991,6 +2416,14 @@ class IFSsimulator(): def applyRadiationDamage(self): """ + + + Returns + ------- + None. + + """ + """ Applies CDM03 radiation model to the image being constructed. .. seealso:: Class :`CDM03` @@ -2011,10 +2444,18 @@ class IFSsimulator(): 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. """ @@ -2029,13 +2470,22 @@ class IFSsimulator(): 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. + 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') # blue zone 1 @@ -2061,7 +2511,6 @@ class IFSsimulator(): self.image_b[0:1024+overscan, 2418*3:2418*3+2048+prescan+overscan] += np.random.normal( loc=0.0, scale=self.information['rn2_b'], size=(1344, 2418)) - # 3 self.log.info('readnoise added in blue channel') @@ -2089,6 +2538,14 @@ class IFSsimulator(): def appFrameTransferEffect(self): """ + + + Returns + ------- + None. + + """ + """ apply frame transfer effect to the four part images """ prescan = 50 @@ -2136,10 +2593,18 @@ class IFSsimulator(): self.image_r[0:1536+overscan, prescan+3442*3:3442*3 + 3072+prescan] += self.frame_r_2[0:1536+overscan, 0:3072] - # 3 + ########################################################################## def electrons2ADU(self): """ + + + Returns + ------- + None. + + """ + """ Convert from electrons to ADUs using the value read from the configuration file. """ self.log.info( @@ -2177,6 +2642,14 @@ class IFSsimulator(): 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 @@ -2199,7 +2672,7 @@ class IFSsimulator(): self.log.info('Bias counts were added to the blue image') - ############################################################################ + ####################################################################### # red zone 4 self.image_r[0:1856, 0:3442] += self.information['bias4_r'] @@ -2213,14 +2686,21 @@ class IFSsimulator(): # zone 2 self.image_r[0:1856, 3442*3:3442*4] += self.information['bias2_r'] - ########################################################################## + ####################################################################### self.log.info('Bias counts were added to the red image') - # 3 - +############################################################################### 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. @@ -2265,13 +2745,23 @@ class IFSsimulator(): self.log.error( 'Cannot include pre- and overscan because of an unknown quadrant!') -#################################################################################### +############################################################################### def applyBleeding_yan(self): """ - Apply bleeding along the CCD columns if the number of electrons in a pixel exceeds the full-well capacity. + + + 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 + Bleeding is modelled in the parallel direction only, because the + CCD273s are assumed not to bleed in serial direction. :return: None @@ -2355,12 +2845,25 @@ class IFSsimulator(): 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. @@ -2390,6 +2893,14 @@ class IFSsimulator(): ################################################################################################## def applyImageShift(self): + """ + + + Returns + ------- + None. + + """ np.random.seed(9*self.simnumber) ud = np.random.random() # Choose a random rotation @@ -2409,10 +2920,17 @@ class IFSsimulator(): self.information['dec'] = dy*self.information['pixel_size'] -# 33 - +############################################################################## 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'] @@ -2436,6 +2954,14 @@ class IFSsimulator(): ############################################################################### def CCDreadout(self): + """ + + + Returns + ------- + None. + + """ imgb = self.image_b.copy() prescan = int(self.information['prescan']) @@ -2449,7 +2975,7 @@ class IFSsimulator(): y1 = 0+prescan y2 = y1+2048 temp[x1:x2, y1:y2] = imgb[0:1024, 0:2048] - ########## zone 3, OSG , left to right ################# + # zone 3, OSG , left to right ################# # np.fliplr(b2) ## left to right # np.flipud(b3) ## down to up x1 = 0 @@ -2458,7 +2984,7 @@ class IFSsimulator(): y1 = 2418+prescan y2 = y1+2048 temp[x1:x2, y1:y2] = np.fliplr(imgb[0:1024, 2048:4096]) - ########## zone 1, OSE,down to up ################### + # zone 1, OSE,down to up ################### x1 = 0 x2 = x1+1024 @@ -2475,7 +3001,7 @@ class IFSsimulator(): self.image_b = temp - ############################################################################## + ####################################################################### imgr = self.image_r.copy() temp = np.zeros((1856, 13768)) @@ -2517,6 +3043,19 @@ class IFSsimulator(): 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. @@ -2760,8 +3299,6 @@ class IFSsimulator(): ########## ########## finish header for 0 layer ###################### ############################################################# - - # 3 # header # new image HDU, blue channel, layer 1 if HeaderTest == 'yes': @@ -2794,8 +3331,8 @@ class IFSsimulator(): hdu_b.header['BZERO'] = (np.float64(32768), '') hdu_b.header['BUNIT'] = ('ADU', 'physical unit of array values') - ########################################################################### - ######### instrument information ######################################## + ####################################################################### + ######### instrument information ################################### if self.source == 'SCI' or self.source == 'COMP': hdu_b.header['CMIRRPOS'] = ( @@ -2824,9 +3361,9 @@ class IFSsimulator(): hdu_b.header['IFSSTAT'] = ( np.int32(0), 'IFS components status parameter') - ############################################################################## + ###################################################################### - #################### detector information############################# + #################### detector information########################## hdu_b.header['CAMERA'] = ('Blue', 'camera of IFS') hdu_b.header['DETSN'] = ('CCD231-c4-00', 'detector serial number') @@ -2888,7 +3425,7 @@ class IFSsimulator(): hdu_b.header['DETTEMP1'] = (np.float32( 0.0), 'detector temperature at EXPT1_'+str(k+1)+' (K)') - ###################################################### end revised on 2024.2.27 ### + #######################end revised on 2024.2.27 ### hdu_b.header['BIN_X'] = ( np.int16(1), 'bin number in X (wavelength)') @@ -2944,7 +3481,7 @@ class IFSsimulator(): 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 @@ -2980,7 +3517,7 @@ class IFSsimulator(): ofd_r.header['OBJECT'] = ( self.information['name_obj'][:30], 'object name') ofd_r.header['TARGET'] = ( - (self.information['target']), 'target name (hhmmss.s+ddmmss)') + (self.information['target']), 'target name (hhmmss.s+ddmmss)') ofd_r.header['OBSID'] = (str(obsid), 'observation ID') @@ -3142,7 +3679,6 @@ class IFSsimulator(): ########## finish header for 0 layer ###################### ############################################################# - # 3 # header # new image HDU, red channel, layer 1 if HeaderTest == 'yes': @@ -3175,7 +3711,7 @@ class IFSsimulator(): hdu_r.header['BUNIT'] = ('ADU', 'physical unit of array values') ######### instrument information ###### - ##hdu_r.header['PMIRRPOS']=(bool(False), 'FSM pointing,T: to MCI, F: not to MCI') + if self.source == 'SCI' or self.source == 'COMP': hdu_r.header['CMIRRPOS'] = ( @@ -3204,7 +3740,7 @@ class IFSsimulator(): hdu_r.header['IFSSTAT'] = ( np.int32(0), 'IFS components status parameter') - ################### detector information############################# + ################### detector information############################ hdu_r.header['CAMERA'] = ('Red', 'camera of IFS') hdu_r.header['DETSN'] = ('CCD231-c6-00', 'detector serial number') @@ -3231,7 +3767,7 @@ class IFSsimulator(): 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 @@ -3343,10 +3879,7 @@ class IFSsimulator(): hdu1.header.add_comment( '========================================================================', after='EPOCH') - ######################################################################################################################### - hdu2.header.add_comment( - '========================================================================', after='FITSSWV') - hdu2.header.add_comment('OBJECT INFORMATION', after='FITSSWV') + ###################################################################### hdu2.header.add_comment( '========================================================================', after='FITSSWV') @@ -3383,7 +3916,7 @@ class IFSsimulator(): '========================================================================', 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( @@ -3406,7 +3939,7 @@ class IFSsimulator(): hdulist_r.writeto(self.file_r, output_verify='ignore', checksum=True) - ################################################################################## + ###################################################################### def earthshine(self, theta): """ @@ -3472,7 +4005,22 @@ class IFSsimulator(): ################################################################################### ################################################################################## + def CalskyNoise(self, lam): + """ + + + Parameters + ---------- + lam : TYPE + DESCRIPTION. + + Returns + ------- + Ns_skynoise : TYPE + DESCRIPTION. + + """ # calculate the reference flux # calculate sky noise; @@ -3505,6 +4053,19 @@ class IFSsimulator(): ############## 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 ################### @@ -3722,14 +4283,14 @@ class IFSsimulator(): img3 = conv - ######################## get subimage ##################### + ######################## 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 ##################### + ####### get photons from sub-image ##################### subimage.scale = 0.01 # here pixelscale=0.01arcsec subimage.setOrigin(0, 0) @@ -3737,14 +4298,14 @@ class IFSsimulator(): 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 @@ -3832,6 +4393,33 @@ class IFSsimulator(): ################################################################################## 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 @@ -4108,8 +4696,6 @@ class IFSsimulator(): 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) @@ -4394,6 +4980,21 @@ class IFSsimulator(): ################################################################################################# def sim_calibration_img(self, exposuretime, source): + """ + + + Parameters + ---------- + exposuretime : TYPE + DESCRIPTION. + source : TYPE + DESCRIPTION. + + Returns + ------- + None. + + """ sizeout = 100 ################# load the hole mask file ################### @@ -4522,7 +5123,7 @@ class IFSsimulator(): slice_image = ns*self.telarea*Width_lambda*10*exptime -################################################################################################## +############################################################################ if slice_image.max() < 0.1: continue @@ -4530,7 +5131,7 @@ class IFSsimulator(): 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 ; @@ -4735,7 +5336,7 @@ class IFSsimulator(): hdu2 = fits.ImageHDU(Wave) - hdu2.header.append(('WaveUnit', 'nm')) + hdu2.header.append(('WaveUnit', 'nm')) newd = fits.HDUList([hdu1, hdu2]) @@ -4758,6 +5359,21 @@ class IFSsimulator(): ################################################################################################# 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. @@ -4789,9 +5405,10 @@ class IFSsimulator(): if simnumber <= 200 and simnumber > 150: self.information['exptime'] = 1200 - self.skyfilepath = self.information['dir_path']+self.information['sky_fitsin'] + self.skyfilepath = self.information['dir_path'] + \ + self.information['sky_fitsin'] - print('self.skyfilepath = ',self.skyfilepath) + print('self.skyfilepath = ', self.skyfilepath) # self.earthshine_theta=30.0 # in degree @@ -5030,7 +5647,29 @@ class IFSsimulator(): ############################################################################ ############################################################################ + + def runIFSsim(sourcein, configfile, iLoop, applyhole='no'): + """ + + + 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)