Skip to content
GitLab
Commit 5070c190 authored by Yan Zhaojun's avatar Yan Zhaojun
Browse files

update

parent 8165a134
Pipeline #7134 failed with stage
in 0 seconds
Hide whitespace changes
Inline Side-by-side
csst_ifs_sim/CTI/CTI.py
View file @ 5070c190
...@@ -15,7 +15,7 @@ parameters in parallel and serial direction. ...@@ -15,7 +15,7 @@ parameters in parallel and serial direction.
import numpy as np import numpy as np
#CDM03bidir # CDM03bidir
class CDM03bidir(): class CDM03bidir():
""" """
Class to run CDM03 CTI model, class Fortran routine to perform the actual CDM03 calculations. Class to run CDM03 CTI model, class Fortran routine to perform the actual CDM03 calculations.
...@@ -27,6 +27,7 @@ class CDM03bidir(): ...@@ -27,6 +27,7 @@ class CDM03bidir():
:param log: instance to Python logging :param log: instance to Python logging
:type log: logging instance :type log: logging instance
""" """
def __init__(self, settings, data, log=None): def __init__(self, settings, data, log=None):
""" """
Class constructor. Class constructor.
...@@ -39,7 +40,8 @@ class CDM03bidir(): ...@@ -39,7 +40,8 @@ class CDM03bidir():
:type log: logging instance :type log: logging instance
""" """
self.data = data self.data = data
self.values = dict(quads=(0,1,2,3), xsize=2048, ysize=2066, dob=0.0, rdose=8.0e9) self.values = dict(quads=(0, 1, 2, 3), xsize=2048,
ysize=2066, dob=0.0, rdose=8.0e9)
self.values.update(settings) self.values.update(settings)
self.log = log self.log = log
self._setupLogger() self._setupLogger()
...@@ -47,44 +49,45 @@ class CDM03bidir(): ...@@ -47,44 +49,45 @@ class CDM03bidir():
# #default CDM03 settings # #default CDM03 settings
# self.params = dict(beta_p=0.6, beta_s=0.6, fwc=200000., vth=1.168e7, vg=6.e-11, t=20.48e-3, # self.params = dict(beta_p=0.6, beta_s=0.6, fwc=200000., vth=1.168e7, vg=6.e-11, t=20.48e-3,
# sfwc=730000., svg=1.0e-10, st=5.0e-6, parallel=1., serial=1.) # sfwc=730000., svg=1.0e-10, st=5.0e-6, parallel=1., serial=1.)
self.params = dict(beta_p=0.6, beta_s=0.6, fwc=100000., vth=1.168e7, vg=6.e-11, t=1.0e-3, self.params = dict(beta_p=0.6, beta_s=0.6, fwc=100000., vth=1.168e7, vg=6.e-11, t=1.0e-3,
sfwc=700000., svg=1.0e-10, st=1.0e-6, parallel=1., serial=0.) sfwc=700000., svg=1.0e-10, st=1.0e-6, parallel=1., serial=0.)
#update with inputs # update with inputs
self.params.update(self.values) self.params.update(self.values)
#read in trap information # read in trap information
trapdata = np.loadtxt(self.values['dir_path']+self.values['paralleltrapfile']) trapdata = np.loadtxt(
self.values['dir_path']+self.values['paralleltrapfile'])
if trapdata.ndim > 1: if trapdata.ndim > 1:
self.nt_p = trapdata[:, 0] self.nt_p = trapdata[:, 0]
self.sigma_p = trapdata[:, 1] self.sigma_p = trapdata[:, 1]
self.taur_p = trapdata[:, 2] self.taur_p = trapdata[:, 2]
else: else:
#only one trap species # only one trap species
self.nt_p = [trapdata[0],] self.nt_p = [trapdata[0],]
self.sigma_p = [trapdata[1],] self.sigma_p = [trapdata[1],]
self.taur_p = [trapdata[2],] self.taur_p = [trapdata[2],]
trapdata = np.loadtxt(self.values['dir_path']+self.values['serialtrapfile']) trapdata = np.loadtxt(
self.values['dir_path']+self.values['serialtrapfile'])
if trapdata.ndim > 1: if trapdata.ndim > 1:
self.nt_s = trapdata[:, 0] self.nt_s = trapdata[:, 0]
self.sigma_s = trapdata[:, 1] self.sigma_s = trapdata[:, 1]
self.taur_s = trapdata[:, 2] self.taur_s = trapdata[:, 2]
else: else:
#only one trap species # only one trap species
self.nt_s = [trapdata[0],] self.nt_s = [trapdata[0],]
self.sigma_s = [trapdata[1],] self.sigma_s = [trapdata[1],]
self.taur_s = [trapdata[2],] self.taur_s = [trapdata[2],]
#scale thibaut's values # scale thibaut's values
if 'thibaut' in self.values['parallelTrapfile']: if 'thibaut' in self.values['parallelTrapfile']:
self.nt_p /= 0.576 #thibaut's values traps / pixel self.nt_p /= 0.576 # thibaut's values traps / pixel
self.sigma_p *= 1.e4 #thibaut's values in m**2 self.sigma_p *= 1.e4 # thibaut's values in m**2
if 'thibaut' in self.values['serialTrapfile']: if 'thibaut' in self.values['serialTrapfile']:
self.nt_s *= 0.576 #thibaut's values traps / pixel #should be division? self.nt_s *= 0.576 # thibaut's values traps / pixel #should be division?
self.sigma_s *= 1.e4 #thibaut's values in m**2 self.sigma_s *= 1.e4 # thibaut's values in m**2
def _setupLogger(self): def _setupLogger(self):
""" """
...@@ -94,7 +97,6 @@ class CDM03bidir(): ...@@ -94,7 +97,6 @@ class CDM03bidir():
# if self.log is None: # if self.log is None:
# self.logger = False # self.logger = False
def applyRadiationDamage(self, data, iquadrant=0): def applyRadiationDamage(self, data, iquadrant=0):
""" """
Apply radian damage based on FORTRAN CDM03 model. The method assumes that Apply radian damage based on FORTRAN CDM03 model. The method assumes that
...@@ -133,8 +135,8 @@ class CDM03bidir(): ...@@ -133,8 +135,8 @@ class CDM03bidir():
:return: image that has been run through the CDM03 model :return: image that has been run through the CDM03 model
:rtype: ndarray :rtype: ndarray
""""" """""
#return data # return data
iflip = iquadrant / 2 iflip = iquadrant / 2
jflip = iquadrant % 2 jflip = iquadrant % 2
...@@ -158,27 +160,22 @@ class CDM03bidir(): ...@@ -158,27 +160,22 @@ class CDM03bidir():
self.log.info('jflip=%i' % jflip) self.log.info('jflip=%i' % jflip)
################################################################################# #################################################################################
###modify # modify
#sys.path.append('../so') # sys.path.append('../so')
from ifs_so import cdm03bidir from ifs_so import cdm03bidir
# from ifs_so.cdm03.cpython-38-x86_64-linux-gnu import cdm03bidir # from ifs_so.cdm03.cpython-38-x86_64-linux-gnu import cdm03bidir
# import cdm03bidir # import cdm03bidir
CTIed = cdm03bidir.cdm03(np.asfortranarray(data), CTIed = cdm03bidir.cdm03(np.asfortranarray(data),
jflip, iflip, jflip, iflip,
self.values['dob'], self.values['rdose'], self.values['dob'], self.values['rdose'],
self.nt_p, self.sigma_p, self.taur_p, self.nt_p, self.sigma_p, self.taur_p,
self.nt_s, self.sigma_s, self.taur_s, self.nt_s, self.sigma_s, self.taur_s,
params, params,
[data.shape[0], data.shape[1], len(self.nt_p), len(self.nt_s), len(self.params)]) [data.shape[0], data.shape[1], len(self.nt_p), len(self.nt_s), len(self.params)])
return np.asanyarray(CTIed) return np.asanyarray(CTIed)
#################################################################################################################
#################################################################################################################
################################################################################################################# #################################################################################################################
csst_ifs_sim/csst_ifs_sim.py
View file @ 5070c190
...@@ -5,6 +5,9 @@ Created on Thu Apr 11 15:18:57 2024 ...@@ -5,6 +5,9 @@ Created on Thu Apr 11 15:18:57 2024
@author: yan @author: yan
""" """
from scipy.interpolate import InterpolatedUnivariateSpline
import logging.handlers
import logging
from astropy.utils.iers import conf from astropy.utils.iers import conf
import galsim import galsim
from scipy.integrate import simps from scipy.integrate import simps
...@@ -28,7 +31,7 @@ import astropy.coordinates as coord ...@@ -28,7 +31,7 @@ import astropy.coordinates as coord
import ctypes import ctypes
import sys import sys
##sys.path.append('./csst_ifs_sim') # sys.path.append('./csst_ifs_sim')
conf.auto_max_age = None conf.auto_max_age = None
""" """
...@@ -104,6 +107,7 @@ class CDM03bidir(): ...@@ -104,6 +107,7 @@ class CDM03bidir():
:param log: instance to Python logging :param log: instance to Python logging
:type log: logging instance :type log: logging instance
""" """
def __init__(self, settings, data, log=None): def __init__(self, settings, data, log=None):
""" """
Class constructor. Class constructor.
...@@ -116,7 +120,8 @@ class CDM03bidir(): ...@@ -116,7 +120,8 @@ class CDM03bidir():
:type log: logging instance :type log: logging instance
""" """
self.data = data self.data = data
self.values = dict(quads=(0,1,2,3), xsize=2048, ysize=2066, dob=0.0, rdose=8.0e9) self.values = dict(quads=(0, 1, 2, 3), xsize=2048,
ysize=2066, dob=0.0, rdose=8.0e9)
self.values.update(settings) self.values.update(settings)
self.log = log self.log = log
self._setupLogger() self._setupLogger()
...@@ -124,44 +129,45 @@ class CDM03bidir(): ...@@ -124,44 +129,45 @@ class CDM03bidir():
# #default CDM03 settings # #default CDM03 settings
# self.params = dict(beta_p=0.6, beta_s=0.6, fwc=200000., vth=1.168e7, vg=6.e-11, t=20.48e-3, # self.params = dict(beta_p=0.6, beta_s=0.6, fwc=200000., vth=1.168e7, vg=6.e-11, t=20.48e-3,
# sfwc=730000., svg=1.0e-10, st=5.0e-6, parallel=1., serial=1.) # sfwc=730000., svg=1.0e-10, st=5.0e-6, parallel=1., serial=1.)
self.params = dict(beta_p=0.6, beta_s=0.6, fwc=100000., vth=1.168e7, vg=6.e-11, t=1.0e-3, self.params = dict(beta_p=0.6, beta_s=0.6, fwc=100000., vth=1.168e7, vg=6.e-11, t=1.0e-3,
sfwc=700000., svg=1.0e-10, st=1.0e-6, parallel=1., serial=0.) sfwc=700000., svg=1.0e-10, st=1.0e-6, parallel=1., serial=0.)
#update with inputs # update with inputs
self.params.update(self.values) self.params.update(self.values)
#read in trap information # read in trap information
trapdata = np.loadtxt(self.values['dir_path']+self.values['paralleltrapfile']) trapdata = np.loadtxt(
self.values['dir_path']+self.values['paralleltrapfile'])
if trapdata.ndim > 1: if trapdata.ndim > 1:
self.nt_p = trapdata[:, 0] self.nt_p = trapdata[:, 0]
self.sigma_p = trapdata[:, 1] self.sigma_p = trapdata[:, 1]
self.taur_p = trapdata[:, 2] self.taur_p = trapdata[:, 2]
else: else:
#only one trap species # only one trap species
self.nt_p = [trapdata[0],] self.nt_p = [trapdata[0],]
self.sigma_p = [trapdata[1],] self.sigma_p = [trapdata[1],]
self.taur_p = [trapdata[2],] self.taur_p = [trapdata[2],]
trapdata = np.loadtxt(self.values['dir_path']+self.values['serialtrapfile']) trapdata = np.loadtxt(
self.values['dir_path']+self.values['serialtrapfile'])
if trapdata.ndim > 1: if trapdata.ndim > 1:
self.nt_s = trapdata[:, 0] self.nt_s = trapdata[:, 0]
self.sigma_s = trapdata[:, 1] self.sigma_s = trapdata[:, 1]
self.taur_s = trapdata[:, 2] self.taur_s = trapdata[:, 2]
else: else:
#only one trap species # only one trap species
self.nt_s = [trapdata[0],] self.nt_s = [trapdata[0],]
self.sigma_s = [trapdata[1],] self.sigma_s = [trapdata[1],]
self.taur_s = [trapdata[2],] self.taur_s = [trapdata[2],]
#scale thibaut's values # scale thibaut's values
if 'thibaut' in self.values['parallelTrapfile']: if 'thibaut' in self.values['parallelTrapfile']:
self.nt_p /= 0.576 #thibaut's values traps / pixel self.nt_p /= 0.576 # thibaut's values traps / pixel
self.sigma_p *= 1.e4 #thibaut's values in m**2 self.sigma_p *= 1.e4 # thibaut's values in m**2
if 'thibaut' in self.values['serialTrapfile']: if 'thibaut' in self.values['serialTrapfile']:
self.nt_s *= 0.576 #thibaut's values traps / pixel #should be division? self.nt_s *= 0.576 # thibaut's values traps / pixel #should be division?
self.sigma_s *= 1.e4 #thibaut's values in m**2 self.sigma_s *= 1.e4 # thibaut's values in m**2
def _setupLogger(self): def _setupLogger(self):
""" """
...@@ -169,7 +175,6 @@ class CDM03bidir(): ...@@ -169,7 +175,6 @@ class CDM03bidir():
""" """
self.logger = True self.logger = True
def applyRadiationDamage(self, data, iquadrant=0): def applyRadiationDamage(self, data, iquadrant=0):
""" """
Apply radian damage based on FORTRAN CDM03 model. The method assumes that Apply radian damage based on FORTRAN CDM03 model. The method assumes that
...@@ -208,8 +213,8 @@ class CDM03bidir(): ...@@ -208,8 +213,8 @@ class CDM03bidir():
:return: image that has been run through the CDM03 model :return: image that has been run through the CDM03 model
:rtype: ndarray :rtype: ndarray
""""" """""
#return data # return data
iflip = iquadrant / 2 iflip = iquadrant / 2
jflip = iquadrant % 2 jflip = iquadrant % 2
...@@ -235,10 +240,10 @@ class CDM03bidir(): ...@@ -235,10 +240,10 @@ class CDM03bidir():
############################################ ############################################
################################################################################# #################################################################################
###modify # modify
#sys.path.append('../so') # sys.path.append('../so')
#from .ifs_so import cdm03bidir #from .ifs_so import cdm03bidir
try: try:
from ifs_so import cdm03bidir from ifs_so import cdm03bidir
except: except:
from .ifs_so import cdm03bidir from .ifs_so import cdm03bidir
...@@ -246,20 +251,20 @@ class CDM03bidir(): ...@@ -246,20 +251,20 @@ class CDM03bidir():
# from ifs_so.cdm03.cpython-38-x86_64-linux-gnu import cdm03bidir # from ifs_so.cdm03.cpython-38-x86_64-linux-gnu import cdm03bidir
# import cdm03bidir # import cdm03bidir
CTIed = cdm03bidir.cdm03(np.asfortranarray(data), CTIed = cdm03bidir.cdm03(np.asfortranarray(data),
jflip, iflip, jflip, iflip,
self.values['dob'], self.values['rdose'], self.values['dob'], self.values['rdose'],
self.nt_p, self.sigma_p, self.taur_p, self.nt_p, self.sigma_p, self.taur_p,
self.nt_s, self.sigma_s, self.taur_s, self.nt_s, self.sigma_s, self.taur_s,
params, params,
[data.shape[0], data.shape[1], len(self.nt_p), len(self.nt_s), len(self.params)]) [data.shape[0], data.shape[1], len(self.nt_p), len(self.nt_s), len(self.params)])
return np.asanyarray(CTIed) return np.asanyarray(CTIed)
################################################################################################################# #################################################################################################################
################################################################################################################# #################################################################################################################
""" """
These functions can be used for logging information. These functions can be used for logging information.
...@@ -267,8 +272,6 @@ These functions can be used for logging information. ...@@ -267,8 +272,6 @@ These functions can be used for logging information.
:version: 0.3 :version: 0.3
""" """
import logging
import logging.handlers
def lg(log_filename, loggername='logger'): def lg(log_filename, loggername='logger'):
...@@ -285,16 +288,17 @@ def lg(log_filename, loggername='logger'): ...@@ -285,16 +288,17 @@ def lg(log_filename, loggername='logger'):
logger.setLevel(logging.DEBUG) logger.setLevel(logging.DEBUG)
# Add the log message handler to the logger # Add the log message handler to the logger
handler = logging.handlers.RotatingFileHandler(log_filename) handler = logging.handlers.RotatingFileHandler(log_filename)
#maxBytes=20, backupCount=5) # maxBytes=20, backupCount=5)
# create formatter # create formatter
formatter = logging.Formatter('%(asctime)s - %(module)s - %(funcName)s - %(levelname)s - %(message)s') formatter = logging.Formatter(
'%(asctime)s - %(module)s - %(funcName)s - %(levelname)s - %(message)s')
# add formatter to ch # add formatter to ch
handler.setFormatter(formatter) handler.setFormatter(formatter)
# add handler to logger # add handler to logger
if (logger.hasHandlers()): if (logger.hasHandlers()):
logger.handlers.clear() logger.handlers.clear()
logger.addHandler(handler) logger.addHandler(handler)
return logger return logger
...@@ -309,7 +313,7 @@ def IFSinformation(): ...@@ -309,7 +313,7 @@ def IFSinformation():
The file provides a function that returns IFS related information such as pixel The file provides a function that returns IFS related information such as pixel
size, dark current, gain... size, dark current, gain...
Returns a dictionary describing VIS. The following information is provided Returns a dictionary describing VIS. The following information is provided
(id: value - reference):: (id: value - reference)::
...@@ -319,11 +323,11 @@ def IFSinformation(): ...@@ -319,11 +323,11 @@ def IFSinformation():
######################################################################################################### #########################################################################################################
#out = dict(readnoise=4, pixel_size=0.1, dark=0.0008333, fullwellcapacity=90000, bluesize=4000, redsize=6000, readtime=300.) #out = dict(readnoise=4, pixel_size=0.1, dark=0.0008333, fullwellcapacity=90000, bluesize=4000, redsize=6000, readtime=300.)
out=dict() out = dict()
out.update({'dob' : 0, 'rdose' : 8.0e9, out.update({'dob': 0, 'rdose': 8.0e9,
'parallelTrapfile' : 'cdm_euclid_parallel.dat', 'serialTrapfile' : 'cdm_euclid_serial.dat', 'parallelTrapfile': 'cdm_euclid_parallel.dat', 'serialTrapfile': 'cdm_euclid_serial.dat',
'beta_s' : 0.6, 'beta_p': 0.6, 'fwc' : 90000, 'vth' : 1.168e7, 't' : 20.48e-3, 'vg' : 6.e-11, 'beta_s': 0.6, 'beta_p': 0.6, 'fwc': 90000, 'vth': 1.168e7, 't': 20.48e-3, 'vg': 6.e-11,
'st' : 5.0e-6, 'sfwc' : 730000., 'svg' : 1.0e-10}) 'st': 5.0e-6, 'sfwc': 730000., 'svg': 1.0e-10})
return out return out
...@@ -338,16 +342,15 @@ def CCDnonLinearityModel(data, beta=6e-7): ...@@ -338,16 +342,15 @@ def CCDnonLinearityModel(data, beta=6e-7):
:rtype: float or ndarray :rtype: float or ndarray
""" """
out = data-beta*data**2 out = data-beta*data**2
return out return out
############################################################################# #############################################################################
from scipy.interpolate import InterpolatedUnivariateSpline
class cosmicrays(): class cosmicrays():
""" """
Cosmic Rays Cosmic Rays
=========== ===========
...@@ -369,6 +372,7 @@ class cosmicrays(): ...@@ -369,6 +372,7 @@ class cosmicrays():
:param information: cosmic ray track information (file containing track length and energy information) and :param information: cosmic ray track information (file containing track length and energy information) and
exposure time. exposure time.
""" """
def __init__(self, log, image, crInfo=None, information=None): def __init__(self, log, image, crInfo=None, information=None):
""" """
Cosmic ray generation class. Can either draw events from distributions or Cosmic ray generation class. Can either draw events from distributions or
...@@ -380,14 +384,14 @@ class cosmicrays(): ...@@ -380,14 +384,14 @@ class cosmicrays():
:param information: cosmic ray track information (file containing track length and energy information) and :param information: cosmic ray track information (file containing track length and energy information) and
exposure time. exposure time.
""" """
#setup logger # setup logger
self.log = log self.log = log
#image and size #image and size
self.image = image.copy() self.image = image.copy()
self.ysize, self.xsize = self.image.shape self.ysize, self.xsize = self.image.shape
#set up the information dictionary, first with defaults and then overwrite with inputs if given # set up the information dictionary, first with defaults and then overwrite with inputs if given
self.information = (dict(cosmicraylengths='/home/yan/csst-master/data/cdf_cr_length.dat', self.information = (dict(cosmicraylengths='/home/yan/csst-master/data/cdf_cr_length.dat',
cosmicraydistance='/home/yan/csst-master/data/cdf_cr_total.dat', cosmicraydistance='/home/yan/csst-master/data/cdf_cr_total.dat',
exptime=565)) exptime=565))
...@@ -413,10 +417,10 @@ class cosmicrays(): ...@@ -413,10 +417,10 @@ class cosmicrays():
:return: cosmic ray map (image) :return: cosmic ray map (image)
:rtype: nd-array :rtype: nd-array
""" """
#create empty array # create empty array
crImage = np.zeros((self.ysize, self.xsize), dtype=np.float64) crImage = np.zeros((self.ysize, self.xsize), dtype=np.float64)
#x and y shifts # x and y shifts
dx = l * np.cos(phi) / 2. dx = l * np.cos(phi) / 2.
dy = l * np.sin(phi) / 2. dy = l * np.sin(phi) / 2.
mskdx = np.abs(dx) < 1e-8 mskdx = np.abs(dx) < 1e-8
...@@ -446,7 +450,7 @@ class cosmicrays(): ...@@ -446,7 +450,7 @@ class cosmicrays():
jhi[msk] = self.ysize jhi[msk] = self.ysize
jhi = jhi.astype(int) jhi = jhi.astype(int)
#loop over the individual events # loop over the individual events
for i, luminosity in enumerate(lum): for i, luminosity in enumerate(lum):
n = 0 # count the intercepts n = 0 # count the intercepts
...@@ -454,7 +458,7 @@ class cosmicrays(): ...@@ -454,7 +458,7 @@ class cosmicrays():
x = [] x = []
y = [] y = []
#Compute X intercepts on the pixel grid # Compute X intercepts on the pixel grid
if ilo[i] < ihi[i]: if ilo[i] < ihi[i]:
for xcoord in range(ilo[i], ihi[i]): for xcoord in range(ilo[i], ihi[i]):
ok = (xcoord - x0[i]) / dx[i] ok = (xcoord - x0[i]) / dx[i]
...@@ -472,7 +476,7 @@ class cosmicrays(): ...@@ -472,7 +476,7 @@ class cosmicrays():
x.append(xcoord) x.append(xcoord)
y.append(y0[i] + ok * dy[i]) y.append(y0[i] + ok * dy[i])
#Compute Y intercepts on the pixel grid # Compute Y intercepts on the pixel grid
if jlo[i] < jhi[i]: if jlo[i] < jhi[i]:
for ycoord in range(jlo[i], jhi[i]): for ycoord in range(jlo[i], jhi[i]):
ok = (ycoord - y0[i]) / dy[i] ok = (ycoord - y0[i]) / dy[i]
...@@ -490,13 +494,13 @@ class cosmicrays(): ...@@ -490,13 +494,13 @@ class cosmicrays():
x.append(x0[i] + ok * dx[i]) x.append(x0[i] + ok * dx[i])
y.append(ycoord) y.append(ycoord)
#check if no intercepts were found # check if no intercepts were found
if n < 1: if n < 1:
xc = int(np.floor(x0[i])) xc = int(np.floor(x0[i]))
yc = int(np.floor(y0[i])) yc = int(np.floor(y0[i]))
crImage[yc, xc] += luminosity crImage[yc, xc] += luminosity
#Find the arguments that sort the intersections along the track # Find the arguments that sort the intersections along the track
u = np.asarray(u) u = np.asarray(u)
x = np.asarray(x) x = np.asarray(x)
y = np.asarray(y) y = np.asarray(y)
...@@ -507,7 +511,7 @@ class cosmicrays(): ...@@ -507,7 +511,7 @@ class cosmicrays():
x = x[args] x = x[args]
y = y[args] y = y[args]
#Decide which cell each interval traverses, and the path length # Decide which cell each interval traverses, and the path length
for i in range(1, n - 1): for i in range(1, n - 1):
w = u[i + 1] - u[i] w = u[i + 1] - u[i]
cx = int(1 + np.floor((x[i + 1] + x[i]) / 2.)) cx = int(1 + np.floor((x[i + 1] + x[i]) / 2.))
...@@ -518,8 +522,6 @@ class cosmicrays(): ...@@ -518,8 +522,6 @@ class cosmicrays():
return crImage return crImage
def _drawEventsToCoveringFactor(self, coveringFraction=3.0, limit=1000, verbose=False): def _drawEventsToCoveringFactor(self, coveringFraction=3.0, limit=1000, verbose=False):
""" """
Generate cosmic ray events up to a covering fraction and include it to a cosmic ray map (self.cosmicrayMap). Generate cosmic ray events up to a covering fraction and include it to a cosmic ray map (self.cosmicrayMap).
...@@ -534,55 +536,57 @@ class cosmicrays(): ...@@ -534,55 +536,57 @@ class cosmicrays():
:return: None :return: None
""" """
if coveringFraction>1 or coveringFraction<0.1: if coveringFraction > 1 or coveringFraction < 0.1:
self.log.error( self.log.error(
'coveringFraction error, it shoub be in [0.1, 1]!') 'coveringFraction error, it shoub be in [0.1, 1]!')
raise ValueError( raise ValueError(
'coveringFraction error, it shoub be in [0.1, 1]!') 'coveringFraction error, it shoub be in [0.1, 1]!')
self.cosmicrayMap = np.zeros((self.ysize, self.xsize)) self.cosmicrayMap = np.zeros((self.ysize, self.xsize))
#how many events to draw at once, too large number leads to exceeding the covering fraction # how many events to draw at once, too large number leads to exceeding the covering fraction
cr_n = int(295 * self.information['exptime'] / 565. * coveringFraction / 1.4) cr_n = int(
295 * self.information['exptime'] / 565. * coveringFraction / 1.4)
covering = 0.0 covering = 0.0
while covering < coveringFraction: while covering < coveringFraction:
#pseudo-random numbers taken from a uniform distribution between 0 and 1 # pseudo-random numbers taken from a uniform distribution between 0 and 1
np.random.seed() np.random.seed()
luck = np.random.rand(cr_n) luck = np.random.rand(cr_n)
#draw the length of the tracks # draw the length of the tracks
ius = InterpolatedUnivariateSpline(self.cr['cr_cdf'], self.cr['cr_u']) ius = InterpolatedUnivariateSpline(
self.cr['cr_cdf'], self.cr['cr_u'])
self.cr['cr_l'] = ius(luck) self.cr['cr_l'] = ius(luck)
if limit is None: if limit is None:
ius = InterpolatedUnivariateSpline(self.cr['cr_cde'], self.cr['cr_v']) ius = InterpolatedUnivariateSpline(
self.cr['cr_cde'], self.cr['cr_v'])
self.cr['cr_e'] = ius(luck) self.cr['cr_e'] = ius(luck)
else: else:
#set the energy directly to the limit # set the energy directly to the limit
self.cr['cr_e'] = np.asarray([limit,]) self.cr['cr_e'] = np.asarray([limit,])
#Choose the properties such as positions and an angle from a random Uniform dist # Choose the properties such as positions and an angle from a random Uniform dist
np.random.seed() np.random.seed()
cr_x = self.xsize * np.random.rand(int(np.floor(cr_n))) cr_x = self.xsize * np.random.rand(int(np.floor(cr_n)))
np.random.seed() np.random.seed()
cr_y = self.ysize * np.random.rand(int(np.floor(cr_n))) cr_y = self.ysize * np.random.rand(int(np.floor(cr_n)))
np.random.seed() np.random.seed()
cr_phi = np.pi * np.random.rand(int(np.floor(cr_n))) cr_phi = np.pi * np.random.rand(int(np.floor(cr_n)))
#find the intercepts # find the intercepts
self.cosmicrayMap += self._cosmicRayIntercepts(self.cr['cr_e'], cr_x, cr_y, self.cr['cr_l'], cr_phi) self.cosmicrayMap += self._cosmicRayIntercepts(
self.cr['cr_e'], cr_x, cr_y, self.cr['cr_l'], cr_phi)
#count the covering factor # count the covering factor
area_cr = np.count_nonzero(self.cosmicrayMap) area_cr = np.count_nonzero(self.cosmicrayMap)
covering = 100.*area_cr / (self.xsize*self.ysize) covering = 100.*area_cr / (self.xsize*self.ysize)
def addUpToFraction(self, coveringFraction, limit=None, verbose=False): def addUpToFraction(self, coveringFraction, limit=None, verbose=False):
""" """
Add cosmic ray events up to the covering Fraction. Add cosmic ray events up to the covering Fraction.
...@@ -597,9 +601,10 @@ class cosmicrays(): ...@@ -597,9 +601,10 @@ class cosmicrays():
:return: image with cosmic rays :return: image with cosmic rays
:rtype: ndarray :rtype: ndarray
""" """
self._drawEventsToCoveringFactor(coveringFraction, limit=limit, verbose=verbose) self._drawEventsToCoveringFactor(
coveringFraction, limit=limit, verbose=verbose)
#paste cosmic rays # paste cosmic rays
self.image += self.cosmicrayMap self.image += self.cosmicrayMap
return self.image return self.image
...@@ -1148,6 +1153,8 @@ def beta_angle(x_sat, y_sat, z_sat, vx_sat, vy_sat, vz_sat, ra_obj, dec_obj): ...@@ -1148,6 +1153,8 @@ def beta_angle(x_sat, y_sat, z_sat, vx_sat, vy_sat, vz_sat, ra_obj, dec_obj):
return angle return angle
############################################################################### ###############################################################################
def find_min(arr): def find_min(arr):
min_val = arr[0] min_val = arr[0]
min_index = 0 min_index = 0
...@@ -1158,6 +1165,8 @@ def find_min(arr): ...@@ -1158,6 +1165,8 @@ def find_min(arr):
return min_val, min_index return min_val, min_index
################################# #################################
def find_max(arr): def find_max(arr):
max_val = arr[0] max_val = arr[0]
max_index = 0 max_index = 0
...@@ -1557,11 +1566,12 @@ def get_dx_dy_blue(wave): ...@@ -1557,11 +1566,12 @@ def get_dx_dy_blue(wave):
# # 色散方向 # # 色散方向
# dxdl = 0.2*np.array([-9.1519, -1.00000000e-06, 3.50000000e-08, -5.00000000e-09, # 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]) # -1.70000000e-11, 4.00949787e-12, -6.16873452e-15])
##### update @2024.10.16 # update @2024.10.16
dydl=np.array([ 2447.9, -1/0.141, 0.0000075, 0.00000078, -0.000000000007] ) ; #色散方向 dydl = np.array([2447.9, -1/0.141, 0.0000075,
dxdl=np.array([ 5.46, -1.5e-02, 3.5e-08, -5.0e-09]) #垂直方向 0.00000078, -0.000000000007]) # 色散方向
dxdl = np.array([5.46, -1.5e-02, 3.5e-08, -5.0e-09]) # 垂直方向
dx = 0.0 dx = 0.0
dy = 0.0 dy = 0.0
for i in range(len(dxdl)): for i in range(len(dxdl)):
...@@ -1593,10 +1603,11 @@ def get_dx_dy_red(wave): ...@@ -1593,10 +1603,11 @@ def get_dx_dy_red(wave):
# 0.0000028, -0.0000000000007, 0.0, 0.0]) # 色散方向 # 0.0000028, -0.0000000000007, 0.0, 0.0]) # 色散方向
# dxdl = 0.00325*np.array([-1638.8, 4.0e-2, 5.500e-3, - # dxdl = 0.00325*np.array([-1638.8, 4.0e-2, 5.500e-3, -
# 5.2e-10, 1.7000e-10, 7.1e-13, -5.16e-15]) # 垂直方向 # 5.2e-10, 1.7000e-10, 7.1e-13, -5.16e-15]) # 垂直方向
## update @2014.10.17 # update @2014.10.17
dydl=np.array([3519.78622, -1/0.1555, 0.0000048, 0.00000028, -0.0000000000007] ) #色散方向 dydl = np.array([3519.78622, -1/0.1555, 0.0000048,
dxdl=np.array([-5.6305, 1.0e-2, 5.500e-7, -5.2e-10]) #垂直方向 0.00000028, -0.0000000000007]) # 色散方向
dxdl = np.array([-5.6305, 1.0e-2, 5.500e-7, -5.2e-10]) # 垂直方向
dx = 0.0 dx = 0.0
dy = 0.0 dy = 0.0
for i in range(len(dxdl)): for i in range(len(dxdl)):
...@@ -1778,7 +1789,7 @@ class IFSsimulator(): ...@@ -1778,7 +1789,7 @@ class IFSsimulator():
self.readoutNoise = self.config.getboolean( self.readoutNoise = self.config.getboolean(
self.section, 'readoutnoise') self.section, 'readoutnoise')
self.appbianpai = self.config.getboolean( self.appbianpai = self.config.getboolean(
self.section, 'appbianpai') self.section, 'appbianpai')
...@@ -2002,16 +2013,16 @@ class IFSsimulator(): ...@@ -2002,16 +2013,16 @@ class IFSsimulator():
self.processConfigs() self.processConfigs()
self._createEmpty() self._createEmpty()
self.information['dir_path']=dir_path self.information['dir_path'] = dir_path
self.information['result_path']=result_path self.information['result_path'] = result_path
self.information['debug'] = debug self.information['debug'] = debug
self.source=sourcein self.source = sourcein
self.simnumber=simnumber self.simnumber = simnumber
############################################################ ############################################################
...@@ -2019,7 +2030,6 @@ class IFSsimulator(): ...@@ -2019,7 +2030,6 @@ class IFSsimulator():
result_day = now.strftime("%Y-%m-%d") result_day = now.strftime("%Y-%m-%d")
if self.source == 'LAMP': if self.source == 'LAMP':
if applyhole == 'yes': if applyhole == 'yes':
ss = '_with_hole_' ss = '_with_hole_'
...@@ -2041,10 +2051,9 @@ class IFSsimulator(): ...@@ -2041,10 +2051,9 @@ class IFSsimulator():
# self.result_path = '../IFS_simData/'+self.source+ss+result_day # self.result_path = '../IFS_simData/'+self.source+ss+result_day
# else: # else:
# self.result_path = '/data/ifspip/CCD_ima/'+self.source+ss+result_day # self.result_path = '/data/ifspip/CCD_ima/'+self.source+ss+result_day
self.result_path = self.information['result_path'] + \
'/'+self.source+ss+result_day
self.result_path= self.information['result_path']+'/'+self.source+ss+result_day
print(self.information['result_path']) print(self.information['result_path'])
if os.path.isdir(self.result_path) == False: if os.path.isdir(self.result_path) == False:
...@@ -2063,7 +2072,7 @@ class IFSsimulator(): ...@@ -2063,7 +2072,7 @@ class IFSsimulator():
############################################################ ############################################################
self.log = lg(self.result_path+'/log_file/IFS_' + self.log = lg(self.result_path+'/log_file/IFS_' +
self.source+'_'+data_time+'_Num_'+str(simnumber)+'.log') self.source+'_'+data_time+'_Num_'+str(simnumber)+'.log')
self.log.info('STARTING A NEW SIMULATION') self.log.info('STARTING A NEW SIMULATION')
...@@ -2163,14 +2172,14 @@ class IFSsimulator(): ...@@ -2163,14 +2172,14 @@ class IFSsimulator():
slice_red['py'][i] = 50+250+randRedpos[i]*4 slice_red['py'][i] = 50+250+randRedpos[i]*4
slice_red['px'][i] = 3.55/0.015*(i-16)+1190.0+118 slice_red['px'][i] = 3.55/0.015*(i-16)+1190.0+118
####### #######
###### flip the fringe up to down,down to up@2024.10.16 # flip the fringe up to down,down to up@2024.10.16
self.slice_blue = dict() self.slice_blue = dict()
self.slice_red = dict() self.slice_red = dict()
self.slice_blue['py'] =2000-slice_blue['py'] self.slice_blue['py'] = 2000-slice_blue['py']
self.slice_blue['px'] = slice_blue['px'] self.slice_blue['px'] = slice_blue['px']
self.slice_red['py'] =3000-slice_red['py'] self.slice_red['py'] = 3000-slice_red['py']
self.slice_red['px'] =slice_red['px'] self.slice_red['px'] = slice_red['px']
####################################################################### #######################################################################
maskSlice = dict() maskSlice = dict()
...@@ -2419,57 +2428,56 @@ class IFSsimulator(): ...@@ -2419,57 +2428,56 @@ class IFSsimulator():
""" """
self.log.info('Added dark current to bule and red channel') self.log.info('Added dark current to bule and red channel')
if self.information['dark1_b'] > 0.001 or self.information['dark1_b'] > 0.001:
if self.information['dark1_b']>0.001 or self.information['dark1_b']>0.001:
self.log.error( self.log.error(
'dark1_b value error, it shoub be in [0.0001, 0.001]!') 'dark1_b value error, it shoub be in [0.0001, 0.001]!')
raise ValueError( raise ValueError(
'dark1_b value error, it shoub be in [0.0001, 0.001]!') 'dark1_b value error, it shoub be in [0.0001, 0.001]!')
if self.information['dark2_b']>0.001 or self.information['dark2_b']>0.001: if self.information['dark2_b'] > 0.001 or self.information['dark2_b'] > 0.001:
self.log.error( self.log.error(
'dark2_b value error, it shoub be in [0.0001, 0.001]!') 'dark2_b value error, it shoub be in [0.0001, 0.001]!')
raise ValueError( raise ValueError(
'dark2_b value error, it shoub be in [0.0001, 0.001]!') 'dark2_b value error, it shoub be in [0.0001, 0.001]!')
if self.information['dark3_b']>0.001 or self.information['dark3_b']>0.001: if self.information['dark3_b'] > 0.001 or self.information['dark3_b'] > 0.001:
self.log.error( self.log.error(
'dark3_b value error, it shoub be in [0.0001, 0.001]!') 'dark3_b value error, it shoub be in [0.0001, 0.001]!')
raise ValueError( raise ValueError(
'dark3_b value error, it shoub be in [0.0001, 0.001]!') 'dark3_b value error, it shoub be in [0.0001, 0.001]!')
if self.information['dark4_b']>0.001 or self.information['dark4_b']>0.001: if self.information['dark4_b'] > 0.001 or self.information['dark4_b'] > 0.001:
self.log.error( self.log.error(
'dark4_b value error, it shoub be in [0.0001, 0.001]!') 'dark4_b value error, it shoub be in [0.0001, 0.001]!')
raise ValueError( raise ValueError(
'dark4_b value error, it shoub be in [0.0001, 0.001]!') 'dark4_b value error, it shoub be in [0.0001, 0.001]!')
#################################################################### ####################################################################
if self.information['dark1_r']>0.001 or self.information['dark1_r']>0.001: if self.information['dark1_r'] > 0.001 or self.information['dark1_r'] > 0.001:
self.log.error( self.log.error(
'dark1_r value error, it shoub be in [0.0001, 0.001]!') 'dark1_r value error, it shoub be in [0.0001, 0.001]!')
raise ValueError( raise ValueError(
'dark1_r value error, it shoub be in [0.0001, 0.001]!') 'dark1_r value error, it shoub be in [0.0001, 0.001]!')
if self.information['dark2_r']>0.001 or self.information['dark2_r']>0.001: if self.information['dark2_r'] > 0.001 or self.information['dark2_r'] > 0.001:
self.log.error( self.log.error(
'dark2_r value error, it shoub be in [0.0001, 0.001]!') 'dark2_r value error, it shoub be in [0.0001, 0.001]!')
raise ValueError( raise ValueError(
'dark2_r value error, it shoub be in [0.0001, 0.001]!') 'dark2_r value error, it shoub be in [0.0001, 0.001]!')
if self.information['dark3_r']>0.001 or self.information['dark3_r']>0.001: if self.information['dark3_r'] > 0.001 or self.information['dark3_r'] > 0.001:
self.log.error( self.log.error(
'dark3_r value error, it shoub be in [0.0001, 0.001]!') 'dark3_r value error, it shoub be in [0.0001, 0.001]!')
raise ValueError( raise ValueError(
'dark3_r value error, it shoub be in [0.0001, 0.001]!') 'dark3_r value error, it shoub be in [0.0001, 0.001]!')
if self.information['dark4_r']>0.001 or self.information['dark4_r']>0.001: if self.information['dark4_r'] > 0.001 or self.information['dark4_r'] > 0.001:
self.log.error( self.log.error(
'dark4_r value error, it shoub be in [0.0001, 0.001]!') 'dark4_r value error, it shoub be in [0.0001, 0.001]!')
raise ValueError( raise ValueError(
'dark4_r value error, it shoub be in [0.0001, 0.001]!') 'dark4_r value error, it shoub be in [0.0001, 0.001]!')
########################################################################### ###########################################################################
# blue zone 4 # blue zone 4
self.image_b[0:1024, 0:2048] += self.information['exptime'] * \ self.image_b[0:1024, 0:2048] += self.information['exptime'] * \
...@@ -2690,59 +2698,57 @@ class IFSsimulator(): ...@@ -2690,59 +2698,57 @@ class IFSsimulator():
average=0.0 and std=readout noise. average=0.0 and std=readout noise.
""" """
self.log.info('readnoise added in blue channel') self.log.info('readnoise added in blue channel')
if self.information['rn1_b']>10 or self.information['rn1_b']<3: if self.information['rn1_b'] > 10 or self.information['rn1_b'] < 3:
self.log.error( self.log.error(
'rn1_b value error, it shoub be in [3, 10]!') 'rn1_b value error, it shoub be in [3, 10]!')
raise ValueError( raise ValueError(
'rn1_b value error, it shoub be in [3, 10]!') 'rn1_b value error, it shoub be in [3, 10]!')
if self.information['rn2_b']>10 or self.information['rn2_b']<3: if self.information['rn2_b'] > 10 or self.information['rn2_b'] < 3:
self.log.error( self.log.error(
'rn2_b value error, it shoub be in [3, 10]!') 'rn2_b value error, it shoub be in [3, 10]!')
raise ValueError( raise ValueError(
'rn2_b value error, it shoub be in [3, 10]!') 'rn2_b value error, it shoub be in [3, 10]!')
if self.information['rn3_b']>10 or self.information['rn3_b']<3: if self.information['rn3_b'] > 10 or self.information['rn3_b'] < 3:
self.log.error( self.log.error(
'rn3_b value error, it shoub be in [3, 10]!') 'rn3_b value error, it shoub be in [3, 10]!')
raise ValueError( raise ValueError(
'rn3_b value error, it shoub be in [3, 10]!') 'rn3_b value error, it shoub be in [3, 10]!')
if self.information['rn4_b']>10 or self.information['rn4_b']<3: if self.information['rn4_b'] > 10 or self.information['rn4_b'] < 3:
self.log.error( self.log.error(
'rn4_b value error, it shoub be in [3, 10]!') 'rn4_b value error, it shoub be in [3, 10]!')
raise ValueError( raise ValueError(
'rn4_b value error, it shoub be in [3, 10]!') 'rn4_b value error, it shoub be in [3, 10]!')
#################################################################### ####################################################################
if self.information['rn1_r']>10 or self.information['rn1_r']<3: if self.information['rn1_r'] > 10 or self.information['rn1_r'] < 3:
self.log.error( self.log.error(
'rn1_r value error, it shoub be in [3, 10]!') 'rn1_r value error, it shoub be in [3, 10]!')
raise ValueError( raise ValueError(
'rn1_r value error, it shoub be in [3, 10]!') 'rn1_r value error, it shoub be in [3, 10]!')
if self.information['rn2_r']>10 or self.information['rn2_r']<3: if self.information['rn2_r'] > 10 or self.information['rn2_r'] < 3:
self.log.error( self.log.error(
'rn2_r value error, it shoub be in [3, 10]!') 'rn2_r value error, it shoub be in [3, 10]!')
raise ValueError( raise ValueError(
'rn2_r value error, it shoub be in [3, 10]!') 'rn2_r value error, it shoub be in [3, 10]!')
if self.information['rn3_r']>10 or self.information['rn3_r']<3: if self.information['rn3_r'] > 10 or self.information['rn3_r'] < 3:
self.log.error( self.log.error(
'rn3_r value error, it shoub be in [3, 10]!') 'rn3_r value error, it shoub be in [3, 10]!')
raise ValueError( raise ValueError(
'rn3_r value error, it shoub be in [3, 10]!') 'rn3_r value error, it shoub be in [3, 10]!')
if self.information['rn4_r']>10 or self.information['rn4_r']<3: if self.information['rn4_r'] > 10 or self.information['rn4_r'] < 3:
self.log.error( self.log.error(
'rn4_r value error, it shoub be in [3, 10]!') 'rn4_r value error, it shoub be in [3, 10]!')
raise ValueError( raise ValueError(
'rn4_r value error, it shoub be in [3, 10]!') 'rn4_r value error, it shoub be in [3, 10]!')
########################################################33 # 33
# blue zone 1 # blue zone 1
np.random.seed() np.random.seed()
prescan = 50 prescan = 50
...@@ -2861,56 +2867,55 @@ class IFSsimulator(): ...@@ -2861,56 +2867,55 @@ class IFSsimulator():
""" """
Convert from electrons to ADUs using the value read from the configuration file. Convert from electrons to ADUs using the value read from the configuration file.
""" """
if self.information['gain1_b']>2 or self.information['gain1_b']<1: if self.information['gain1_b'] > 2 or self.information['gain1_b'] < 1:
self.log.error( self.log.error(
'gain1_b value error, it shoub be in [1, 2]!') 'gain1_b value error, it shoub be in [1, 2]!')
raise ValueError( raise ValueError(
'gain1_b value error, it shoub be in [1, 2]!') 'gain1_b value error, it shoub be in [1, 2]!')
if self.information['gain2_b']>2 or self.information['gain2_b']<1: if self.information['gain2_b'] > 2 or self.information['gain2_b'] < 1:
self.log.error( self.log.error(
'gain2_b value error, it shoub be in [1, 2]!') 'gain2_b value error, it shoub be in [1, 2]!')
raise ValueError( raise ValueError(
'gain2_b value error, it shoub be in [1, 2]!') 'gain2_b value error, it shoub be in [1, 2]!')
if self.information['gain3_b']>2 or self.information['gain3_b']<1: if self.information['gain3_b'] > 2 or self.information['gain3_b'] < 1:
self.log.error( self.log.error(
'gain3_b value error, it shoub be in [1, 2]!') 'gain3_b value error, it shoub be in [1, 2]!')
raise ValueError( raise ValueError(
'gain3_b value error, it shoub be in [1, 2]!') 'gain3_b value error, it shoub be in [1, 2]!')
if self.information['gain4_b']>2 or self.information['gain4_b']<1: if self.information['gain4_b'] > 2 or self.information['gain4_b'] < 1:
self.log.error( self.log.error(
'gain4_b value error, it shoub be in [1, 2]!') 'gain4_b value error, it shoub be in [1, 2]!')
raise ValueError( raise ValueError(
'gain4_b value error, it shoub be in [1, 2]!') 'gain4_b value error, it shoub be in [1, 2]!')
#################################################################### ####################################################################
if self.information['gain1_r']>2 or self.information['gain1_r']<1: if self.information['gain1_r'] > 2 or self.information['gain1_r'] < 1:
self.log.error( self.log.error(
'gain1_r value error, it shoub be in [1, 2]!') 'gain1_r value error, it shoub be in [1, 2]!')
raise ValueError( raise ValueError(
'gain1_r value error, it shoub be in [1, 2]!') 'gain1_r value error, it shoub be in [1, 2]!')
if self.information['gain2_r']>2 or self.information['gain2_r']<1: if self.information['gain2_r'] > 2 or self.information['gain2_r'] < 1:
self.log.error( self.log.error(
'gain2_r value error, it shoub be in [1, 2]!') 'gain2_r value error, it shoub be in [1, 2]!')
raise ValueError( raise ValueError(
'gain2_r value error, it shoub be in [1, 2]!') 'gain2_r value error, it shoub be in [1, 2]!')
if self.information['gain3_r']>2 or self.information['gain3_r']<1: if self.information['gain3_r'] > 2 or self.information['gain3_r'] < 1:
self.log.error( self.log.error(
'gain3_r value error, it shoub be in [1, 2]!') 'gain3_r value error, it shoub be in [1, 2]!')
raise ValueError( raise ValueError(
'gain3_r value error, it shoub be in [1, 2]!') 'gain3_r value error, it shoub be in [1, 2]!')
if self.information['gain4_r']>2 or self.information['gain4_r']<1: if self.information['gain4_r'] > 2 or self.information['gain4_r'] < 1:
self.log.error( self.log.error(
'gain4_r value error, it shoub be in [1, 2]!') 'gain4_r value error, it shoub be in [1, 2]!')
raise ValueError('gain4_r value error, it shoub be in [1, 2]!') raise ValueError('gain4_r value error, it shoub be in [1, 2]!')
#################################################################### ####################################################################
self.log.info( self.log.info(
'Converting from electrons to ADUs using a factor of gain') 'Converting from electrons to ADUs using a factor of gain')
...@@ -2936,7 +2941,7 @@ class IFSsimulator(): ...@@ -2936,7 +2941,7 @@ class IFSsimulator():
##########third part, means old zone 1 ################### ##########third part, means old zone 1 ###################
self.image_b[0:1344, 2418*2:2418*3] /= self.information['gain3_b'] self.image_b[0:1344, 2418*2:2418*3] /= self.information['gain3_b']
#### fourth part, means old zone 2 # fourth part, means old zone 2
self.image_b[0:1344, 2418*3:2418*4] /= self.information['gain4_b'] self.image_b[0:1344, 2418*3:2418*4] /= self.information['gain4_b']
############################################################################ ############################################################################
...@@ -2974,50 +2979,50 @@ class IFSsimulator(): ...@@ -2974,50 +2979,50 @@ class IFSsimulator():
""" """
if self.information['bias1_b']>2000 or self.information['bias1_b']<100: if self.information['bias1_b'] > 2000 or self.information['bias1_b'] < 100:
self.log.error( self.log.error(
'bias1_b value error, it shoub be in [100, 2000]!') 'bias1_b value error, it shoub be in [100, 2000]!')
raise ValueError( raise ValueError(
'bias1_b value error, it shoub be in [100, 2000]!') 'bias1_b value error, it shoub be in [100, 2000]!')
if self.information['bias2_b']>2000 or self.information['bias2_b']<100: if self.information['bias2_b'] > 2000 or self.information['bias2_b'] < 100:
self.log.error( self.log.error(
'bias2_b value error, it shoub be in [100, 2000]!') 'bias2_b value error, it shoub be in [100, 2000]!')
raise ValueError( raise ValueError(
'bias2_b value error, it shoub be in [100, 2000]!') 'bias2_b value error, it shoub be in [100, 2000]!')
if self.information['bias3_b']>2000 or self.information['bias3_b']<100: if self.information['bias3_b'] > 2000 or self.information['bias3_b'] < 100:
self.log.error( self.log.error(
'bias3_b value error, it shoub be in [100, 2000]!') 'bias3_b value error, it shoub be in [100, 2000]!')
raise ValueError( raise ValueError(
'bias3_b value error, it shoub be in [100, 2000]!') 'bias3_b value error, it shoub be in [100, 2000]!')
if self.information['bias4_b']>2000 or self.information['bias4_b']<100: if self.information['bias4_b'] > 2000 or self.information['bias4_b'] < 100:
self.log.error( self.log.error(
'bias4_b value error, it shoub be in [100, 2000]!') 'bias4_b value error, it shoub be in [100, 2000]!')
raise ValueError( raise ValueError(
'bias4_b value error, it shoub be in [100, 2000]!') 'bias4_b value error, it shoub be in [100, 2000]!')
#################################################################### ####################################################################
if self.information['bias1_r']>2000 or self.information['bias1_r']<100: if self.information['bias1_r'] > 2000 or self.information['bias1_r'] < 100:
self.log.error( self.log.error(
'bias1_r value error, it shoub be in [100, 2000]!') 'bias1_r value error, it shoub be in [100, 2000]!')
raise ValueError( raise ValueError(
'bias1_r value error, it shoub be in [100, 2000]!') 'bias1_r value error, it shoub be in [100, 2000]!')
if self.information['bias2_r']>2000 or self.information['bias2_r']<100: if self.information['bias2_r'] > 2000 or self.information['bias2_r'] < 100:
self.log.error( self.log.error(
'bias2_r value error, it shoub be in [100, 2000]!') 'bias2_r value error, it shoub be in [100, 2000]!')
raise ValueError( raise ValueError(
'bias2_r value error, it shoub be in [100, 2000]!') 'bias2_r value error, it shoub be in [100, 2000]!')
if self.information['bias3_r']>2000 or self.information['bias3_r']<100: if self.information['bias3_r'] > 2000 or self.information['bias3_r'] < 100:
self.log.error( self.log.error(
'bias3_r value error, it shoub be in [100, 2000]!') 'bias3_r value error, it shoub be in [100, 2000]!')
raise ValueError( raise ValueError(
'bias3_r value error, it shoub be in [100, 2000]!') 'bias3_r value error, it shoub be in [100, 2000]!')
if self.information['bias4_r']>2000 or self.information['bias4_r']<100: if self.information['bias4_r'] > 2000 or self.information['bias4_r'] < 100:
self.log.error( self.log.error(
'bias4_r value error, it shoub be in [100, 2000]!') 'bias4_r value error, it shoub be in [100, 2000]!')
raise ValueError( raise ValueError(
...@@ -3264,7 +3269,6 @@ class IFSsimulator(): ...@@ -3264,7 +3269,6 @@ class IFSsimulator():
# def applyImageShift(self): # def applyImageShift(self):
# """ # """
# Returns # Returns
# ------- # -------
# None. # None.
...@@ -3367,8 +3371,8 @@ class IFSsimulator(): ...@@ -3367,8 +3371,8 @@ class IFSsimulator():
# y1 = 2418*3+prescan # y1 = 2418*3+prescan
# y2 = y1+2048 # y2 = y1+2048
# temp[x1:x2, y1:y2] = np.flipud(np.fliplr(imgb[1024:2048, 2048:4096])) # temp[x1:x2, y1:y2] = np.flipud(np.fliplr(imgb[1024:2048, 2048:4096]))
### update 2024.10.18 # update 2024.10.18
# first part, old OSG part ,## shift: left to right # first part, old OSG part ,## shift: left to right
x1 = 0 x1 = 0
x2 = x1+1024 x2 = x1+1024
...@@ -3399,7 +3403,6 @@ class IFSsimulator(): ...@@ -3399,7 +3403,6 @@ class IFSsimulator():
y1 = 2418*3+prescan y1 = 2418*3+prescan
y2 = y1+2048 y2 = y1+2048
temp[x1:x2, y1:y2] = np.flipud(np.fliplr(imgb[1024:2048, 2048:4096])) temp[x1:x2, y1:y2] = np.flipud(np.fliplr(imgb[1024:2048, 2048:4096]))
self.image_b = temp self.image_b = temp
...@@ -3437,15 +3440,15 @@ class IFSsimulator(): ...@@ -3437,15 +3440,15 @@ class IFSsimulator():
# y1 = 3442*3+prescan # y1 = 3442*3+prescan
# y2 = y1+3072 # y2 = y1+3072
# temp[x1:x2, y1:y2] = np.flipud(np.fliplr(imgr[1536:3072, 3072:6144])) # temp[x1:x2, y1:y2] = np.flipud(np.fliplr(imgr[1536:3072, 3072:6144]))
###### update @2024.10.18 # update @2024.10.18
# readout image ,first part, old OSG, shift: left to right # readout image ,first part, old OSG, shift: left to right
x1 = 0 x1 = 0
x2 = x1+1536 x2 = x1+1536
y1 = 0+prescan y1 = 0+prescan
y2 = y1+3072 y2 = y1+3072
temp[x1:x2, y1:y2] = np.fliplr(imgr[0:1536, 3072:6144]) temp[x1:x2, y1:y2] = np.fliplr(imgr[0:1536, 3072:6144])
########## second part, old OSH ,no change; ################# ########## second part, old OSH ,no change; #################
# np.fliplr(b2) ## left to right # np.fliplr(b2) ## left to right
# np.flipud(b3) ## down to up # np.flipud(b3) ## down to up
...@@ -3469,7 +3472,6 @@ class IFSsimulator(): ...@@ -3469,7 +3472,6 @@ class IFSsimulator():
y1 = 3442*3+prescan y1 = 3442*3+prescan
y2 = y1+3072 y2 = y1+3072
temp[x1:x2, y1:y2] = np.flipud(np.fliplr(imgr[1536:3072, 3072:6144])) temp[x1:x2, y1:y2] = np.flipud(np.fliplr(imgr[1536:3072, 3072:6144]))
self.image_r = temp self.image_r = temp
...@@ -3496,17 +3498,15 @@ class IFSsimulator(): ...@@ -3496,17 +3498,15 @@ class IFSsimulator():
Updates header with the input values and flags used during simulation. Updates header with the input values and flags used during simulation.
""" """
###### ######
###### add random pointing error to telescope, the telescope poingt parameter in # add random pointing error to telescope, the telescope poingt parameter in
###### fits header have random pointing error # fits header have random pointing error
ud_ra = np.random.random() # Choose a random shift in arcsec ud_ra = np.random.random() # Choose a random shift in arcsec
self.information['ra_pnt0']=self.information['dec_pnt0']+0.01*ud_ra self.information['ra_pnt0'] = self.information['dec_pnt0']+0.01*ud_ra
ud_dec= np.random.random()
self.information['dec_pnt0']=self.information['dec_pnt0']+0.01*ud_dec
ud_dec = np.random.random()
self.information['dec_pnt0'] = self.information['dec_pnt0']+0.01*ud_dec
HeaderTest = 'no' HeaderTest = 'no'
sim_ver = str(self.information['sim_ver']) sim_ver = str(self.information['sim_ver'])
...@@ -4677,7 +4677,6 @@ class IFSsimulator(): ...@@ -4677,7 +4677,6 @@ class IFSsimulator():
####################################################### #######################################################
# CCD quantum efficiency # CCD quantum efficiency
CCD_Qe_lam = np.interp(lam, self.CCD_Qe[:, 0], self.CCD_Qe[:, 1]) CCD_Qe_lam = np.interp(lam, self.CCD_Qe[:, 0], self.CCD_Qe[:, 1])
...@@ -4769,7 +4768,6 @@ class IFSsimulator(): ...@@ -4769,7 +4768,6 @@ class IFSsimulator():
photons_blue.addTo(blue_img) photons_blue.addTo(blue_img)
# fits.writeto('blueImg.fits',blue_img.array,overwrite=True) # fits.writeto('blueImg.fits',blue_img.array,overwrite=True)
if (lam >= 560.0) & (lam <= 1000.0): if (lam >= 560.0) & (lam <= 1000.0):
...@@ -4885,32 +4883,35 @@ class IFSsimulator(): ...@@ -4885,32 +4883,35 @@ class IFSsimulator():
self.information['ra_obj'] = a[0].header['RA'] # in degree self.information['ra_obj'] = a[0].header['RA'] # in degree
self.information['dec_obj'] = a[0].header['DEC'] # in degree self.information['dec_obj'] = a[0].header['DEC'] # in degree
disRa = (telRa-skyRa)/3600.0 # convert unit of arcsec to degree disRa = (telRa-skyRa)/3600.0 # convert unit of arcsec to degree
disDec = (telDec-skyDec)/3600.0 # convert unit of arcsec to degree disDec = (telDec-skyDec)/3600.0 # convert unit of arcsec to degree
# RA_PNT0 =OBJ_RA + DIS_RA/cos(OBJ_DEC), DEC_PNT0= OBJ_DEC DIS_DEC # RA_PNT0 =OBJ_RA + DIS_RA/cos(OBJ_DEC), DEC_PNT0= OBJ_DEC DIS_DEC
self.information['ra_pnt0'] = a[0].header['RA'] + \ self.information['ra_pnt0'] = a[0].header['RA'] + \
disRa/np.cos(a[0].header['DEC']/180.0*np.pi) disRa/np.cos(a[0].header['DEC']/180.0*np.pi)
self.information['dec_pnt0'] = a[0].header['DEC'] + disDec self.information['dec_pnt0'] = a[0].header['DEC'] + disDec
############### calculate the earthshine and zodiacal noise ,new code 2023.11.1 ############ ############### calculate the earthshine and zodiacal noise ,new code 2023.11.1 ############
############### ###############
self.log.info('Real telescope pointing in Ra = %f, Dec = %f' % (self.information['ra_pnt0'], self.information['dec_pnt0'])) self.log.info('Real telescope pointing in Ra = %f, Dec = %f' % (
self.information['ra_pnt0'], self.information['dec_pnt0']))
ra = self.information['ra_pnt0'] ra = self.information['ra_pnt0']
dec = self.information['dec_pnt0'] dec = self.information['dec_pnt0']
time_jd = time2jd(self.dt) time_jd = time2jd(self.dt)
if self.appbianpai: if self.appbianpai:
sn = self.simnumber-1
sn=self.simnumber-1; x_sat = float(self.bianpai_data['x_sat']
x_sat = float(self.bianpai_data['x_sat'][sn*5+self.exptime_start_index]) [sn*5+self.exptime_start_index])
y_sat = float(self.bianpai_data['y_sat'][sn*5+self.exptime_start_index]) y_sat = float(self.bianpai_data['y_sat']
z_sat = float(self.bianpai_data['z_sat'][sn*5+self.exptime_start_index]) [sn*5+self.exptime_start_index])
z_sat = float(self.bianpai_data['z_sat']
[sn*5+self.exptime_start_index])
### ###
else: else:
...@@ -4996,18 +4997,24 @@ class IFSsimulator(): ...@@ -4996,18 +4997,24 @@ class IFSsimulator():
################################ ################################
############## doppler effect to photons.wavelength ############# ############## doppler effect to photons.wavelength #############
if self.appbianpai: if self.appbianpai:
sn=self.simnumber-1; sn = self.simnumber-1
x_sat = float(self.bianpai_data['x_sat'][sn*5+self.exptime_start_index]) x_sat = float(self.bianpai_data['x_sat']
y_sat = float(self.bianpai_data['y_sat'][sn*5+self.exptime_start_index]) [sn*5+self.exptime_start_index])
z_sat = float(self.bianpai_data['z_sat'][sn*5+self.exptime_start_index]) y_sat = float(self.bianpai_data['y_sat']
[sn*5+self.exptime_start_index])
vx_sat = float(self.bianpai_data['vx_sat'][sn*5+self.exptime_start_index]) z_sat = float(self.bianpai_data['z_sat']
vy_sat = float(self.bianpai_data['vy_sat'][sn*5+self.exptime_start_index]) [sn*5+self.exptime_start_index])
vz_sat = float(self.bianpai_data['vz_sat'][sn*5+self.exptime_start_index])
vx_sat = float(
self.bianpai_data['vx_sat'][sn*5+self.exptime_start_index])
vy_sat = float(
self.bianpai_data['vy_sat'][sn*5+self.exptime_start_index])
vz_sat = float(
self.bianpai_data['vz_sat'][sn*5+self.exptime_start_index])
else: else:
# self.orbit_pars # self.orbit_pars
x_sat = float(self.orbit_pars[self.orbit_exp_num, 1]) x_sat = float(self.orbit_pars[self.orbit_exp_num, 1])
y_sat = float(self.orbit_pars[self.orbit_exp_num, 2]) y_sat = float(self.orbit_pars[self.orbit_exp_num, 2])
...@@ -5033,35 +5040,42 @@ class IFSsimulator(): ...@@ -5033,35 +5040,42 @@ class IFSsimulator():
################################################# #################################################
if self.appbianpai: if self.appbianpai:
sn=self.simnumber-1; sn = self.simnumber-1
p1x = float(self.bianpai_data['x_sat'][sn*5+self.exptime_end_index]) p1x = float(self.bianpai_data['x_sat']
p1y = float(self.bianpai_data['y_sat'][sn*5+self.exptime_end_index]) [sn*5+self.exptime_end_index])
p1z = float(self.bianpai_data['z_sat'][sn*5+self.exptime_end_index]) p1y = float(self.bianpai_data['y_sat']
[sn*5+self.exptime_end_index])
p1vx = float(self.bianpai_data['vx_sat'][sn*5+self.exptime_end_index]) p1z = float(self.bianpai_data['z_sat']
p1vy = float(self.bianpai_data['vy_sat'][sn*5+self.exptime_end_index]) [sn*5+self.exptime_end_index])
p1vz = float(self.bianpai_data['vz_sat'][sn*5+self.exptime_end_index])
p1vx = float(self.bianpai_data['vx_sat']
[sn*5+self.exptime_end_index])
p1vy = float(self.bianpai_data['vy_sat']
[sn*5+self.exptime_end_index])
p1vz = float(self.bianpai_data['vz_sat']
[sn*5+self.exptime_end_index])
else: else:
# exposure end time is t2 ; # exposure end time is t2 ;
t2 = self.dt+timedelta(seconds=self.information['exptime']) t2 = self.dt+timedelta(seconds=self.information['exptime'])
### data read time is the exposure end time plus readouttime ### ### data read time is the exposure end time plus readouttime ###
# data read end time is t3 # data read end time is t3
t2jd = time2jd(t2) t2jd = time2jd(t2)
if self.orbit_pars[-1, 0] < t2jd: # orbit parameters are not in currenct txt file # orbit parameters are not in currenct txt file
if self.orbit_pars[-1, 0] < t2jd:
self.orbit_file_num = self.orbit_file_num+1 self.orbit_file_num = self.orbit_file_num+1
fn = self.information['dir_path'] + \ fn = self.information['dir_path'] + \
'IFS_inputdata/TianCe/orbit20160925/' + \ 'IFS_inputdata/TianCe/orbit20160925/' + \
str(self.orbit_file_num)+'.txt' str(self.orbit_file_num)+'.txt'
self.orbit_pars = np.loadtxt(fn) self.orbit_pars = np.loadtxt(fn)
self.orbit_exp_num = 0 self.orbit_exp_num = 0
for k in range(self.orbit_exp_num, len(self.orbit_pars), 1): for k in range(self.orbit_exp_num, len(self.orbit_pars), 1):
if t2jd-self.orbit_pars[k, 0] < 0: if t2jd-self.orbit_pars[k, 0] < 0:
break break
if k == 0: if k == 0:
...@@ -5072,14 +5086,14 @@ class IFSsimulator(): ...@@ -5072,14 +5086,14 @@ class IFSsimulator():
self.orbit_pars[k, 2])*deltaT.seconds/120 self.orbit_pars[k, 2])*deltaT.seconds/120
p1z = self.orbit_pars[k, 3]-(self.orbit_pars[k+1, 3] - p1z = self.orbit_pars[k, 3]-(self.orbit_pars[k+1, 3] -
self.orbit_pars[k, 3])*deltaT.seconds/120 self.orbit_pars[k, 3])*deltaT.seconds/120
p1vx = self.orbit_pars[k, 4]-(self.orbit_pars[k+1, 4] - p1vx = self.orbit_pars[k, 4]-(self.orbit_pars[k+1, 4] -
self.orbit_pars[k, 4])*deltaT.seconds/120 self.orbit_pars[k, 4])*deltaT.seconds/120
p1vy = self.orbit_pars[k, 5]-(self.orbit_pars[k+1, 5] - p1vy = self.orbit_pars[k, 5]-(self.orbit_pars[k+1, 5] -
self.orbit_pars[k, 5])*deltaT.seconds/120 self.orbit_pars[k, 5])*deltaT.seconds/120
p1vz = self.orbit_pars[k, 6]-(self.orbit_pars[k+1, 6] - p1vz = self.orbit_pars[k, 6]-(self.orbit_pars[k+1, 6] -
self.orbit_pars[k, 6])*deltaT.seconds/120 self.orbit_pars[k, 6])*deltaT.seconds/120
else: else:
deltaT = jd2time(self.orbit_pars[k, 0])-t2 deltaT = jd2time(self.orbit_pars[k, 0])-t2
p1x = self.orbit_pars[k-1, 1]+(self.orbit_pars[k, 1] - p1x = self.orbit_pars[k-1, 1]+(self.orbit_pars[k, 1] -
...@@ -5088,7 +5102,7 @@ class IFSsimulator(): ...@@ -5088,7 +5102,7 @@ class IFSsimulator():
self.orbit_pars[k-1, 2])*deltaT.seconds/120 self.orbit_pars[k-1, 2])*deltaT.seconds/120
p1z = self.orbit_pars[k-1, 3]+(self.orbit_pars[k, 3] - p1z = self.orbit_pars[k-1, 3]+(self.orbit_pars[k, 3] -
self.orbit_pars[k-1, 3])*deltaT.seconds/120 self.orbit_pars[k-1, 3])*deltaT.seconds/120
p1vx = self.orbit_pars[k-1, 4]+(self.orbit_pars[k, 4] - p1vx = self.orbit_pars[k-1, 4]+(self.orbit_pars[k, 4] -
self.orbit_pars[k-1, 4])*deltaT.seconds/120 self.orbit_pars[k-1, 4])*deltaT.seconds/120
p1vy = self.orbit_pars[k-1, 5]+(self.orbit_pars[k, 5] - p1vy = self.orbit_pars[k-1, 5]+(self.orbit_pars[k, 5] -
...@@ -5394,7 +5408,6 @@ class IFSsimulator(): ...@@ -5394,7 +5408,6 @@ class IFSsimulator():
energy_blue = energy_blue+sum(photons_blue.flux) energy_blue = energy_blue+sum(photons_blue.flux)
################# #################
# fits.writeto('blueImg.fits',blue_img.array,overwrite=True) # fits.writeto('blueImg.fits',blue_img.array,overwrite=True)
...@@ -5851,15 +5864,12 @@ class IFSsimulator(): ...@@ -5851,15 +5864,12 @@ class IFSsimulator():
self.source = sourcein self.source = sourcein
self.simnumber = simnumber self.simnumber = simnumber
#self.configure(simnumber,dir_path,result_path) # print the configfile name and path; # self.configure(simnumber,dir_path,result_path) # print the configfile name and path;
self.debug = self.information['debug'] self.debug = self.information['debug']
if self.information['exptime'] > 2000 or self.information['exptime'] < 0:
if self.information['exptime']>2000 or self.information['exptime']<0:
self.log.error( self.log.error(
'exptime value error, it shoub be in [0, 2000]!') 'exptime value error, it shoub be in [0, 2000]!')
exit(2) exit(2)
...@@ -5882,33 +5892,34 @@ class IFSsimulator(): ...@@ -5882,33 +5892,34 @@ class IFSsimulator():
self.information['sky_fitsin'] self.information['sky_fitsin']
print('self.skyfilepath = ', self.skyfilepath) print('self.skyfilepath = ', self.skyfilepath)
################################################################### ###################################################################
if self.appbianpai: if self.appbianpai:
### load yunxingbianpai csv file # load yunxingbianpai csv file
############ load star data catlog ##################### ############ load star data catlog #####################
starcat=self.information['bianpai_file'] starcat = self.information['bianpai_file']
##starcat='selection_20230517_concat.fits' # starcat='selection_20230517_concat.fits'
################################################### ###################################################
self.log.info('Stat catlog file name is %s' % (starcat)) self.log.info('Stat catlog file name is %s' % (starcat))
########################################## ##########################################
df=pd.read_csv(self.information['dir_path']+'IFS_inputdata/TianCe/'+starcat) df = pd.read_csv(
self.information['dir_path']+'IFS_inputdata/TianCe/'+starcat)
################################################################### ###################################################################
sn=self.simnumber-1; sn = self.simnumber-1
arr=np.array(df['time'][sn*5:sn*5+5]); arr = np.array(df['time'][sn*5:sn*5+5])
self.exptime_start_jd,self.exptime_start_index=find_min(arr); self.exptime_start_jd, self.exptime_start_index = find_min(arr)
self.exptime_end_jd, self.exptime_end_index =find_max(arr); self.exptime_end_jd, self.exptime_end_index = find_max(arr)
###self.earthshine_theta=df['earth_angle'][sn*5+index] # in degree # self.earthshine_theta=df['earth_angle'][sn*5+index] # in degree
self.dt = jd2time(self.exptime_start_jd); self.dt = jd2time(self.exptime_start_jd)
self.bianpai_data=df; self.bianpai_data = df
################################################################## ##################################################################
else: else:
#### load orbit parameters ##### #### load orbit parameters #####
flag = 0 flag = 0
self.dt = datetime.utcnow() self.dt = datetime.utcnow()
...@@ -5916,40 +5927,38 @@ class IFSsimulator(): ...@@ -5916,40 +5927,38 @@ class IFSsimulator():
self.simnumber*(self.information['exptime']+self.information['readouttime']+125)/120) self.simnumber*(self.information['exptime']+self.information['readouttime']+125)/120)
now_dt = datetime.utcnow() now_dt = datetime.utcnow()
now_jd = time2jd(now_dt) now_jd = time2jd(now_dt)
for k in range(1, 50, 1): for k in range(1, 50, 1):
# fn=father_path+'/IFS_inputdata/TianCe/orbit20160925/'+str(k)+'.txt'; # fn=father_path+'/IFS_inputdata/TianCe/orbit20160925/'+str(k)+'.txt';
fn = self.information['dir_path'] + \ fn = self.information['dir_path'] + \
'IFS_inputdata/TianCe/orbit20160925/'+str(k)+'.txt' 'IFS_inputdata/TianCe/orbit20160925/'+str(k)+'.txt'
d = np.loadtxt(fn) d = np.loadtxt(fn)
for kk in range(len(d[:, 0])): for kk in range(len(d[:, 0])):
if now_jd-d[kk, 0] <= 0: if now_jd-d[kk, 0] <= 0:
flag = 1 flag = 1
break break
if flag == 1: if flag == 1:
break break
################################################################## ##################################################################
if kk + self.dt_num < len(d[:, 0]): if kk + self.dt_num < len(d[:, 0]):
exptime_start_jd = d[kk+self.dt_num, 0] exptime_start_jd = d[kk+self.dt_num, 0]
self.orbit_pars = d self.orbit_pars = d
self.orbit_file_num = k self.orbit_file_num = k
self.orbit_exp_num = kk self.orbit_exp_num = kk
else: else:
fn = self.information['dir_path'] + \ fn = self.information['dir_path'] + \
'IFS_inputdata/TianCe/orbit20160925/'+str(k+1)+'.txt' 'IFS_inputdata/TianCe/orbit20160925/'+str(k+1)+'.txt'
d = np.loadtxt(fn) d = np.loadtxt(fn)
self.orbit_pars = d self.orbit_pars = d
self.orbit_file_num = k+1 self.orbit_file_num = k+1
self.orbit_exp_num = self.dt_num self.orbit_exp_num = self.dt_num
...@@ -5958,7 +5967,6 @@ class IFSsimulator(): ...@@ -5958,7 +5967,6 @@ class IFSsimulator():
################################################################## ##################################################################
################################################################## ##################################################################
self.TianCe_day = self.dt.strftime("%Y-%m-%d") self.TianCe_day = self.dt.strftime("%Y-%m-%d")
self.TianCe_exp_start = dt2hmd(self.dt) self.TianCe_exp_start = dt2hmd(self.dt)
...@@ -5975,7 +5983,7 @@ class IFSsimulator(): ...@@ -5975,7 +5983,7 @@ class IFSsimulator():
sky_rot = 2 * (ud-0.5) * 5 sky_rot = 2 * (ud-0.5) * 5
dsmax = np.floor(50*np.cos(sky_rot/180*np.pi))-34 dsmax = np.floor(50*np.cos(sky_rot/180*np.pi))-34
np.random.seed(11*self.simnumber) np.random.seed(11*self.simnumber)
ud = np.random.random() # Choose a random shift in arcsec ud = np.random.random() # Choose a random shift in arcsec
telRa = 2 * (ud-0.5) * dsmax * 0.1 # in arcsec telRa = 2 * (ud-0.5) * dsmax * 0.1 # in arcsec
...@@ -6181,9 +6189,9 @@ def runIFSsim(sourcein, configfile, dir_path, result_path, iLoop, debug, applyho ...@@ -6181,9 +6189,9 @@ def runIFSsim(sourcein, configfile, dir_path, result_path, iLoop, debug, applyho
simulate = dict() simulate = dict()
simulate[iLoop] = IFSsimulator(configfile) simulate[iLoop] = IFSsimulator(configfile)
simulate[iLoop].configure(sourcein, dir_path, result_path, iLoop, debug, applyhole) # load the configfile; simulate[iLoop].configure(sourcein, dir_path, result_path,
iLoop, debug, applyhole) # load the configfile;
if applyhole == 'yes' and sourcein == 'LAMP': if applyhole == 'yes' and sourcein == 'LAMP':
simulate[iLoop].information['holemask'] = 'yes' simulate[iLoop].information['holemask'] = 'yes'
...@@ -6194,9 +6202,9 @@ def runIFSsim(sourcein, configfile, dir_path, result_path, iLoop, debug, applyho ...@@ -6194,9 +6202,9 @@ def runIFSsim(sourcein, configfile, dir_path, result_path, iLoop, debug, applyho
# simulate[iLoop].information['dir_path'] = dir_path # simulate[iLoop].information['dir_path'] = dir_path
# simulate[iLoop].information['debug'] = debug # simulate[iLoop].information['debug'] = debug
# simulate[iLoop].information['result_path'] = result_path # simulate[iLoop].information['result_path'] = result_path
simulate[iLoop].simulate(sourcein, iLoop) simulate[iLoop].simulate(sourcein, iLoop)
return 1 return 1
......
tests/test_ifs_sim.py
View file @ 5070c190
...@@ -15,7 +15,7 @@ from csst_ifs_sim import csst_ifs_sim ...@@ -15,7 +15,7 @@ from csst_ifs_sim import csst_ifs_sim
import sys import sys
### sys.path.append('IFS_git/csst_ifs_sim/csst_ifs_sim') # sys.path.append('IFS_git/csst_ifs_sim/csst_ifs_sim')
class TestDemoFunction(unittest.TestCase): class TestDemoFunction(unittest.TestCase):
...@@ -35,27 +35,28 @@ class TestDemoFunction(unittest.TestCase): ...@@ -35,27 +35,28 @@ class TestDemoFunction(unittest.TestCase):
This case aims to test whether the demo function returns `1` if input is `None`. This case aims to test whether the demo function returns `1` if input is `None`.
""" """
# demo function test # demo function test
dir_path = os.path.join(os.environ['UNIT_TEST_DATA_ROOT'],'ifs_sim/') dir_path = os.path.join(os.environ['UNIT_TEST_DATA_ROOT'], 'ifs_sim/')
print(dir_path) print(dir_path)
print(sys.version ) print(sys.version)
###configfile = dir_path+'IFS_inputdata/configData/IFS_sim_C90.config' ###configfile = dir_path+'IFS_inputdata/configData/IFS_sim_C90.config'
configfile = './csst_ifs_sim/ifs_data/IFS_sim_C90.config' configfile = './csst_ifs_sim/ifs_data/IFS_sim_C90.config'
sourcein = 'SCI' sourcein = 'SCI'
print(configfile) print(configfile)
debug=True debug = True
result_path=dir_path+'ifs_sim_result' result_path = dir_path+'ifs_sim_result'
csst_ifs_sim.runIFSsim(sourcein, configfile, dir_path, result_path, 1, debug,'no') csst_ifs_sim.runIFSsim(sourcein, configfile,
dir_path, result_path, 1, debug, 'no')
self.assertEqual( self.assertEqual(
1 , 1, 1, 1,
"case 1: SCI sim passes.", "case 1: SCI sim passes.",
) )
############################################################## ##############################################################
def test_ifs_sim_2(self): def test_ifs_sim_2(self):
""" """
...@@ -73,26 +74,27 @@ class TestDemoFunction(unittest.TestCase): ...@@ -73,26 +74,27 @@ class TestDemoFunction(unittest.TestCase):
This case aims to test whether the demo function returns `1` if input is `None`. This case aims to test whether the demo function returns `1` if input is `None`.
""" """
# demo function test # demo function test
dir_path = os.path.join(os.environ['UNIT_TEST_DATA_ROOT'],'ifs_sim/') dir_path = os.path.join(os.environ['UNIT_TEST_DATA_ROOT'], 'ifs_sim/')
print(dir_path) print(dir_path)
print(sys.version ) print(sys.version)
###configfile = dir_path+'IFS_inputdata/configData/IFS_sim_C90.config' ###configfile = dir_path+'IFS_inputdata/configData/IFS_sim_C90.config'
configfile = './csst_ifs_sim/ifs_data/IFS_sim_C90.config' configfile = './csst_ifs_sim/ifs_data/IFS_sim_C90.config'
sourcein = 'BIAS' sourcein = 'BIAS'
print(configfile) print(configfile)
debug=True debug = True
result_path=dir_path+'ifs_sim_result' result_path = dir_path+'ifs_sim_result'
csst_ifs_sim.runIFSsim(sourcein, configfile, dir_path, result_path, 1, debug,'no') csst_ifs_sim.runIFSsim(sourcein, configfile,
dir_path, result_path, 1, debug, 'no')
self.assertEqual( self.assertEqual(
1 , 1, 1, 1,
"case 2: sim passes.", "case 2: sim passes.",
) )
################################################################### ###################################################################
def test_ifs_sim_3(self): def test_ifs_sim_3(self):
""" """
...@@ -110,26 +112,27 @@ class TestDemoFunction(unittest.TestCase): ...@@ -110,26 +112,27 @@ class TestDemoFunction(unittest.TestCase):
This case aims to test whether the demo function returns `1` if input is `None`. This case aims to test whether the demo function returns `1` if input is `None`.
""" """
# demo function test # demo function test
dir_path = os.path.join(os.environ['UNIT_TEST_DATA_ROOT'],'ifs_sim/') dir_path = os.path.join(os.environ['UNIT_TEST_DATA_ROOT'], 'ifs_sim/')
print(dir_path) print(dir_path)
print(sys.version ) print(sys.version)
###configfile = dir_path+'IFS_inputdata/configData/IFS_sim_C90.config' ###configfile = dir_path+'IFS_inputdata/configData/IFS_sim_C90.config'
configfile = './csst_ifs_sim/ifs_data/IFS_sim_C90.config' configfile = './csst_ifs_sim/ifs_data/IFS_sim_C90.config'
sourcein = 'DARK' sourcein = 'DARK'
print(configfile) print(configfile)
debug=True debug = True
result_path=dir_path+'ifs_sim_result' result_path = dir_path+'ifs_sim_result'
csst_ifs_sim.runIFSsim(sourcein, configfile, dir_path, result_path, 1, debug,'no') csst_ifs_sim.runIFSsim(sourcein, configfile,
dir_path, result_path, 1, debug, 'no')
self.assertEqual( self.assertEqual(
1 , 1, 1, 1,
"case 3: sim passes.", "case 3: sim passes.",
) )
################################################################### ###################################################################
def test_ifs_sim_4(self): def test_ifs_sim_4(self):
""" """
...@@ -147,27 +150,28 @@ class TestDemoFunction(unittest.TestCase): ...@@ -147,27 +150,28 @@ class TestDemoFunction(unittest.TestCase):
This case aims to test whether the demo function returns `1` if input is `None`. This case aims to test whether the demo function returns `1` if input is `None`.
""" """
# demo function test # demo function test
dir_path = os.path.join(os.environ['UNIT_TEST_DATA_ROOT'],'ifs_sim/') dir_path = os.path.join(os.environ['UNIT_TEST_DATA_ROOT'], 'ifs_sim/')
print(dir_path) print(dir_path)
print(sys.version ) print(sys.version)
###configfile = dir_path+'IFS_inputdata/configData/IFS_sim_C90.config' ###configfile = dir_path+'IFS_inputdata/configData/IFS_sim_C90.config'
configfile = './csst_ifs_sim/ifs_data/IFS_sim_C90.config' configfile = './csst_ifs_sim/ifs_data/IFS_sim_C90.config'
sourcein = 'LAMP' sourcein = 'LAMP'
print(configfile) print(configfile)
debug=True debug = True
result_path=dir_path+'ifs_sim_result' result_path = dir_path+'ifs_sim_result'
csst_ifs_sim.runIFSsim(sourcein, configfile, dir_path, result_path, 1, debug,'no') csst_ifs_sim.runIFSsim(sourcein, configfile,
dir_path, result_path, 1, debug, 'no')
self.assertEqual( self.assertEqual(
1 , 1, 1, 1,
"case 4: sim passes.", "case 4: sim passes.",
) )
################################################################### ###################################################################
def test_ifs_sim_5(self): def test_ifs_sim_5(self):
""" """
...@@ -185,25 +189,26 @@ class TestDemoFunction(unittest.TestCase): ...@@ -185,25 +189,26 @@ class TestDemoFunction(unittest.TestCase):
This case aims to test whether the demo function returns `1` if input is `None`. This case aims to test whether the demo function returns `1` if input is `None`.
""" """
# demo function test # demo function test
dir_path = os.path.join(os.environ['UNIT_TEST_DATA_ROOT'],'ifs_sim/') dir_path = os.path.join(os.environ['UNIT_TEST_DATA_ROOT'], 'ifs_sim/')
print(dir_path) print(dir_path)
print(sys.version ) print(sys.version)
###configfile = dir_path+'IFS_inputdata/configData/IFS_sim_C90.config' ###configfile = dir_path+'IFS_inputdata/configData/IFS_sim_C90.config'
configfile = './csst_ifs_sim/ifs_data/IFS_sim_C90.config' configfile = './csst_ifs_sim/ifs_data/IFS_sim_C90.config'
sourcein = 'LAMP' sourcein = 'LAMP'
print(configfile) print(configfile)
debug=True debug = True
result_path=dir_path+'ifs_sim_result' result_path = dir_path+'ifs_sim_result'
csst_ifs_sim.runIFSsim(sourcein, configfile, dir_path, result_path, 1, debug,'yes') csst_ifs_sim.runIFSsim(sourcein, configfile,
dir_path, result_path, 1, debug, 'yes')
self.assertEqual( self.assertEqual(
1 , 1, 1, 1,
"case 5: sim passes.", "case 5: sim passes.",
) )
################################################################### ###################################################################
def test_ifs_sim_6(self): def test_ifs_sim_6(self):
""" """
...@@ -221,24 +226,23 @@ class TestDemoFunction(unittest.TestCase): ...@@ -221,24 +226,23 @@ class TestDemoFunction(unittest.TestCase):
This case aims to test whether the demo function returns `1` if input is `None`. This case aims to test whether the demo function returns `1` if input is `None`.
""" """
# demo function test # demo function test
dir_path = os.path.join(os.environ['UNIT_TEST_DATA_ROOT'],'ifs_sim/') dir_path = os.path.join(os.environ['UNIT_TEST_DATA_ROOT'], 'ifs_sim/')
print(dir_path) print(dir_path)
print(sys.version ) print(sys.version)
###configfile = dir_path+'IFS_inputdata/configData/IFS_sim_C90.config' ###configfile = dir_path+'IFS_inputdata/configData/IFS_sim_C90.config'
configfile = './csst_ifs_sim/ifs_data/IFS_sim_C90.config' configfile = './csst_ifs_sim/ifs_data/IFS_sim_C90.config'
sourcein = 'FLAT' sourcein = 'FLAT'
print(configfile) print(configfile)
debug=True debug = True
result_path=dir_path+'ifs_sim_result' result_path = dir_path+'ifs_sim_result'
csst_ifs_sim.runIFSsim(sourcein, configfile, dir_path, result_path, 1, debug, 'no') csst_ifs_sim.runIFSsim(sourcein, configfile,
dir_path, result_path, 1, debug, 'no')
self.assertEqual( self.assertEqual(
1 , 1, 1, 1,
"case 6: sim passes.", "case 6: sim passes.",
) )
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment