Skip to content
GitLab
Commit 1b761b45 authored by Chen Wei's avatar Chen Wei
Browse files

Merge branch 'master' into 'main' 

build

See merge request shaosim/ifs!1
parents 15040d3f d2cef509
Show whitespace changes
Inline Side-by-side
.idea/inspectionProfiles/profiles_settings.xml 0 → 100644
View file @ 1b761b45
<component name="InspectionProjectProfileManager">
<settings>
<option name="USE_PROJECT_PROFILE" value="false" />
<version value="1.0" />
</settings>
</component>
\ No newline at end of file
LICENSE 0 → 100644
View file @ 1b761b45
MIT License
Copyright (c) 2022 CSST-L1
Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:
The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.
README.md
View file @ 1b761b45
# Ifs
# csst_ifs_common
......@@ -15,14 +15,14 @@ Already a pro? Just edit this README.md and make it your own. Want to make it ea
```
cd existing_repo
git remote add origin https://csst-tb.bao.ac.cn/code/shaosim/ifs.git
git remote add origin https://csst-tb.bao.ac.cn/code/csst-l1/ifs/csst_ifs_common.git
git branch -M main
git push -uf origin main
```
## Integrate with your tools
- [ ] [Set up project integrations](https://csst-tb.bao.ac.cn/code/shaosim/ifs/-/settings/integrations)
- [ ] [Set up project integrations](http://10.3.10.28/code/csst-l1/ifs/csst_ifs_common/-/settings/integrations)
## Collaborate with your team
......
csst_ifs_sim/CTI/CTI.py 0 → 100644
View file @ 1b761b45
"""
Charge Transfer Inefficiency
============================
This file contains a simple class to run a CDM03 CTI model developed by Alex Short (ESA).
This now contains both the official CDM03 and a new version that allows different trap
parameters in parallel and serial direction.
:requires: NumPy
:requires: CDM03 (FORTRAN code, f2py -c -m cdm03bidir cdm03bidir.f90)
:version: 0.35
"""
import numpy as np
# try:
# import cdm03bidir
# #import cdm03bidirTest as cdm03bidir #for testing purposes only
# except ImportError:
# print('import CTI module')
# #print ('No CDM03bidir module available, please compile it: f2py -c -m cdm03bidir cdm03bidir.f90')
#CDM03bidir
class CDM03bidir():
"""
Class to run CDM03 CTI model, class Fortran routine to perform the actual CDM03 calculations.
:param settings: input parameters
:type settings: dict
:param data: input data to be radiated
:type data: ndarray
:param log: instance to Python logging
:type log: logging instance
"""
def __init__(self, settings, data, log=None):
"""
Class constructor.
:param settings: input parameters
:type settings: dict
:param data: input data to be radiated
:type data: ndarray
:param log: instance to Python logging
:type log: logging instance
"""
self.data = data
self.values = dict(quads=(0,1,2,3), xsize=2048, ysize=2066, dob=0.0, rdose=8.0e9)
self.values.update(settings)
self.log = log
self._setupLogger()
#default CDM03 settings
self.params = dict(beta_p=0.6, beta_s=0.6, fwc=200000., vth=1.168e7, vg=6.e-11, t=20.48e-3,
sfwc=730000., svg=1.0e-10, st=5.0e-6, parallel=1., serial=1.)
#update with inputs
self.params.update(self.values)
#read in trap information
trapdata = np.loadtxt(self.values['parallelTrapfile'])
if trapdata.ndim > 1:
self.nt_p = trapdata[:, 0]
self.sigma_p = trapdata[:, 1]
self.taur_p = trapdata[:, 2]
else:
#only one trap species
self.nt_p = [trapdata[0],]
self.sigma_p = [trapdata[1],]
self.taur_p = [trapdata[2],]
trapdata = np.loadtxt(self.values['serialTrapfile'])
if trapdata.ndim > 1:
self.nt_s = trapdata[:, 0]
self.sigma_s = trapdata[:, 1]
self.taur_s = trapdata[:, 2]
else:
#only one trap species
self.nt_s = [trapdata[0],]
self.sigma_s = [trapdata[1],]
self.taur_s = [trapdata[2],]
#scale thibaut's values
if 'thibaut' in self.values['parallelTrapfile']:
self.nt_p /= 0.576 #thibaut's values traps / pixel
self.sigma_p *= 1.e4 #thibaut's values in m**2
if 'thibaut' in self.values['serialTrapfile']:
self.nt_s *= 0.576 #thibaut's values traps / pixel #should be division?
self.sigma_s *= 1.e4 #thibaut's values in m**2
def _setupLogger(self):
"""
Set up the logger.
"""
self.logger = True
if self.log is None:
self.logger = False
def radiateFullCCD(self):
"""
This routine allows the whole CCD to be run through a radiation damage mode.
The routine takes into account the fact that the amplifiers are in the corners
of the CCD. The routine assumes that the CCD is using four amplifiers.
There is an excess of .copy() calls, which should probably be cleaned up. However,
given that I had problem with the Fortran code, I have kept the calls. If memory
becomes an issue then this should be cleaned.
:return: radiation damaged image
:rtype: ndarray
"""
ydim, xdim = self.data.shape
out = np.zeros((xdim, ydim))
#transpose the data, because Python has different convention than Fortran
data = self.data.transpose().copy()
for quad in self.values['quads']:
if self.logger:
self.log.info('Adding CTI to Q%i' % quad)
if quad == 0:
d = data[0:self.values['xsize'], 0:self.values['ysize']].copy()
tmp = self.applyRadiationDamage(d, iquadrant=quad).copy()
out[0:self.values['xsize'], 0:self.values['ysize']] = tmp
elif quad == 1:
d = data[self.values['xsize']:, :self.values['ysize']].copy()
tmp = self.applyRadiationDamage(d, iquadrant=quad).copy()
out[self.values['xsize']:, :self.values['ysize']] = tmp
elif quad == 2:
d = data[:self.values['xsize'], self.values['ysize']:].copy()
tmp = self.applyRadiationDamage(d, iquadrant=quad).copy()
out[:self.values['xsize'], self.values['ysize']:] = tmp
elif quad == 3:
d = data[self.values['xsize']:, self.values['ysize']:].copy()
tmp = self.applyRadiationDamage(d, iquadrant=quad).copy()
out[self.values['xsize']:, self.values['ysize']:] = tmp
else:
print( 'ERROR -- too many quadrants!!' )
self.log.error('Too many quadrants! This method allows only four quadrants.')
return out.transpose()
def radiateFullCCD2(self):
"""
This routine allows the whole CCD to be run through a radiation damage mode.
The routine takes into account the fact that the amplifiers are in the corners
of the CCD. The routine assumes that the CCD is using four amplifiers.
There is an excess of .copy() calls, which should probably be cleaned up. However,
given that I had problem with the Fortran code, I have kept the calls. If memory
becomes an issue then this should be cleaned.
:return: radiation damaged image
:rtype: ndarray
"""
ydim, xdim = self.data.shape
out = np.empty((ydim, xdim))
#transpose the data, because Python has different convention than Fortran
data = self.data.copy()
for quad in self.values['quads']:
if self.logger:
self.log.info('Adding CTI to Q%i' % quad)
if quad == 0:
d = data[:self.values['ysize'], :self.values['xsize']].copy()
tmp = self.applyRadiationDamage(d, iquadrant=quad).copy()
out[:self.values['ysize'], :self.values['xsize']] = tmp
elif quad == 1:
d = data[:self.values['ysize'], self.values['xsize']:].copy()
tmp = self.applyRadiationDamage(d, iquadrant=quad).copy()
out[:self.values['ysize'], self.values['xsize']:] = tmp
elif quad == 2:
d = data[self.values['ysize']:, :self.values['xsize']].copy()
tmp = self.applyRadiationDamage(d, iquadrant=quad).copy()
out[self.values['ysize']:, :self.values['xsize']] = tmp
elif quad == 3:
d = data[self.values['ysize']:, self.values['xsize']:].copy()
tmp = self.applyRadiationDamage(d, iquadrant=quad).copy()
out[self.values['ysize']:, self.values['xsize']:] = tmp
else:
print( 'ERROR -- too many quadrants!!')
self.log.error('Too many quadrants! This method allows only four quadrants.')
return out
def applyRadiationDamage(self, data, iquadrant=0):
"""
Apply radian damage based on FORTRAN CDM03 model. The method assumes that
input data covers only a single quadrant defined by the iquadrant integer.
:param data: imaging data to which the CDM03 model will be applied to.
:type data: ndarray
:param iquandrant: number of the quadrant to process
:type iquandrant: int
cdm03 - Function signature::
sout = cdm03(sinp,iflip,jflip,dob,rdose,in_nt,in_sigma,in_tr,[xdim,ydim,zdim])
Required arguments:
sinp : input rank-2 array('d') with bounds (xdim,ydim)
iflip : input int
jflip : input int
dob : input float
rdose : input float
in_nt : input rank-1 array('d') with bounds (zdim)
in_sigma : input rank-1 array('d') with bounds (zdim)
in_tr : input rank-1 array('d') with bounds (zdim)
Optional arguments:
xdim := shape(sinp,0) input int
ydim := shape(sinp,1) input int
zdim := len(in_nt) input int
Return objects:
sout : rank-2 array('d') with bounds (xdim,ydim)
.. Note:: Because Python/NumPy arrays are different row/column based, one needs
to be extra careful here. NumPy.asfortranarray will be called to get
an array laid out in Fortran order in memory. Before returning the
array will be laid out in memory in C-style (row-major order).
:return: image that has been run through the CDM03 model
:rtype: ndarray
"""""
#return data
iflip = iquadrant / 2
jflip = iquadrant % 2
params = [self.params['beta_p'], self.params['beta_s'], self.params['fwc'], self.params['vth'],
self.params['vg'], self.params['t'], self.params['sfwc'], self.params['svg'], self.params['st'],
self.params['parallel'], self.params['serial']]
if self.logger:
self.log.info('nt_p=' + str(self.nt_p))
self.log.info('nt_s=' + str(self.nt_s))
self.log.info('sigma_p= ' + str(self.sigma_p))
self.log.info('sigma_s= ' + str(self.sigma_s))
self.log.info('taur_p= ' + str(self.taur_p))
self.log.info('taur_s= ' + str(self.taur_s))
self.log.info('dob=%f' % self.values['dob'])
self.log.info('rdose=%e' % self.values['rdose'])
self.log.info('xsize=%i' % data.shape[1])
self.log.info('ysize=%i' % data.shape[0])
self.log.info('quadrant=%i' % iquadrant)
self.log.info('iflip=%i' % iflip)
self.log.info('jflip=%i' % jflip)
#################################################################################
###modify
import sys
#sys.path.append('../so')
from ifs_so import cdm03bidir
# from ifs_so.cdm03.cpython-38-x86_64-linux-gnu import cdm03bidir
# import cdm03bidir
CTIed = cdm03bidir.cdm03(np.asfortranarray(data),
jflip, iflip,
self.values['dob'], self.values['rdose'],
self.nt_p, self.sigma_p, self.taur_p,
self.nt_s, self.sigma_s, self.taur_s,
params,
[data.shape[0], data.shape[1], len(self.nt_p), len(self.nt_s), len(self.params)])
return np.asanyarray(CTIed)
#################################################################################################################
class CDM03():
"""
Class to run CDM03 CTI model, class Fortran routine to perform the actual CDM03 calculations.
:param data: input data to be radiated
:type data: ndarray
:param input: input parameters
:type input: dictionary
:param log: instance to Python logging
:type log: logging instance
"""
def __init__(self, input, data, log=None):
"""
Class constructor.
:param data: input data to be radiated
:type data: ndarray
:param input: input parameters
:type input: dictionary
:param log: instance to Python logging
:type log: logging instance
"""
try:
import cdm03
except ImportError:
print( 'No CDM03 module available, please compile it: f2py -c -m cdm03 cdm03.f90')
self.data = data
self.values = dict(quads=(0,1,2,3), xsize=2048, ysize=2066, dob=0.0, rdose=8.0e9)
self.values.update(input)
self.log = log
self._setupLogger()
def _setupLogger(self):
"""
Set up the logger.
"""
self.logger = True
if self.log is None:
self.logger = False
def radiateFullCCD(self):
"""
This routine allows the whole CCD to be run through a radiation damage mode.
The routine takes into account the fact that the amplifiers are in the corners
of the CCD. The routine assumes that the CCD is using four amplifiers.
There is an excess of .copy() calls, which should probably be cleaned up. However,
given that I had problem with the Fortran code, I have kept the calls. If memory
becomes an issue then this should be cleaned.
:return: radiation damaged image
:rtype: ndarray
"""
ydim, xdim = self.data.shape
out = np.zeros((xdim, ydim))
#transpose the data, because Python has different convention than Fortran
data = self.data.transpose().copy()
for quad in self.values['quads']:
if self.logger:
self.log.info('Adding CTI to Q%i' % quad)
if quad == 0:
d = data[0:self.values['xsize'], 0:self.values['ysize']].copy()
tmp = self.applyRadiationDamage(d, iquadrant=quad).copy()
out[0:self.values['xsize'], 0:self.values['ysize']] = tmp
elif quad == 1:
d = data[self.values['xsize']:, :self.values['ysize']].copy()
tmp = self.applyRadiationDamage(d, iquadrant=quad).copy()
out[self.values['xsize']:, :self.values['ysize']] = tmp
elif quad == 2:
d = data[:self.values['xsize'], self.values['ysize']:].copy()
tmp = self.applyRadiationDamage(d, iquadrant=quad).copy()
out[:self.values['xsize'], self.values['ysize']:] = tmp
elif quad == 3:
d = data[self.values['xsize']:, self.values['ysize']:].copy()
tmp = self.applyRadiationDamage(d, iquadrant=quad).copy()
out[self.values['xsize']:, self.values['ysize']:] = tmp
else:
print ('ERROR -- too many quadrants!!')
self.log.error('Too many quadrants! This method allows only four quadrants.')
return out.transpose()
def radiateFullCCD2(self):
"""
This routine allows the whole CCD to be run through a radiation damage mode.
The routine takes into account the fact that the amplifiers are in the corners
of the CCD. The routine assumes that the CCD is using four amplifiers.
There is an excess of .copy() calls, which should probably be cleaned up. However,
given that I had problem with the Fortran code, I have kept the calls. If memory
becomes an issue then this should be cleaned.
:return: radiation damaged image
:rtype: ndarray
"""
ydim, xdim = self.data.shape
out = np.empty((ydim, xdim))
#transpose the data, because Python has different convention than Fortran
data = self.data.copy()
for quad in self.values['quads']:
if self.logger:
self.log.info('Adding CTI to Q%i' % quad)
if quad == 0:
d = data[:self.values['ysize'], :self.values['xsize']].copy()
tmp = self.applyRadiationDamage(d, iquadrant=quad).copy()
out[:self.values['ysize'], :self.values['xsize']] = tmp
elif quad == 1:
d = data[:self.values['ysize'], self.values['xsize']:].copy()
tmp = self.applyRadiationDamage(d, iquadrant=quad).copy()
out[:self.values['ysize'], self.values['xsize']:] = tmp
elif quad == 2:
d = data[self.values['ysize']:, :self.values['xsize']].copy()
tmp = self.applyRadiationDamage(d, iquadrant=quad).copy()
out[self.values['ysize']:, :self.values['xsize']] = tmp
elif quad == 3:
d = data[self.values['ysize']:, self.values['xsize']:].copy()
tmp = self.applyRadiationDamage(d, iquadrant=quad).copy()
out[self.values['ysize']:, self.values['xsize']:] = tmp
else:
print ('ERROR -- too many quadrants!!')
self.log.error('Too many quadrants! This method allows only four quadrants.')
return out
def applyRadiationDamage(self, data, iquadrant=0):
"""
Apply radian damage based on FORTRAN CDM03 model. The method assumes that
input data covers only a single quadrant defined by the iquadrant integer.
:param data: imaging data to which the CDM03 model will be applied to.
:type data: ndarray
:param iquandrant: number of the quadrant to process
:type iquandrant: int
cdm03 - Function signature::
sout = cdm03(sinp,iflip,jflip,dob,rdose,in_nt,in_sigma,in_tr,[xdim,ydim,zdim])
Required arguments:
sinp : input rank-2 array('d') with bounds (xdim,ydim)
iflip : input int
jflip : input int
dob : input float
rdose : input float
in_nt : input rank-1 array('d') with bounds (zdim)
in_sigma : input rank-1 array('d') with bounds (zdim)
in_tr : input rank-1 array('d') with bounds (zdim)
Optional arguments:
xdim := shape(sinp,0) input int
ydim := shape(sinp,1) input int
zdim := len(in_nt) input int
Return objects:
sout : rank-2 array('d') with bounds (xdim,ydim)
.. Note:: Because Python/NumPy arrays are different row/column based, one needs
to be extra careful here. NumPy.asfortranarray will be called to get
an array laid out in Fortran order in memory. Before returning the
array will be laid out in memory in C-style (row-major order).
:return: image that has been run through the CDM03 model
:rtype: ndarray
"""
#read in trap information
trapdata = np.loadtxt(self.values['trapfile'])
nt = trapdata[:, 0]
sigma = trapdata[:, 1]
taur = trapdata[:, 2]
iflip = iquadrant / 2
jflip = iquadrant % 2
if self.logger:
self.log.info('nt=' + str(nt))
self.log.info('sigma= ' + str(sigma))
self.log.info('taur= ' + str(taur))
self.log.info('dob=%f' % self.values['dob'])
self.log.info('rdose=%e' % self.values['rdose'])
self.log.info('xsize=%i' % data.shape[1])
self.log.info('ysize=%i' % data.shape[0])
self.log.info('quadrant=%i' % iquadrant)
self.log.info('iflip=%i' % iflip)
self.log.info('jflip=%i' % jflip)
# #call Fortran routine
# CTIed = cdm03.cdm03(np.asfortranarray(data),
# iflip, jflip,
# self.values['dob'], self.values['rdose'],
# nt, sigma, taur)
###modify
import sys
sys.path.append('../CTI')
import cdm03
#################################################################################
CTIed = cdm03.cdm03(np.asfortranarray(data),
jflip, iflip,
self.values['dob'], self.values['rdose'],
nt,sigma,taur )
return np.asanyarray(CTIed)
#################################################################################################################
csst_ifs_sim/csst_ifs_sim.py 0 → 100644
View file @ 1b761b45
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
@author: yan
"""
"""
The CSST IFS Image Simulator
=============================================
This file contains an image simulator for the CSST IFS.
The approximate sequence of events in the simulator is as follows:
#. Read in a configuration file, which defines for example,
detector characteristics (bias, dark and readout noise, gain,
plate scale and pixel scale, oversampling factor, exposure time etc.).
#. Read in another file containing charge trap definitions (for CTI modelling).
#. Read in a file defining the cosmic rays (trail lengths and cumulative distributions).
#. Load the wavefront aberration data used to calculate PSF with defined wavelength and field of view.
#. Apply calibration unit flux to mimic flat field exposures [optional].
#. Apply a multiplicative flat-field map to emulate pixel-to-pixel non-uniformity [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 photon (Poisson) noise [optional]
#. Add cosmetic defects from an input file [optional].
#. Add overscan regions in the serial direction [optional].
#. Apply the CDM03 radiation damage model [optional].
#. Apply 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).
#. Finally the simulated image is converted to a FITS file.
:version: 2023.04
Future Work
-----------
.. todo::
#.
#.
#. ......
Contact Information: zhaojunyan@shao.ac.cn
-------
"""
from scipy.integrate import simps
import pandas as pd
from datetime import datetime, timedelta
import julian
from astropy.coordinates import get_sun
from astropy.time import Time
from astropy import units as u
from astropy.coordinates import SkyCoord, EarthLocation
from scipy import interpolate
import numpy as np
import scipy.io as sio
from astropy.io import fits
from scipy.signal import fftconvolve
import galsim
from scipy import ndimage
import cmath
import os, sys
import configparser as ConfigParser
from optparse import OptionParser
sys.path.append('../')
from CTI import CTI
from support import logger as lg
from support import cosmicrays
from support import IFSinstrumentModel
###############################################################################
def str2time(strTime):
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 to mjd
def time2mjd(dateT):
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):
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++2400000.5
#mjd to datetime
def mjd2time(mjd):
t0=datetime(1858,11,17,0,0,0,0)#
return t0+timedelta(days=mjd)
def jd2time(jd):
mjd=jd2mjd(jd)
return mjd2time(mjd)
#mjd to jd
def mjd2jd(mjd):
return mjd+2400000.5
def jd2mjd(jd):
return jd-2400000.5
def dt2hmd(dt):
## dt is datetime
hour=dt.hour
minute=dt.minute
second=dt.second
if hour<10:
str_h='0'+str(hour)
else:
str_h=str(hour)
if minute<10:
str_m='0'+str(minute)
else:
str_m=str(minute)
if second<10:
str_d='0'+str(second+dt.microsecond*1e-6)
else:
str_d=str(second+dt.microsecond*1e-6)
return str_h+':'+str_m+':'+str_d
###############################################################################
###############################################################################
def deg2HMS(ra, dec, rou=False):
'''convert deg to ra's HMS or dec's DHS'''
RA, DEC, rs, ds = '000000', '000000', '', ''
if dec:
if str(dec)[0] == '-':
ds, dec = '-', abs(dec)
deg = int(dec)
decM = abs(int((dec-deg)*60))
if rou:
decS = round((abs((dec-deg)*60)-decM)*60,1)
else:
decS = int((abs((dec-deg)*60)-decM)*60)
if deg ==0:
deg ="00"
elif deg <10:
deg = "0%s"%deg
if decM ==0:
decM ="00"
elif decM <10:
decM = "0%s"%decM
if decS ==0:
decS ="00"
elif decS <10:
decS = "0%s"%decS
DEC = '{0}{1}{2}{3}'.format(ds, deg, decM, decS)
if ra:
if str(ra)[0] == '-':
rs, ra = '-', abs(ra)
raH = int(ra/15)
raM = int(((ra/15)-raH)*60)
if rou:
raS = round(((((ra/15)-raH)*60)-raM)*60,1)
else:
raS = int(((((ra/15)-raH)*60)-raM)*60)
if raH ==0:
raH = "00"
elif raH <10:
raH = "0%s"%raH
if raM ==0:
raM = "00"
elif raM <10:
raM = "0%s"%raM
if raS ==0:
raS = "00"
elif raS <10:
raS = "0%s"%raS
RA = '{0}{1}{2}{3}'.format(rs, raH, raM, raS)
if ds=='-':
return RA+DEC
else:
return RA+'+'+DEC
###############################################################################
###########################################################
def beta_angle( x_sat, y_sat, z_sat, vx_sat, vy_sat, vz_sat, ra_obj, dec_obj):
# get the vector for next motion
x_sat2 = x_sat + vx_sat
y_sat2 = y_sat + vy_sat
z_sat2 = z_sat + vz_sat
vector1 = np.stack((x_sat, y_sat, z_sat), axis=0)
vector2 = np.stack((x_sat2, y_sat2, z_sat2), axis=0)
vector_normal = np.cross(vector1, vector2)
location_normal = EarthLocation.from_geocentric(vector_normal[0],
vector_normal[1],
vector_normal[2],
unit=u.km)
radec_normal = SkyCoord(ra=location_normal.lon,
dec=location_normal.lat,
frame='gcrs')
lb_normal = radec_normal.transform_to('geocentrictrueecliptic')
radec_sun = SkyCoord(ra=ra_obj*u.degree, dec=dec_obj*u.degree, frame='gcrs')
lb_sun = radec_sun.transform_to('geocentrictrueecliptic')
# get the angle between normal and solar
angle = 90 - lb_normal.separation(lb_sun).degree
return angle
##########################################################
##########################################################
def LSR_velocity(ra,dec,velocity,Obstime):
local = EarthLocation.from_geodetic(lat=22.349368*u.deg, lon=113.584068*u.deg, height=10*u.m)
### convert ra and dec to
source = SkyCoord(ra*u.deg, dec*u.deg,frame='icrs', unit=(u.hourangle,u.deg) )
l = source.galactic.l.deg
b = source.galactic.b.deg
c = SkyCoord(l=l*u.degree, b=b*u.degree, frame='galactic')
c_icrs = c.transform_to('icrs')
barycorr = c_icrs.radial_velocity_correction(obstime=Time(Obstime), location=local)
velocity = velocity + barycorr.value/1000
#print(barycorr.value/1000)
l = l * np.pi / 180
b = b * np.pi / 180
return velocity + 9 * np.cos(l) * np.cos(b) + 12 * np.sin(l) * np.cos(b) + 7 * np.sin(b)
###############################################################################
def rotation_yan(x0,y0,x1,y1,angle):
#% point A (x0,y0)
#% point B(x1,y1)
#% roate angle ,point B rotate with point A
alpha=angle/180*3.1415926 #% in radians
x2=(x1-x0)*np.cos(alpha)-(y1-y0)*np.sin(alpha)+x0
y2=(x1-x0)*np.sin(alpha)+(y1-y0)*np.cos(alpha)+y0
return x2, y2
######################################################################
def centroid(data):
h,w = np.shape(data)
x = np.arange(0,w)
y = np.arange(0,h)
x1=np.ones((1,h))
y1=np.ones((w,1))
cx=(np.dot(np.dot(x1,data),y)) /(np.dot(np.dot(x1,data),y1))
cy=(np.dot(np.dot(x ,data),y1)) /(np.dot(np.dot(x1,data),y1))
return np.float64(cx),np.float64(cy)
####################################################################
def SpecCube2photon(data,wave):
# 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;
planckh= 6.62620*10**-27 # % erg s;
cc=2.99792458*10**17 # in nm/s
# fov2=0.01 # in arcsec^2
# telarea=3.1415926*100*100 # in cm^2，
fluxlam=1e-17*data # convert original unit to unit of erg/s/A/cm^2
lam=wave # wave in nm ;;
ephoton=planckh*cc/lam # in erg/photon,
Nphoton =fluxlam/ephoton # in unit of photons/cm2/s/A
return Nphoton
######################################################################
######################################################################
def zero_pad(image, shape, position='corner'):
"""
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
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
def anySampledPSF(wavefront, pupil , Q , sizeout):
'''
written on 2020.12.04 by Yan @Shao
% wavefront sampled in the united circle;
% pupil function samped as wavefront;
% sample ratio Q=lambda/D/pixel;
% lambda is wavelength;
% D is diameter;
% pixel is in arcsec;
% or sample ratio Q=lambda*f/D/pixelsize
% f is focal length;
% pixelsize is the actural size of the detector;
% make sure all the varia have the same unit;
% the returned PSF has sum value of 1
'''
m,n = np.shape(wavefront)
phase=2.0*np.pi*wavefront #% the phase of the input wavefront;
####generalized pupil function of channel1;
if Q>=1:
pk1=zero_pad(pupil*np.exp(cmath.sqrt(-1)*(phase)),( int(np.floor(m*Q)),int(np.floor(m*Q))), position='corner')
psf=abs(np.fft.fftshift(np.fft.fft2(pk1)))**2
#psf=psf/psf.sum()
else: #%%% case: Q<1
# %fft method
if Q>=0.5: #% in this case Q<1 and Q>=0.5
pk1=zero_pad(pupil*np.exp(cmath.sqrt(-1)*(phase)), (int(2*np.floor(m*Q)),int(2*np.floor(m*Q))), position='corner')
mypsf=np.fft.fft2(pk1);
t=mypsf[0::2,0::2];
psf=abs(np.fft.fftshift(t))**2
#psf=mypsf3/mypsf3.sum();
else:
if Q>=0.25:
pk1=zero_pad(pupil*np.exp(cmath.sqrt(-1)*(phase)),( int(4*np.floor(m*Q)),int(4*np.floor(m*Q))), position='corner')
mypsf=np.fft.fft2(pk1);
t=mypsf[0::4,0::4];
psf=abs(np.fft.fftshift(t))**2
#psf=mypsf3/mypsf3.sum();
else:
print('---- Q<0.25 , stop running')
sys.exit(0)
########################################################################################
mm,nn=np.shape(psf)
if np.mod(sizeout,2)==0:
Nx=sizeout/2-1
else:
Nx=(sizeout+1)/2-1
if max(mm,nn)<sizeout:
psfout=np.zeros((sizeout,sizeout))
[cx,cy]=np.where(psf==np.max(psf))
s=np.ceil(min(mm,nn)/2)-abs(cx-cy)-1
Nx=int(Nx)
cx=int(cx)
cy=int(cy)
s=int(s)
psfout[Nx-s:Nx+s-1,Nx-s:Nx+s-1]=psf[cx-s:cx+s-1,cy-s:cy+s-1]
else:
psfout=np.zeros((sizeout,sizeout));
[cx,cy]=np.where(psf==np.max(psf))
s1=int(np.ceil(min(mm,nn)/2)-abs(cx-cy)-1)
s2=int(sizeout/2)
s=min(s1,s2)
s=s-1
Nx=int(Nx)
cx=int(cx)
cy=int(cy)
psfout[Nx-s:Nx+s-1,Nx-s:Nx+s-1]=psf[cx-s:cx+s-1,cy-s:cy+s-1]
####
cx,cy=centroid(psfout)
Nt=sizeout
#psf=ndimage.shift(psfout,[Nt/2-cy, Nt/2-cx],order=1, mode='nearest' ) ## for my fft method
psf=ndimage.shift(psfout,[Nt/2-cy-1, Nt/2-cx-1],order=1, mode='nearest' ) #for convolve method
psf=psf/psf.sum()
return psf
##############################################################################
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]) #垂直方向
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):
#wave is the wavelength in nm;
dydl=np.array([-640.0239901372472, 0.0018, 0.00048, 0.0000028, -0.0000000000007, 0.0 , 0.0] ) #色散方向
dxdl=0.00325*np.array([-1638.8, 4.0e-2, 5.500e-3, -5.2e-10, 1.7000e-10, 7.1e-13, -5.16e-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 getSpectrum(Spectrum0,lam,sigma):
#
wave=Spectrum0[:,0] ## wavelength；
Qt =Spectrum0[:,1] ##intensity；
d=abs(lam-wave);
if min(d)>3*sigma:
SpmOut=1.0/20000
else:
column=np.where(d==min(d))
SpmOut=Qt[column[0]]*np.exp(-(lam-wave[column[0]])**2/2/sigma**2)
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():
"""
CSST IFS Simulator
The image that is being build is in::
self.image
self.image_b blue channel
self.image_r red channel
:param opts: OptionParser instance
:type opts: OptionParser instance
"""
def __init__(self, opts):
"""
Class Constructor.
:param opts: OptionParser instance
:type opts: OptionParser instance
"""
####################################
self.configfile = opts.configfile
if opts.section is None:
####self.section = 'DEFAULT'
self.section = 'TEST' #####simulation section;
else:
self.section = opts.section
#load instrument model, these values are also stored in the FITS header
self.information = IFSinstrumentModel.IFSinformation()
#update settings with defaults
self.information.update(dict(quadrant=int(0),
ccdx=int(0),
ccdy=int(0),
psfoversampling=1.0,
ccdxgap=1.643,
ccdygap=8.116,
fullwellcapacity=90000,
readouttime=4,
rdose=8.0e9,
coveringfraction=0.5 ))
######################################################################
self.configure() #print the configfile name and path;
self.information.update(dict(
cosmicraylengths =self.information['indata_path']+'/cdf_cr_length.dat',
cosmicraydistance=self.information['indata_path']+'/cdf_cr_total.dat',
parallelTrapfile =self.information['indata_path']+'/cdm_euclid_parallel.dat',
serialTrapfile =self.information['indata_path']+'/cdm_euclid_serial.dat',
cosmeticsfile_b =self.information['indata_path']+'/Cosmetics_b.txt',
cosmeticsfile_r =self.information['indata_path']+'/Cosmetics_r.txt' ))
# public system information
self.pixelscale=0.1 # the pixel size is 0.1 arcsec
self.pixelsize=15 # the pixel physical size is 15 micron
self.Fnum=15.469875 # the f number= focal_length/D;
self.telD=2 # tht telescope size is 2 meter, diamter size;
self.telarea=3.1415926*(self.telD/2)**2*100*100 # telescope square in cm^2
self.fov2=0.01 # the pixel square
self.planckh= 6.62620*1e-27 # % erg s;
self.cc=2.99792458*1e17 # in nm/s
self.light_FWHM=0.175
self.HgArsigma=self.light_FWHM/2.35; ## sigma value of the light source;
#### load system optical and CCD efficiency data
matfn0=self.information['indata_path']+'/TotalQ200923.mat'
da0=sio.loadmat(matfn0)
self.optical_blue_Q=da0['opticalQ_1'] # optical efficiency of blue channel
self.optical_red_Q =da0['opticalQ_2'] # optical efficiency of red channel
self.CCD_Qe=da0['ccdQE'] # red channel
############################################################# load all useful data;
############### load wavefront data;
matfn2=self.information['indata_path']+'/opd_638nm.mat'
da2=sio.loadmat(matfn2)
opd0=da2['opd'] # opd unit is in meter
self.opd0=opd0*1e9 # convert opd to nm
self.pupil=abs(opd0)>0.0
####################
da=fits.open(self.information['indata_path']+'/opd1.fits')
self.opd1=da[0].data
da=fits.open(self.information['indata_path']+'/opd2.fits')
self.opd2=da[0].data
da=fits.open(self.information['indata_path']+'/opd3.fits')
self.opd3=da[0].data
################################################################################
## slice information; the slice is put in vertial direction;
# the silce long size is 64 pixels
slice_blue=dict()
slice_red =dict()
randRedpos=np.array([0.52835362, 1.1105926 , -0.81667794, 0.88860884, -2.78092636,
-0.15810022, -1.56726852, -0.71601855, -1.31768647, 1.73107581,
0.4933876 , 2.83673377, 0.22226286, -0.02634998, 0.35539383,
-0.91989574, -2.44856212, 0.91020484, -3.03097852, -1.11638071,
1.28360669, -0.12521128, -0.3907698 , 0.70183575, 1.00578099,
1.67339662, 0.18067182, -0.56303075, 0.40959616, 1.45518379,
-0.93194744, 0.41492972])
randBluepos=np.array([ 0.97449008, -0.21371406, -1.62513338, -3.06938604, 1.72615283,
0.73333374, 0.80923919, -0.9418576 , -0.16806578, -1.04416631,
2.20155068, -0.0900156 , 0.07597706, 0.76208373, 0.29426592,
-0.89434703, 0.34017826, 1.16936499, 0.10977829, -1.31179539,
-0.50859787, -1.01891651, -0.95791581, -1.53018041, 0.88358827,
0.07837641, -0.86157157, -1.18070438, 0.53970599, 1.4913733 ,
2.10938775, 1.23213412])
#####################
for i in range(32):
if i==0:
slice_blue['py']=np.zeros(32)
slice_blue['px']=np.zeros(32)
slice_red['py']=np.zeros(32)
slice_red['px']=np.zeros(32)
if i<16:
slice_blue['py'][i]=50+randBluepos[i]*4
slice_blue['px'][i]=3.55/0.015*i+166.0
slice_red['py'][i]=50+randBluepos[i]*4
slice_red['px'][i]=3.55/0.015*i+1190.0
else:
slice_blue['py'][i]=50+250+randBluepos[i]*4
slice_blue['px'][i]=3.55/0.015*(i-16)+166.0+118
slice_red['py'][i]=50+250+randBluepos[i]*4
slice_red['px'][i]=3.55/0.015*(i-16)+1190.0+118
#######
self.slice_blue=slice_blue
self.slice_red=slice_red
###############################################################################
maskSlice=dict()
maskSlit=dict()
sizeout=100
for k in range(32):
maskSlice[str(k)]=np.zeros((sizeout,sizeout))
maskSlice[str(k)][2*k+18:2*k+20, int(sizeout/2)-32:int(sizeout/2)+32]=1
maskSlit[str(k)]=np.zeros((sizeout,sizeout))
maskSlit[str(k)][2*k+18:2*k+20, int(sizeout/2)-37:int(sizeout/2)+37]=1
self.maskSlice=maskSlice
self.maskSlit =maskSlit
################################################################################################
############################################################################3
def readConfigs(self):
"""
Reads the config file information using configParser and sets up a logger.
"""
self.config = ConfigParser.RawConfigParser()
self.config.read_file(open(self.configfile))
###############################################################################
def processConfigs(self):
"""
Processes configuration information and save the information to a dictionary self.information.
The configuration file may look as follows::
[TEST]
For explanation of each field, see /data/test.config. Note that if an input field does not exist,
then the values are taken from the default instrument model as described in
support.IFSinstrumentModel.VISinformation(). Any of the defaults can be overwritten by providing
a config file with a correct field name.
"""
#parse options and update the information dictionary
options = self.config.options(self.section)
settings = {}
for option in options:
try:
settings[option] = self.config.getint(self.section, option)
except ValueError:
try:
settings[option] = self.config.getfloat(self.section, option)
except ValueError:
settings[option] = self.config.get(self.section, option)
self.information.update(settings)
self.cosmicRays = self.config.getboolean(self.section, 'cosmicRays')
self.darknoise = self.config.getboolean(self.section, 'darknoise')
self.cosmetics = self.config.getboolean(self.section, 'cosmetics')
self.radiationDamage = self.config.getboolean(self.section, 'radiationDamage')
self.bleeding = self.config.getboolean(self.section, 'bleeding')
self.skyback = self.config.getboolean(self.section, 'skyback')
#these don't need to be in the config file
try:
self.nonlinearity = self.config.getboolean(self.section, 'nonlinearity')
except:
self.nonlinearity = False
try:
self.flatfieldM = self.config.getboolean(self.section, 'flatfieldM')
except:
self.flatfieldM = False
try:
self.readoutNoise = self.config.getboolean(self.section, 'readoutNoise')
except:
self.readoutNoise = True
####################################################################
self.booleans = dict(nonlinearity=self.nonlinearity,
flatfieldM=self.flatfieldM,
cosmicRays=self.cosmicRays,
darknoise=self.darknoise,
cosmetics=self.cosmetics,
radiationDamage=self.radiationDamage,
bleeding=self.bleeding,
skyback =self.skyback)
#####################################################################
now=datetime.now()
result_day=now.strftime("%Y-%m-%d")
self.result_path=self.information['result_path']+'/'+result_day
if os.path.isdir(self.result_path)==False:
os.mkdir(self.result_path)
os.mkdir(self.result_path+'/log_file')
os.mkdir(self.result_path+'/sky_Data')
now=datetime.now()
data_time=now.strftime("%Y-%m-%d-%H-%M-%S")
self.log = lg.setUpLogger(self.result_path+'/log_file/IFS_'+'_'+data_time+'.log')
self.log.info('STARTING A NEW SIMULATION')
self.log.info(self.information)
return
#########################################################################################
#######################################################################
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']))
crLengths = np.loadtxt(self.information['cosmicraylengths'])
crDists = np.loadtxt(self.information['cosmicraydistance'])
self.cr = dict(cr_u=crLengths[:, 0], cr_cdf=crLengths[:, 1], cr_cdfn=np.shape(crLengths)[0],
cr_v=crDists[:, 0], cr_cde=crDists[:, 1], cr_cden=np.shape(crDists)[0])
##############################################################################################
def configure(self):
"""
Configures the simulator with input information and creates and empty array to which the final image will
be build on.
"""
self.readConfigs()
self.processConfigs()
self.log.info('Read in the configuration files and created an empty array')
#################################################################################################################
def MakeFlatMatrix(self, img, seed):
####
ysize, xsize=img.shape
np.random.seed(seed)
r1,r2,r3,r4 = np.random.random(4)
a1 = -0.5 + 0.2*r1
a2 = -0.5 + 0.2*r2
a3 = r3+5
a4 = r4+5
xmin,xmax,ymin,ymax = 0, xsize,0, ysize
Flty, Fltx = np.mgrid[ymin:ymax, xmin:xmax]
np.random.seed(seed)
p1,p2,bg=np.random.poisson(1000, 3)
Fltz = 1e-6*(a1 * (Fltx-p1) ** 2 + a2 * (Flty-p2) ** 2 - a3*Fltx - a4*Flty) + bg*20
FlatMat = Fltz/np.mean(Fltz)
return FlatMat
#########################################################################################################
def addLampFlux(self):
"""
Include flux from the calibration source.
"""
self.image_b += fits.getdata(self.information['flatflux'])
self.image_r += fits.getdata(self.information['flatflux'])
self.log.info('Flux from the calibration unit included (%s)' % self.information['flatflux'])
def applyflatfield(self):
"""
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.
"""
###
flat_b=self.MakeFlatMatrix(self.image_b, 100)
flat_r=self.MakeFlatMatrix(self.image_r, 200)
self.image_b *= flat_b
self.image_r *= flat_r
self.log.info('Applied flatfield to images.')
return
########################################################
###############################################################################
def addCosmicRays(self):
"""
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()
self.cr['exptime'] = self.information['exptime'] #to scale the number of cosmics with exposure time
#cosmic ray image
crImage_b = np.zeros((2048, 4096), dtype=np.float64)
crImage_r = np.zeros((3072, 6144), dtype=np.float64)
#cosmic ray instance
cosmics_b = cosmicrays.cosmicrays(self.log, crImage_b, crInfo=self.cr)
cosmics_r = cosmicrays.cosmicrays(self.log, crImage_r, crInfo=self.cr)
#add cosmic rays up to the covering fraction
#CCD_cr = cosmics.addUpToFraction(self.information['coveringFraction'], limit=None)
CCD_cr_b = cosmics_b.addUpToFraction(self.information['coveringfraction'], limit=None)
CCD_cr_r = cosmics_r.addUpToFraction(self.information['coveringfraction'], limit=None)
#paste the information
self.image_b += CCD_cr_b
self.image_r += CCD_cr_r
#count the covering factor
area_cr_b = np.count_nonzero(CCD_cr_b)
area_cr_r = np.count_nonzero(CCD_cr_r)
#self.log.info('The cosmic ray covering factor is %i pixels ' % area_cr)
self.log.info('The cosmic ray in blue channel covering factor is %i pixels ' % area_cr_b)
self.log.info('The cosmic ray in red channel covering factor is %i pixels ' % area_cr_r)
#########################################################
#########################################################################
#########################################################################
def applyDarkCurrent(self):
"""
Apply dark current. Scales the dark with the exposure time.
Additionally saves the image without noise to a FITS file.
"""
self.log.info('Added dark current to bule and red channel' )
########## blue zone 1
self.image_b[0:1024,0:2048] += self.information['exptime'] * self.information['dark1_b']
########## zone 4 #################
self.image_b[1024:2048,0:2048] += self.information['exptime'] * self.information['dark4_b']
########## zone 2 ###################
self.image_b[0:1024,2048:4096] += self.information['exptime'] * self.information['dark2_b']
########## zone 3
self.image_b[1024:2048,2048:4096]+= self.information['exptime'] * self.information['dark3_b']
########## red zone 1
self.image_r[0:1536, 0:3072] += self.information['exptime'] * self.information['dark1_r']
########## zone 4 #################
self.image_r[1536:3712,0:3072] += self.information['exptime'] * self.information['dark4_r']
########## zone 2 ###################
self.image_r[0:1536,3072:6144] += self.information['exptime'] * self.information['dark2_r']
########## zone 3
self.image_r[1536:3072,3072:6144] += self.information['exptime'] * self.information['dark3_r']
#######################################################################################################3
def applyCosmicBackground(self):
"""
Apply dark the cosmic background. Scales the background with the exposure time.
Additionally saves the image without noise to a FITS file.
"""
#add background
bcgr = self.information['exptime'] * self.information['cosmic_bkgd']
#self.image += bcgr
self.image_b += bcgr
self.image_r += bcgr
self.log.info('Added cosmic background = %f' % bcgr)
if self.cosmicRays:
#self.imagenoCR += bcgr
self.imagenoCR_b += bcgr
self.imagenoCR_r += bcgr
##########################################################################################
##############################################################################
def applyPoissonNoise(self):
"""
Add Poisson noise to the image.
"""
rounded = np.rint(self.image_b) ### round to
residual = self.image_b.copy() - rounded #ugly workaround for multiple rounding operations...
rounded[rounded < 0.0] = 0.0
np.random.seed()
self.image_b = np.random.poisson(rounded).astype(np.float64)
self.log.info('Added Poisson noise on channel blue')
self.image_b += residual
rounded = np.rint(self.image_r) ### round to
residual = self.image_r.copy() - rounded #ugly workaround for multiple rounding operations...
rounded[rounded < 0.0] = 0.0
np.random.seed()
self.image_r = np.random.poisson(rounded).astype(np.float64)
self.log.info('Added Poisson noise on channel red')
self.image_r += residual
###################################################################################################################
def applyCosmetics(self):
"""
Apply cosmetic defects described in the input file.
Warning:: This method does not work if the input file has exactly one line.
"""
cosmetics = np.loadtxt(self.information['cosmeticsfile_b'])
x = np.round(cosmetics[:, 0]).astype(int)
y = np.round(cosmetics[:, 1]).astype(int)
value = cosmetics[:, 2]
cosmetics_b=np.zeros((3712,6784))
cosmetics_r=np.zeros((3712,6784))
self.log.info('Adding cosmetic defects to blue channel:' )
for xc, yc, val in zip(x, y, value):
if 0 <= xc <= 6784 and 0 <= yc <= 3712:
#self.image[yc, xc] = val
self.image_b[yc, xc] = val
cosmetics_b[yc,xc]=val
self.log.info('x=%i, y=%i, value=%f' % (xc, yc, val))
######################################################################################################
cosmetics = np.loadtxt(self.information['cosmeticsfile_r'])
x = np.round(cosmetics[:, 0]).astype(int)
y = np.round(cosmetics[:, 1]).astype(int)
value = cosmetics[:, 2]
self.log.info('Adding cosmetic defects to red channel:' )
for xc, yc, val in zip(x, y, value):
if 0 <= xc <= 6784 and 0 <= yc <= 3712:
#self.image[yc, xc] = val
self.image_r[yc, xc] = val
cosmetics_r[yc,xc]=val
self.log.info('x=%i, y=%i, value=%f' % (xc, yc, val))
################################################################################
def applyRadiationDamage(self):
"""
Applies CDM03 radiation model to the image being constructed.
.. seealso:: Class :`CDM03`
"""
self.log.debug('Starting to apply radiation damage model...')
#at this point we can give fake data...
cti = CTI.CDM03bidir(self.information, [], log=self.log)
#here we need the right input data
self.image_b = cti.applyRadiationDamage(self.image_b.copy().transpose(), iquadrant=self.information['quadrant']).transpose()
self.log.info('Radiation damage added.')
self.log.debug('Starting to apply radiation damage model...')
#at this point we can give fake data...
cti = CTI.CDM03bidir(self.information, [], log=self.log)
#here we need the right input data
self.image_r = cti.applyRadiationDamage(self.image_r.copy().transpose(), iquadrant=self.information['quadrant']).transpose()
self.log.info('Radiation damage added.')
##################################################################################
def applyNonlinearity(self):
"""
Applies a CCD273 non-linearity model to the image being constructed.
"""
self.log.debug('Starting to apply non-linearity model...')
self.image_b = IFSinstrumentModel.CCDnonLinearityModel(self.image_b.copy())
self.log.info('Non-linearity effects included.')
self.log.debug('Starting to apply non-linearity model...')
self.image_r = IFSinstrumentModel.CCDnonLinearityModel(self.image_r.copy())
self.log.info('Non-linearity effects included.')
#####################################################################################
def applyReadoutNoise(self):
"""
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.
"""
self.log.info('readnoise added in blue channel')
########## blue zone 1
np.random.seed()
self.image_b[0:1856,0:3392] += np.random.normal(loc=0.0, scale=self.information['rn1_b'], size=(1856,3392))
########## zone 4 #################
np.random.seed()
self.image_b[1856:3712,0:3392] += np.random.normal(loc=0.0, scale=self.information['rn4_b'], size=(1856,3392))
########## zone 2 ###################
np.random.seed()
self.image_b[0:1856,3392:6784] += np.random.normal(loc=0.0, scale=self.information['rn2_b'], size=(1856,3392))
########## zone 3
np.random.seed()
self.image_b[1856:3712,3392:6784]+= np.random.normal(loc=0.0, scale=self.information['rn3_b'], size=(1856,3392))
############################################################################
self.log.info('readnoise added in blue channel')
########## red zone 1
np.random.seed()
self.image_r[0:1856, 0:3392] += np.random.normal(loc=0.0, scale=self.information['rn1_r'], size=(1856,3392))
########## zone 4 #################
np.random.seed()
self.image_r[1856:3712,0:3392] += np.random.normal(loc=0.0, scale=self.information['rn4_r'], size=(1856,3392))
########## zone 2 ###################
np.random.seed()
self.image_r[0:1856,3392:6784] += np.random.normal(loc=0.0, scale=self.information['rn2_r'], size=(1856,3392))
########## zone 3
np.random.seed()
self.image_r[1856:3712,3392:6784]+= np.random.normal(loc=0.0, scale=self.information['rn3_r'], size=(1856,3392))
##########################################################################################
def electrons2ADU(self):
"""
Convert from electrons to ADUs using the value read from the configuration file.
"""
###############################################################
self.log.info('Converting from electrons to ADUs using a factor of gain' )
########## blue zone 1
self.image_b[0:1856,0:3392] /= self.information['gain1_b']
########## zone 4 #################
self.image_b[1856:3712,0:3392] /= self.information['gain4_b']
########## zone 2 ###################
self.image_b[0:1856,3392:6784] /= self.information['gain2_b']
########## zone 3
self.image_b[1856:3712,3392:6784]/= self.information['gain3_b']
############################################################################
########## red zone 1
self.image_r[0:1856, 0:3392] /= self.information['gain1_r']
########## zone 4 #################
self.image_r[1856:3712,0:3392] /= self.information['gain4_r']
########## zone 2 ###################
self.image_r[0:1856,3392:6784] /= self.information['gain2_r']
########## zone 3
self.image_r[1856:3712,3392:6784] /= self.information['gain3_r']
##########################################################################3
####################################################################################
def applyBias(self):
"""
Adds a bias level to the image being constructed.
The value of bias is read from the configure file and stored
in the information dictionary (key bias).
"""
########## blue zone 1
self.image_b[0:1856,0:3392] += self.information['bias1_b']
########## zone 4 #################
self.image_b[1856:3712,0:3392] += self.information['bias4_b']
########## zone 2 ###################
self.image_b[0:1856,3392:6784] += self.information['bias2_b']
########## zone 3
self.image_b[1856:3712,3392:6784] += self.information['bias3_b']
self.log.info('Bias counts were added to the blue image' )
############################################################################
########## red zone 1
self.image_r[0:1856, 0:3392] += self.information['bias1_r']
########## zone 4 #################
self.image_r[1856:3712,0:3392] += self.information['bias4_r']
########## zone 2 ###################
self.image_r[0:1856,3392:6784] += self.information['bias2_r']
########## zone 3
self.image_r[1856:3712,3392:6784] += self.information['bias3_r']
##########################################################################
self.log.info('Bias counts were added to the red image' )
###############################################################################
###############################################################################
def applyBleeding_yan(self):
"""
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
serial direction.
:return: None
"""
if self.image_b.max()>self.information['fullwellcapacity']:
self.log.info('Applying column bleeding to blue CCD image...')
#loop over each column, as bleeding is modelled column-wise
for i, column in enumerate(self.image_b.T):
sum = 0.
for j, value in enumerate(column):
#first round - from bottom to top (need to half the bleeding)
overload = value - self.information['fullwellcapacity']
if overload > 0.:
overload /= 2.
#self.image[j, i] -= overload
self.image_b[j, i] -= overload
sum += overload
elif sum > 0.:
if -overload > sum:
overload = -sum
self.image_b[j, i] -= overload
sum += overload
for i, column in enumerate(self.image_b.T):
sum = 0.
for j, value in enumerate(column[::-1]):
#second round - from top to bottom (bleeding was half'd already, so now full)
overload = value - self.information['fullwellcapacity']
if overload > 0.:
#self.image[-j-1, i] -= overload
self.image_b[-j-1, i] -= overload
sum += overload
elif sum > 0.:
if -overload > sum:
overload = -sum
#self.image[-j-1, i] -= overload
self.image_b[-j-1, i] -= overload
sum += overload
print('Applying column bleeding to blue image finished.......')
######################################################################
if self.image_r.max()>self.information['fullwellcapacity']:
self.log.info('Applying column bleeding to red CCD image...')
for i, column in enumerate(self.image_r.T):
sum = 0.
for j, value in enumerate(column):
#first round - from bottom to top (need to half the bleeding)
overload = value - self.information['fullwellcapacity']
if overload > 0.:
overload /= 2.
#self.image[j, i] -= overload
self.image_r[j, i] -= overload
sum += overload
elif sum > 0.:
if -overload > sum:
overload = -sum
#self.image[j, i] -= overload
self.image_r[j, i] -= overload
sum += overload
for i, column in enumerate(self.image_r.T):
sum = 0.
for j, value in enumerate(column[::-1]):
#second round - from top to bottom (bleeding was half'd already, so now full)
overload = value - self.information['fullwellcapacity']
if overload > 0.:
#self.image[-j-1, i] -= overload
self.image_r[-j-1, i] -= overload
sum += overload
elif sum > 0.:
if -overload > sum:
overload = -sum
#self.image[-j-1, i] -= overload
self.image_r[-j-1, i] -= overload
sum += overload
print('Applying column bleeding to red image finished.......')
############################################################################
############################################################################
def discretise(self, max=2**16-1):
"""
Converts a floating point image array (self.image) to an integer array with max values
defined by the argument max.
:param max: maximum value the the integer array may contain [default 65k]
:type max: float
:return: None
"""
#avoid negative numbers in case bias level was not added
self.image_b[self.image_b < 0.0] = 0.
#cut of the values larger than max
self.image_b[self.image_b > max] = max
self.image_b = np.rint(self.image_b).astype(int)
self.log.info('Maximum and total values of the image are %i and %i, respectively' % (np.max(self.image_b),
np.sum(self.image_b)))
#avoid negative numbers in case bias level was not added
self.image_r[self.image_r < 0.0] = 0.
#cut of the values larger than max
self.image_r[self.image_r > max] = max
self.image_r = np.rint(self.image_r).astype(int)
self.log.info('Maximum and total values of the image are %i and %i, respectively' % (np.max(self.image_r),
np.sum(self.image_r)))
##################################################################################################
def applyImageShift(self):
np.random.seed()
ud= np.random.random() # Choose a random rotation
dx= 2* (ud-0.5) * self.information['shiftmax']
np.random.seed()
ud= np.random.random() # Choose a random rotation
dy= 2* (ud-0.5) * self.information['shiftmax']
self.image_b= ndimage.shift(self.image_b.copy(), [dy+self.information['shift_b_y'] , dx+self.information['shift_b_x']], order=0, mode='nearest')
self.image_r= ndimage.shift(self.image_r.copy(), [dy+self.information['shift_r_y'] , dx+self.information['shift_r_x']], order=0, mode='nearest')
self.log.info('Applied image shifting to g r i channels.')
self.information['ra'] = dx*self.information['pixel_size']
self.information['dec'] = dy*self.information['pixel_size']
######################################################################################################33
def applyImageRotate(self ):
np.random.seed()
ud = np.random.random() # Choose a random rotation
angle = 2* (ud-0.5) * self.information['tel_rotmax']
inputimg=self.image_b.copy()
rotimg=ndimage.rotate(inputimg,angle+self.information['rotate_b'],order=1,reshape=False) # here we choose reshape=False, the rotated image will
self.image_b=rotimg
inputimg=self.image_r.copy()
rotimg=ndimage.rotate(inputimg,angle+self.information['rotate_r'],order=1,reshape=False) # here we choose reshape=False, the rotated image will
self.image_r=rotimg
self.information['Tel_rot']=angle
self.log.info('Applied telescope rotation with angle (in degree)= %f.', angle)
###############################################################################
def CCDreadout(self):
imgb=self.image_b.copy()
temp=np.zeros((3712,6784))
########## zone 1
x1=0
x2=x1+1024
y1=0
y2=y1+2048
temp[x1:x2,y1:y2]=imgb[0:1024,0:2048]
########## zone 4 #################
x1=2688
x2=x1+1024
y1=0
y2=y1+2048
temp[x1:x2,y1:y2]=imgb[1024:2048,0:2048]
########## zone 2 ###################
x1=0
x2=x1+1024
y1=6784-2048
y2=y1+2048
temp[x1:x2,y1:y2]=imgb[0:1024,2048:4096]
########## zone 3
x1=2688
x2=x1+1024
y1=6784-2048
y2=y1+2048
temp[x1:x2,y1:y2]=imgb[1024:2048,2048:4096]
self.image_b=temp
##############################################################################
imgr=self.image_r.copy()
temp=np.zeros((3712,6784))
########## zone 1
x1=0
x2=x1+1536
y1=0
y2=y1+3072
temp[x1:x2,y1:y2]=imgr[0:1536, 0:3072]
########## zone 4 #################
x1=2176
x2=x1+1536
y1=0
y2=y1+3072
temp[x1:x2,y1:y2]=imgr[1536:3712,0:3072]
########## zone 2 ###################
x1=0
x2=x1+1536
y1=6784-3072
y2=y1+3072
temp[x1:x2,y1:y2]=imgr[0:1536,3072:6144]
########## zone 3
x1=2176
x2=x1+1536
y1=6784-3072
y2=y1+3072
temp[x1:x2,y1:y2]=imgr[1536:3072,3072:6144]
self.image_r=temp
return
##############################################################################
def writeOutputs(self):
"""
Writes out a FITS file using PyFITS and converts the image array to 16bit unsigned integer as
appropriate for VIS.
Updates header with the input values and flags used during simulation.
"""
## Readout information
self.source='sci'
now=datetime.now()
data_time=now.strftime("%Y-%m-%d %H:%M:%S")
exp_endtime=now.strftime("%Y%m%d%H%M%S")
start=now-timedelta(seconds=self.information['exptime'])
exp_starttime=start.strftime("%Y%m%d%H%M%S")
#write the actual file
obsid=300000000+1
data_time=self.dt.strftime("%Y-%m-%d %H:%M:%S")
exp_starttime=self.dt.strftime("%Y%m%d%H%M%S")
### exposure end time is t2 ;
t2=self.dt+timedelta(seconds=self.information['exptime'])
exp_endtime=t2.strftime("%Y%m%d%H%M%S")
t3=self.dt+timedelta(seconds=self.information['exptime'])+timedelta(seconds=self.information['readouttime'])
filename_b='CSST_IFS_B_'+self.source+'_'+exp_starttime+'_'+exp_endtime+'_'+str(obsid)+'_X_L0_VER_'+self.information['img_ver']+'.fits'
file_b=self.result_path+'/sky_Data/'+filename_b
filename_r='CSST_IFS_R_'+self.source+'_'+exp_starttime+'_'+exp_endtime+'_'+str(obsid)+'_X_L0_VER_'+self.information['img_ver']+'.fits'
file_r=self.result_path+'/sky_Data/'+filename_r
#create a new FITS file, using HDUList instance
ofd_b = fits.PrimaryHDU()
ofd_b.header['GROUPS']=( bool(False), 'always F')
ofd_b.header['DATE'] =( data_time, 'date this file was written' )
ofd_b.header['FILENAME']=(filename_b, ' file name C48 ')
ofd_b.header['OBSTYPE'] =( self.source, 'observation type raw,flt, mask, bias, dark, sci')
ofd_b.header['TELESCOP']=('CSST', 'always CSST')
ofd_b.header['INSTRUME']=( 'IFS', ' ')
ofd_b.header['RADECSYS']=('ICRS', ' always ICRS ')
ofd_b.header['EQUINOX'] =( float(2000.0), 'always 2000.0')
ofd_b.header['FITSCREA']=( '4.2.1', 'FITS create software version')
######### Object information #############
ofd_b.header['OBJECT']=( self.information['name_obj'], 'object name')
ofd_b.header['TARGET']=( (self.information['target']), 'target name, hhmmss+ddmmss')
ofd_b.header['OBJ_RA'] =(np.float64(self.information['ra_obj']) , 'RA of the object in deg')
ofd_b.header['OBJ_DEC']=(np.float64(self.information['dec_obj']) , 'DEC of the object in deg')
ofd_b.header['RA_PNT0']=( np.float64(self.information['ra_pnt0']) , 'RA of the pointing (degrees) at EXPSTART')
ofd_b.header['DEC_PNT0']=(np.float64(self.information['dec_pnt0']) , 'DEC of the pointing (degrees) at EXPSTART')
##############
ofd_b.header['OBSID'] =(str(obsid) , 'observation ID, 3+8bit')
######## Telescope information ###############
# ofd_b.header['COMMENT'] ='=========================================================='
# ofd_b.header['COMMENT'] ='Telescope information'
# ofd_b.header['COMMENT'] ='=========================================================='
ofd_b.header['REFFRAME']=('CSSTGSC-1.0' , 'guide star catalog version')
ofd_b.header['DATE-OBS']=(data_time , 'date of the observation (yyyy-mm-dd hh:mm:ss)')
ofd_b.header['EXPSTART']=(np.float64(time2jd(self.dt)), 'exposure start time')
ofd_b.header['SUNANGL0']=(np.float32(0.0) , 'angle between sun and optical axis at EXPSTART')
ofd_b.header['MOONANG0']=(np.float32(0.0) , 'angle between moon and optical axis at EXPSTART')
ofd_b.header['POS_ANG0']=(np.float64(0.0), 'angle between optical axis and the North Pole at EXPSTART in arcsec')
ofd_b.header['TEL_ALT0']=(np.float64(0.0), 'angle between optical axis and the ground- piston at EXPSTART in deg')
ofd_b.header['HOODSTA0']=(np.float32(0.0) , 'lens hood altitude at EXPSTART')
ofd_b.header['HOODANG0']=(np.float32(0.0), 'lens hood azimuth at EXPSTART')
ofd_b.header['POSI0_X'] =(np.float64(self.information['POSI0_X']) , 'the orbital position of CSST in X direction at EXPSTART')
ofd_b.header['POSI0_Y'] =(np.float64(self.information['POSI0_Y']) , 'the orbital position of CSST in Y direction at EXPSTART')
ofd_b.header['POSI0_Z'] =(np.float64(self.information['POSI0_Z']) , 'the orbital position of CSST in Z direction at EXPSTART')
ofd_b.header['VELO0_X']=( np.float64(self.information['VELO0_X']) , 'the orbital velocity of CSST in X direction at EXPSTART')
ofd_b.header['VELO0_Y']=( np.float64(self.information['VELO0_Y']) , 'the orbital velocity of CSST in Y direction at EXPSTART')
ofd_b.header['VELO0_Z']=( np.float64(self.information['VELO0_Z']) , 'the orbital velocity of CSST in Z direction at EXPSTART')
ofd_b.header['Euler0_1']=( np.float64(0.0), 'Euler angle 1 at EXPSTART')
ofd_b.header['Euler0_2']=( np.float64(0.0), 'Euler angle 2 at EXPSTART')
ofd_b.header['Euler0_3']=( np.float64(0.0), 'Euler angle 3 at EXPSTART')
ofd_b.header['EXPEND'] =(np.float64(time2jd(t2)) , 'exposure end time')
ofd_b.header['SUNANGL1']=(np.float32(0.0), 'angle between sun and optical axis at EXPEND')
ofd_b.header['MOONANG1']=(np.float32(0.0) , 'angle between moon and optical axis at EXPEND ')
ofd_b.header['POS_ANG1']=(np.float64(0.0) , 'angle between optical axis and the North Pole at EXPEND in arcsec')
ofd_b.header['TEL_ALT1']=(np.float64(0.0) , 'angle between optical axis and the ground- piston at EXPEND in deg ')
ofd_b.header['HOODSTA1']=(np.float32(0.0), 'lens hood altitude at EXPEND ')
ofd_b.header['HOODANG1']=(np.float32(0.0), 'lens hood azimuth at EXPEND ')
ofd_b.header['POSI1_X'] =(np.float64(self.information['POSI1_X']) , 'the orbital position of CSST in X direction at EXPEND')
ofd_b.header['POSI1_Y'] =(np.float64(self.information['POSI1_Y']) , 'the orbital position of CSST in Y direction at EXPEND')
ofd_b.header['POSI1_Z'] =(np.float64(self.information['POSI1_Z']) , 'the orbital position of CSST in Z direction at EXPEND')
ofd_b.header['VELO1_X']=(np.float64(self.information['VELO1_X']) , 'the orbital velocity of CSST in X direction at EXPEND')
ofd_b.header['VELO1_Y']=(np.float64(self.information['VELO1_Y']) , 'the orbital velocity of CSST in Y direction at EXPEND')
ofd_b.header['VELO1_Z']=(np.float64(self.information['VELO1_Z']) , 'the orbital velocity of CSST in Z direction at EXPEND')
ofd_b.header['Euler1_1']=( np.float64(0.0), 'Euler angle 1 at EXPEND')
ofd_b.header['Euler1_2']=( np.float64(0.0), 'Euler angle 2 at EXPEND')
ofd_b.header['Euler1_3']=( np.float64(0.0), 'Euler angle 3 at EXPEND')
ofd_b.header['RA_PNT1']=(np.float64(ofd_b.header['RA_PNT0']), 'RA of the pointing (degrees) at EXPEND in deg')
ofd_b.header['DEC_PNT1']=(np.float64(ofd_b.header['DEC_PNT0']), 'DEC of the pointing (degrees) at EXPEND in deg')
ofd_b.header['EXPTIME']=(self.information['exptime'], 'exposure duration')
ofd_b.header['EPOCH'] =(np.float32(0.0), 'coordinate epoch')
ofd_b.header['CHECKSUM']=( 0 , 'hdu-checksum')
########## finish header for 0 layer
#############################################################################3
##### header
b1= self.image_b[1856:3712,0:3392] #b4
b2= self.image_b[1856:3712,3392:6784] #b3
b3= self.image_b[0:1856,0:3392] #b1
b4= self.image_b[0:1856,3392:6784] #b2
####### do Flip the b2 b2 and b4 array in the up/down or left/right direction.
b2=np.fliplr(b2) ## left to right
b3=np.flipud(b3) ## down to up
b4=np.fliplr(b4) ## left to right and down to up
b4=np.flipud(b4)
bb=np.hstack((b1,b2,b3,b4))
#new image HDU, blue channel, layer 1
hdu_b =fits.ImageHDU(data=np.uint16(bb))
######### instrument information ######
#####
hdu_b.header['PMIRRPOS']=(bool(False), 'FSM pointing,T: to MCI, F: not to MCI')
if self.source =='sci':
hdu_b.header['CMIRRPOS']=(bool(False), 'position of calibration switch mirror,T: for calibration, F: not')
else:
hdu_b.header['CMIRRPOS']=(bool(True), 'position of calibration switch mirror,T: for calibration, F: not')
if self.source=='flat':
hdu_b.header['FLAMP'] =(int(1), 'status of flat lamp,0: off, 1: , 2: ')
else:
hdu_b.header['FLAMP'] =(int(0), 'status of flat lamp,0: off, 1: , 2: ')
if self.source=='lamp':
hdu_b.header['ALAMP'] =(int(1),'status of atomic emission line lamp,0: off, 1: , 2: ')
else:
hdu_b.header['ALAMP'] =(int(0),'status of atomic emission line lamp,0: off, 1: , 2: ')
#############
hdu_b.header['IFSMODE'] =(int(0), 'IFS working mode')
hdu_b.header['IFSTEMP'] =(float(0.0), 'IFS components temperature in degC')
hdu_b.header['IFSSTAT'] =(int(0), 'IFS components status parameter')
##############################################################################
################### detector information#############################
# hdu_b.header['COMMENT'] ='=========================================================='
# hdu_b.header['COMMENT'] ='Detector information'
# hdu_b.header['COMMENT'] ='=========================================================='
hdu_b.header['CAMERA'] =('Blue','camera of IFS')
hdu_b.header['DETNAM'] =('CCD231-c4','detector name')
hdu_b.header['DETSIZE'] =('', 'detector size')
hdu_b.header['DATASEC'] =('', 'data section')
hdu_b.header['PIXSCAL1']=(1856, 'pixel scale for axis 1')
hdu_b.header['PIXSCAL2']=(3392*4, 'pixel scale for axis 2')
hdu_b.header['PIXSIZE1']=(15, 'pixel size in um')
hdu_b.header['PIXSIZE2']=(15, 'pixel size in um')
hdu_b.header['NCHAN'] =(4, 'number of readout channels')
hdu_b.header['NCHAN1'] =(2, 'number of horizontal channels')
hdu_b.header['NCHAN2'] =(2, 'number of verticalchannels')
hdu_b.header['PSCAN1'] =(0, 'horizontal prescan width, per readout channel')
hdu_b.header['PSCAN2'] =(0, 'vertical prescan width, per readout channel')
hdu_b.header['OSCAN1'] =(0,' horizontal overscan width, per readout channel')
hdu_b.header['OSCAN2'] =(0, 'vertical overscan width, per readout channel')
## Readout information
# hdu_b.header['COMMENT'] ='============================================================='
# hdu_b.header['COMMENT'] ='Readout information'
# hdu_b.header['COMMENT'] ='============================================================='
hdu_b.header['READT0'] =(np.float64(time2jd(t2)),'read start time (UTC)')
hdu_b.header['READT1'] =(np.float64(time2jd(t3)), 'read end time (UTC)')
hdu_b.header['DETTEMP0']=(np.float32(0.0), 'detector temperature at READT0')
hdu_b.header['DETTEMP1']=(np.float32(0.0), 'detector temperature at READT1')
hdu_b.header['BIN_X'] =(0, 'bin number in X (wavelength)')
hdu_b.header['BIN_Y'] =(0, 'bin number in Y (spatial)')
hdu_b.header['GAIN1'] =(self.information['gain4_b'],'CCD gain (channel 1)')
hdu_b.header['GAIN2'] =(self.information['gain3_b'],'CCD gain (channel 2)')
hdu_b.header['GAIN3'] =(self.information['gain1_b'],'CCD gain (channel 3)')
hdu_b.header['GAIN4'] =(self.information['gain2_b'],'CCD gain (channel 4)')
hdu_b.header['DARK1'] =(self.information['dark4_b'],'CCD dark (channel 1)')
hdu_b.header['DARK2'] =(self.information['dark3_b'],'CCD dark (channel 2)')
hdu_b.header['DARK3'] =(self.information['dark1_b'],'CCD dark (channel 3)')
hdu_b.header['DARK4'] =(self.information['dark2_b'],'CCD dark (channel 4)')
hdu_b.header['RDNOIS1'] =(self.information['rn4_b'],'read noise (channel 1')
hdu_b.header['RDNOIS2'] =(self.information['rn3_b'],'read noise (channel 2')
hdu_b.header['RDNOIS3'] =(self.information['rn1_b'],'read noise (channel 3')
hdu_b.header['RDNOIS4'] =(self.information['rn2_b'],'read noise (channel 4')
hdu_b.header['DETBIA1'] =(self.information['bias4_b'],'amplifier bias voltage (channel1)')
hdu_b.header['DETBIA2'] =(self.information['bias3_b'],'amplifier bias voltage (channel2)')
hdu_b.header['DETBIA3'] =(self.information['bias1_b'],'amplifier bias voltage (channel3)')
hdu_b.header['DETBIA4'] =(self.information['bias2_b'],'amplifier bias voltage (channel4)')
hdu_b.header['RDSPEED'] =(100,'read speed (in MHz)')
hdu_b.header['EXPTIME'] =(self.information['exptime'],'exposure time in seconds')
hdu_b.header['Img_Ver'] =(self.information['img_ver'], 'IFS CCD image Version')
hdu_b.header['sky_obj'] =(self.skyfilepath,'input sky fits filepath')
##########################################################
#################### red camera ######################
#create a new FITS file, using HDUList instance
ofd_r = fits.PrimaryHDU()
ofd_r.header['GROUPS']=( bool(False), 'always F')
ofd_r.header['DATE'] =( data_time, 'date this file was written' )
ofd_r.header['FILENAME']=(filename_r, ' file name C48 ')
ofd_r.header['OBSTYPE'] =( self.source, 'observation type raw,flt, mask, bias, dark, sci')
ofd_r.header['TELESCOP']=('CSST' , 'always CSST')
ofd_r.header['INSTRUME']=( 'IFS', ' ')
ofd_r.header['RADECSYS']=('ICRS', ' always ICRS ')
ofd_r.header['EQUINOX'] =( float(2000.0), 'always 2000.0')
ofd_r.header['FITSCREA']=( '4.2.1' , 'FITS create software version')
######### Object information #############
######### Object information #############
# ofd_r.header['COMMENT']='======================================================================='
# ofd_r.header['COMMENT']='Object information'
# ofd_r.header['COMMENT']='======================================================================='
ofd_r.header['OBJECT']=( self.information['name_obj'], 'object name')
ofd_r.header['TARGET']=( (self.information['target']), 'target name, hhmmss+ddmmss')
ofd_r.header['OBJ_RA'] =(np.float64(self.information['ra_obj']) , 'RA of the object in deg')
ofd_r.header['OBJ_DEC']=(np.float64(self.information['dec_obj']) , 'DEC of the object in deg')
ofd_r.header['RA_PNT0']=( np.float64(self.information['ra_pnt0']) , 'RA of the pointing (degrees) at EXPSTART')
ofd_r.header['DEC_PNT0']=(np.float64(self.information['dec_pnt0']) , 'DEC of the pointing (degrees) at EXPSTART')
ofd_r.header['OBSID'] =(str(obsid) , 'observation ID, 3+8bit ')
######## Telescope information ###############
# ofd_r.header['COMMENT']='======================================================================='
# ofd_r.header['COMMENT']='Telescope information'
# ofd_r.header['COMMENT']='======================================================================='
ofd_r.header['REFFRAME']=('CSSTGSC-1.0' , 'guide star catalog version ')
ofd_r.header['DATE-OBS']=(data_time , 'date of the observation (yyyy-mm-dd hh:mm:ss)')
ofd_r.header['EXPSTART']=(np.float64(exp_starttime), 'exposure start time ')
ofd_r.header['SUNANGL0']=(np.float32(0.0) , 'angle between sun and optical axis at EXPSTART ')
ofd_r.header['MOONANG0']=(np.float32(0.0) , 'angle between moon and optical axis at EXPSTART ')
ofd_r.header['POS_ANG0']=(np.float64(0.0), 'angle between optical axis and the North Pole at EXPSTART in arcsec ')
ofd_r.header['TEL_ALT0']=(np.float64(0.0), 'angle between optical axis and the ground- piston at EXPSTART in deg ')
ofd_r.header['HOODSTA0']=(np.float32(0.0) , 'lens hood altitude at EXPSTART ')
ofd_r.header['HOODANG0']=(np.float32(0.0), 'lens hood azimuth at EXPSTART ')
ofd_r.header['POSI0_X'] =(np.float64(self.information['POSI0_X']) , 'the orbital position of CSST in X direction at EXPSTART')
ofd_r.header['POSI0_Y'] =(np.float64(self.information['POSI0_Y']) , 'the orbital position of CSST in Y direction at EXPSTART')
ofd_r.header['POSI0_Z'] =(np.float64(self.information['POSI0_Z']) , 'the orbital position of CSST in Z direction at EXPSTART')
ofd_r.header['VELO0_X']=( np.float64(self.information['VELO0_X']) , 'the orbital velocity of CSST in X direction at EXPSTART')
ofd_r.header['VELO0_Y']=( np.float64(self.information['VELO0_Y']) , 'the orbital velocity of CSST in Y direction at EXPSTART')
ofd_r.header['VELO0_Z']=( np.float64(self.information['VELO0_Z']) , 'the orbital velocity of CSST in Z direction at EXPSTART')
ofd_r.header['Euler0_1']=( np.float64(0.0), 'Euler angle 1 at EXPSTART')
ofd_r.header['Euler0_2']=( np.float64(0.0), 'Euler angle 2 at EXPSTART')
ofd_r.header['Euler0_3']=( np.float64(0.0), 'Euler angle 3 at EXPSTART')
ofd_r.header['EXPEND'] =(np.float64(exp_endtime) , 'exposure end time')
ofd_r.header['SUNANGL1']=(np.float32(0.0), 'angle between sun and optical axis at EXPEND')
ofd_r.header['MOONANG1']=(np.float32(0.0) , 'angle between moon and optical axis at EXPEND ')
ofd_r.header['POS_ANG1']=(np.float64(0.0) , 'angle between optical axis and the North Pole at EXPEND in arcsec')
ofd_r.header['TEL_ALT1']=(np.float64(0.0) , 'angle between optical axis and the ground- piston at EXPEND in deg ')
ofd_r.header['HOODSTA1']=(np.float32(0.0), 'lens hood altitude at EXPEND ')
ofd_r.header['HOODANG1']=(np.float32(0.0), 'lens hood azimuth at EXPEND ')
ofd_r.header['POSI1_X'] =(np.float64(self.information['POSI1_X']) , 'the orbital position of CSST in X direction at EXPEND')
ofd_r.header['POSI1_Y'] =(np.float64(self.information['POSI1_Y']) , 'the orbital position of CSST in Y direction at EXPEND')
ofd_r.header['POSI1_Z'] =(np.float64(self.information['POSI1_Z']) , 'the orbital position of CSST in Z direction at EXPEND')
ofd_r.header['VELO1_X']=(np.float64(self.information['VELO1_X']) , 'the orbital velocity of CSST in X direction at EXPEND')
ofd_r.header['VELO1_Y']=(np.float64(self.information['VELO1_Y']) , 'the orbital velocity of CSST in Y direction at EXPEND')
ofd_r.header['VELO1_Z']=(np.float64(self.information['VELO1_Z']) , 'the orbital velocity of CSST in Z direction at EXPEND')
ofd_r.header['Euler1_1']=( np.float64(0.0), 'Euler angle 1 at EXPEND')
ofd_r.header['Euler1_2']=( np.float64(0.0), 'Euler angle 2 at EXPEND')
ofd_r.header['Euler1_3']=( np.float64(0.0), 'Euler angle 3 at EXPEND')
ofd_r.header['RA_PNT1']=(np.float64(ofd_r.header['RA_PNT0']), 'RA of the pointing (degrees) at EXPEND in deg')
ofd_r.header['DEC_PNT1']=(np.float64(ofd_r.header['DEC_PNT0']), 'DEC of the pointing (degrees) at EXPEND in deg')
ofd_r.header['EXPTIME']=(self.information['exptime'], 'exposure duration')
ofd_r.header['EPOCH'] =(np.float32(0.0), 'coordinate epoch')
ofd_r.header['CHECKSUM']=( 0 , 'hdu-checksum')
### finish 0 layer header
########## finish header for 0 layer
# ########## blue zone 1--to--3
# self.image_r[0:1856,0:3392] += self.information['bias1_r']
# ########## zone 4 --to---1 #################
# self.image_r[1856:3712,0:3392] += self.information['bias4_r']
# ########## zone 2 ----to----4 ###################
# self.image_r[0:1856,3392:6784] += self.information['bias2_r']
# ########## zone 3 ---to------2
# self.image_r[1856:3712,3392:6784] += self.information['bias3_r']
#############################################################################3
##### header
b1= self.image_r[1856:3712,0:3392]
b2= self.image_r[1856:3712,3392:6784]
b3= self.image_r[0:1856,0:3392]
b4= self.image_r[0:1856,3392:6784]
####### do Flip the b2 b2 and b4 array in the up/down or left/right direction.
b2=np.fliplr(b2) ## left to right
b3=np.flipud(b3) ## down to up
b4=np.fliplr(b4) ## left to right and down to up
b4=np.flipud(b4)
rr=np.hstack((b1,b2,b3,b4))
#new image HDU, blue channel, layer 1
hdu_r =fits.ImageHDU(data=np.uint16(rr))
#########################################
######### instrument information ######
hdu_r.header['PMIRRPOS']=(bool(False), 'FSM pointing,T: to MCI, F: not to MCI')
hdu_r.header['CMIRRPOS']=(bool(False), 'position of calibration switch mirror,T: for calibration, F: not')
hdu_r.header['FLAMP'] =(int(0), 'status of flat lamp,0: off, 1: , 2: ')
hdu_r.header['ALAMP'] =(int(0),'status of atomic emission line lamp,0: off, 1: , 2: ')
hdu_r.header['IFSMODE'] =(int(0), 'IFS working mode')
hdu_r.header['IFSTEMP'] =(float(0.0), 'IFS components temperature in degC')
hdu_r.header['IFSSTAT'] =(int(0), 'IFS components status parameter')
################### detector information#############################
# hdu_r.header['COMMENT']='======================================================================='
# hdu_r.header['COMMENT']='Detector information'
# hdu_r.header['COMMENT']='======================================================================='
hdu_r.header['CAMERA'] =('Red','camera of IFS')
hdu_r.header['DETNAM'] =('CCD231-c4','detector name')
hdu_r.header['DETSIZE'] =('', 'detector size')
hdu_r.header['DATASEC'] =('', 'data section')
hdu_r.header['PIXSCAL1']=(1856, 'pixel scale for axis 1')
hdu_r.header['PIXSCAL2']=(3392*4, 'pixel scale for axis 2')
hdu_r.header['PIXSIZE1']=(15, 'pixel size in um')
hdu_r.header['PIXSIZE2']=(15, 'pixel size in um')
hdu_r.header['NCHAN'] =(4, 'number of readout channels')
hdu_r.header['NCHAN1'] =(2, 'number of horizontal channels')
hdu_r.header['NCHAN2'] =(2, 'number of verticalchannels')
hdu_r.header['PSCAN1'] =(0, 'horizontal prescan width, per readout channel')
hdu_r.header['PSCAN2'] =(0, 'vertical prescan width, per readout channel')
hdu_r.header['OSCAN1'] =(0,' horizontal overscan width, per readout channel')
hdu_r.header['OSCAN2'] =(0, 'vertical overscan width, per readout channel')
#####################################################################################################
## Readout information
# hdu_r.header['COMMENT']='======================================================================='
# hdu_r.header['COMMENT']='Readout information'
# hdu_r.header['COMMENT']='======================================================================='
hdu_r.header['READT0'] =(np.float64(time2jd(t2)),'read start time (UTC)')
hdu_r.header['READT1'] =(np.float64(time2jd(t3)), 'read end time (UTC)')
hdu_r.header['DETTEMP0']=(np.float32(0.0), 'detector temperature at READT0')
hdu_r.header['DETTEMP1']=(np.float32(0.0), 'detector temperature at READT1')
hdu_r.header['BIN_X'] =(0, 'bin number in X (wavelength)')
hdu_r.header['BIN_Y'] =(0, 'bin number in Y (spatial)')
hdu_r.header['GAIN1'] =(self.information['gain4_r'],'CCD gain (channel 1)')
hdu_r.header['GAIN2'] =(self.information['gain3_r'],'CCD gain (channel 2)')
hdu_r.header['GAIN3'] =(self.information['gain1_r'],'CCD gain (channel 3)')
hdu_r.header['GAIN4'] =(self.information['gain2_r'],'CCD gain (channel 4)')
hdu_r.header['DARK1'] =(self.information['dark4_r'],'CCD dark (channel 1)')
hdu_r.header['DARK2'] =(self.information['dark3_r'],'CCD dark (channel 2)')
hdu_r.header['DARK3'] =(self.information['dark1_r'],'CCD dark (channel 3)')
hdu_r.header['DARK4'] =(self.information['dark2_r'],'CCD dark (channel 4)')
hdu_r.header['RDNOIS1'] =(self.information['rn4_r'],'read noise (channel 1')
hdu_r.header['RDNOIS2'] =(self.information['rn3_r'],'read noise (channel 2')
hdu_r.header['RDNOIS3'] =(self.information['rn1_r'],'read noise (channel 3')
hdu_r.header['RDNOIS4'] =(self.information['rn2_r'],'read noise (channel 4')
hdu_r.header['DETBIA1'] =(self.information['bias4_r'],'amplifier bias voltage (channel1)')
hdu_r.header['DETBIA2'] =(self.information['bias3_r'],'amplifier bias voltage (channel2)')
hdu_r.header['DETBIA3'] =(self.information['bias1_r'],'amplifier bias voltage (channel3)')
hdu_r.header['DETBIA4'] =(self.information['bias2_r'],'amplifier bias voltage (channel4)')
hdu_r.header['RDSPEED'] =(100,'read speed (in MHz)')
hdu_r.header['EXPTIME'] =(self.information['exptime'],'exposure time in seconds')
hdu_r.header['Img_Ver'] =(self.information['img_ver'], 'IFS CCD image Version')
hdu_r.header['sky_obj'] =(self.skyfilepath,'input sky fits filename')
hdulist_b=fits.HDUList([ofd_b, hdu_b])
hdulist_b.writeto(file_b, overwrite=True)
#print('IFS_b.fits is created ')
hdulist_r=fits.HDUList([ofd_r, hdu_r])
hdulist_r.writeto(file_r, overwrite=True)
#print('IFS_r.fits is created ')
##################################################################################
def earthshine(self,theta):
"""
For given theta angle, return the earth-shine spectrum.
:param theta: angle (in degree) from the target to earth limb.
:return: the scaled solar spectrum
template_wave: unit in A
template_flux: unit in erg/s/cm^2/A/arcsec^2
"""
# read solar template
solar_template = pd.read_csv(self.information['indata_path']+'/refs/solar_spec.dat', sep='\s+',
header=None, comment='#')
template_wave = solar_template[0].values
template_flux = solar_template[1].values
# read earth shine surface brightness
earthshine_curve = pd.read_csv(self.information['indata_path']+'/refs/earthshine.dat',
header=None, comment='#')
angle = earthshine_curve[0].values
surface_brightness = earthshine_curve[1].values
# read V-band throughtput
cat_filter_V = pd.read_csv(self.information['indata_path']+'/refs/filter_Bessell_V.dat', sep='\s+',
header=None, comment='#')
filter_wave = cat_filter_V[0].values
filter_response = cat_filter_V[1].values
# interplate to the target wavelength in V-band
ind_filter = (template_wave >= np.min(filter_wave)) & (template_wave <= np.max(filter_wave))
filter_wave_interp = template_wave[ind_filter]
filter_response_interp = np.interp(filter_wave_interp, filter_wave, filter_response)
filter_constant = simps(filter_response_interp * filter_wave_interp, filter_wave_interp)
template_constant = simps(filter_response_interp * template_wave[ind_filter] * template_flux[ind_filter],
template_wave[ind_filter])
dwave = filter_wave_interp[1:] - filter_wave_interp[:-1]
wave_eff = np.nansum(dwave * filter_wave_interp[1:] * filter_response_interp[1:]) / \
np.nansum(dwave * filter_response_interp[1:])
# get the normalized value at theta.
u0 = np.interp(theta, angle, surface_brightness) # mag/arcsec^2
u0 = 10**((u0 + 48.6)/(-2.5)) # target flux in erg/s/cm^2/Hz unit
u0 = u0 * 3e18 / wave_eff**2 # erg/s/cm^2/A/arcsec^2
factor = u0 * filter_constant / template_constant
norm_flux = template_flux * factor # erg/s/cm^2/A/arcsec^2
self.earthshine_wave=template_wave # A
self.earthshine_flux=norm_flux
return
########################################################################################################################################################################################################################################################
def zodiacal(self, ra, dec, time):
"""
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
:param time: the specified string that in ISO format i.e., yyyy-mm-dd.
:return:
wave_A: wavelength of the zodical spectrum
spec_mjy: flux of the zodical spectrum, in unit of MJy/sr
spec_erg: flux of the zodical spectrum, in unit of erg/s/cm^2/A/sr
"""
# get solar position
dt = datetime.fromisoformat(time)
jd = julian.to_jd(dt, fmt='jd')
t = Time(jd, format='jd', scale='utc')
astro_sun = get_sun(t)
ra_sun, dec_sun = astro_sun.gcrs.ra.deg, astro_sun.gcrs.dec.deg
radec_sun = SkyCoord(ra=ra_sun*u.degree, dec=dec_sun*u.degree, frame='gcrs')
lb_sun = radec_sun.transform_to('geocentrictrueecliptic')
# get offsets between the target and sun.
radec_obj = SkyCoord(ra=ra*u.degree, dec=dec*u.degree, frame='icrs')
lb_obj = radec_obj.transform_to('geocentrictrueecliptic')
beta = abs(lb_obj.lat.degree)
lamda = abs(lb_obj.lon.degree - lb_sun.lon.degree)
# interpolated zodical surface brightness at 0.5 um
zodi = pd.read_csv(self.information['indata_path']+'/refs/zodi_map.dat', sep='\s+', header=None, comment='#')
beta_angle = np.array([0, 5, 10, 15, 20, 25, 30, 45, 60, 75])
lamda_angle = np.array([0, 5, 10, 15, 20, 25, 30, 35, 40, 45,
60, 75, 90, 105, 120, 135, 150, 165, 180])
xx, yy = np.meshgrid(beta_angle, lamda_angle)
f = interpolate.interp2d(xx, yy, zodi, kind='linear')
zodi_obj = f(beta, lamda) # 10^−8 W m−2 sr−1 um−1
# read the zodical spectrum in the ecliptic
cat_spec = pd.read_csv(self.information['indata_path']+'/refs/solar_spec.dat', sep='\s+', header=None, comment='#')
wave = cat_spec[0].values # A
spec0 = cat_spec[1].values # 10^-8 W m^−2 sr^−1 μm^−1
zodi_norm = 252 # 10^-8 W m^−2 sr^−1 μm^−1
spec = spec0 * (zodi_obj / zodi_norm) * 1e-8 # W m^−2 sr^−1 μm^−1
# convert to the commonly used unit of MJy/sr, erg/s/cm^2/A/sr
wave_A = wave # A
spec_mjy = spec * 0.1 * wave_A**2 / 3e18 * 1e23 * 1e-6 # MJy/sr
spec_erg = spec * 0.1 # erg/s/cm^2/A/sr
spec_erg2 = spec_erg / 4.25452e10 # erg/s/cm^2/A/arcsec^2
self.zodiacal_wave=wave_A # in A
self.zodiacal_flux=spec_erg2
return
###################################################################################
##################################################################################
def CalskyNoise(self, lam):
# calculate sky noise;
planckh= 6.62620*10**-27 # erg s;
cc=2.99792458*10**17 # in nm/s
fov2=0.01 # arcsec^2
# lam is input wavelength in nm
##########################################
self.earthshine_wave # A
self.earthshine_flux # erg/s/cm^2/A/arcsec^2
earthshine_flux=np.interp(lam*10.0, self.earthshine_wave,self.earthshine_flux) # flux from zodiacal
self.zodiacal_wave # in A
self.zodiacal_flux # erg/s/cm^2/A/arcsec^2
zodiacal_flux=np.interp(lam*10.0, self.zodiacal_wave,self.zodiacal_flux) # flux from zodiacal
fluxlam_sky=(earthshine_flux+zodiacal_flux)*fov2 # erg/s/cm2/A
###############
ephoton=planckh*cc/lam # in erg/photon, cc与lambda单位需要一致；
Ns_skynoise=fluxlam_sky/ephoton # in unit of photons/cm2/s/A
return Ns_skynoise
#################################################################################
def sim_sky_img(self,skyfitsfilename, skyRa_shift, skyDec_shift, sky_rot, exposuretime):
############################################################################
### load fits file
indatafile=skyfitsfilename
a=fits.open(indatafile)
#######################################
self.information['name_obj']=a[0].header['OBJECT']
self.information['ra_obj'] =a[0].header['RA'] ### in degree
self.information['dec_obj'] =a[0].header['DEC'] ### in degree
disRa =(skyRa_shift)/3600.0 ##convert unit of degree to arcsec
disDec=(skyDec_shift)/3600.0 ##convert unit of degree to arcsec
self.information['ra_pnt0'] =a[0].header['RA'] + disRa/np.cos(a[0].header['DEC']/180.0*np.pi)
self.information['dec_pnt0']=a[0].header['DEC']+ disDec
self.earthshine(self.earthshine_theta)
self.zodiacal(self.information['ra_obj'], self.information['dec_obj'], self.zodiacal_time)
self.information['target']=deg2HMS(self.information['ra_obj'], self.information['dec_obj'])
### main input data
SpecCube=a[1].data ## spectrum data cube;
Wave=0.1*a[2].data # the relatived wavelength which is converted from Unit A to nm
#print('Wave data header', hdr)
######################################################################################
exptime=self.information['exptime'] #exposure time
dis_dx=disRa # image shift Ra in arcsec
dis_dy=disDec # image shift Dec in arcsec
sizeout= len(SpecCube[:,0,0])
blue_img=galsim.Image(np.zeros((2048,4096)),copy=True)
blue_img.scale=self.pixelscale
blue_img.setOrigin=(0,0)
red_img=galsim.Image(np.zeros((3072,6144)),copy=True)
red_img.scale=self.pixelscale
red_img.setOrigin=(0,0)
blue_sensor=galsim.Sensor()
red_sensor=galsim.Sensor()
deltalam=np.mean(np.diff(Wave))
energy=0.0
energy_blue=0.0
energy_red=0.0
width_blue=0
################################
############## doppler effect to photons.wavelength #############
#self.orbit_pars
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])
vx_sat=float(self.orbit_pars[self.orbit_exp_num,4])
vy_sat=float(self.orbit_pars[self.orbit_exp_num,5])
vz_sat=float(self.orbit_pars[self.orbit_exp_num,6])
self.information['POSI0_X']=x_sat
self.information['POSI0_Y']=y_sat
self.information['POSI0_Z']=z_sat
self.information['VELO0_X']=vx_sat
self.information['VELO0_Y']=vy_sat
self.information['VELO0_Z']=vz_sat
theta1=beta_angle( x_sat,y_sat,z_sat,vx_sat,vy_sat,vz_sat, self.information['ra_obj'], self.information['dec_obj'])
v1=np.sqrt(vx_sat**2+vy_sat**2+vz_sat**2)*np.cos(theta1/180.0*np.pi) # velocity at stat exposure time
vv1=LSR_velocity(self.information['ra_obj'],self.information['dec_obj'],v1, self.TianCe_day)
#################################################
### exposure end time is t2 ;
t2=self.dt+timedelta(seconds=self.information['exptime'])
t2jd=time2jd(t2)
if self.orbit_pars[-1,0]<t2jd: ## orbit parameters are not in currenct txt file
self.orbit_file_num=self.orbit_file_num+1
fn=self.information['indata_path']+'/refs/orbit20160925/'+str(self.orbit_file_num)+'.txt'
self.orbit_pars=np.loadtxt(fn)
self.orbit_exp_num=0
for k in range(self.orbit_exp_num, len(self.orbit_pars),1):
if t2jd-self.orbit_pars[k,0]<0:
break
if k==0:
deltaT=jd2time(self.orbit_pars[k,0])-t2
p1x=self.orbit_pars[k,1]-(self.orbit_pars[k+1,1]-self.orbit_pars[k,1])*deltaT.seconds/120
p1y=self.orbit_pars[k,2]-(self.orbit_pars[k+1,2]-self.orbit_pars[k,2])*deltaT.seconds/120
p1z=self.orbit_pars[k,3]-(self.orbit_pars[k+1,3]-self.orbit_pars[k,3])*deltaT.seconds/120
p1vx=self.orbit_pars[k,4]-(self.orbit_pars[k+1,4]-self.orbit_pars[k,4])*deltaT.seconds/120
p1vx=self.orbit_pars[k,5]-(self.orbit_pars[k+1,5]-self.orbit_pars[k,5])*deltaT.seconds/120
p1vx=self.orbit_pars[k,6]-(self.orbit_pars[k+1,6]-self.orbit_pars[k,6])*deltaT.seconds/120
else:
deltaT=jd2time(self.orbit_pars[k,0])-t2
p1x=self.orbit_pars[k-1,1]+(self.orbit_pars[k,1]-self.orbit_pars[k-1,1])*deltaT.seconds/120
p1y=self.orbit_pars[k-1,2]+(self.orbit_pars[k,2]-self.orbit_pars[k-1,2])*deltaT.seconds/120
p1z=self.orbit_pars[k-1,3]+(self.orbit_pars[k,3]-self.orbit_pars[k-1,3])*deltaT.seconds/120
p1vx=self.orbit_pars[k-1,4]+(self.orbit_pars[k,4]-self.orbit_pars[k-1,4])*deltaT.seconds/120
p1vy=self.orbit_pars[k-1,5]+(self.orbit_pars[k,5]-self.orbit_pars[k-1,5])*deltaT.seconds/120
p1vz=self.orbit_pars[k-1,6]+(self.orbit_pars[k,6]-self.orbit_pars[k-1,6])*deltaT.seconds/120
#######
self.information['POSI1_X']=p1x
self.information['POSI1_Y']=p1y
self.information['POSI1_Z']=p1z
self.information['VELO1_X']=p1vx
self.information['VELO1_Y']=p1vy
self.information['VELO1_Z']=p1vz
theta2=beta_angle(p1x, p1y, p1z, p1vx, p1vy, p1vz,self.information['ra_obj'], self.information['dec_obj'])
v2=np.sqrt(p1vx**2+p1vy**2+p1vz**2)*np.cos(theta2/180.0*np.pi) # velocity at end exposure time
vv2=LSR_velocity(self.information['ra_obj'],self.information['dec_obj'], v2, self.TianCe_day)
#### get slice and slit mask
########################################################################################
for ilam in range(len(Wave)):
#print('ilam=', ilam)
if ilam%500==0:
self.log.info('ilam = %i' % ilam)
lam=Wave[ilam] # the wavelength of the i-th frame data
if lam<350:
continue
###############################################
Specimg=SpecCube[:,:,ilam] # get the i-th frame of the input SpecCube data;
Nspecimg=SpecCube2photon(Specimg,lam) # convert to photons/cm2/s/A
Nskynoise=self.CalskyNoise(lam)
### add sky noise
if self.skyback:
Nimg=Nspecimg+Nskynoise
else:
Nimg=Nspecimg
# multipe the tel area, exposure time, and bandwidh;
Nimg=Nimg*self.telarea*exptime*deltalam*10.0 # photons/cm2/s/A to photons, here Nimg size is 100*100
##########################################################################
### shift image with photons position and rotate them round the image true center
img=galsim.Image(Nimg,copy=True)
img.scale = self.pixelscale
img.setOrigin(0,0)
photons=galsim.PhotonArray.makeFromImage(img, max_flux=max(Nimg.max()/1000.0, 1.0))
### now shift and rotated photons
self.information['shift_dx']=dis_dx
self.information['shift_dy']=dis_dy
photons.x=photons.x-dis_dx ### apply shift to photons position in ra direction
photons.y=photons.y-dis_dy ### apply shift to photons position in dec direction
x2,y2=rotation_yan(img.true_center.x, img.true_center.y, photons.x/img.scale, photons.y/img.scale, sky_rot)
#print('rotation time=', t2-t1)
photons2=galsim.PhotonArray(N=len(x2), x=x2*img.scale, y=y2*img.scale, flux = photons.flux)
photons=photons2
## convert the photon image to galsim photons;
###################################################################
image0 = galsim.Image(sizeout,sizeout)
image0.scale = self.pixelscale
image0.setOrigin(0,0)
rotphotons=galsim.PhotonArray(N=len(photons.x), x=x2, y=y2, flux = photons2.flux)
sensor = galsim.Sensor()
sensor.accumulate(rotphotons, image0)
#####################################
#####################################################################
### do convolve image0 with PSF0 from primay CSST ###
### calculate the PSF0 at this wavelength
Q=lam*1e-3*self.Fnum/self.pixelsize
wavefront0=self.opd0/lam
wavefront=ndimage.rotate(wavefront0, -sky_rot, order=1, reshape=False) # here we choose reshape=False, rotate the wavefront
psf0=anySampledPSF(wavefront, self.pupil , Q , 64)
conv = fftconvolve(image0.array, psf0, mode='same')
conv[conv < 0.0] = 0.0
img0=conv
################ opd1 rms =0.075 @632.8nm
wavefront=self.opd1/lam
wavefront=ndimage.rotate(wavefront, -sky_rot, order=1, reshape=False) # here we choose reshape=False, rotate the wavefront
psf1=anySampledPSF(wavefront, self.pupil , Q , 32)
##### do convolve with psf1
conv = fftconvolve(img0, psf1, mode='same')
conv[conv < 0.0] = 0.0
img0=conv
#######################################################
energy=energy+img0[50-32:50+32,50-32:50+32].sum() ## calculate the slice image energy;
CCD_Qe_lam=np.interp(lam,self.CCD_Qe[:,0], self.CCD_Qe[:,1]) ## CCD quantum efficiency
Qe_blue = np.interp(lam, self.optical_blue_Q[:,0], self.optical_blue_Q[:,1]) # optical efficiency , convert the wavelength to A
Qe_red = np.interp(lam, self.optical_red_Q[:,0], self.optical_red_Q[:,1]) # optical efficiency, convert the wavelength to A
Qe_blue=Qe_blue*CCD_Qe_lam
Qe_red =Qe_red *CCD_Qe_lam
lam1=lam*(1+vv1/(3.0*1e5))
lam2=lam*(1+vv2/(3.0*1e5))
######## consider the slice optical efficiency in different slicer channel
da=fits.open(self.information['indata_path']+'/slicer_QE.fits')
slicer_Qe=da[0].data
img0=img0*slicer_Qe
########## do the slice effect ###################
for k in range(32):
#### do slice effect to get slice image
img1=img0*self.maskSlice[str(k)]
############ get opd2 and PSF2 ######################
wavefront2=self.opd2/lam
psf2=anySampledPSF(wavefront2, self.pupil , Q , 64)
psf2=ndimage.rotate(psf2, -sky_rot, order=1, reshape=False) # here we choose reshape=False, the rotated image will
##### do convolve
##################
conv = fftconvolve(img1, psf2, mode='same')
#suppress negative numbers
conv[conv < 0.0] = 0.0
img2=conv
##############do Slit Mask ###########################
img2=img2*self.maskSlit[str(k)]
######### get opd3 and PSF3 ##########################
wavefront3=self.opd3/lam
psf3=anySampledPSF(wavefront3, self.pupil , Q , 64)
psf3=ndimage.rotate(psf3, -sky_rot, order=1, reshape=False) # here we choose reshape=False, the rotated image will
##### do convolve
########################################
conv = fftconvolve(img2, psf3, mode='same')
#suppress negative numbers
conv[conv < 0.0] = 0.0
img3=conv
######################## get subimage #####################
subimage = galsim.Image(80,80)
subimage.array[:,:]=img3[int(sizeout/2)-40:int(sizeout/2)+40,int(sizeout/2)-40:int(sizeout/2)+40]
######################## get photons from sub-image #####################
subimage.scale = self.pixelscale
subimage.setOrigin(0,0)
photons=galsim.PhotonArray.makeFromImage(subimage, max_flux=max(img3.max()/1000.0,1.0))
#############################################################################
###### 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
np.random.seed()
wavesample=lam1+(lam2-lam1)*np.random.rand(p_num)
if (lam>=350.0 and lam<=650.0):
## bulue channel
photons_blue=galsim.PhotonArray(N=p_num, x=photons.x[idx0], y=photons.y[idx0], flux = Qe_blue*photons.flux[idx0],wavelength=wavesample)
dx_blue, dy_blue = get_dx_dy_blue(wavesample)
photons_blue.x=photons_blue.x/self.pixelscale+dx_blue+self.slice_blue['px'][k]
photons_blue.y=photons_blue.y/self.pixelscale+dy_blue+self.slice_blue['py'][k]
blue_sensor.accumulate(photons_blue, blue_img)
energy_blue=energy_blue+sum(photons_blue.flux)
width_blue=width_blue+deltalam/32.0
if (lam>=560.0) & (lam<=1000.0 ):
## red channel
photons_red=galsim.PhotonArray(N=p_num,x=photons.x[idx0],y=photons.y[idx0], flux = Qe_red*photons.flux[idx0],wavelength=wavesample)
dx_red, dy_red = get_dx_dy_red(wavesample)
photons_red.x=photons_red.x/self.pixelscale+dx_red+self.slice_red['px'][k]
photons_red.y=photons_red.y/self.pixelscale+dy_red+self.slice_red['py'][k]
red_sensor.accumulate(photons_red, red_img)
energy_red=energy_red+sum(photons_red.flux)
####################################################################################
## stray light will cover 2% of input total light;
blue_img.array[:,:]=blue_img.array[:,:]+0.01*energy/2048/4096
red_img.array[:,:] =red_img.array[:,:]+ 0.01*energy/3072/6144
self.image_b=blue_img.array
self.image_r=red_img.array
return
#################################################################################################
#################################################################################################
def simulate(self, skyfitsin, skyRa_shift, skyDec_shift,sky_rot, exptime):
"""
Create a single simulated image of a quadrant defined by the configuration file.
"""
#self.configure() #print the configfile name and path;
self.dt=datetime.now()
self.information['exptime']=exptime
self.skyfilepath=skyfitsin
np.random.seed()
ud = np.random.random() # Choose a random
self.earthshine_theta=ud * 60 # in degree
##################################################################
#### load orbit parameters #####
flag=0
for k in range(1,50,1):
fn=self.information['indata_path']+'/refs/orbit20160925/'+str(k)+'.txt';
d=np.loadtxt(fn);
self.dt_num=int((self.information['exptime']+self.information['readouttime']+125)/120)
now_dt=datetime.utcnow()
now_jd=time2jd(now_dt)
for kk in range(len(d[:,0])):
if now_jd-d[kk,0]<=0:
flag=1
break
if flag==1:
break
#####################end for
self.orbit_pars=d
self.orbit_file_num=k
self.orbit_exp_num =kk
exptime_start_jd=d[kk,0] #### jd time, utc format
self.dt=julian.from_jd(exptime_start_jd, fmt='jd')
self.TianCe_day=self.dt.strftime("%Y-%m-%d") ###str(self.dt.year)+'-'+str(self.dt.month)+'-'+str(self.dt.day)
self.TianCe_exp_start=dt2hmd(self.dt)
self.zodiacal_time=self.TianCe_day
#######################################################################
self.sim_sky_img(skyfitsin, skyRa_shift, skyDec_shift, sky_rot, self.information['exptime'])
self.information['sky_rot']=sky_rot
self.information['skyRa_shift']=skyRa_shift
self.information['skyRa_shift']=skyDec_shift
###############################################################################
############ add some effect to images #######################################
if self.flatfieldM:
self.applyflatfield()
#print('Applying flatfieldM finished.......')
if self.darknoise:
self.applyDarkCurrent()
if self.cosmicRays:
self.addCosmicRays()
#print('Applying cosmicRays finished.......')
if self.bleeding:
self.applyBleeding_yan()
#print('Applying bleeding finished.......')
self.applyPoissonNoise()
if self.nonlinearity:
self.applyNonlinearity()
#print('Applying nonlinearity finished.......')
if self.radiationDamage:
self.applyRadiationDamage()
#print('Applying radiationDamage finished.......')
##### cut original CCD image to four parts by four read out channels and zones
self.CCDreadout()
if self.readoutNoise:
self.applyReadoutNoise()
#print('Applying readoutNoise finished.......')
self.electrons2ADU()
self.applyBias()
if self.cosmetics:
self.applyCosmetics()
#print('Applying cosmetics finished.......')
self.discretise()
self.writeOutputs()
self.log.info('Using the following input values:')
for key, value in self.information.items():
self.log.info('%s = %s' % (key, value))
self.log.info('Using the following booleans:')
for key, value in self.booleans.items():
self.log.info('%s = %s' % (key, value))
self.log.info('Finished the simulation.')
##############################################################################################
##############################################################################################
def runIFSsim(configfile):
opts, args = processArgs()
opts.configfile=configfile
simulate= IFSsimulator(opts)
skyfitsin=simulate.information['skyfitsin']
exptime=simulate.information['exptime']
sky_ra_dis=simulate.information['sky_ra_dis']
sky_dec_dis=simulate.information['sky_dec_dis']
sky_angle_dis=simulate.information['sky_angle_dis']
simulate.simulate(skyfitsin, sky_ra_dis, sky_dec_dis, sky_angle_dis, exptime)
########################## begin main fucntion #######################################
##############################################################################################
##############################################################################################
if __name__ == "__main__":
if len(sys.argv[:]) <2:
configfile='./ifs_data/ifs_sim_example.config'
###########################################################################################
if len(sys.argv[:]) >=2:
configfile=sys.argv[1]
if not os.path.exists(configfile):
print('The given input configfile path is wrong......')
sys.exit(1)
################################################
runIFSsim(configfile)
print('---The CSST-IFS simulation is successful!---')
#'/home/yan/IFS_FabuCode/InputData/IFS_sim_Fabu.config'
############################################################
csst_ifs_sim/support/IFSinstrumentModel.py 0 → 100644
View file @ 1b761b45
"""
VIS Instrument Model
====================
The file provides a function that returns VIS related information such as pixel
size, dark current, gain, zeropoint, and sky background.
:requires: NumPy
:requires: numexpr
:version: 0.7
"""
# import datetime, math
# import numpy as np
# import numexpr as ne
def IFSinformation():
"""
Returns a dictionary describing VIS. The following information is provided (id: value - reference)::
:return: instrument model parameters
:rtype: dict
"""
#########################################################################################################
#out = dict(readnoise=4, pixel_size=0.1, dark=0.0008333, fullwellcapacity=90000, bluesize=4000, redsize=6000, readtime=300.)
out=dict()
out.update({'dob' : 0, 'rdose' : 8.0e9,
'parallelTrapfile' : 'cdm_euclid_parallel.dat', 'serialTrapfile' : 'cdm_euclid_serial.dat',
'beta_s' : 0.6, 'beta_p': 0.6, 'fwc' : 90000, 'vth' : 1.168e7, 't' : 20.48e-3, 'vg' : 6.e-11,
'st' : 5.0e-6, 'sfwc' : 730000., 'svg' : 1.0e-10})
return out
def CCDnonLinearityModel(data, beta=6e-7):
"""
The non-linearity is modelled based on the results presented.
:param data: data to which the non-linearity model is being applied to
:type data: ndarray
:return: input data after conversion with the non-linearity model
:rtype: float or ndarray
"""
out = data-beta*data**2
return out
if __name__ == '__main__':
print('IFSinstrument')
csst_ifs_sim/support/cosmicrays.py 0 → 100644
View file @ 1b761b45
"""
Cosmic Rays
===========
This simple class can be used to include cosmic ray events to an image.
By default the cosmic ray events are drawn from distributions describing
the length and energy of the events. Such distributions can be generated
for example using Stardust code (http://ieeexplore.ieee.org/stamp/stamp.jsp?arnumber=04636917).
The energy of the cosmic ray events can also be set to constant for
testing purposes. The class can be used to draw a single cosmic ray
event or up to a covering fraction.
:requires: NumPy
:requires: SciPy
:version: 0.2
"""
import numpy as np
from scipy.interpolate import InterpolatedUnivariateSpline
class cosmicrays():
"""
Cosmic ray generation class. Can either draw events from distributions or
set the energy of the events to a constant.
:param log: logger instance
:param image: image to which cosmic rays are added to (a copy is made not to change the original numpy array)
:param crInfo: column information (cosmic ray file)
:param information: cosmic ray track information (file containing track length and energy information) and
exposure time.
"""
def __init__(self, log, image, crInfo=None, information=None):
"""
Cosmic ray generation class. Can either draw events from distributions or
set the energy of the events to a constant.
:param log: logger instance
:param image: image to which cosmic rays are added to (a copy is made not to change the original numpy array)
:param crInfo: column information (cosmic ray file)
:param information: cosmic ray track information (file containing track length and energy information) and
exposure time.
"""
#setup logger
self.log = log
#image and size
self.image = image.copy()
self.ysize, self.xsize = self.image.shape
#set up the information dictionary, first with defaults and then overwrite with inputs if given
self.information = (dict(cosmicraylengths='/home/yan/csst-master/data/cdf_cr_length.dat',
cosmicraydistance='/home/yan/csst-master/data/cdf_cr_total.dat',
exptime=565))
if information is not None:
self.information.update(information)
if crInfo is not None:
self.cr = crInfo
else:
self._readCosmicrayInformation()
def _readCosmicrayInformation(self):
self.log.info('Reading in cosmic ray information from %s and %s' % (self.information['cosmicraylengths'],
self.information['cosmicraydistance']))
#read in the information from the files
crLengths = np.loadtxt(self.information['cosmicraylengths'])
crDists = np.loadtxt(self.information['cosmicraydistance'])
#set up the cosmic ray information dictionary
self.cr = dict(cr_u=crLengths[:, 0], cr_cdf=crLengths[:, 1], cr_cdfn=np.shape(crLengths)[0],
cr_v=crDists[:, 0], cr_cde=crDists[:, 1], cr_cden=np.shape(crDists)[0])
return self.cr
def _cosmicRayIntercepts(self, lum, x0, y0, l, phi):
"""
Derive cosmic ray streak intercept points.
:param lum: luminosities of the cosmic ray tracks
:param x0: central positions of the cosmic ray tracks in x-direction
:param y0: central positions of the cosmic ray tracks in y-direction
:param l: lengths of the cosmic ray tracks
:param phi: orientation angles of the cosmic ray tracks
:return: cosmic ray map (image)
:rtype: nd-array
"""
#create empty array
crImage = np.zeros((self.ysize, self.xsize), dtype=np.float64)
#x and y shifts
dx = l * np.cos(phi) / 2.
dy = l * np.sin(phi) / 2.
mskdx = np.abs(dx) < 1e-8
mskdy = np.abs(dy) < 1e-8
dx[mskdx] = 0.
dy[mskdy] = 0.
#pixels in x-direction
ilo = np.round(x0.copy() - dx)
msk = ilo < 0.
ilo[msk] = 0
ilo = ilo.astype(np.int)
ihi = 1 + np.round(x0.copy() + dx)
msk = ihi > self.xsize
ihi[msk] = self.xsize
ihi = ihi.astype(np.int)
#pixels in y-directions
jlo = np.round(y0.copy() - dy)
msk = jlo < 0.
jlo[msk] = 0
jlo = jlo.astype(np.int)
jhi = 1 + np.round(y0.copy() + dy)
msk = jhi > self.ysize
jhi[msk] = self.ysize
jhi = jhi.astype(np.int)
#loop over the individual events
for i, luminosity in enumerate(lum):
n = 0 # count the intercepts
u = []
x = []
y = []
#Compute X intercepts on the pixel grid
if ilo[i] < ihi[i]:
for xcoord in range(ilo[i], ihi[i]):
ok = (xcoord - x0[i]) / dx[i]
if np.abs(ok) <= 0.5:
n += 1
u.append(ok)
x.append(xcoord)
y.append(y0[i] + ok * dy[i])
else:
for xcoord in range(ihi[i], ilo[i]):
ok = (xcoord - x0[i]) / dx[i]
if np.abs(ok) <= 0.5:
n += 1
u.append(ok)
x.append(xcoord)
y.append(y0[i] + ok * dy[i])
#Compute Y intercepts on the pixel grid
if jlo[i] < jhi[i]:
for ycoord in range(jlo[i], jhi[i]):
ok = (ycoord - y0[i]) / dy[i]
if np.abs(ok) <= 0.5:
n += 1
u.append(ok)
x.append(x0[i] + ok * dx[i])
y.append(ycoord)
else:
for ycoord in range(jhi[i], jlo[i]):
ok = (ycoord - y0[i]) / dy[i]
if np.abs(ok) <= 0.5:
n += 1
u.append(ok)
x.append(x0[i] + ok * dx[i])
y.append(ycoord)
#check if no intercepts were found
if n < 1:
xc = int(np.floor(x0[i]))
yc = int(np.floor(y0[i]))
crImage[yc, xc] += luminosity
#Find the arguments that sort the intersections along the track
u = np.asarray(u)
x = np.asarray(x)
y = np.asarray(y)
args = np.argsort(u)
u = u[args]
x = x[args]
y = y[args]
#Decide which cell each interval traverses, and the path length
for i in range(1, n - 1):
w = u[i + 1] - u[i]
cx = int(1 + np.floor((x[i + 1] + x[i]) / 2.))
cy = int(1 + np.floor((y[i + 1] + y[i]) / 2.))
if 0 <= cx < self.xsize and 0 <= cy < self.ysize:
crImage[cy, cx] += (w * luminosity)
return crImage
def _drawCosmicRays(self, limit=None):
"""
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.
"""
#estimate the number of cosmics
cr_n = self.xsize * self.ysize * 0.014 / 43.263316 * 2.
#scale with exposure time, the above numbers are for the nominal 565s exposure
cr_n *= (self.information['exptime'] / 565.0)
#assume a power-law intensity distribution for tracks
fit = dict(cr_lo=1.0e3, cr_hi=1.0e5, cr_q=2.0e0)
fit['q1'] = 1.0e0 - fit['cr_q']
fit['en1'] = fit['cr_lo'] ** fit['q1']
fit['en2'] = fit['cr_hi'] ** fit['q1']
#pseudo-random numbers taken from a uniform distribution between 0 and 1
np.random.seed()
luck = np.random.rand(int(np.floor(cr_n)))
#draw the length of the tracks
if self.cr['cr_cdfn'] > 1:
ius = InterpolatedUnivariateSpline(self.cr['cr_cdf'], self.cr['cr_u'])
self.cr['cr_l'] = ius(luck)
else:
self.cr['cr_l'] = np.sqrt(1.0 - luck ** 2) / luck
#draw the energy of the tracks
if self.cr['cr_cden'] > 1:
ius = InterpolatedUnivariateSpline(self.cr['cr_cde'], self.cr['cr_v'])
self.cr['cr_e'] = ius(luck)
else:
np.random.seed()
self.cr['cr_e'] = (fit['en1'] + (fit['en2'] - fit['en1']) *
np.random.rand(int(np.floor(cr_n)))) ** (1.0 / fit['q1'])
#Choose the properties such as positions and an angle from a random Uniform dist
np.random.seed()
cr_x = self.xsize * np.random.rand(int(np.floor(cr_n)))
np.random.seed()
cr_y = self.ysize * np.random.rand(int(np.floor(cr_n)))
np.random.seed()
cr_phi = np.pi * np.random.rand(int(np.floor(cr_n)))
#find the intercepts
if limit is None:
self.cosmicrayMap = self._cosmicRayIntercepts(self.cr['cr_e'], cr_x, cr_y, self.cr['cr_l'], cr_phi)
print ('Number of cosmic ray events:', len(self.cr['cr_e']))
else:
#limit to electron levels < limit
msk = self.cr['cr_e'] < limit
print ('Number of cosmic ray events: %i / %i' % (len(self.cr['cr_e'][msk]), int(np.floor(cr_n))))
self.cosmicrayMap = self._cosmicRayIntercepts(self.cr['cr_e'][msk], cr_x[msk], cr_y[msk],
self.cr['cr_l'][msk], cr_phi[msk])
#count the covering factor
area_cr = np.count_nonzero(self.cosmicrayMap)
text = 'The cosmic ray covering factor is %i pixels i.e. %.3f per cent' \
% (area_cr, 100.*area_cr / (self.xsize*self.ysize))
self.log.info(text)
print (text)
def _drawSingleEvent(self, limit=1000, cr_n=1):
"""
Generate a single cosmic ray event and include it to a cosmic ray map (self.cosmicrayMap).
:param limit: limiting energy for the cosmic ray event
:type limit: float
:param cr_n: number of cosmic ray events to include
:type cr_n: int
:return: None
"""
#pseudo-random numbers taken from a uniform distribution between 0 and 1
np.random.seed()
luck = np.random.rand(cr_n)
#draw the length of the tracks
ius = InterpolatedUnivariateSpline(self.cr['cr_cdf'], self.cr['cr_u'])
self.cr['cr_l'] = ius(luck)
#set the energy directly to the limit
self.cr['cr_e'] = np.asarray([limit, ]*cr_n)
#Choose the properties such as positions and an angle from a random Uniform dist
np.random.seed()
cr_x = self.xsize * np.random.rand(int(np.floor(cr_n)))
np.random.seed()
cr_y = self.ysize * np.random.rand(int(np.floor(cr_n)))
np.random.seed()
cr_phi = np.pi * np.random.rand(int(np.floor(cr_n)))
#find the intercepts
self.cosmicrayMap = self._cosmicRayIntercepts(self.cr['cr_e'], cr_x, cr_y, self.cr['cr_l'], cr_phi)
#count the covering factor
area_cr = np.count_nonzero(self.cosmicrayMap)
text = 'The cosmic ray covering factor is %i pixels i.e. %.3f per cent' \
% (area_cr, 100.*area_cr / (self.xsize*self.ysize))
self.log.info(text)
print( text)
def _drawEventsToCoveringFactor(self, coveringFraction=3.0, limit=1000, verbose=False):
"""
Generate cosmic ray events up to a covering fraction and include it to a cosmic ray map (self.cosmicrayMap).
:param coveringFraction: covering fraction of cosmic rya events in per cent of total number of pixels
:type coveringFraction: float
:param limit: limiting energy for the cosmic ray event [None = draw from distribution]
:type limit: None or float
:param verbose: print out information to stdout
:type verbose: bool
:return: None
"""
self.cosmicrayMap = np.zeros((self.ysize, self.xsize))
#how many events to draw at once, too large number leads to exceeding the covering fraction
cr_n = int(295 * self.information['exptime'] / 565. * coveringFraction / 1.4)
covering = 0.0
while covering < coveringFraction:
#pseudo-random numbers taken from a uniform distribution between 0 and 1
np.random.seed()
luck = np.random.rand(cr_n)
#draw the length of the tracks
ius = InterpolatedUnivariateSpline(self.cr['cr_cdf'], self.cr['cr_u'])
self.cr['cr_l'] = ius(luck)
if limit is None:
ius = InterpolatedUnivariateSpline(self.cr['cr_cde'], self.cr['cr_v'])
self.cr['cr_e'] = ius(luck)
else:
#set the energy directly to the limit
self.cr['cr_e'] = np.asarray([limit,])
#Choose the properties such as positions and an angle from a random Uniform dist
np.random.seed()
cr_x = self.xsize * np.random.rand(int(np.floor(cr_n)))
np.random.seed()
cr_y = self.ysize * np.random.rand(int(np.floor(cr_n)))
np.random.seed()
cr_phi = np.pi * np.random.rand(int(np.floor(cr_n)))
#find the intercepts
self.cosmicrayMap += self._cosmicRayIntercepts(self.cr['cr_e'], cr_x, cr_y, self.cr['cr_l'], cr_phi)
#count the covering factor
area_cr = np.count_nonzero(self.cosmicrayMap)
covering = 100.*area_cr / (self.xsize*self.ysize)
text = 'The cosmic ray covering factor is %i pixels i.e. %.3f per cent' % (area_cr, covering)
self.log.info(text)
if verbose:
print( text)
def addCosmicRays(self, limit=None):
"""
Include cosmic rays to the image given.
:return: image with cosmic rays
:rtype: ndarray
"""
self._drawCosmicRays(limit=limit)
#paste cosmic rays
self.image += self.cosmicrayMap
return self.image
def addSingleEvent(self, limit=None):
"""
Include a single cosmic ray event to the image given.
:return: image with cosmic rays
:rtype: ndarray
"""
self._drawSingleEvent(limit=limit)
#paste cosmic rays
self.image += self.cosmicrayMap
return self.image
def addUpToFraction(self, coveringFraction, limit=None, verbose=False):
"""
Add cosmic ray events up to the covering Fraction.
:param coveringFraction: covering fraction of cosmic rya events in per cent of total number of pixels
:type coveringFraction: float
:param limit: limiting energy for the cosmic ray event [None = draw from distribution]
:type limit: None or float
:param verbose: print out information to stdout
:type verbose: bool
:return: image with cosmic rays
:rtype: ndarray
"""
self._drawEventsToCoveringFactor(coveringFraction, limit=limit, verbose=verbose)
#paste cosmic rays
self.image += self.cosmicrayMap
return self.image
if __name__ == "__main__":
print()
csst_ifs_sim/support/logger.py 0 → 100644
View file @ 1b761b45
"""
These functions can be used for logging information.
.. Warning:: logger is not multiprocessing safe.
:version: 0.3
"""
import logging
import logging.handlers
def setUpLogger(log_filename, loggername='logger'):
"""
Sets up a logger.
:param: log_filename: name of the file to save the log.
:param: loggername: name of the logger
:return: logger instance
"""
# create logger
logger = logging.getLogger(loggername)
logger.setLevel(logging.DEBUG)
# Add the log message handler to the logger
handler = logging.handlers.RotatingFileHandler(log_filename)
#maxBytes=20, backupCount=5)
# create formatter
formatter = logging.Formatter('%(asctime)s - %(module)s - %(funcName)s - %(levelname)s - %(message)s')
# add formatter to ch
handler.setFormatter(formatter)
# add handler to logger
if (logger.hasHandlers()):
logger.handlers.clear()
logger.addHandler(handler)
return logger
setup.py 0 → 100644
View file @ 1b761b45
import setuptools
with open("README.md", "r") as fh:
long_description = fh.read()
setuptools.setup(
name='csst_ifs_sim',
version='2.0.0-IFS1.0.0',
author='CSST Team',
author_email='zhaojunyan@shao.ac.cn',
description='The CSST - ifs - sim', # short description
long_description=long_description,
long_description_content_type="text/markdown",
url='https://csst-tb.bao.ac.cn/',
# project_urls={
# 'Source': 'https://csst-tb.bao.ac.cn/code/csst-l1/ifs/csst_ifs_sim',
# },
packages=setuptools.find_packages(),
license='MIT',
classifiers=["Development Status :: 5 - Production/Stable",
"Intended Audience :: Science/Research",
"License :: OSI Approved :: MIT License",
"Operating System :: OS Independent",
"Programming Language :: Python :: 3.8",
"Topic :: Scientific/Engineering :: Physics",
"Topic :: Scientific/Engineering :: Astronomy"],
package_dir={'csst_ifs_sim': 'csst_ifs_sim'},
# include_package_data=True,
package_data={"": ["LICENSE", "README.md"],
"csst_ifs_sim": ["ifs_so/*", "ifs_data/*", "ifs_data/refs/*", "ifs_data/refs/orbit20160925/*"]},
# install_requires=['sphinx',
# 'numpy',
# 'scipy', 'matplotlib',
# 'astropy', 'healpy', 'ccdproc', 'deepCR', 'photutils'],
python_requires='>=3.8',
)
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