Skip to content
GitLab
Commit 9e7cecd1 authored by Yan Zhaojun's avatar Yan Zhaojun
Browse files

debug

parent 2b7c8ba3
Pipeline #4479 passed with stage
in 0 seconds
Show whitespace changes
Inline Side-by-side
csst_mci_sim/csst_mci_sim.py
View file @ 9e7cecd1
...@@ -55,79 +55,36 @@ Contact Information: zhaojunyan@shao.ac.cn ...@@ -55,79 +55,36 @@ Contact Information: zhaojunyan@shao.ac.cn
""" """
import galsim import galsim
import pandas as pd import pandas as pd
from scipy.integrate import simps #from scipy.integrate import simps
from datetime import datetime, timedelta from datetime import datetime, timedelta
from astropy.time import Time from astropy.time import Time
from astropy.coordinates import get_sun from astropy.coordinates import get_sun
# import scipy.io as sio
# import cmath
import numpy as np import numpy as np
from tqdm import tqdm from tqdm import tqdm
from astropy.table import Table from astropy.table import Table
# from astropy.wcs import WCS as WCS
from astropy.io import fits from astropy.io import fits
from astropy import units as u from astropy import units as u
import os, sys, math import os, sys, math
import configparser as ConfigParser import configparser as ConfigParser
from matplotlib import pyplot as plt
# from optparse import OptionParser
# from matplotlib import pyplot as plt
from scipy import ndimage from scipy import ndimage
sys.path.append('./csst_mci_sim') sys.path.append('./csst_mci_sim')
##sys.path.append('..')
from CTI import CTI from CTI import CTI
from support import logger as lg from support import logger as lg
from support import cosmicrays from support import cosmicrays
from support import shao from support import shao
from support import sed from support import sed
#from shao import onOrbitObsPosition
from support import MCIinstrumentModel from support import MCIinstrumentModel
from joblib import Parallel, delayed from joblib import Parallel, delayed
#import multiprocessing
from astropy.coordinates import SkyCoord from astropy.coordinates import SkyCoord
from scipy import interpolate from scipy import interpolate
from scipy.signal import fftconvolve from scipy.signal import fftconvolve
import pandas import pandas
# sys.path.append('../MCI_inputData/TianCe')
# sys.path.append('../MCI_inputData/SED_Code')
# set the folder
FOLDER ='../'
#####################################################################################
###############################################################################
import ctypes import ctypes
import astropy.coordinates as coord import astropy.coordinates as coord
from scipy.interpolate import interp1d from scipy.interpolate import interp1d
########################### functions #########################
def transRaDec2D(ra, dec): def transRaDec2D(ra, dec):
# radec转为竞天程序里的ob, 赤道坐标系下的笛卡尔三维坐标xyz. # radec转为竞天程序里的ob, 赤道坐标系下的笛卡尔三维坐标xyz.
x1 = np.cos(dec / 57.2957795) * np.cos(ra / 57.2957795) x1 = np.cos(dec / 57.2957795) * np.cos(ra / 57.2957795)
...@@ -191,41 +148,6 @@ def earth_angle(time_jd, x_sat, y_sat, z_sat, ra_obj, dec_obj): ...@@ -191,41 +148,6 @@ def earth_angle(time_jd, x_sat, y_sat, z_sat, ra_obj, dec_obj):
angle = 180 - angle_a - angle_b.degree angle = 180 - angle_a - angle_b.degree
return angle return angle
#################################################3
# def MCIinformation():
# """
# Returns a dictionary describing MCI CCD. The following information is provided (id: value - reference)::
# dob: 0 - CDM03 (Short et al. 2010)
# fwc: 90000 - CCD spec EUCL-EST-RS-6-002 (for CDM03)
# rdose: 30000000000.0 - derived (above the PLM requirement)
# sfwc: 730000.0 - CDM03 (Short et al. 2010), see also the CCD spec EUCL-EST-RS-6-002
# st: 5e-06 - CDM03 (Short et al. 2010)
# svg: 1e-10 - CDM03 (Short et al. 2010)
# t: 0.01024 - CDM03 (Short et al. 2010)
# trapfile: cdm_euclid.dat - CDM03 (derived, refitted to CCD204 data)
# vg: 6e-11 - CDM03 (Short et al. 2010)
# vth: 11680000.0 - CDM03 (Short et al. 2010)
# :return: instrument model parameters
# :rtype: dict
# """
# #########################################################################################################
# out=dict()
# out.update({'dob' : 0, 'rdose' : 8.0e9,
# 'parallelTrapfile' : 'cdm_euclid_parallel.dat', 'serialTrapfile' : 'cdm_euclid_serial.dat',
# 'beta_s' : 0.6, 'beta_p': 0.6, 'fwc' : 90000, 'vth' : 1.168e7, 't' : 20.48e-3, 'vg' : 6.e-11,
# 'st' : 5.0e-6, 'sfwc' : 730000., 'svg' : 1.0e-10})
# return out
####################################### #######################################
def CCDnonLinearityModel(data, beta=6e-7): def CCDnonLinearityModel(data, beta=6e-7):
...@@ -346,13 +268,7 @@ def make_c_coor(fov, step): ...@@ -346,13 +268,7 @@ def make_c_coor(fov, step):
############################################################################## ##############################################################################
############################################################################### ##############################################################################
# def str2time(strTime):
# if len(strTime)>20:#
# msec=int(float('0.'+strTime[20:])*1000000) # microsecond
# str2=strTime[0:19]+' '+str(msec)
# return datetime.strptime(str2,'%Y %m %d %H %M %S %f')
#datetime to mjd #datetime to mjd
def time2mjd(dateT): def time2mjd(dateT):
...@@ -444,18 +360,10 @@ def deg2HMS(ra0, dec0): ...@@ -444,18 +360,10 @@ def deg2HMS(ra0, dec0):
hhmmss=ss[0:2]+ss[3:5]+ra_ss hhmmss=ss[0:2]+ss[3:5]+ra_ss
return hhmmss+dms_d+dms_m+dms_s return hhmmss+dms_d+dms_m+dms_s
################################################################################ ###############################################################################
# def cut_radius(rcutstar, mstar, mag):
# return rcutstar * 10**(0.4*2*(mstar-mag)/3.125)
# def core_radius(rcorestar, mstar, mag):
# return rcorestar * 10**(0.4*(mstar-mag)/2)
# def v_disp(sigstar, mstar, mag):
# return sigstar * 10**(0.4*(mstar-mag)/3.57)
############################################################################### ###############################################################################
##################################################################################### #############################################################################
#### distortion function #### distortion function
def distortField(ra, dec, ch): def distortField(ra, dec, ch):
...@@ -586,18 +494,6 @@ def cal_pos(center_ra, center_dec, rotTel, rotSky, sRa, sDec): ...@@ -586,18 +494,6 @@ def cal_pos(center_ra, center_dec, rotTel, rotSky, sRa, sDec):
# return a.reshape(sh).sum(-1).sum(1) # return a.reshape(sh).sum(-1).sum(1)
##################################################################### #####################################################################
# def float2char(a, z):
# ### a is input float value
# ### transfer float a to chars and save z bis
# b=str(a)
# n=len(b)
# if z>=n:
# return b
# else:
# return b[:z]
######################################################### #########################################################
def centroid(data): def centroid(data):
...@@ -621,327 +517,8 @@ def centroid(data): ...@@ -621,327 +517,8 @@ def centroid(data):
############################################################################### ###############################################################################
# def fftrange(n, dtype=float):
# """FFT-aligned coordinate grid for n samples."""
# return np.arange(-n//2, -n//2+n, dtype=dtype)
###############################################################################
# def dft2(ary, Q, samples, shift=None):
# """Compute the two dimensional Discrete Fourier Transform of a matrix.
# Parameters
# ----------
# ary : `numpy.ndarray`
# an array, 2D, real or complex. Not fftshifted.
# Q : `float`
# oversampling / padding factor to mimic an FFT. If Q=2, Nyquist sampled
# samples : `int` or `Iterable`
# number of samples in the output plane.
# If an int, used for both dimensions. If an iterable, used for each dim
# shift : `float`, optional
# shift of the output domain, as a frequency. Same broadcast
# rules apply as with samples.
# Returns
# -------
# `numpy.ndarray`
# 2D array containing the shifted transform.
# Equivalent to ifftshift(fft2(fftshift(ary))) modulo output
# sampling/grid differences
# """
# # this is for dtype stabilization
# Q = float(Q) ## Q=lambda*f/D/pixelsize
# n, m = ary.shape ###ary maybe is the pupil function
# N, M = samples,samples
# X, Y, U, V = (fftrange(n) for n in (m, n, M, N))
# a = 1 / Q
# ###################################################################
# Eout_fwd = np.exp(-1j * 2 * np.pi * a / n * np.outer(Y, V).T)
# Ein_fwd = np.exp(-1j * 2 * np.pi * a / m * np.outer(X, U))
# #############################################################################
# out = Eout_fwd @ ary @ Ein_fwd #### ary is the input pupil function
# return out
##############################################################################
# def idft2(ary, Q, samples, shift=None):
# """Compute the two dimensional inverse Discrete Fourier Transform of a matrix.
# Parameters
# ----------
# ary : `numpy.ndarray`
# an array, 2D, real or complex. Not fftshifted.
# Q : `float`
# oversampling / padding factor to mimic an FFT. If Q=2, Nyquist sampled
# samples : `int` or `Iterable`
# number of samples in the output plane.
# If an int, used for both dimensions. If an iterable, used for each dim
# shift : `float`, optional
# shift of the output domain, as a frequency. Same broadcast
# rules apply as with samples.
# Returns
# -------
# `numpy.ndarray`
# 2D array containing the shifted transform.
# Equivalent to ifftshift(ifft2(fftshift(ary))) modulo output
# sampling/grid differences
# """
# # this is for dtype stabilization
# Q = float(Q) ## Q=lambda*f/D/pixelsize
# n, m = ary.shape ###ary maybe is the pupil function
# N, M = samples
# X, Y, U, V = (fftrange(n) for n in (m, n, M, N))
# ###############################################################################
# nm = n*m
# NM = N*M
# r = NM/nm
# a = 1 / Q
# ###################################################################
# # Eout_fwd = np.exp(-1j * 2 * np.pi * a / n * np.outer(Y, V).T)
# # Ein_fwd = np.exp(-1j * 2 * np.pi * a / m * np.outer(X, U))
# Eout_rev = np.exp(1j * 2 * np.pi * a / n * np.outer(Y, V).T) * (1/r)
# Ein_rev = np.exp(1j * 2 * np.pi * a / m * np.outer(X, U)) * (1/nm)
# ###########################################################################
# out = Eout_rev @ ary @ Ein_rev
# return out
###############################################################################
############################################################################### ###############################################################################
# def opd2psf(opd,cwave,oversampling):
# '''
# calculate psf from opd data with defined wavelength
# Parameters
# ----------
# opd : TYPE
# the input opd data.
# cwave : TYPE
# the defined wavelength.
# '''
# D=2;
# focLength=41.253 # % MCI focal length in meters
# pixel=10/oversampling # % pixel size in microns
# opd=opd*1e9; ## convert opd unit of meter to nm
# zoomopd=ndimage.zoom(opd, 2, order=1)
# m,n=np.shape(zoomopd);
# clambda=cwave;
# wfe=zoomopd/clambda; #%% the main telescope aberration ,in wavelength;
# Q=clambda*1e-3*focLength/D/pixel
# pupil=np.zeros((m,m));
# pupil[abs(zoomopd)>0]=1;
# #idx=pupil>0;
# Nt=512
# phase=2*np.pi*wfe; #% wavefront phase in radians
# #%generalized pupil function of channel1;
# pk=pupil*np.exp(cmath.sqrt(-1)*(phase))
# pf=dft2(pk, Q, Nt)
# pf2 = abs(pf)**2
# psfout=pf2/pf2.sum()
# cx,cy=centroid(psfout)
# psf=ndimage.shift(psfout,[Nt/2-1-cy, Nt/2-1-cx],order=1, mode='nearest' )
# return psf
########################################################################
# def cal_PSF_new(channel,wfetimes,oversampling):
# # filed postion: Ra=ys1, Dec=ys2
# ############ use opd cal PSF on the defined field;
# #% psf interpolation test 20210207
# # xfmin=-0.07 # % the filed minmum value in the x direction in degrees
# # xfmax= 0.07 #% the filed maxmum value in the x direction in degrees
# # yfmin=-0.35 #% the filed minmum value in the y direction in degrees
# # yfmax=-0.49 # % the filed maxmum value in the y direction in degrees
# cfx= 0.2*3600; # % field center in x derection in arcsec;
# cfy= -0.45*3600; # % field center in y derection in arcsec;
# #wavelength=[255,337,419,501,583,665,747,829,911,1000]
# cwave=dict()
# if channel == 'g':
# cwave[channel]=475
# waven=4
# if channel == 'r':
# cwave[channel]=625
# waven=6
# if channel == 'i':
# cwave[channel]=776
# waven=7
# n=np.arange(1,101,1)
# ffx= 0.1+ 0.02222222*((n-1)%10)
# ffy=-0.35-0.02222222*np.floor((n-0.1)/10)
# psf=dict()
# fx=dict()
# fy=dict()
# fx[channel]=ffx*3600-cfx
# fy[channel]=ffy*3600-cfy
# Nt=512
# psf[channel]=np.zeros((len(n),Nt,Nt))
# for i in range(len(n)):
# #waven
# # waven=4
# # fieldn=10
# file=self.information['dir_path']+'MCI_inputData/MCI_wavefront/wave_'+str(waven+1)+'/wavefront/opd_'+str(i+1)+'.mat'
# data=sio.loadmat(file)
# opd=data['opd'] ## opd data;
# psf[channel][i,:,:]=opd2psf(wfetimes*opd, cwave[channel],oversampling)
# hdu1=fits.PrimaryHDU(psf[channel])
# hdu1.header.append(('channel', channel, 'MCI PSF channel'))
# hdu1.header['pixsize']=(0.05/oversampling, 'PSF pixel size in arcsec')
# hdu1.header['wfetimes']=(wfetimes, 'wavefront error magnification to CSST primary wavefront ')
# dtime=datetime.datetime.utcnow().strftime('%Y -%m -%d %H: %M: %S')
# hdu1.header.add_history('PSF is generated on :'+dtime)
# hdu2=fits.ImageHDU(fx[channel])
# hdu2.header.append(('field_X', 'arcsec'))
# hdu3=fits.ImageHDU(fy[channel])
# hdu3.header.append(('field_Y', 'arcsec'))
# newd=fits.HDUList([hdu1,hdu2,hdu3])
# PSFfilename=self.information['dir_path']+'MCI_inputData/PSF/PSF_'+channel+'.fits'
# newd.writeto(PSFfilename,overwrite=True)
# return
#############################################################################
# def cal_Filter_PSF(wfetimes):
# # filed postion: Ra=ys1, Dec=ys2
# ############ use opd cal PSF on the defined field;
# #% psf interpolation test 20210207
# # xfmin=-0.07 # % the filed minmum value in the x direction in degrees
# # xfmax= 0.07 #% the filed maxmum value in the x direction in degrees
# # yfmin=-0.35 #% the filed minmum value in the y direction in degrees
# # yfmax=-0.49 # % the filed maxmum value in the y direction in degrees
# oversampling=2
# cfx= 0.2*3600; # % field center in x derection in arcsec;
# cfy= -0.45*3600; # % field center in y derection in arcsec;
# wavelist =np.array([255, 337,419,501,583,665,747,829,911,1000])
# filterP=np.load(self.information['dir_path']+'MCI_inputData/MCI_filters/mci_filterPWTC.npy',allow_pickle=True).item()
# fn=np.arange(1,101,1) ### PSF field point
# ffx= 0.1+ 0.02222222*((fn-1)%10)
# ffy=-0.35-0.02222222*np.floor((fn-0.1)/10)
# psf=dict()
# fx=dict()
# fy=dict()
# fx=ffx*3600-cfx
# fy=ffy*3600-cfy
# Nt=512
# for filterk in filterP.keys():
# filtername=filterk
# ## print(filtername)
# filterwaveC=filterP[filterk]['Clambda']
# filterFwhm=filterP[filterk]['FWHM']
# fdis=abs(filterwaveC-wavelist)
# findex=np.argsort(fdis)
# waven=findex[0]
# psf[filtername]=dict()
# psf[filtername]['psf_field_X'] =fx
# psf[filtername]['psf_field_Y'] =fy
# psf[filtername]['psf_scale'] =0.05/oversampling
# psf[filtername]['wfetimes'] =wfetimes
# psf[filtername]['psf_mat'] =np.zeros( (Nt,Nt,7,len(fn)) )
# psf[filtername]['filter_name'] =filtername
# psf[filtername]['filter_channel']=filterP[filterk]['channel']
# psf[filtername]['psf_iwave'] =np.zeros(7)
# for ii in range(len(fn)):
# #waven
# # waven=4
# # fieldn=10
# file=self.information['dir_path']+'MCI_input/MCI_wavefront/wave_'+str(waven+1)+'/wavefront/opd_'+str(ii+1)+'.mat'
# data=sio.loadmat(file)
# opd=data['opd'] ## opd data;
# for kk in range(7):
# ###
# wavek=filterwaveC+filterFwhm*(kk-3)*0.25
# psfd=opd2psf(wfetimes*opd, wavek, oversampling)
# psf[filtername]['psf_mat'][:,:,kk,ii]=psfd[:,:]
# if ii==0:
# psf[filtername]['psf_iwave'][kk]=wavek
# np.save(self.information['dir_path']+'mci_sim_result/'+filtername+'_PSF.npy', psf[filtername])
# return
#########################################################################
''' '''
PSF interpolation for MCI_sim PSF interpolation for MCI_sim
...@@ -1145,8 +722,6 @@ class MCIsimulator(): ...@@ -1145,8 +722,6 @@ class MCIsimulator():
#ghost ratio can be in engineering format, so getfloat does not capture it... #ghost ratio can be in engineering format, so getfloat does not capture it...
self.information['ghostratio'] = float(self.config.get(self.section, 'ghostRatio')) self.information['ghostratio'] = float(self.config.get(self.section, 'ghostRatio'))
################################################################################################### ###################################################################################################
################################################################################################## ##################################################################################################
...@@ -1201,27 +776,12 @@ class MCIsimulator(): ...@@ -1201,27 +776,12 @@ class MCIsimulator():
save_cosmicrays =self.save_cosmicrays ) save_cosmicrays =self.save_cosmicrays )
##################################################################### #####################################################################
# self.log.info('Using the following input values:')
# for key, value in self.information.items():
# self.log.info('%s = %s' % (key, value))
# self.log.info('Using the following booleans:')
# for key, value in self.booleans.items():
# self.log.info('%s = %s' % (key, value))
return return
######################################################################################### ##############################################################################
# def make_c_coor(self, bs, nc): ###############################################################################
# '''
# Draw the mesh grids for a bs*bs box with nc*nc pixels
# '''
# ds=bs/nc
# xx01 = np.linspace(-bs/2.0,bs/2.0-ds,nc)+0.5*ds
# xx02 = np.linspace(-bs/2.0,bs/2.0-ds,nc)+0.5*ds
# xg2,xg1 = np.meshgrid(xx01,xx02)
# return xg1,xg2
##########################################################################################
def _createEmpty(self): def _createEmpty(self):
""" """
...@@ -1236,32 +796,8 @@ class MCIsimulator(): ...@@ -1236,32 +796,8 @@ class MCIsimulator():
self.image_i=np.zeros((self.information['ysize'], self.information['xsize']), dtype=float) self.image_i=np.zeros((self.information['ysize'], self.information['xsize']), dtype=float)
return return
######################################################################################################################### ###############################################################################
##############################################################################
#########################################################################################################
# def smoothingWithChargeDiffusion(self, image, sigma=(0.32, 0.32)):
# """
# Smooths a given image with a gaussian kernel with widths given as sigmas.
# This smoothing can be used to mimic charge diffusion within the CCD.
# The default values are from Table 8-2 of CCD_273_Euclid_secification_1.0.130812.pdf converted
# to sigmas (FWHM / (2sqrt(2ln2)) and rounded up to the second decimal.
# .. Note:: This method should not be called for the full image if the charge spreading
# has already been taken into account in the system PSF to avoid float counting.
# :param image: image array which is smoothed with the kernel
# :type image: ndarray
# :param sigma: widths of the gaussian kernel that approximates the charge diffusion [0.32, 0.32].
# :param sigma: tuple
# :return: smoothed image array
# :rtype: ndarray
# """
# return ndimage.filters.gaussian_filter(image, sigma)
def _loadGhostModel(self): def _loadGhostModel(self):
""" """
Reads in a ghost model from a FITS file and stores the data to self.ghostModel. Reads in a ghost model from a FITS file and stores the data to self.ghostModel.
...@@ -1316,16 +852,12 @@ class MCIsimulator(): ...@@ -1316,16 +852,12 @@ class MCIsimulator():
##print('print information:', self.information) ##print('print information:', self.information)
########################################################################### ###########################################################################
now=datetime.utcnow() now=datetime.utcnow()
#data_time=now.strftime("%Y-%m-%d-%H-%M-%S") #data_time=now.strftime("%Y-%m-%d-%H-%M-%S")
result_day=now.strftime("%Y-%m-%d") result_day=now.strftime("%Y-%m-%d")
if self.information['dir_path'] =='/nfsdata/share/simulation-unittest/mci_sim/': if self.information['dir_path'] =='/nfsdata/share/simulation-unittest/mci_sim/':
self.result_path = self.information['dir_path'] +'mci_sim_result/'+result_day self.result_path = self.information['dir_path'] +'mci_sim_result/'+result_day
else: else:
...@@ -1339,7 +871,6 @@ class MCIsimulator(): ...@@ -1339,7 +871,6 @@ class MCIsimulator():
self.result_path = '/data/mcisimdata/'+result_day self.result_path = '/data/mcisimdata/'+result_day
if os.path.isdir(self.result_path)==False: if os.path.isdir(self.result_path)==False:
os.mkdir(self.result_path) os.mkdir(self.result_path)
os.mkdir(self.result_path+'/cali_Data') os.mkdir(self.result_path+'/cali_Data')
...@@ -1359,7 +890,6 @@ class MCIsimulator(): ...@@ -1359,7 +890,6 @@ class MCIsimulator():
self.log = lg.setUpLogger(self.result_path+'/log_Data/MCIsim_'+self.source+'_'+data_time+'_Num_'+str(simnumber)+'.log') self.log = lg.setUpLogger(self.result_path+'/log_Data/MCIsim_'+self.source+'_'+data_time+'_Num_'+str(simnumber)+'.log')
self.log.info('-------STARTING A NEW SIMULATION------------') self.log.info('-------STARTING A NEW SIMULATION------------')
for key, value in self.information.items(): for key, value in self.information.items():
...@@ -1376,16 +906,11 @@ class MCIsimulator(): ...@@ -1376,16 +906,11 @@ class MCIsimulator():
#### load telescope efficiency data #### load telescope efficiency data
self.tel_eff=np.load(self.information['dir_path']+'MCI_inputData/tel_eff/tel_eff.npy',allow_pickle=True).item() self.tel_eff=np.load(self.information['dir_path']+'MCI_inputData/tel_eff/tel_eff.npy',allow_pickle=True).item()
#### load MCI filter data #### load MCI filter data
self.filterP=np.load(self.information['dir_path']+'MCI_inputData/MCI_filters/mci_filterPWTC.npy',allow_pickle=True).item() self.filterP=np.load(self.information['dir_path']+'MCI_inputData/MCI_filters/mci_filterPWTC.npy',allow_pickle=True).item()
return return
################################################################################ ################################################################################
...@@ -1521,14 +1046,6 @@ class MCIsimulator(): ...@@ -1521,14 +1046,6 @@ class MCIsimulator():
wave = np.linspace(2500,11000, 8501) # set the wavelenth for MCI from 250nm to 1100nm wave = np.linspace(2500,11000, 8501) # set the wavelenth for MCI from 250nm to 1100nm
# wavearr =self.earthshine_wave # A
# fluxarr =self.earthshine_flux # erg/s/cm^2/A/arcsec^2
# earthshinefluxarr=np.interp(wave, wavearr, fluxarr) # erg/s/cm^2/A/arcsec^2
# wavearr =self.zodiacal_wave # A
# fluxarr =self.zodiacal_flux # erg/s/cm^2/A/arcsec^2
# zodiacalfluxarr=np.interp(wave, wavearr, fluxarr) # erg/s/cm^2/A/arcsec^2
skynoise_flux=self.earthshine_flux+self.zodiacal_flux # erg/s/cm^2/A/arcsec^2 skynoise_flux=self.earthshine_flux+self.zodiacal_flux # erg/s/cm^2/A/arcsec^2
...@@ -1607,13 +1124,6 @@ class MCIsimulator(): ...@@ -1607,13 +1124,6 @@ class MCIsimulator():
self.log.info('Stat catlog file name is %s' % (starcat)) self.log.info('Stat catlog file name is %s' % (starcat))
######################################################################
# da=fits.open(self.information['dir_path']+'MCI_inputData/star_input/'+starcat)
# df2=da[1].data
# del da
########################################## ##########################################
df2=pandas.read_csv(self.information['dir_path']+'MCI_inputData/star_input/'+starcat) df2=pandas.read_csv(self.information['dir_path']+'MCI_inputData/star_input/'+starcat)
...@@ -1753,8 +1263,6 @@ class MCIsimulator(): ...@@ -1753,8 +1263,6 @@ class MCIsimulator():
########## finish save PSF fits ########################################## ########## finish save PSF fits ##########################################
############################################################################ ############################################################################
###### add binary stars to the star catlog #######
nsrcs=len(self.star['ra_gaia']) nsrcs=len(self.star['ra_gaia'])
...@@ -1790,8 +1298,6 @@ class MCIsimulator(): ...@@ -1790,8 +1298,6 @@ class MCIsimulator():
####################################################################### #######################################################################
nlayccd = 0 nlayccd = 0
############################################ ############################################
ra = self.information['ra_pnt0'] ra = self.information['ra_pnt0']
...@@ -1803,9 +1309,6 @@ class MCIsimulator(): ...@@ -1803,9 +1309,6 @@ class MCIsimulator():
y_sat=float(self.orbit_pars[self.orbit_exp_num,2]) y_sat=float(self.orbit_pars[self.orbit_exp_num,2])
z_sat=float(self.orbit_pars[self.orbit_exp_num,3]) 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 wave0, zodi0 = self.zodiacal(ra, dec, self.TianCe_day) # erg/s/cm^2/A/arcsec^2
self.log.info('dir_path:') self.log.info('dir_path:')
...@@ -1899,9 +1402,6 @@ class MCIsimulator(): ...@@ -1899,9 +1402,6 @@ class MCIsimulator():
#################### generate star image ########## #################### generate star image ##########
for j in tqdm(range(nsrcs)): ### nsrcs length for j in tqdm(range(nsrcs)): ### nsrcs length
if j==0: if j==0:
...@@ -1913,9 +1413,7 @@ class MCIsimulator(): ...@@ -1913,9 +1413,7 @@ class MCIsimulator():
starRa = 3600.0*( newRa[j] ) # ra of star, arcsecond starRa = 3600.0*( newRa[j] ) # ra of star, arcsecond
starDec = 3600.0*( newDec[j] ) # dec of star, arcsecond starDec = 3600.0*( newDec[j] ) # dec of star, arcsecond
################################################################### ###################################################################
################################################################### ###################################################################
fsx,fsy=cal_pos(center_ra,center_dec, rotTelPos, rotSkyPos, starRa, starDec) fsx,fsy=cal_pos(center_ra,center_dec, rotTelPos, rotSkyPos, starRa, starDec)
row= fsy/pixelscale+self.information['ysize']/2-0.5 ### row number on CCD image row= fsy/pixelscale+self.information['ysize']/2-0.5 ### row number on CCD image
...@@ -1928,7 +1426,6 @@ class MCIsimulator(): ...@@ -1928,7 +1426,6 @@ class MCIsimulator():
nlayccd=nlayccd+1 nlayccd=nlayccd+1
if self.debug: if self.debug:
if nlayccd>1: if nlayccd>1:
break break
...@@ -1976,8 +1473,6 @@ class MCIsimulator(): ...@@ -1976,8 +1473,6 @@ class MCIsimulator():
####### use SED_code to generate star SED ####### ####### use SED_code to generate star SED #######
# SED of j-th star # SED of j-th star
umag = self.star['umag'][j] umag = self.star['umag'][j]
gmag = self.star['gmag'][j] gmag = self.star['gmag'][j]
rmag = self.star['rmag'][j] rmag = self.star['rmag'][j]
...@@ -1987,7 +1482,6 @@ class MCIsimulator(): ...@@ -1987,7 +1482,6 @@ class MCIsimulator():
# Input redshift # Input redshift
redshift = 0 redshift = 0
# Loading galaxy SED template # Loading galaxy SED template
t=sed.Gal_Temp(self.information['dir_path']) t=sed.Gal_Temp(self.information['dir_path'])
# Calculating the magnitude (u, g, r, i, z) of each template # Calculating the magnitude (u, g, r, i, z) of each template
...@@ -1996,8 +1490,6 @@ class MCIsimulator(): ...@@ -1996,8 +1490,6 @@ class MCIsimulator():
ugriz = np.array([umag, gmag, rmag, imag, zmag]) ugriz = np.array([umag, gmag, rmag, imag, zmag])
star_flux = sed.Model_Galaxy_SED(wave, ugriz, redshift, t, self.information['dir_path']) star_flux = sed.Model_Galaxy_SED(wave, ugriz, redshift, t, self.information['dir_path'])
##cal_SED_photons(self, flux_arr): ##cal_SED_photons(self, flux_arr):
intscales,PSF_eff_weight=self.cal_SED_photons(star_flux) intscales,PSF_eff_weight=self.cal_SED_photons(star_flux)
...@@ -2020,8 +1512,6 @@ class MCIsimulator(): ...@@ -2020,8 +1512,6 @@ class MCIsimulator():
if nlayccd %10 == 0: if nlayccd %10 == 0:
self.log.info('Star number on CCD is = %i' % (nlayccd)) self.log.info('Star number on CCD is = %i' % (nlayccd))
###### ######
st_phot_C1.append(intscales['g']) st_phot_C1.append(intscales['g'])
st_phot_C2.append(intscales['r']) st_phot_C2.append(intscales['r'])
...@@ -2115,8 +1605,6 @@ class MCIsimulator(): ...@@ -2115,8 +1605,6 @@ class MCIsimulator():
cx0,cy0=centroid(image.array) cx0,cy0=centroid(image.array)
#self.log.info('finish appfat ........')
############################################ ############################################
############################################ ############################################
photons.x=photons.x-cx0*0.025 photons.x=photons.x-cx0*0.025
...@@ -2176,8 +1664,6 @@ class MCIsimulator(): ...@@ -2176,8 +1664,6 @@ class MCIsimulator():
tab['st_magi']=np.asarray(st_magi) tab['st_magi']=np.asarray(st_magi)
tab['st_magz']=np.asarray(st_magz) tab['st_magz']=np.asarray(st_magz)
###################################33 ###################################33
...@@ -2196,63 +1682,6 @@ class MCIsimulator(): ...@@ -2196,63 +1682,6 @@ class MCIsimulator():
################################################################################ ################################################################################
################################################################################ ################################################################################
#################################################################################
########################################################################
# def earthshine(self, theta):
# """
# For given theta angle, return the earth-shine spectrum.
# :param theta: angle (in degree) from the target to earth limb.
# :return: the scaled solar spectrum
# template_wave: unit in A
# template_flux: unit in erg/s/cm^2/A/arcsec^2
# """
# # read solar template
# solar_template = pd.read_csv(self.information['dir_path']+'MCI_inputData/refs/solar_spec.dat', sep='\s+',
# header=None, comment='#')
# template_wave = solar_template[0].values
# template_flux = solar_template[1].values
# # read earth shine surface brightness
# earthshine_curve = pd.read_csv(self.information['dir_path']+'MCI_inputData/refs/earthshine.dat',
# header=None, comment='#')
# angle = earthshine_curve[0].values
# surface_brightness = earthshine_curve[1].values
# # read V-band throughtput
# cat_filter_V = pd.read_csv(self.information['dir_path']+'MCI_inputData/refs/filter_Bessell_V.dat', sep='\s+',
# header=None, comment='#')
# filter_wave = cat_filter_V[0].values
# filter_response = cat_filter_V[1].values
# # interplate to the target wavelength in V-band
# ind_filter = (template_wave >= np.min(filter_wave)) & (template_wave <= np.max(filter_wave))
# filter_wave_interp = template_wave[ind_filter]
# filter_response_interp = np.interp(filter_wave_interp, filter_wave, filter_response)
# filter_constant = simps(filter_response_interp * filter_wave_interp, filter_wave_interp)
# template_constant = simps(filter_response_interp * template_wave[ind_filter] * template_flux[ind_filter],
# template_wave[ind_filter])
# dwave = filter_wave_interp[1:] - filter_wave_interp[:-1]
# wave_eff = np.nansum(dwave * filter_wave_interp[1:] * filter_response_interp[1:]) / \
# np.nansum(dwave * filter_response_interp[1:])
# # get the normalized value at theta.
# u0 = np.interp(theta, angle, surface_brightness) # mag/arcsec^2
# u0 = 10**((u0 + 48.6)/(-2.5)) # target flux in erg/s/cm^2/Hz unit
# u0 = u0 * 3e18 / wave_eff**2 # erg/s/cm^2/A/arcsec^2
# factor = u0 * filter_constant / template_constant
# norm_flux = template_flux * factor # erg/s/cm^2/A/arcsec^2
# self.earthshine_wave=template_wave # A
# self.earthshine_flux=norm_flux
# return
# ########################################################################################################################################################################################################################################################
def zodiacal(self, ra, dec, time): def zodiacal(self, ra, dec, time):
""" """
...@@ -2821,8 +2250,6 @@ class MCIsimulator(): ...@@ -2821,8 +2250,6 @@ class MCIsimulator():
############################################################################### ###############################################################################
###############################################################################
############################################################################### ###############################################################################
def generatePRNU(self, ave=1.0, sigma=0.01): def generatePRNU(self, ave=1.0, sigma=0.01):
...@@ -2924,9 +2351,6 @@ class MCIsimulator(): ...@@ -2924,9 +2351,6 @@ class MCIsimulator():
cosmics_r = cosmicrays.cosmicrays(self.log, crImage, self.information['exptime'], crInfo=self.cr) cosmics_r = cosmicrays.cosmicrays(self.log, crImage, self.information['exptime'], crInfo=self.cr)
cosmics_i = cosmicrays.cosmicrays(self.log, crImage, self.information['exptime'], crInfo=self.cr) cosmics_i = cosmicrays.cosmicrays(self.log, crImage, self.information['exptime'], crInfo=self.cr)
#add cosmic rays up to the covering fraction
#CCD_cr = cosmics.addUpToFraction(self.information['coveringFraction'], limit=None)
CCD_cr_g = cosmics_g.addUpToFraction(self.information['coveringfraction']*self.information['exptime']/300.0, limit=None) CCD_cr_g = cosmics_g.addUpToFraction(self.information['coveringfraction']*self.information['exptime']/300.0, limit=None)
CCD_cr_r = cosmics_r.addUpToFraction(self.information['coveringfraction']*self.information['exptime']/300.0, limit=None) CCD_cr_r = cosmics_r.addUpToFraction(self.information['coveringfraction']*self.information['exptime']/300.0, limit=None)
CCD_cr_i = cosmics_i.addUpToFraction(self.information['coveringfraction']*self.information['exptime']/300.0, limit=None) CCD_cr_i = cosmics_i.addUpToFraction(self.information['coveringfraction']*self.information['exptime']/300.0, limit=None)
...@@ -2937,8 +2361,6 @@ class MCIsimulator(): ...@@ -2937,8 +2361,6 @@ class MCIsimulator():
self.image_r += CCD_cr_r self.image_r += CCD_cr_r
self.image_i += CCD_cr_i self.image_i += CCD_cr_i
#save cosmic ray image map
#self.cosmicMap = CCD_cr
if self.save_cosmicrays and self.cosmicRays: if self.save_cosmicrays and self.cosmicRays:
...@@ -3120,8 +2542,6 @@ class MCIsimulator(): ...@@ -3120,8 +2542,6 @@ class MCIsimulator():
y = np.round(cosmetics[:, 1]).astype(int) ## col number y = np.round(cosmetics[:, 1]).astype(int) ## col number
value = np.round(cosmetics[:, 2]).astype(int) value = np.round(cosmetics[:, 2]).astype(int)
#cosmetics_r=np.zeros((4636,235526))
self.log.info('Adding cosmetic defects to C2 channel:' ) self.log.info('Adding cosmetic defects to C2 channel:' )
for xc, yc, val in zip(x, y, value): for xc, yc, val in zip(x, y, value):
if 0 < xc < 4616 and 27 < (yc % 1499) <1179: if 0 < xc < 4616 and 27 < (yc % 1499) <1179:
...@@ -3817,10 +3237,6 @@ class MCIsimulator(): ...@@ -3817,10 +3237,6 @@ class MCIsimulator():
ofd_g.header['EPOCH'] =(np.float32(time2jd(t3.utcnow())), 'equinox of pointing RA and Dec') ofd_g.header['EPOCH'] =(np.float32(time2jd(t3.utcnow())), 'equinox of pointing RA and Dec')
ofd_g.header['EXPTIME']=(np.float32(self.information['exptime']), 'exposure time (s)') ofd_g.header['EXPTIME']=(np.float32(self.information['exptime']), 'exposure time (s)')
# ss1='HDU checksum'
# ss2='data unit checksum'
# ofd_g.header['CHECKSUM']=( data_time[:19] , ss1)
# ofd_g.header['DATASUM'] =( data_time[:19] , ss2)
ofd_g.header['CHECKSUM'] =(0, '0') ofd_g.header['CHECKSUM'] =(0, '0')
################################################ ################################################
########## finish header for 0 layer #################### ########## finish header for 0 layer ####################
...@@ -4039,10 +3455,7 @@ class MCIsimulator(): ...@@ -4039,10 +3455,7 @@ class MCIsimulator():
hdu_g.header['EPOCH'] =(np.float32(time2jd(t3.utcnow())), 'equinox of pointing RA and Dec') hdu_g.header['EPOCH'] =(np.float32(time2jd(t3.utcnow())), 'equinox of pointing RA and Dec')
hdu_g.header['CHECKSUM'] =(0, '0') hdu_g.header['CHECKSUM'] =(0, '0')
############################################################## ##############################################################
#####
# hdu_g.header['CHECKSUM'] =('','HDU checksum updated yyyy-mm-ddTHH:MM:SS')
# hdu_g.header['DATASUM'] =('','data unit checksum updated yyyy-mm-ddTHH:MM:SS')
################ finish MCI header for image layer in G channel ################ finish MCI header for image layer in G channel
##############################################################################################333 ##############################################################################################333
...@@ -4199,10 +3612,6 @@ class MCIsimulator(): ...@@ -4199,10 +3612,6 @@ class MCIsimulator():
ofd_r.header['EPOCH'] =(np.float32(time2jd(t3.utcnow())), 'equinox of pointing RA and Dec') ofd_r.header['EPOCH'] =(np.float32(time2jd(t3.utcnow())), 'equinox of pointing RA and Dec')
ofd_r.header['EXPTIME']=(np.float32(self.information['exptime']), 'exposure time (s)') ofd_r.header['EXPTIME']=(np.float32(self.information['exptime']), 'exposure time (s)')
# ss1='HDU checksum'
# ss2='data unit checksum'
# ofd_r.header['CHECKSUM']=( data_time[:19] , ss1)
# ofd_r.header['DATASUM'] =( data_time[:19] , ss2)
ofd_r.header['CHECKSUM'] =(0, '0') ofd_r.header['CHECKSUM'] =(0, '0')
################################################ ################################################
########## finish header for 0 layer #################### ########## finish header for 0 layer ####################
...@@ -4419,10 +3828,6 @@ class MCIsimulator(): ...@@ -4419,10 +3828,6 @@ class MCIsimulator():
hdu_r.header['EPOCH'] =(np.float32(time2jd(t3.utcnow())), 'equinox of pointing RA and Dec') hdu_r.header['EPOCH'] =(np.float32(time2jd(t3.utcnow())), 'equinox of pointing RA and Dec')
hdu_r.header['CHECKSUM'] =(0, '0') hdu_r.header['CHECKSUM'] =(0, '0')
############################################################## ##############################################################
#####
# hdu_r.header['CHECKSUM'] =('','HDU checksum updated yyyy-mm-ddTHH:MM:SS')
# hdu_r.header['DATASUM'] =('','data unit checksum updated yyyy-mm-ddTHH:MM:SS')
################ finish MCI header for image layer in G channel ################ finish MCI header for image layer in G channel
##############################################################################################333 ##############################################################################################333
...@@ -4580,10 +3985,6 @@ class MCIsimulator(): ...@@ -4580,10 +3985,6 @@ class MCIsimulator():
ofd_i.header['EPOCH'] =(np.float32(time2jd(t3.utcnow())), 'equinox of pointing RA and Dec') ofd_i.header['EPOCH'] =(np.float32(time2jd(t3.utcnow())), 'equinox of pointing RA and Dec')
ofd_i.header['EXPTIME']=(np.float32(self.information['exptime']), 'exposure time (s)') ofd_i.header['EXPTIME']=(np.float32(self.information['exptime']), 'exposure time (s)')
# ss1='HDU checksum'
# ss2='data unit checksum'
# ofd_i.header['CHECKSUM']=( data_time[:19] , ss1)
# ofd_i.header['DATASUM'] =( data_time[:19] , ss2)
ofd_i.header['CHECKSUM'] =(0, '0') ofd_i.header['CHECKSUM'] =(0, '0')
################################################ ################################################
########## finish header for 0 layer #################### ########## finish header for 0 layer ####################
...@@ -4805,8 +4206,6 @@ class MCIsimulator(): ...@@ -4805,8 +4206,6 @@ class MCIsimulator():
hdu_i.header['CHECKSUM'] =(0, '0') hdu_i.header['CHECKSUM'] =(0, '0')
############################################################## ##############################################################
#####
################ finish MCI header for image layer in G channel ################ finish MCI header for image layer in G channel
##############################################################################################333 ##############################################################################################333
...@@ -5336,7 +4735,7 @@ class MCIsimulator(): ...@@ -5336,7 +4735,7 @@ class MCIsimulator():
def runMCIsim(sourcein,configfile,dir_path, debug, iLoop): def runMCIsim(sourcein,configfile,dir_path, debug, iLoop):
print('路径Test：dir_path', dir_path) print('Path Test：dir_path', dir_path)
sim= dict() sim= dict()
sim[iLoop] = MCIsimulator(configfile) sim[iLoop] = MCIsimulator(configfile)
......
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