Skip to content
GitLab
Commit 4617177b authored by zjyan's avatar zjyan
Browse files

update pep8 ifs_sim

parent 6326829d
Pipeline #3977 failed with stage
in 0 seconds
Hide whitespace changes
Inline Side-by-side
.DS_Store 0 → 100644
View file @ 4617177b
File added
csst_ifs_sim/csst_ifs_sim.py
View file @ 4617177b
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Mon Mar 7 14:53:41 2022
Created on Thu Apr 11 15:18:57 2024
@author: yan
"""
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
# import sys
import os
from scipy.interpolate import interp1d
import astropy.coordinates as coord
import ctypes
from csst_ifs_sim.support import IFSinstrumentModel
from csst_ifs_sim.support import cosmicrays
from csst_ifs_sim.support import logger as lg
from csst_ifs_sim.CTI import CTI
# from optparse import OptionParser
import configparser as ConfigParser
import cmath
from scipy import ndimage
from scipy.signal import fftconvolve
from astropy.io import fits
import scipy.io as sio
import numpy as np
from scipy import interpolate
from astropy.coordinates import SkyCoord, EarthLocation
from astropy import units as u
from astropy.time import Time
from astropy.coordinates import get_sun
from datetime import datetime, timedelta
import pandas as pd
from scipy.integrate import simps
import galsim
# sys.path.append('../')
from joblib import Parallel, delayed
from astropy.utils.iers import conf
conf.auto_max_age = None
"""
"""
"""
Contact Information: zhaojunyan@shao.ac.cn
The CSST IFS Image Simulator
=============================================
......@@ -40,11 +77,14 @@ a star or an extended source (i.e. a galaxy).
(first step), and finally overlay onto the detector according to its position.
#. Apply calibration unit flux to mimic flat field exposures [optional].
#. Apply a multiplicative flat-field map to emulate pixel-to-pixel non-uniformity [optional].
#. Apply a multiplicative flat-field map to emulate pixel-to-pixel
non-uniformity [optional].
#. Add a charge injection line (horizontal and/or vertical) [optional].
#. Add cosmic ray tracks onto the CCD with random positions but known distribution [optional].
#. Add cosmic ray tracks onto the CCD with random positions but known
distribution [optional].
#. Apply detector charge bleeding in column direction [optional].
#. Add constant dark current and background light from Zodiacal light [optional].
#. Add constant dark current and background light from Zodiacal
light [optional].
#. Include spatially uniform scattered light to the pixel grid [optional].
#. Add photon (Poisson) noise [optional]
#. Add cosmetic defects from an input file [optional].
......@@ -53,135 +93,52 @@ a star or an extended source (i.e. a galaxy).
#. Apply CCD273 non-linearity model to the pixel data [optional].
#. Add readout noise selected from a Gaussian distribution [optional].
#. Convert from electrons to ADUs using a given gain factor.
#. Add a given bias level and discretise the counts (the output is going to be in 16bit unsigned integers).
#. Add a given bias level and discretise the counts (the output is going to be
in 16bit unsigned integers).
#. Finally the simulated image is converted to a FITS file, a WCS is assigned
and the output is saved to the current working directory.
.. Warning:: The code is still work in progress and new features are being added.
The code has been tested, but nevertheless bugs may be lurking in corners, so
please report any weird or inconsistent simulations to the author.
Dependencies
------------
This script depends on the following packages:
:requires: PyFITS (tested with 3.0.6 and 3.3)
:requires: NumPy (tested with 1.6.1, 1.7.1, 1.8.0, 1.9.1, and 1.9.2)
:requires: numexpr (tested with 2.0.1 and 2.3.1)
:requires: SciPy (tested with 0.10.1, 0.12, 0.13, 0.14, and 0.15.1)
:requires: vissim-python package
.. Note:: This class is Python 3 compatible.
.. Note:: CUDA acceleration requires an NVIDIA GPU that supports CUDA and PyFFT and PyCUDA packages.
Note that the CUDA acceleration does not speed up the computations unless oversampled PSF
is being used. If > 2GB of memory is available on the GPU, speed up up to a factor of 50 is
possible.
Testing
-------
the input confiefile as follows:
config_file=['IFS_test_210412.config' ]
you can modify the configfile to change the simulaiton parameters.
please run the as follows::
python IFS_sim_Pro.py
This will run the unittest and compare the result to a previously calculated image.
Please inspect the standard output for results.
Running the test will produce an image representing IFS lower left (0th) quadrant of the CCD (x, y) = (0, 0). Because
noise and cosmic rays are randomised one cannot directly compare the science
outputs but we must rely on the outputs that are free from random effects. The data subdirectory
contains a file named "NoNoiseNoCr_IFS_b.fits, NoNoiseNoCr_IFS_r.fits, NoNoiseNoCr_IFS_i.fits", which is the comparison image without
any noise or cosmic rays.
Warning:: The code is still work in progress and new features are being added.
The code has been tested, but nevertheless bugs may be lurking in corners, so
please report any weird or inconsistent simulations to the author.
:version: 2022.12.15
Note:: This class is Python 3 compatible.
@230613 ,add new lamp hole simulation
Contact Information: zhaojunyan@shao.ac.cn
-------
--23.4.10 --- update the fits header as defined order 0 data format
2023.12.26 add date frame transfer effect and software version from the config file
2024.1.19... add multiple exposure mode function;
2023.06.13 , add new lamp hole simulation
2023.12.26 add date frame transfer effect and software version
2024.1.19 add multiple exposure mode function;
2024.3.20 --- update the fits header as defined order 0 data format
"""
import sys
import os
from scipy.interpolate import interp1d
import astropy.coordinates as coord
import ctypes
from csst_ifs_sim.support import IFSinstrumentModel
from csst_ifs_sim.support import cosmicrays
from csst_ifs_sim.support import logger as lg
from csst_ifs_sim.CTI import CTI
from optparse import OptionParser
import configparser as ConfigParser
import cmath
from scipy import ndimage
from scipy.signal import fftconvolve
from astropy.io import fits
import scipy.io as sio
import numpy as np
from scipy import interpolate
from astropy.coordinates import SkyCoord, EarthLocation
from astropy import units as u
from astropy.time import Time
from astropy.coordinates import get_sun
from datetime import datetime, timedelta
import pandas as pd
from scipy.integrate import simps
from astropy.utils.iers import conf
conf.auto_max_age = None
import galsim
###import julian
sys.path.append('../')
#from joblib import Parallel, delayed
# from astropy.table import Table
###############################################################################
#from sky_bkg import sky_bkg
filterPivotWave = {'nuv': 2875.5, 'u': 3629.6, 'g': 4808.4,
'r': 6178.2, 'i': 7609.0, 'z': 9012.9, 'y': 9627.9}
filterIndex = {'nuv': 0, 'u': 1, 'g': 2, 'r': 3, 'i': 4, 'z': 5, 'y': 6}
# filterCCD = {'nuv':'UV0','u':'UV0','g':'Astro_MB','r':'Astro_MB', 'i':'Basic_NIR', 'z':'Basic_NIR','y':'Basic_NIR'}
# bandRange = {'nuv':[2504.0,3230.0],'u':[3190.0,4039.0],'g':[3989.0,5498.0],'r':[5438.0,6956.0], 'i':[6886.0,8469.0],
# 'z':[8379.0,10855.0],'y':[9217.0, 10900.0], 'GU':[2550, 4000],'GV':[4000, 6200],'GI':[6200,10000]}
# Instrument_dir = '/Users/linlin/Data/csst/straylightsim-master/Instrument/'
# SpecOrder = ['-2','-1','0','1','2']
#
# filterMirrorEff = {'nuv':0.54,'u':0.68,'g':0.8,'r':0.8, 'i':0.8, 'z':0.8,'y':0.8}
def transRaDec2D(ra, dec):
# radec转为竞天程序里的ob, 赤道坐标系下的笛卡尔三维坐标xyz.
"""
ra,dec转为竞天程序里的ob, 赤道坐标系下的笛卡尔三维坐标xyz.
Parameters
----------
ra : float
ra in deg.
dec : float
dec in deg.
Returns
-------
TYPE
DESCRIPTION.
"""
x1 = np.cos(dec / 57.2957795) * np.cos(ra / 57.2957795)
y1 = np.cos(dec / 57.2957795) * np.sin(ra / 57.2957795)
z1 = np.sin(dec / 57.2957795)
......@@ -190,6 +147,22 @@ def transRaDec2D(ra, dec):
def flux2ill(wave, flux):
"""
Parameters
----------
wave : TYPE
DESCRIPTION.
flux : TYPE
DESCRIPTION.
Returns
-------
E : TYPE
DESCRIPTION.
"""
# erg/s/cm^2/A/arcsec^2 to W/m^2
......@@ -208,22 +181,6 @@ def flux2ill(wave, flux):
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)
......@@ -232,14 +189,30 @@ def flux2ill(wave, flux):
delta_lamba = 0.1 # nm
E = np.sum(flux3_interp * delta_lamba)
# pdb.set_trace()
return E
################################################################
def ill2flux(path, E):
"""
Parameters
----------
path : TYPE
DESCRIPTION.
E : TYPE
DESCRIPTION.
Returns
-------
wave0 : TYPE
DESCRIPTION.
spec_scaled : TYPE
DESCRIPTION.
"""
# use template from sky_bkg (background_spec_hst.dat)
filename = path+'IFS_inputdata/refs/background_spec_hst.dat'
......@@ -270,10 +243,33 @@ def ill2flux(path, E):
##############################################################
##########################################################################################
def earth_angle(time_jd, x_sat, y_sat, z_sat, ra_obj, dec_obj):
"""
#
Parameters
----------
time_jd : TYPE
DESCRIPTION.
x_sat : TYPE
DESCRIPTION.
y_sat : TYPE
DESCRIPTION.
z_sat : TYPE
DESCRIPTION.
ra_obj : TYPE
DESCRIPTION.
dec_obj : TYPE
DESCRIPTION.
Returns
-------
angle : TYPE
DESCRIPTION.
"""
ra_sat = np.arctan2(y_sat, x_sat) / np.pi * 180
dec_sat = np.arctan2(z_sat, np.sqrt(x_sat**2+y_sat**2)) / np.pi * 180
......@@ -299,8 +295,7 @@ def earth_angle(time_jd, x_sat, y_sat, z_sat, ra_obj, dec_obj):
return angle
# 3
#####################################################################
class StrayLight(object):
......@@ -344,8 +339,22 @@ class StrayLight(object):
self.slcdll.Init(str.encode(self.deFn), str.encode(
self.PSTFn), str.encode(self.RFn), str.encode(self.ZolFn))
############################################################################
def caculateStarLightFilter(self, filter='i'):
"""
Parameters
----------
filter : TYPE, optional
DESCRIPTION. The default is 'i'.
Returns
-------
TYPE
DESCRIPTION.
"""
sat = (ctypes.c_double*3)()
sat[:] = self.sat
......@@ -367,8 +376,22 @@ class StrayLight(object):
band_star_e2 = star_e2[:][filterIndex[filter]]
return max(band_star_e1, band_star_e2)
###############################################################################
def caculateEarthShineFilter(self, filter='i'):
"""
#
Parameters
----------
filter : TYPE, optional
DESCRIPTION. The default is 'i'.
Returns
-------
TYPE
DESCRIPTION.
"""
sat = (ctypes.c_double*3)()
sat[:] = self.sat
ob = (ctypes.c_double*3)()
......@@ -379,8 +402,8 @@ class StrayLight(object):
self.slcdll.ComposeY(ob, py1, py2)
earth_e1 = (ctypes.c_double*7)()
self.slcdll.EarthShine(self.jtime, sat, ob, py1,
earth_e1) # e[7]代表7个波段的照度
self.slcdll.EarthShine(self.jtime, sat, ob, py1, earth_e1)
# e[7]代表7个波段的照度
earth_e2 = (ctypes.c_double*7)()
self.slcdll.EarthShine(self.jtime, sat, ob, py2, earth_e2)
......@@ -395,23 +418,65 @@ class StrayLight(object):
def str2time(strTime):
"""
Parameters
----------
strTime : TYPE
DESCRIPTION.
Returns
-------
TYPE
DESCRIPTION.
"""
if len(strTime) > 20: # 暂时未用到
msec = int(float('0.'+strTime[20:])*1000000) # 微秒
str2 = strTime[0:19]+' '+str(msec)
return datetime.strptime(str2, '%Y %m %d %H %M %S %f')
# datetime类转mjd
##########################################################################
def time2mjd(dateT):
"""
Parameters
----------
dateT : TYPE
DESCRIPTION.
Returns
-------
TYPE
DESCRIPTION.
"""
t0 = datetime(1858, 11, 17, 0, 0, 0, 0) # 简化儒略日起始日
mjd = (dateT-t0).days
mjd_s = dateT.hour*3600.0+dateT.minute*60.0 + \
dateT.second+dateT.microsecond/1000000.0
return mjd+mjd_s/86400.0
###########################################################################
def time2jd(dateT):
"""
Parameters
----------
dateT : TYPE
DESCRIPTION.
Returns
-------
TYPE
DESCRIPTION.
"""
t0 = datetime(1858, 11, 17, 0, 0, 0, 0) # 简化儒略日起始日
mjd = (dateT-t0).days
mjd_s = dateT.hour*3600.0+dateT.minute*60.0 + \
......@@ -419,29 +484,99 @@ def time2jd(dateT):
return mjd+mjd_s/86400.0++2400000.5
# mjd转datetime类
#############################################################################
def mjd2time(mjd):
"""
Parameters
----------
mjd : TYPE
DESCRIPTION.
Returns
-------
TYPE
DESCRIPTION.
"""
t0 = datetime(1858, 11, 17, 0, 0, 0, 0) # 简化儒略日起始日
return t0+timedelta(days=mjd)
#############################################################################
def jd2time(jd):
"""
Parameters
----------
jd : TYPE
DESCRIPTION.
Returns
-------
TYPE
DESCRIPTION.
"""
mjd = jd2mjd(jd)
return mjd2time(mjd)
# mjd和jd互转
############################################################################
def mjd2jd(mjd):
return mjd+2400000.5
"""
Parameters
----------
mjd : TYPE
DESCRIPTION.
Returns
-------
TYPE
DESCRIPTION.
"""
return mjd+2400000.5
###########################################################################
def jd2mjd(jd):
return jd-2400000.5
"""
Parameters
----------
jd : TYPE
DESCRIPTION.
Returns
-------
TYPE
DESCRIPTION.
"""
return jd-2400000.5
##########################################################################
def dt2hmd(dt):
"""
Parameters
----------
dt : TYPE
DESCRIPTION.
Returns
-------
TYPE
DESCRIPTION.
"""
## dt is datetime
hour = dt.hour
minute = dt.minute
......@@ -464,9 +599,23 @@ def dt2hmd(dt):
return str_h+':'+str_m+':'+str_d
###############################################################################
def float2char(a, z):
"""
def float2char(a, z):
Parameters
----------
a : TYPE
DESCRIPTION.
z : TYPE
DESCRIPTION.
Returns
-------
TYPE
DESCRIPTION.
"""
# a is input float value
# transfer float a to chars and save z bis
b = str(a)
......@@ -481,6 +630,22 @@ def float2char(a, z):
###############################################################################
def deg2HMS(ra0, dec0):
"""
Parameters
----------
ra0 : TYPE
DESCRIPTION.
dec0 : TYPE
DESCRIPTION.
Returns
-------
TYPE
DESCRIPTION.
"""
'''convert deg to ra's HMS and dec's DMS'''
c = SkyCoord(ra=ra0*u.degree, dec=dec0*u.degree)
ss = c.to_string('hmsdms')
......@@ -517,9 +682,35 @@ def deg2HMS(ra0, dec0):
############################################################################
###########################################################
def beta_angle(x_sat, y_sat, z_sat, vx_sat, vy_sat, vz_sat, ra_obj, dec_obj):
"""
Parameters
----------
x_sat : TYPE
DESCRIPTION.
y_sat : TYPE
DESCRIPTION.
z_sat : TYPE
DESCRIPTION.
vx_sat : TYPE
DESCRIPTION.
vy_sat : TYPE
DESCRIPTION.
vz_sat : TYPE
DESCRIPTION.
ra_obj : TYPE
DESCRIPTION.
dec_obj : TYPE
DESCRIPTION.
Returns
-------
angle : TYPE
DESCRIPTION.
"""
# get the vector for next motion
x_sat2 = x_sat + vx_sat
......@@ -549,10 +740,28 @@ def beta_angle(x_sat, y_sat, z_sat, vx_sat, vy_sat, vz_sat, ra_obj, dec_obj):
return angle
###############################################################################
##########################################################
def LSR_velocity(ra, dec, velocity, Obstime):
"""
Parameters
----------
ra : TYPE
DESCRIPTION.
dec : TYPE
DESCRIPTION.
velocity : TYPE
DESCRIPTION.
Obstime : TYPE
DESCRIPTION.
Returns
-------
TYPE
DESCRIPTION.
"""
# '''
# 输入参数为ra，dec 原始视向速度in km/s,观测时间
# local为观测地位置，c为利用输入坐标建立的坐标系对象。
......@@ -580,6 +789,30 @@ def LSR_velocity(ra, dec, velocity, Obstime):
def rotation_yan(x0, y0, x1, y1, angle):
"""
Parameters
----------
x0 : TYPE
DESCRIPTION.
y0 : TYPE
DESCRIPTION.
x1 : TYPE
DESCRIPTION.
y1 : TYPE
DESCRIPTION.
angle : TYPE
DESCRIPTION.
Returns
-------
x2 : TYPE
DESCRIPTION.
y2 : TYPE
DESCRIPTION.
"""
# % point A (x0,y0)
# % point B(x1,y1)
# % roate angle ,point B rotate with point A
......@@ -591,6 +824,20 @@ def rotation_yan(x0, y0, x1, y1, angle):
def centroid(data):
"""
Parameters
----------
data : TYPE
DESCRIPTION.
Returns
-------
TYPE
DESCRIPTION.
"""
h, w = np.shape(data)
x = np.arange(0, w)
y = np.arange(0, h)
......@@ -600,10 +847,24 @@ def centroid(data):
cy = (np.dot(np.dot(x, data), y1)) / (np.dot(np.dot(x1, data), y1))
return float(cx), float(cy)
# 33
###############################################################################
def centroidN(data):
"""
def centroidN(data):
Parameters
----------
data : TYPE
DESCRIPTION.
Returns
-------
cx : TYPE
DESCRIPTION.
cy : TYPE
DESCRIPTION.
"""
'''
calculate the centroid of the input two-dimentional image data
......@@ -624,8 +885,23 @@ def centroidN(data):
####################################################################
def SpecCube2photon(data, wave):
"""
Parameters
----------
data : TYPE
DESCRIPTION.
wave : TYPE
DESCRIPTION.
Returns
-------
Nphoton : TYPE
DESCRIPTION.
"""
# calcutle photons from original spec cube data,
# data: input data, two-dimentional, unit : 1e-17 erg/s/A/cm^2
# wave: the relative wavefront, unit in nm;
......@@ -640,11 +916,26 @@ def SpecCube2photon(data, wave):
return Nphoton
######################################################################
######################################################################
def mag2flux(starmag, lam):
"""
####################################################################################
Parameters
----------
starmag : TYPE
DESCRIPTION.
lam : TYPE
DESCRIPTION.
Returns
-------
fluxlam : TYPE
DESCRIPTION.
photonu : TYPE
DESCRIPTION.
def mag2flux(starmag, lam):
"""
# starmag: input AB mag
# lambda: wavelength in nm
......@@ -662,8 +953,23 @@ def mag2flux(starmag, lam):
###############################################################################
def flux2photons(flux, lam):
"""
Parameters
----------
flux : TYPE
DESCRIPTION.
lam : TYPE
DESCRIPTION.
Returns
-------
photonu : TYPE
DESCRIPTION.
"""
# flux: # erg/s/cm2/A/
# lam: # wavelength in nm
......@@ -674,17 +980,48 @@ def flux2photons(flux, lam):
return photonu
###############################################################################
##################################################
def fftrange(n, dtype=float):
"""
def fftrange(n, dtype=float):
Parameters
----------
n : TYPE
DESCRIPTION.
dtype : TYPE, optional
DESCRIPTION. The default is float.
Returns
-------
TYPE
DESCRIPTION.
"""
"""FFT-aligned coordinate grid for n samples."""
return np.arange(-n//2, -n//2+n, dtype=dtype)
###############################################################################
def dft2(ary, Q, samples, shift=None):
"""
def dft2(ary, Q, samples, shift=None):
Parameters
----------
ary : TYPE
DESCRIPTION.
Q : TYPE
DESCRIPTION.
samples : TYPE
DESCRIPTION.
shift : TYPE, optional
DESCRIPTION. The default is None.
Returns
-------
out : TYPE
DESCRIPTION.
"""
"""Compute the two dimensional Discrete Fourier Transform of a matrix.
Parameters
......@@ -723,59 +1060,35 @@ def dft2(ary, Q, samples, shift=None):
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 anySampledPSFnew(wavefront, pupil, Q, sizeout):
"""
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.
wavefront : TYPE
DESCRIPTION.
pupil : TYPE
DESCRIPTION.
Q : TYPE
DESCRIPTION.
sizeout : TYPE
DESCRIPTION.
Returns
-------
`numpy.ndarray`
2D array containing the shifted transform.
Equivalent to ifftshift(ifft2(fftshift(ary))) modulo output
sampling/grid differences
psf : TYPE
DESCRIPTION.
"""
# 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_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 anySampledPSFnew(wavefront, pupil, Q, sizeout):
'''
'''
written on 2020.12.04 by Yan @Shao
% wavefront sampled in the united circle;
% pupil function samped as wavefront;
......@@ -804,7 +1117,7 @@ def anySampledPSFnew(wavefront, pupil, Q, sizeout):
cx, cy = centroid(psfout)
Nt = sizeout
# psf=ndimage.shift(psfout,[Nt/2-cy, Nt/2-cx],order=1, mode='nearest' ) ## for my fft method
# for convolve method
psf = ndimage.shift(
psfout, [Nt/2-cy-1, Nt/2-cx-1], order=1, mode='nearest')
......@@ -812,79 +1125,60 @@ def anySampledPSFnew(wavefront, pupil, Q, sizeout):
return psf
###########################################################################################
##############################################################################
##############################################################################
def zero_pad(image, shape, position='corner'):
def get_dx_dy_blue(wave):
"""
Extends image to a certain size with zeros
Parameters
----------
image: real 2d `np.ndarray`
Input image
shape: tuple of int
Desired output shape of the image
position : str, optional
The position of the input image in the output one:
* 'corner'
top-left corner (default)
* 'center'
centered
wave : TYPE
DESCRIPTION.
Returns
-------
padded_img: real `np.ndarray`
The zero-padded image
"""
shape = np.asarray(shape, dtype=int)
imshape = np.asarray(image.shape, dtype=int)
if np.alltrue(imshape == shape):
return image
if np.any(shape <= 0):
raise ValueError("ZERO_PAD: null or negative shape given")
dshape = shape - imshape
if np.any(dshape < 0):
raise ValueError("ZERO_PAD: target size smaller than source one")
pad_img = np.zeros(shape, dtype=image.dtype)
idx, idy = np.indices(imshape)
if position == 'center':
if np.any(dshape % 2 != 0):
raise ValueError("ZERO_PAD: source and target shapes "
"have different parity.")
offx, offy = dshape // 2
else:
offx, offy = (0, 0)
pad_img[idx + offx, idy + offy] = image
return pad_img
################################################################################
##############################################################################
dx : TYPE
DESCRIPTION.
dy : TYPE
DESCRIPTION.
def get_dx_dy_blue(wave):
"""
# wave is the wavelength in nm;
# dx is in dispersion direction, dy is in vertical direction;
dydl = np.array([-423.256, 0.001, 0.00075,
0.0000078, -0.0000000000007, 0.0, 0.0]) # 色散方向
dxdl = 0.2*np.array([-9.1519, -1.00000000e-06, 3.50000000e-08, -5.00000000e-09, -
1.70000000e-11, 4.00949787e-12, -6.16873452e-15]) # 垂直方向
0.0000078, -0.0000000000007, 0.0, 0.0])
# 色散方向
dxdl = 0.2*np.array([-9.1519, -1.00000000e-06, 3.50000000e-08, -5.00000000e-09,
-1.70000000e-11, 4.00949787e-12, -6.16873452e-15])
# 垂直方向
dx = 0.0
dy = 0.0
for i in range(len(dxdl)):
dx = dx+dxdl[i]*wave**(i)
dy = dy+dydl[i]*wave**(i)
return dx, dy
##############################################################################
def get_dx_dy_red(wave):
"""
Parameters
----------
wave : TYPE
DESCRIPTION.
Returns
-------
dx : TYPE
DESCRIPTION.
dy : TYPE
DESCRIPTION.
"""
# wave is the wavelength in nm;
dydl = np.array([-640.0239901372472, 0.0018, 0.00048,
0.0000028, -0.0000000000007, 0.0, 0.0]) # 色散方向
......@@ -900,8 +1194,25 @@ def get_dx_dy_red(wave):
##############################################################################################
def getSpectrum(Spectrum0, lam, sigma):
"""
Parameters
----------
Spectrum0 : TYPE
DESCRIPTION.
lam : TYPE
DESCRIPTION.
sigma : TYPE
DESCRIPTION.
Returns
-------
SpmOut : TYPE
DESCRIPTION.
"""
# %获得输入波长lambda的光谱强度
wave = Spectrum0[:, 0] # %原始光谱数据给定的波长；
......@@ -918,39 +1229,8 @@ def getSpectrum(Spectrum0, lam, sigma):
return SpmOut
#####################################################################################
####################################################################################################################
def processArgs(printHelp=False):
"""
Processes command line arguments.
"""
parser = OptionParser()
parser.add_option('-c', '--configfile', dest='configfile',
help="Name of the configuration file", metavar="string")
parser.add_option('-s', '--section', dest='section',
help="Name of the section of the config file [SCIENCE]", metavar="string")
parser.add_option('-q', '--quadrant', dest='quadrant', help='CCD quadrant to simulate [0, 1, 2, 3]',
metavar='int')
parser.add_option('-x', '--xCCD', dest='xCCD', help='CCD number in X-direction within the FPA matrix',
metavar='int')
parser.add_option('-y', '--yCCD', dest='yCCD', help='CCD number in Y-direction within the FPA matrix',
metavar='int')
parser.add_option('-d', '--debug', dest='debug', action='store_true',
help='Debugging mode on')
parser.add_option('-t', '--test', dest='test', action='store_true',
help='Run unittest')
parser.add_option('-f', '--fixed', dest='fixed', help='Use a fixed seed for the random number generators',
metavar='int')
if printHelp:
parser.print_help()
else:
return parser.parse_args()
##############################################################################################
##############################################################################
##############################################################################
class IFSsimulator():
"""
......@@ -978,27 +1258,9 @@ class IFSsimulator():
####################################
self.configfile = configfile
self.section = 'TEST' # simulation section;
# if opts.section is None:
# ####self.section = 'DEFAULT'
# self.section = 'TEST' # simulation section;
# else:
# self.section = opts.section
self.section = 'TEST' # simulation section;
# if opts.debug is None:
# self.debug = False
# else:
# self.debug = opts.debug
# try:
# self.random = opts.testing
# except:
# self.random = False
# load instrument model, these values are also stored in the FITS header
# load instrument model
self.information = IFSinstrumentModel.IFSinformation()
# update settings with defaults
......@@ -1026,20 +1288,29 @@ class IFSsimulator():
ra=0.0,
dec=0.0,
injection=80000.0,
coveringfraction=0.5, # CR: 3.0% is for 565s exposure
###################################################
# cosmicraylengths=self.information['dir_path']+'IFS_inputdata/cdf_cr_length.dat',
# cosmicraydistance=self.information['dir_path']+'IFS_inputdata/cdf_cr_total.dat',
# parallelTrapfile=self.information['dir_path']+'IFS_inputdata/cdm_euclid_parallel.dat',
# serialTrapfile=self.information['dir_path']+'IFS_inputdata/cdm_euclid_serial.dat',
# cosmeticsFile='../IFS_inputdata/cosmetics.dat',
coveringfraction=0.5,
# CR: 3.0% is for 565s exposure
#########################################
mode='same'))
#####################################
# 3
def readConfigs(self, simnumber):
"""
Parameters
----------
simnumber : TYPE
DESCRIPTION.
Returns
-------
None.
"""
"""
Reads the config file information using configParser and sets up a logger.
"""
self.config = ConfigParser.RawConfigParser()
......@@ -1048,9 +1319,9 @@ class IFSsimulator():
# self.log.info(self.information)
###############################################################################
def processConfigs(self):
"""
Processes configuration information and save the information to a dictionary self.information.
......@@ -1083,15 +1354,6 @@ class IFSsimulator():
# force gain to be float
self.information['e_adu'] = float(self.information['e_adu'])
# #name of the output file, include quadrants and CCDs
###################################################################################################
#self.information['output'] = self.config.get(self.section, 'output')
##################################################################################################
# #booleans to control the flow
# self.chargeInjectionx = self.config.getboolean(self.section, 'chargeInjectionx')
# self.chargeInjectiony = self.config.getboolean(self.section, 'chargeInjectiony')
self.cosmicRays = self.config.getboolean(self.section, 'cosmicRays')
self.darknoise = self.config.getboolean(self.section, 'darknoise')
self.cosmetics = self.config.getboolean(self.section, 'cosmetics')
......@@ -1126,7 +1388,7 @@ class IFSsimulator():
except:
self.intscale = True
##############################################################################
######################################################################
# 3
......@@ -1146,7 +1408,7 @@ class IFSsimulator():
#####################################################################
#########################################################################################
############################################################################
def _createEmpty(self):
"""
......@@ -1163,12 +1425,32 @@ class IFSsimulator():
self.pixel = 0.1 # arcsec, pixel scale size;
#########################################################################################################
##########################################################
##############################################################################
##############################################################################
def zodiacal(self, ra, dec, time):
"""
For given RA, DEC and TIME, return the interpolated zodical spectrum in Leinert-1998.
Parameters
----------
ra : TYPE
DESCRIPTION.
dec : TYPE
DESCRIPTION.
time : TYPE
DESCRIPTION.
Returns
-------
wave_A : TYPE
DESCRIPTION.
spec_erg2 : TYPE
DESCRIPTION.
"""
"""
For given RA, DEC and TIME, return the interpolated zodical
spectrum in Leinert-1998.
:param ra: RA in unit of degree, ICRS frame
:param dec: DEC in unit of degree, ICRS frame
......@@ -1183,7 +1465,7 @@ class IFSsimulator():
# get solar position
dt = datetime.fromisoformat(time)
jd = time2jd(dt)
##jd = julian.to_jd(dt, fmt='jd')
## jd = julian.to_jd(dt, fmt='jd')
t = Time(jd, format='jd', scale='utc')
......@@ -1236,10 +1518,26 @@ class IFSsimulator():
# self.zodiacal_flux=spec_erg2
return wave_A, spec_erg2
###################################################################################
##########################################################################
def smoothingWithChargeDiffusion(self, image, sigma=(0.32, 0.32)):
"""
Parameters
----------
image : TYPE
DESCRIPTION.
sigma : TYPE, optional
DESCRIPTION. The default is (0.32, 0.32).
Returns
-------
TYPE
DESCRIPTION.
"""
"""
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.
......@@ -1258,21 +1556,28 @@ class IFSsimulator():
:rtype: ndarray
"""
return ndimage.filters.gaussian_filter(image, sigma)
#################################################################################
###############################################################################
def readCosmicRayInformation(self):
"""
Returns
-------
None.
def readCosmicRayInformation(self):
"""
"""
Reads in the cosmic ray track information from two input files.
Stores the information to a dictionary called cr.
"""
# self.log.info('Reading in cosmic ray information from %s and %s' % (self.information['cosmicraylengths'],
# self.information['cosmicraydistance']))
self.log.info('Reading in cosmic ray information from %s and %s' % (self.information['cosmicraylengths'],
self.information['cosmicraydistance']))
print(self.information)
crLengths = np.loadtxt(
self.information['dir_path']+self.information['cosmicraylengths'])
self.information['dir_path']+self.information['cosmicraylengths'])
crDists = np.loadtxt(
self.information['dir_path']+self.information['cosmicraydistance'])
......@@ -1280,92 +1585,106 @@ class IFSsimulator():
cr_v=crDists[:, 0], cr_cde=crDists[:, 1], cr_cden=np.shape(crDists)[0])
###############################################################################
def _writeFITSfile(self, image, filename):
"""
:param image: image array to save
:type image: ndarray
:param filename: name of the output file, e.g. file.fits
:type filename: str
# def _writeFITSfile(self, image, filename):
# """
# :param image: image array to save
# :type image: ndarray
# :param filename: name of the output file, e.g. file.fits
# :type filename: str
:return: None
"""
if os.path.isfile(filename):
os.remove(filename)
# :return: None
# """
# if os.path.isfile(filename):
# os.remove(filename)
# create a new FITS file, using HDUList instance
ofd = fits.HDUList(fits.PrimaryHDU())
# # create a new FITS file, using HDUList instance
# ofd = fits.HDUList(fits.PrimaryHDU())
# new image HDU
hdu = fits.ImageHDU(data=image)
# # new image HDU
# hdu = fits.ImageHDU(data=image)
# update and verify the header
hdu.header.add_history('Created by IFSsim at %s' %
datetime.isoformat(datetime.utcnow()))
hdu.verify('fix')
# # update and verify the header
# hdu.header.add_history('Created by IFSsim at %s' %
# datetime.isoformat(datetime.utcnow()))
# hdu.verify('fix')
ofd.append(hdu)
# ofd.append(hdu)
# write the actual file
ofd.writeto(filename)
# # write the actual file
# ofd.writeto(filename)
###############################################################################
def writeFITSfile(self, data, filename, unsigned16bit=False):
"""
Writes out a simple FITS file.
:param data: data to be written
:type data: ndarray
:param filename: name of the output file
:type filename: str
:param unsigned16bit: whether to scale the data using bzero=32768
:type unsigned16bit: bool
# def writeFITSfile(self, data, filename, unsigned16bit=False):
# """
# Writes out a simple FITS file.
# :param data: data to be written
# :type data: ndarray
# :param filename: name of the output file
# :type filename: str
# :param unsigned16bit: whether to scale the data using bzero=32768
# :type unsigned16bit: bool
# :return: None
# """
# if os.path.isfile(filename):
# os.remove(filename)
# # create a new FITS file, using HDUList instance
# ofd = fits.HDUList(fits.PrimaryHDU())
# # new image HDU
# hdu = fits.ImageHDU(data=data)
# # add input keywords to the header
# for key, value in self.information.items():
# # truncate long keys
# if len(key) > 8:
# key = key[:7]
# try:
# hdu.header.set(key.upper(), value)
# except:
# try:
# hdu.header.set(key.upper(), str(value))
# except:
# pass
# # write booleans
# for key, value in self.booleans.items():
# # truncate long keys
# if len(key) > 8:
# key = key[:7]
# hdu.header.set(key.upper(), str(value), 'Boolean Flags')
# # update and verify the header
# hdu.header.add_history(
# 'This is an itermediate data product no the final output!')
# hdu.verify('fix')
# ofd.append(hdu)
# # write the actual file
# ofd.writeto(filename)
:return: None
###############################################################################
def configure(self, simnumber):
"""
if os.path.isfile(filename):
os.remove(filename)
# create a new FITS file, using HDUList instance
ofd = fits.HDUList(fits.PrimaryHDU())
# new image HDU
hdu = fits.ImageHDU(data=data)
# add input keywords to the header
for key, value in self.information.items():
# truncate long keys
if len(key) > 8:
key = key[:7]
try:
hdu.header.set(key.upper(), value)
except:
try:
hdu.header.set(key.upper(), str(value))
except:
pass
# write booleans
for key, value in self.booleans.items():
# truncate long keys
if len(key) > 8:
key = key[:7]
hdu.header.set(key.upper(), str(value), 'Boolean Flags')
# update and verify the header
hdu.header.add_history(
'This is an itermediate data product no the final output!')
hdu.verify('fix')
Parameters
----------
simnumber : TYPE
DESCRIPTION.
ofd.append(hdu)
# write the actual file
ofd.writeto(filename)
Returns
-------
None.
###############################################################################
def configure(self, simnumber):
"""
Configures the simulator with input information and creates and empty array to which the final image will
"""
Configures the simulator with input information and creates and
empty array to which the final image will
be build on.
"""
self.readConfigs(simnumber)
......@@ -1391,8 +1710,9 @@ class IFSsimulator():
else:
ss = '_'
if self.information['dir_path']=='/nfsdata/share/simulation-unittest/ifs_sim/':
self.result_path = self.information['dir_path']+'ifs_sim_result/'+self.source+ss+result_day
if self.information['dir_path'] == '/nfsdata/share/simulation-unittest/ifs_sim/':
self.result_path = self.information['dir_path'] + \
'ifs_sim_result/'+self.source+ss+result_day
else:
home_path = os.environ['HOME']
......@@ -1402,9 +1722,6 @@ class IFSsimulator():
self.result_path = '../IFS_simData_'+self.source+ss+result_day
else:
self.result_path = '/data/ifspip/CCD_ima/IFS_simData_'+self.source+ss+result_day
if os.path.isdir(self.result_path) == False:
os.mkdir(self.result_path)
......@@ -1526,7 +1843,7 @@ class IFSsimulator():
self.slice_blue = slice_blue
self.slice_red = slice_red
###############################################################################
#######################################################################
maskSlice = dict()
maskSlit = dict()
sizeout = 100
......@@ -1541,19 +1858,37 @@ class IFSsimulator():
self.maskSlice = maskSlice
self.maskSlit = maskSlit
################################################################################################
######################################################################
#################################################################################################################
##############################################################################
def generateflat(self, ave=1.0, sigma=0.01):
"""
Parameters
----------
ave : TYPE, optional
DESCRIPTION. The default is 1.0.
sigma : TYPE, optional
DESCRIPTION. The default is 0.01.
Returns
-------
TYPE
DESCRIPTION.
TYPE
DESCRIPTION.
"""
"""
Creates a flat field image with given properties.
:return: flat field image
:rtype: ndarray
"""
self.log.info('Generating a flat field...')
self.log.info(
'The flat field has mean value of 1 and a given fluctuations, usually either 1 or 2 percent defined by sigma= %d...' % sigma)
self.log.info('The flat field has mean value of 1 and a given fluctuations,
usually either 1 or 2 percent defined by sigma= %d...' % sigma)
np.random.seed(5*self.simnumber)
self.flat_b = np.random.normal(loc=ave, scale=sigma, size=(2048, 4096))
......@@ -1589,10 +1924,18 @@ class IFSsimulator():
return self.flat_b, self.flat_r
#########################################################################################################
##########################################################################
def addLampFlux(self):
"""
Returns
-------
None.
"""
"""
Include flux from the calibration source.
"""
......@@ -1604,6 +1947,22 @@ class IFSsimulator():
#############################################################################
def MakeFlatMatrix(self, img, seed):
"""
Parameters
----------
img : TYPE
DESCRIPTION.
seed : TYPE
DESCRIPTION.
Returns
-------
FlatMat : TYPE
DESCRIPTION.
"""
####
ysize, xsize = img.shape
np.random.seed(seed)
......@@ -1625,6 +1984,14 @@ class IFSsimulator():
def getFrameTransferImg(self):
"""
Returns
-------
None.
"""
"""
get frame transfer image from original image of self.image_b and self.image_r
"""
self.frame_b_4 = np.zeros((1024+320, 2048), dtype=float)
......@@ -1670,11 +2037,24 @@ class IFSsimulator():
###############################################################################
def applyflatfield(self, simnumber):
"""
Parameters
----------
simnumber : TYPE
DESCRIPTION.
Returns
-------
None.
"""
"""
Applies multiplicative flat field to emulate pixel-to-pixel non-uniformity.
Because the pixel-to-pixel non-uniformity effect (i.e. multiplicative) flat fielding takes place
before CTI and other effects, the flat field file must be the same size as the pixels that see
the sky.
the sky.
"""
# apply Flat field to image
###
......@@ -1690,31 +2070,30 @@ class IFSsimulator():
return
#################################################################################
###############################################################################
def addChargeInjection(self):
"""
Add either horizontal or vertical charge injection line to the image.
###############################################################################
def addCosmicRays(self, idk):
"""
if self.chargeInjectionx:
#self.image[1500:1511, :] = self.information['injection']
self.image_b[1500:1511, :] = self.information['injection']
self.image_r[1500:1511, :] = self.information['injection']
self.log.info('Adding vertical charge injection line')
Parameters
----------
idk : TYPE
DESCRIPTION.
if self.chargeInjectiony:
#self.image[:, 1500:1511] = self.information['injection']
self.image_b[:, 1500:1511] = self.information['injection']
self.image_r[:, 1500:1511] = self.information['injection']
self.log.info('Adding horizontal charge injection line')
Returns
-------
None.
###############################################################################
def addCosmicRays(self, idk):
"""
Add cosmic rays to the arrays based on a power-law intensity distribution for tracks.
Cosmic ray properties (such as location and angle) are chosen from random Uniform distribution.
For details, see the documentation for the cosmicrays class in the support package.
"""
Add cosmic rays to the arrays based on a power-law intensity
distribution for tracks.
Cosmic ray properties (such as location and angle) are chosen from
random Uniform distribution.
For details, see the documentation for the cosmicrays class in the
support package.
"""
self.readCosmicRayInformation()
# to scale the number of cosmics with exposure time
......@@ -1761,6 +2140,14 @@ class IFSsimulator():
def addReadoutTrails(self):
"""
Returns
-------
None.
"""
"""
Add readout trails resulting from reading out the shutter open.
Quadrants assumed to be numbered:
......@@ -1791,7 +2178,7 @@ class IFSsimulator():
else:
self.image_b = data
# 3
flux_ratio = self.information['readouttime'] / float(
self.information['redsize']) / self.information['exptime']
# make a copy, this will be updated
......@@ -1815,10 +2202,18 @@ class IFSsimulator():
else:
self.image_r = data
##############################################################################################################
##############################################################################
def applyDarkCurrent(self):
"""
Returns
-------
None.
"""
"""
Apply dark current. Scales the dark with the exposure time.
Additionally saves the image without noise to a FITS file.
......@@ -1859,11 +2254,18 @@ class IFSsimulator():
self.information['dark2_r']
# 3
##############################################################################
def applyCosmicBackground(self):
"""
Returns
-------
None.
"""
"""
Apply dark the cosmic background. Scales the background with the exposure time.
Additionally saves the image without noise to a FITS file.
......@@ -1882,10 +2284,18 @@ class IFSsimulator():
self.imagenoCR_b += bcgr
self.imagenoCR_r += bcgr
##########################################################################################
##########################################################################
def applyScatteredLight(self):
"""
Returns
-------
None.
"""
"""
Adds spatially uniform scattered light to the image.
"""
sl = self.information['exptime'] * self.information['scattered_light']
......@@ -1898,6 +2308,14 @@ class IFSsimulator():
##############################################################################
def applyPoissonNoise(self):
"""
Returns
-------
None.
"""
"""
Add Poisson noise to the image.
"""
......@@ -1923,9 +2341,16 @@ class IFSsimulator():
###################################################################################################################
def applyCosmetics(self):
"""
Returns
-------
None.
"""
"""
Apply cosmetic defects described in the input file.
#Number of hot and dead pixels from MSSL/Euclid/TR/12003 Issue 2 Draft b
......@@ -1952,7 +2377,7 @@ class IFSsimulator():
self.image_b[xc, yc] = val
# cosmetics_b[xc,yc]=val
self.log.info('x=%i, y=%i, value=%f' % (xc, yc, val))
######################################################################################################
#############################################################################
cosmetics = np.loadtxt(
self.information['dir_path']+self.information['cosmeticsfile_r'])
......@@ -1991,6 +2416,14 @@ class IFSsimulator():
def applyRadiationDamage(self):
"""
Returns
-------
None.
"""
"""
Applies CDM03 radiation model to the image being constructed.
.. seealso:: Class :`CDM03`
......@@ -2011,10 +2444,18 @@ class IFSsimulator():
self.image_r = cti.applyRadiationDamage(self.image_r.copy(
).transpose(), iquadrant=self.information['quadrant']).transpose()
self.log.info('Radiation damage added.')
##################################################################################
##############################################################################
def applyNonlinearity(self):
"""
Returns
-------
None.
"""
"""
Applies a CCD273 non-linearity model to the image being constructed.
"""
......@@ -2029,13 +2470,22 @@ class IFSsimulator():
self.image_r.copy())
self.log.info('Non-linearity effects included.')
########################################################################################
######################################################################
def applyReadoutNoise(self):
"""
Returns
-------
None.
"""
"""
Applies readout noise to the image being constructed.
The noise is drawn from a Normal (Gaussian) distribution with average=0.0 and std=readout noise.
The noise is drawn from a Normal (Gaussian) distribution with
average=0.0 and std=readout noise.
"""
self.log.info('readnoise added in blue channel')
# blue zone 1
......@@ -2061,7 +2511,6 @@ class IFSsimulator():
self.image_b[0:1024+overscan, 2418*3:2418*3+2048+prescan+overscan] += np.random.normal(
loc=0.0, scale=self.information['rn2_b'], size=(1344, 2418))
# 3
self.log.info('readnoise added in blue channel')
......@@ -2089,6 +2538,14 @@ class IFSsimulator():
def appFrameTransferEffect(self):
"""
Returns
-------
None.
"""
"""
apply frame transfer effect to the four part images
"""
prescan = 50
......@@ -2136,10 +2593,18 @@ class IFSsimulator():
self.image_r[0:1536+overscan, prescan+3442*3:3442*3 +
3072+prescan] += self.frame_r_2[0:1536+overscan, 0:3072]
# 3
##########################################################################
def electrons2ADU(self):
"""
Returns
-------
None.
"""
"""
Convert from electrons to ADUs using the value read from the configuration file.
"""
self.log.info(
......@@ -2177,6 +2642,14 @@ class IFSsimulator():
def applyBias(self):
"""
Returns
-------
None.
"""
"""
Adds a bias level to the image being constructed.
The value of bias is read from the configure file and stored
......@@ -2199,7 +2672,7 @@ class IFSsimulator():
self.log.info('Bias counts were added to the blue image')
############################################################################
#######################################################################
# red zone 4
self.image_r[0:1856, 0:3442] += self.information['bias4_r']
......@@ -2213,14 +2686,21 @@ class IFSsimulator():
# zone 2
self.image_r[0:1856, 3442*3:3442*4] += self.information['bias2_r']
##########################################################################
#######################################################################
self.log.info('Bias counts were added to the red image')
# 3
###############################################################################
def addPreOverScans(self):
"""
Returns
-------
None.
"""
"""
Add pre- and overscan regions to the self.image. These areas are added only in the serial direction.
Because the 1st and 3rd quadrant are read out in to a different serial direction than the nominal
orientation, in these images the regions are mirrored.
......@@ -2265,13 +2745,23 @@ class IFSsimulator():
self.log.error(
'Cannot include pre- and overscan because of an unknown quadrant!')
####################################################################################
###############################################################################
def applyBleeding_yan(self):
"""
Apply bleeding along the CCD columns if the number of electrons in a pixel exceeds the full-well capacity.
Returns
-------
None.
"""
"""
Apply bleeding along the CCD columns if the number of electrons in
a pixel exceeds the full-well capacity.
Bleeding is modelled in the parallel direction only, because the CCD273s are assumed not to bleed in
Bleeding is modelled in the parallel direction only, because the
CCD273s are assumed not to bleed in
serial direction.
:return: None
......@@ -2355,12 +2845,25 @@ class IFSsimulator():
sum += overload
print('Applying column bleeding to red image finished.......')
############################################################################
###########################################################################
############################################################################
def discretise(self, max=2**16-1):
"""
Parameters
----------
max : TYPE, optional
DESCRIPTION. The default is 2**16-1.
Returns
-------
None.
"""
"""
Converts a floating point image array (self.image) to an integer array with max values
defined by the argument max.
......@@ -2390,6 +2893,14 @@ class IFSsimulator():
##################################################################################################
def applyImageShift(self):
"""
Returns
-------
None.
"""
np.random.seed(9*self.simnumber)
ud = np.random.random() # Choose a random rotation
......@@ -2409,10 +2920,17 @@ class IFSsimulator():
self.information['dec'] = dy*self.information['pixel_size']
# 33
##############################################################################
def applyImageRotate(self):
"""
Returns
-------
None.
"""
np.random.seed(10*self.simnumber)
ud = np.random.random() # Choose a random rotation
angle = 2 * (ud-0.5) * self.information['tel_rotmax']
......@@ -2436,6 +2954,14 @@ class IFSsimulator():
###############################################################################
def CCDreadout(self):
"""
Returns
-------
None.
"""
imgb = self.image_b.copy()
prescan = int(self.information['prescan'])
......@@ -2449,7 +2975,7 @@ class IFSsimulator():
y1 = 0+prescan
y2 = y1+2048
temp[x1:x2, y1:y2] = imgb[0:1024, 0:2048]
########## zone 3, OSG , left to right #################
# zone 3, OSG , left to right #################
# np.fliplr(b2) ## left to right
# np.flipud(b3) ## down to up
x1 = 0
......@@ -2458,7 +2984,7 @@ class IFSsimulator():
y1 = 2418+prescan
y2 = y1+2048
temp[x1:x2, y1:y2] = np.fliplr(imgb[0:1024, 2048:4096])
########## zone 1, OSE,down to up ###################
# zone 1, OSE,down to up ###################
x1 = 0
x2 = x1+1024
......@@ -2475,7 +3001,7 @@ class IFSsimulator():
self.image_b = temp
##############################################################################
#######################################################################
imgr = self.image_r.copy()
temp = np.zeros((1856, 13768))
......@@ -2517,6 +3043,19 @@ class IFSsimulator():
def writeOutputs(self, obnum):
"""
Parameters
----------
obnum : TYPE
DESCRIPTION.
Returns
-------
None.
"""
"""
Writes out a FITS file using PyFITS and converts the image array to 16bit unsigned integer as
appropriate for VIS.
......@@ -2760,8 +3299,6 @@ class IFSsimulator():
##########
########## finish header for 0 layer ######################
#############################################################
# 3
# header
# new image HDU, blue channel, layer 1
if HeaderTest == 'yes':
......@@ -2794,8 +3331,8 @@ class IFSsimulator():
hdu_b.header['BZERO'] = (np.float64(32768), '')
hdu_b.header['BUNIT'] = ('ADU', 'physical unit of array values')
###########################################################################
######### instrument information ########################################
#######################################################################
######### instrument information ###################################
if self.source == 'SCI' or self.source == 'COMP':
hdu_b.header['CMIRRPOS'] = (
......@@ -2824,9 +3361,9 @@ class IFSsimulator():
hdu_b.header['IFSSTAT'] = (
np.int32(0), 'IFS components status parameter')
##############################################################################
######################################################################
#################### detector information#############################
#################### detector information##########################
hdu_b.header['CAMERA'] = ('Blue', 'camera of IFS')
hdu_b.header['DETSN'] = ('CCD231-c4-00', 'detector serial number')
......@@ -2888,7 +3425,7 @@ class IFSsimulator():
hdu_b.header['DETTEMP1'] = (np.float32(
0.0), 'detector temperature at EXPT1_'+str(k+1)+' (K)')
###################################################### end revised on 2024.2.27 ###
#######################end revised on 2024.2.27 ###
hdu_b.header['BIN_X'] = (
np.int16(1), 'bin number in X (wavelength)')
......@@ -2944,7 +3481,7 @@ class IFSsimulator():
if self.source == 'LAMP' and self.information['holemask'] == 'yes':
hdu_b.header['Hole'] = ('yes', 'apply hole to LAMP')
######################################################################################################
#####################################################################
#################### red camera ######################
# create a new FITS file, using HDUList instance
......@@ -2980,7 +3517,7 @@ class IFSsimulator():
ofd_r.header['OBJECT'] = (
self.information['name_obj'][:30], 'object name')
ofd_r.header['TARGET'] = (
(self.information['target']), 'target name (hhmmss.s+ddmmss)')
(self.information['target']), 'target name (hhmmss.s+ddmmss)')
ofd_r.header['OBSID'] = (str(obsid), 'observation ID')
......@@ -3142,7 +3679,6 @@ class IFSsimulator():
########## finish header for 0 layer ######################
#############################################################
# 3
# header
# new image HDU, red channel, layer 1
if HeaderTest == 'yes':
......@@ -3175,7 +3711,7 @@ class IFSsimulator():
hdu_r.header['BUNIT'] = ('ADU', 'physical unit of array values')
######### instrument information ######
##hdu_r.header['PMIRRPOS']=(bool(False), 'FSM pointing,T: to MCI, F: not to MCI')
if self.source == 'SCI' or self.source == 'COMP':
hdu_r.header['CMIRRPOS'] = (
......@@ -3204,7 +3740,7 @@ class IFSsimulator():
hdu_r.header['IFSSTAT'] = (
np.int32(0), 'IFS components status parameter')
################### detector information#############################
################### detector information############################
hdu_r.header['CAMERA'] = ('Red', 'camera of IFS')
hdu_r.header['DETSN'] = ('CCD231-c6-00', 'detector serial number')
......@@ -3231,7 +3767,7 @@ class IFSsimulator():
hdu_r.header['OSCAN2'] = (
np.int32(320), 'vertical overscan height, per readout channel')
#####################################################################################################
##############################################################################
# Readout information
frame_time_r = 0.13824 # data frame transfer time in red camera
......@@ -3343,10 +3879,7 @@ class IFSsimulator():
hdu1.header.add_comment(
'========================================================================', after='EPOCH')
#########################################################################################################################
hdu2.header.add_comment(
'========================================================================', after='FITSSWV')
hdu2.header.add_comment('OBJECT INFORMATION', after='FITSSWV')
######################################################################
hdu2.header.add_comment(
'========================================================================', after='FITSSWV')
......@@ -3383,7 +3916,7 @@ class IFSsimulator():
'========================================================================', before='EXPT0_1')
hdulist_b.writeto(self.file_b, output_verify='ignore', checksum=True)
#####################################################################################################
#########################################################################
hdulist_r = fits.HDUList([hdu2, hdu_r])
hdu_r.header.add_comment(
......@@ -3406,7 +3939,7 @@ class IFSsimulator():
hdulist_r.writeto(self.file_r, output_verify='ignore', checksum=True)
##################################################################################
######################################################################
def earthshine(self, theta):
"""
......@@ -3472,7 +4005,22 @@ class IFSsimulator():
###################################################################################
##################################################################################
def CalskyNoise(self, lam):
"""
Parameters
----------
lam : TYPE
DESCRIPTION.
Returns
-------
Ns_skynoise : TYPE
DESCRIPTION.
"""
# calculate the reference flux
# calculate sky noise;
......@@ -3505,6 +4053,19 @@ class IFSsimulator():
##############
def sim_lamp_hole_img(self, exposuretime):
"""
Parameters
----------
exposuretime : TYPE
DESCRIPTION.
Returns
-------
None.
"""
sizeout = 1000 # here image has the pixelscale =0.01 arcsec
################# create the slicer and slit file for high resolution 0.01arcsec simulation ###################
......@@ -3722,14 +4283,14 @@ class IFSsimulator():
img3 = conv
######################## get subimage #####################
######################## get subimag #####################
subimage = galsim.Image(800, 800)
subimage.array[:, :] = img3[int(
sizeout/2)-400:int(sizeout/2)+400, int(sizeout/2)-400:int(sizeout/2)+400]
# fits.writeto('subimage.fits',subimage.array,overwrite=True)
######################## get photons from sub-image #####################
####### get photons from sub-image #####################
subimage.scale = 0.01 # here pixelscale=0.01arcsec
subimage.setOrigin(0, 0)
......@@ -3737,14 +4298,14 @@ class IFSsimulator():
photons = galsim.PhotonArray.makeFromImage(
subimage, max_flux=max(img3.max()/10000.0, 1))
#############################################################################
##############################################################
# do something for each photons;
##############################################################################
#############################################################
idx0 = np.where(photons.flux > 1e-3)
# energy=energy+sum(photons.flux[idx0]) ### totla energy for slice image
p_num = len(idx0[0])
###############################################################################################
##############################################################
# find photons for blue channel, and make the flux multiple the optical and CCD efficiency
......@@ -3832,6 +4393,33 @@ class IFSsimulator():
##################################################################################
def sim_sky_img(self, skyfitsfilename, skyRa, skyDec, sky_rot, telRa, telDec, exposuretime, simnumber):
"""
Parameters
----------
skyfitsfilename : TYPE
DESCRIPTION.
skyRa : TYPE
DESCRIPTION.
skyDec : TYPE
DESCRIPTION.
sky_rot : TYPE
DESCRIPTION.
telRa : TYPE
DESCRIPTION.
telDec : TYPE
DESCRIPTION.
exposuretime : TYPE
DESCRIPTION.
simnumber : TYPE
DESCRIPTION.
Returns
-------
None.
"""
# load the input original data cube simulated by Fengshuai
......@@ -4108,8 +4696,6 @@ class IFSsimulator():
self.log.error('Error in datacube, negative values !!!!!!!')
continue
##############################################################################
# if ilam==100:
# fits.writeto('/media/yan//IFSsim/IFSsim Data/original_Img.fits', Nimg, overwrite=True)
......@@ -4394,6 +4980,21 @@ class IFSsimulator():
#################################################################################################
def sim_calibration_img(self, exposuretime, source):
"""
Parameters
----------
exposuretime : TYPE
DESCRIPTION.
source : TYPE
DESCRIPTION.
Returns
-------
None.
"""
sizeout = 100
################# load the hole mask file ###################
......@@ -4522,7 +5123,7 @@ class IFSsimulator():
slice_image = ns*self.telarea*Width_lambda*10*exptime
##################################################################################################
############################################################################
if slice_image.max() < 0.1:
continue
......@@ -4530,7 +5131,7 @@ class IFSsimulator():
if slice_image.min() < 0:
self.log.error('Error in datacube, negative values !!!!!!!')
continue
##################################################################################################
#############################################################################
# for calibration , there is no primay CSST system, but diffraction effect should consider ;
......@@ -4735,7 +5336,7 @@ class IFSsimulator():
hdu2 = fits.ImageHDU(Wave)
hdu2.header.append(('WaveUnit', 'nm'))
hdu2.header.append(('WaveUnit', 'nm'))
newd = fits.HDUList([hdu1, hdu2])
......@@ -4758,6 +5359,21 @@ class IFSsimulator():
#################################################################################################
def simulate(self, sourcein, simnumber):
"""
Parameters
----------
sourcein : TYPE
DESCRIPTION.
simnumber : TYPE
DESCRIPTION.
Returns
-------
None.
"""
# def simulate(self, skyfitsfilename, skyRa, skyDec, sky_rot, telRa, telDec, exposuretime,simnumber):
"""
Create a single simulated image of a quadrant defined by the configuration file.
......@@ -4789,9 +5405,10 @@ class IFSsimulator():
if simnumber <= 200 and simnumber > 150:
self.information['exptime'] = 1200
self.skyfilepath = self.information['dir_path']+self.information['sky_fitsin']
self.skyfilepath = self.information['dir_path'] + \
self.information['sky_fitsin']
print('self.skyfilepath = ',self.skyfilepath)
print('self.skyfilepath = ', self.skyfilepath)
# self.earthshine_theta=30.0 # in degree
......@@ -5030,7 +5647,29 @@ class IFSsimulator():
############################################################################
############################################################################
def runIFSsim(sourcein, configfile, iLoop, applyhole='no'):
"""
Parameters
----------
sourcein : TYPE
DESCRIPTION.
configfile : TYPE
DESCRIPTION.
iLoop : TYPE
DESCRIPTION.
applyhole : TYPE, optional
DESCRIPTION. The default is 'no'.
Returns
-------
int
DESCRIPTION.
"""
simulate = dict()
simulate[iLoop] = IFSsimulator(configfile)
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment