From 9e7cecd10ad1f770a1bf4a6102f92070168dff3c Mon Sep 17 00:00:00 2001 From: yan Date: Sat, 11 May 2024 14:40:25 +0800 Subject: [PATCH] debug --- csst_mci_sim/csst_mci_sim.py | 633 +---------------------------------- 1 file changed, 16 insertions(+), 617 deletions(-) diff --git a/csst_mci_sim/csst_mci_sim.py b/csst_mci_sim/csst_mci_sim.py index b214047..69ca23f 100644 --- a/csst_mci_sim/csst_mci_sim.py +++ b/csst_mci_sim/csst_mci_sim.py @@ -55,79 +55,36 @@ Contact Information: zhaojunyan@shao.ac.cn """ import galsim import pandas as pd -from scipy.integrate import simps - +#from scipy.integrate import simps from datetime import datetime, timedelta from astropy.time import Time from astropy.coordinates import get_sun - -# import scipy.io as sio -# import cmath - import numpy as np - from tqdm import tqdm - from astropy.table import Table - -# from astropy.wcs import WCS as WCS - from astropy.io import fits - from astropy import units as u - import os, sys, math - import configparser as ConfigParser - -# from optparse import OptionParser - -# from matplotlib import pyplot as plt - +from matplotlib import pyplot as plt from scipy import ndimage - sys.path.append('./csst_mci_sim') -##sys.path.append('..') - from CTI import CTI from support import logger as lg from support import cosmicrays - from support import shao - from support import sed -#from shao import onOrbitObsPosition - from support import MCIinstrumentModel - from joblib import Parallel, delayed -#import multiprocessing - from astropy.coordinates import SkyCoord - from scipy import interpolate - from scipy.signal import fftconvolve - import pandas - -# sys.path.append('../MCI_inputData/TianCe') -# sys.path.append('../MCI_inputData/SED_Code') - - -# set the folder -FOLDER ='../' - - -##################################################################################### -############################################################################### import ctypes - import astropy.coordinates as coord - from scipy.interpolate import interp1d - +########################### functions ######################### def transRaDec2D(ra, dec): # radec转为竞天程序里的ob, 赤道坐标系下的笛卡尔三维坐标xyz. 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): angle = 180 - angle_a - angle_b.degree 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): @@ -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 def time2mjd(dateT): @@ -444,18 +360,10 @@ def deg2HMS(ra0, dec0): hhmmss=ss[0:2]+ss[3:5]+ra_ss 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 def distortField(ra, dec, ch): @@ -586,18 +494,6 @@ def cal_pos(center_ra, center_dec, rotTel, rotSky, sRa, sDec): # 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): @@ -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 @@ -1145,8 +722,6 @@ class MCIsimulator(): #ghost ratio can be in engineering format, so getfloat does not capture it... self.information['ghostratio'] = float(self.config.get(self.section, 'ghostRatio')) - - ################################################################################################### ################################################################################################## @@ -1201,27 +776,12 @@ class MCIsimulator(): 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 - ######################################################################################### - # 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): """ @@ -1236,32 +796,8 @@ class MCIsimulator(): self.image_i=np.zeros((self.information['ysize'], self.information['xsize']), dtype=float) 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): """ Reads in a ghost model from a FITS file and stores the data to self.ghostModel. @@ -1316,16 +852,12 @@ class MCIsimulator(): ##print('print information:', self.information) - ########################################################################### now=datetime.utcnow() #data_time=now.strftime("%Y-%m-%d-%H-%M-%S") result_day=now.strftime("%Y-%m-%d") - - - if self.information['dir_path'] =='/nfsdata/share/simulation-unittest/mci_sim/': self.result_path = self.information['dir_path'] +'mci_sim_result/'+result_day else: @@ -1338,7 +870,6 @@ class MCIsimulator(): else: self.result_path = '/data/mcisimdata/'+result_day - if os.path.isdir(self.result_path)==False: os.mkdir(self.result_path) @@ -1358,8 +889,7 @@ class MCIsimulator(): 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------------') for key, value in self.information.items(): @@ -1376,16 +906,11 @@ class MCIsimulator(): #### load telescope efficiency data - self.tel_eff=np.load(self.information['dir_path']+'MCI_inputData/tel_eff/tel_eff.npy',allow_pickle=True).item() #### load MCI filter data self.filterP=np.load(self.information['dir_path']+'MCI_inputData/MCI_filters/mci_filterPWTC.npy',allow_pickle=True).item() - - - - return ################################################################################ @@ -1521,14 +1046,6 @@ class MCIsimulator(): 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 @@ -1607,13 +1124,6 @@ class MCIsimulator(): 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) @@ -1753,9 +1263,7 @@ class MCIsimulator(): ########## finish save PSF fits ########################################## ############################################################################ - ###### add binary stars to the star catlog ####### - nsrcs=len(self.star['ra_gaia']) self.log.info('load star catlog successfully') @@ -1790,8 +1298,6 @@ class MCIsimulator(): ####################################################################### nlayccd = 0 - - ############################################ ra = self.information['ra_pnt0'] @@ -1803,9 +1309,6 @@ class MCIsimulator(): 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 self.log.info('dir_path:') @@ -1898,10 +1401,7 @@ class MCIsimulator(): st_magz=[] #################### generate star image ########## - - - - + for j in tqdm(range(nsrcs)): ### nsrcs length if j==0: @@ -1913,9 +1413,7 @@ class MCIsimulator(): starRa = 3600.0*( newRa[j] ) # ra of star, arcsecond starDec = 3600.0*( newDec[j] ) # dec of star, arcsecond - ################################################################### - ################################################################### 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 @@ -1928,7 +1426,6 @@ class MCIsimulator(): nlayccd=nlayccd+1 - if self.debug: if nlayccd>1: break @@ -1976,8 +1473,6 @@ class MCIsimulator(): ####### use SED_code to generate star SED ####### # SED of j-th star - - umag = self.star['umag'][j] gmag = self.star['gmag'][j] rmag = self.star['rmag'][j] @@ -1987,7 +1482,6 @@ class MCIsimulator(): # Input redshift redshift = 0 - # Loading galaxy SED template t=sed.Gal_Temp(self.information['dir_path']) # Calculating the magnitude (u, g, r, i, z) of each template @@ -1996,8 +1490,6 @@ class MCIsimulator(): ugriz = np.array([umag, gmag, rmag, imag, zmag]) star_flux = sed.Model_Galaxy_SED(wave, ugriz, redshift, t, self.information['dir_path']) - - ##cal_SED_photons(self, flux_arr): intscales,PSF_eff_weight=self.cal_SED_photons(star_flux) @@ -2019,9 +1511,7 @@ class MCIsimulator(): if nlayccd %10 == 0: self.log.info('Star number on CCD is = %i' % (nlayccd)) - - - + ###### st_phot_C1.append(intscales['g']) st_phot_C2.append(intscales['r']) @@ -2114,9 +1604,7 @@ class MCIsimulator(): cx0,cy0=centroid(image.array) - - #self.log.info('finish appfat ........') - + ############################################ ############################################ photons.x=photons.x-cx0*0.025 @@ -2176,8 +1664,6 @@ class MCIsimulator(): tab['st_magi']=np.asarray(st_magi) tab['st_magz']=np.asarray(st_magz) - - ###################################33 @@ -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): """ @@ -2820,8 +2249,6 @@ class MCIsimulator(): ############################################################################### - -############################################################################### ############################################################################### @@ -2924,9 +2351,6 @@ class MCIsimulator(): 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) - #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_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) @@ -2937,8 +2361,6 @@ class MCIsimulator(): self.image_r += CCD_cr_r self.image_i += CCD_cr_i - #save cosmic ray image map - #self.cosmicMap = CCD_cr if self.save_cosmicrays and self.cosmicRays: @@ -3120,8 +2542,6 @@ class MCIsimulator(): y = np.round(cosmetics[:, 1]).astype(int) ## col number value = np.round(cosmetics[:, 2]).astype(int) - #cosmetics_r=np.zeros((4636,235526)) - self.log.info('Adding cosmetic defects to C2 channel:' ) for xc, yc, val in zip(x, y, value): if 0 < xc < 4616 and 27 < (yc % 1499) <1179: @@ -3817,10 +3237,6 @@ class MCIsimulator(): 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)') - # 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') ################################################ ########## finish header for 0 layer #################### @@ -4039,10 +3455,7 @@ class MCIsimulator(): 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'] =('','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 ##############################################################################################333 @@ -4199,10 +3612,6 @@ class MCIsimulator(): 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)') - # 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') ################################################ ########## finish header for 0 layer #################### @@ -4419,10 +3828,6 @@ class MCIsimulator(): 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'] =('','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 ##############################################################################################333 @@ -4580,10 +3985,6 @@ class MCIsimulator(): 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)') - # 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') ################################################ ########## finish header for 0 layer #################### @@ -4805,8 +4206,6 @@ class MCIsimulator(): hdu_i.header['CHECKSUM'] =(0, '0') ############################################################## - ##### - ################ finish MCI header for image layer in G channel ##############################################################################################333 @@ -5336,7 +4735,7 @@ class MCIsimulator(): def runMCIsim(sourcein,configfile,dir_path, debug, iLoop): - print('路径Test：dir_path', dir_path) + print('Path Test：dir_path', dir_path) sim= dict() sim[iLoop] = MCIsimulator(configfile) -- GitLab