Skip to content
GitLab
Commit 8165a134 authored by Yan Zhaojun's avatar Yan Zhaojun
Browse files

update

parent 07fc48c3
Pipeline #7092 failed with stage
in 0 seconds
Hide whitespace changes
Inline Side-by-side
csst_ifs_sim/csst_ifs_sim.py
View file @ 8165a134
......@@ -28,7 +28,6 @@ import astropy.coordinates as coord
import ctypes
import sys
##sys.path.append('./csst_ifs_sim')
conf.auto_max_age = None
......@@ -1149,6 +1148,25 @@ def beta_angle(x_sat, y_sat, z_sat, vx_sat, vy_sat, vz_sat, ra_obj, dec_obj):
return angle
###############################################################################
def find_min(arr):
min_val = arr[0]
min_index = 0
for i in range(1, len(arr)):
if arr[i] < min_val:
min_val = arr[i]
min_index = i
return min_val, min_index
#################################
def find_max(arr):
max_val = arr[0]
max_index = 0
for i in range(1, len(arr)):
if arr[i] > max_val:
max_val = arr[i]
max_index = i
return max_val, max_index
##########################################################
......@@ -1534,12 +1552,16 @@ def get_dx_dy_blue(wave):
"""
# wave is the wavelength in nm;
# dx is in dispersion direction, dy is in vertical direction;
dydl = np.array([-423.256, 0.001, 0.00075,
0.0000078, -0.0000000000007, 0.0, 0.0])
# 色散方向
dxdl = 0.2*np.array([-9.1519, -1.00000000e-06, 3.50000000e-08, -5.00000000e-09,
-1.70000000e-11, 4.00949787e-12, -6.16873452e-15])
# 垂直方向
# dydl = np.array([-423.256, 0.001, 0.00075,
# 0.0000078, -0.0000000000007, 0.0, 0.0])
# # 色散方向
# dxdl = 0.2*np.array([-9.1519, -1.00000000e-06, 3.50000000e-08, -5.00000000e-09,
# -1.70000000e-11, 4.00949787e-12, -6.16873452e-15])
##### update @2024.10.16
dydl=np.array([ 2447.9, -1/0.141, 0.0000075, 0.00000078, -0.000000000007] ) ; #色散方向
dxdl=np.array([ 5.46, -1.5e-02, 3.5e-08, -5.0e-09]) #垂直方向
dx = 0.0
dy = 0.0
for i in range(len(dxdl)):
......@@ -1567,11 +1589,14 @@ def get_dx_dy_red(wave):
"""
# wave is the wavelength in nm;
dydl = np.array([-640.0239901372472, 0.0018, 0.00048,
0.0000028, -0.0000000000007, 0.0, 0.0]) # 色散方向
dxdl = 0.00325*np.array([-1638.8, 4.0e-2, 5.500e-3, -
5.2e-10, 1.7000e-10, 7.1e-13, -5.16e-15]) # 垂直方向
# dydl = np.array([-640.0239901372472, 0.0018, 0.00048,
# 0.0000028, -0.0000000000007, 0.0, 0.0]) # 色散方向
# dxdl = 0.00325*np.array([-1638.8, 4.0e-2, 5.500e-3, -
# 5.2e-10, 1.7000e-10, 7.1e-13, -5.16e-15]) # 垂直方向
## update @2014.10.17
dydl=np.array([3519.78622, -1/0.1555, 0.0000048, 0.00000028, -0.0000000000007] ) #色散方向
dxdl=np.array([-5.6305, 1.0e-2, 5.500e-7, -5.2e-10]) #垂直方向
dx = 0.0
dy = 0.0
for i in range(len(dxdl)):
......@@ -1753,6 +1778,9 @@ class IFSsimulator():
self.readoutNoise = self.config.getboolean(
self.section, 'readoutnoise')
self.appbianpai = self.config.getboolean(
self.section, 'appbianpai')
self.intscale = True
######################################################################
......@@ -1769,7 +1797,8 @@ class IFSsimulator():
bleeding=self.bleeding,
sky_noise=self.sky_noise,
intscale=self.intscale,
save_cosmicrays=self.save_cosmicrays)
save_cosmicrays=self.save_cosmicrays,
appbianpai=self.appbianpai)
############################################################################
......@@ -2134,9 +2163,14 @@ class IFSsimulator():
slice_red['py'][i] = 50+250+randRedpos[i]*4
slice_red['px'][i] = 3.55/0.015*(i-16)+1190.0+118
#######
self.slice_blue = slice_blue
self.slice_red = slice_red
###### flip the fringe up to down,down to up@2024.10.16
self.slice_blue = dict()
self.slice_red = dict()
self.slice_blue['py'] =2000-slice_blue['py']
self.slice_blue['px'] = slice_blue['px']
self.slice_red['py'] =3000-slice_red['py']
self.slice_red['px'] =slice_red['px']
#######################################################################
maskSlice = dict()
......@@ -2223,7 +2257,7 @@ class IFSsimulator():
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
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)
......@@ -2439,35 +2473,35 @@ class IFSsimulator():
# blue zone 4
self.image_b[0:1024, 0:2048] += self.information['exptime'] * \
self.information['dark4_b']
self.information['dark2_b']
########## zone 1 #################
self.image_b[1024:2048, 0:2048] += self.information['exptime'] * \
self.information['dark1_b']
self.information['dark3_b']
########## zone 3 ###################
self.image_b[0:1024, 2048:4096] += self.information['exptime'] * \
self.information['dark3_b']
self.information['dark1_b']
# zone 2
self.image_b[1024:2048, 2048:4096] += self.information['exptime'] * \
self.information['dark2_b']
self.information['dark4_b']
# red zone 4
self.image_r[0:1536, 0:3072] += self.information['exptime'] * \
self.information['dark4_r']
self.information['dark2_r']
########## zone 1 #################
self.image_r[1536:3712, 0:3072] += self.information['exptime'] * \
self.information['dark1_r']
self.information['dark3_r']
########## zone 3 ###################
self.image_r[0:1536, 3072:6144] += self.information['exptime'] * \
self.information['dark3_r']
self.information['dark1_r']
# zone 2
self.image_r[1536:3072, 3072:6144] += self.information['exptime'] * \
self.information['dark2_r']
self.information['dark4_r']
##############################################################################
......@@ -2715,44 +2749,44 @@ class IFSsimulator():
overscan = 320
# zone 4 , OSH #### 1024, 2048
self.image_b[0:1024+overscan, 0:2048+prescan+overscan] += np.random.normal(
loc=0.0, scale=self.information['rn4_b'], size=(1344, 2418))
loc=0.0, scale=self.information['rn2_b'], size=(1344, 2418))
########## zone 3, OSG #################
np.random.seed()
self.image_b[0:1024+overscan, 2418:2418+2048+prescan+overscan] += np.random.normal(
loc=0.0, scale=self.information['rn3_b'], size=(1344, 2418))
loc=0.0, scale=self.information['rn1_b'], size=(1344, 2418))
########## zone 1, OSE ###################
np.random.seed()
self.image_b[0:1024+overscan, 2418*2:2418*2+2048+prescan+overscan] += np.random.normal(
loc=0.0, scale=self.information['rn1_b'], size=(1344, 2418))
loc=0.0, scale=self.information['rn3_b'], size=(1344, 2418))
########## zone 2, OSF ###############
np.random.seed()
self.image_b[0:1024+overscan, 2418*3:2418*3+2048+prescan+overscan] += np.random.normal(
loc=0.0, scale=self.information['rn2_b'], size=(1344, 2418))
loc=0.0, scale=self.information['rn4_b'], size=(1344, 2418))
self.log.info('readnoise added in blue channel')
# red zone 4 , OSH### 1536* 3072
np.random.seed()
self.image_r[0:1536+overscan, 0:3072+prescan+overscan] += np.random.normal(
loc=0.0, scale=self.information['rn4_r'], size=(1856, 3442))
loc=0.0, scale=self.information['rn2_r'], size=(1856, 3442))
########## zone 3 ,OSG #################
np.random.seed()
self.image_r[0:1536+overscan, 3442:3442+3072+prescan+overscan] += np.random.normal(
loc=0.0, scale=self.information['rn3_r'], size=(1856, 3442))
loc=0.0, scale=self.information['rn1_r'], size=(1856, 3442))
########## zone 1 ,OSE ###################
np.random.seed()
self.image_r[0:1536+overscan, 3442*2:3442*2+3072+prescan+overscan] += np.random.normal(
loc=0.0, scale=self.information['rn1_r'], size=(1856, 3442))
loc=0.0, scale=self.information['rn3_r'], size=(1856, 3442))
########## zone 2,OSF ###########
np.random.seed()
self.image_r[0:1536+overscan, 3442*3:3442*3+3072+prescan+overscan] += np.random.normal(
loc=0.0, scale=self.information['rn2_r'], size=(1856, 3442))
loc=0.0, scale=self.information['rn4_r'], size=(1856, 3442))
##########################################################################################
......@@ -2880,31 +2914,44 @@ class IFSsimulator():
self.log.info(
'Converting from electrons to ADUs using a factor of gain')
# blue zone 4
self.image_b[0:1344, 0:2418] /= self.information['gain4_b']
# # blue zone 4
# self.image_b[0:1344, 0:2418] /= self.information['gain4_b']
########## zone 3 #################
self.image_b[0:1344, 2418:2418*2] /= self.information['gain3_b']
# ########## zone 3 #################
# self.image_b[0:1344, 2418:2418*2] /= self.information['gain3_b']
########## zone 1 ###################
self.image_b[0:1344, 2418*2:2418*3] /= self.information['gain1_b']
# ########## zone 1 ###################
# self.image_b[0:1344, 2418*2:2418*3] /= self.information['gain1_b']
# zone 2
self.image_b[0:1344, 2418*3:2418*4] /= self.information['gain2_b']
# # zone 2
# self.image_b[0:1344, 2418*3:2418*4] /= self.information['gain2_b']
################# update @2024.10.18 ###
# first part, menas old blue zone 4
self.image_b[0:1344, 0:2418] /= self.information['gain1_b']
##########second part, means old zone 3 #################
self.image_b[0:1344, 2418:2418*2] /= self.information['gain2_b']
##########third part, means old zone 1 ###################
self.image_b[0:1344, 2418*2:2418*3] /= self.information['gain3_b']
#### fourth part, means old zone 2
self.image_b[0:1344, 2418*3:2418*4] /= self.information['gain4_b']
############################################################################
##########################
# red zone 4
self.image_r[0:1856, 0:3442] /= self.information['gain4_r']
self.image_r[0:1856, 0:3442] /= self.information['gain1_r']
########## zone 3 #################
self.image_r[0:1856, 3442:3442*2] /= self.information['gain3_r']
self.image_r[0:1856, 3442:3442*2] /= self.information['gain2_r']
########## zone 1 ###################
self.image_r[0:1856, 3442*2:3442*3] /= self.information['gain1_r']
self.image_r[0:1856, 3442*2:3442*3] /= self.information['gain3_r']
# zone 2
self.image_r[0:1856, 3442*3:3442*4] /= self.information['gain2_r']
self.image_r[0:1856, 3442*3:3442*4] /= self.information['gain4_r']
# 3
......@@ -2979,16 +3026,16 @@ class IFSsimulator():
########################################################################
# blue zone 4
self.image_b[0:1344, 0:2418] += self.information['bias4_b']
self.image_b[0:1344, 0:2418] += self.information['bias2_b']
########## zone 3 #################
self.image_b[0:1344, 2418:2418*2] += self.information['bias3_b']
self.image_b[0:1344, 2418:2418*2] += self.information['bias1_b']
########## zone 1 ###################
self.image_b[0:1344, 2418*2:2418*3] += self.information['bias1_b']
self.image_b[0:1344, 2418*2:2418*3] += self.information['bias3_b']
# zone 2
self.image_b[0:1344, 2418*3:2418*4] += self.information['bias2_b']
self.image_b[0:1344, 2418*3:2418*4] += self.information['bias4_b']
############################################################################
......@@ -2997,16 +3044,16 @@ class IFSsimulator():
#######################################################################
# red zone 4
self.image_r[0:1856, 0:3442] += self.information['bias4_r']
self.image_r[0:1856, 0:3442] += self.information['bias2_r']
########## zone 3 #################
self.image_r[0:1856, 3442:3442*2] += self.information['bias3_r']
self.image_r[0:1856, 3442:3442*2] += self.information['bias1_r']
########## zone 1 ###################
self.image_r[0:1856, 3442*2:3442*3] += self.information['bias1_r']
self.image_r[0:1856, 3442*2:3442*3] += self.information['bias3_r']
# zone 2
self.image_r[0:1856, 3442*3:3442*4] += self.information['bias2_r']
self.image_r[0:1856, 3442*3:3442*4] += self.information['bias4_r']
#######################################################################
......@@ -3290,14 +3337,46 @@ class IFSsimulator():
overscan = int(self.information['overscan'])
temp = np.zeros((1344, 9672))
# zone 4, OSH ######0:1024, 50:2048+50
# # zone 4, OSH ######0:1024, 50:2048+50
# x1 = 0
# x2 = x1+1024
# y1 = 0+prescan
# y2 = y1+2048
# temp[x1:x2, y1:y2] = imgb[0:1024, 0:2048]
# # zone 3, OSG , left to right #################
# # np.fliplr(b2) ## left to right
# # np.flipud(b3) ## down to up
# x1 = 0
# x2 = x1+1024
# y1 = 2418+prescan
# y2 = y1+2048
# temp[x1:x2, y1:y2] = np.fliplr(imgb[0:1024, 2048:4096])
# # zone 1, OSE,down to up ###################
# x1 = 0
# x2 = x1+1024
# y1 = 2418*2+prescan
# y2 = y1+2048
# temp[x1:x2, y1:y2] = np.flipud(imgb[1024:2048, 0:2048])
# ########## zone 2, OSF down to yp ,left to right #######
# x1 = 0
# x2 = x1+1024
# y1 = 2418*3+prescan
# y2 = y1+2048
# temp[x1:x2, y1:y2] = np.flipud(np.fliplr(imgb[1024:2048, 2048:4096]))
### update 2024.10.18
# first part, old OSG part ,## shift: left to right
x1 = 0
x2 = x1+1024
y1 = 0+prescan
y2 = y1+2048
temp[x1:x2, y1:y2] = imgb[0:1024, 0:2048]
# zone 3, OSG , left to right #################
temp[x1:x2, y1:y2] = np.fliplr(imgb[0:1024, 2048:4096])
# second part, old OSH, no shift #################
# np.fliplr(b2) ## left to right
# np.flipud(b3) ## down to up
x1 = 0
......@@ -3305,21 +3384,22 @@ class IFSsimulator():
y1 = 2418+prescan
y2 = y1+2048
temp[x1:x2, y1:y2] = np.fliplr(imgb[0:1024, 2048:4096])
# zone 1, OSE,down to up ###################
temp[x1:x2, y1:y2] = imgb[0:1024, 0:2048]
### third part, old OSE, down to up ###################
x1 = 0
x2 = x1+1024
y1 = 2418*2+prescan
y2 = y1+2048
temp[x1:x2, y1:y2] = np.flipud(imgb[1024:2048, 0:2048])
########## zone 2, OSF down to yp ,left to right #######
########## fourth part, old OSF part; down to yp ,left to right #######
x1 = 0
x2 = x1+1024
y1 = 2418*3+prescan
y2 = y1+2048
temp[x1:x2, y1:y2] = np.flipud(np.fliplr(imgb[1024:2048, 2048:4096]))
self.image_b = temp
......@@ -3327,14 +3407,46 @@ class IFSsimulator():
imgr = self.image_r.copy()
temp = np.zeros((1856, 13768))
# zone 4, OSH ######0:1024, 50:2048+50
# # zone 4, OSH ######0:1024, 50:2048+50
# x1 = 0
# x2 = x1+1536
# y1 = 0+prescan
# y2 = y1+3072
# temp[x1:x2, y1:y2] = imgr[0:1536, 0:3072]
# ########## zone 3, OSG , left to right #################
# # np.fliplr(b2) ## left to right
# # np.flipud(b3) ## down to up
# x1 = 0
# x2 = x1+1536
# y1 = 3442+prescan
# y2 = y1+3072
# temp[x1:x2, y1:y2] = np.fliplr(imgr[0:1536, 3072:6144])
# ########## zone 1, OSE,down to up ############################
# x1 = 0
# x2 = x1+1536
# y1 = 3442*2+prescan
# y2 = y1+3072
# temp[x1:x2, y1:y2] = np.flipud(imgr[1536:3072, 0:3072])
# ########## zone 2, OSF down to yp ,left to right ################
# x1 = 0
# x2 = x1+1536
# y1 = 3442*3+prescan
# y2 = y1+3072
# temp[x1:x2, y1:y2] = np.flipud(np.fliplr(imgr[1536:3072, 3072:6144]))
###### update @2024.10.18
# readout image ,first part, old OSG, shift: left to right
x1 = 0
x2 = x1+1536
y1 = 0+prescan
y2 = y1+3072
temp[x1:x2, y1:y2] = imgr[0:1536, 0:3072]
########## zone 3, OSG , left to right #################
temp[x1:x2, y1:y2] = np.fliplr(imgr[0:1536, 3072:6144])
########## second part, old OSH ,no change; #################
# np.fliplr(b2) ## left to right
# np.flipud(b3) ## down to up
x1 = 0
......@@ -3342,21 +3454,22 @@ class IFSsimulator():
y1 = 3442+prescan
y2 = y1+3072
temp[x1:x2, y1:y2] = np.fliplr(imgr[0:1536, 3072:6144])
########## zone 1, OSE,down to up ############################
temp[x1:x2, y1:y2] = imgr[0:1536, 0:3072]
########## third part , old OSE , down to up ############################
x1 = 0
x2 = x1+1536
y1 = 3442*2+prescan
y2 = y1+3072
temp[x1:x2, y1:y2] = np.flipud(imgr[1536:3072, 0:3072])
########## zone 2, OSF down to yp ,left to right ################
########## fourth part, old OSF, down to up ,left to right ################
x1 = 0
x2 = x1+1536
y1 = 3442*3+prescan
y2 = y1+3072
temp[x1:x2, y1:y2] = np.flipud(np.fliplr(imgr[1536:3072, 3072:6144]))
self.image_r = temp
......@@ -3383,6 +3496,18 @@ class IFSsimulator():
Updates header with the input values and flags used during simulation.
"""
######
###### add random pointing error to telescope, the telescope poingt parameter in
###### fits header have random pointing error
ud_ra = np.random.random() # Choose a random shift in arcsec
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
HeaderTest = 'no'
sim_ver = str(self.information['sim_ver'])
......@@ -4760,24 +4885,38 @@ class IFSsimulator():
self.information['ra_obj'] = a[0].header['RA'] # in degree
self.information['dec_obj'] = a[0].header['DEC'] # in degree
disRa = (telRa-skyRa)/3600.0 # convert unit of arcsec to degree
disRa = (telRa-skyRa)/3600.0 # convert unit of arcsec to degree
disDec = (telDec-skyDec)/3600.0 # convert unit of arcsec to degree
# RA_PNT0 =OBJ_RA + DIS_RA/cos(OBJ_DEC), DEC_PNT0= OBJ_DEC DIS_DEC
self.information['ra_pnt0'] = a[0].header['RA'] + \
disRa/np.cos(a[0].header['DEC']/180.0*np.pi)
self.information['dec_pnt0'] = a[0].header['DEC'] + disDec
############### calculate the earthshine and zodiacal noise ,new code 2023.11.1 ############
###############
self.log.info('Real telescope pointing in Ra = %f, Dec = %f' % (self.information['ra_pnt0'], self.information['dec_pnt0']))
ra = self.information['ra_pnt0']
dec = self.information['dec_pnt0']
time_jd = time2jd(self.dt)
if self.appbianpai:
sn=self.simnumber-1;
x_sat = float(self.bianpai_data['x_sat'][sn*5+self.exptime_start_index])
y_sat = float(self.bianpai_data['y_sat'][sn*5+self.exptime_start_index])
z_sat = float(self.bianpai_data['z_sat'][sn*5+self.exptime_start_index])
###
else:
x_sat = float(self.orbit_pars[self.orbit_exp_num, 1])
y_sat = float(self.orbit_pars[self.orbit_exp_num, 2])
z_sat = float(self.orbit_pars[self.orbit_exp_num, 3])
x_sat = float(self.orbit_pars[self.orbit_exp_num, 1])
y_sat = float(self.orbit_pars[self.orbit_exp_num, 2])
z_sat = float(self.orbit_pars[self.orbit_exp_num, 3])
wave0, zodi0 = self.zodiacal(
ra, dec, self.TianCe_day) # erg/s/cm^2/A/arcsec^2
......@@ -4856,14 +4995,26 @@ class IFSsimulator():
width_blue = 0
################################
############## doppler effect to photons.wavelength #############
# self.orbit_pars
x_sat = float(self.orbit_pars[self.orbit_exp_num, 1])
y_sat = float(self.orbit_pars[self.orbit_exp_num, 2])
z_sat = float(self.orbit_pars[self.orbit_exp_num, 3])
vx_sat = float(self.orbit_pars[self.orbit_exp_num, 4])
vy_sat = float(self.orbit_pars[self.orbit_exp_num, 5])
vz_sat = float(self.orbit_pars[self.orbit_exp_num, 6])
if self.appbianpai:
sn=self.simnumber-1;
x_sat = float(self.bianpai_data['x_sat'][sn*5+self.exptime_start_index])
y_sat = float(self.bianpai_data['y_sat'][sn*5+self.exptime_start_index])
z_sat = float(self.bianpai_data['z_sat'][sn*5+self.exptime_start_index])
vx_sat = float(self.bianpai_data['vx_sat'][sn*5+self.exptime_start_index])
vy_sat = float(self.bianpai_data['vy_sat'][sn*5+self.exptime_start_index])
vz_sat = float(self.bianpai_data['vz_sat'][sn*5+self.exptime_start_index])
else:
# self.orbit_pars
x_sat = float(self.orbit_pars[self.orbit_exp_num, 1])
y_sat = float(self.orbit_pars[self.orbit_exp_num, 2])
z_sat = float(self.orbit_pars[self.orbit_exp_num, 3])
vx_sat = float(self.orbit_pars[self.orbit_exp_num, 4])
vy_sat = float(self.orbit_pars[self.orbit_exp_num, 5])
vz_sat = float(self.orbit_pars[self.orbit_exp_num, 6])
self.information['POSI0_X'] = x_sat
self.information['POSI0_Y'] = y_sat
......@@ -4881,56 +5032,69 @@ class IFSsimulator():
self.information['ra_obj'], self.information['dec_obj'], v1, self.TianCe_day)
#################################################
# exposure end time is t2 ;
t2 = self.dt+timedelta(seconds=self.information['exptime'])
### data read time is the exposure end time plus readouttime ###
# data read end time is t3
t2jd = time2jd(t2)
if self.orbit_pars[-1, 0] < t2jd: # orbit parameters are not in currenct txt file
self.orbit_file_num = self.orbit_file_num+1
fn = self.information['dir_path'] + \
'IFS_inputdata/TianCe/orbit20160925/' + \
str(self.orbit_file_num)+'.txt'
self.orbit_pars = np.loadtxt(fn)
self.orbit_exp_num = 0
for k in range(self.orbit_exp_num, len(self.orbit_pars), 1):
if t2jd-self.orbit_pars[k, 0] < 0:
break
if k == 0:
deltaT = jd2time(self.orbit_pars[k, 0])-t2
p1x = self.orbit_pars[k, 1]-(self.orbit_pars[k+1, 1] -
self.orbit_pars[k, 1])*deltaT.seconds/120
p1y = self.orbit_pars[k, 2]-(self.orbit_pars[k+1, 2] -
self.orbit_pars[k, 2])*deltaT.seconds/120
p1z = self.orbit_pars[k, 3]-(self.orbit_pars[k+1, 3] -
self.orbit_pars[k, 3])*deltaT.seconds/120
p1vx = self.orbit_pars[k, 4]-(self.orbit_pars[k+1, 4] -
self.orbit_pars[k, 4])*deltaT.seconds/120
p1vx = self.orbit_pars[k, 5]-(self.orbit_pars[k+1, 5] -
self.orbit_pars[k, 5])*deltaT.seconds/120
p1vx = self.orbit_pars[k, 6]-(self.orbit_pars[k+1, 6] -
self.orbit_pars[k, 6])*deltaT.seconds/120
if self.appbianpai:
sn=self.simnumber-1;
p1x = float(self.bianpai_data['x_sat'][sn*5+self.exptime_end_index])
p1y = float(self.bianpai_data['y_sat'][sn*5+self.exptime_end_index])
p1z = float(self.bianpai_data['z_sat'][sn*5+self.exptime_end_index])
p1vx = float(self.bianpai_data['vx_sat'][sn*5+self.exptime_end_index])
p1vy = float(self.bianpai_data['vy_sat'][sn*5+self.exptime_end_index])
p1vz = float(self.bianpai_data['vz_sat'][sn*5+self.exptime_end_index])
else:
deltaT = jd2time(self.orbit_pars[k, 0])-t2
p1x = self.orbit_pars[k-1, 1]+(self.orbit_pars[k, 1] -
self.orbit_pars[k-1, 1])*deltaT.seconds/120
p1y = self.orbit_pars[k-1, 2]+(self.orbit_pars[k, 2] -
self.orbit_pars[k-1, 2])*deltaT.seconds/120
p1z = self.orbit_pars[k-1, 3]+(self.orbit_pars[k, 3] -
self.orbit_pars[k-1, 3])*deltaT.seconds/120
p1vx = self.orbit_pars[k-1, 4]+(self.orbit_pars[k, 4] -
self.orbit_pars[k-1, 4])*deltaT.seconds/120
p1vy = self.orbit_pars[k-1, 5]+(self.orbit_pars[k, 5] -
self.orbit_pars[k-1, 5])*deltaT.seconds/120
p1vz = self.orbit_pars[k-1, 6]+(self.orbit_pars[k, 6] -
self.orbit_pars[k-1, 6])*deltaT.seconds/120
# exposure end time is t2 ;
t2 = self.dt+timedelta(seconds=self.information['exptime'])
### data read time is the exposure end time plus readouttime ###
# data read end time is t3
t2jd = time2jd(t2)
if self.orbit_pars[-1, 0] < t2jd: # orbit parameters are not in currenct txt file
self.orbit_file_num = self.orbit_file_num+1
fn = self.information['dir_path'] + \
'IFS_inputdata/TianCe/orbit20160925/' + \
str(self.orbit_file_num)+'.txt'
self.orbit_pars = np.loadtxt(fn)
self.orbit_exp_num = 0
for k in range(self.orbit_exp_num, len(self.orbit_pars), 1):
if t2jd-self.orbit_pars[k, 0] < 0:
break
if k == 0:
deltaT = jd2time(self.orbit_pars[k, 0])-t2
p1x = self.orbit_pars[k, 1]-(self.orbit_pars[k+1, 1] -
self.orbit_pars[k, 1])*deltaT.seconds/120
p1y = self.orbit_pars[k, 2]-(self.orbit_pars[k+1, 2] -
self.orbit_pars[k, 2])*deltaT.seconds/120
p1z = self.orbit_pars[k, 3]-(self.orbit_pars[k+1, 3] -
self.orbit_pars[k, 3])*deltaT.seconds/120
p1vx = self.orbit_pars[k, 4]-(self.orbit_pars[k+1, 4] -
self.orbit_pars[k, 4])*deltaT.seconds/120
p1vy = self.orbit_pars[k, 5]-(self.orbit_pars[k+1, 5] -
self.orbit_pars[k, 5])*deltaT.seconds/120
p1vz = self.orbit_pars[k, 6]-(self.orbit_pars[k+1, 6] -
self.orbit_pars[k, 6])*deltaT.seconds/120
else:
deltaT = jd2time(self.orbit_pars[k, 0])-t2
p1x = self.orbit_pars[k-1, 1]+(self.orbit_pars[k, 1] -
self.orbit_pars[k-1, 1])*deltaT.seconds/120
p1y = self.orbit_pars[k-1, 2]+(self.orbit_pars[k, 2] -
self.orbit_pars[k-1, 2])*deltaT.seconds/120
p1z = self.orbit_pars[k-1, 3]+(self.orbit_pars[k, 3] -
self.orbit_pars[k-1, 3])*deltaT.seconds/120
p1vx = self.orbit_pars[k-1, 4]+(self.orbit_pars[k, 4] -
self.orbit_pars[k-1, 4])*deltaT.seconds/120
p1vy = self.orbit_pars[k-1, 5]+(self.orbit_pars[k, 5] -
self.orbit_pars[k-1, 5])*deltaT.seconds/120
p1vz = self.orbit_pars[k-1, 6]+(self.orbit_pars[k, 6] -
self.orbit_pars[k-1, 6])*deltaT.seconds/120
#######
self.information['POSI1_X'] = p1x
......@@ -4961,7 +5125,7 @@ class IFSsimulator():
ilam = klam*alfa
if ilam >= 6000:
if ilam >= 4000:
break
# print(ilam)
......@@ -5692,7 +5856,7 @@ class IFSsimulator():
self.debug = self.information['debug']
self.dt = datetime.utcnow()
if self.information['exptime']>2000 or self.information['exptime']<0:
......@@ -5718,65 +5882,89 @@ class IFSsimulator():
self.information['sky_fitsin']
print('self.skyfilepath = ', self.skyfilepath)
# self.earthshine_theta=30.0 # in degree
##################################################################
#### load orbit parameters #####
flag = 0
self.dt_num = int(
self.simnumber*(self.information['exptime']+self.information['readouttime']+125)/120)
now_dt = datetime.utcnow()
now_jd = time2jd(now_dt)
for k in range(1, 50, 1):
# fn=father_path+'/IFS_inputdata/TianCe/orbit20160925/'+str(k)+'.txt';
fn = self.information['dir_path'] + \
'IFS_inputdata/TianCe/orbit20160925/'+str(k)+'.txt'
d = np.loadtxt(fn)
for kk in range(len(d[:, 0])):
if now_jd-d[kk, 0] <= 0:
flag = 1
break
if flag == 1:
break
################################
if kk + self.dt_num < len(d[:, 0]):
exptime_start_jd = d[kk+self.dt_num, 0]
self.orbit_pars = d
self.orbit_file_num = k
self.orbit_exp_num = kk
###################################################################
if self.appbianpai:
### load yunxingbianpai csv file
############ load star data catlog #####################
starcat=self.information['bianpai_file']
##starcat='selection_20230517_concat.fits'
###################################################
self.log.info('Stat catlog file name is %s' % (starcat))
##########################################
df=pd.read_csv(self.information['dir_path']+'IFS_inputdata/TianCe/'+starcat)
###################################################################
sn=self.simnumber-1;
arr=np.array(df['time'][sn*5:sn*5+5]);
self.exptime_start_jd,self.exptime_start_index=find_min(arr);
self.exptime_end_jd, self.exptime_end_index =find_max(arr);
###self.earthshine_theta=df['earth_angle'][sn*5+index] # in degree
self.dt = jd2time(self.exptime_start_jd);
self.bianpai_data=df;
##################################################################
else:
#### load orbit parameters #####
flag = 0
self.dt = datetime.utcnow()
self.dt_num = int(
self.simnumber*(self.information['exptime']+self.information['readouttime']+125)/120)
now_dt = datetime.utcnow()
now_jd = time2jd(now_dt)
for k in range(1, 50, 1):
# fn=father_path+'/IFS_inputdata/TianCe/orbit20160925/'+str(k)+'.txt';
fn = self.information['dir_path'] + \
'IFS_inputdata/TianCe/orbit20160925/'+str(k)+'.txt'
d = np.loadtxt(fn)
for kk in range(len(d[:, 0])):
if now_jd-d[kk, 0] <= 0:
flag = 1
break
if flag == 1:
break
##################################################################
if kk + self.dt_num < len(d[:, 0]):
exptime_start_jd = d[kk+self.dt_num, 0]
self.orbit_pars = d
self.orbit_file_num = k
self.orbit_exp_num = kk
else:
fn = self.information['dir_path'] + \
'IFS_inputdata/TianCe/orbit20160925/'+str(k+1)+'.txt'
d = np.loadtxt(fn)
self.orbit_pars = d
self.orbit_file_num = k+1
self.orbit_exp_num = self.dt_num
self.dt = jd2time(exptime_start_jd)
fn = self.information['dir_path'] + \
'IFS_inputdata/TianCe/orbit20160925/'+str(k+1)+'.txt'
d = np.loadtxt(fn)
self.orbit_pars = d
self.orbit_file_num = k+1
self.orbit_exp_num = self.dt_num
##########################################
########################################################################
#self.dt=julian.from_jd(exptime_start_jd, fmt='jd')
self.dt = jd2time(exptime_start_jd)
##################################################################
##################################################################
# str(self.dt.year)+'-'+str(self.dt.month)+'-'+str(self.dt.day)
self.TianCe_day = self.dt.strftime("%Y-%m-%d")
self.TianCe_exp_start = dt2hmd(self.dt)
self.zodiacal_time = self.TianCe_day
###################################################################
skyRa = 0
skyDec = 0
......@@ -5787,6 +5975,7 @@ class IFSsimulator():
sky_rot = 2 * (ud-0.5) * 5
dsmax = np.floor(50*np.cos(sky_rot/180*np.pi))-34
np.random.seed(11*self.simnumber)
ud = np.random.random() # Choose a random shift in arcsec
telRa = 2 * (ud-0.5) * dsmax * 0.1 # in arcsec
......@@ -5813,6 +6002,7 @@ class IFSsimulator():
#######################################################################
elif self.source == 'FLAT':
self.dt = datetime.utcnow()
self.information['exptime'] = 200+simnumber*100
self.sim_calibration_img(self.information['exptime'], 'FLAT')
self.information['sky_rot'] = 0
......@@ -5820,6 +6010,7 @@ class IFSsimulator():
########
elif self.source == 'LAMP':
self.dt = datetime.utcnow()
self.information['exptime'] = 200+simnumber*100
self.information['sky_rot'] = 0
......@@ -5831,12 +6022,14 @@ class IFSsimulator():
#########
elif self.source == 'DARK':
self.dt = datetime.utcnow()
self.information['sky_rot'] = 0
self.information['exptime'] = simnumber*3600*24
self.image_b = np.zeros((2048, 4096))
self.image_r = np.zeros((3072, 6144))
elif self.source == 'BIAS':
self.dt = datetime.utcnow()
self.information['sky_rot'] = 0
self.information['exptime'] = 0
self.image_b = np.zeros((2048, 4096))
......
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