diff --git a/csst_mci_sim/__pycache__/csst_mci_sim.cpython-311.pyc b/csst_mci_sim/__pycache__/csst_mci_sim.cpython-311.pyc index 12d60db044ecf6362e20dd28e52394913b24d0ed..6c1399c8b032d263a2aa6d974a1c89c0e29d82f2 100644 Binary files a/csst_mci_sim/__pycache__/csst_mci_sim.cpython-311.pyc and b/csst_mci_sim/__pycache__/csst_mci_sim.cpython-311.pyc differ diff --git a/csst_mci_sim/csst_mci_sim.py b/csst_mci_sim/csst_mci_sim.py index 777ab42a6352de14d40d711673156c07f14d05a4..72d50182aa6f94e0f565d8977be5f142064e1401 100644 --- a/csst_mci_sim/csst_mci_sim.py +++ b/csst_mci_sim/csst_mci_sim.py @@ -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: diff --git a/csst_mci_sim/mci_data/mci_all_9K.config b/csst_mci_sim/mci_data/mci_all_9K.config index 885a24df0e8a0bf826742c24f122ff9f65e30b6e..5b45b7e45c12281f373e41247756195832c2029a 100755 --- a/csst_mci_sim/mci_data/mci_all_9K.config +++ b/csst_mci_sim/mci_data/mci_all_9K.config @@ -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 ############################################## ##############################################