diff --git a/csst_ifs_sim/csst_ifs_sim.py b/csst_ifs_sim/csst_ifs_sim.py index d809ab871d4c4a650392953ead3889746c16be7b..f59d022cde448cc7c5a59ce2040b9b974a482b65 100644 --- a/csst_ifs_sim/csst_ifs_sim.py +++ b/csst_ifs_sim/csst_ifs_sim.py @@ -5,6 +5,27 @@ 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 +from support import IFSinstrumentModel +from support import cosmicrays +from support import logger as lg +from CTI import CTI import os from scipy.interpolate import interp1d import astropy.coordinates as coord @@ -15,27 +36,6 @@ import sys # from csst_ifs_sim.support import logger as lg # from csst_ifs_sim.CTI import CTI sys.path.append('./csst_ifs_sim') -from CTI import CTI -from support import logger as lg -from support import cosmicrays -from support import IFSinstrumentModel -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 -from astropy.utils.iers import conf conf.auto_max_age = None """ @@ -56,10 +56,7 @@ The approximate sequence of events in the simulator is as follows: #. 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 with defined - wavelength and field of view. - -#. Apply calibration unit flux to mimic flat field exposures [optional]. +#. 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 @@ -67,12 +64,11 @@ distribution. #. Apply detector charge bleeding in column direction. #. Add constant dark current and background light from Zodiacal light. -#. Include spatially uniform scattered light to the pixel grid. #. 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 [optional]. -#. Apply CCD273 non-linearity model to the pixel data. +#. 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 @@ -82,17 +78,17 @@ 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 +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 ##### + def transRaDec2D(ra, dec): """ @@ -116,52 +112,6 @@ def transRaDec2D(ra, dec): return np.array([x1, y1, z1]) ############################################################################### - -# 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 - -# # 1 W/m^2/sr/μm = 0.10 erg/cm^2/s/sr/A -# # 1 sr = 1 rad^2 = 4.25452e10 arcsec^2 -# # 1 J/s = 1 W -# # 1 J = 10^7 erg - -# # convert erg/s/cm^2/A/arcsec^2 to erg/s/cm^2/A/sr -# flux1 = flux / (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 -# flux3 = flux2 * np.pi * D**2 / 4 / f**2 / 10**3 - -# # 对波长积分 -# f = interp1d(wave, flux3) -# wave_interp = np.arange(3800, 7800) -# flux3_interp = f(wave_interp) -# # 输出单位 W/m^2 -# delta_lamba = 0.1 # nm -# E = np.sum(flux3_interp * delta_lamba) - -# return E - ################################################################ @@ -214,7 +164,8 @@ def ill2flux(path, E): ############################################################## -########################################################################################## +############################################################################## + def earth_angle(time_jd, x_sat, y_sat, z_sat, ra_obj, dec_obj): """ @@ -268,6 +219,7 @@ def earth_angle(time_jd, x_sat, y_sat, z_sat, ra_obj, dec_obj): ##################################################################### + class StrayLight(object): def __init__(self, path, jtime=2460843., sat=np.array([0, 0, 0]), radec=np.array([0, 0])): @@ -311,6 +263,7 @@ 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'): """ @@ -326,7 +279,8 @@ class StrayLight(object): DESCRIPTION. """ - filterIndex = {'nuv': 0, 'u': 1, 'g': 2, 'r': 3, 'i': 4, 'z': 5, 'y': 6} + filterIndex = {'nuv': 0, 'u': 1, 'g': 2, + 'r': 3, 'i': 4, 'z': 5, 'y': 6} sat = (ctypes.c_double*3)() sat[:] = self.sat @@ -349,6 +303,7 @@ class StrayLight(object): return max(band_star_e1, band_star_e2) ############################################################################### + def caculateEarthShineFilter(self, filter='i'): """ # @@ -364,7 +319,8 @@ class StrayLight(object): DESCRIPTION. """ - filterIndex = {'nuv': 0, 'u': 1, 'g': 2, 'r': 3, 'i': 4, 'z': 5, 'y': 6} + filterIndex = {'nuv': 0, 'u': 1, 'g': 2, + 'r': 3, 'i': 4, 'z': 5, 'y': 6} sat = (ctypes.c_double*3)() sat[:] = self.sat @@ -377,7 +333,7 @@ class StrayLight(object): earth_e1 = (ctypes.c_double*7)() self.slcdll.EarthShine(self.jtime, sat, ob, py1, earth_e1) - # e[7]代表7个波段的照度 + # e[7]代表7个波段的照度 earth_e2 = (ctypes.c_double*7)() self.slcdll.EarthShine(self.jtime, sat, ob, py2, earth_e2) @@ -391,29 +347,6 @@ 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): """ @@ -436,6 +369,8 @@ def time2mjd(dateT): return mjd+mjd_s/86400.0 ########################################################################### + + def time2jd(dateT): """ @@ -460,6 +395,7 @@ def time2jd(dateT): # mjd转datetime类 ############################################################################# + def mjd2time(mjd): """ @@ -479,6 +415,8 @@ def mjd2time(mjd): return t0+timedelta(days=mjd) ############################################################################# + + def jd2time(jd): """ @@ -500,6 +438,8 @@ def jd2time(jd): # mjd和jd互转 ############################################################################ + + def mjd2jd(mjd): """ @@ -518,6 +458,8 @@ def mjd2jd(mjd): return mjd+2400000.5 ########################################################################### + + def jd2mjd(jd): """ @@ -536,6 +478,8 @@ def jd2mjd(jd): return jd-2400000.5 ########################################################################## + + def dt2hmd(dt): """ @@ -573,36 +517,9 @@ def dt2hmd(dt): return str_h+':'+str_m+':'+str_d ############################################################################### -# 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) -# n = len(b) - -# if z >= n: -# return b -# else: -# return b[:z] - - ############################################################################### + def deg2HMS(ra0, dec0): """ @@ -715,6 +632,8 @@ def beta_angle(x_sat, y_sat, z_sat, vx_sat, vy_sat, vz_sat, ra_obj, dec_obj): ############################################################################### ########################################################## + + def LSR_velocity(ra, dec, velocity, Obstime): """ @@ -822,42 +741,7 @@ def centroid(data): return float(cx), float(cy) ############################################################################### -# def centroidN(data): -# """ - - -# Parameters -# ---------- -# data : TYPE -# DESCRIPTION. - -# Returns -# ------- -# cx : TYPE -# DESCRIPTION. -# cy : TYPE -# DESCRIPTION. - -# """ -# ''' -# calculate the centroid of the input two-dimentional image data - -# Parameters -# ---------- -# data : input image. - -# Returns -# ------- -# cx: the centroid column number, in horizontal direction definet in python image show -# cy: the centroid row number , in vertical direction - -# ''' -# ### -# from scipy import ndimage -# cy, cx = ndimage.center_of_mass(data) -# return cx, cy -#################################################################### def SpecCube2photon(data, wave): """ @@ -891,6 +775,8 @@ def SpecCube2photon(data, wave): ###################################################################### ###################################################################### + + def mag2flux(starmag, lam): """ @@ -927,6 +813,7 @@ def mag2flux(starmag, lam): ############################################################################### + def flux2photons(flux, lam): """ @@ -954,6 +841,8 @@ def flux2photons(flux, lam): return photonu ############################################################################### + + def fftrange(n, dtype=float): """ @@ -975,6 +864,8 @@ def fftrange(n, dtype=float): return np.arange(-n//2, -n//2+n, dtype=dtype) ############################################################################### + + def dft2(ary, Q, samples, shift=None): """ @@ -1041,6 +932,8 @@ def dft2(ary, Q, samples, shift=None): return out ############################################################################## ############################################################################### + + def anySampledPSFnew(wavefront, pupil, Q, sizeout): """ @@ -1091,7 +984,7 @@ def anySampledPSFnew(wavefront, pupil, Q, sizeout): 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') @@ -1103,6 +996,7 @@ def anySampledPSFnew(wavefront, pupil, Q, sizeout): ############################################################################## + def get_dx_dy_blue(wave): """ @@ -1136,6 +1030,7 @@ def get_dx_dy_blue(wave): return dx, dy ############################################################################## + def get_dx_dy_red(wave): """ @@ -1168,6 +1063,7 @@ def get_dx_dy_red(wave): ############################################################################################## + def getSpectrum(Spectrum0, lam, sigma): """ @@ -1241,7 +1137,7 @@ class IFSsimulator(): self.information.update(dict(quadrant=int(0), ccdx=int(0), ccdy=int(0), - + ccdxgap=1.643, ccdygap=8.116, bluesize=4000, @@ -1286,11 +1182,10 @@ class IFSsimulator(): self.config.read_file(open(self.configfile)) - # self.log.info(self.information) + self.log.info(self.information) ############################################################################### - def processConfigs(self): """ Processes configuration information and save the information to a dictionary self.information. @@ -1334,30 +1229,18 @@ class IFSsimulator(): 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.intscale = True - ###self.config.getboolean(self.section, 'intscale') - - + self.intscale = True ###################################################################### - # 3 - self.information['variablePSF'] = False self.booleans = dict(nonlinearity=self.nonlinearity, @@ -1372,8 +1255,6 @@ class IFSsimulator(): intscale=self.intscale, save_cosmicrays=self.save_cosmicrays) - ##################################################################### - ############################################################################ def _createEmpty(self): @@ -1393,6 +1274,7 @@ class IFSsimulator(): ############################################################################## ############################################################################## + def zodiacal(self, ra, dec, time): """ @@ -1489,7 +1371,6 @@ class IFSsimulator(): # def smoothingWithChargeDiffusion(self, image, sigma=(0.32, 0.32)): # """ - # Parameters # ---------- # image : TYPE @@ -1540,7 +1421,7 @@ class IFSsimulator(): """ self.log.info('Reading in cosmic ray information from %s and %s' % (self.information['cosmicraylengths'], self.information['cosmicraydistance'])) - #### print(self.information) + # print(self.information) crLengths = np.loadtxt( self.information['dir_path']+self.information['cosmicraylengths']) @@ -1551,88 +1432,6 @@ 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 - - # :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=image) - - # # update and verify the header - # hdu.header.add_history('Created by IFSsim at %s' % - # datetime.isoformat(datetime.utcnow())) - # hdu.verify('fix') - - # ofd.append(hdu) - - # # 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 - - # :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) - ############################################################################### def configure(self, simnumber): """ @@ -1826,94 +1625,11 @@ class IFSsimulator(): 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) - - # np.random.seed(5*self.simnumber) - # self.flat_b = np.random.normal(loc=ave, scale=sigma, size=(2048, 4096)) - - # np.random.seed(55*self.simnumber) - # self.flat_r = np.random.normal(loc=ave, scale=sigma, size=(3072, 6144)) - - # s1 = self.flat_b - - # hdu1 = fits.PrimaryHDU(s1) - - # hdu1.header.set('sigma', sigma) - - # dtime = datetime.utcnow().strftime('%Y -%m -%d %H: %M: %S') - # hdu1.header.add_history( - # 'flat image of blue channel is generated on :'+dtime) - - # f1 = '../flat_Blue_'+str(sigma)+'.fits' - # fits.writeto(f1, s1, header=hdu1.header, overwrite=True) - - # s2 = self.flat_r - - # hdu1 = fits.PrimaryHDU(s2) - - # hdu1.header.set('sigma', sigma) - - # dtime = datetime.utcnow().strftime('%Y -%m -%d %H: %M: %S') - # hdu1.header.add_history( - # 'flat image of red channel is generated on :'+dtime) - - # f2 = '../flat_Red_'+str(sigma)+'.fits' - # fits.writeto(f2, s2, header=hdu1.header, overwrite=True) - - # return self.flat_b, self.flat_r - - ########################################################################## - - # def addLampFlux(self): - # """ - - - # Returns - # ------- - # None. - - # """ - # """ - # Include flux from the calibration source. - # """ - - # self.image_b += fits.getdata(self.information['flatflux']) - # self.image_r += fits.getdata(self.information['flatflux']) - - # self.log.info('Flux from the calibration unit included (%s)' % - # self.information['flatflux']) +############################################################################# - ############################################################################# def MakeFlatMatrix(self, img, seed): """ - + Parameters ---------- @@ -1949,7 +1665,7 @@ class IFSsimulator(): def getFrameTransferImg(self): """ - + Returns ------- @@ -1971,17 +1687,18 @@ class IFSsimulator(): # blue channle b_4 = np.sum(self.image_b[0:1024, 0:2048], axis=0) - - ###### flip right to left of part3 and part2 - temp_b_3=np.fliplr(self.image_b[0:1024, 2048:4096]) ## left to right + + # flip right to left of part3 and part2 + temp_b_3 = np.fliplr(self.image_b[0:1024, 2048:4096]) # left to right b_3 = np.sum(temp_b_3, axis=0) - temp_b_2=np.fliplr(self.image_b[1024:1024+1024, 2048:2048+2048]) ## left to right + temp_b_2 = np.fliplr( + self.image_b[1024:1024+1024, 2048:2048+2048]) # left to right b_2 = np.sum(temp_b_2, axis=0) - + b_1 = np.sum(self.image_b[1024:1024+1024, 0:2048], axis=0) #### - + for k in range(1024+320): # part self.frame_b_4[k, :] += b_4*0.09216/self.information['exptime'] @@ -1992,21 +1709,19 @@ class IFSsimulator(): ################### # red channle r_4 = np.sum(self.image_r[0:1536, 0:3072], axis=0) - - - ###### flip right to left of part3 and part2 - temp_r_3=np.fliplr(self.image_r[0:1536, 3072:6144]) ## left to right + + # flip right to left of part3 and part2 + temp_r_3 = np.fliplr(self.image_r[0:1536, 3072:6144]) # left to right r_3 = np.sum(temp_r_3, axis=0) - temp_r_2=np.fliplr(self.image_r[1536:1536+1536, 3072:6144]) ## left to right + temp_r_2 = np.fliplr( + self.image_r[1536:1536+1536, 3072:6144]) # left to right r_2 = np.sum(temp_r_2, axis=0) - - + #r_3 = np.sum(self.image_r[0:1536, 3072:6144], axis=0) #r_2 = np.sum(self.image_r[1536:1536+1536, 3072:6144], axis=0) - - + r_1 = np.sum(self.image_r[1536:1536+1536, 0:3072], axis=0) for k in range(1536+320): @@ -2020,7 +1735,7 @@ class IFSsimulator(): ############################################################################### def applyflatfield(self, simnumber): """ - + Parameters ---------- @@ -2058,7 +1773,7 @@ class IFSsimulator(): ############################################################################### def addCosmicRays(self, idk): """ - + Parameters ---------- @@ -2121,75 +1836,12 @@ class IFSsimulator(): ########################################################## ######################################################################### - # def addReadoutTrails(self): - # """ - - - # Returns - # ------- - # None. - - # """ - # """ - # Add readout trails resulting from reading out the shutter open. - - # Quadrants assumed to be numbered: - # 2 3 - # 0 1 - # """ - # flux_ratio = self.information['readouttime'] / float( - # self.information['bluesize']) / self.information['exptime'] - - # # make a copy, this will be updated - # data = self.image_b.copy() - # # Amplifier at different positions depending on the quadrant number ! - # # left side is 0, 2 and right side is 1, 3 starting from bottom i.e. - # # going clock wise from lower left we have 0, 2, 3, and 1 quadrants. - # if self.information['quadrant'] in (2, 3): - # data = data[::-1, :] - - # data_shift = data.copy() * flux_ratio - # size1, size2 = data.shape - - # for i in range(1, size2, 1): - # data_shift2 = np.roll(data_shift, i, axis=0) - # data_shift2[:i, :] = 0.0 - # data += data_shift2 - - # if self.information['quadrant'] in (2, 3): - # self.image_b = data[::-1, :] - # else: - # self.image_b = data - - - # flux_ratio = self.information['readouttime'] / float( - # self.information['redsize']) / self.information['exptime'] - # # make a copy, this will be updated - # data = self.image_r.copy() - # # Amplifier at different positions depending on the quadrant number ! - # # left side is 0, 2 and right side is 1, 3 starting from bottom i.e. - # # going clock wise from lower left we have 0, 2, 3, and 1 quadrants. - # if self.information['quadrant'] in (2, 3): - # data = data[::-1, :] - - # data_shift = data.copy() * flux_ratio - # size1, size2 = data.shape - - # for i in range(1, size2, 1): - # data_shift2 = np.roll(data_shift, i, axis=0) - # data_shift2[:i, :] = 0.0 - # data += data_shift2 - - # if self.information['quadrant'] in (2, 3): - # self.image_r = data[::-1, :] - # else: - # self.image_r = data ############################################################################## def applyDarkCurrent(self): """ - + Returns ------- @@ -2238,60 +1890,11 @@ class IFSsimulator(): ############################################################################## - - # 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. - # """ - - # # add background - # bcgr = self.information['exptime'] * self.information['cosmic_bkgd'] - - # self.image_b += bcgr - # self.image_r += bcgr - - # self.log.info('Added cosmic background = %f' % bcgr) - - # if self.cosmicRays: - - # 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'] - - # self.image_b += sl - # self.image_r += sl - - # self.log.info('Added scattered light = %f' % sl) - ############################################################################## + def applyPoissonNoise(self): """ - + Returns ------- @@ -2322,11 +1925,11 @@ class IFSsimulator(): self.image_r += residual -################################################################################################################### +############################################################################## def applyCosmetics(self): """ - + Returns ------- @@ -2399,7 +2002,7 @@ class IFSsimulator(): def applyRadiationDamage(self): """ - + Returns ------- @@ -2415,13 +2018,11 @@ class IFSsimulator(): self.log.debug('Starting to apply radiation damage model...') # at this point we can give fake data... cti = CTI.CDM03bidir(self.information, [], log=self.log) - - ### here we need the right input data + + # 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... @@ -2434,7 +2035,7 @@ class IFSsimulator(): def applyNonlinearity(self): """ - + Returns ------- @@ -2460,7 +2061,7 @@ class IFSsimulator(): def applyReadoutNoise(self): """ - + Returns ------- @@ -2497,7 +2098,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)) - self.log.info('readnoise added in blue channel') # red zone 4 , OSH### 1536* 3072 @@ -2524,7 +2124,7 @@ class IFSsimulator(): def appFrameTransferEffect(self): """ - + Returns ------- @@ -2583,7 +2183,7 @@ class IFSsimulator(): def electrons2ADU(self): """ - + Returns ------- @@ -2628,7 +2228,7 @@ class IFSsimulator(): def applyBias(self): """ - + Returns ------- @@ -2679,7 +2279,7 @@ class IFSsimulator(): ############################################################################### def addPreOverScans(self): """ - + Returns ------- @@ -2735,7 +2335,7 @@ class IFSsimulator(): def applyBleeding_yan(self): """ - + Returns ------- @@ -2837,7 +2437,7 @@ class IFSsimulator(): def discretise(self, max=2**16-1): """ - + Parameters ---------- @@ -2878,39 +2478,38 @@ class IFSsimulator(): ################################################################################################## - def applyImageShift(self): - """ - + # def applyImageShift(self): + # """ - Returns - ------- - None. - """ + # 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'] + # np.random.seed(9*self.simnumber) + # ud = np.random.random() # Choose a random rotation + # dx = 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') + # 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'] + # 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 # ------- @@ -2939,9 +2538,10 @@ class IFSsimulator(): # 'Applied telescope rotation with angle (in degree)= %f.', angle) ############################################################################### + def CCDreadout(self): """ - + Returns ------- @@ -3029,7 +2629,7 @@ class IFSsimulator(): def writeOutputs(self, obnum): """ - + Parameters ---------- @@ -3301,8 +2901,8 @@ class IFSsimulator(): 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: + + if self.information['exposuretimes'] > 1: hdu_b.header['NAXIS3'] = ( np.int32(self.information['exposuretimes']), 'length of third data axis') @@ -3681,8 +3281,8 @@ class IFSsimulator(): hdu_r.header['NAXIS2'] = ( np.int32(13768), 'length of second data axis') - - if self.information['exposuretimes']>1: + + if self.information['exposuretimes'] > 1: hdu_r.header['NAXIS3'] = ( np.int32(self.information['exposuretimes']), 'length of third data axis') @@ -3993,10 +3593,9 @@ class IFSsimulator(): ################################################################################### ################################################################################## - def CalskyNoise(self, lam): """ - + Parameters ---------- @@ -4042,7 +3641,7 @@ class IFSsimulator(): def sim_lamp_hole_img(self, exposuretime): """ - + Parameters ---------- @@ -4217,8 +3816,7 @@ class IFSsimulator(): ####################################################### - # energy=energy+img0[500-302:500+320,500-320:500+320].sum() ## calculate the slice image energy; - + # CCD quantum efficiency CCD_Qe_lam = np.interp(lam, self.CCD_Qe[:, 0], self.CCD_Qe[:, 1]) @@ -4310,11 +3908,6 @@ class IFSsimulator(): photons_blue.addTo(blue_img) - # zblue=blue_img.array - - # 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) @@ -4380,7 +3973,7 @@ class IFSsimulator(): def sim_sky_img(self, skyfitsfilename, skyRa, skyDec, sky_rot, telRa, telDec, exposuretime, simnumber): """ - + Parameters ---------- @@ -4901,8 +4494,7 @@ class IFSsimulator(): energy_blue = energy_blue+sum(photons_blue.flux) - width_blue = width_blue+deltalam/32.0 - + ################# # fits.writeto('blueImg.fits',blue_img.array,overwrite=True) @@ -4965,7 +4557,7 @@ class IFSsimulator(): ################################################################################################# def sim_calibration_img(self, exposuretime, source): """ - + Parameters ---------- @@ -5144,18 +4736,6 @@ class IFSsimulator(): conv[conv < 0.0] = 0.0 img0 = conv - # temp=img0*maskSlice['31']*maskSlit['31'] - # print(image0.array.sum()) - # print(img1.sum()) - # fits.writeto('img.fits',image0.array,overwrite=True) - # fits.writeto('convimg.fits',img0,overwrite=True) - - ####### - - ####################################################### - - # 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]) @@ -5263,9 +4843,8 @@ class IFSsimulator(): photons_blue.y = photons_blue.y/self.pixelscale + \ dy_blue+self.slice_blue['py'][k] - - - #photons_blue.y =2048-photons_blue.y + + #photons_blue.y =2048-photons_blue.y blue_sensor.accumulate(photons_blue, blue_img) @@ -5287,8 +4866,7 @@ class IFSsimulator(): 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) @@ -5348,7 +4926,7 @@ class IFSsimulator(): def simulate(self, sourcein, simnumber): """ - + Parameters ---------- @@ -5370,19 +4948,18 @@ class IFSsimulator(): :return: None """ - self.source = sourcein self.simnumber = simnumber self.configure(simnumber) # print the configfile name and path; - - self.debug=self.information['debug'] + + self.debug = self.information['debug'] self.dt = datetime.utcnow() if self.source == 'SCI' or self.source == 'COMP': - if simnumber <= 50 and simnumber>20: + if simnumber <= 50 and simnumber > 20: self.information['exptime'] = 300 if simnumber <= 100 and simnumber > 50: @@ -5536,20 +5113,18 @@ class IFSsimulator(): exptimes = self.information['exposuretimes'] else: exptimes = 1 - self.information['exposuretimes'] =1 - - + self.information['exposuretimes'] = 1 + # if self.debug: # exptimes=1 - if 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): @@ -5617,16 +5192,14 @@ class IFSsimulator(): print('Applying cosmetics finished.......') self.discretise() - - - if exptimes>1: + + 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) @@ -5650,7 +5223,7 @@ class IFSsimulator(): def runIFSsim(sourcein, configfile, dir_path, iLoop, debug, applyhole='no'): """ - + Parameters ---------- @@ -5680,11 +5253,11 @@ def runIFSsim(sourcein, configfile, dir_path, iLoop, debug, applyhole='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['debug'] = debug + simulate[iLoop].simulate(sourcein, iLoop) return 1 -############################## end ########################################## \ No newline at end of file +############################## end ##########################################