Skip to content
GitLab
Commit 0484bf5c authored by Yan Zhaojun's avatar Yan Zhaojun
Browse files

debug

parent f4148616
Pipeline #4084 passed with stage
in 0 seconds
Hide whitespace changes
Inline Side-by-side
csst_mci_sim/csst_mci_sim.py
View file @ 0484bf5c
......@@ -163,52 +163,52 @@ def transRaDec2D(ra, dec):
return np.array([x1, y1, z1])
###############################################################################
def flux2ill(wave, flux):
# erg/s/cm^2/A/arcsec^2 to W/m^2
# 1 W/m^2/sr/μm = 0.10 erg/cm^2/s/sr/A
# 1 sr = 1 rad^2 = 4.25452e10 arcsec^2
# 1 J/s = 1 W
# 1 J = 10^7 erg
# convert erg/s/cm^2/A/arcsec^2 to erg/s/cm^2/A/sr
flux1 = flux / (1/4.25452e10)
# convert erg/s/cm^2/A/sr to W/m^2/sr/um
flux2 = flux1 * 10
# 对接收面积积分，输出单位 W/m^2/nm
D = 2 # meter
f = 28 # meter
flux3 = flux2 * np.pi * D**2 / 4 / f**2 / 10**3
# # # u_lambda: W/m^2/nm
# # obj0:
# V = 73 * 1e-8 # W/m^2/nm
# Es = V * np.pi * D**2 / 4 / f**2 / 10**3
# c1 = 3.7418e-16
# c2 = 1.44e-2
# t = 5700
# # wave需要代入meter.
# wave0 = np.arange(380, 780) # nm
# delta_lamba0 = 1 # nm
# u_lambda = c1 / (wave0*1e-9)**5 / (np.exp(c2 / (wave0*1e-9) / t) - 1)
# f_lambda = (u_lambda / u_lambda[wave0 == 500]) * Es
# E0 = np.sum(f_lambda * 1)
# # plt.plot(wave, u_lambda)
# # plt.show()
# 对波长积分
f = interp1d(wave, flux3)
wave_interp = np.arange(3800, 7800)
flux3_interp = f(wave_interp)
# 输出单位 W/m^2
delta_lamba = 0.1 # nm
E = np.sum(flux3_interp * delta_lamba)
# pdb.set_trace()
return E
# def flux2ill(wave, flux):
# # erg/s/cm^2/A/arcsec^2 to W/m^2
# # 1 W/m^2/sr/μm = 0.10 erg/cm^2/s/sr/A
# # 1 sr = 1 rad^2 = 4.25452e10 arcsec^2
# # 1 J/s = 1 W
# # 1 J = 10^7 erg
# # convert erg/s/cm^2/A/arcsec^2 to erg/s/cm^2/A/sr
# flux1 = flux / (1/4.25452e10)
# # convert erg/s/cm^2/A/sr to W/m^2/sr/um
# flux2 = flux1 * 10
# # 对接收面积积分，输出单位 W/m^2/nm
# D = 2 # meter
# f = 28 # meter
# flux3 = flux2 * np.pi * D**2 / 4 / f**2 / 10**3
# # # # u_lambda: W/m^2/nm
# # # obj0:
# # V = 73 * 1e-8 # W/m^2/nm
# # Es = V * np.pi * D**2 / 4 / f**2 / 10**3
# # c1 = 3.7418e-16
# # c2 = 1.44e-2
# # t = 5700
# # # wave需要代入meter.
# # wave0 = np.arange(380, 780) # nm
# # delta_lamba0 = 1 # nm
# # u_lambda = c1 / (wave0*1e-9)**5 / (np.exp(c2 / (wave0*1e-9) / t) - 1)
# # f_lambda = (u_lambda / u_lambda[wave0 == 500]) * Es
# # E0 = np.sum(f_lambda * 1)
# # # plt.plot(wave, u_lambda)
# # # plt.show()
# # 对波长积分
# f = interp1d(wave, flux3)
# wave_interp = np.arange(3800, 7800)
# flux3_interp = f(wave_interp)
# # 输出单位 W/m^2
# delta_lamba = 0.1 # nm
# E = np.sum(flux3_interp * delta_lamba)
# # pdb.set_trace()
# return E
################################################################
def ill2flux(E,path):
......@@ -270,36 +270,36 @@ def earth_angle(time_jd, x_sat, y_sat, z_sat, ra_obj, dec_obj):
#################################################3
def MCIinformation():
"""
Returns a dictionary describing MCI CCD. The following information is provided (id: value - reference)::
# 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)
# 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
"""
# :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})
# #########################################################################################################
# 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
# return out
#######################################
def CCDnonLinearityModel(data, beta=6e-7):
......@@ -400,29 +400,29 @@ class StrayLight(object):
###############################################################################
###############################################################################
def make_c_coor(fov, step):
"""
Draw the mesh grids for a fov * fov box with given step .
"""
# def make_c_coor(fov, step):
# """
# Draw the mesh grids for a fov * fov box with given step .
# """
nc=int(fov/step)
# nc=int(fov/step)
bs=fov
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
# bs=fov
# 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 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')
# 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):
......@@ -515,14 +515,14 @@ def deg2HMS(ra0, dec0):
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 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 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)
# def v_disp(sigstar, mstar, mag):
# return sigstar * 10**(0.4*(mstar-mag)/3.57)
###############################################################################
#####################################################################################
......@@ -561,56 +561,56 @@ def distortField(ra, dec, ch):
return distortRa, distortDec
#####################################################################################
def world_to_pixel(sra,
sdec,
rotsky,
tra,
tdec,
x_center_pixel,
y_center_pixel,
pixelsize=0.05):
"""_summary_
Parameters
----------
sra : np.array
star ra such as np.array([22,23,24]),unit:deg
sdec : np.array
stat dec such as np.array([10,20,30]),unit:deg
rotsky : float
rotation angel of the telescope,unit:deg
tra : float
telescope pointing ra,such as 222.2 deg
rdec : float
telescope pointing dec,such as 33.3 deg
x_center_pixel : float
x refer point of MCI,usually the center of the CCD,which start from (0,0)
y_center_pixel : float
y refer point of MCI,usually the center of the CCD,which start from(0,0)
pixelsize : float
pixelsize for MCI ccd, default :0.05 arcsec / pixel
Returns
-------
pixel_position:list with two np.array
such as [array([124.16937605, 99. , 99. ]),
array([ 97.52378661, 99. , 100.50014483])]
"""
theta_r = rotsky / 180 * np.pi
w = WCS(naxis=2)
w.wcs.crpix = [x_center_pixel + 1, y_center_pixel + 1] # pixel from (1, 1)
w.wcs.cd = np.array([[
-np.cos(-theta_r) * pixelsize / 3600,
np.sin(-theta_r) * pixelsize / 3600
],
[
np.sin(-theta_r) * pixelsize / 3600,
np.cos(-theta_r) * pixelsize / 3600
]])
w.wcs.crval = [tra, tdec]
w.wcs.ctype = ["RA---TAN", "DEC--TAN"]
pixel_position = w.all_world2pix(sra, sdec, 0)
return pixel_position
# def world_to_pixel(sra,
# sdec,
# rotsky,
# tra,
# tdec,
# x_center_pixel,
# y_center_pixel,
# pixelsize=0.05):
# """_summary_
# Parameters
# ----------
# sra : np.array
# star ra such as np.array([22,23,24]),unit:deg
# sdec : np.array
# stat dec such as np.array([10,20,30]),unit:deg
# rotsky : float
# rotation angel of the telescope,unit:deg
# tra : float
# telescope pointing ra,such as 222.2 deg
# rdec : float
# telescope pointing dec,such as 33.3 deg
# x_center_pixel : float
# x refer point of MCI,usually the center of the CCD,which start from (0,0)
# y_center_pixel : float
# y refer point of MCI,usually the center of the CCD,which start from(0,0)
# pixelsize : float
# pixelsize for MCI ccd, default :0.05 arcsec / pixel
# Returns
# -------
# pixel_position:list with two np.array
# such as [array([124.16937605, 99. , 99. ]),
# array([ 97.52378661, 99. , 100.50014483])]
# """
# theta_r = rotsky / 180 * np.pi
# w = WCS(naxis=2)
# w.wcs.crpix = [x_center_pixel + 1, y_center_pixel + 1] # pixel from (1, 1)
# w.wcs.cd = np.array([[
# -np.cos(-theta_r) * pixelsize / 3600,
# np.sin(-theta_r) * pixelsize / 3600
# ],
# [
# np.sin(-theta_r) * pixelsize / 3600,
# np.cos(-theta_r) * pixelsize / 3600
# ]])
# w.wcs.crval = [tra, tdec]
# w.wcs.ctype = ["RA---TAN", "DEC--TAN"]
# pixel_position = w.all_world2pix(sra, sdec, 0)
# return pixel_position
###############################################################################
......@@ -635,39 +635,39 @@ def cal_pos(center_ra, center_dec, rotTel, rotSky, sRa, sDec):
#####################################################################################
def krebin(a, sample):
""" Fast Rebinning with flux conservation
New shape must be an integer divisor of the current shape.
This algorithm is much faster than rebin_array
Parameters
----------
a : array_like
input array
sample : rebinned sample 2 or 3 or 4 or other int value
"""
# Klaus P's fastrebin from web
size=a.shape
shape=np.zeros(2,dtype=int)
shape[0]=int(size[0]/sample)
shape[1]=int(size[1]/sample)
sh = shape[0], a.shape[0] // shape[0], shape[1], a.shape[1] // shape[1]
return a.reshape(sh).sum(-1).sum(1)
# def krebin(a, sample):
# """ Fast Rebinning with flux conservation
# New shape must be an integer divisor of the current shape.
# This algorithm is much faster than rebin_array
# Parameters
# ----------
# a : array_like
# input array
# sample : rebinned sample 2 or 3 or 4 or other int value
# """
# # Klaus P's fastrebin from web
# size=a.shape
# shape=np.zeros(2,dtype=int)
# shape[0]=int(size[0]/sample)
# shape[1]=int(size[1]/sample)
# sh = shape[0], a.shape[0] // shape[0], shape[1], a.shape[1] // shape[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)
# 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]
# if z>=n:
# return b
# else:
# return b[:z]
#########################################################
def centroid(data):
......@@ -744,270 +744,270 @@ def dft2(ary, Q, samples, shift=None):
return out
##############################################################################
def idft2(ary, Q, samples, shift=None):
"""Compute the two dimensional inverse Discrete Fourier Transform of a matrix.
# 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.
# 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
# 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
# """
# # 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
# 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))
# 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
# ###############################################################################
# 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)
# ###################################################################
# # 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
# out = Eout_rev @ ary @ Ein_rev
return out
# return out
###############################################################################
###############################################################################
def opd2psf(opd,cwave,oversampling):
'''
calculate psf from opd data with defined wavelength
# 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.
# 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
# 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
# opd=opd*1e9; ## convert opd unit of meter to nm
zoomopd=ndimage.zoom(opd, 2, order=1)
# zoomopd=ndimage.zoom(opd, 2, order=1)
m,n=np.shape(zoomopd);
clambda=cwave;
# m,n=np.shape(zoomopd);
# clambda=cwave;
wfe=zoomopd/clambda; #%% the main telescope aberration ,in wavelength;
# wfe=zoomopd/clambda; #%% the main telescope aberration ,in wavelength;
Q=clambda*1e-3*focLength/D/pixel
# Q=clambda*1e-3*focLength/D/pixel
pupil=np.zeros((m,m));
# pupil=np.zeros((m,m));
pupil[abs(zoomopd)>0]=1;
#idx=pupil>0;
# pupil[abs(zoomopd)>0]=1;
# #idx=pupil>0;
Nt=512
# Nt=512
phase=2*np.pi*wfe; #% wavefront phase in radians
# phase=2*np.pi*wfe; #% wavefront phase in radians
#%generalized pupil function of channel1;
pk=pupil*np.exp(cmath.sqrt(-1)*(phase))
# #%generalized pupil function of channel1;
# pk=pupil*np.exp(cmath.sqrt(-1)*(phase))
pf=dft2(pk, Q, Nt)
# pf=dft2(pk, Q, Nt)
pf2 = abs(pf)**2
psfout=pf2/pf2.sum()
# pf2 = abs(pf)**2
# psfout=pf2/pf2.sum()
cx,cy=centroid(psfout)
# cx,cy=centroid(psfout)
psf=ndimage.shift(psfout,[Nt/2-1-cy, Nt/2-1-cx],order=1, mode='nearest' )
# psf=ndimage.shift(psfout,[Nt/2-1-cy, Nt/2-1-cx],order=1, mode='nearest' )
return psf
# return psf
########################################################################
def cal_PSF_new(channel,wfetimes,oversampling):
# filed postion: Ra=ys1, Dec=ys2
# 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
# ############ 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;
# 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]
# #wavelength=[255,337,419,501,583,665,747,829,911,1000]
cwave=dict()
# 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
# 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)
# 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)
# 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
# psf=dict()
# fx=dict()
# fy=dict()
# fx[channel]=ffx*3600-cfx
fy[channel]=ffy*3600-cfy
# fy[channel]=ffy*3600-cfy
Nt=512
# Nt=512
psf[channel]=np.zeros((len(n),Nt,Nt))
# 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)
# 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=fits.PrimaryHDU(psf[channel])
hdu1.header.append(('channel', channel, 'MCI PSF channel'))
# hdu1.header.append(('channel', channel, 'MCI PSF channel'))
hdu1.header['pixsize']=(0.05/oversampling, 'PSF pixel size in arcsec')
# hdu1.header['pixsize']=(0.05/oversampling, 'PSF pixel size in arcsec')
hdu1.header['wfetimes']=(wfetimes, 'wavefront error magnification to CSST primary wavefront ')
# 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)
# 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=fits.ImageHDU(fx[channel])
hdu2.header.append(('field_X', 'arcsec'))
# hdu2.header.append(('field_X', 'arcsec'))
hdu3=fits.ImageHDU(fy[channel])
# hdu3=fits.ImageHDU(fy[channel])
hdu3.header.append(('field_Y', 'arcsec'))
# hdu3.header.append(('field_Y', 'arcsec'))
newd=fits.HDUList([hdu1,hdu2,hdu3])
# newd=fits.HDUList([hdu1,hdu2,hdu3])
PSFfilename=self.information['dir_path']+'MCI_inputData/PSF/PSF_'+channel+'.fits'
# PSFfilename=self.information['dir_path']+'MCI_inputData/PSF/PSF_'+channel+'.fits'
newd.writeto(PSFfilename,overwrite=True)
# newd.writeto(PSFfilename,overwrite=True)
return
# return
#############################################################################
def cal_Filter_PSF(wfetimes):
# filed postion: Ra=ys1, Dec=ys2
# 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
# ############ 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;
# 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])
# 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()
# 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)
# 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
# psf=dict()
# fx=dict()
# fy=dict()
# fx=ffx*3600-cfx
# fy=ffy*3600-cfy
# Nt=512
for filterk in filterP.keys():
filtername=filterk
# for filterk in filterP.keys():
# filtername=filterk
## print(filtername)
# ## print(filtername)
filterwaveC=filterP[filterk]['Clambda']
filterFwhm=filterP[filterk]['FWHM']
# filterwaveC=filterP[filterk]['Clambda']
# filterFwhm=filterP[filterk]['FWHM']
fdis=abs(filterwaveC-wavelist)
# fdis=abs(filterwaveC-wavelist)
findex=np.argsort(fdis)
# findex=np.argsort(fdis)
waven=findex[0]
# waven=findex[0]
psf[filtername]=dict()
# 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)
# 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
# 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)
# psfd=opd2psf(wfetimes*opd, wavek, oversampling)
psf[filtername]['psf_mat'][:,:,kk,ii]=psfd[:,:]
# psf[filtername]['psf_mat'][:,:,kk,ii]=psfd[:,:]
if ii==0:
psf[filtername]['psf_iwave'][kk]=wavek
# if ii==0:
# psf[filtername]['psf_iwave'][kk]=wavek
np.save(self.information['dir_path']+'mci_sim_result/'+filtername+'_PSF.npy', psf[filtername])
# np.save(self.information['dir_path']+'mci_sim_result/'+filtername+'_PSF.npy', psf[filtername])
return
# return
......@@ -1169,13 +1169,9 @@ class MCIsimulator():
# parallelTrapfile =self.information['dir_path']+'MCI_inputData/data/cdm_euclid_parallel.dat',
# serialTrapfile =self.information['dir_path']+'MCI_inputData/data/cdm_euclid_serial.dat',
mode='same'))
####################################################
####################################################
###############################################################################
def readConfigs(self,simnumber,source):
"""
......@@ -1192,7 +1188,7 @@ class MCIsimulator():
################################################################################################3
###########################################################################
def processConfigs(self):
"""
......@@ -1287,16 +1283,16 @@ class MCIsimulator():
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 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):
......@@ -1316,26 +1312,26 @@ class MCIsimulator():
#########################################################################################################
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.
# 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.
# 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.
# .. 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
# :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)
# :return: smoothed image array
# :rtype: ndarray
# """
# return ndimage.filters.gaussian_filter(image, sigma)
def _loadGhostModel(self):
......@@ -1398,7 +1394,7 @@ class MCIsimulator():
now=datetime.utcnow()
#data_time=now.strftime("%Y-%m-%d-%H-%M-%S")
result_day=now.strftime("%Y-%m-%d")
#self.result_path='../MCI_simData_'+result_day
......@@ -1722,109 +1718,109 @@ class MCIsimulator():
#################### cal_PSF_array #################################
# if self.save_starpsf==100:
if self.save_starpsf==100:
# self.log.info('calculate and save star PSF data......')
self.log.info('calculate and save star PSF data......')
# primary_g=fits.PrimaryHDU()
# PSF_g=fits.HDUList([primary_g])
primary_g=fits.PrimaryHDU()
PSF_g=fits.HDUList([primary_g])
# primary_r=fits.PrimaryHDU()
# PSF_r=fits.HDUList([primary_r])
primary_r=fits.PrimaryHDU()
PSF_r=fits.HDUList([primary_r])
# primary_i=fits.PrimaryHDU()
# PSF_i=fits.HDUList([primary_i])
primary_i=fits.PrimaryHDU()
PSF_i=fits.HDUList([primary_i])
# fov=0.05*min(self.information['xsize'],self.information['ysize'])
fov=0.05*min(self.information['xsize'],self.information['ysize'])
# dec_arr,ra_arr=make_c_coor(fov,25)
dec_arr,ra_arr=make_c_coor(fov,25)
# k=0;
k=0;
# #######################################################################
# for ii in range(len(ra_arr)):
# for jj in range(len(ra_arr)):
#######################################################################
for ii in range(len(ra_arr)):
for jj in range(len(ra_arr)):
# ##################################################################
##################################################################
# galRa = ra_arr[ii,jj] +center_ra*3600 # ra of PFS, arcsecond
# galDec = dec_arr[ii,jj]+center_dec*3600 # dec of PSF, arcsecond
galRa = ra_arr[ii,jj] +center_ra*3600 # ra of PFS, arcsecond
galDec = dec_arr[ii,jj]+center_dec*3600 # dec of PSF, arcsecond
# fsx,fsy=cal_pos(center_ra,center_dec, rotTelPos, rotSkyPos, galRa, galDec)
fsx,fsy=cal_pos(center_ra,center_dec, rotTelPos, rotSkyPos, galRa, galDec)
# ############# do field distottion ##########
############# do field distottion ##########
# if self.distortion:
if self.distortion:
# for i in range(3):
# ch=channel[i]
# fpx,fpy=distortField(fsx, fsy, ch)
for i in range(3):
ch=channel[i]
fpx,fpy=distortField(fsx, fsy, ch)
# else:
else:
# fpx=fsx
# fpy=fsy
fpx=fsx
fpy=fsy
# ###################################################################
###################################################################
# psf=dict()
psf=dict()
# for i in range(3):
for i in range(3):
# ch=channel[i]
ch=channel[i]
# psfmat=self.get_PSF(fpx, fpy, ch)
psfmat=self.get_PSF(fpx, fpy, ch)
# temp=0
temp=0
# for iwave in range(7):
# temp=temp+1/7.0*psfmat[:,:,iwave]
for iwave in range(7):
temp=temp+1/7.0*psfmat[:,:,iwave]
# temp=temp/temp.sum()
temp=temp/temp.sum()
# ####rotate the PSF data
# if abs(theta.deg)>0:
# psf[ch]=ndimage.rotate(temp, theta.deg, order=1, reshape=False) # here we choose reshape=False, the rotated image will
# else:
# psf[ch]=temp
# ###################################
####rotate the PSF data
if abs(theta.deg)>0:
psf[ch]=ndimage.rotate(temp, theta.deg, order=1, reshape=False) # here we choose reshape=False, the rotated image will
else:
psf[ch]=temp
###################################
# hdu_g=fits.ImageHDU(psf['g'])
# hdu_g.header['ra'] = ra_arr[ii,jj]
# hdu_g.header['dec'] = dec_arr[ii,jj]
# hdu_g.header['rot_deg']=-theta.deg
# PSF_g.append(hdu_g)
# ###################################
# hdu_r=fits.ImageHDU(psf['r'])
# hdu_r.header['ra'] = ra_arr[ii,jj]
# hdu_r.header['dec'] = dec_arr[ii,jj]
# hdu_r.header['rot_deg']=-theta.deg
# PSF_r.append(hdu_r)
# ###################################
# hdu_i=fits.ImageHDU(psf['i'])
# hdu_i.header['ra'] = ra_arr[ii,jj]
# hdu_i.header['dec']= dec_arr[ii,jj]
# hdu_i.header['rot_deg']=-theta.deg
# PSF_i.append(hdu_i)
# ###################################
# del hdu_g
# del hdu_r
# del hdu_i
# k=k+1
hdu_g=fits.ImageHDU(psf['g'])
hdu_g.header['ra'] = ra_arr[ii,jj]
hdu_g.header['dec'] = dec_arr[ii,jj]
hdu_g.header['rot_deg']=-theta.deg
PSF_g.append(hdu_g)
###################################
hdu_r=fits.ImageHDU(psf['r'])
hdu_r.header['ra'] = ra_arr[ii,jj]
hdu_r.header['dec'] = dec_arr[ii,jj]
hdu_r.header['rot_deg']=-theta.deg
PSF_r.append(hdu_r)
###################################
hdu_i=fits.ImageHDU(psf['i'])
hdu_i.header['ra'] = ra_arr[ii,jj]
hdu_i.header['dec']= dec_arr[ii,jj]
hdu_i.header['rot_deg']=-theta.deg
PSF_i.append(hdu_i)
###################################
del hdu_g
del hdu_r
del hdu_i
k=k+1
# ############################################
# file_g=self.result_path+'/PSF_Data/'+'PSF_sim_No.'+str(self.information['simnumber'])+'_C1.fits'
# PSF_g.writeto(file_g,overwrite=True)
############################################
file_g=self.result_path+'/PSF_Data/'+'PSF_sim_No.'+str(self.information['simnumber'])+'_C1.fits'
PSF_g.writeto(file_g,overwrite=True)
# file_r=self.result_path+'/PSF_Data/'+'PSF_sim_No.'+str(self.information['simnumber'])+'_C2.fits'
# PSF_r.writeto(file_r,overwrite=True)
file_r=self.result_path+'/PSF_Data/'+'PSF_sim_No.'+str(self.information['simnumber'])+'_C2.fits'
PSF_r.writeto(file_r,overwrite=True)
# file_i=self.result_path+'/PSF_Data/'+'PSF_sim_No.'+str(self.information['simnumber'])+'_C3.fits'
# PSF_i.writeto(file_i,overwrite=True)
file_i=self.result_path+'/PSF_Data/'+'PSF_sim_No.'+str(self.information['simnumber'])+'_C3.fits'
PSF_i.writeto(file_i,overwrite=True)
# del PSF_g
# del PSF_r
# del PSF_i
del PSF_g
del PSF_r
del PSF_i
########## finish save PSF fits ##########################################
############################################################################
......@@ -1919,9 +1915,6 @@ class MCIsimulator():
f2 = interp1d(earthshine_wave0, earthshine_flux0)
earthshine_mci = f2(wave_mci)
# df = pd.DataFrame({'wave': wave_ifs,
# 'zodi': zodi_ifs,
# 'earthshine': earthshine_ifs})
self.zodiacal_wave = wave_mci # in A
self.zodiacal_flux = zodi_mci
......@@ -1929,9 +1922,7 @@ class MCIsimulator():
self.earthshine_wave = wave_mci # A
self.earthshine_flux = earthshine_mci
########################################################################################
#######################################################################################
#########################################################################################
self.cal_sky_noise()
......@@ -1952,11 +1943,8 @@ class MCIsimulator():
self.information['CD2_1']= np.sin(-theta)*self.information['pixel_size']/3600.0 ####
self.information['CD2_2']= np.cos(-theta)*self.information['pixel_size']/3600.0 ####
#######################################################################
if self.TianceEffect:
#######################################################################
if self.TianceEffect:
ra_list = self.star['ra_gaia'].tolist()
dec_list = self.star['dec_gaia'].tolist()
pmra_list = self.star['pmra_gaia'].tolist()
......@@ -1993,8 +1981,6 @@ class MCIsimulator():
st_magi=[]
st_magz=[]
#################### generate star image ##########
if self.debug:
......@@ -2116,12 +2102,8 @@ class MCIsimulator():
ugriz = np.array([umag, gmag, rmag, imag, zmag])
star_flux = sed.Model_Galaxy_SED(wave, ugriz, redshift, t, self.information['dir_path'])
else:
umag = binary_star[j-nsrcs, 2]
gmag = binary_star[j-nsrcs, 3]
rmag = binary_star[j-nsrcs, 4]
......@@ -2143,9 +2125,7 @@ class MCIsimulator():
##cal_SED_photons(self, flux_arr):
intscales,PSF_eff_weight=self.cal_SED_photons(star_flux)
gx=dict()
gy=dict()
......@@ -2180,7 +2160,6 @@ class MCIsimulator():
st_magi.append(imag)
st_magz.append(zmag)
## self.log.info('begin cal PSF')
######################################################################
psf=dict()
......@@ -2227,9 +2206,6 @@ class MCIsimulator():
cx0,cy0=centroid(stamp_img.array)
if self.appFatt :
### apply treering and bright fatter and diffusion;
......@@ -2266,8 +2242,6 @@ class MCIsimulator():
photons.addTo(final_image[ch])
###############################
###### debug code
# fi= galsim.ImageF(int(self.information['xsize']), int(self.information['ysize']))
......@@ -2482,21 +2456,16 @@ class MCIsimulator():
self.information['ra_obj'] =self.information['star_ra']
self.information['dec_obj'] =self.information['star_dec']
#####################################################################
# self.earthshine(self.earthshine_theta)
# self.zodiacal(self.information['star_ra'], self.information['star_dec'], self.dt.strftime("%Y-%m-%d"))
####################################################################
##################################################################################
ra = self.information['ra_pnt0']
dec = self.information['dec_pnt0']
time_jd=time2jd(self.dt)
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
wave0, zodi0 = self.zodiacal(self.information['ra_pnt0'], self.information['dec_pnt0'], self.TianCe_day) # erg/s/cm^2/A/arcsec^2
# EarthShine from straylight
sl = StrayLight(self.information['dir_path'], jtime=time_jd, sat=np.array([x_sat, y_sat, z_sat]),
......@@ -2632,7 +2601,7 @@ class MCIsimulator():
print('k2=',k2)
#
filename=self.information['dir_path']+'MCI_inputData/galaxy_Input/Lens_SED_IMG_0.025_230626/Lens_img_cut_IMG_'+str(k2+1)+'.fits'
filename=self.information['dir_path']+'MCI_inputData/galaxy_Input/noLens_SED_IMG_0.025_230626/Lens_img_cut_IMG_'+str(k2+1)+'.fits'
self.log.info('galaxy_Input image path is: %s' %(filename))
......@@ -2643,7 +2612,7 @@ class MCIsimulator():
srcs_cat=fits.open(filename)
#### load galaxy SED fitsfile ###
filename=self.information['dir_path']+'MCI_inputData/galaxy_Input/Lens_SED_IMG_0.025_230626/Lens_img_cut_SED_'+str(k2+1)+'.fits'
filename=self.information['dir_path']+'MCI_inputData/galaxy_Input/noLens_SED_IMG_0.025_230626/Lens_img_cut_SED_'+str(k2+1)+'.fits'
srcs_sed=fits.open(filename)
self.log.info('galaxy_Input SED path is: %s' %(filename))
......@@ -2655,10 +2624,10 @@ class MCIsimulator():
for kkk in range(1,len(srcs_cat)):
t1=srcs_cat[kkk].header['new_ra'] -self.information['gal_ra'] +self.information['star_ra']
ra_list.append(t1)
ra_list.append(float(t1))
t2=srcs_cat[kkk].header['new_dec'] -self.information['gal_dec'] +self.information['star_dec']
dec_list.append(t2)
dec_list.append(float(t2))
if self.TianceEffect:
......@@ -2669,7 +2638,7 @@ class MCIsimulator():
rv_list = [0.0 for i in range(len(ra_list))]
################################################
newRa, newDec = shao.onOrbitObsPosition(ra_list, dec_list, pmra_list, \
newRa, newDec = shao.onOrbitObsPosition(self.information['dir_path'],ra_list, dec_list, pmra_list, \
pmdec_list, rv_list, parallax_list, len(ra_list), \
self.information['pos_x'], self.information['pos_y'], self.information['pos_z'], self.information['velocity_x'],self.information['velocity_y'],self.information['velocity_z'], "J2000", self.TianCe_day, self.TianCe_exp_start)
else:
......
csst_mci_sim/mci_data/mci_all_9K.config
View file @ 0484bf5c
......@@ -62,7 +62,7 @@ TianceEffect = yes
intscale = yes
ghosts = no
ghosts = yes
shutterEffect = yes
......@@ -70,7 +70,7 @@ flatfieldM = yes
PRNUeffect = yes
appFatt = no
appFatt = yes
sky_shift_rot = yes
......@@ -80,9 +80,9 @@ sim_star = yes
sim_galaxy = yes
save_starpsf = no
save_starpsf = yes
save_cosmicrays = no
save_cosmicrays = yes
##############################################
##############################################
......
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