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

Merge branch 'master' into 'main' 

build

See merge request shaosim/mci!1
parents 3a19183e f957f062
Pipeline #567 failed with stages
in 1 minute and 8 seconds
Hide whitespace changes
Inline Side-by-side
LICENSE 0 → 100644
View file @ aecb88aa
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 @ aecb88aa
# mci
# 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/mci.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/mci/-/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_mci_sim/CTI/CTI.py 0 → 100644
View file @ aecb88aa
"""
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_mci_sim/csst_mci_sim.py 0 → 100644
View file @ aecb88aa
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
@author: yan, zhaojun"""
"""
The CSST MCI Image Simulator
=============================================
This file contains an image simulator for the CSST MCI.
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, 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).
#. Loop over the number of exposures to co-add and for each object in the galaxy fits file:
#. 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 [optional].
#. Add photon (Poisson) noise [optional]
#. Add cosmetic defects from an input file [optional].
#. Add prescan and 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 [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, a WCS is assigned
and the output is saved to the current working directory.
.. Warning:: The code is still work in progress and new features are being added.
The code has been tested, but nevertheless bugs may be lurking in corners, so
please report any weird or inconsistent simulations to the author.
Contact Information: zhaojunyan@shao.ac.cn
-------
"""
import galsim
import pandas as pd
from scipy.integrate import simps
import julian
from datetime import datetime, timedelta
from astropy.time import Time
from astropy.coordinates import get_sun
import numpy as np
from tqdm import tqdm
from astropy.wcs import WCS as WCS
from astropy.io import fits
from astropy import units as u
import os, sys, math
import configparser as ConfigParser
from optparse import OptionParser
from scipy import ndimage
from astropy.coordinates import SkyCoord
from scipy import interpolate
from scipy.signal import fftconvolve
sys.path.append('..')
from CTI import CTI
from support import logger as lg
from support import cosmicrays
from support import MCIinstrumentModel
from support import shao
from support import sed
# set the folder
FOLDER ='../'
#####################################################################################
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()
###############################################################################
def make_c_coor(fov, step):
"""
Draw the mesh grids for a fov * fov box with given step .
"""
nc=int(fov/step)
bs=fov
ds = bs / nc
xx01 = np.linspace(-bs / 2.0, bs / 2.0 - ds, nc) + 0.5 * ds
xx02 = np.linspace(-bs / 2.0, bs / 2.0 - ds, nc) + 0.5 * ds
xg2, xg1 = np.meshgrid(xx01, xx02)
return xg1, xg2
##############################################################################
###############################################################################
def str2time(strTime):
if len(strTime)>20:#
msec=int(float('0.'+strTime[20:])*1000000) # microsecond
str2=strTime[0:19]+' '+str(msec)
return datetime.strptime(str2,'%Y %m %d %H %M %S %f')
#datetime to mjd
def time2mjd(dateT):
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 cut_radius(rcutstar, mstar, mag):
return rcutstar * 10**(0.4*2*(mstar-mag)/3.125)
def core_radius(rcorestar, mstar, mag):
return rcorestar * 10**(0.4*(mstar-mag)/2)
def v_disp(sigstar, mstar, mag):
return sigstar * 10**(0.4*(mstar-mag)/3.57)
###############################################################################
#####################################################################################
#### distortion function
def distortField(ra, dec, ch):
###% ra ,dec are the idea position in arcsec , and the center position is ra=0, dec=0
'''MCI geometry distortion effect in channel G R and I'''
distortA=dict()
distortB=dict()
distortA['g']=np.array([0.000586224851462036,0.999778507355332,-0.000377391397200834,-4.40403573019393e-06,0.000147444528630966,-1.93281825050268e-06,5.73520238714447e-07,-1.05194725549731e-08,7.55124671251098e-07,-3.85623216175251e-09,-1.36254798567168e-11,1.36983654024573e-09,-4.53104347730230e-11,1.50641840169747e-09,-1.05226890727366e-11,1.06471556810228e-12,-1.42299485385165e-16,5.90008722878028e-12,2.36749977009538e-17,4.11061448730916e-12,3.39287277469805e-17])
distortB['g']=np.array([1.00113878750680,-0.000495107770623867,0.999162436209357,9.67138672907941e-05,-3.88808001621361e-06,0.000267940195644359,-3.03504629416955e-09,7.59508802931395e-07,-1.13952684052500e-08,9.58672893217047e-07,2.96397490595616e-10,-3.01506961685930e-11,2.22570315950441e-09,-4.26077028191289e-11,2.07336026666512e-09,1.71450079597168e-16,2.94455108664049e-12,8.18246797239659e-17,8.32235684990184e-12,1.05203957047210e-17,5.58541337682320e-12])
distortA['r']=np.array([0.000530712822898294,0.999814397089242,-0.000366783811357609,-1.08650906916023e-05,0.000146500108063480,-4.48994291741769e-06,5.66992328511032e-07,-2.10826363791224e-08,7.51050898843171e-07,-7.74459106063614e-09,-5.78712436084704e-11,1.36389366137393e-09,-9.50057131576629e-11,1.49666815592076e-09,-2.16074939825761e-11,2.07124539578673e-12,-1.32787329448528e-13,5.78542850850562e-12,-5.82328830750271e-14,4.04881815088026e-12,2.46095156355282e-15])
distortB['r']=np.array([1.00110972320988,-0.000483253233767592,0.999181377618578,9.60316928622784e-05,-9.06574382424422e-06,0.000266147076486786,-6.64513480754565e-09,7.53794491639207e-07,-2.25340150955077e-08,9.53032549996219e-07,2.93844965469408e-10,-6.80521155195455e-11,2.21272371816576e-09,-9.13097728847483e-11,2.06225316601308e-09,9.41096324034730e-14,3.00253083795091e-12,-1.40935245750903e-13,8.27122551593984e-12,-2.41500808188601e-13,5.52074803508046e-12])
distortA['i']=np.array([0.000532645905256854,0.999813271238129,-0.000366606839338804,-1.06782246147986e-05,0.000146529811926289,-4.54488847969004e-06,5.68028930846942e-07,-2.11055358580630e-08,7.51187628288423e-07,-7.93885390157909e-09,-5.63104333653064e-11,1.36559362569725e-09,-9.64353839395137e-11,1.50231201562648e-09,-2.04159956020294e-11,1.91408535007684e-12,-7.46369323634455e-14,5.86071470639700e-12,-1.65328163262914e-13,4.07648238374705e-12,3.40880674871652e-14])
distortB['i']=np.array([1.00111276400037,-0.000484310769918440,0.999180286636823,9.61240720951134e-05,-8.76526511019577e-06,0.000266192433565878,-6.59462433108999e-09,7.54391611102465e-07,-2.30957980169170e-08,9.53612393886641e-07,2.95558631849358e-10,-7.08307092409029e-11,2.21952307845185e-09,-9.03816603147003e-11,2.06450141393973e-09,-3.77836686311826e-14,3.13358046393416e-12,-5.20417845240954e-14,8.24487152802918e-12,-1.17846469609206e-13,5.35830020655191e-12])
x=ra/0.05*10/1000 # convert ra in arcsec to mm
y=dec/0.05*10/1000 # convert ra in arcsec to mm
dd=np.array([1, x, y, x*x, x*y, y*y, x**3, x**2*y, x*y**2, y**3, x**4, x**3*y, x**2*y**2, x*y**3, y**4, x**5, x**4*y, x**3*y**2, x**2*y**3 , x*y**4, y**5])
xc= np.dot(distortA[ch],dd)
yc= np.dot(distortB[ch],dd)
distortRa =xc*1000/10*0.05-0.00265
distortDec=yc*1000/10*0.05-5.0055
return distortRa, distortDec
#####################################################################################
def world_to_pixel(sra,
sdec,
rotsky,
tra,
tdec,
x_center_pixel,
y_center_pixel,
pixelsize=0.05):
"""_summary_
Parameters
----------
sra : np.array
star ra such as np.array([22,23,24]),unit:deg
sdec : np.array
stat dec such as np.array([10,20,30]),unit:deg
rotsky : float
rotation angel of the telescope,unit:deg
tra : float
telescope pointing ra,such as 222.2 deg
rdec : float
telescope pointing dec,such as 33.3 deg
x_center_pixel : float
x refer point of MCI,usually the center of the CCD,which start from (0,0)
y_center_pixel : float
y refer point of MCI,usually the center of the CCD,which start from(0,0)
pixelsize : float
pixelsize for MCI ccd, default :0.05 arcsec / pixel
Returns
-------
pixel_position:list with two np.array
such as [array([124.16937605, 99. , 99. ]),
array([ 97.52378661, 99. , 100.50014483])]
"""
theta_r = rotsky / 180 * np.pi
w = WCS(naxis=2)
w.wcs.crpix = [x_center_pixel + 1, y_center_pixel + 1] # pixel from (1, 1)
w.wcs.cd = np.array([[
-np.cos(-theta_r) * pixelsize / 3600,
np.sin(-theta_r) * pixelsize / 3600
],
[
np.sin(-theta_r) * pixelsize / 3600,
np.cos(-theta_r) * pixelsize / 3600
]])
w.wcs.crval = [tra, tdec]
w.wcs.ctype = ["RA---TAN", "DEC--TAN"]
pixel_position = w.all_world2pix(sra, sdec, 0)
return pixel_position
###############################################################################
def cal_pos(center_ra, center_dec, rotTel, rotSky, sRa, sDec):
# ''' center_ra: telescope pointing in ra direction, unit: degree
# center_dec: telescope pointing in dec direction, unit: degree
# rotTel: telescope totation in degree
# rotSky: telescope rotation in degree
# sRa: source ra in arcsec
# sDec: source dec in arcsec
boresight = galsim.CelestialCoord(ra=-center_ra*galsim.degrees, dec=center_dec*galsim.degrees)
q = rotTel - rotSky
cq, sq = np.cos(q), np.sin(q)
jac = [cq, -sq, sq, cq]
affine = galsim.AffineTransform(*jac)
wcs = galsim.TanWCS(affine, boresight, units=galsim.arcsec)
objCoord = galsim.CelestialCoord(ra =-sRa*galsim.arcsec, dec=sDec*galsim.arcsec)
field_pos = wcs.toImage(objCoord)
return -field_pos.x, field_pos.y ## return the field position
#####################################################################################
def krebin(a, sample):
""" Fast Rebinning with flux conservation
New shape must be an integer divisor of the current shape.
This algorithm is much faster than rebin_array
Parameters
----------
a : array_like
input array
sample : rebinned sample 2 or 3 or 4 or other int value
"""
# Klaus P's fastrebin from web
size=a.shape
shape=np.zeros(2,dtype=int)
shape[0]=int(size[0]/sample)
shape[1]=int(size[1]/sample)
sh = shape[0], a.shape[0] // shape[0], shape[1], a.shape[1] // shape[1]
return a.reshape(sh).sum(-1).sum(1)
#####################################################################
#########################################################
def centroid(data):
'''
calculate the centroid of the input two-dimentional image data
Parameters
----------
data : input image.
Returns
-------
cx: the centroid column number, in horizontal direction definet in python image show
cy: the centroid row number , in vertical direction
'''
###
from scipy import ndimage
cy,cx=ndimage.center_of_mass(data)
return cx, cy
###############################################################################
###############################################################################
#########################################################################
'''
PSF interpolation for MCI_sim
'''
import scipy.spatial as spatial
###find neighbors-KDtree ###
def findNeighbors(tx, ty, px, py, dn=5):
"""
find nearest neighbors by 2D-KDTree
Parameters:
tx, ty (float, float): a given position
px, py (numpy.array, numpy.array): position data for tree
dr (float-optional): distance in arcsec
dn (int-optional): nearest-N
OnlyDistance (bool-optional): only use distance to find neighbors. Default: True
Returns:
indexq (numpy.array): index
"""
datax = px
datay = py
tree = spatial.KDTree(list(zip(datax.ravel(), datay.ravel())))
indexq=[]
dist, indexq = tree.query([tx, ty], dn)
return indexq
###############################################################################
###PSF-IDW###
def psfMaker_IDW(px, py, PSFMat, cen_col, cen_row, IDWindex=5, OnlyNeighbors=True):
"""
psf interpolation by IDW
Parameters:
px, py (float, float): position of the target
PSFMat (numpy.array): PSF matrix array at given position and wavelength
cen_col, cen_row (numpy.array, numpy.array): potions of the psf centers in arcsec
IDWindex (int-optional): the power index of IDW
OnlyNeighbors (bool-optional): only neighbors are used for psf interpolation
Returns:
psfMaker (numpy.array)
"""
minimum_psf_weight = 1e-8
ref_col = px ### target position in x
ref_row = py ### target position in y
ngy, ngx = PSFMat[:, :, 0, 0].shape ### PSF data size
npsf = PSFMat[:, :, :,:].shape[3] ### total psf number in one wavelength
psfWeight = np.zeros([npsf])
if OnlyNeighbors == True:
neigh = findNeighbors(px, py, cen_col, cen_row, dn=5)
# if hoc is not None:
# neigh = findNeighbors_hoclist(cen_col, cen_row, tx=px,ty=py, dn=4, hoc=hoc, hoclist=hoclist)
neighFlag = np.zeros(npsf)
neighFlag[neigh] = 1
for ipsf in range(npsf):
if OnlyNeighbors == True:
if neighFlag[ipsf] != 1:
continue
dist = np.sqrt((ref_col - cen_col[ipsf])**2 + (ref_row - cen_row[ipsf])**2)
if IDWindex == 1:
psfWeight[ipsf] = dist
if IDWindex == 2:
psfWeight[ipsf] = dist**2
if IDWindex == 3:
psfWeight[ipsf] = dist**3
if IDWindex == 4:
psfWeight[ipsf] = dist**4
psfWeight[ipsf] = max(psfWeight[ipsf], minimum_psf_weight)
psfWeight[ipsf] = 1./psfWeight[ipsf]
psfWeight /= np.sum(psfWeight)
psfMaker = np.zeros([ngy, ngx, 7 ], dtype=np.float32)
for waven in range(7):
for ipsf in range(npsf):
if OnlyNeighbors == True:
if neighFlag[ipsf] != 1:
continue
iPSFMat = PSFMat[:, :, waven,ipsf].copy()
ipsfWeight = psfWeight[ipsf]
psfMaker[:, :, waven] += iPSFMat * ipsfWeight
psfMaker[:, :, waven] /= np.nansum(psfMaker[:, :, waven])
return psfMaker
###############################################################################
###############################################################################
###############################################################################
class MCIsimulator():
"""
CSST MCI Simulator
The image that is being build is in::
self.image
self.image_g g channel
self.image_r r channel
self.image_i i 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
###########################################################################
self.information = MCIinstrumentModel.MCIinformation()
#update settings with defaults
self.information.update(dict(quadrant=int(0),
ccdx=int(0),
ccdy=int(0),
ccdxgap=1.643,
ccdygap=8.116,
xsize=2000,
ysize=2000,
fullwellcapacity=90000,
dark=0.001,
exptime=300.0,
readouttime=4.,
rdose=8.0e9,
ghostCutoff=22.0,
ghostRatio=5.e-5,
coveringfraction=1.0,
pixel_size=0.05))
self.configure(1) #print the configfile name and path;
self.information.update(dict(
cosmicraylengths =self.information['indata_path']+'/data/cdf_cr_length.dat',
cosmicraydistance=self.information['indata_path']+'/data/cdf_cr_total.dat',
parallelTrapfile =self.information['indata_path']+'/data/cdm_euclid_parallel.dat',
serialTrapfile =self.information['indata_path']+'/data/cdm_euclid_serial.dat',
cosmeticsfile_g =self.information['indata_path']+'/data/Cosmetics_g.txt',
cosmeticsfile_r =self.information['indata_path']+'/data/Cosmetics_r.txt' ,
cosmeticsfile_i =self.information['indata_path']+'/data/Cosmetics_i.txt' ))
###############################################################################
def readConfigs(self,simnumber):
"""
Reads the config file information using configParser and sets up a logger.
"""
self.config = ConfigParser.RawConfigParser()
#self.config.readfp(open(self.configfile))
self.config.read_file(open(self.configfile))
################################################################################################3
def processConfigs(self):
"""
Processes configuration information and save the information to a dictionary self.information.
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.MCIinstrumentModel.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)
#ghost ratio can be in engineering format, so getfloat does not capture it...
try:
self.information['ghostRatio'] = float(self.config.get(self.section, 'ghostRatio'))
except:
pass
#########
#booleans to control the flow
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.overscans = self.config.getboolean(self.section, 'overscans')
self.nonlinearity = self.config.getboolean(self.section, 'nonlinearity')
self.readoutNoise = self.config.getboolean(self.section, 'readoutNoise')
self.skyback = self.config.getboolean(self.section, 'skyback')
self.TianceEffect = self.config.getboolean(self.section, 'TianceEffect')
self.intscale = self.config.getboolean(self.section, 'intscale')
self.ghosts = self.config.getboolean(self.section, 'ghosts')
self.shutterEffect =self.config.getboolean(self.section, 'shutterEffect')
self.flatfieldM =self.config.getboolean(self.section, 'flatfieldm')
self.PRNUeffect =self.config.getboolean(self.section, 'PRNUeffect')
self.appFatt =self.config.getboolean(self.section, 'appFatt')
self.sky_shift_rot =self.config.getboolean(self.section, 'sky_shift_rot')
self.distortion =self.config.getboolean(self.section, 'distortion')
self.sim_star =self.config.getboolean(self.section, 'sim_star')
self.sim_galaxy =self.config.getboolean(self.section, 'sim_galaxy')
###############################################3#################################
self.booleans = dict(cosmicRays =self.cosmicRays,
darknoise =self.darknoise,
cosmetics =self.cosmetics,
radiationDamage=self.radiationDamage,
bleeding =self.bleeding,
overscans =self.overscans,
nonlinearity =self.nonlinearity,
readoutNoise =self.readoutNoise,
skyback =self.skyback,
TianceEffect =self.TianceEffect,
intscale =self.intscale,
ghosts =self.ghosts ,
shutterEffect =self.shutterEffect,
flatfieldM =self.flatfieldM,
PRNUeffect =self.PRNUeffect,
appFatt =self.appFatt ,
sky_shift_rot =self.sky_shift_rot,
distortion =self.distortion,
sim_star =self.sim_star,
sim_galaxy =self.sim_galaxy)
#####################################################################
now=datetime.now()
#data_time=now.strftime("%Y-%m-%d-%H-%M-%S")
result_day=now.strftime("%Y-%m-%d")
self.result_path=self.information['result_path']+'/'+result_day ### CSST1
if os.path.isdir(self.result_path)==False:
os.mkdir(self.result_path)
os.mkdir(self.result_path+'/cali_Data')
os.mkdir(self.result_path+'/log_Data')
os.mkdir(self.result_path+'/sky_Data')
#############################################################
data_time=now.strftime("%Y-%m-%d-%H-%M-%S")
self.log = lg.setUpLogger(self.result_path+'/log_Data/MCIsim_'+'_'+data_time+'.log')
self.log.info('-------STARTING A NEW SIMULATION------------')
self.log.info(self.information)
##############################################
#load instrument model, these values are also stored in the FITS header
self.information['G_filters']=["F275W", "F280N","NUV", "WU", "CBU", "F343N", "u", "F373N", "F395N", "F336W"]
self.information['R_filters']=["F487N", "F502N", "CBV", "r", "F656N", "F658N", "F467M", "F555W", "F606W", "F673N"]
self.information['I_filters']=["z", "y", "F815N", "CBI", "F925N", "F960M", "F968N", "F845M" ,"F850LP" ,"F814W"]
#### load telescope efficiency data
self.tel_eff=np.load(self.information['indata_path']+'/tel_eff/tel_eff.npy',allow_pickle=True).item()
#### load MCI filter data
self.filterP=np.load(self.information['indata_path']+'/MCI_filters/mci_filterPWTC.npy',allow_pickle=True).item()
#####################################################################
self.log.info('Using the following input values:')
for key, value in self.information.items():
self.log.info('%s = %s' % (key, value))
self.log.info('Using the following booleans:')
for key, value in self.booleans.items():
self.log.info('%s = %s' % (key, value))
return
#########################################################################################
def make_c_coor(self, bs, nc):
'''
Draw the mesh grids for a bs*bs box with nc*nc pixels
'''
ds=bs/nc
xx01 = np.linspace(-bs/2.0,bs/2.0-ds,nc)+0.5*ds
xx02 = np.linspace(-bs/2.0,bs/2.0-ds,nc)+0.5*ds
xg2,xg1 = np.meshgrid(xx01,xx02)
return xg1,xg2
##########################################################################################
def _createEmpty(self):
"""
Creates and empty array of a given x and y size full of zeros.
add g r i channel images;
Creates lensing parameters;
"""
self.image_g=np.zeros((self.information['ysize'], self.information['xsize']), dtype=np.float64)
self.image_r=np.zeros((self.information['ysize'], self.information['xsize']), dtype=np.float64)
self.image_i=np.zeros((self.information['ysize'], self.information['xsize']), dtype=np.float64)
return
#########################################################################################################################
#########################################################################################################
def _loadGhostModel(self):
"""
Reads in a ghost model from a FITS file and stores the data to self.ghostModel.
Currently assumes that the ghost model has already been properly scaled and that the pixel
scale of the input data corresponds to the nominal VIS pixel scale. Futhermore, assumes that the
distance to the ghost from y=0 is appropriate (given current knowledge, about 750 VIS pixels).
"""
self.ghostOffsetX=self.information['ghostoffsetx']
self.ghostOffsetY=self.information['ghostoffsety']
self.ghostMax = self.information['ghostratio']
self.log.info('Maximum in the ghost model %e' % self.ghostMax)
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])
return
###############################################################################
###############################################################################
def configure(self,simnumber):
"""
Configures the simulator with input information and creates and empty array to which the final image will
be build on.
"""
self.readConfigs(simnumber)
self.processConfigs()
self._createEmpty()
self.log.info('Read in the configuration files and created an empty array')
return
################################################################################
def load_filter_PSF(self):
##### load filter PSF in three channels ,G , R and I
###
# load filter name in three channels from information
self.filter_g=self.information['filter_g']
self.filter_r=self.information['filter_r']
self.filter_i=self.information['filter_i']
self.filter_psf=dict()
filtername=self.filter_g
self.filter_psf['g']=np.load(self.information['indata_path']+'/PSF/'+filtername+'_PSF.npy', allow_pickle=True).item()
filtername=self.filter_r
self.filter_psf['r']=np.load(self.information['indata_path']+'/PSF/'+filtername+'_PSF.npy', allow_pickle=True).item()
filtername=self.filter_i
self.filter_psf['i']=np.load(self.information['indata_path']+'/PSF/'+filtername+'_PSF.npy', allow_pickle=True).item()
return
###############################################################################################
def get_PSF(self, px, py, ch):
"""
Get the PSF at a given image position
Parameters:
px,py : target position
filternmae : the selected filter name
ch : the MCI channle , g or r or i
Returns:
PSF: PSF array at given 7 wavelength and at target position
"""
PSFMat = self.filter_psf[ch]['psf_mat']
cen_col= self.filter_psf[ch]['psf_field_X']
cen_row= self.filter_psf[ch]['psf_field_Y']
imPSF = psfMaker_IDW(px, py, PSFMat, cen_col, cen_row, IDWindex=2)
return imPSF
###############################################################################
###############################################################################
def cal_SED_photons(self, flux):
## flux_arr, input flux array with unit of 1e-17 erg/s/A/cm^2
wave = np.linspace(2500,10000,8501) # set the wavelenth for MCI from 250nm to 1100nm
wave=wave/10
# calcutle photons from original spec cube data,
# data: input data, two-dimentional, unit : 1e-17 erg/s/A/cm^2
planckh= 6.62620*10**-27 # % erg s;
cc=2.99792458*10**17 # in nm/s
telarea=3.1415926*100*100 # in cm^2，望远镜聚光面积�?
fluxlam=1e-17*(flux) # convert original unit to unit of erg/s/A/cm^2
# wave in nm ;;
ephoton=planckh*cc/wave # single photon energy in erg/photon; cc与lambda单位需要一致；
Nphoton =fluxlam/ephoton*telarea*self.information['exptime'] # in unit of photons
#### load telescope efficiency data
teleff=dict()
teleff['g']=self.tel_eff['G'] ### fisrt colunmm is wavelength in nm
teleff['r']=self.tel_eff['R'] ### second colunmm is efficiency
teleff['i']=self.tel_eff['I']
#####
ft=dict()
ft['g']=self.filter_g
ft['r']=self.filter_r
ft['i']=self.filter_i
# load filter name in three channels
ft_wave=dict()
ft_wave['g']=self.filterP[ft['g']]['wave_nm']
ft_wave['r']=self.filterP[ft['r']]['wave_nm']
ft_wave['i']=self.filterP[ft['i']]['wave_nm']
ft_eff=dict()
ft_eff['g']=self.filterP[ft['g']]['throughout']
ft_eff['r']=self.filterP[ft['r']]['throughout']
ft_eff['i']=self.filterP[ft['i']]['throughout']
chlist=['g','r','i']
Sumphoton=dict()
PSF_eff_weight=dict()
for k in range(3):
ch=chlist[k]
tel_wavearr=teleff[ch][:,0]
tel_effarr =teleff[ch][:,1]
tel_effnew=np.interp(wave, tel_wavearr, tel_effarr)
ft_wavearr=ft_wave[ch]
ft_effarr =ft_eff[ch]
ft_effnew =np.interp(wave, ft_wavearr, ft_effarr)
#### SED flux*filter_flux*tel_eff
eff_arr=Nphoton*ft_effnew*tel_effnew
Sumphoton[ch]=np.sum(eff_arr)
iwave=self.filter_psf[ch]['psf_iwave'][:] ### seven wavelength for the seleced filter
psf_coeff =np.interp(iwave, wave, eff_arr)
if psf_coeff.sum()==0:
psf_coeff=1.0/7*np.ones(7)
else:
psf_coeff=psf_coeff/psf_coeff.sum()
PSF_eff_weight[ch]=psf_coeff
return Sumphoton, PSF_eff_weight
##################################################
###############################################################################
def cal_sky_noise(self):
####### add earthshine #######
####### add earthshine #######
#self.earthshine_wave # A
#self.earthshine_flux # erg/s/cm^2/A/arcsec^2
wave = np.linspace(2500,11000, 8501) # set the wavelenth for MCI from 250nm to 1100nm
wavearr =self.earthshine_wave # A
fluxarr =self.earthshine_flux # erg/s/cm^2/A/arcsec^2
earthshinefluxarr=np.interp(wave, wavearr, fluxarr) # erg/s/cm^2/A/arcsec^2
wavearr =self.zodiacal_wave # A
fluxarr =self.zodiacal_flux # erg/s/cm^2/A/arcsec^2
zodiacalfluxarr=np.interp(wave, wavearr, fluxarr) # erg/s/cm^2/A/arcsec^2
skynoise_flux=earthshinefluxarr+zodiacalfluxarr # erg/s/cm^2/A/arcsec^2
wave = np.linspace(2500,10000,8501) # set the wavelenth for MCI from 250nm to 1100nm
wave=wave/10
# data: input data, two-dimentional, unit : 1e-17 erg/s/A/cm^2
planckh= 6.62620*10**-27 # % erg s;
cc=2.99792458*10**17 # in nm/s
telarea=3.1415926*100*100 # in cm^2，
ephoton=planckh*cc/wave # in erg/photon,
Nphoton_flux =skynoise_flux/ephoton*telarea*0.05*0.05*self.information['exptime'] # each pixel will have the photons in unit of photons/s/A
#### load telescope efficiency data
teleff=dict()
teleff['g']=self.tel_eff['G'] ### fisrt colunmm is wavelength in nm
teleff['r']=self.tel_eff['R'] ### second colunmm is efficiency
teleff['i']=self.tel_eff['I']
#####
ft=dict()
ft['g']=self.filter_g
ft['r']=self.filter_r
ft['i']=self.filter_i
# load filter name in three channels
ft_wave=dict()
ft_wave['g']=self.filterP[ft['g']]['wave_nm']
ft_wave['r']=self.filterP[ft['r']]['wave_nm']
ft_wave['i']=self.filterP[ft['i']]['wave_nm']
ft_eff=dict()
ft_eff['g']=self.filterP[ft['g']]['throughout']
ft_eff['r']=self.filterP[ft['r']]['throughout']
ft_eff['i']=self.filterP[ft['i']]['throughout']
chlist=['g','r','i']
self.Sky_Noise=dict()
for k in range(3):
ch=chlist[k]
tel_wavearr=teleff[ch][:,0]
tel_effarr =teleff[ch][:,1]
tel_effnew=np.interp(wave, tel_wavearr, tel_effarr)
ft_wavearr=ft_wave[ch]
ft_effarr =ft_eff[ch]
ft_effnew =np.interp(wave, ft_wavearr, ft_effarr)
self.Sky_Noise[ch]=np.sum(Nphoton_flux*ft_effnew*tel_effnew)
return self.Sky_Noise
##################################################
##############################################################################
def cal_StarSED_img(self):
#####################
pixelscale=self.information['pixel_size'] ## arcsec, pixel scale size
ghostOffsetX =self.information['ghostoffsetx']
ghostOffsetY =self.information['ghostoffsety']
ghostMx =self.information['ghostratio']
rotTelPos=self.information['rotTelPos']
rotSkyPos=self.information['rotSkyPos']
theta = rotTelPos - rotSkyPos
############ load star data catlog #####################
df0=pd.read_csv(self.information['indata_path']+'/star_input/GaiaSource_675688-675713.csv')
####
df=df0.fillna(1000)
df2=df[ (df['phot_rp_mean_mag']>16)&(df['phot_bp_mean_mag']>16) & (df['phot_bp_mean_mag']<31) & (df['pmra']<100) & (df['phot_rp_mean_mag']<31)]
df2.index = range(len(df2))
## limit the filed of view to 600*600 arcsec^2 square zone;
df3=df2[ (abs(df2['ra']-df2['ra'].mean())<300/3600.0) & (abs(df2['dec']-df2['dec'].mean())<300/3600.0) ]
df3.index = range(len(df3))
self.star=df3
del df0,df,df2,df3
self.information['ra_obj'] = self.star['ra'].mean()
self.information['dec_obj'] = self.star['dec'].mean()
center_ra =self.information['T_disRa'] +self.information['ra_obj']
center_dec=self.information['T_disDec'] +self.information['dec_obj']
self.information['ra_pnt0'] =center_ra
self.information['dec_pnt0']=center_dec
channel=['g','r','i']
#######################################################################
nsrcs=len(self.star)
##################################################################
obj=dict()
obj['g']=np.zeros(3)
obj['r']=np.zeros(3)
obj['i']=np.zeros(3)
fullimg=dict()
star_input=dict()
star_input['g']=np.zeros((nsrcs,3))
star_input['r']=np.zeros((nsrcs,3))
star_input['i']=np.zeros((nsrcs,3))
star_output=dict()
star_output['g']=np.zeros((nsrcs,3))
star_output['r']=np.zeros((nsrcs,3))
star_output['i']=np.zeros((nsrcs,3))
final_image=dict()
final_image['g'] = galsim.ImageF(int(self.information['xsize']), int(self.information['ysize']))
final_image['r'] = galsim.ImageF(int(self.information['xsize']), int(self.information['ysize']))
final_image['i'] = galsim.ImageF(int(self.information['xsize']), int(self.information['ysize']))
final_image['g'].setOrigin(0,0)
final_image['r'].setOrigin(0,0)
final_image['i'].setOrigin(0,0)
#######################################################################
nlayccd = 0
#### calculate sky noise #####
self.earthshine(self.earthshine_theta)
self.zodiacal(self.information['ra_obj'], self.information['dec_obj'], self.dt.strftime("%Y-%m-%d"))
self.cal_sky_noise()
#################################
self.information['target'] =deg2HMS(self.information['ra_obj'], self.information['dec_obj'])
self.information['CRVAL1']=self.information['ra_pnt0']
self.information['CRVAL2']=self.information['dec_pnt0']
self.information['CRPIX1']=self.information['xsize']/2-0.5 ####
self.information['CRPIX2']=self.information['ysize']/2-0.5 ####
self.information['CD1_1']=-np.cos(-theta)*self.information['pixel_size']/3600.0 ####
self.information['CD1_2']= np.sin(-theta)*self.information['pixel_size']/3600.0 ####
self.information['CD2_1']= np.sin(-theta)*self.information['pixel_size']/3600.0 ####
self.information['CD2_2']= np.cos(-theta)*self.information['pixel_size']/3600.0 ####
#######################################################################
if self.TianceEffect:
ra_list = self.star['ra'].tolist()
dec_list = self.star['dec'].tolist()
pmra_list = self.star['pmra'].tolist()
pmdec_list = self.star['pmdec'].tolist()
parallax_list = self.star['parallax'].tolist()
rv_list = [0.0 for i in range(len(ra_list))]
################################################
newRa, newDec = shao.onOrbitObsPosition(ra_list, dec_list, pmra_list, \
pmdec_list, rv_list, parallax_list, len(ra_list), \
self.information['pos_x'], self.information['pos_y'], self.information['pos_z'], self.information['velocity_x'],self.information['velocity_y'],self.information['velocity_z'], "J2000", self.TianCe_day, self.TianCe_exp_start)
else:
newRa =self.star['ra']
newDec =self.star['dec']
######################################################
#################### generate star image ##########
for j in tqdm(range(nsrcs)): ### nsrcs length
starRa = 3600.0*( newRa[j] ) # ra of star, arcsecond
starDec = 3600.0*( newDec[j] ) # dec of star, arcsecond
###################################################################
fsx,fsy=cal_pos(center_ra,center_dec, rotTelPos, rotSkyPos, starRa, starDec)
row= fsy/pixelscale+self.information['ysize']/2-0.5 ### row number on CCD image
col= fsx/pixelscale+self.information['xsize']/2-0.5 ### col number on CCD image
if row>=self.information['ysize']+100 or row<=-100 or col>=self.information['xsize']+100 or col <=-100:
continue
nlayccd=nlayccd+1
################################################################
for i in range(3):
###
ch=channel[i]
star_input[ch][int(nlayccd)-1, 0] =fsx #ra
star_input[ch][int(nlayccd)-1, 1] =fsy #dec
############# do field distottion ##########
if self.distortion:
for i in range(3):
ch=channel[i]
fpx,fpy=distortField(fsx, fsy, ch)
obj[ch][0]=fpx
obj[ch][1]=fpy
star_output[ch][int(nlayccd)-1, 0] =fpx #ra
star_output[ch][int(nlayccd)-1, 1] =fpy #dec
else:
fpx=fsx
fpy=fsy
for i in range(3):
ch=channel[i]
obj[ch][0]=fpx
obj[ch][1]=fpy
star_output[ch][int(nlayccd)-1, 0] =fpx #ra
star_output[ch][int(nlayccd)-1, 1] =fpy #dec
#######################################################################
####### use SED_code to generate star SED #######
# SED of j-th star
if j<nsrcs-8:
bp = self.star['phot_bp_mean_mag'][j]
rp = self.star['phot_rp_mean_mag'][j]
else:
bp = binary_star[j-nsrcs, 2]
rp = bp
wave = np.linspace(2500, 11000, 8501)
# Loading stellar SED template
t=sed.Star_Temp(self.information['indata_path'])
# Calculating the magnitude (Bp, Rp) of each template
t.toMag()
# Calculating flux
star_flux = sed.Model_Stellar_SED(wave, bp, rp, t, self.information['indata_path']) # unit of 10-17 erg/s/A/cm2
#########################################
##cal_SED_photons(self, flux_arr):
intscales,PSF_eff_weight=self.cal_SED_photons(star_flux)
if j<nsrcs-8:
self.log.info('star number=%i, original Ra(in deg)=%f, original Dec(in deg)=%f, bp=%f,rp=%f, TotPhotons_G=%f, TotPhotons_R=%f,TotPhotons_I=%f' % (nlayccd, newRa[j], newDec[j],self.star['phot_bp_mean_mag'][j], self.star['phot_rp_mean_mag'][j],intscales['g'],intscales['r'],intscales['i']))
gx=dict()
gy=dict()
for i in range(3):
ch=channel[i]
gx[ch]= obj[ch][1]/pixelscale+self.information['ysize']/2-0.5 ### row number on CCD image
gy[ch]= obj[ch][0]/pixelscale+self.information['xsize']/2-0.5 ### col number on CCD image
self.log.info('Channel in =%s, PosX(in pixel)=%f, PosY(in pixel)=%f' % (channel[i], gx[ch], gy[ch]))
######################################################################
psf=dict()
for i in range(3):
ch=channel[i]
#######
psfmat=self.get_PSF(fsx, fsy, ch)
for iwave in range(7):
if iwave==0:
temp=PSF_eff_weight[ch][iwave]*psfmat[:,:,iwave]
else:
temp=temp+PSF_eff_weight[ch][iwave]*psfmat[:,:,iwave]
temp=temp/temp.sum()
####rotate the PSF data
if abs(theta.deg>0):
psf[ch]=ndimage.rotate(temp, theta.deg, order=1, reshape=True) # here we choose reshape=False, the rotated image will
else:
psf[ch]=temp
conv = psf[ch]
conv=conv/conv.sum()
conv=conv*intscales[ch]
#################################################
stamp_img = galsim.Image(conv, copy=True)
stamp_img.setOrigin(0,0)
stamp_img.scale=0.025
#################################################
photons=galsim.PhotonArray.makeFromImage(stamp_img, max_flux=max(1.0, stamp_img.array.max()/10000.0))
cx0,cy0=centroid(stamp_img.array)
if self.appFatt :
### apply treering and bright fatter and diffusion;
SimpleTreeRing = galsim.SiliconSensor().simple_treerings(amplitude=self.information['treering'])
cx, cy = centroid(conv)
pos = galsim.PositionD(x=int(cx), y=int(cy)) ### set subimge center as the pos
siliconsensor = galsim.SiliconSensor(strength=self.information['fatter'], diffusion_factor=self.information['difusion'], treering_func=SimpleTreeRing, treering_center=pos)
image = galsim.Image(conv.shape[0], conv.shape[1])
image.scale = 0.025
image.setOrigin(0,0)
photons.x=photons.x/image.scale
photons.y=photons.y/image.scale
siliconsensor.accumulate(photons, image)
photons=galsim.PhotonArray.makeFromImage(image, max_flux=1.0)
cx0,cy0=centroid(image.array)
############################################
############################################
photons.x=photons.x-cx0*0.025
photons.y=photons.y-cy0*0.025
##########################################################
#### add ghost image ####
if self.ghosts:
ghostphotons=galsim.PhotonArray(len(photons.x),x=photons.x, y=photons.y, flux=photons.flux)
ghostphotons.flux=ghostMx*ghostphotons.flux
ghostphotons.x=photons.x/0.05+ghostOffsetY
ghostphotons.y=photons.y/0.05+ghostOffsetX
if ghostphotons.flux.max()>0.1/1500:
ghostphotons.addTo(final_image[ch])
#########################################################
photons.x=photons.x/0.05+gy[ch]
photons.y=photons.y/0.05+gx[ch]
photons.addTo(final_image[ch])
###############################################################
###############################################################
#print('total star nummber on CCD is: ',nlayccd )
#################################################################
#######################################################################
fullimg['g']=final_image['g'].array
fullimg['r']=final_image['r'].array
fullimg['i']=final_image['i'].array
return fullimg
################################################################################
################################################################################
#################################################################################
########################################################################
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)
#xx, yy = np.meshgrid(beta_angle, lamda_angle,indexing='ij', sparse=True)
f = interpolate.interp2d(xx, yy, zodi, kind='linear')
#f = interpolate.RegularGridInterpolator((xx, yy), zodi, method='linear')
zodi_obj = f(beta, lamda) #
# 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 #
zodi_norm = 252 #
spec = spec0 * (zodi_obj / zodi_norm) * 1e-8 #
# 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 cal_Lensing_galaxy_img(self):
#####################
self.dsx_arc = self.information['pixel_size'] # arcsec, pixel scale size;
self.bsz_arc = max(self.information['ysize'] ,self.information['xsize'])*self.dsx_arc # arcsec, field size;
self.bsz_deg = self.bsz_arc/3600.
#############################################################
### galaxy catlog information
ra = 74.0563449685417
dec = -34.155241212932374
# halo_id = 229600100382
self.information['star_ra'] = 249.70093286177536
self.information['star_dec'] =-44.60485406119162
self.information['gal_ra'] = ra
self.information['gal_dec']= dec
self.information['ra_obj'] =self.information['star_ra']
self.information['dec_obj'] =self.information['star_dec']
#####################################################################
self.earthshine(self.earthshine_theta)
self.zodiacal(self.information['star_ra'], self.information['star_dec'], self.dt.strftime("%Y-%m-%d"))
self.cal_sky_noise()
############ load galaxy data with SED ############################
losscale =self.information['pixel_size']
rotTelPos=self.information['rotTelPos']
rotSkyPos=self.information['rotSkyPos']
theta = rotTelPos - rotSkyPos
center_ra =self.information['T_disRa'] +self.information['star_ra']
center_dec =self.information['T_disDec'] +self.information['star_dec']
self.information['ra_pnt0'] =center_ra
self.information['dec_pnt0']=center_dec
self.information['target'] =deg2HMS(self.information['star_ra'], self.information['star_dec'])
self.information['CRVAL1']=self.information['ra_pnt0']
self.information['CRVAL2']=self.information['dec_pnt0']
self.information['CRPIX1']=self.information['xsize']/2-0.5 ####
self.information['CRPIX2']=self.information['ysize']/2-0.5 ####
self.information['CD1_1'] =-np.cos(-theta)*self.information['pixel_size']/3600.0 ####
self.information['CD1_2'] = np.sin(-theta)*self.information['pixel_size']/3600.0 ####
self.information['CD2_1'] = np.sin(-theta)*self.information['pixel_size']/3600.0 ####
self.information['CD2_2'] = np.cos(-theta)*self.information['pixel_size']/3600.0 ####
obj=dict()
obj['g']=np.zeros(3)
obj['r']=np.zeros(3)
obj['i']=np.zeros(3)
fullimg=dict()
galaxy_input=dict()
nsrcs = 72000
galaxy_input['g']=np.zeros((nsrcs,3))
galaxy_input['r']=np.zeros((nsrcs,3))
galaxy_input['i']=np.zeros((nsrcs,3))
galaxy_output=dict()
galaxy_output['g']=np.zeros((nsrcs,3))
galaxy_output['r']=np.zeros((nsrcs,3))
galaxy_output['i']=np.zeros((nsrcs,3))
channel=['g','r','i']
final_image=dict()
final_image['g'] = galsim.ImageF(int(self.information['xsize']), int(self.information['ysize']))
final_image['r'] = galsim.ImageF(int(self.information['xsize']), int(self.information['ysize']))
final_image['i'] = galsim.ImageF(int(self.information['xsize']), int(self.information['ysize']))
final_image['g'].setOrigin(0,0)
final_image['r'].setOrigin(0,0)
final_image['i'].setOrigin(0,0)
####################################################################
nlayccd = 0
####################generate loc image
for k2 in (range(13)): ###
print('k2=',k2)
#
filename=self.information['indata_path']+'/galaxy_Input/Lens_SED_IMG_V3_0.025/Lens_img_cut_IMG_'+str(k2+1)+'.fits'
if not os.path.exists(filename):
#print('finish load all the input galaxy image fits files already')
break
srcs_cat=fits.open(filename)
#### load galaxy SED fitsfile ###
filename=self.information['indata_path']+'/galaxy_Input/Lens_SED_IMG_V3_0.025/Lens_img_cut_SED_'+str(k2+1)+'.fits'
srcs_sed=fits.open(filename)
###################################################################
for k1 in tqdm(range(1,len(srcs_cat))):
galRa = 3600*(srcs_cat[k1].header['new_ra'] -self.information['gal_ra'] +self.information['star_ra']) # ra of galaxies, arcsecond
galDec = 3600*(srcs_cat[k1].header['new_dec']-self.information['gal_dec']+self.information['star_dec']) # dec of galaxies, arcsecond
#################################################################
fsx,fsy=cal_pos(center_ra,center_dec, rotTelPos, rotSkyPos, galRa, galDec)
row= fsy/losscale+self.information['ysize']/2-0.5 ### row number on CCD image, dec direction ,y direction
col= fsx/losscale+self.information['xsize']/2-0.5 ### col number on CCD image, ra direction, x direction
if row>=self.information['ysize']+100 or row<=-100 or col>=self.information['xsize']+100 or col <=-100:
continue
# # SED of j-th galaxy ,# unit of 10-17 erg/s/A/cm2
gal_flux=srcs_sed[k1].data
################################
### rotate the lensed_images_g ###
if abs(theta.deg)>0:
lensed_images_g=ndimage.rotate(srcs_cat[k1].data, -theta.deg, order=1, reshape=True) # here we choose reshape=False, the rotated image will
else:
lensed_images_g=srcs_cat[k1].data
if lensed_images_g.sum()<=0:
continue
nlayccd=nlayccd+1
################################################################
for i in range(3):
###
ch=channel[i]
galaxy_input[ch][int(nlayccd)-1, 0] =fsx #ra
galaxy_input[ch][int(nlayccd)-1, 1] =fsy #dec
#
############# do field distottion ##########
if self.distortion:
for i in range(3):
ch=channel[i]
fpx,fpy=distortField(fsx, fsy, ch)
obj[ch][0]=fpx
obj[ch][1]=fpy
galaxy_output[ch][int(nlayccd)-1, 0] =fpx #ra
galaxy_output[ch][int(nlayccd)-1, 1] =fpy #dec
else:
fpx=fsx
fpy=fsy
for i in range(3):
ch=channel[i]
obj[ch][0]=fpx
obj[ch][1]=fpy
galaxy_output[ch][int(nlayccd)-1, 0] =fpx #ra
galaxy_output[ch][int(nlayccd)-1, 1] =fpy #dec
#######################################################################
self.log.info('Galaxy number=%i, Ra(in arcsec)=%f, Dec(in arcsec)=%f' % (nlayccd, fpx, fpy))
##cal_SED_photons(self, flux_arr):
intscales,PSF_eff_weight=self.cal_SED_photons(gal_flux)
img=dict()
lensed_images_g=lensed_images_g/lensed_images_g.sum()
img['g']=lensed_images_g*intscales['g']
img['r']=lensed_images_g*intscales['r']
img['i']=lensed_images_g*intscales['i']
################################################
gx=dict()
gy=dict()
for i in range(3):
ch=channel[i]
gy[ch]= obj[ch][1]/losscale+self.information['ysize']/2-0.5 ### row number on CCD image
gx[ch]= obj[ch][0]/losscale+self.information['xsize']/2-0.5 ### col number on CCD image
######################################################################
psf=dict()
for i in range(3):
ch=channel[i]
psfmat=self.get_PSF(fsx, fsy, ch)
temp=0
for iwave in range(7):
temp=temp+PSF_eff_weight[ch][iwave]*psfmat[:,:,iwave]
temp=temp/temp.sum()
##### rotate the PSF data
if abs(theta.deg)>0:
psf[ch]=ndimage.rotate(temp, theta.deg, order=1, reshape=False) # here we choose reshape=False, the rotated image will
else:
psf[ch]=temp
conv = fftconvolve(img[ch], psf[ch], mode='full')
#suppress negative numbers
conv[conv < 0.0] = 0.0
#################################################
stamp_img = galsim.Image(conv, copy=True)
stamp_img.setOrigin(0,0)
stamp_img.scale=0.025
photons=galsim.PhotonArray.makeFromImage(stamp_img, max_flux=1.0)
cx0,cy0=centroid(stamp_img.array)
########### apply fat effect #####
if self.appFatt =='yes':
### apply treering and bright fatter and diffusion;
SimpleTreeRing = galsim.SiliconSensor().simple_treerings(amplitude=self.information['treering'])
cx, cy = centroid(conv)
pos = galsim.PositionD(x=int(cx), y=int(cy)) ### set subimge center as the pos
siliconsensor = galsim.SiliconSensor(strength=self.information['fatter'], diffusion_factor=self.information['difusion'], treering_func=SimpleTreeRing, treering_center=pos)
image = galsim.Image(conv.shape[1],conv.shape[0])
image.scale = 0.025
image.setOrigin(0,0)
photons.x=photons.x/image.scale
photons.y=photons.y/image.scale
siliconsensor.accumulate(photons, image)
photons=galsim.PhotonArray.makeFromImage(image, max_flux=1.0)
cx0,cy0=centroid(image.array)
##############################################################
##############################################################
##################################################################
##################################################################
photons.x=photons.x-cx0*0.025
photons.y=photons.y-cy0*0.025
photons.x=photons.x/0.05+gx[ch]
photons.y=photons.y/0.05+gy[ch]
photons.addTo(final_image[ch])
####################################################
########################################################################
##print('total los nummber on CCD is: ',nlayccd )
#################################################################
#################################################################
fullimg['g']=final_image['g'].array
fullimg['r']=final_image['r'].array
fullimg['i']=final_image['i'].array
return fullimg
###############################################################################
###############################################################################
def generatePRNU(self, ave=1.0, sigma=0.01):
"""
Creates a PRNU field image with given properties.
:return: PRNU field image
:rtype: ndarray
"""
self.log.info('Generating a flat field...')
self.log.info('The flat field has mean value of 1 and a given fluctuations, usually either 1 or 2 percent defined by sigma= %d...' % sigma)
self.PRNU=dict()
np.random.seed(self.filterP[self.filter_g]['seed'])
self.PRNU['g'] = np.random.normal(loc=ave, scale=sigma, size=(self.information['ysize'], self.information['xsize']))
np.random.seed(self.filterP[self.filter_r]['seed'])
self.PRNU['r'] = np.random.normal(loc=ave, scale=sigma, size=(self.information['ysize'], self.information['xsize']))
np.random.seed(self.filterP[self.filter_i]['seed'])
self.PRNU['i'] = np.random.normal(loc=ave, scale=sigma, size=(self.information['ysize'], self.information['xsize']))
return self.PRNU
###############################################################################
###############################################################################
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 applyFlat(self):
#### apply Flat field to image
self.FlatMat=dict()
###
self.FlatMat['g']=self.MakeFlatMatrix(self.image_g, self.filterP[self.filter_g]['seed'])
self.image_g*=self.FlatMat['g']
###
self.FlatMat['r']=self.MakeFlatMatrix(self.image_r, self.filterP[self.filter_r]['seed'])
self.image_r*=self.FlatMat['r']
###
self.FlatMat['i']=self.MakeFlatMatrix(self.image_i, self.filterP[self.filter_i]['seed'])
self.image_g*=self.FlatMat['i']
return
######################################################################################################
def applyPRNUeffect(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.
"""
### generate flatfiledfile, with
prnu_sigma=self.information['prnu_sigma']
self.generatePRNU(ave=1.0, sigma=prnu_sigma)
self.image_g *= self.PRNU['g']
self.image_r *= self.PRNU['r']
self.image_i *= self.PRNU['i']
self.log.info('Applied flatfield to images.')
##################################################################################################
##################################################################################################
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 = np.zeros((self.information['ysize'], self.information['xsize']), dtype=np.float64)
#cosmic ray instance
cosmics_g = cosmicrays.cosmicrays(self.log, crImage, crInfo=self.cr)
cosmics_r = cosmicrays.cosmicrays(self.log, crImage, crInfo=self.cr)
cosmics_i = cosmicrays.cosmicrays(self.log, crImage, crInfo=self.cr)
#add cosmic rays up to the covering fraction
CCD_cr_g = cosmics_g.addUpToFraction(self.information['coveringfraction'], limit=None)
CCD_cr_r = cosmics_r.addUpToFraction(self.information['coveringfraction'], limit=None)
CCD_cr_i = cosmics_i.addUpToFraction(self.information['coveringfraction'], limit=None)
#paste the information
self.image_g += CCD_cr_g
self.image_r += CCD_cr_r
self.image_i += CCD_cr_i
#save cosmic ray image map
self.cosmicMap_g = CCD_cr_g
self.cosmicMap_r = CCD_cr_r
self.cosmicMap_i = CCD_cr_i
#count the covering factor
area_cr_g = np.count_nonzero(self.cosmicMap_g)
area_cr_r = np.count_nonzero(self.cosmicMap_r)
area_cr_i = np.count_nonzero(self.cosmicMap_i)
#self.log.info('The cosmic ray covering factor is %i pixels ' % area_cr)
self.log.info('The cosmic ray in G channel covering factor is %i pixels ' % area_cr_g)
self.log.info('The cosmic ray in R channel covering factor is %i pixels ' % area_cr_r)
self.log.info('The cosmic ray in I channel covering factor is %i pixels ' % area_cr_i)
########################################################################################
def ShutterEffectMat(self, img, t_exp, t_shutter=1.3, dist_bearing=4000, dt=0.001):
# Generate Shutter-Effect normalized image
# t_shutter: time of shutter movement
# dist_bearing: distance between two bearings of shutter leaves
# dt: delta_t of sampling
SampleNumb = int(t_shutter/dt+1)
DistHalf = dist_bearing/2
t = np.arange(SampleNumb)*dt
a_arr = 5.84*np.sin(2*math.pi/t_shutter*t)
v = np.zeros(SampleNumb)
theta = np.zeros(SampleNumb)
x = np.arange(SampleNumb)/(SampleNumb-1)*dist_bearing
s = np.zeros(SampleNumb)
s1 = np.zeros(SampleNumb)
s2 = np.zeros(SampleNumb)
brt = np.zeros(SampleNumb)
idx = np.arange(SampleNumb)
sidx = np.zeros(SampleNumb)
s1idx = np.zeros(SampleNumb)
s2idx = np.zeros(SampleNumb)
v[0] = 0
theta[0] = 0
for i in range(SampleNumb-1):
v[i+1] = v[i]+a_arr[i]*dt
theta[i+1] = theta[i]+v[i]*dt
s1[i] = DistHalf*np.cos(theta[i])
s2[i] = dist_bearing-DistHalf*np.cos(theta[i])
s1idx[i] = int(s1[i]/dist_bearing*(SampleNumb))
s2idx[i] = int(s2[i]/dist_bearing*(SampleNumb))
brt[(idx>s1idx[i]) & (idx<s2idx[i])] += dt
if t_exp>t_shutter*2:
brt = brt*2+(t_exp-t_shutter*2)
else:
brt = brt*2
x = (x-dist_bearing/2)*100
intp = interpolate.splrep(x, brt, s=0)
xmin = 0 #
xmax = img.shape[1] # row
ymin = 0 #
ymax = img.shape[0] # column
if xmin<np.min(x) or xmax>np.max(x):
raise LookupError("Out of focal-plane bounds in X-direction.")
if ymin<0 or ymax>25331:
raise LookupError("Out of focal-plane bounds in Y-direction.")
sizex = xmax-xmin
sizey = ymax-ymin
xnewgrid = np.mgrid[xmin:(xmin+sizex)]
expeffect = interpolate.splev(xnewgrid, intp, der=0)
expeffect /= np.max(expeffect)
exparrnormal = np.tile(expeffect, (sizey,1))
# Image *= exparrnormal
return exparrnormal
########################################################################################
def addshuttereffect(self):
''' apply shutter effect to image'''
#
self.shutterMat=self.ShutterEffectMat(self.image_g, self.information['exptime'])
self.image_g *= self.shutterMat
self.image_r *= self.shutterMat
self.image_i *= self.shutterMat
return
###############################################################################
###############################################################################
def applyDarkCurrent(self):
"""
Apply dark current. Scales the dark with the exposure time.
Additionally saves the image without noise to a FITS file.
"""
#add dark
dark = self.information['exptime'] * self.information['dark']
self.image_g += dark
self.image_r += dark
self.image_i += dark
self.log.info('Added dark current = %f' % dark)
return
###########################################################################
###########################################################################
def applyPoissonSkyNoise(self):
"""
Add Poisson sky noise to the image.
"""
skynoise_g=self.Sky_Noise['g']*np.ones_like(self.image_g)
np.random.seed()
self.image_g =self.image_g+ np.random.poisson(lam=skynoise_g)
self.log.info('Added Poisson sky noise on channel g')
skynoise_r=self.Sky_Noise['r']*np.ones_like(self.image_r)
np.random.seed()
self.image_r =self.image_r+ np.random.poisson(lam=skynoise_r)
self.log.info('Added Poisson sky noise on channel r')
skynoise_i=self.Sky_Noise['i']*np.ones_like(self.image_i)
np.random.seed()
self.image_i =self.image_i+ np.random.poisson(lam=skynoise_i)
self.log.info('Added Poisson sky noise on channel i')
return
###############################################################################
###############################################################################
def applyCosmetics(self):
"""
Apply cosmetic defects described in the input file.
#Number of hot and dead pixels from MSSL/Euclid/TR/12003 Issue 2 Draft b
.. Warning:: This method does not work if the input file has exactly one line.
"""
#######################################################################
cosmetics = np.loadtxt(self.information['indata_path']+'/data/Cosmetics_g.txt')
x = np.round(cosmetics[:, 0]).astype(int) ## row number
y = np.round(cosmetics[:, 1]).astype(int) ## col number
value = np.round(cosmetics[:, 1]).astype(int)
#cosmetics_g=np.zeros((4636,235526))
self.log.info('Adding cosmetic defects to G channel:' )
for xc, yc, val in zip(x, y, value):
if 0 < xc < 4616 and 27 < (yc % 1499) <1179:
self.image_g[xc, yc] = val
#cosmetics_g[yc,xc]=val
self.log.info('x=%i, y=%i, value=%f' % (xc, yc, val))
######################################################################
#######################################################################
cosmetics = np.loadtxt(self.information['indata_path']+'/data/Cosmetics_r.txt')
x = np.round(cosmetics[:, 0]).astype(int) ## row number
y = np.round(cosmetics[:, 1]).astype(int) ## col number
value = np.round(cosmetics[:, 1]).astype(int)
#cosmetics_r=np.zeros((4636,235526))
self.log.info('Adding cosmetic defects to R channel:' )
for xc, yc, val in zip(x, y, value):
if 0 < xc < 4616 and 27 < (yc % 1499) <1179:
self.image_r[xc, yc] = val
#cosmetics_r[yc,xc]=val
self.log.info('x=%i, y=%i, value=%f' % (xc, yc, val))
##############################################################################
#######################################################################
cosmetics = np.loadtxt(self.information['indata_path']+'/data/Cosmetics_i.txt')
x = np.round(cosmetics[:, 0]).astype(int) ## row number
y = np.round(cosmetics[:, 1]).astype(int) ## col number
value = np.round(cosmetics[:, 1]).astype(int)
#cosmetics_i=np.zeros((4636,235526))
self.log.info('Adding cosmetic defects to I channel:' )
for xc, yc, val in zip(x, y, value):
if 0 < xc < 4616 and 27 < (yc % 1499) <1179:
self.image_i[xc, yc] = val
#cosmetics_i[yc,xc]=val
self.log.info('x=%i, y=%i, value=%f' % (xc, yc, val))
return
#############################################################################################################################
def applyRadiationDamage(self):
"""
Applies CDM03 radiation model to the image being constructed.
.. seealso:: Class :`CDM03`
"""
#save image without CTI
self.noCTI_g = self.image_g.copy()
self.log.info('Save NoCti Nonoise image fits file to %s ' % ('noctinonoise' + 'mci_g.fits'))
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_g = cti.applyRadiationDamage(self.image_g.copy().transpose(), iquadrant=self.information['quadrant']).transpose()
self.log.info('Radiation damage added.')
#################################################
self.noCTI_r = self.image_r.copy()
self.log.info('Save NoCti Nonoise image fits file to %s ' % ('noctinonoise' + 'mci_r.fits'))
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.')
##################################################
self.noCTI_i = self.image_i.copy()
self.log.info('Save NoCti Nonoise image fits file to %s ' % ('noctinonoise' + 'mci_i.fits'))
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_i = cti.applyRadiationDamage(self.image_i.copy().transpose(), iquadrant=self.information['quadrant']).transpose()
self.log.info('Radiation damage added.')
return
###########################################################################
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_g = MCIinstrumentModel.CCDnonLinearityModel(self.image_g.copy())
self.log.info('Non-linearity effects included.')
########################################################################
self.log.debug('Starting to apply non-linearity model...')
self.image_r = MCIinstrumentModel.CCDnonLinearityModel(self.image_r.copy())
self.log.info('Non-linearity effects included.')
########################################################################
self.log.debug('Starting to apply non-linearity model...')
self.image_i = MCIinstrumentModel.CCDnonLinearityModel(self.image_i.copy())
self.log.info('Non-linearity effects included.')
return
###############################################################################
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.
"""
if self.information['xsize']==9216 and self.information['ysize']==9232 :
xmin=dict()
xmax=dict()
ymin=dict()
ymax=dict()
############## calculate the non-over subimg of the final matrix of the image
for k in range(1,17):
xmin[k]=0
xmax[k]=4616+self.overscan
ymin[k]=(1152+self.prescan+self.overscan)*(k-1)
ymax[k]=(1152+self.prescan+self.overscan)*(k)
###############
np.random.seed()
noise_g = np.random.normal(loc=0.0, scale=self.information['g_rdnois'+str(k)], size=(4616+self.overscan,1152+self.overscan+self.prescan))
self.image_g[xmin[k]:xmax[k], ymin[k]:ymax[k]]+=noise_g
self.log.info('Bias of %i counts were added to the g band image' % self.information['g_rdnois'+str(k)])
np.random.seed()
noise_r = np.random.normal(loc=0.0, scale=self.information['r_rdnois'+str(k)], size=(4616+self.overscan,1152+self.overscan+self.prescan))
self.image_r[xmin[k]:xmax[k], ymin[k]:ymax[k]]+=noise_r
self.log.info('Bias of %i counts were added to the r band image' % self.information['r_rdnois'+str(k)])
np.random.seed()
noise_i = np.random.normal(loc=0.0, scale=self.information['i_rdnois'+str(k)], size=(4616+self.overscan,1152+self.overscan+self.prescan))
self.image_i[xmin[k]:xmax[k], ymin[k]:ymax[k]]+=noise_i
self.log.info('Bias of %i counts were added to the i band image' % self.information['i_rdnois'+str(k)])
### end for
else:
np.random.seed()
noise = np.random.normal(loc=0.0, scale=4, size=self.image_g.shape)
self.log.info('Sum of readnoise in g channel = %f' % np.sum(noise))
#add to the image
self.image_g += noise
#######################################################
np.random.seed()
noise = np.random.normal(loc=0.0, scale=4, size=self.image_r.shape)
self.log.info('Sum of readnoise in r channel= %f' % np.sum(noise))
#add to the image
self.image_r += noise
########################################
np.random.seed()
noise = np.random.normal(loc=0.0, scale=4, size=self.image_i.shape)
self.log.info('Sum of readnoise in i channel= %f' % np.sum(noise))
#add to the image
self.image_i += noise
return
##########################################################################################################
def electrons2ADU(self):
"""
Convert from electrons to ADUs using the value read from the configuration file.
"""
#####
if self.overscans:
######################
xmin=dict()
xmax=dict()
ymin=dict()
ymax=dict()
############## calculate the non-over subimg of the final matrix of the image
for k in range(1,17):
xmin[k]=0
xmax[k]=4616+self.overscan
ymin[k]=(1152+self.overscan+self.prescan)*(k-1)
ymax[k]=(1152+self.overscan+self.prescan)*(k)
###############
self.image_g[xmin[k]:xmax[k], ymin[k]:ymax[k]]/=self.information['g_gain'+str(k)]
self.image_r[xmin[k]:xmax[k], ymin[k]:ymax[k]]/=self.information['r_gain'+str(k)]
self.image_i[xmin[k]:xmax[k], ymin[k]:ymax[k]]/=self.information['i_gain'+str(k)]
### end for
return
####################################################################################
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).
"""
xmin=dict()
xmax=dict()
ymin=dict()
ymax=dict()
############## calculate the non-over subimg of the final matrix of the image
for k in range(1,17):
xmin[k]=0
xmax[k]=4616+self.overscan
ymin[k]=(1152+self.prescan+self.overscan)*(k-1)
ymax[k]=(1152+self.prescan+self.overscan)*(k)
###############
self.image_g[xmin[k]:xmax[k], ymin[k]:ymax[k]]+=self.information['g_detbia'+str(k)]
self.log.info('Bias of %i counts were added to the g band image' % self.information['g_detbia'+str(k)])
self.image_r[xmin[k]:xmax[k], ymin[k]:ymax[k]]+=self.information['r_detbia'+str(k)]
self.log.info('Bias of %i counts were added to the g band image' % self.information['r_detbia'+str(k)])
self.image_i[xmin[k]:xmax[k], ymin[k]:ymax[k]]+=self.information['i_detbia'+str(k)]
self.log.info('Bias of %i counts were added to the g band image' % self.information['i_detbia'+str(k)])
### end for
return
##############################################################################
####################################################################################
def applyBleeding(self, img, direction='horizon'):
"""
Apply bleeding along the CCD readout direction if the number of electrons in a pixel exceeds the full-well capacity.
direction: this is the CCD readout direction, and default value is horizon direction
:return: the bleeding image
"""
data=img.copy()
if data.max()<self.information['fullwellcapacity']:
return img
else:
if direction=='horizon':
#loop over each column, as bleeding is modelled column-wise
for i, column in enumerate(data): ### select one solumnn
if column.max() <=self.information['fullwellcapacity']:
continue
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
data[i, j] -= overload
sum += overload
elif sum > 0.:
if -overload > sum:
overload = -sum
#self.image[j, i] -= overload
data[i,j] -= overload
sum += overload
######################################################
for i, column in enumerate(data):
if column.max() <=self.information['fullwellcapacity']:
continue
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
data[ i,-j-1] -= overload
sum += overload
elif sum > 0.:
if -overload > sum:
overload = -sum
#self.image[-j-1, i] -= overload
data[i,-j-1,] -= overload
sum += overload
else:
#loop over each column, as bleeding is modelled column-wise
for i, column in enumerate(data.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
data[j, i] -= overload
sum += overload
elif sum > 0.:
if -overload > sum:
overload = -sum
#self.image[j, i] -= overload
data[j, i] -= overload
sum += overload
################################
for i, column in enumerate(data.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
data[-j-1, i] -= overload
sum += overload
elif sum > 0.:
if -overload > sum:
overload = -sum
#self.image[-j-1, i] -= overload
data[-j-1, i] -= overload
sum += overload
######print('Applying column bleeding finished.......')
return data
############################################################################
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_g[self.image_g < 0.0] = 0.
#cut of the values larger than max
self.image_g[self.image_g > max] = max
self.image_g = np.rint(self.image_g).astype(int)
self.log.info('Maximum and total values of the image are %i and %i, respectively' % (np.max(self.image_g),
np.sum(self.image_g)))
#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)))
###############################################################################
#avoid negative numbers in case bias level was not added
self.image_i[self.image_i < 0.0] = 0.
#cut of the values larger than max
self.image_i[self.image_i > max] = max
self.image_i = np.rint(self.image_i).astype(int)
self.log.info('Maximum and total values of the image are %i and %i, respectively' % (np.max(self.image_i),
np.sum(self.image_i)))
return
##################################################################################################
def addScanEffect(self, img):
#### read CCD image to matrix as defined data format
'''
function: add prescan and overscan empty zone to the CCD physical image
return: img
'''
### xsize is column, ysize is row
self.prescan =self.information['prescan'] ## horizon direction , columnn
self.overscan=self.information['overscan'] ## vertical direction, row
row = self.information['ysize']
col = self.information['xsize']
if row !=9232 or col !=9216: #####
print('image size is not correct.......')
sys.exit(1)
####
###########################################
xmin=dict()
xmax=dict()
ymin=dict()
ymax=dict()
#### physical CCD zone 1-16
for n in range(1,17):
#### physical zone index
if n<9:
xmin[n]=0 # row number
xmax[n]=4616
ymin[n]=(n-1)*1152 # col number
ymax[n]=ymin[n]+1152
else:
xmin[n]=4616
xmax[n]=4616*2
ymin[n]=(16-n)*1152
ymax[n]=ymin[n]+1152
###########################################
temp=np.zeros((int(4616+self.overscan), int((1152+self.overscan+self.prescan)*16)))
matxmin=dict()
matxmax=dict()
matymin=dict()
matymax=dict()
############## get the non-over subimg of the final matrix of the image
for k in range(1,17):
### big new image include prescan and overscan;
matxmin[k]=0 ### row number
matxmax[k]=4616 ### row number
matymin[k]=(1152+self.overscan+self.prescan)*(k-1)+self.prescan ## col number start
matymax[k]= matymin[k]+1152 ## col number end
################################
for k in range(1,17):
if k<=4: ### zos1 zos2 zos3 zos4
subimg=img[xmin[k]:xmax[k],ymin[k]:ymax[k]] # get k-th subimg from input img
subimg=np.flipud(subimg) ## down to up, up to down
if k>=5 and k<=8: ### zos5 zos6 zos7 zos8
subimg=img[xmin[k]:xmax[k],ymin[k]:ymax[k]] # get k-th subimg from input img
subimg=np.flipud(subimg) ## down to up, up to down
subimg=np.fliplr(subimg) ## left to right ,right to left
if k>=9 and k<=12: ### zos5 zos6 zos7 zos8
subimg=img[xmin[k]:xmax[k],ymin[k]:ymax[k]] # get k-th subimg from input img
subimg=np.fliplr(subimg) ## left to right ,right to left
if k>=13 and k<=16: ### zos5 zos6 zos7 zos8
subimg=img[xmin[k]:xmax[k],ymin[k]:ymax[k]] # get k-th subimg from input img
######################
temp[matxmin[k]:matxmax[k],matymin[k]:matymax[k]]=subimg[:,:]
#### end for
return temp
###########################################################################
def applyPoissonNoise(self):
"""
Add Poisson noise to the image.
"""
rounded = np.rint(self.image_g) ### round to
residual = self.image_g.copy() - rounded #ugly workaround for multiple rounding operations...
rounded[rounded < 0.0] = 0.0
self.image_g = np.random.poisson(rounded).astype(np.float64)
self.log.info('Added Poisson noise on channel g')
self.image_g += 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
self.image_r = np.random.poisson(rounded).astype(np.float64)
self.log.info('Added Poisson noise on channel r')
self.image_r += residual
rounded = np.rint(self.image_i) ### round to
residual = self.image_i.copy() - rounded #ugly workaround for multiple rounding operations...
rounded[rounded < 0.0] = 0.0
self.image_i = np.random.poisson(rounded).astype(np.float64)
self.log.info('Added Poisson noise on channel i')
self.image_i += residual
#########################################################################################
def writeOutputs(self, obnum=1):
"""
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
##########################################################################
obsid=200000000+obnum
#create a new FITS file, using HDUList instance
ofd_g = fits.PrimaryHDU()
#### add primary header
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")
### data read time is the exposure end time ###
### data read end time is t3
t3=self.dt+timedelta(seconds=self.information['exptime'])+timedelta(seconds=self.information['readouttime'])
filename_g='CSST_MCI_C1_'+self.source+'_'+exp_starttime+'_'+exp_endtime+'_'+str(obsid)+'_'+ \
str(self.filterP[self.filter_g]['number'])+'_L0_'+self.information['ver']+'.fits'
ofd_g.header['GROUPS']=( bool(False), 'always F')
ofd_g.header['DATE'] =( t3.strftime("%Y-%m-%d %H:%M:%S"), 'date this file was written' )
ofd_g.header['FILENAME']=(filename_g, ' file name C48 ')
ofd_g.header['OBSTYPE'] =( self.source, 'observation type raw,flt, mask, bias, dark, sci')
ofd_g.header['TELESCOP']=('CSST', 'always CSST')
ofd_g.header['INSTRUME']=( 'MCI', ' ')
ofd_g.header['RADECSYS']=('ICRS', ' always ICRS ')
ofd_g.header['EQUINOX'] =( float(2000.0), 'always 2000.0')
ofd_g.header['FITSCREA']=( '4.2.1', 'FITS create software version')
######### Object information #############
ofd_g.header['OBJECT']=( self.information['name_obj'], 'object name')
ofd_g.header['TARGET']=( self.information['target'], 'target name, hhmmss+ddmmss')
ofd_g.header['OBJ_RA'] =(np.float64(self.information['ra_obj']) , 'RA of the object in deg')
ofd_g.header['OBJ_DEC']=(np.float64(self.information['dec_obj']) , 'DEC of the object in deg')
ofd_g.header['RA_PNT0']=( np.float64(self.information['ra_pnt0']) , 'RA of the pointing (degrees) at EXPSTART')
ofd_g.header['DEC_PNT0']=(np.float64(self.information['dec_pnt0']) , 'DEC of the pointing (degrees) at EXPSTART')
##############
ofd_g.header['OBSID'] =(str(obsid) , 'observation ID, 3+8bit')
######## Telescope information ###############
# ofd_g.header['COMMENT'] ='=========================================================='
# ofd_g.header['COMMENT'] ='Telescope information'
# ofd_g.header['COMMENT'] ='=========================================================='
ofd_g.header['REFFRAME']=('CSSTGSC-1.0' , 'guide star catalog version')
ofd_g.header['DATE-OBS']=(self.dt.strftime("%Y-%m-%d %H:%M:%S") , 'date of the observation (yyyy-mm-dd hh:mm:ss)')
ofd_g.header['EXPSTART']=(time2jd(self.dt), 'exposure start time, UTC')
ofd_g.header['SUNANGL0']=(np.float32(0.0) , 'angle between sun and optical axis at EXPSTART')
ofd_g.header['MOONANG0']=(np.float32(0.0) , 'angle between moon and optical axis at EXPSTART')
ofd_g.header['POS_ANG0']=(np.float64(0.0), 'angle between optical axis and the North Pole at EXPSTART in arcsec')
ofd_g.header['TEL_ALT0']=(np.float64(0.0), 'angle between optical axis and the ground- piston at EXPSTART in deg')
ofd_g.header['HOODSTA0']=(np.float32(0.0) , 'lens hood altitude at EXPSTART')
ofd_g.header['HOODANG0']=(np.float32(0.0), 'lens hood azimuth at EXPSTART')
ofd_g.header['POSI0_X'] =(np.float64(self.information['pos_x']) , 'the orbital position of CSST in X direction at EXPSTART')
ofd_g.header['POSI0_Y'] =(np.float64(self.information['pos_y']) , 'the orbital position of CSST in Y direction at EXPSTART')
ofd_g.header['POSI0_Z'] =(np.float64(self.information['pos_z']) , 'the orbital position of CSST in Z direction at EXPSTART')
ofd_g.header['VELO0_X']=( np.float64(self.information['velocity_x']) , 'the orbital velocity of CSST in X direction at EXPSTART')
ofd_g.header['VELO0_Y']=( np.float64(self.information['velocity_y']) , 'the orbital velocity of CSST in Y direction at EXPSTART')
ofd_g.header['VELO0_Z']=( np.float64(self.information['velocity_z']) , 'the orbital velocity of CSST in Z direction at EXPSTART')
ofd_g.header['Euler0_1']=( np.float64(0.0), 'Euler angle 1 at EXPSTART')
ofd_g.header['Euler0_2']=( np.float64(0.0), 'Euler angle 2 at EXPSTART')
ofd_g.header['Euler0_3']=( np.float64(0.0), 'Euler angle 3 at EXPSTART')
ofd_g.header['EXPEND'] =(time2jd(t2) , 'exposure end time,UTC')
ofd_g.header['SUNANGL1']=(np.float32(0.0), 'angle between sun and optical axis at EXPEND')
ofd_g.header['MOONANG1']=(np.float32(0.0) , 'angle between moon and optical axis at EXPEND ')
ofd_g.header['POS_ANG1']=(np.float64(0.0) , 'angle between optical axis and the North Pole at EXPEND in arcsec')
ofd_g.header['TEL_ALT1']=(np.float64(0.0) , 'angle between optical axis and the ground- piston at EXPEND in deg ')
ofd_g.header['HOODSTA1']=(np.float32(0.0), 'lens hood altitude at EXPEND ')
ofd_g.header['HOODANG1']=(np.float32(0.0), 'lens hood azimuth at EXPEND ')
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']+'/TianCe/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
ofd_g.header['POSI1_X'] =(np.float64(p1x) , 'the orbital position of CSST in X direction at EXPEND')
ofd_g.header['POSI1_Y'] =(np.float64(p1y) , 'the orbital position of CSST in Y direction at EXPEND')
ofd_g.header['POSI1_Z'] =(np.float64(p1z) , 'the orbital position of CSST in Z direction at EXPEND')
ofd_g.header['VELO1_X']=(np.float64(p1vx) , 'the orbital velocity of CSST in X direction at EXPEND')
ofd_g.header['VELO1_Y']=(np.float64(p1vy) , 'the orbital velocity of CSST in Y direction at EXPEND')
ofd_g.header['VELO1_Z']=(np.float64(p1vz) , 'the orbital velocity of CSST in Z direction at EXPEND')
ofd_g.header['Euler1_1']=( np.float64(0.0), 'Euler angle 1 at EXPEND')
ofd_g.header['Euler1_2']=( np.float64(0.0), 'Euler angle 2 at EXPEND')
ofd_g.header['Euler1_3']=( np.float64(0.0), 'Euler angle 3 at EXPEND')
ofd_g.header['RA_PNT1']=(np.float64(0.0), 'RA of the pointing (degrees) at EXPEND in deg')
ofd_g.header['DEC_PNT1']=(np.float64(0.0), 'DEC of the pointing (degrees) at EXPEND in deg')
ofd_g.header['EXPTIME']=(self.information['exptime'], 'exposure duration')
ofd_g.header['EPOCH'] =(np.float32(0.0), 'coordinate epoch')
ofd_g.header['CHECKSUM']=( 0 , 'hdu-checksum')
########## finish header for 0 layer ####################
#new image HDU
hdu_g = fits.ImageHDU(data=np.uint16(self.image_g))
### MCI header for immage layer ###
#### Instrument information (at EXPSTART time)
hdu_g.header['PMIRRPOS']=( bool(True) , 'FSM pointing,T: point to MCI, F: not point to MCI ')
hdu_g.header['DIFFUSER']=( bool(True) , 'insert diffuser for flat calibration,T: insert diffuser, F: not insert diffuser ')
hdu_g.header['FLAMP'] =( np.int16(0) , 'status of flat lamp,0: off, 1: 1, 2: 2, 3: 3')
hdu_g.header['MCITEMP'] =( np.float32(0.0), 'MCI components temperature')
hdu_g.header['MCISTAT'] =( np.int16(0) , 'MCI components status parameters')
##### Filter information #####
hdu_g.header['FILTER'] =(self.filter_g , 'filter band name')
##### Detector information ####
hdu_g.header['DETNAM'] =('C1' , 'detector name for channle 1')
hdu_g.header['DETTEMP'] =(np.float32(0) , 'detector temperature')
hdu_g.header['DETSIZE'] =(str(self.image_g.shape[0])+'*'+str(self.image_g.shape[1]), 'detector size')
hdu_g.header['DATASEC'] =(str(self.information['ysize'])+'*'+str( self.information['xsize']), 'data section')
hdu_g.header['PIXSCAL1']=(np.float32(0.05) , 'pixel scale for axis 1')
hdu_g.header['PIXSCAL2']=(np.float32(0.05) , 'pixel scale for axis 2')
hdu_g.header['PIXSIZE1']=(np.float32(10.0) , 'pixel size for axis 1')
hdu_g.header['PIXSIZE2']=(np.float32(10.0) , 'pixel size for axis 2')
hdu_g.header['NCHAN'] =(np.int16(16) , 'number of readout channels')
hdu_g.header['NCHAN1'] =(np.int16(8) , 'number of horizontal channels')
hdu_g.header['NCHAN2'] =(np.int16(2) , 'number of vertical channels')
hdu_g.header['PSCAN1'] =(np.int32(27) , 'horizontal prescan width, per readout channel')
hdu_g.header['PSCAN2'] =(np.int32(8) , 'vertical prescan height, per readout channel')
hdu_g.header['OSCAN1'] =(np.int32(16) , 'horizontal overscan width, per readout channel')
hdu_g.header['OSCAN2'] =(np.int32(8) , 'vertical overscan width, per readout channel')
hdu_g.header['BIN_X'] =(np.int16(0) , 'bin number in X (wavelength)')
hdu_g.header['BIN_Y'] =(np.int16(0) , 'bin number in Y (spatial)')
#### World coordinate system and related parameters #####
hdu_g.header['WCSAXES']=(np.int16(2) , 'number of World Coordinate System axes')
hdu_g.header['CRPIX1'] =( np.float64(self.information['CRPIX1']) , 'x-coordinate of reference pixel')
hdu_g.header['CRPIX2'] =( np.float64(self.information['CRPIX2']) , 'y-coordinate of reference pixel')
hdu_g.header['CRVAL1'] =( np.float64(self.information['CRVAL1']) , 'first axis value at reference pixel')
hdu_g.header['CRVAL2'] =( np.float64(self.information['CRVAL2']) , 'second axis value at reference pixel')
hdu_g.header['CTYPE1'] =( 'RA---TAN' , 'the coordinate type for the first axis')
hdu_g.header['CTYPE2'] =( 'DEC---TAN' , 'the coordinate type for the second axis')
hdu_g.header['CD1_1'] =( np.float64(self.information['CD1_1']) , 'partial of first axis coordinate w.r.t. x')
hdu_g.header['CD1_2'] =( np.float64(self.information['CD1_2']) , 'partial of first axis coordinate w.r.t. y')
hdu_g.header['CD2_1'] =( np.float64(self.information['CD2_1']) , 'partial of second axis coordinate w.r.t. x')
hdu_g.header['CD2_2'] =( np.float64(self.information['CD2_2']) , 'partial of second axis coordinate w.r.t. y')
hdu_g.header['others'] =( 0.0 , 'other specific parameters')
###### Readout information #######
hdu_g.header['GAIN1'] =( np.float32(self.information['g_gain1']) , 'CCD gain (channel 1)')
hdu_g.header['GAIN2'] =( np.float32(self.information['g_gain2']) , 'CCD gain (channel 2)')
hdu_g.header['GAIN3'] =( np.float32(self.information['g_gain3']) , 'CCD gain (channel 3)')
hdu_g.header['GAIN4'] =( np.float32(self.information['g_gain4']) , 'CCD gain (channel 4)')
hdu_g.header['GAIN5'] =( np.float32(self.information['g_gain5']) , 'CCD gain (channel 5)')
hdu_g.header['GAIN6'] =( np.float32(self.information['g_gain6']) , 'CCD gain (channel 6)')
hdu_g.header['GAIN7'] =( np.float32(self.information['g_gain7']) , 'CCD gain (channel 7)')
hdu_g.header['GAIN8'] =( np.float32(self.information['g_gain8']) , 'CCD gain (channel 8)')
hdu_g.header['GAIN9'] =( np.float32(self.information['g_gain9']) , 'CCD gain (channel 9)')
hdu_g.header['GAIN10']=( np.float32(self.information['g_gain10']) , 'CCD gain (channel 10)')
hdu_g.header['GAIN11']=( np.float32(self.information['g_gain11']) , 'CCD gain (channel 11)')
hdu_g.header['GAIN12']=( np.float32(self.information['g_gain12']) , 'CCD gain (channel 12)')
hdu_g.header['GAIN13']=( np.float32(self.information['g_gain13']) , 'CCD gain (channel 13)')
hdu_g.header['GAIN14']=( np.float32(self.information['g_gain14']) , 'CCD gain (channel 14)')
hdu_g.header['GAIN15']=( np.float32(self.information['g_gain15']) , 'CCD gain (channel 15)')
hdu_g.header['GAIN16']=( np.float32(self.information['g_gain16']) , 'CCD gain (channel 16)')
#######
hdu_g.header['RDNOIS1'] =( np.float32(self.information['g_rdnois1']) , 'read noise (channel 1)')
hdu_g.header['RDNOIS2'] =( np.float32(self.information['g_rdnois2']) , 'read noise (channel 2)')
hdu_g.header['RDNOIS3'] =( np.float32(self.information['g_rdnois3']) , 'read noise (channel 3)')
hdu_g.header['RDNOIS4'] =( np.float32(self.information['g_rdnois4']) , 'read noise (channel 4)')
hdu_g.header['RDNOIS5'] =( np.float32(self.information['g_rdnois5']) , 'read noise (channel 5)')
hdu_g.header['RDNOIS6'] =( np.float32(self.information['g_rdnois6']) , 'read noise (channel 6)')
hdu_g.header['RDNOIS7'] =( np.float32(self.information['g_rdnois7']) , 'read noise (channel 7)')
hdu_g.header['RDNOIS8'] =( np.float32(self.information['g_rdnois8']) , 'read noise (channel 8)')
hdu_g.header['RDNOIS9'] =( np.float32(self.information['g_rdnois9']) , 'read noise (channel 9)')
hdu_g.header['RDNOIS10']=( np.float32(self.information['g_rdnois10']) , 'read noise (channel 10)')
hdu_g.header['RDNOIS11']=( np.float32(self.information['g_rdnois11']) , 'read noise (channel 11)')
hdu_g.header['RDNOIS12']=( np.float32(self.information['g_rdnois12']) , 'read noise (channel 12)')
hdu_g.header['RDNOIS13']=( np.float32(self.information['g_rdnois13']) , 'read noise (channel 13)')
hdu_g.header['RDNOIS14']=( np.float32(self.information['g_rdnois14']) , 'read noise (channel 14)')
hdu_g.header['RDNOIS15']=( np.float32(self.information['g_rdnois15']) , 'read noise (channel 15)')
hdu_g.header['RDNOIS16']=( np.float32(self.information['g_rdnois16']) , 'read noise (channel 16)')
#######
hdu_g.header['DETBIA1']=( np.float32(self.information['g_detbia1']) , 'amplifier bias voltage (channel 1)')
hdu_g.header['DETBIA2']=( np.float32(self.information['g_detbia2']) , 'amplifier bias voltage (channel 2)')
hdu_g.header['DETBIA3']=( np.float32(self.information['g_detbia3']) , 'amplifier bias voltage (channel 3)')
hdu_g.header['DETBIA4']=( np.float32(self.information['g_detbia4']) , 'amplifier bias voltage (channel 4)')
hdu_g.header['DETBIA5']=( np.float32(self.information['g_detbia5']) , 'amplifier bias voltage (channel 5)')
hdu_g.header['DETBIA6']=( np.float32(self.information['g_detbia6']) , 'amplifier bias voltage (channel 6)')
hdu_g.header['DETBIA7']=( np.float32(self.information['g_detbia7']) , 'amplifier bias voltage (channel 7)')
hdu_g.header['DETBIA8']=( np.float32(self.information['g_detbia8']) , 'amplifier bias voltage (channel 8)')
hdu_g.header['DETBIA9']=( np.float32(self.information['g_detbia9']) , 'amplifier bias voltage (channel 9)')
hdu_g.header['DETBIA10']=( np.float32(self.information['g_detbia10']) , 'amplifier bias voltage (channel 10)')
hdu_g.header['DETBIA11']=( np.float32(self.information['g_detbia11']) , 'amplifier bias voltage (channel 11)')
hdu_g.header['DETBIA12']=( np.float32(self.information['g_detbia12']) , 'amplifier bias voltage (channel 12)')
hdu_g.header['DETBIA13']=( np.float32(self.information['g_detbia13']) , 'amplifier bias voltage (channel 13)')
hdu_g.header['DETBIA14']=( np.float32(self.information['g_detbia14']) , 'amplifier bias voltage (channel 14)')
hdu_g.header['DETBIA15']=( np.float32(self.information['g_detbia15']) , 'amplifier bias voltage (channel 15)')
hdu_g.header['DETBIA16']=( np.float32(self.information['g_detbia16']) , 'amplifier bias voltage (channel 16)')
######
hdu_g.header['READT0'] =( time2jd(t2) , 'read start time (UTC)')
hdu_g.header['READT1'] =( time2jd(t3) , 'read end time (UTC)')
hdu_g.header['DETTEMP0']=( np.float32(0) , 'detector temperature at READT0')
hdu_g.header['DETTEMP1']=( np.float32(0) , 'detector temperature at READT1')
hdu_g.header['RDSPEED'] =( np.float32(3) , 'read speed (in MHz)')
##### Shutter information #####
hdu_g.header['SHUTHWV'] =( 'XXXXXXXXXXXXXXXX' , 'shutter hardware version')
hdu_g.header['SHUTSWV'] =( 'XXXXXXXXXXXXXXXX' , 'shutter software')
hdu_g.header['SHUTSTAT']=( bool(True) , 'shutter status,T: open, F: close')
hdu_g.header['SHTOPEN0']=( np.float64(0) , 'shutter open time (begin)')
hdu_g.header['SHTOPEN1']=( np.float64(0) , 'shutter open time (end)')
hdu_g.header['SHTCLOS0']=( np.float64(0) , 'shutter close time (begin)')
hdu_g.header['SHTCLOS1']=( np.float64(0) , 'shutter close time (end)')
################ finish MCI header for image layer
##############################################################################################333
hdulist_g=fits.HDUList([ofd_g, hdu_g])
if self.source=='DARK' or self.source=='FLAT' or self.source=='BIAS':
file_g=self.result_path+'/cali_Data/'+filename_g
else:
file_g=self.result_path+'/sky_Data/'+filename_g
hdulist_g.writeto(file_g, overwrite=True)
###########################################################################
### r band
#create a new FITS file, using HDUList instance
ofd_r = fits.PrimaryHDU()
filename_r='CSST_MCI_C2_'+self.source+'_'+exp_starttime+'_'+exp_endtime+'_'+str(obsid)+'_'+ \
str(self.filterP[self.filter_r]['number'])+'_L0_'+self.information['ver']+'.fits'
######### add primary header ##########
ofd_r.header['GROUPS']=( bool(False), 'always F')
ofd_r.header['DATE'] =( ofd_g.header['DATE'], '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']=( 'MCI', ' ')
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 #############
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']=(ofd_g.header['DATE-OBS'] , 'date of the observation (yyyy-mm-dd hh:mm:ss)')
ofd_r.header['EXPSTART']=(ofd_g.header['EXPSTART'], '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['pos_x']) , 'the orbital position of CSST in X direction at EXPSTART')
ofd_r.header['POSI0_Y'] =(np.float64(self.information['pos_y']) , 'the orbital position of CSST in Y direction at EXPSTART')
ofd_r.header['POSI0_Z'] =(np.float64(self.information['pos_z']) , 'the orbital position of CSST in Z direction at EXPSTART')
ofd_r.header['VELO0_X']=( np.float64(self.information['velocity_x']) , 'the orbital velocity of CSST in X direction at EXPSTART')
ofd_r.header['VELO0_Y']=( np.float64(self.information['velocity_y']) , 'the orbital velocity of CSST in Y direction at EXPSTART')
ofd_r.header['VELO0_Z']=( np.float64(self.information['velocity_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'] =(ofd_g.header['EXPEND'] , '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(p1x) , 'the orbital position of CSST in X direction at EXPEND')
ofd_r.header['POSI1_Y'] =(np.float64(p1y) , 'the orbital position of CSST in Y direction at EXPEND')
ofd_r.header['POSI1_Z'] =(np.float64(p1z) , 'the orbital position of CSST in Z direction at EXPEND')
ofd_r.header['VELO1_X']=(np.float64(p1vx) , 'the orbital velocity of CSST in X direction at EXPEND')
ofd_r.header['VELO1_Y']=(np.float64(p1vy) , 'the orbital velocity of CSST in Y direction at EXPEND')
ofd_r.header['VELO1_Z']=(np.float64(p1vz) , '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(0.0), 'RA of the pointing (degrees) at EXPEND in deg')
ofd_r.header['DEC_PNT1']=(np.float64(0.0), '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 header for 0 layer
#new image HDU
hdu_r = fits.ImageHDU(data=np.uint16(self.image_r))
### MCI header for immage layer ###
#### Instrument information (at EXPSTART time)
hdu_r.header['PMIRRPOS']=( bool(True) , 'FSM pointing,T: point to MCI, F: not point to MCI ')
hdu_r.header['DIFFUSER']=( bool(True) , 'insert diffuser for flat calibration,T: insert diffuser, F: not insert diffuser ')
hdu_r.header['FLAMP'] =( np.int16(0) , 'status of flat lamp,0: off, 1: 1, 2: 2, 3: 3')
hdu_r.header['MCITEMP'] =( np.float32(0.0), 'MCI components temperature')
hdu_r.header['MCISTAT'] =( np.int16(0) , 'MCI components status parameters')
##### Filter information #####
hdu_r.header['FILTER'] =(self.filter_r , 'filter band name')
##### Detector information ####
hdu_r.header['DETNAM'] =('C2' , 'detector name for channle 2')
hdu_r.header['DETTEMP'] =(np.float32(0) , 'detector temperature')
hdu_r.header['DETSIZE'] =(str(self.image_r.shape[0])+'*'+str(self.image_r.shape[1]), 'detector size')
hdu_r.header['DATASEC'] =(str(self.information['ysize'])+'*'+str( self.information['xsize']), 'data section')
hdu_r.header['PIXSCAL1']=(np.float32(0.05) , 'pixel scale for axis 1')
hdu_r.header['PIXSCAL2']=(np.float32(0.05) , 'pixel scale for axis 2')
hdu_r.header['PIXSIZE1']=(np.float32(10.0) , 'pixel size for axis 1')
hdu_r.header['PIXSIZE2']=(np.float32(10.0) , 'pixel size for axis 2')
hdu_r.header['NCHAN'] =(np.int16(16) , 'number of readout channels')
hdu_r.header['NCHAN1'] =(np.int16(8) , 'number of horizontal channels')
hdu_r.header['NCHAN2'] =(np.int16(2) , 'number of vertical channels')
hdu_r.header['PSCAN1'] =(np.int32(27) , 'horizontal prescan width, per readout channel')
hdu_r.header['PSCAN2'] =(np.int32(8) , 'vertical prescan height, per readout channel')
hdu_r.header['OSCAN1'] =(np.int32(16) , 'horizontal overscan width, per readout channel')
hdu_r.header['OSCAN2'] =(np.int32(8) , 'vertical overscan width, per readout channel')
hdu_r.header['BIN_X'] =(np.int16(0) , 'bin number in X (wavelength)')
hdu_r.header['BIN_Y'] =(np.int16(0) , 'bin number in Y (spatial)')
#### World coordinate system and related parameters #####
hdu_r.header['WCSAXES']=( np.int16(2) , 'number of World Coordinate System axes')
hdu_r.header['CRPIX1'] =( np.float64(self.information['CRPIX1']) , 'x-coordinate of reference pixel')
hdu_r.header['CRPIX2'] =( np.float64(self.information['CRPIX2']) , 'y-coordinate of reference pixel')
hdu_r.header['CRVAL1'] =( np.float64(self.information['CRVAL1']) , 'first axis value at reference pixel')
hdu_r.header['CRVAL2'] =( np.float64(self.information['CRVAL2']) , 'second axis value at reference pixel')
hdu_r.header['CTYPE1'] =( 'RA---TAN' , 'the coordinate type for the first axis')
hdu_r.header['CTYPE2'] =( 'DEC---TAN' , 'the coordinate type for the second axis')
hdu_r.header['CD1_1'] =( np.float64(self.information['CD1_1']) , 'partial of first axis coordinate w.r.t. x')
hdu_r.header['CD1_2'] =( np.float64(self.information['CD1_2']) , 'partial of first axis coordinate w.r.t. y')
hdu_r.header['CD2_1'] =( np.float64(self.information['CD2_1']) , 'partial of second axis coordinate w.r.t. x')
hdu_r.header['CD2_2'] =( np.float64(self.information['CD2_2']) , 'partial of second axis coordinate w.r.t. y')
hdu_r.header['others'] =( 0.0 , 'other specific parameters')
###### Readout information #######
hdu_r.header['GAIN1'] =( np.float32(self.information['r_gain1']) , 'CCD gain (channel 1)')
hdu_r.header['GAIN2'] =( np.float32(self.information['r_gain2']) , 'CCD gain (channel 2)')
hdu_r.header['GAIN3'] =( np.float32(self.information['r_gain3']) , 'CCD gain (channel 3)')
hdu_r.header['GAIN4'] =( np.float32(self.information['r_gain4']) , 'CCD gain (channel 4)')
hdu_r.header['GAIN5'] =( np.float32(self.information['r_gain5']) , 'CCD gain (channel 5)')
hdu_r.header['GAIN6'] =( np.float32(self.information['r_gain6']) , 'CCD gain (channel 6)')
hdu_r.header['GAIN7'] =( np.float32(self.information['r_gain7']) , 'CCD gain (channel 7)')
hdu_r.header['GAIN8'] =( np.float32(self.information['r_gain8']) , 'CCD gain (channel 8)')
hdu_r.header['GAIN9'] =( np.float32(self.information['r_gain9']) , 'CCD gain (channel 9)')
hdu_r.header['GAIN10']=( np.float32(self.information['r_gain10']) , 'CCD gain (channel 10)')
hdu_r.header['GAIN11']=( np.float32(self.information['r_gain11']) , 'CCD gain (channel 11)')
hdu_r.header['GAIN12']=( np.float32(self.information['r_gain12']) , 'CCD gain (channel 12)')
hdu_r.header['GAIN13']=( np.float32(self.information['r_gain13']) , 'CCD gain (channel 13)')
hdu_r.header['GAIN14']=( np.float32(self.information['r_gain14']) , 'CCD gain (channel 14)')
hdu_r.header['GAIN15']=( np.float32(self.information['r_gain15']) , 'CCD gain (channel 15)')
hdu_r.header['GAIN16']=( np.float32(self.information['r_gain16']) , 'CCD gain (channel 16)')
#######
hdu_r.header['RDNOIS1'] =( np.float32(self.information['r_rdnois1']) , 'read noise (channel 1)')
hdu_r.header['RDNOIS2'] =( np.float32(self.information['r_rdnois2']) , 'read noise (channel 2)')
hdu_r.header['RDNOIS3'] =( np.float32(self.information['r_rdnois3']) , 'read noise (channel 3)')
hdu_r.header['RDNOIS4'] =( np.float32(self.information['r_rdnois4']) , 'read noise (channel 4)')
hdu_r.header['RDNOIS5'] =( np.float32(self.information['r_rdnois5']) , 'read noise (channel 5)')
hdu_r.header['RDNOIS6'] =( np.float32(self.information['r_rdnois6']) , 'read noise (channel 6)')
hdu_r.header['RDNOIS7'] =( np.float32(self.information['r_rdnois7']) , 'read noise (channel 7)')
hdu_r.header['RDNOIS8'] =( np.float32(self.information['r_rdnois8']) , 'read noise (channel 8)')
hdu_r.header['RDNOIS9'] =( np.float32(self.information['r_rdnois9']) , 'read noise (channel 9)')
hdu_r.header['RDNOIS10']=( np.float32(self.information['r_rdnois10']) , 'read noise (channel 10)')
hdu_r.header['RDNOIS11']=( np.float32(self.information['r_rdnois11']) , 'read noise (channel 11)')
hdu_r.header['RDNOIS12']=( np.float32(self.information['r_rdnois12']) , 'read noise (channel 12)')
hdu_r.header['RDNOIS13']=( np.float32(self.information['r_rdnois13']) , 'read noise (channel 13)')
hdu_r.header['RDNOIS14']=( np.float32(self.information['r_rdnois14']) , 'read noise (channel 14)')
hdu_r.header['RDNOIS15']=( np.float32(self.information['r_rdnois15']) , 'read noise (channel 15)')
hdu_r.header['RDNOIS16']=( np.float32(self.information['r_rdnois16']) , 'read noise (channel 16)')
#######
hdu_r.header['DETBIA1']=( np.float32(self.information['r_detbia1']) , 'amplifier bias voltage (channel 1)')
hdu_r.header['DETBIA2']=( np.float32(self.information['r_detbia2']) , 'amplifier bias voltage (channel 2)')
hdu_r.header['DETBIA3']=( np.float32(self.information['r_detbia3']) , 'amplifier bias voltage (channel 3)')
hdu_r.header['DETBIA4']=( np.float32(self.information['r_detbia4']) , 'amplifier bias voltage (channel 4)')
hdu_r.header['DETBIA5']=( np.float32(self.information['r_detbia5']) , 'amplifier bias voltage (channel 5)')
hdu_r.header['DETBIA6']=( np.float32(self.information['r_detbia6']) , 'amplifier bias voltage (channel 6)')
hdu_r.header['DETBIA7']=( np.float32(self.information['r_detbia7']) , 'amplifier bias voltage (channel 7)')
hdu_r.header['DETBIA8']=( np.float32(self.information['r_detbia8']) , 'amplifier bias voltage (channel 8)')
hdu_r.header['DETBIA9']=( np.float32(self.information['r_detbia9']) , 'amplifier bias voltage (channel 9)')
hdu_r.header['DETBIA10']=( np.float32(self.information['r_detbia10']) , 'amplifier bias voltage (channel 10)')
hdu_r.header['DETBIA11']=( np.float32(self.information['r_detbia11']) , 'amplifier bias voltage (channel 11)')
hdu_r.header['DETBIA12']=( np.float32(self.information['r_detbia12']) , 'amplifier bias voltage (channel 12)')
hdu_r.header['DETBIA13']=( np.float32(self.information['r_detbia13']) , 'amplifier bias voltage (channel 13)')
hdu_r.header['DETBIA14']=( np.float32(self.information['r_detbia14']) , 'amplifier bias voltage (channel 14)')
hdu_r.header['DETBIA15']=( np.float32(self.information['r_detbia15']) , 'amplifier bias voltage (channel 15)')
hdu_r.header['DETBIA16']=( np.float32(self.information['r_detbia16']) , 'amplifier bias voltage (channel 16)')
######
######
hdu_r.header['READT0'] =( hdu_g.header['READT0'] , 'read start time (UTC)')
hdu_r.header['READT1'] =( hdu_g.header['READT1'] , 'read end time (UTC)')
hdu_r.header['DETTEMP0']=( np.float32(0) , 'detector temperature at READT0')
hdu_r.header['DETTEMP1']=( np.float32(0) , 'detector temperature at READT1')
hdu_r.header['RDSPEED'] =( np.float32(3) , 'read speed (in MHz)')
##### Shutter information #####
hdu_r.header['SHUTHWV'] =( 'XXXXXXXXXXXXXXXX' , 'shutter hardware version')
hdu_r.header['SHUTSWV'] =( 'XXXXXXXXXXXXXXXX' , 'shutter software')
hdu_r.header['SHUTSTAT']=( bool(True) , 'shutter status,T: open, F: close')
hdu_r.header['SHTOPEN0']=( np.float64(0) , 'shutter open time (begin)')
hdu_r.header['SHTOPEN1']=( np.float64(0) , 'shutter open time (end)')
hdu_r.header['SHTCLOS0']=( np.float64(0) , 'shutter close time (begin)')
hdu_r.header['SHTCLOS1']=( np.float64(0) , 'shutter close time (end)')
################ finish MCI header for image layer
#write the actual file
hdulist_r=fits.HDUList([ofd_r, hdu_r])
#file_r=self.result_path+'/sky_Data/'+filename_r
if self.source=='DARK' or self.source=='FLAT' or self.source=='BIAS':
file_r=self.result_path+'/cali_Data/'+filename_r
else:
file_r=self.result_path+'/sky_Data/'+filename_r
hdulist_r.writeto(file_r, overwrite=True)
#################################################################################
### i band
#create a new FITS file, using HDUList instance
ofd_i = fits.PrimaryHDU()
filename_i='CSST_MCI_C3_'+self.source+'_'+exp_starttime+'_'+exp_endtime+'_'+str(obsid)+'_'+ \
str(self.filterP[self.filter_i]['number'])+'_L0_'+self.information['ver']+'.fits'
#### add primary header
ofd_i.header['GROUPS']=( bool(False), 'always F')
ofd_i.header['DATE'] =( ofd_g.header['DATE'], 'date this file was written' )
ofd_i.header['FILENAME']=(filename_i, ' file name C48 ')
ofd_i.header['OBSTYPE'] =( self.source, 'observation type raw,flt, mask, bias, dark, sci')
ofd_i.header['TELESCOP']=('CSST', 'always CSST')
ofd_i.header['INSTRUME']=( 'MCI', ' ')
ofd_i.header['RADECSYS']=('ICRS', ' always ICRS ')
ofd_i.header['EQUINOX'] =( float(2000.0), 'always 2000.0')
ofd_i.header['FITSCREA']=( '4.2.1', 'FITS create software version')
######### Object information #############
ofd_i.header['OBJECT']=( self.information['name_obj'], 'object name')
ofd_i.header['TARGET']=( self.information['target'], 'target name, hhmmss+ddmmss')
ofd_i.header['OBJ_RA'] =(np.float64(self.information['ra_obj']) , 'RA of the object in deg')
ofd_i.header['OBJ_DEC']=(np.float64(self.information['dec_obj']) , 'DEC of the object in deg')
ofd_i.header['RA_PNT0'] =( np.float64(self.information['ra_pnt0']) , 'RA of the pointing (degrees) at EXPSTART')
ofd_i.header['DEC_PNT0']=(np.float64(self.information['dec_pnt0']) , 'DEC of the pointing (degrees) at EXPSTART')
##############
ofd_i.header['OBSID'] =(str(obsid) , 'observation ID, 3+8bit')
######## Telescope information ###############
# ofd_i.header['COMMENT'] ='=========================================================='
# ofd_i.header['COMMENT'] ='Telescope information'
# ofd_i.header['COMMENT'] ='=========================================================='
ofd_i.header['REFFRAME']=('CSSTGSC-1.0' , 'guide star catalog version')
ofd_i.header['DATE-OBS']=(ofd_g.header['DATE-OBS'] , 'date of the observation (yyyy-mm-dd hh:mm:ss)')
ofd_i.header['EXPSTART']=(ofd_g.header['EXPSTART'], 'exposure start time')
ofd_i.header['SUNANGL0']=(np.float32(0.0) , 'angle between sun and optical axis at EXPSTART')
ofd_i.header['MOONANG0']=(np.float32(0.0) , 'angle between moon and optical axis at EXPSTART')
ofd_i.header['POS_ANG0']=(np.float64(0.0), 'angle between optical axis and the North Pole at EXPSTART in arcsec')
ofd_i.header['TEL_ALT0']=(np.float64(0.0), 'angle between optical axis and the ground- piston at EXPSTART in deg')
ofd_i.header['HOODSTA0']=(np.float32(0.0) , 'lens hood altitude at EXPSTART')
ofd_i.header['HOODANG0']=(np.float32(0.0), 'lens hood azimuth at EXPSTART')
ofd_i.header['POSI0_X'] =(np.float64(self.information['pos_x']) , 'the orbital position of CSST in X direction at EXPSTART')
ofd_i.header['POSI0_Y'] =(np.float64(self.information['pos_y']) , 'the orbital position of CSST in Y direction at EXPSTART')
ofd_i.header['POSI0_Z'] =(np.float64(self.information['pos_z']) , 'the orbital position of CSST in Z direction at EXPSTART')
ofd_i.header['VELO0_X']=( np.float64(self.information['velocity_x']) , 'the orbital velocity of CSST in X direction at EXPSTART')
ofd_i.header['VELO0_Y']=( np.float64(self.information['velocity_y']) , 'the orbital velocity of CSST in Y direction at EXPSTART')
ofd_i.header['VELO0_Z']=( np.float64(self.information['velocity_z']) , 'the orbital velocity of CSST in Z direction at EXPSTART')
ofd_i.header['Euler0_1']=( np.float64(0.0), 'Euler angle 1 at EXPSTART')
ofd_i.header['Euler0_2']=( np.float64(0.0), 'Euler angle 2 at EXPSTART')
ofd_i.header['Euler0_3']=( np.float64(0.0), 'Euler angle 3 at EXPSTART')
ofd_i.header['EXPEND'] =(ofd_g.header['EXPEND'] , 'exposure end time')
ofd_i.header['SUNANGL1']=(np.float32(0.0), 'angle between sun and optical axis at EXPEND')
ofd_i.header['MOONANG1']=(np.float32(0.0) , 'angle between moon and optical axis at EXPEND ')
ofd_i.header['POS_ANG1']=(np.float64(0.0) , 'angle between optical axis and the North Pole at EXPEND in arcsec')
ofd_i.header['TEL_ALT1']=(np.float64(0.0) , 'angle between optical axis and the ground- piston at EXPEND in deg ')
ofd_i.header['HOODSTA1']=(np.float32(0.0), 'lens hood altitude at EXPEND ')
ofd_i.header['HOODANG1']=(np.float32(0.0), 'lens hood azimuth at EXPEND ')
ofd_i.header['POSI1_X'] =(np.float64(p1x) , 'the orbital position of CSST in X direction at EXPEND')
ofd_i.header['POSI1_Y'] =(np.float64(p1y) , 'the orbital position of CSST in Y direction at EXPEND')
ofd_i.header['POSI1_Z'] =(np.float64(p1z) , 'the orbital position of CSST in Z direction at EXPEND')
ofd_i.header['VELO1_X']=(np.float64(p1vx) , 'the orbital velocity of CSST in X direction at EXPEND')
ofd_i.header['VELO1_Y']=(np.float64(p1vy) , 'the orbital velocity of CSST in Y direction at EXPEND')
ofd_i.header['VELO1_Z']=(np.float64(p1vz) , 'the orbital velocity of CSST in Z direction at EXPEND')
ofd_i.header['Euler1_1']=( np.float64(0.0), 'Euler angle 1 at EXPEND')
ofd_i.header['Euler1_2']=( np.float64(0.0), 'Euler angle 2 at EXPEND')
ofd_i.header['Euler1_3']=( np.float64(0.0), 'Euler angle 3 at EXPEND')
ofd_i.header['RA_PNT1']=(np.float64(0.0), 'RA of the pointing (degrees) at EXPEND in deg')
ofd_i.header['DEC_PNT1']=(np.float64(0.0), 'DEC of the pointing (degrees) at EXPEND in deg')
ofd_i.header['EXPTIME']=(self.information['exptime'], 'exposure duration')
ofd_i.header['EPOCH'] =(np.float32(0.0), 'coordinate epoch')
ofd_i.header['CHECKSUM']=( 0 , 'hdu-checksum')
########## finish header for 0 layer
#new image HDU
hdu_i = fits.ImageHDU(data=np.uint16(self.image_i))
### MCI header for immage layer ###
#### Instrument information (at EXPSTART time)
hdu_i.header['PMIRRPOS']=( bool(True) , 'FSM pointing,T: point to MCI, F: not point to MCI ')
hdu_i.header['DIFFUSER']=( bool(True) , 'insert diffuser for flat calibration,T: insert diffuser, F: not insert diffuser ')
hdu_i.header['FLAMP'] =( np.int16(0) , 'status of flat lamp,0: off, 1: 1, 2: 2, 3: 3')
hdu_i.header['MCITEMP'] =( np.float32(0.0), 'MCI components temperature')
hdu_i.header['MCISTAT'] =( np.int16(0) , 'MCI components status parameters')
##### Filter information #####
hdu_i.header['FILTER'] =(self.filter_i , 'filter band name')
##### Detector information ####
hdu_i.header['DETNAM'] =('C3' , 'detector name for channle 3')
hdu_i.header['DETTEMP'] =(np.float32(0) , 'detector temperature')
hdu_i.header['DETSIZE'] =(str(self.image_i.shape[0])+'*'+str(self.image_i.shape[1]), 'detector size')
hdu_i.header['DATASEC'] =(str(self.information['ysize'])+'*'+str( self.information['xsize']), 'data section')
hdu_i.header['PIXSCAL1']=(np.float32(0.05) , 'pixel scale for axis 1')
hdu_i.header['PIXSCAL2']=(np.float32(0.05) , 'pixel scale for axis 2')
hdu_i.header['PIXSIZE1']=(np.float32(10.0) , 'pixel size for axis 1')
hdu_i.header['PIXSIZE2']=(np.float32(10.0) , 'pixel size for axis 2')
hdu_i.header['NCHAN'] =(np.int16(16) , 'number of readout channels')
hdu_i.header['NCHAN1'] =(np.int16(8) , 'number of horizontal channels')
hdu_i.header['NCHAN2'] =(np.int16(2) , 'number of vertical channels')
hdu_i.header['PSCAN1'] =(np.int32(27) , 'horizontal prescan width, per readout channel')
hdu_i.header['PSCAN2'] =(np.int32(8) , 'vertical prescan height, per readout channel')
hdu_i.header['OSCAN1'] =(np.int32(16) , 'horizontal overscan width, per readout channel')
hdu_i.header['OSCAN2'] =(np.int32(8) , 'vertical overscan width, per readout channel')
hdu_i.header['BIN_X'] =(np.int16(0) , 'bin number in X (wavelength)')
hdu_i.header['BIN_Y'] =(np.int16(0) , 'bin number in Y (spatial)')
#### World coordinate system and related parameters #####
hdu_i.header['WCSAXES']=( np.int16(2) , 'number of World Coordinate System axes')
hdu_i.header['CRPIX1'] =( np.float64(self.information['CRPIX1']) , 'x-coordinate of reference pixel')
hdu_i.header['CRPIX2'] =( np.float64(self.information['CRPIX2']) , 'y-coordinate of reference pixel')
hdu_i.header['CRVAL1'] =( np.float64(self.information['CRVAL1']) , 'first axis value at reference pixel')
hdu_i.header['CRVAL2'] =( np.float64(self.information['CRVAL2']) , 'second axis value at reference pixel')
hdu_i.header['CTYPE1'] =( 'RA---TAN' , 'the coordinate type for the first axis')
hdu_i.header['CTYPE2'] =( 'DEC---TAN' , 'the coordinate type for the second axis')
hdu_i.header['CD1_1'] =( np.float64(self.information['CD1_1']) , 'partial of first axis coordinate w.r.t. x')
hdu_i.header['CD1_2'] =( np.float64(self.information['CD1_2']) , 'partial of first axis coordinate w.r.t. y')
hdu_i.header['CD2_1'] =( np.float64(self.information['CD2_1']) , 'partial of second axis coordinate w.r.t. x')
hdu_i.header['CD2_2'] =( np.float64(self.information['CD2_2']) , 'partial of second axis coordinate w.r.t. y')
hdu_i.header['others'] =( 0.0 , 'other specific parameters')
###### Readout information #######
hdu_i.header['GAIN1'] =( np.float32(self.information['i_gain1']) , 'CCD gain (channel 1)')
hdu_i.header['GAIN2'] =( np.float32(self.information['i_gain2']) , 'CCD gain (channel 2)')
hdu_i.header['GAIN3'] =( np.float32(self.information['i_gain3']) , 'CCD gain (channel 3)')
hdu_i.header['GAIN4'] =( np.float32(self.information['i_gain4']) , 'CCD gain (channel 4)')
hdu_i.header['GAIN5'] =( np.float32(self.information['i_gain5']) , 'CCD gain (channel 5)')
hdu_i.header['GAIN6'] =( np.float32(self.information['i_gain6']) , 'CCD gain (channel 6)')
hdu_i.header['GAIN7'] =( np.float32(self.information['i_gain7']) , 'CCD gain (channel 7)')
hdu_i.header['GAIN8'] =( np.float32(self.information['i_gain8']) , 'CCD gain (channel 8)')
hdu_i.header['GAIN9'] =( np.float32(self.information['i_gain9']) , 'CCD gain (channel 9)')
hdu_i.header['GAIN10']=( np.float32(self.information['i_gain10']) , 'CCD gain (channel 10)')
hdu_i.header['GAIN11']=( np.float32(self.information['i_gain11']) , 'CCD gain (channel 11)')
hdu_i.header['GAIN12']=( np.float32(self.information['i_gain12']) , 'CCD gain (channel 12)')
hdu_i.header['GAIN13']=( np.float32(self.information['i_gain13']) , 'CCD gain (channel 13)')
hdu_i.header['GAIN14']=( np.float32(self.information['i_gain14']) , 'CCD gain (channel 14)')
hdu_i.header['GAIN15']=( np.float32(self.information['i_gain15']) , 'CCD gain (channel 15)')
hdu_i.header['GAIN16']=( np.float32(self.information['i_gain16']) , 'CCD gain (channel 16)')
#######
hdu_i.header['RDNOIS1'] =( np.float32(self.information['i_rdnois1']) , 'read noise (channel 1)')
hdu_i.header['RDNOIS2'] =( np.float32(self.information['i_rdnois2']) , 'read noise (channel 2)')
hdu_i.header['RDNOIS3'] =( np.float32(self.information['i_rdnois3']) , 'read noise (channel 3)')
hdu_i.header['RDNOIS4'] =( np.float32(self.information['i_rdnois4']) , 'read noise (channel 4)')
hdu_i.header['RDNOIS5'] =( np.float32(self.information['i_rdnois5']) , 'read noise (channel 5)')
hdu_i.header['RDNOIS6'] =( np.float32(self.information['i_rdnois6']) , 'read noise (channel 6)')
hdu_i.header['RDNOIS7'] =( np.float32(self.information['i_rdnois7']) , 'read noise (channel 7)')
hdu_i.header['RDNOIS8'] =( np.float32(self.information['i_rdnois8']) , 'read noise (channel 8)')
hdu_i.header['RDNOIS9'] =( np.float32(self.information['i_rdnois9']) , 'read noise (channel 9)')
hdu_i.header['RDNOIS10']=( np.float32(self.information['i_rdnois10']) , 'read noise (channel 10)')
hdu_i.header['RDNOIS11']=( np.float32(self.information['i_rdnois11']) , 'read noise (channel 11)')
hdu_i.header['RDNOIS12']=( np.float32(self.information['i_rdnois12']) , 'read noise (channel 12)')
hdu_i.header['RDNOIS13']=( np.float32(self.information['i_rdnois13']) , 'read noise (channel 13)')
hdu_i.header['RDNOIS14']=( np.float32(self.information['i_rdnois14']) , 'read noise (channel 14)')
hdu_i.header['RDNOIS15']=( np.float32(self.information['i_rdnois15']) , 'read noise (channel 15)')
hdu_i.header['RDNOIS16']=( np.float32(self.information['i_rdnois16']) , 'read noise (channel 16)')
#######
hdu_i.header['DETBIA1']=( np.float32(self.information['i_detbia1']) , 'amplifier bias voltage (channel 1)')
hdu_i.header['DETBIA2']=( np.float32(self.information['i_detbia2']) , 'amplifier bias voltage (channel 2)')
hdu_i.header['DETBIA3']=( np.float32(self.information['i_detbia3']) , 'amplifier bias voltage (channel 3)')
hdu_i.header['DETBIA4']=( np.float32(self.information['i_detbia4']) , 'amplifier bias voltage (channel 4)')
hdu_i.header['DETBIA5']=( np.float32(self.information['i_detbia5']) , 'amplifier bias voltage (channel 5)')
hdu_i.header['DETBIA6']=( np.float32(self.information['i_detbia6']) , 'amplifier bias voltage (channel 6)')
hdu_i.header['DETBIA7']=( np.float32(self.information['i_detbia7']) , 'amplifier bias voltage (channel 7)')
hdu_i.header['DETBIA8']=( np.float32(self.information['i_detbia8']) , 'amplifier bias voltage (channel 8)')
hdu_i.header['DETBIA9']=( np.float32(self.information['i_detbia9']) , 'amplifier bias voltage (channel 9)')
hdu_i.header['DETBIA10']=( np.float32(self.information['i_detbia10']) , 'amplifier bias voltage (channel 10)')
hdu_i.header['DETBIA11']=( np.float32(self.information['i_detbia11']) , 'amplifier bias voltage (channel 11)')
hdu_i.header['DETBIA12']=( np.float32(self.information['i_detbia12']) , 'amplifier bias voltage (channel 12)')
hdu_i.header['DETBIA13']=( np.float32(self.information['i_detbia13']) , 'amplifier bias voltage (channel 13)')
hdu_i.header['DETBIA14']=( np.float32(self.information['i_detbia14']) , 'amplifier bias voltage (channel 14)')
hdu_i.header['DETBIA15']=( np.float32(self.information['i_detbia15']) , 'amplifier bias voltage (channel 15)')
hdu_i.header['DETBIA16']=( np.float32(self.information['i_detbia16']) , 'amplifier bias voltage (channel 16)')
######
######
hdu_i.header['READT0'] =( hdu_g.header['READT0'] , 'read start time (UTC)')
hdu_i.header['READT1'] =( hdu_g.header['READT1'] , 'read end time (UTC)')
hdu_i.header['DETTEMP0']=( np.float32(0) , 'detector temperature at READT0')
hdu_i.header['DETTEMP1']=( np.float32(0) , 'detector temperature at READT1')
hdu_i.header['RDSPEED'] =( np.float32(3) , 'read speed (in MHz)')
##### Shutter information #####
hdu_i.header['SHUTHWV'] =( 'XXXXXXXXXXXXXXXX' , 'shutter hardware version')
hdu_i.header['SHUTSWV'] =( 'XXXXXXXXXXXXXXXX' , 'shutter software')
hdu_i.header['SHUTSTAT']=( bool(True) , 'shutter status,T: open, F: close')
hdu_i.header['SHTOPEN0']=( np.float64(0) , 'shutter open time (begin)')
hdu_i.header['SHTOPEN1']=( np.float64(0) , 'shutter open time (end)')
hdu_i.header['SHTCLOS0']=( np.float64(0) , 'shutter close time (begin)')
hdu_i.header['SHTCLOS1']=( np.float64(0) , 'shutter close time (end)')
################ finish MCI header for image layer
#write the actual file
hdulist_i=fits.HDUList([ofd_i, hdu_i])
#file_i=self.result_path+'/sky_Data/'+filename_i
if self.source=='DARK' or self.source=='FLAT' or self.source=='BIAS':
file_i=self.result_path+'/cali_Data/'+filename_i
else:
file_i=self.result_path+'/sky_Data/'+filename_i
hdulist_i.writeto(file_i, overwrite=True)
#################################################################################################
def simulate(self, simnumber):
"""
Create a single simulated image of a quadrant defined by the configuration file.
Will do all steps defined in the config file sequentially.
:return: None
"""
##############################################################################################
##############################################################################################
sourcelist=['XDF','CS','PI','TO']
self.source=self.information['sourcein']
self.information['ver']='1.0'
self._loadGhostModel()
self.load_filter_PSF()
self.information['simnumber']=simnumber
############################################################################
flag=0
for k in range(1,50,1):
fn=self.information['indata_path']+'/TianCe/orbit20160925/'+str(k)+'.txt';
d=np.loadtxt(fn);
self.dt_num=int(self.information['simnumber']*(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
self.orbit_pars=d
self.orbit_file_num=k
self.orbit_exp_num =kk
self.information['pos_x']=float(self.orbit_pars[self.orbit_exp_num,1])
self.information['pos_y']=float(self.orbit_pars[self.orbit_exp_num,2])
self.information['pos_z']=float(self.orbit_pars[self.orbit_exp_num,3])
self.information['velocity_x']=float(self.orbit_pars[self.orbit_exp_num,4])
self.information['velocity_y']=float(self.orbit_pars[self.orbit_exp_num,5])
self.information['velocity_z']=float(self.orbit_pars[self.orbit_exp_num,6])
exptime_start_jd=d[kk,0] #### jd time, utc format
self.dt=julian.from_jd(exptime_start_jd, fmt='jd')
#############################################################################
#############################################################################
if self.source in sourcelist:
self.earthshine_theta=30.0 # in degree
self.TianCe_day=str(self.dt.year)+self.dt.strftime("-%m-%d") ####str(self.dt.year)+'-'+str(self.dt.month)+'-'+str(self.dt.day)
self.TianCe_exp_start=dt2hmd(self.dt)
###################################
if self.sky_shift_rot:
np.random.seed()
ud = np.random.random() # Choose a random pointing ra in degree
dis_ra = 2*(ud-1)*10/3600.0
np.random.seed()
ud = np.random.random() # Choose a random pointing dec in degree
dis_dec = 2*(ud-1)*10/3600.0 #
np.random.seed()
ud = np.random.random() # Choose a random telescope rotation in degree
rotTelPos = 2*(ud-1)*5.0*galsim.degrees # the image will rotate in clockwise rotaion between -10-10 degree
rotSkyPos = 0.0*galsim.degrees
################################
else:
dis_ra = 0 #
dis_dec = 0 # T
rotTelPos = 0.0*galsim.degrees # the image will rotate in clockwise rotaion between -10-10 degree
rotSkyPos = 0.0*galsim.degrees
#####################################################
self.information['T_disRa'] =dis_ra # telescope pointing Ra displacement
self.information['T_disDec'] =dis_dec # telescope pointing Dec displacement
self.information['rotTelPos']=rotTelPos # telescope rotation in degree
self.information['rotSkyPos']=rotSkyPos # sky rotation in degree
self.information['name_obj'] ='MCI_obj'
theta = self.information['rotTelPos'] - self.information['rotSkyPos']
self.information['rotangle']=theta.deg
self.log.info('dis_ra(in pixel)=%f, dis_dec(in pixel)=%f, sky_rot(in deg)=%f' % (dis_ra*3600/0.05, dis_dec*3600/0.05, theta.deg))
#######################################################################################################
#############################################################################
########## simulate star images from star SED data ##################
if self.sim_star :
############################
starimg=self.cal_StarSED_img()
##########################################################################
# # ########################################################################
# ###### generate galaxy SED image without Lensing effect ##########
if self.sim_galaxy :
losimg=self.cal_Lensing_galaxy_img()
if self.sim_star ==False and self.sim_galaxy == False:
print('error: No image will be simulated...... ')
sys.exit(1)
#### stack the image
if 'losimg' in dir():
self.image_g+=losimg['g']
self.image_r+=losimg['r']
self.image_i+=losimg['i']
if 'starimg' in dir():
self.image_g+=starimg['g']
self.image_r+=starimg['r']
self.image_i+=starimg['i']
###########################################################################################################
#########################################################################################################
####### end if source =='SKY'
elif self.source=='FLAT':
theta=0
self.information['name_obj']='FLAT'
self.information['rotangle']=0
self.information['CRVAL1']=0
self.information['CRVAL2']=0
self.information['CRPIX1']=self.information['xsize']/2-0.5 ####
self.information['CRPIX2']=self.information['ysize']/2-0.5 ####
self.information['tel_Pra']=0
self.information['tel_Pdec']=0
self.information['CD1_1']=-np.cos(theta)*self.information['pixel_size'] ####
self.information['CD1_2']= np.sin(theta)*self.information['pixel_size'] ####
self.information['CD2_1']= np.sin(theta)*self.information['pixel_size'] ####
self.information['CD2_2']= np.cos(theta)*self.information['pixel_size'] ####
self.information['ra_obj'] =0
self.information['dec_obj'] =0
self.information['ra_pnt0'] =0
self.information['dec_pnt0']=0
self.information['target'] =0
self.image_g +=100*self.information['exptime']
self.image_r +=100*self.information['exptime']
self.image_i +=100*self.information['exptime']
elif self.source=='DARK':
theta=0
self.information['name_obj']='DARK'
self.information['rotangle']=0
self.information['CRVAL1']=0
self.information['CRVAL2']=0
self.information['tel_Pra']=0
self.information['tel_Pdec']=0
self.information['CRPIX1']=self.information['xsize']/2-0.5 ####
self.information['CRPIX2']=self.information['ysize']/2-0.5 ####
self.information['CD1_1']=-np.cos(theta)*self.information['pixel_size'] ####
self.information['CD1_2']= np.sin(theta)*self.information['pixel_size'] ####
self.information['CD2_1']= np.sin(theta)*self.information['pixel_size'] ####
self.information['CD2_2']= np.cos(theta)*self.information['pixel_size'] ####
self.information['ra_obj'] =0
self.information['dec_obj'] =0
self.information['ra_pnt0'] =0
self.information['dec_pnt0']=0
self.information['target'] =0
self.image_g=np.zeros((self.information['ysize'], self.information['xsize']), dtype=np.float64)
self.image_r=np.zeros((self.information['ysize'], self.information['xsize']), dtype=np.float64)
self.image_i=np.zeros((self.information['ysize'], self.information['xsize']), dtype=np.float64)
elif self.source=='BIAS':
theta=0
self.information['exptime']=0.0
self.information['name_obj']='BIAS'
self.information['rotangle']=0
self.information['CRVAL1']=0
self.information['CRVAL2']=0
self.information['tel_Pra']=0
self.information['tel_Pdec']=0
self.information['CRPIX1']=self.information['xsize']/2-0.5 ####
self.information['CRPIX2']=self.information['ysize']/2-0.5 ####
self.information['CD1_1']=-np.cos(theta)*self.information['pixel_size'] ####
self.information['CD1_2']= np.sin(theta)*self.information['pixel_size'] ####
self.information['CD2_1']= np.sin(theta)*self.information['pixel_size'] ####
self.information['CD2_2']= np.cos(theta)*self.information['pixel_size'] ####
self.information['ra_obj'] =0
self.information['dec_obj'] =0
self.information['ra_pnt0'] =0
self.information['dec_pnt0']=0
self.information['target'] =0
self.image_g=np.zeros((self.information['ysize'], self.information['xsize']), dtype=np.float64)
self.image_r=np.zeros((self.information['ysize'], self.information['xsize']), dtype=np.float64)
self.image_i=np.zeros((self.information['ysize'], self.information['xsize']), dtype=np.float64)
else:
print('Souce is not correct and programe will return')
sys.exit(1)
########################################################################################################################
########################################################################################################################
#Apply flat-field large scale structure for one chip
if self.flatfieldM:
self.applyFlat()
#print('applyFlat')
########################################################
if self.PRNUeffect:
self.applyPRNUeffect()
#print('applyPRNUeffect')
####################################################################
if self.source in sourcelist:
if self.cosmicRays:
self.addCosmicRays()
#print('addCosmicRays finisth')
##################################################
if self.skyback and self.source in sourcelist:
self.applyPoissonSkyNoise()
#print('apply Poisson Sky Noise')
#################################################
if self.darknoise:
self.applyDarkCurrent()
#print('applyDarkCurrent finisth')
#################################################
if self.shutterEffect:
self.addshuttereffect()
#print('add shutter effect ')
#################################################
if self.source not in sourcelist:
self.applyPoissonNoise()
###################################################
if self.nonlinearity:
self.applyNonlinearity()
#print('applyNonlinearity')
# ################################################
if self.bleeding:
self.image_g=self.applyBleeding(self.image_g.copy())
self.image_r=self.applyBleeding(self.image_r.copy())
self.image_i=self.applyBleeding(self.image_i.copy())
#print('apply bleeding effect finisth')
################################################
#### read CCD image to output matrix as defined
if self.overscans and self.information['xsize']==9216 and self.information['ysize']==9232 :
self.image_g=self.addScanEffect(self.image_g.copy())
self.image_r=self.addScanEffect(self.image_r.copy())
self.image_i=self.addScanEffect(self.image_i.copy())
#print('apply overscan zone 1-16')
#################################################################
if self.radiationDamage:
self.applyRadiationDamage()
#print('applyRadiationDamage()')
######################################################################
###### apply readoutNoise ######
if self.readoutNoise:
self.applyReadoutNoise()
#print('apply readout noise')
######################################################################
if self.information['xsize']==9216 and self.information['ysize']==9232 :
self.electrons2ADU() ### apply gain
#print('apply 1-16 zone by different gain ')
else:
self.image_g=self.image_g/1.5
self.image_r=self.image_r/1.5
self.image_i=self.image_i/1.5
#print('apply gain finish')
######################################################################
#size of the output image array, xsize is column, ysize is row, xsize = 9216,ysize = 9232
if self.information['xsize']==9216 and self.information['ysize']==9232 :
self.applyBias() ###
#print('apply bias finish')
else:
self.image_g=self.image_g+500.0
self.image_r=self.image_r+500.0
self.image_i=self.image_i+500.0
###################################################################
if self.cosmetics and self.information['xsize']==9216 and self.information['ysize']==9232 :
self.applyCosmetics()
#print('applyCosmetics')
##################################################
if self.intscale:
self.discretise()
##################################################
self.writeOutputs(simnumber)
self.log.info('Finished the simulation.')
#print('The iLoop= % d simlaiton finished. ' )
################################################################################################
def mcisim(iLoop,configfile):
opts, args = processArgs()
opts.configfile=configfile
simulate = MCIsimulator(opts)
simulate.configure(1)
simulate.simulate(1)
return
################################################################################################
################################################################################################
################################################################################################
################################### main function #############################################
#############################################################################
if __name__ == "__main__":
if len(sys.argv[:]) <2:
configfile='./mci_data/mci_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)
mcisim(1, configfile)
print('MCI simulation is successful!')
\ No newline at end of file
csst_mci_sim/mci_data/mci_sim_example.config 0 → 100644
View file @ aecb88aa
[TEST]
### user should define file path below
inData_path =/home/yan/MCI_sim_Fabu/MCI_inputData
result_path =/home/yan/MCI_sim_Fabu/mci_sim_result
#size of the output image array, xsize is column, ysize is row, xsize = 9216,ysize = 9232
xsize = 1000
ysize = 1000
##prescan and overscan , do not change!!!
prescan = 27
overscan = 320
###sourcein = XDF, DARK, FLAT, BIAS
sourcein=XDF
##################################################
#### ####
####Control flags (can be yes/no) ####
##################################################
cosmicRays = yes
darknoise = yes
cosmetics = yes
radiationDamage= no
bleeding = yes
overscans = yes
nonlinearity = yes
readoutNoise = yes
skyback = yes
TianceEffect = yes
intscale = yes
ghosts = yes
shutterEffect = yes
flatfieldM = yes
PRNUeffect = yes
appFatt = no
sky_shift_rot = yes
distortion = no
sim_star = yes
sim_galaxy = yes
##############################################
##############################################
#CCD properties
fullwellcapacity = 90000
#dark noise in electrons per second
dark = 0.001
#exposure to simulate, exposure time
exptime = 300.0
###PNRU matrix sigma
prnu_sigma=0.001
### cosmicray coveringFraction
coveringFraction=1.0
#offset from the object, note that at the moment this is fixed, but in reality a focal plane position dependent.
ghostOffsetX = 50
ghostOffsetY = 50
ghostRatio = 1e-04
####treering effect,bright fatter effect and difusion effect
treering=0.1
fatter=1.0
difusion=0.1
### the choosen Filters in three CCDs
filter_g=u
filter_r=F555W
filter_i=F814W
##['G_filters']=["F275W", "F280N","NUV", "WU", "CBU", "F343N", "u", "F373N", "F395N"," F336W"]
##['R_filters']=["F487N", "F502N", "CBV", "r", "F656N", "F658N", "F467M", "F555W", "F606W", "F673N"]
##['I_filters']=["z", "y", "F815N", "CBI", "F925N", "F960M", "F968N", "F845M" ,"F850LP" ,"F814W"]
####################
g_gain1=1.53
g_gain2=1.54
g_gain3=1.55
g_gain4=1.53
g_gain5=1.51
g_gain6=1.56
g_gain7=1.58
g_gain8=1.53
g_gain9=1.54
g_gain10=1.57
g_gain11=1.51
g_gain12=1.53
g_gain13=1.55
g_gain14=1.57
g_gain15=1.53
g_gain16=1.52
#######################
g_rdnois1=4.53
g_rdnois2=4.54
g_rdnois3=4.55
g_rdnois4=4.53
g_rdnois5=4.51
g_rdnois6=4.56
g_rdnois7=4.58
g_rdnois8=4.53
g_rdnois9=4.54
g_rdnois10=4.57
g_rdnois11=4.51
g_rdnois12=4.53
g_rdnois13=4.55
g_rdnois14=4.57
g_rdnois15=4.53
g_rdnois16=4.52
############################
g_detbia1=500
g_detbia2=510
g_detbia3=514
g_detbia4=520
g_detbia5=524
g_detbia6=540
g_detbia7=530
g_detbia8=532
g_detbia9=534
g_detbia10=526
g_detbia11=532
g_detbia12=516
g_detbia13=540
g_detbia14=560
g_detbia15=536
g_detbia16=528
############################
r_gain1=1.52
r_gain2=1.54
r_gain3=1.52
r_gain4=1.55
r_gain5=1.56
r_gain6=1.52
r_gain7=1.54
r_gain8=1.6
r_gain9=1.54
r_gain10=1.52
r_gain11=1.54
r_gain12=1.57
r_gain13=1.52
r_gain14=1.55
r_gain15=1.52
r_gain16=1.55
################################
r_rdnois1=4.23
r_rdnois2=4.24
r_rdnois3=4.25
r_rdnois4=4.23
r_rdnois5=4.21
r_rdnois6=4.26
r_rdnois7=4.28
r_rdnois8=4.23
r_rdnois9=4.24
r_rdnois10=4.27
r_rdnois11=4.21
r_rdnois12=4.23
r_rdnois13=4.25
r_rdnois14=4.27
r_rdnois15=4.23
r_rdnois16=4.22
###################################
r_detbia1=600
r_detbia2=610
r_detbia3=614
r_detbia4=620
r_detbia5=624
r_detbia6=640
r_detbia7=630
r_detbia8=632
r_detbia9=634
r_detbia10=626
r_detbia11=632
r_detbia12=616
r_detbia13=640
r_detbia14=660
r_detbia15=636
r_detbia16=628
###################################
i_gain1=1.62
i_gain2=1.64
i_gain3=1.62
i_gain4=1.65
i_gain5=1.66
i_gain6=1.62
i_gain7=1.64
i_gain8=1.6
i_gain9=1.64
i_gain10=1.62
i_gain11=1.64
i_gain12=1.67
i_gain13=1.62
i_gain14=1.65
i_gain15=1.62
i_gain16=1.65
###################################
i_rdnois1=4.63
i_rdnois2=4.64
i_rdnois3=4.65
i_rdnois4=4.63
i_rdnois5=4.61
i_rdnois6=4.66
i_rdnois7=4.68
i_rdnois8=4.63
i_rdnois9=4.64
i_rdnois10=4.67
i_rdnois11=4.61
i_rdnois12=4.63
i_rdnois13=4.65
i_rdnois14=4.67
i_rdnois15=4.63
i_rdnois16=4.62
#######################################
i_detbia1=400
i_detbia2=410
i_detbia3=414
i_detbia4=420
i_detbia5=424
i_detbia6=440
i_detbia7=430
i_detbia8=432
i_detbia9=434
i_detbia10=426
i_detbia11=432
i_detbia12=416
i_detbia13=440
i_detbia14=460
i_detbia15=436
i_detbia16=428
####### end #######################
csst_mci_sim/support/MCIinstrumentModel.py 0 → 100644
View file @ aecb88aa
"""
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 matplotlib
# import matplotlib.pyplot as plt
# import datetime, math
# import numpy as np
# import numexpr as ne
def MCIinformation():
"""
Returns a dictionary describing MCI CCD. The following information is provided (id: value - reference)::
dob: 0 - CDM03 (Short et al. 2010)
fwc: 90000 - CCD spec EUCL-EST-RS-6-002 (for CDM03)
rdose: 30000000000.0 - derived (above the PLM requirement)
sfwc: 730000.0 - CDM03 (Short et al. 2010), see also the CCD spec EUCL-EST-RS-6-002
st: 5e-06 - CDM03 (Short et al. 2010)
svg: 1e-10 - CDM03 (Short et al. 2010)
t: 0.01024 - CDM03 (Short et al. 2010)
trapfile: cdm_euclid.dat - CDM03 (derived, refitted to CCD204 data)
vg: 6e-11 - CDM03 (Short et al. 2010)
vth: 11680000.0 - CDM03 (Short et al. 2010)
:return: instrument model parameters
:rtype: dict
"""
#########################################################################################################
out=dict()
out.update({'dob' : 0, 'rdose' : 8.0e9,
'parallelTrapfile' : 'cdm_euclid_parallel.dat', 'serialTrapfile' : 'cdm_euclid_serial.dat',
'beta_s' : 0.6, 'beta_p': 0.6, 'fwc' : 90000, 'vth' : 1.168e7, 't' : 20.48e-3, 'vg' : 6.e-11,
'st' : 5.0e-6, 'sfwc' : 730000., 'svg' : 1.0e-10})
return out
def CCDnonLinearityModel(data, beta=6e-7):
"""
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()
csst_mci_sim/support/cosmicrays.py 0 → 100644
View file @ aecb88aa
"""
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__":
import sys
sys.path.append('/home/yan/csst-master/')
from support import logger as lg
from support import files as fileIO
from scipy import ndimage
from astropy.io import fits
#set up logger
log = lg.setUpLogger('VISsim.log')
#test section
crImage = np.zeros((2066, 2048), dtype=np.float64)
#cosmic ray instance
cosmics = cosmicrays(log, crImage)
#add cosmic rays up to the covering fraction
CCD_cr = cosmics.addUpToFraction(1.4, limit=None, verbose=True)
print(CCD_cr)
print(type(CCD_cr))
effected = np.count_nonzero(CCD_cr)
print (effected, effected*100./(CCD_cr.shape[0]*CCD_cr.shape[1]))
#save to FITS
fits.writeto(CCD_cr, 'cosmicrayTest.fits', overwrite=True)
#smooth with a charge diffusion kernel
smooth = ndimage.filters.gaussian_filter(CCD_cr, (0.32, 0.32))
fits.writeto(smooth, 'cosmicrayTestSmoothed.fits')
csst_mci_sim/support/logger.py 0 → 100644
View file @ aecb88aa
"""
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
class SimpleLogger(object):
"""
A simple class to create a log file or print the information on screen.
"""
def __init__(self, filename, verbose=False):
self.file = open(filename, 'w')
self.verbose = verbose
def write(self, text):
"""
Writes text either to file or screen.
"""
print >> self.file, text
if self.verbose: print( text)
csst_mci_sim/support/sed.py 0 → 100644
View file @ aecb88aa
#!/usr/bin/env python
# -*- coding:utf-8 -*-
# @Author: Shuai Feng (hebtu.edu.cn)
# @Time: 2022-09-25
import numpy as np
from scipy.interpolate import interp1d
from astropy.io import fits
from astropy.io import ascii
from astropy import units as u
# ----------------
# Magnitude Module
from scipy.interpolate import interp1d
def Calzetti_Law(wave, Rv = 4.05):
"""Dust Extinction Curve by Calzetti et al. (2000)
Args:
wave (float): Wavelength
Rv (float, optional): Extinction curve. Defaults to 4.05.
Returns:
float: Extinction value E(B-V)
"""
wave_number = 1./(wave * 1e-4)
reddening_curve = np.zeros(len(wave))
idx = np.logical_and(wave >= 1200, wave <= 6300)
reddening_curve[idx] = 2.659 * ( -2.156 + 1.509 * wave_number[idx] - 0.198 * \
(wave_number[idx] ** 2)) + 0.011 * (wave_number[idx] **3 ) + Rv
idx = np.logical_and(wave >= 6300, wave <= 22000)
reddening_curve[idx] = 2.659 * ( -1.857 + 1.040 * wave_number[idx]) + Rv
return reddening_curve
def reddening(wave, flux, ebv = 0.0, law = 'calzetti', Rv = 4.05):
"""
Reddening an input spectra through a given reddening curve.
Args:
wave (float): Wavelength of input spectra
flux (float): Flux of input spectra
ebv (float, optional): Extinction value. Defaults to 0.
law (str, optional): Extinction curve. Defaults to 'calzetti'.
Rv (float, optional): _description_. Defaults to 4.05.
Returns:
float: Flux of spectra after reddening.
"""
if law == 'calzetti':
curve = Calzetti_Law(wave, Rv = Rv)
fluxNew = flux / (10. ** (0.4 * ebv * curve))
return fluxNew
def flux_to_mag(wave, flux, filepath, band='GAIA_bp'):
"""Convert flux of given spectra to magnitude
Args:
wave (float): Wavelength
flux (float): Flux of spectra
band (str, optional): Filter band name. Defaults to 'GAIA_bp'.
Returns:
float: value of magnitude
"""
# /home/yan/MCI_sim/MCI_input/SED_Code/data
import os
#parent = os.path.dirname(os.path.realpath(__file__))
parent=filepath
band = ascii.read(parent+'/SED_Code/seddata/' + band + '.dat')
wave0= band['col1']
curv0= band['col2']
# Setting the response
func = interp1d(wave0, curv0)
response = np.copy(wave)
ind_extra = (wave > max(wave0)) | (wave < min(wave0))
response[ind_extra] = 0
ind_inside = (wave < max(wave0)) & (wave > min(wave0))
response[ind_inside] = func(wave[ind_inside])
# Total Flux
Tflux = np.trapz(flux * response, wave) / np.trapz(response, wave)
return -2.5 * np.log10(Tflux)
def calibrate(wave, flux, mag, filepath,band='GAIA_bp'):
"""
Calibrate the spectra according to the magnitude.
Args:
wave (float): Wavelength
flux (float): Flux of spectra
mag (float): Input magnitude.
band (str, optional): Filter band name. Defaults to 'GAIA_bp'.
Returns:
float: Flux of calibrated spectra. Units: 1e-17 erg/s/A/cm^2
"""
inst_mag = flux_to_mag(wave, flux,filepath ,band = band)
instflux = 10 ** (-0.4 * inst_mag)
realflux = (mag * u.STmag).to(u.erg/u.s/u.cm**2/u.AA).value
# Normalization
flux_ratio = realflux / instflux
flux_calibrate = flux * flux_ratio * 1e17 # Units: 10^-17 erg/s/A/cm^2
return flux_calibrate
# ------------
# SED Template
class Gal_Temp():
"""
Template of Galaxy SED
"""
def __init__(self,parent):
import os
self.parent = parent
hdulist = fits.open(self.parent+'/SED_Code/seddata/galaxy_temp.fits')
self.wave = hdulist[1].data['wave']
self.flux = hdulist[2].data
self.age_grid = hdulist[3].data['logAge']
self.feh_grid = hdulist[3].data['FeH']
def toMag(self, redshift = 0):
"""Calculating magnitude
Args:
redshift (float, optional): redshift of spectra. Defaults to 0.
"""
wave = self.wave * (1 + redshift)
self.umag = flux_to_mag(wave, self.flux, self.parent,band='SDSS_u')
self.gmag = flux_to_mag(wave, self.flux, self.parent,band='SDSS_g')
self.rmag = flux_to_mag(wave, self.flux, self.parent,band='SDSS_r')
self.imag = flux_to_mag(wave, self.flux, self.parent,band='SDSS_i')
self.zmag = flux_to_mag(wave, self.flux, self.parent,band='SDSS_z')
class Star_Temp():
"""
Template of Stellar SED
"""
def __init__(self, parent):
#import os
#parent = os.path.dirname(os.path.realpath(__file__))
##/home/yan/MCI_sim_Fabu/MCI_inputData/SED_Code
self.parent = parent
hdulist = fits.open(parent+'/SED_Code/seddata/stellar_temp.fits')
self.wave = hdulist[1].data['wave']
self.flux = hdulist[2].data
self.Teff_grid = hdulist[3].data['Teff']
self.FeH_grid = hdulist[3].data['FeH']
def toMag(self):
wave = self.wave
self.bpmag = flux_to_mag(wave, self.flux, self.parent,band='GAIA_bp')
self.rpmag = flux_to_mag(wave, self.flux, self.parent,band='GAIA_rp')
# -------------
# SED Modelling
def Model_Stellar_SED(wave, bp, rp, temp, filepath):
"""Modelling stellar SED based on bp, rp magnitude
Args:
wave (float): Wavelength
bp (float): Magnitude of GAIA BP band
rp (float): Magnitude of GAIA RP band
temp (class): Class of stellar template
Returns:
float array: Spectral energy distribution of stellar SED,
which have the same length to the input wave
"""
color0 = bp - rp
colors = temp.bpmag - temp.rpmag
idx = np.argmin(np.abs(colors - color0))
flux0 = temp.flux[idx]
flux1 = np.interp(wave, temp.wave, flux0)
flux = calibrate(wave, flux1, rp, filepath, band = 'GAIA_rp')
return flux
def Model_Galaxy_SED(wave, ugriz, z, temp):
"""Modelling galaxy SED based on u,g,r,i,z magnitude
Args:
wave (float): Wavelength
ugriz (float, array): The array of magnitude of SDSS ugriz band
z (float): Redshift
temp (class): Class of gaalxy template
Returns:
float array: Spectral energy distribution of stellar SED,
which have the same length to the input wave
"""
sed = 10. ** (-0.4 * ugriz)
sed = sed / sed[2]
ntemp = len(temp.rmag)
dmag = np.zeros(ntemp)
for j in range(ntemp):
ugriz0 = np.array([temp.umag[j], temp.gmag[j], temp.rmag[j],
temp.imag[j], temp.zmag[j]])
sed0 = 10. ** (-0.4 * ugriz0)
sed0 = sed0 / sed0[2]
dmag[j] = np.sum(np.abs(sed - sed0))
idx = np.argmin(dmag)
flux0 = temp.flux[idx]
# Effect of E(B-V)
ri0 = ugriz[2] - ugriz[3]
ri = temp.rmag - temp.imag
dri = ri0 - ri[idx]
Alambda = Calzetti_Law(np.array([6213 / (1 + z), 7625 / (1 + z)]))
eri0 = (Alambda[0] - Alambda[1])
ebv = dri/eri0
if ebv<0:
ebv=0
if ebv>0.5:
ebv=0.5
flux1 = reddening(temp.wave, flux0, ebv = ebv)
flux2 = np.interp(wave, temp.wave * (1 + z), flux1)
flux = calibrate(wave, flux2, ugriz[2], band = 'SDSS_r')
return flux
\ No newline at end of file
csst_mci_sim/support/shao.py 0 → 100644
View file @ aecb88aa
from ctypes import *
def checkInputList(input_list, n):
if type(input_list) != type([1, 2, 3]):
raise TypeError("Input type is not list!", input_list)
for i in input_list:
if type(i) != type(1.1):
if type(i) != type(1):
raise TypeError("Input list's element is not float or int!", input_list)
if len(input_list) != n:
raise RuntimeError("Length of input list is not equal to stars' number!", input_list)
def onOrbitObsPosition(input_ra_list, input_dec_list, input_pmra_list, input_pmdec_list, input_rv_list, \
input_parallax_list, input_nstars, input_x, input_y, input_z, input_vx, input_vy, \
input_vz, input_epoch, input_date_str, input_time_str):
#Check input parameters
if type(input_nstars) != type(1):
raise TypeError("Parameter 7 is not int!", input_nstars)
checkInputList(input_ra_list, input_nstars)
checkInputList(input_dec_list, input_nstars)
checkInputList(input_pmra_list, input_nstars)
checkInputList(input_pmdec_list, input_nstars)
checkInputList(input_rv_list, input_nstars)
checkInputList(input_parallax_list, input_nstars)
if type(input_x) != type(1.1):
raise TypeError("Parameter 8 is not double!", input_x)
if type(input_y) != type(1.1):
raise TypeError("Parameter 9 is not double!", input_y)
if type(input_z) != type(1.1):
raise TypeError("Parameter 10 is not double!", input_z)
if type(input_vx) != type(1.1):
raise TypeError("Parameter 11 is not double!", input_vx)
if type(input_vy) != type(1.1):
raise TypeError("Parameter 12 is not double!", input_vy)
if type(input_vz) != type(1.1):
raise TypeError("Parameter 13 is not double!", input_vz)
#Convert km -> m
input_x = input_x*1000.0
input_y = input_y*1000.0
input_z = input_z*1000.0
input_vx = input_vx*1000.0
input_vy = input_vy*1000.0
input_vz = input_vz*1000.0
if type(input_date_str) != type("2025-03-05"):
raise TypeError("Parameter 15 is not string!", input_date_str)
else:
input_date_str = input_date_str.strip()
if not (input_date_str[4]=="-" and input_date_str[7]=="-"):
raise TypeError("Parameter 15 format error (1)!", input_date_str)
else:
tmp = input_date_str.split("-")
if len(tmp) != 3:
raise TypeError("Parameter 15 format error (2)!", input_date_str)
input_year = int(tmp[0])
input_month = int(tmp[1])
input_day = int(tmp[2])
if not (input_year>=1900 and input_year<=2100):
raise TypeError("Parameter 15 year range error [1900 ~ 2100]!", input_year)
if not (input_month>=1 and input_month<=12):
raise TypeError("Parameter 15 month range error [1 ~ 12]!", input_month)
if not (input_day>=1 and input_day<=31):
raise TypeError("Parameter 15 day range error [1 ~ 31]!", input_day)
if type(input_time_str) != type("20:15:15.15"):
raise TypeError("Parameter 16 is not string!", input_time_str)
else:
input_time_str = input_time_str.strip()
if not (input_time_str[2]==":" and input_time_str[5]==":"):
raise TypeError("Parameter 16 format error (1)!", input_time_str)
else:
tmp = input_time_str.split(":")
if len(tmp) != 3:
raise TypeError("Parameter 16 format error (2)!", input_time_str)
input_hour = int(tmp[0])
input_minute = int(tmp[1])
input_second = float(tmp[2])
if not (input_hour>=0 and input_hour<=23):
raise TypeError("Parameter 16 hour range error [0 ~ 23]!", input_hour)
if not (input_minute>=0 and input_minute<=59):
raise TypeError("Parameter 16 minute range error [0 ~ 59]!", input_minute)
if not (input_second>=0 and input_second<60.0):
raise TypeError("Parameter 16 second range error [0 ~ 60)!", input_second)
#Inital dynamic lib
# import os
# currfile=os.getcwd()
# print(currfile)
shao = cdll.LoadLibrary('./mci_so/libshao.so')
shao.onOrbitObs.restype = c_int
d3 = c_double * 3
shao.onOrbitObs.argtypes = [c_double, c_double, c_double, c_double, c_double, c_double, \
c_int, c_int, c_int, c_int, c_int, c_double, \
c_double, POINTER(d3), POINTER(d3), \
c_int, c_int, c_int, c_int, c_int, c_double, \
POINTER(c_double), POINTER(c_double) ]
output_ra_list = list()
output_dec_list = list()
for i in range(input_nstars):
input_ra = c_double(input_ra_list[i])
input_dec = c_double(input_dec_list[i])
input_pmra = c_double(input_pmra_list[i])
input_pmdec = c_double(input_pmdec_list[i])
input_rv = c_double(input_rv_list[i])
input_parallax = c_double(input_parallax_list[i])
p3 = d3(input_x, input_y, input_z)
v3 = d3(input_vx, input_vy, input_vz)
input_year_c = c_int(input_year)
input_month_c = c_int(input_month)
input_day_c = c_int(input_day)
input_hour_c = c_int(input_hour)
input_minute_c = c_int(input_minute)
input_second_c = c_double(input_second)
DAT = c_double(37.0)
output_ra = c_double(0.0)
output_dec = c_double(0.0)
rs = shao.onOrbitObs(input_ra, input_dec, input_parallax, input_pmra, input_pmdec, input_rv, \
input_year_c, input_month_c, input_day_c, input_hour_c, input_minute_c, input_second_c, \
DAT, byref(p3), byref(v3), \
input_year_c, input_month_c, input_day_c, input_hour_c, input_minute_c, input_second_c, \
byref(output_ra), byref(output_dec))
if rs != 0:
raise RuntimeError("Calculate error!")
output_ra_list.append(output_ra.value)
output_dec_list.append(output_dec.value)
return output_ra_list, output_dec_list
\ No newline at end of file
setup.py 0 → 100644
View file @ aecb88aa
import setuptools
with open("README.md", "r") as fh:
long_description = fh.read()
setuptools.setup(
name='csst_mci_sim',
version='2.0.0-MCI1.0.0',
author='CSST Team',
author_email='zhaojunyan@shao.ac.cn',
description='The CSST - mci - 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_mci_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_mci_sim': 'csst_mci_sim'},
# include_package_data=True,
package_data={"": ["LICENSE", "README.md"],
"csst_mci_sim": ["mci_so/*", "mci_data/*"]},
# 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