Skip to content
GitLab
Commit d6274e12 authored by Yan Zhaojun's avatar Yan Zhaojun
Browse files

debug

parent 91e3a59d
Pipeline #4461 failed with stage
in 0 seconds
Show whitespace changes
Inline Side-by-side
csst_ifs_sim/csst_ifs_sim.py
View file @ d6274e12
...@@ -5,6 +5,27 @@ Created on Thu Apr 11 15:18:57 2024 ...@@ -5,6 +5,27 @@ Created on Thu Apr 11 15:18:57 2024
@author: yan @author: yan
""" """
from astropy.utils.iers import conf
import galsim
from scipy.integrate import simps
import pandas as pd
from datetime import datetime, timedelta
from astropy.coordinates import get_sun
from astropy.time import Time
from astropy import units as u
from astropy.coordinates import SkyCoord, EarthLocation
from scipy import interpolate
import numpy as np
import scipy.io as sio
from astropy.io import fits
from scipy.signal import fftconvolve
from scipy import ndimage
import cmath
import configparser as ConfigParser
from support import IFSinstrumentModel
from support import cosmicrays
from support import logger as lg
from CTI import CTI
import os import os
from scipy.interpolate import interp1d from scipy.interpolate import interp1d
import astropy.coordinates as coord import astropy.coordinates as coord
...@@ -15,27 +36,6 @@ import sys ...@@ -15,27 +36,6 @@ import sys
# from csst_ifs_sim.support import logger as lg # from csst_ifs_sim.support import logger as lg
# from csst_ifs_sim.CTI import CTI # from csst_ifs_sim.CTI import CTI
sys.path.append('./csst_ifs_sim') sys.path.append('./csst_ifs_sim')
from CTI import CTI
from support import logger as lg
from support import cosmicrays
from support import IFSinstrumentModel
import configparser as ConfigParser
import cmath
from scipy import ndimage
from scipy.signal import fftconvolve
from astropy.io import fits
import scipy.io as sio
import numpy as np
from scipy import interpolate
from astropy.coordinates import SkyCoord, EarthLocation
from astropy import units as u
from astropy.time import Time
from astropy.coordinates import get_sun
from datetime import datetime, timedelta
import pandas as pd
from scipy.integrate import simps
import galsim
from astropy.utils.iers import conf
conf.auto_max_age = None conf.auto_max_age = None
""" """
...@@ -56,10 +56,7 @@ The approximate sequence of events in the simulator is as follows: ...@@ -56,10 +56,7 @@ The approximate sequence of events in the simulator is as follows:
#. Read in another file containing charge trap definitions (for CTI modelling). #. Read in another file containing charge trap definitions (for CTI modelling).
#. Read in a file defining the cosmic rays (trail lengths and #. Read in a file defining the cosmic rays (trail lengths and
cumulative distributions). cumulative distributions).
#. Load the wavefront aberration data used to calculate PSF with defined #. Load the wavefront aberration data used to calculate PSF .
wavelength and field of view.
#. Apply calibration unit flux to mimic flat field exposures [optional].
#. Apply a multiplicative flat-field map to emulate pixel-to-pixel #. Apply a multiplicative flat-field map to emulate pixel-to-pixel
non-uniformity. non-uniformity.
#. Add cosmic ray tracks onto the CCD with random positions but known #. Add cosmic ray tracks onto the CCD with random positions but known
...@@ -67,12 +64,11 @@ distribution. ...@@ -67,12 +64,11 @@ distribution.
#. Apply detector charge bleeding in column direction. #. Apply detector charge bleeding in column direction.
#. Add constant dark current and background light from Zodiacal #. Add constant dark current and background light from Zodiacal
light. light.
#. Include spatially uniform scattered light to the pixel grid.
#. Add photon (Poisson) noise #. Add photon (Poisson) noise
#. Add cosmetic defects from an input file. #. Add cosmetic defects from an input file.
#. Add pre- and overscan regions in the serial direction. #. Add pre- and overscan regions in the serial direction.
#. Apply the CDM03 radiation damage model [optional]. #. Apply the CDM03 radiation damage model.
#. Apply CCD273 non-linearity model to the pixel data. #. Apply non-linearity model to the pixel data.
#. Add readout noise selected from a Gaussian distribution. #. Add readout noise selected from a Gaussian distribution.
#. Convert from electrons to ADUs using a given gain factor. #. Convert from electrons to ADUs using a given gain factor.
#. Add a given bias level and discretise the counts #. Add a given bias level and discretise the counts
...@@ -82,17 +78,17 @@ Warning:: The code is still work in progress and new features are being added. ...@@ -82,17 +78,17 @@ 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 The code has been tested, but nevertheless bugs may be lurking in corners, so
please report any weird or inconsistent simulations to the author. please report any weird or inconsistent simulations to the author.
Note:: This class is Python 3 compatible. Note:: This class is Python 3 compatible.
2023.06.13 , add new lamp hole simulation 2023.06.13 , add new lamp hole simulation
2023.12.26 add date frame transfer effect and software version 2023.12.26 add date frame transfer effect and software version
2024.1.19 add multiple exposure mode function; 2024.1.19 add multiple exposure mode function;
2024.3.20 --- update the fits header as defined order 0 data format 2024.3.20 --- update the fits header as defined order 0 data format
2024.5.10 ---updata and correct the bug of frame transfer effect simulation
""" """
#### functions definition ##### #### functions definition #####
def transRaDec2D(ra, dec): def transRaDec2D(ra, dec):
""" """
...@@ -116,52 +112,6 @@ def transRaDec2D(ra, dec): ...@@ -116,52 +112,6 @@ def transRaDec2D(ra, dec):
return np.array([x1, y1, z1]) return np.array([x1, y1, z1])
############################################################################### ###############################################################################
# def flux2ill(wave, flux):
# """
# Parameters
# ----------
# wave : TYPE
# DESCRIPTION.
# flux : TYPE
# DESCRIPTION.
# Returns
# -------
# E : TYPE
# DESCRIPTION.
# """
# # erg/s/cm^2/A/arcsec^2 to W/m^2
# # 1 W/m^2/sr/μm = 0.10 erg/cm^2/s/sr/A
# # 1 sr = 1 rad^2 = 4.25452e10 arcsec^2
# # 1 J/s = 1 W
# # 1 J = 10^7 erg
# # convert erg/s/cm^2/A/arcsec^2 to erg/s/cm^2/A/sr
# flux1 = flux / (1/4.25452e10)
# # convert erg/s/cm^2/A/sr to W/m^2/sr/um
# flux2 = flux1 * 10
# # 对接收面积积分，输出单位 W/m^2/nm
# D = 2 # meter
# f = 28 # meter
# flux3 = flux2 * np.pi * D**2 / 4 / f**2 / 10**3
# # 对波长积分
# f = interp1d(wave, flux3)
# wave_interp = np.arange(3800, 7800)
# flux3_interp = f(wave_interp)
# # 输出单位 W/m^2
# delta_lamba = 0.1 # nm
# E = np.sum(flux3_interp * delta_lamba)
# return E
################################################################ ################################################################
...@@ -214,7 +164,8 @@ def ill2flux(path, E): ...@@ -214,7 +164,8 @@ def ill2flux(path, E):
############################################################## ##############################################################
########################################################################################## ##############################################################################
def earth_angle(time_jd, x_sat, y_sat, z_sat, ra_obj, dec_obj): def earth_angle(time_jd, x_sat, y_sat, z_sat, ra_obj, dec_obj):
""" """
...@@ -268,6 +219,7 @@ def earth_angle(time_jd, x_sat, y_sat, z_sat, ra_obj, dec_obj): ...@@ -268,6 +219,7 @@ def earth_angle(time_jd, x_sat, y_sat, z_sat, ra_obj, dec_obj):
##################################################################### #####################################################################
class StrayLight(object): class StrayLight(object):
def __init__(self, path, jtime=2460843., sat=np.array([0, 0, 0]), radec=np.array([0, 0])): def __init__(self, path, jtime=2460843., sat=np.array([0, 0, 0]), radec=np.array([0, 0])):
...@@ -311,6 +263,7 @@ class StrayLight(object): ...@@ -311,6 +263,7 @@ class StrayLight(object):
self.slcdll.Init(str.encode(self.deFn), str.encode( self.slcdll.Init(str.encode(self.deFn), str.encode(
self.PSTFn), str.encode(self.RFn), str.encode(self.ZolFn)) self.PSTFn), str.encode(self.RFn), str.encode(self.ZolFn))
############################################################################ ############################################################################
def caculateStarLightFilter(self, filter='i'): def caculateStarLightFilter(self, filter='i'):
""" """
...@@ -326,7 +279,8 @@ class StrayLight(object): ...@@ -326,7 +279,8 @@ class StrayLight(object):
DESCRIPTION. DESCRIPTION.
""" """
filterIndex = {'nuv': 0, 'u': 1, 'g': 2, 'r': 3, 'i': 4, 'z': 5, 'y': 6} filterIndex = {'nuv': 0, 'u': 1, 'g': 2,
'r': 3, 'i': 4, 'z': 5, 'y': 6}
sat = (ctypes.c_double*3)() sat = (ctypes.c_double*3)()
sat[:] = self.sat sat[:] = self.sat
...@@ -349,6 +303,7 @@ class StrayLight(object): ...@@ -349,6 +303,7 @@ class StrayLight(object):
return max(band_star_e1, band_star_e2) return max(band_star_e1, band_star_e2)
############################################################################### ###############################################################################
def caculateEarthShineFilter(self, filter='i'): def caculateEarthShineFilter(self, filter='i'):
""" """
# #
...@@ -364,7 +319,8 @@ class StrayLight(object): ...@@ -364,7 +319,8 @@ class StrayLight(object):
DESCRIPTION. DESCRIPTION.
""" """
filterIndex = {'nuv': 0, 'u': 1, 'g': 2, 'r': 3, 'i': 4, 'z': 5, 'y': 6} filterIndex = {'nuv': 0, 'u': 1, 'g': 2,
'r': 3, 'i': 4, 'z': 5, 'y': 6}
sat = (ctypes.c_double*3)() sat = (ctypes.c_double*3)()
sat[:] = self.sat sat[:] = self.sat
...@@ -391,29 +347,6 @@ class StrayLight(object): ...@@ -391,29 +347,6 @@ class StrayLight(object):
############################################################################### ###############################################################################
# def str2time(strTime):
# """
# Parameters
# ----------
# strTime : TYPE
# DESCRIPTION.
# Returns
# -------
# TYPE
# DESCRIPTION.
# """
# if len(strTime) > 20: # 暂时未用到
# msec = int(float('0.'+strTime[20:])*1000000) # 微秒
# str2 = strTime[0:19]+' '+str(msec)
# return datetime.strptime(str2, '%Y %m %d %H %M %S %f')
# datetime类转mjd
##########################################################################
def time2mjd(dateT): def time2mjd(dateT):
""" """
...@@ -436,6 +369,8 @@ def time2mjd(dateT): ...@@ -436,6 +369,8 @@ def time2mjd(dateT):
return mjd+mjd_s/86400.0 return mjd+mjd_s/86400.0
########################################################################### ###########################################################################
def time2jd(dateT): def time2jd(dateT):
""" """
...@@ -460,6 +395,7 @@ def time2jd(dateT): ...@@ -460,6 +395,7 @@ def time2jd(dateT):
# mjd转datetime类 # mjd转datetime类
############################################################################# #############################################################################
def mjd2time(mjd): def mjd2time(mjd):
""" """
...@@ -479,6 +415,8 @@ def mjd2time(mjd): ...@@ -479,6 +415,8 @@ def mjd2time(mjd):
return t0+timedelta(days=mjd) return t0+timedelta(days=mjd)
############################################################################# #############################################################################
def jd2time(jd): def jd2time(jd):
""" """
...@@ -500,6 +438,8 @@ def jd2time(jd): ...@@ -500,6 +438,8 @@ def jd2time(jd):
# mjd和jd互转 # mjd和jd互转
############################################################################ ############################################################################
def mjd2jd(mjd): def mjd2jd(mjd):
""" """
...@@ -518,6 +458,8 @@ def mjd2jd(mjd): ...@@ -518,6 +458,8 @@ def mjd2jd(mjd):
return mjd+2400000.5 return mjd+2400000.5
########################################################################### ###########################################################################
def jd2mjd(jd): def jd2mjd(jd):
""" """
...@@ -536,6 +478,8 @@ def jd2mjd(jd): ...@@ -536,6 +478,8 @@ def jd2mjd(jd):
return jd-2400000.5 return jd-2400000.5
########################################################################## ##########################################################################
def dt2hmd(dt): def dt2hmd(dt):
""" """
...@@ -573,36 +517,9 @@ def dt2hmd(dt): ...@@ -573,36 +517,9 @@ def dt2hmd(dt):
return str_h+':'+str_m+':'+str_d return str_h+':'+str_m+':'+str_d
############################################################################### ###############################################################################
# def float2char(a, z):
# """
# Parameters
# ----------
# a : TYPE
# DESCRIPTION.
# z : TYPE
# DESCRIPTION.
# Returns
# -------
# TYPE
# DESCRIPTION.
# """
# # a is input float value
# # transfer float a to chars and save z bis
# b = str(a)
# n = len(b)
# if z >= n:
# return b
# else:
# return b[:z]
############################################################################### ###############################################################################
def deg2HMS(ra0, dec0): def deg2HMS(ra0, dec0):
""" """
...@@ -715,6 +632,8 @@ def beta_angle(x_sat, y_sat, z_sat, vx_sat, vy_sat, vz_sat, ra_obj, dec_obj): ...@@ -715,6 +632,8 @@ def beta_angle(x_sat, y_sat, z_sat, vx_sat, vy_sat, vz_sat, ra_obj, dec_obj):
############################################################################### ###############################################################################
########################################################## ##########################################################
def LSR_velocity(ra, dec, velocity, Obstime): def LSR_velocity(ra, dec, velocity, Obstime):
""" """
...@@ -822,42 +741,7 @@ def centroid(data): ...@@ -822,42 +741,7 @@ def centroid(data):
return float(cx), float(cy) return float(cx), float(cy)
############################################################################### ###############################################################################
# def centroidN(data):
# """
# Parameters
# ----------
# data : TYPE
# DESCRIPTION.
# Returns
# -------
# cx : TYPE
# DESCRIPTION.
# cy : TYPE
# DESCRIPTION.
# """
# '''
# 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
####################################################################
def SpecCube2photon(data, wave): def SpecCube2photon(data, wave):
""" """
...@@ -891,6 +775,8 @@ def SpecCube2photon(data, wave): ...@@ -891,6 +775,8 @@ def SpecCube2photon(data, wave):
###################################################################### ######################################################################
###################################################################### ######################################################################
def mag2flux(starmag, lam): def mag2flux(starmag, lam):
""" """
...@@ -927,6 +813,7 @@ def mag2flux(starmag, lam): ...@@ -927,6 +813,7 @@ def mag2flux(starmag, lam):
############################################################################### ###############################################################################
def flux2photons(flux, lam): def flux2photons(flux, lam):
""" """
...@@ -954,6 +841,8 @@ def flux2photons(flux, lam): ...@@ -954,6 +841,8 @@ def flux2photons(flux, lam):
return photonu return photonu
############################################################################### ###############################################################################
def fftrange(n, dtype=float): def fftrange(n, dtype=float):
""" """
...@@ -975,6 +864,8 @@ def fftrange(n, dtype=float): ...@@ -975,6 +864,8 @@ def fftrange(n, dtype=float):
return np.arange(-n//2, -n//2+n, dtype=dtype) return np.arange(-n//2, -n//2+n, dtype=dtype)
############################################################################### ###############################################################################
def dft2(ary, Q, samples, shift=None): def dft2(ary, Q, samples, shift=None):
""" """
...@@ -1041,6 +932,8 @@ def dft2(ary, Q, samples, shift=None): ...@@ -1041,6 +932,8 @@ def dft2(ary, Q, samples, shift=None):
return out return out
############################################################################## ##############################################################################
############################################################################### ###############################################################################
def anySampledPSFnew(wavefront, pupil, Q, sizeout): def anySampledPSFnew(wavefront, pupil, Q, sizeout):
""" """
...@@ -1103,6 +996,7 @@ def anySampledPSFnew(wavefront, pupil, Q, sizeout): ...@@ -1103,6 +996,7 @@ def anySampledPSFnew(wavefront, pupil, Q, sizeout):
############################################################################## ##############################################################################
def get_dx_dy_blue(wave): def get_dx_dy_blue(wave):
""" """
...@@ -1136,6 +1030,7 @@ def get_dx_dy_blue(wave): ...@@ -1136,6 +1030,7 @@ def get_dx_dy_blue(wave):
return dx, dy return dx, dy
############################################################################## ##############################################################################
def get_dx_dy_red(wave): def get_dx_dy_red(wave):
""" """
...@@ -1168,6 +1063,7 @@ def get_dx_dy_red(wave): ...@@ -1168,6 +1063,7 @@ def get_dx_dy_red(wave):
############################################################################################## ##############################################################################################
def getSpectrum(Spectrum0, lam, sigma): def getSpectrum(Spectrum0, lam, sigma):
""" """
...@@ -1286,11 +1182,10 @@ class IFSsimulator(): ...@@ -1286,11 +1182,10 @@ class IFSsimulator():
self.config.read_file(open(self.configfile)) self.config.read_file(open(self.configfile))
# self.log.info(self.information) self.log.info(self.information)
############################################################################### ###############################################################################
def processConfigs(self): def processConfigs(self):
""" """
Processes configuration information and save the information to a dictionary self.information. Processes configuration information and save the information to a dictionary self.information.
...@@ -1334,30 +1229,18 @@ class IFSsimulator(): ...@@ -1334,30 +1229,18 @@ class IFSsimulator():
self.section, 'save_cosmicrays') self.section, 'save_cosmicrays')
self.sky_noise = self.config.getboolean(self.section, 'sky_noise') self.sky_noise = self.config.getboolean(self.section, 'sky_noise')
self.nonlinearity = self.config.getboolean( self.nonlinearity = self.config.getboolean(
self.section, 'nonlinearity') self.section, 'nonlinearity')
self.flatfieldM = self.config.getboolean( self.flatfieldM = self.config.getboolean(
self.section, 'flatfieldM') self.section, 'flatfieldM')
self.readoutNoise = self.config.getboolean( self.readoutNoise = self.config.getboolean(
self.section, 'readoutnoise') self.section, 'readoutnoise')
self.intscale = True self.intscale = True
###self.config.getboolean(self.section, 'intscale')
###################################################################### ######################################################################
# 3
self.information['variablePSF'] = False self.information['variablePSF'] = False
self.booleans = dict(nonlinearity=self.nonlinearity, self.booleans = dict(nonlinearity=self.nonlinearity,
...@@ -1372,8 +1255,6 @@ class IFSsimulator(): ...@@ -1372,8 +1255,6 @@ class IFSsimulator():
intscale=self.intscale, intscale=self.intscale,
save_cosmicrays=self.save_cosmicrays) save_cosmicrays=self.save_cosmicrays)
#####################################################################
############################################################################ ############################################################################
def _createEmpty(self): def _createEmpty(self):
...@@ -1393,6 +1274,7 @@ class IFSsimulator(): ...@@ -1393,6 +1274,7 @@ class IFSsimulator():
############################################################################## ##############################################################################
############################################################################## ##############################################################################
def zodiacal(self, ra, dec, time): def zodiacal(self, ra, dec, time):
""" """
...@@ -1489,7 +1371,6 @@ class IFSsimulator(): ...@@ -1489,7 +1371,6 @@ class IFSsimulator():
# def smoothingWithChargeDiffusion(self, image, sigma=(0.32, 0.32)): # def smoothingWithChargeDiffusion(self, image, sigma=(0.32, 0.32)):
# """ # """
# Parameters # Parameters
# ---------- # ----------
# image : TYPE # image : TYPE
...@@ -1540,7 +1421,7 @@ class IFSsimulator(): ...@@ -1540,7 +1421,7 @@ class IFSsimulator():
""" """
self.log.info('Reading in cosmic ray information from %s and %s' % (self.information['cosmicraylengths'], self.log.info('Reading in cosmic ray information from %s and %s' % (self.information['cosmicraylengths'],
self.information['cosmicraydistance'])) self.information['cosmicraydistance']))
#### print(self.information) # print(self.information)
crLengths = np.loadtxt( crLengths = np.loadtxt(
self.information['dir_path']+self.information['cosmicraylengths']) self.information['dir_path']+self.information['cosmicraylengths'])
...@@ -1551,88 +1432,6 @@ class IFSsimulator(): ...@@ -1551,88 +1432,6 @@ class IFSsimulator():
cr_v=crDists[:, 0], cr_cde=crDists[:, 1], cr_cden=np.shape(crDists)[0]) cr_v=crDists[:, 0], cr_cde=crDists[:, 1], cr_cden=np.shape(crDists)[0])
############################################################################### ###############################################################################
# def _writeFITSfile(self, image, filename):
# """
# :param image: image array to save
# :type image: ndarray
# :param filename: name of the output file, e.g. file.fits
# :type filename: str
# :return: None
# """
# if os.path.isfile(filename):
# os.remove(filename)
# # create a new FITS file, using HDUList instance
# ofd = fits.HDUList(fits.PrimaryHDU())
# # new image HDU
# hdu = fits.ImageHDU(data=image)
# # update and verify the header
# hdu.header.add_history('Created by IFSsim at %s' %
# datetime.isoformat(datetime.utcnow()))
# hdu.verify('fix')
# ofd.append(hdu)
# # write the actual file
# ofd.writeto(filename)
###############################################################################
# def writeFITSfile(self, data, filename, unsigned16bit=False):
# """
# Writes out a simple FITS file.
# :param data: data to be written
# :type data: ndarray
# :param filename: name of the output file
# :type filename: str
# :param unsigned16bit: whether to scale the data using bzero=32768
# :type unsigned16bit: bool
# :return: None
# """
# if os.path.isfile(filename):
# os.remove(filename)
# # create a new FITS file, using HDUList instance
# ofd = fits.HDUList(fits.PrimaryHDU())
# # new image HDU
# hdu = fits.ImageHDU(data=data)
# # add input keywords to the header
# for key, value in self.information.items():
# # truncate long keys
# if len(key) > 8:
# key = key[:7]
# try:
# hdu.header.set(key.upper(), value)
# except:
# try:
# hdu.header.set(key.upper(), str(value))
# except:
# pass
# # write booleans
# for key, value in self.booleans.items():
# # truncate long keys
# if len(key) > 8:
# key = key[:7]
# hdu.header.set(key.upper(), str(value), 'Boolean Flags')
# # update and verify the header
# hdu.header.add_history(
# 'This is an itermediate data product no the final output!')
# hdu.verify('fix')
# ofd.append(hdu)
# # write the actual file
# ofd.writeto(filename)
############################################################################### ###############################################################################
def configure(self, simnumber): def configure(self, simnumber):
""" """
...@@ -1826,91 +1625,8 @@ class IFSsimulator(): ...@@ -1826,91 +1625,8 @@ class IFSsimulator():
self.maskSlit = maskSlit self.maskSlit = maskSlit
###################################################################### ######################################################################
############################################################################## #############################################################################
# def generateflat(self, ave=1.0, sigma=0.01):
# """
# Parameters
# ----------
# ave : TYPE, optional
# DESCRIPTION. The default is 1.0.
# sigma : TYPE, optional
# DESCRIPTION. The default is 0.01.
# Returns
# -------
# TYPE
# DESCRIPTION.
# TYPE
# DESCRIPTION.
# """
# """
# Creates a flat field image with given properties.
# :return: flat field image
# :rtype: ndarray
# """
# self.log.info('Generating a flat field...')
# self.log.info('The flat field has mean value of 1 and a given fluctuations, usually either 1 or 2 percent defined by sigma= %d...' % sigma)
# np.random.seed(5*self.simnumber)
# self.flat_b = np.random.normal(loc=ave, scale=sigma, size=(2048, 4096))
# np.random.seed(55*self.simnumber)
# self.flat_r = np.random.normal(loc=ave, scale=sigma, size=(3072, 6144))
# s1 = self.flat_b
# hdu1 = fits.PrimaryHDU(s1)
# hdu1.header.set('sigma', sigma)
# dtime = datetime.utcnow().strftime('%Y -%m -%d %H: %M: %S')
# hdu1.header.add_history(
# 'flat image of blue channel is generated on :'+dtime)
# f1 = '../flat_Blue_'+str(sigma)+'.fits'
# fits.writeto(f1, s1, header=hdu1.header, overwrite=True)
# s2 = self.flat_r
# hdu1 = fits.PrimaryHDU(s2)
# hdu1.header.set('sigma', sigma)
# dtime = datetime.utcnow().strftime('%Y -%m -%d %H: %M: %S')
# hdu1.header.add_history(
# 'flat image of red channel is generated on :'+dtime)
# f2 = '../flat_Red_'+str(sigma)+'.fits'
# fits.writeto(f2, s2, header=hdu1.header, overwrite=True)
# return self.flat_b, self.flat_r
##########################################################################
# def addLampFlux(self):
# """
# Returns
# -------
# None.
# """
# """
# Include flux from the calibration source.
# """
# self.image_b += fits.getdata(self.information['flatflux'])
# self.image_r += fits.getdata(self.information['flatflux'])
# self.log.info('Flux from the calibration unit included (%s)' %
# self.information['flatflux'])
#############################################################################
def MakeFlatMatrix(self, img, seed): def MakeFlatMatrix(self, img, seed):
""" """
...@@ -1972,11 +1688,12 @@ class IFSsimulator(): ...@@ -1972,11 +1688,12 @@ class IFSsimulator():
# blue channle # blue channle
b_4 = np.sum(self.image_b[0:1024, 0:2048], axis=0) b_4 = np.sum(self.image_b[0:1024, 0:2048], axis=0)
###### flip right to left of part3 and part2 # flip right to left of part3 and part2
temp_b_3=np.fliplr(self.image_b[0:1024, 2048:4096]) ## left to right temp_b_3 = np.fliplr(self.image_b[0:1024, 2048:4096]) # left to right
b_3 = np.sum(temp_b_3, axis=0) b_3 = np.sum(temp_b_3, axis=0)
temp_b_2=np.fliplr(self.image_b[1024:1024+1024, 2048:2048+2048]) ## left to right temp_b_2 = np.fliplr(
self.image_b[1024:1024+1024, 2048:2048+2048]) # left to right
b_2 = np.sum(temp_b_2, axis=0) b_2 = np.sum(temp_b_2, axis=0)
b_1 = np.sum(self.image_b[1024:1024+1024, 0:2048], axis=0) b_1 = np.sum(self.image_b[1024:1024+1024, 0:2048], axis=0)
...@@ -1993,20 +1710,18 @@ class IFSsimulator(): ...@@ -1993,20 +1710,18 @@ class IFSsimulator():
# red channle # red channle
r_4 = np.sum(self.image_r[0:1536, 0:3072], axis=0) r_4 = np.sum(self.image_r[0:1536, 0:3072], axis=0)
# flip right to left of part3 and part2
###### flip right to left of part3 and part2 temp_r_3 = np.fliplr(self.image_r[0:1536, 3072:6144]) # left to right
temp_r_3=np.fliplr(self.image_r[0:1536, 3072:6144]) ## left to right
r_3 = np.sum(temp_r_3, axis=0) r_3 = np.sum(temp_r_3, axis=0)
temp_r_2=np.fliplr(self.image_r[1536:1536+1536, 3072:6144]) ## left to right temp_r_2 = np.fliplr(
self.image_r[1536:1536+1536, 3072:6144]) # left to right
r_2 = np.sum(temp_r_2, axis=0) r_2 = np.sum(temp_r_2, axis=0)
#r_3 = np.sum(self.image_r[0:1536, 3072:6144], axis=0) #r_3 = np.sum(self.image_r[0:1536, 3072:6144], axis=0)
#r_2 = np.sum(self.image_r[1536:1536+1536, 3072:6144], axis=0) #r_2 = np.sum(self.image_r[1536:1536+1536, 3072:6144], axis=0)
r_1 = np.sum(self.image_r[1536:1536+1536, 0:3072], axis=0) r_1 = np.sum(self.image_r[1536:1536+1536, 0:3072], axis=0)
for k in range(1536+320): for k in range(1536+320):
...@@ -2121,69 +1836,6 @@ class IFSsimulator(): ...@@ -2121,69 +1836,6 @@ class IFSsimulator():
########################################################## ##########################################################
######################################################################### #########################################################################
# def addReadoutTrails(self):
# """
# Returns
# -------
# None.
# """
# """
# Add readout trails resulting from reading out the shutter open.
# Quadrants assumed to be numbered:
# 2 3
# 0 1
# """
# flux_ratio = self.information['readouttime'] / float(
# self.information['bluesize']) / self.information['exptime']
# # make a copy, this will be updated
# data = self.image_b.copy()
# # Amplifier at different positions depending on the quadrant number !
# # left side is 0, 2 and right side is 1, 3 starting from bottom i.e.
# # going clock wise from lower left we have 0, 2, 3, and 1 quadrants.
# if self.information['quadrant'] in (2, 3):
# data = data[::-1, :]
# data_shift = data.copy() * flux_ratio
# size1, size2 = data.shape
# for i in range(1, size2, 1):
# data_shift2 = np.roll(data_shift, i, axis=0)
# data_shift2[:i, :] = 0.0
# data += data_shift2
# if self.information['quadrant'] in (2, 3):
# self.image_b = data[::-1, :]
# else:
# self.image_b = data
# flux_ratio = self.information['readouttime'] / float(
# self.information['redsize']) / self.information['exptime']
# # make a copy, this will be updated
# data = self.image_r.copy()
# # Amplifier at different positions depending on the quadrant number !
# # left side is 0, 2 and right side is 1, 3 starting from bottom i.e.
# # going clock wise from lower left we have 0, 2, 3, and 1 quadrants.
# if self.information['quadrant'] in (2, 3):
# data = data[::-1, :]
# data_shift = data.copy() * flux_ratio
# size1, size2 = data.shape
# for i in range(1, size2, 1):
# data_shift2 = np.roll(data_shift, i, axis=0)
# data_shift2[:i, :] = 0.0
# data += data_shift2
# if self.information['quadrant'] in (2, 3):
# self.image_r = data[::-1, :]
# else:
# self.image_r = data
############################################################################## ##############################################################################
...@@ -2238,57 +1890,8 @@ class IFSsimulator(): ...@@ -2238,57 +1890,8 @@ class IFSsimulator():
############################################################################## ##############################################################################
# def applyCosmicBackground(self):
# """
# Returns
# -------
# None.
# """
# """
# Apply dark the cosmic background. Scales the background with the exposure time.
# Additionally saves the image without noise to a FITS file.
# """
# # add background
# bcgr = self.information['exptime'] * self.information['cosmic_bkgd']
# self.image_b += bcgr
# self.image_r += bcgr
# self.log.info('Added cosmic background = %f' % bcgr)
# if self.cosmicRays:
# self.imagenoCR_b += bcgr
# self.imagenoCR_r += bcgr
##########################################################################
# def applyScatteredLight(self):
# """
# Returns
# -------
# None.
# """
# """
# Adds spatially uniform scattered light to the image.
# """
# sl = self.information['exptime'] * self.information['scattered_light']
# self.image_b += sl
# self.image_r += sl
# self.log.info('Added scattered light = %f' % sl)
############################################################################## ##############################################################################
def applyPoissonNoise(self): def applyPoissonNoise(self):
""" """
...@@ -2322,7 +1925,7 @@ class IFSsimulator(): ...@@ -2322,7 +1925,7 @@ class IFSsimulator():
self.image_r += residual self.image_r += residual
################################################################################################################### ##############################################################################
def applyCosmetics(self): def applyCosmetics(self):
""" """
...@@ -2416,13 +2019,11 @@ class IFSsimulator(): ...@@ -2416,13 +2019,11 @@ class IFSsimulator():
# at this point we can give fake data... # at this point we can give fake data...
cti = CTI.CDM03bidir(self.information, [], log=self.log) cti = CTI.CDM03bidir(self.information, [], log=self.log)
### here we need the right input data # here we need the right input data
self.image_b = cti.applyRadiationDamage(self.image_b.copy( self.image_b = cti.applyRadiationDamage(self.image_b.copy(
).transpose(), iquadrant=self.information['quadrant']).transpose() ).transpose(), iquadrant=self.information['quadrant']).transpose()
self.log.info('Radiation damage added.') self.log.info('Radiation damage added.')
self.log.debug('Starting to apply radiation damage model...') self.log.debug('Starting to apply radiation damage model...')
# at this point we can give fake data... # at this point we can give fake data...
cti = CTI.CDM03bidir(self.information, [], log=self.log) cti = CTI.CDM03bidir(self.information, [], log=self.log)
...@@ -2497,7 +2098,6 @@ class IFSsimulator(): ...@@ -2497,7 +2098,6 @@ class IFSsimulator():
self.image_b[0:1024+overscan, 2418*3:2418*3+2048+prescan+overscan] += np.random.normal( self.image_b[0:1024+overscan, 2418*3:2418*3+2048+prescan+overscan] += np.random.normal(
loc=0.0, scale=self.information['rn2_b'], size=(1344, 2418)) loc=0.0, scale=self.information['rn2_b'], size=(1344, 2418))
self.log.info('readnoise added in blue channel') self.log.info('readnoise added in blue channel')
# red zone 4 , OSH### 1536* 3072 # red zone 4 , OSH### 1536* 3072
...@@ -2878,32 +2478,32 @@ class IFSsimulator(): ...@@ -2878,32 +2478,32 @@ class IFSsimulator():
################################################################################################## ##################################################################################################
def applyImageShift(self): # def applyImageShift(self):
""" # """
Returns # Returns
------- # -------
None. # None.
""" # """
np.random.seed(9*self.simnumber) # np.random.seed(9*self.simnumber)
ud = np.random.random() # Choose a random rotation # ud = np.random.random() # Choose a random rotation
dx = 2 * (ud-0.5) * self.information['shiftmax'] # dx = 2 * (ud-0.5) * self.information['shiftmax']
np.random.seed(99*self.simnumber) # np.random.seed(99*self.simnumber)
ud = np.random.random() # Choose a random rotation # ud = np.random.random() # Choose a random rotation
dy = 2 * (ud-0.5) * self.information['shiftmax'] # dy = 2 * (ud-0.5) * self.information['shiftmax']
self.image_b = ndimage.shift(self.image_b.copy(), [ # self.image_b = ndimage.shift(self.image_b.copy(), [
dy+self.information['shift_b_y'], dx+self.information['shift_b_x']], order=0, mode='nearest') # dy+self.information['shift_b_y'], dx+self.information['shift_b_x']], order=0, mode='nearest')
self.image_r = ndimage.shift(self.image_r.copy(), [ # self.image_r = ndimage.shift(self.image_r.copy(), [
dy+self.information['shift_r_y'], dx+self.information['shift_r_x']], order=0, mode='nearest') # dy+self.information['shift_r_y'], dx+self.information['shift_r_x']], order=0, mode='nearest')
self.log.info('Applied image shifting to g r i channels.') # self.log.info('Applied image shifting to g r i channels.')
self.information['ra'] = dx*self.information['pixel_size'] # self.information['ra'] = dx*self.information['pixel_size']
self.information['dec'] = dy*self.information['pixel_size'] # self.information['dec'] = dy*self.information['pixel_size']
############################################################################## ##############################################################################
...@@ -2911,7 +2511,6 @@ class IFSsimulator(): ...@@ -2911,7 +2511,6 @@ class IFSsimulator():
# def applyImageRotate(self): # def applyImageRotate(self):
# """ # """
# Returns # Returns
# ------- # -------
# None. # None.
...@@ -2939,6 +2538,7 @@ class IFSsimulator(): ...@@ -2939,6 +2538,7 @@ class IFSsimulator():
# 'Applied telescope rotation with angle (in degree)= %f.', angle) # 'Applied telescope rotation with angle (in degree)= %f.', angle)
############################################################################### ###############################################################################
def CCDreadout(self): def CCDreadout(self):
""" """
...@@ -3302,7 +2902,7 @@ class IFSsimulator(): ...@@ -3302,7 +2902,7 @@ class IFSsimulator():
hdu_b.header['NAXIS2'] = (np.int32(9672), 'length of second data axis') hdu_b.header['NAXIS2'] = (np.int32(9672), 'length of second data axis')
if self.information['exposuretimes']>1: if self.information['exposuretimes'] > 1:
hdu_b.header['NAXIS3'] = ( hdu_b.header['NAXIS3'] = (
np.int32(self.information['exposuretimes']), 'length of third data axis') np.int32(self.information['exposuretimes']), 'length of third data axis')
...@@ -3682,7 +3282,7 @@ class IFSsimulator(): ...@@ -3682,7 +3282,7 @@ class IFSsimulator():
hdu_r.header['NAXIS2'] = ( hdu_r.header['NAXIS2'] = (
np.int32(13768), 'length of second data axis') np.int32(13768), 'length of second data axis')
if self.information['exposuretimes']>1: if self.information['exposuretimes'] > 1:
hdu_r.header['NAXIS3'] = ( hdu_r.header['NAXIS3'] = (
np.int32(self.information['exposuretimes']), 'length of third data axis') np.int32(self.information['exposuretimes']), 'length of third data axis')
...@@ -3993,7 +3593,6 @@ class IFSsimulator(): ...@@ -3993,7 +3593,6 @@ class IFSsimulator():
################################################################################### ###################################################################################
################################################################################## ##################################################################################
def CalskyNoise(self, lam): def CalskyNoise(self, lam):
""" """
...@@ -4217,7 +3816,6 @@ class IFSsimulator(): ...@@ -4217,7 +3816,6 @@ class IFSsimulator():
####################################################### #######################################################
# energy=energy+img0[500-302:500+320,500-320:500+320].sum() ## calculate the slice image energy;
# CCD quantum efficiency # CCD quantum efficiency
CCD_Qe_lam = np.interp(lam, self.CCD_Qe[:, 0], self.CCD_Qe[:, 1]) CCD_Qe_lam = np.interp(lam, self.CCD_Qe[:, 0], self.CCD_Qe[:, 1])
...@@ -4310,11 +3908,6 @@ class IFSsimulator(): ...@@ -4310,11 +3908,6 @@ class IFSsimulator():
photons_blue.addTo(blue_img) photons_blue.addTo(blue_img)
# zblue=blue_img.array
# energy_blue=energy_blue+sum(photons_blue.flux)
width_blue = width_blue+Width_lambda/32.0
# fits.writeto('blueImg.fits',blue_img.array,overwrite=True) # fits.writeto('blueImg.fits',blue_img.array,overwrite=True)
...@@ -4901,7 +4494,6 @@ class IFSsimulator(): ...@@ -4901,7 +4494,6 @@ class IFSsimulator():
energy_blue = energy_blue+sum(photons_blue.flux) energy_blue = energy_blue+sum(photons_blue.flux)
width_blue = width_blue+deltalam/32.0
################# #################
...@@ -5144,18 +4736,6 @@ class IFSsimulator(): ...@@ -5144,18 +4736,6 @@ class IFSsimulator():
conv[conv < 0.0] = 0.0 conv[conv < 0.0] = 0.0
img0 = conv img0 = conv
# temp=img0*maskSlice['31']*maskSlit['31']
# print(image0.array.sum())
# print(img1.sum())
# fits.writeto('img.fits',image0.array,overwrite=True)
# fits.writeto('convimg.fits',img0,overwrite=True)
#######
#######################################################
# energy=energy+img0[50-32:50+32,50-32:50+32].sum() ## calculate the slice image energy;
# CCD quantum efficiency # CCD quantum efficiency
CCD_Qe_lam = np.interp(lam, self.CCD_Qe[:, 0], self.CCD_Qe[:, 1]) CCD_Qe_lam = np.interp(lam, self.CCD_Qe[:, 0], self.CCD_Qe[:, 1])
...@@ -5264,7 +4844,6 @@ class IFSsimulator(): ...@@ -5264,7 +4844,6 @@ class IFSsimulator():
photons_blue.y = photons_blue.y/self.pixelscale + \ photons_blue.y = photons_blue.y/self.pixelscale + \
dy_blue+self.slice_blue['py'][k] dy_blue+self.slice_blue['py'][k]
#photons_blue.y =2048-photons_blue.y #photons_blue.y =2048-photons_blue.y
blue_sensor.accumulate(photons_blue, blue_img) blue_sensor.accumulate(photons_blue, blue_img)
...@@ -5288,7 +4867,6 @@ class IFSsimulator(): ...@@ -5288,7 +4867,6 @@ class IFSsimulator():
photons_red.y = photons_red.y/self.pixelscale + \ photons_red.y = photons_red.y/self.pixelscale + \
dy_red+self.slice_red['py'][k] dy_red+self.slice_red['py'][k]
#photons_red.y =3072-photons_red.y #photons_red.y =3072-photons_red.y
red_sensor.accumulate(photons_red, red_img) red_sensor.accumulate(photons_red, red_img)
...@@ -5370,19 +4948,18 @@ class IFSsimulator(): ...@@ -5370,19 +4948,18 @@ class IFSsimulator():
:return: None :return: None
""" """
self.source = sourcein self.source = sourcein
self.simnumber = simnumber self.simnumber = simnumber
self.configure(simnumber) # print the configfile name and path; self.configure(simnumber) # print the configfile name and path;
self.debug=self.information['debug'] self.debug = self.information['debug']
self.dt = datetime.utcnow() self.dt = datetime.utcnow()
if self.source == 'SCI' or self.source == 'COMP': if self.source == 'SCI' or self.source == 'COMP':
if simnumber <= 50 and simnumber>20: if simnumber <= 50 and simnumber > 20:
self.information['exptime'] = 300 self.information['exptime'] = 300
if simnumber <= 100 and simnumber > 50: if simnumber <= 100 and simnumber > 50:
...@@ -5536,12 +5113,11 @@ class IFSsimulator(): ...@@ -5536,12 +5113,11 @@ class IFSsimulator():
exptimes = self.information['exposuretimes'] exptimes = self.information['exposuretimes']
else: else:
exptimes = 1 exptimes = 1
self.information['exposuretimes'] =1 self.information['exposuretimes'] = 1
# if self.debug: # if self.debug:
# exptimes=1 # exptimes=1
if exptimes>1: if exptimes > 1:
self.image_b_N = np.zeros((1344, 9672, exptimes)) self.image_b_N = np.zeros((1344, 9672, exptimes))
self.image_r_N = np.zeros((1856, 13768, exptimes)) self.image_r_N = np.zeros((1856, 13768, exptimes))
...@@ -5550,7 +5126,6 @@ class IFSsimulator(): ...@@ -5550,7 +5126,6 @@ class IFSsimulator():
self.image_r_N = np.zeros((1856, 13768)) self.image_r_N = np.zeros((1856, 13768))
for idk in range(exptimes): for idk in range(exptimes):
############ add some effect to images ####################################### ############ add some effect to images #######################################
...@@ -5618,8 +5193,7 @@ class IFSsimulator(): ...@@ -5618,8 +5193,7 @@ class IFSsimulator():
self.discretise() self.discretise()
if exptimes > 1:
if exptimes>1:
self.image_b_N[:, :, idk] = self.image_b[:, :].copy() self.image_b_N[:, :, idk] = self.image_b[:, :].copy()
self.image_r_N[:, :, idk] = self.image_r[:, :].copy() self.image_r_N[:, :, idk] = self.image_r[:, :].copy()
...@@ -5627,7 +5201,6 @@ class IFSsimulator(): ...@@ -5627,7 +5201,6 @@ class IFSsimulator():
self.image_b_N[:, :] = self.image_b[:, :].copy() self.image_b_N[:, :] = self.image_b[:, :].copy()
self.image_r_N[:, :] = self.image_r[:, :].copy() self.image_r_N[:, :] = self.image_r[:, :].copy()
self.writeOutputs(simnumber) self.writeOutputs(simnumber)
################################################# #################################################
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment