diff --git a/LICENSE b/LICENSE new file mode 100644 index 0000000000000000000000000000000000000000..a9b5e744c0d05f2592bc82552323c2044a5bfbce --- /dev/null +++ b/LICENSE @@ -0,0 +1,21 @@ +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. diff --git a/README.md b/README.md index d2bc7fa0298dea340be95866e6430c0ea6fe9a25..aed373c9603307c3bddd511cc45de92fabdf830a 100644 --- a/README.md +++ b/README.md @@ -1,4 +1,4 @@ -# 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 diff --git a/csst_mci_sim/CTI/CTI.py b/csst_mci_sim/CTI/CTI.py new file mode 100644 index 0000000000000000000000000000000000000000..48c87cb8aee4207a0b4d764e32a0cea44a3d6987 --- /dev/null +++ b/csst_mci_sim/CTI/CTI.py @@ -0,0 +1,495 @@ +""" +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) + + +################################################################################################################# + diff --git a/csst_mci_sim/CTI/__init__.py b/csst_mci_sim/CTI/__init__.py new file mode 100644 index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391 diff --git a/csst_mci_sim/__init__.py b/csst_mci_sim/__init__.py new file mode 100644 index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391 diff --git a/csst_mci_sim/csst_mci_sim.py b/csst_mci_sim/csst_mci_sim.py new file mode 100644 index 0000000000000000000000000000000000000000..78f69aead92b3596119ca099f7b2ec5491ba8a2b --- /dev/null +++ b/csst_mci_sim/csst_mci_sim.py @@ -0,0 +1,3733 @@ +#!/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 j0): + 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]) & (idxt_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 xminnp.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() 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]=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 diff --git a/csst_mci_sim/mci_data/mci_sim_example.config b/csst_mci_sim/mci_data/mci_sim_example.config new file mode 100644 index 0000000000000000000000000000000000000000..4b9e711bab987afd4b6443b3802bfd4e67996f92 --- /dev/null +++ b/csst_mci_sim/mci_data/mci_sim_example.config @@ -0,0 +1,282 @@ +[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 ####################### + + + + + + + + + + + + diff --git a/csst_mci_sim/mci_so/__init__.py b/csst_mci_sim/mci_so/__init__.py new file mode 100644 index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391 diff --git a/csst_mci_sim/mci_so/cdm03.cpython-37m-x86_64-linux-gnu.so b/csst_mci_sim/mci_so/cdm03.cpython-37m-x86_64-linux-gnu.so new file mode 100644 index 0000000000000000000000000000000000000000..85a4cfc31920218e8b634b8d53cf9db62547f63f Binary files /dev/null and b/csst_mci_sim/mci_so/cdm03.cpython-37m-x86_64-linux-gnu.so differ diff --git a/csst_mci_sim/mci_so/cdm03.cpython-38-x86_64-linux-gnu.so b/csst_mci_sim/mci_so/cdm03.cpython-38-x86_64-linux-gnu.so new file mode 100644 index 0000000000000000000000000000000000000000..7bfa02f3f33ab5d5a16e4d35ac622d8bf7a49ed1 Binary files /dev/null and b/csst_mci_sim/mci_so/cdm03.cpython-38-x86_64-linux-gnu.so differ diff --git a/csst_mci_sim/mci_so/cdm03bidir.cpython-37m-x86_64-linux-gnu.so b/csst_mci_sim/mci_so/cdm03bidir.cpython-37m-x86_64-linux-gnu.so new file mode 100644 index 0000000000000000000000000000000000000000..a63494c33d3e3f2c0dd42085e549d38415c682c6 Binary files /dev/null and b/csst_mci_sim/mci_so/cdm03bidir.cpython-37m-x86_64-linux-gnu.so differ diff --git a/csst_mci_sim/mci_so/cdm03bidir.cpython-38-x86_64-linux-gnu.so b/csst_mci_sim/mci_so/cdm03bidir.cpython-38-x86_64-linux-gnu.so new file mode 100644 index 0000000000000000000000000000000000000000..4ddf6dd453c91a1ed8b3448956bc127596d748f8 Binary files /dev/null and b/csst_mci_sim/mci_so/cdm03bidir.cpython-38-x86_64-linux-gnu.so differ diff --git a/csst_mci_sim/mci_so/libshao.so b/csst_mci_sim/mci_so/libshao.so new file mode 100644 index 0000000000000000000000000000000000000000..47d97a44a74c8f011914dc07f52d8957c81eff39 Binary files /dev/null and b/csst_mci_sim/mci_so/libshao.so differ diff --git a/csst_mci_sim/support/MCIinstrumentModel.py b/csst_mci_sim/support/MCIinstrumentModel.py new file mode 100644 index 0000000000000000000000000000000000000000..8777e9e66b1ab47151cc95baf7ef2e57ce846905 --- /dev/null +++ b/csst_mci_sim/support/MCIinstrumentModel.py @@ -0,0 +1,68 @@ +""" +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() diff --git a/csst_mci_sim/support/__init__.py b/csst_mci_sim/support/__init__.py new file mode 100644 index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391 diff --git a/csst_mci_sim/support/cosmicrays.py b/csst_mci_sim/support/cosmicrays.py new file mode 100644 index 0000000000000000000000000000000000000000..3b8e4bafe679feaf540d6dd8afd24ce67797af9f --- /dev/null +++ b/csst_mci_sim/support/cosmicrays.py @@ -0,0 +1,448 @@ +""" +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') diff --git a/csst_mci_sim/support/logger.py b/csst_mci_sim/support/logger.py new file mode 100644 index 0000000000000000000000000000000000000000..68253e7338eeb33aaf334dafa862d96c71fa6a33 --- /dev/null +++ b/csst_mci_sim/support/logger.py @@ -0,0 +1,55 @@ +""" +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) diff --git a/csst_mci_sim/support/sed.py b/csst_mci_sim/support/sed.py new file mode 100644 index 0000000000000000000000000000000000000000..dccb65403796c1f69a76507ade362c717cdf761c --- /dev/null +++ b/csst_mci_sim/support/sed.py @@ -0,0 +1,245 @@ +#!/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 diff --git a/csst_mci_sim/support/shao.py b/csst_mci_sim/support/shao.py new file mode 100644 index 0000000000000000000000000000000000000000..3717c9fb93fd242bd9b1f524bd5cfe919dd4652d --- /dev/null +++ b/csst_mci_sim/support/shao.py @@ -0,0 +1,134 @@ +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 diff --git a/setup.py b/setup.py new file mode 100644 index 0000000000000000000000000000000000000000..84c12698eb53e8e5f82ea89787061a61f5de6687 --- /dev/null +++ b/setup.py @@ -0,0 +1,36 @@ +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', +)