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

Delete csst_ifs_sim.py

parent 1525b803
Hide whitespace changes
Inline Side-by-side
csst_ifs_sim/csst_ifs_sim.py deleted 100644 → 0
View file @ 1525b803
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
@author: yan
"""
"""
The CSST IFS Image Simulator
=============================================
This file contains an image simulator for the CSST IFS.
The approximate sequence of events in the simulator is as follows:
#. Read in a configuration file, which defines for example,
detector characteristics (bias, dark and readout noise, gain,
plate scale and pixel scale, oversampling factor, exposure time etc.).
#. Read in another file containing charge trap definitions (for CTI modelling).
#. Read in a file defining the cosmic rays (trail lengths and cumulative distributions).
#. Load the wavefront aberration data used to calculate PSF with defined wavelength and field of view.
#. Apply calibration unit flux to mimic flat field exposures [optional].
#. Apply a multiplicative flat-field map to emulate pixel-to-pixel non-uniformity [optional].
#. Add cosmic ray tracks onto the CCD with random positions but known distribution [optional].
#. Apply detector charge bleeding in column direction [optional].
#. Add constant dark current and background light from Zodiacal light [optional].
#. Add photon (Poisson) noise [optional]
#. Add cosmetic defects from an input file [optional].
#. Add overscan regions in the serial direction [optional].
#. Apply the CDM03 radiation damage model [optional].
#. Apply non-linearity model to the pixel data [optional].
#. Add readout noise selected from a Gaussian distribution [optional].
#. Convert from electrons to ADUs using a given gain factor.
#. Add a given bias level and discretise the counts (the output is going to be in 16bit unsigned integers).
#. Finally the simulated image is converted to a FITS file.
:version: 2023.04
Future Work
-----------
.. todo::
#.
#.
#. ......
Contact Information: zhaojunyan@shao.ac.cn
-------
"""
from scipy.integrate import simps
import pandas as pd
from datetime import datetime, timedelta
import julian
from astropy.coordinates import get_sun
from astropy.time import Time
from astropy import units as u
from astropy.coordinates import SkyCoord, EarthLocation
from scipy import interpolate
import numpy as np
import scipy.io as sio
from astropy.io import fits
from scipy.signal import fftconvolve
import galsim
from scipy import ndimage
import cmath
import os, sys
import configparser as ConfigParser
from optparse import OptionParser
sys.path.append('../')
from CTI import CTI
from support import logger as lg
from support import cosmicrays
from support import IFSinstrumentModel
###############################################################################
def str2time(strTime):
if len(strTime)>20:#
msec=int(float('0.'+strTime[20:])*1000000) #
str2=strTime[0:19]+' '+str(msec)
return datetime.strptime(str2,'%Y %m %d %H %M %S %f')
#datetime to mjd
def time2mjd(dateT):
t0=datetime(1858,11,17,0,0,0,0)#
mjd=(dateT-t0).days
mjd_s=dateT.hour*3600.0+dateT.minute*60.0+dateT.second+dateT.microsecond/1000000.0
return mjd+mjd_s/86400.0
def time2jd(dateT):
t0=datetime(1858,11,17,0,0,0,0)#
mjd=(dateT-t0).days
mjd_s=dateT.hour*3600.0+dateT.minute*60.0+dateT.second+dateT.microsecond/1000000.0
return mjd+mjd_s/86400.0++2400000.5
#mjd to datetime
def mjd2time(mjd):
t0=datetime(1858,11,17,0,0,0,0)#
return t0+timedelta(days=mjd)
def jd2time(jd):
mjd=jd2mjd(jd)
return mjd2time(mjd)
#mjd to jd
def mjd2jd(mjd):
return mjd+2400000.5
def jd2mjd(jd):
return jd-2400000.5
def dt2hmd(dt):
## dt is datetime
hour=dt.hour
minute=dt.minute
second=dt.second
if hour<10:
str_h='0'+str(hour)
else:
str_h=str(hour)
if minute<10:
str_m='0'+str(minute)
else:
str_m=str(minute)
if second<10:
str_d='0'+str(second+dt.microsecond*1e-6)
else:
str_d=str(second+dt.microsecond*1e-6)
return str_h+':'+str_m+':'+str_d
###############################################################################
###############################################################################
def deg2HMS(ra, dec, rou=False):
'''convert deg to ra's HMS or dec's DHS'''
RA, DEC, rs, ds = '000000', '000000', '', ''
if dec:
if str(dec)[0] == '-':
ds, dec = '-', abs(dec)
deg = int(dec)
decM = abs(int((dec-deg)*60))
if rou:
decS = round((abs((dec-deg)*60)-decM)*60,1)
else:
decS = int((abs((dec-deg)*60)-decM)*60)
if deg ==0:
deg ="00"
elif deg <10:
deg = "0%s"%deg
if decM ==0:
decM ="00"
elif decM <10:
decM = "0%s"%decM
if decS ==0:
decS ="00"
elif decS <10:
decS = "0%s"%decS
DEC = '{0}{1}{2}{3}'.format(ds, deg, decM, decS)
if ra:
if str(ra)[0] == '-':
rs, ra = '-', abs(ra)
raH = int(ra/15)
raM = int(((ra/15)-raH)*60)
if rou:
raS = round(((((ra/15)-raH)*60)-raM)*60,1)
else:
raS = int(((((ra/15)-raH)*60)-raM)*60)
if raH ==0:
raH = "00"
elif raH <10:
raH = "0%s"%raH
if raM ==0:
raM = "00"
elif raM <10:
raM = "0%s"%raM
if raS ==0:
raS = "00"
elif raS <10:
raS = "0%s"%raS
RA = '{0}{1}{2}{3}'.format(rs, raH, raM, raS)
if ds=='-':
return RA+DEC
else:
return RA+'+'+DEC
###############################################################################
###########################################################
def beta_angle( x_sat, y_sat, z_sat, vx_sat, vy_sat, vz_sat, ra_obj, dec_obj):
# get the vector for next motion
x_sat2 = x_sat + vx_sat
y_sat2 = y_sat + vy_sat
z_sat2 = z_sat + vz_sat
vector1 = np.stack((x_sat, y_sat, z_sat), axis=0)
vector2 = np.stack((x_sat2, y_sat2, z_sat2), axis=0)
vector_normal = np.cross(vector1, vector2)
location_normal = EarthLocation.from_geocentric(vector_normal[0],
vector_normal[1],
vector_normal[2],
unit=u.km)
radec_normal = SkyCoord(ra=location_normal.lon,
dec=location_normal.lat,
frame='gcrs')
lb_normal = radec_normal.transform_to('geocentrictrueecliptic')
radec_sun = SkyCoord(ra=ra_obj*u.degree, dec=dec_obj*u.degree, frame='gcrs')
lb_sun = radec_sun.transform_to('geocentrictrueecliptic')
# get the angle between normal and solar
angle = 90 - lb_normal.separation(lb_sun).degree
return angle
##########################################################
##########################################################
def LSR_velocity(ra,dec,velocity,Obstime):
local = EarthLocation.from_geodetic(lat=22.349368*u.deg, lon=113.584068*u.deg, height=10*u.m)
### convert ra and dec to
source = SkyCoord(ra*u.deg, dec*u.deg,frame='icrs', unit=(u.hourangle,u.deg) )
l = source.galactic.l.deg
b = source.galactic.b.deg
c = SkyCoord(l=l*u.degree, b=b*u.degree, frame='galactic')
c_icrs = c.transform_to('icrs')
barycorr = c_icrs.radial_velocity_correction(obstime=Time(Obstime), location=local)
velocity = velocity + barycorr.value/1000
#print(barycorr.value/1000)
l = l * np.pi / 180
b = b * np.pi / 180
return velocity + 9 * np.cos(l) * np.cos(b) + 12 * np.sin(l) * np.cos(b) + 7 * np.sin(b)
###############################################################################
def rotation_yan(x0,y0,x1,y1,angle):
#% point A (x0,y0)
#% point B(x1,y1)
#% roate angle ,point B rotate with point A
alpha=angle/180*3.1415926 #% in radians
x2=(x1-x0)*np.cos(alpha)-(y1-y0)*np.sin(alpha)+x0
y2=(x1-x0)*np.sin(alpha)+(y1-y0)*np.cos(alpha)+y0
return x2, y2
######################################################################
def centroid(data):
h,w = np.shape(data)
x = np.arange(0,w)
y = np.arange(0,h)
x1=np.ones((1,h))
y1=np.ones((w,1))
cx=(np.dot(np.dot(x1,data),y)) /(np.dot(np.dot(x1,data),y1))
cy=(np.dot(np.dot(x ,data),y1)) /(np.dot(np.dot(x1,data),y1))
return np.float64(cx),np.float64(cy)
####################################################################
def SpecCube2photon(data,wave):
# calcutle photons from original spec cube data,
# data: input data, two-dimentional, unit : 1e-17 erg/s/A/cm^2
# wave: the relative wavefront, unit in nm;
planckh= 6.62620*10**-27 # % erg s;
cc=2.99792458*10**17 # in nm/s
# fov2=0.01 # in arcsec^2
# telarea=3.1415926*100*100 # in cm^2，
fluxlam=1e-17*data # convert original unit to unit of erg/s/A/cm^2
lam=wave # wave in nm ;;
ephoton=planckh*cc/lam # in erg/photon,
Nphoton =fluxlam/ephoton # in unit of photons/cm2/s/A
return Nphoton
######################################################################
######################################################################
def zero_pad(image, shape, position='corner'):
"""
Extends image to a certain size with zeros
Parameters
----------
image: real 2d `np.ndarray`
Input image
shape: tuple of int
Desired output shape of the image
position : str, optional
The position of the input image in the output one:
* 'corner'
top-left corner (default)
* 'center'
centered
Returns
-------
padded_img: real `np.ndarray`
The zero-padded image
"""
shape = np.asarray(shape, dtype=int)
imshape = np.asarray(image.shape, dtype=int)
if np.alltrue(imshape == shape):
return image
if np.any(shape <= 0):
raise ValueError("ZERO_PAD: null or negative shape given")
dshape = shape - imshape
if np.any(dshape < 0):
raise ValueError("ZERO_PAD: target size smaller than source one")
pad_img = np.zeros(shape, dtype=image.dtype)
idx, idy = np.indices(imshape)
if position == 'center':
if np.any(dshape % 2 != 0):
raise ValueError("ZERO_PAD: source and target shapes "
"have different parity.")
offx, offy = dshape // 2
else:
offx, offy = (0, 0)
pad_img[idx + offx, idy + offy] = image
return pad_img
def anySampledPSF(wavefront, pupil , Q , sizeout):
'''
written on 2020.12.04 by Yan @Shao
% wavefront sampled in the united circle;
% pupil function samped as wavefront;
% sample ratio Q=lambda/D/pixel;
% lambda is wavelength;
% D is diameter;
% pixel is in arcsec;
% or sample ratio Q=lambda*f/D/pixelsize
% f is focal length;
% pixelsize is the actural size of the detector;
% make sure all the varia have the same unit;
% the returned PSF has sum value of 1
'''
m,n = np.shape(wavefront)
phase=2.0*np.pi*wavefront #% the phase of the input wavefront;
####generalized pupil function of channel1;
if Q>=1:
pk1=zero_pad(pupil*np.exp(cmath.sqrt(-1)*(phase)),( int(np.floor(m*Q)),int(np.floor(m*Q))), position='corner')
psf=abs(np.fft.fftshift(np.fft.fft2(pk1)))**2
#psf=psf/psf.sum()
else: #%%% case: Q<1
# %fft method
if Q>=0.5: #% in this case Q<1 and Q>=0.5
pk1=zero_pad(pupil*np.exp(cmath.sqrt(-1)*(phase)), (int(2*np.floor(m*Q)),int(2*np.floor(m*Q))), position='corner')
mypsf=np.fft.fft2(pk1);
t=mypsf[0::2,0::2];
psf=abs(np.fft.fftshift(t))**2
#psf=mypsf3/mypsf3.sum();
else:
if Q>=0.25:
pk1=zero_pad(pupil*np.exp(cmath.sqrt(-1)*(phase)),( int(4*np.floor(m*Q)),int(4*np.floor(m*Q))), position='corner')
mypsf=np.fft.fft2(pk1);
t=mypsf[0::4,0::4];
psf=abs(np.fft.fftshift(t))**2
#psf=mypsf3/mypsf3.sum();
else:
print('---- Q<0.25 , stop running')
sys.exit(0)
########################################################################################
mm,nn=np.shape(psf)
if np.mod(sizeout,2)==0:
Nx=sizeout/2-1
else:
Nx=(sizeout+1)/2-1
if max(mm,nn)<sizeout:
psfout=np.zeros((sizeout,sizeout))
[cx,cy]=np.where(psf==np.max(psf))
s=np.ceil(min(mm,nn)/2)-abs(cx-cy)-1
Nx=int(Nx)
cx=int(cx)
cy=int(cy)
s=int(s)
psfout[Nx-s:Nx+s-1,Nx-s:Nx+s-1]=psf[cx-s:cx+s-1,cy-s:cy+s-1]
else:
psfout=np.zeros((sizeout,sizeout));
[cx,cy]=np.where(psf==np.max(psf))
s1=int(np.ceil(min(mm,nn)/2)-abs(cx-cy)-1)
s2=int(sizeout/2)
s=min(s1,s2)
s=s-1
Nx=int(Nx)
cx=int(cx)
cy=int(cy)
psfout[Nx-s:Nx+s-1,Nx-s:Nx+s-1]=psf[cx-s:cx+s-1,cy-s:cy+s-1]
####
cx,cy=centroid(psfout)
Nt=sizeout
#psf=ndimage.shift(psfout,[Nt/2-cy, Nt/2-cx],order=1, mode='nearest' ) ## for my fft method
psf=ndimage.shift(psfout,[Nt/2-cy-1, Nt/2-cx-1],order=1, mode='nearest' ) #for convolve method
psf=psf/psf.sum()
return psf
##############################################################################
def get_dx_dy_blue(wave):
#wave is the wavelength in nm;
# dx is in dispersion direction, dy is in vertical direction;
dydl=np.array([-423.256, 0.001, 0.00075, 0.0000078, -0.0000000000007, 0.0, 0.0] ) #色散方向
dxdl=0.2*np.array([-9.1519, -1.00000000e-06, 3.50000000e-08, -5.00000000e-09, -1.70000000e-11, 4.00949787e-12, -6.16873452e-15]) #垂直方向
dx=0.0
dy=0.0
for i in range(len(dxdl)):
dx=dx+dxdl[i]*wave**(i)
dy=dy+dydl[i]*wave**(i)
return dx,dy
def get_dx_dy_red(wave):
#wave is the wavelength in nm;
dydl=np.array([-640.0239901372472, 0.0018, 0.00048, 0.0000028, -0.0000000000007, 0.0 , 0.0] ) #色散方向
dxdl=0.00325*np.array([-1638.8, 4.0e-2, 5.500e-3, -5.2e-10, 1.7000e-10, 7.1e-13, -5.16e-15]) #垂直方向
dx=0.0
dy=0.0
for i in range(len(dxdl)):
dx=dx+dxdl[i]*wave**(i)
dy=dy+dydl[i]*wave**(i)
return dx,dy
##############################################################################################
def getSpectrum(Spectrum0,lam,sigma):
#
wave=Spectrum0[:,0] ## wavelength；
Qt =Spectrum0[:,1] ##intensity；
d=abs(lam-wave);
if min(d)>3*sigma:
SpmOut=1.0/20000
else:
column=np.where(d==min(d))
SpmOut=Qt[column[0]]*np.exp(-(lam-wave[column[0]])**2/2/sigma**2)
return SpmOut
#####################################################################################
####################################################################################################################
def processArgs(printHelp=False):
"""
Processes command line arguments.
"""
parser = OptionParser()
parser.add_option('-c', '--configfile', dest='configfile',
help="Name of the configuration file", metavar="string")
parser.add_option('-s', '--section', dest='section',
help="Name of the section of the config file [SCIENCE]", metavar="string")
parser.add_option('-q', '--quadrant', dest='quadrant', help='CCD quadrant to simulate [0, 1, 2, 3]',
metavar='int')
parser.add_option('-x', '--xCCD', dest='xCCD', help='CCD number in X-direction within the FPA matrix',
metavar='int')
parser.add_option('-y', '--yCCD', dest='yCCD', help='CCD number in Y-direction within the FPA matrix',
metavar='int')
parser.add_option('-d', '--debug', dest='debug', action='store_true',
help='Debugging mode on')
parser.add_option('-t', '--test', dest='test', action='store_true',
help='Run unittest')
parser.add_option('-f', '--fixed', dest='fixed', help='Use a fixed seed for the random number generators',
metavar='int')
if printHelp:
parser.print_help()
else:
return parser.parse_args()
##############################################################################################
class IFSsimulator():
"""
CSST IFS Simulator
The image that is being build is in::
self.image
self.image_b blue channel
self.image_r red channel
:param opts: OptionParser instance
:type opts: OptionParser instance
"""
def __init__(self, opts):
"""
Class Constructor.
:param opts: OptionParser instance
:type opts: OptionParser instance
"""
####################################
self.configfile = opts.configfile
if opts.section is None:
####self.section = 'DEFAULT'
self.section = 'TEST' #####simulation section;
else:
self.section = opts.section
#load instrument model, these values are also stored in the FITS header
self.information = IFSinstrumentModel.IFSinformation()
#update settings with defaults
self.information.update(dict(quadrant=int(0),
ccdx=int(0),
ccdy=int(0),
psfoversampling=1.0,
ccdxgap=1.643,
ccdygap=8.116,
fullwellcapacity=90000,
readouttime=4,
rdose=8.0e9,
coveringfraction=0.5 ))
######################################################################
self.configure() #print the configfile name and path;
self.information.update(dict(
cosmicraylengths =self.information['indata_path']+'/cdf_cr_length.dat',
cosmicraydistance=self.information['indata_path']+'/cdf_cr_total.dat',
parallelTrapfile =self.information['indata_path']+'/cdm_euclid_parallel.dat',
serialTrapfile =self.information['indata_path']+'/cdm_euclid_serial.dat',
cosmeticsfile_b =self.information['indata_path']+'/Cosmetics_b.txt',
cosmeticsfile_r =self.information['indata_path']+'/Cosmetics_r.txt' ))
# public system information
self.pixelscale=0.1 # the pixel size is 0.1 arcsec
self.pixelsize=15 # the pixel physical size is 15 micron
self.Fnum=15.469875 # the f number= focal_length/D;
self.telD=2 # tht telescope size is 2 meter, diamter size;
self.telarea=3.1415926*(self.telD/2)**2*100*100 # telescope square in cm^2
self.fov2=0.01 # the pixel square
self.planckh= 6.62620*1e-27 # % erg s;
self.cc=2.99792458*1e17 # in nm/s
self.light_FWHM=0.175
self.HgArsigma=self.light_FWHM/2.35; ## sigma value of the light source;
#### load system optical and CCD efficiency data
matfn0=self.information['indata_path']+'/TotalQ200923.mat'
da0=sio.loadmat(matfn0)
self.optical_blue_Q=da0['opticalQ_1'] # optical efficiency of blue channel
self.optical_red_Q =da0['opticalQ_2'] # optical efficiency of red channel
self.CCD_Qe=da0['ccdQE'] # red channel
############################################################# load all useful data;
############### load wavefront data;
matfn2=self.information['indata_path']+'/opd_638nm.mat'
da2=sio.loadmat(matfn2)
opd0=da2['opd'] # opd unit is in meter
self.opd0=opd0*1e9 # convert opd to nm
self.pupil=abs(opd0)>0.0
####################
da=fits.open(self.information['indata_path']+'/opd1.fits')
self.opd1=da[0].data
da=fits.open(self.information['indata_path']+'/opd2.fits')
self.opd2=da[0].data
da=fits.open(self.information['indata_path']+'/opd3.fits')
self.opd3=da[0].data
################################################################################
## slice information; the slice is put in vertial direction;
# the silce long size is 64 pixels
slice_blue=dict()
slice_red =dict()
randRedpos=np.array([0.52835362, 1.1105926 , -0.81667794, 0.88860884, -2.78092636,
-0.15810022, -1.56726852, -0.71601855, -1.31768647, 1.73107581,
0.4933876 , 2.83673377, 0.22226286, -0.02634998, 0.35539383,
-0.91989574, -2.44856212, 0.91020484, -3.03097852, -1.11638071,
1.28360669, -0.12521128, -0.3907698 , 0.70183575, 1.00578099,
1.67339662, 0.18067182, -0.56303075, 0.40959616, 1.45518379,
-0.93194744, 0.41492972])
randBluepos=np.array([ 0.97449008, -0.21371406, -1.62513338, -3.06938604, 1.72615283,
0.73333374, 0.80923919, -0.9418576 , -0.16806578, -1.04416631,
2.20155068, -0.0900156 , 0.07597706, 0.76208373, 0.29426592,
-0.89434703, 0.34017826, 1.16936499, 0.10977829, -1.31179539,
-0.50859787, -1.01891651, -0.95791581, -1.53018041, 0.88358827,
0.07837641, -0.86157157, -1.18070438, 0.53970599, 1.4913733 ,
2.10938775, 1.23213412])
#####################
for i in range(32):
if i==0:
slice_blue['py']=np.zeros(32)
slice_blue['px']=np.zeros(32)
slice_red['py']=np.zeros(32)
slice_red['px']=np.zeros(32)
if i<16:
slice_blue['py'][i]=50+randBluepos[i]*4
slice_blue['px'][i]=3.55/0.015*i+166.0
slice_red['py'][i]=50+randBluepos[i]*4
slice_red['px'][i]=3.55/0.015*i+1190.0
else:
slice_blue['py'][i]=50+250+randBluepos[i]*4
slice_blue['px'][i]=3.55/0.015*(i-16)+166.0+118
slice_red['py'][i]=50+250+randBluepos[i]*4
slice_red['px'][i]=3.55/0.015*(i-16)+1190.0+118
#######
self.slice_blue=slice_blue
self.slice_red=slice_red
###############################################################################
maskSlice=dict()
maskSlit=dict()
sizeout=100
for k in range(32):
maskSlice[str(k)]=np.zeros((sizeout,sizeout))
maskSlice[str(k)][2*k+18:2*k+20, int(sizeout/2)-32:int(sizeout/2)+32]=1
maskSlit[str(k)]=np.zeros((sizeout,sizeout))
maskSlit[str(k)][2*k+18:2*k+20, int(sizeout/2)-37:int(sizeout/2)+37]=1
self.maskSlice=maskSlice
self.maskSlit =maskSlit
################################################################################################
############################################################################3
def readConfigs(self):
"""
Reads the config file information using configParser and sets up a logger.
"""
self.config = ConfigParser.RawConfigParser()
self.config.read_file(open(self.configfile))
###############################################################################
def processConfigs(self):
"""
Processes configuration information and save the information to a dictionary self.information.
The configuration file may look as follows::
[TEST]
For explanation of each field, see /data/test.config. Note that if an input field does not exist,
then the values are taken from the default instrument model as described in
support.IFSinstrumentModel.VISinformation(). Any of the defaults can be overwritten by providing
a config file with a correct field name.
"""
#parse options and update the information dictionary
options = self.config.options(self.section)
settings = {}
for option in options:
try:
settings[option] = self.config.getint(self.section, option)
except ValueError:
try:
settings[option] = self.config.getfloat(self.section, option)
except ValueError:
settings[option] = self.config.get(self.section, option)
self.information.update(settings)
self.cosmicRays = self.config.getboolean(self.section, 'cosmicRays')
self.darknoise = self.config.getboolean(self.section, 'darknoise')
self.cosmetics = self.config.getboolean(self.section, 'cosmetics')
self.radiationDamage = self.config.getboolean(self.section, 'radiationDamage')
self.bleeding = self.config.getboolean(self.section, 'bleeding')
self.skyback = self.config.getboolean(self.section, 'skyback')
#these don't need to be in the config file
try:
self.nonlinearity = self.config.getboolean(self.section, 'nonlinearity')
except:
self.nonlinearity = False
try:
self.flatfieldM = self.config.getboolean(self.section, 'flatfieldM')
except:
self.flatfieldM = False
try:
self.readoutNoise = self.config.getboolean(self.section, 'readoutNoise')
except:
self.readoutNoise = True
####################################################################
self.booleans = dict(nonlinearity=self.nonlinearity,
flatfieldM=self.flatfieldM,
cosmicRays=self.cosmicRays,
darknoise=self.darknoise,
cosmetics=self.cosmetics,
radiationDamage=self.radiationDamage,
bleeding=self.bleeding,
skyback =self.skyback)
#####################################################################
now=datetime.now()
result_day=now.strftime("%Y-%m-%d")
self.result_path=self.information['result_path']+'/'+result_day
if os.path.isdir(self.result_path)==False:
os.mkdir(self.result_path)
os.mkdir(self.result_path+'/log_file')
os.mkdir(self.result_path+'/sky_Data')
now=datetime.now()
data_time=now.strftime("%Y-%m-%d-%H-%M-%S")
self.log = lg.setUpLogger(self.result_path+'/log_file/IFS_'+'_'+data_time+'.log')
self.log.info('STARTING A NEW SIMULATION')
self.log.info(self.information)
return
#########################################################################################
#######################################################################
def readCosmicRayInformation(self):
"""
Reads in the cosmic ray track information from two input files.
Stores the information to a dictionary called cr.
"""
self.log.info('Reading in cosmic ray information from %s and %s' % (self.information['cosmicraylengths'],
self.information['cosmicraydistance']))
crLengths = np.loadtxt(self.information['cosmicraylengths'])
crDists = np.loadtxt(self.information['cosmicraydistance'])
self.cr = dict(cr_u=crLengths[:, 0], cr_cdf=crLengths[:, 1], cr_cdfn=np.shape(crLengths)[0],
cr_v=crDists[:, 0], cr_cde=crDists[:, 1], cr_cden=np.shape(crDists)[0])
##############################################################################################
def configure(self):
"""
Configures the simulator with input information and creates and empty array to which the final image will
be build on.
"""
self.readConfigs()
self.processConfigs()
self.log.info('Read in the configuration files and created an empty array')
#################################################################################################################
def MakeFlatMatrix(self, img, seed):
####
ysize, xsize=img.shape
np.random.seed(seed)
r1,r2,r3,r4 = np.random.random(4)
a1 = -0.5 + 0.2*r1
a2 = -0.5 + 0.2*r2
a3 = r3+5
a4 = r4+5
xmin,xmax,ymin,ymax = 0, xsize,0, ysize
Flty, Fltx = np.mgrid[ymin:ymax, xmin:xmax]
np.random.seed(seed)
p1,p2,bg=np.random.poisson(1000, 3)
Fltz = 1e-6*(a1 * (Fltx-p1) ** 2 + a2 * (Flty-p2) ** 2 - a3*Fltx - a4*Flty) + bg*20
FlatMat = Fltz/np.mean(Fltz)
return FlatMat
#########################################################################################################
def addLampFlux(self):
"""
Include flux from the calibration source.
"""
self.image_b += fits.getdata(self.information['flatflux'])
self.image_r += fits.getdata(self.information['flatflux'])
self.log.info('Flux from the calibration unit included (%s)' % self.information['flatflux'])
def applyflatfield(self):
"""
Applies multiplicative flat field to emulate pixel-to-pixel non-uniformity.
Because the pixel-to-pixel non-uniformity effect (i.e. multiplicative) flat fielding takes place
before CTI and other effects, the flat field file must be the same size as the pixels that see
the sky.
"""
###
flat_b=self.MakeFlatMatrix(self.image_b, 100)
flat_r=self.MakeFlatMatrix(self.image_r, 200)
self.image_b *= flat_b
self.image_r *= flat_r
self.log.info('Applied flatfield to images.')
return
########################################################
###############################################################################
def addCosmicRays(self):
"""
Add cosmic rays to the arrays based on a power-law intensity distribution for tracks.
Cosmic ray properties (such as location and angle) are chosen from random Uniform distribution.
For details, see the documentation for the cosmicrays class in the support package.
"""
self.readCosmicRayInformation()
self.cr['exptime'] = self.information['exptime'] #to scale the number of cosmics with exposure time
#cosmic ray image
crImage_b = np.zeros((2048, 4096), dtype=np.float64)
crImage_r = np.zeros((3072, 6144), dtype=np.float64)
#cosmic ray instance
cosmics_b = cosmicrays.cosmicrays(self.log, crImage_b, crInfo=self.cr)
cosmics_r = cosmicrays.cosmicrays(self.log, crImage_r, crInfo=self.cr)
#add cosmic rays up to the covering fraction
#CCD_cr = cosmics.addUpToFraction(self.information['coveringFraction'], limit=None)
CCD_cr_b = cosmics_b.addUpToFraction(self.information['coveringfraction'], limit=None)
CCD_cr_r = cosmics_r.addUpToFraction(self.information['coveringfraction'], limit=None)
#paste the information
self.image_b += CCD_cr_b
self.image_r += CCD_cr_r
#count the covering factor
area_cr_b = np.count_nonzero(CCD_cr_b)
area_cr_r = np.count_nonzero(CCD_cr_r)
#self.log.info('The cosmic ray covering factor is %i pixels ' % area_cr)
self.log.info('The cosmic ray in blue channel covering factor is %i pixels ' % area_cr_b)
self.log.info('The cosmic ray in red channel covering factor is %i pixels ' % area_cr_r)
#########################################################
#########################################################################
#########################################################################
def applyDarkCurrent(self):
"""
Apply dark current. Scales the dark with the exposure time.
Additionally saves the image without noise to a FITS file.
"""
self.log.info('Added dark current to bule and red channel' )
########## blue zone 1
self.image_b[0:1024,0:2048] += self.information['exptime'] * self.information['dark1_b']
########## zone 4 #################
self.image_b[1024:2048,0:2048] += self.information['exptime'] * self.information['dark4_b']
########## zone 2 ###################
self.image_b[0:1024,2048:4096] += self.information['exptime'] * self.information['dark2_b']
########## zone 3
self.image_b[1024:2048,2048:4096]+= self.information['exptime'] * self.information['dark3_b']
########## red zone 1
self.image_r[0:1536, 0:3072] += self.information['exptime'] * self.information['dark1_r']
########## zone 4 #################
self.image_r[1536:3712,0:3072] += self.information['exptime'] * self.information['dark4_r']
########## zone 2 ###################
self.image_r[0:1536,3072:6144] += self.information['exptime'] * self.information['dark2_r']
########## zone 3
self.image_r[1536:3072,3072:6144] += self.information['exptime'] * self.information['dark3_r']
#######################################################################################################3
def applyCosmicBackground(self):
"""
Apply dark the cosmic background. Scales the background with the exposure time.
Additionally saves the image without noise to a FITS file.
"""
#add background
bcgr = self.information['exptime'] * self.information['cosmic_bkgd']
#self.image += bcgr
self.image_b += bcgr
self.image_r += bcgr
self.log.info('Added cosmic background = %f' % bcgr)
if self.cosmicRays:
#self.imagenoCR += bcgr
self.imagenoCR_b += bcgr
self.imagenoCR_r += bcgr
##########################################################################################
##############################################################################
def applyPoissonNoise(self):
"""
Add Poisson noise to the image.
"""
rounded = np.rint(self.image_b) ### round to
residual = self.image_b.copy() - rounded #ugly workaround for multiple rounding operations...
rounded[rounded < 0.0] = 0.0
np.random.seed()
self.image_b = np.random.poisson(rounded).astype(np.float64)
self.log.info('Added Poisson noise on channel blue')
self.image_b += residual
rounded = np.rint(self.image_r) ### round to
residual = self.image_r.copy() - rounded #ugly workaround for multiple rounding operations...
rounded[rounded < 0.0] = 0.0
np.random.seed()
self.image_r = np.random.poisson(rounded).astype(np.float64)
self.log.info('Added Poisson noise on channel red')
self.image_r += residual
###################################################################################################################
def applyCosmetics(self):
"""
Apply cosmetic defects described in the input file.
Warning:: This method does not work if the input file has exactly one line.
"""
cosmetics = np.loadtxt(self.information['cosmeticsfile_b'])
x = np.round(cosmetics[:, 0]).astype(int)
y = np.round(cosmetics[:, 1]).astype(int)
value = cosmetics[:, 2]
cosmetics_b=np.zeros((3712,6784))
cosmetics_r=np.zeros((3712,6784))
self.log.info('Adding cosmetic defects to blue channel:' )
for xc, yc, val in zip(x, y, value):
if 0 <= xc <= 6784 and 0 <= yc <= 3712:
#self.image[yc, xc] = val
self.image_b[yc, xc] = val
cosmetics_b[yc,xc]=val
self.log.info('x=%i, y=%i, value=%f' % (xc, yc, val))
######################################################################################################
cosmetics = np.loadtxt(self.information['cosmeticsfile_r'])
x = np.round(cosmetics[:, 0]).astype(int)
y = np.round(cosmetics[:, 1]).astype(int)
value = cosmetics[:, 2]
self.log.info('Adding cosmetic defects to red channel:' )
for xc, yc, val in zip(x, y, value):
if 0 <= xc <= 6784 and 0 <= yc <= 3712:
#self.image[yc, xc] = val
self.image_r[yc, xc] = val
cosmetics_r[yc,xc]=val
self.log.info('x=%i, y=%i, value=%f' % (xc, yc, val))
################################################################################
def applyRadiationDamage(self):
"""
Applies CDM03 radiation model to the image being constructed.
.. seealso:: Class :`CDM03`
"""
self.log.debug('Starting to apply radiation damage model...')
#at this point we can give fake data...
cti = CTI.CDM03bidir(self.information, [], log=self.log)
#here we need the right input data
self.image_b = cti.applyRadiationDamage(self.image_b.copy().transpose(), iquadrant=self.information['quadrant']).transpose()
self.log.info('Radiation damage added.')
self.log.debug('Starting to apply radiation damage model...')
#at this point we can give fake data...
cti = CTI.CDM03bidir(self.information, [], log=self.log)
#here we need the right input data
self.image_r = cti.applyRadiationDamage(self.image_r.copy().transpose(), iquadrant=self.information['quadrant']).transpose()
self.log.info('Radiation damage added.')
##################################################################################
def applyNonlinearity(self):
"""
Applies a CCD273 non-linearity model to the image being constructed.
"""
self.log.debug('Starting to apply non-linearity model...')
self.image_b = IFSinstrumentModel.CCDnonLinearityModel(self.image_b.copy())
self.log.info('Non-linearity effects included.')
self.log.debug('Starting to apply non-linearity model...')
self.image_r = IFSinstrumentModel.CCDnonLinearityModel(self.image_r.copy())
self.log.info('Non-linearity effects included.')
#####################################################################################
def applyReadoutNoise(self):
"""
Applies readout noise to the image being constructed.
The noise is drawn from a Normal (Gaussian) distribution with average=0.0 and std=readout noise.
"""
self.log.info('readnoise added in blue channel')
########## blue zone 1
np.random.seed()
self.image_b[0:1856,0:3392] += np.random.normal(loc=0.0, scale=self.information['rn1_b'], size=(1856,3392))
########## zone 4 #################
np.random.seed()
self.image_b[1856:3712,0:3392] += np.random.normal(loc=0.0, scale=self.information['rn4_b'], size=(1856,3392))
########## zone 2 ###################
np.random.seed()
self.image_b[0:1856,3392:6784] += np.random.normal(loc=0.0, scale=self.information['rn2_b'], size=(1856,3392))
########## zone 3
np.random.seed()
self.image_b[1856:3712,3392:6784]+= np.random.normal(loc=0.0, scale=self.information['rn3_b'], size=(1856,3392))
############################################################################
self.log.info('readnoise added in blue channel')
########## red zone 1
np.random.seed()
self.image_r[0:1856, 0:3392] += np.random.normal(loc=0.0, scale=self.information['rn1_r'], size=(1856,3392))
########## zone 4 #################
np.random.seed()
self.image_r[1856:3712,0:3392] += np.random.normal(loc=0.0, scale=self.information['rn4_r'], size=(1856,3392))
########## zone 2 ###################
np.random.seed()
self.image_r[0:1856,3392:6784] += np.random.normal(loc=0.0, scale=self.information['rn2_r'], size=(1856,3392))
########## zone 3
np.random.seed()
self.image_r[1856:3712,3392:6784]+= np.random.normal(loc=0.0, scale=self.information['rn3_r'], size=(1856,3392))
##########################################################################################
def electrons2ADU(self):
"""
Convert from electrons to ADUs using the value read from the configuration file.
"""
###############################################################
self.log.info('Converting from electrons to ADUs using a factor of gain' )
########## blue zone 1
self.image_b[0:1856,0:3392] /= self.information['gain1_b']
########## zone 4 #################
self.image_b[1856:3712,0:3392] /= self.information['gain4_b']
########## zone 2 ###################
self.image_b[0:1856,3392:6784] /= self.information['gain2_b']
########## zone 3
self.image_b[1856:3712,3392:6784]/= self.information['gain3_b']
############################################################################
########## red zone 1
self.image_r[0:1856, 0:3392] /= self.information['gain1_r']
########## zone 4 #################
self.image_r[1856:3712,0:3392] /= self.information['gain4_r']
########## zone 2 ###################
self.image_r[0:1856,3392:6784] /= self.information['gain2_r']
########## zone 3
self.image_r[1856:3712,3392:6784] /= self.information['gain3_r']
##########################################################################3
####################################################################################
def applyBias(self):
"""
Adds a bias level to the image being constructed.
The value of bias is read from the configure file and stored
in the information dictionary (key bias).
"""
########## blue zone 1
self.image_b[0:1856,0:3392] += self.information['bias1_b']
########## zone 4 #################
self.image_b[1856:3712,0:3392] += self.information['bias4_b']
########## zone 2 ###################
self.image_b[0:1856,3392:6784] += self.information['bias2_b']
########## zone 3
self.image_b[1856:3712,3392:6784] += self.information['bias3_b']
self.log.info('Bias counts were added to the blue image' )
############################################################################
########## red zone 1
self.image_r[0:1856, 0:3392] += self.information['bias1_r']
########## zone 4 #################
self.image_r[1856:3712,0:3392] += self.information['bias4_r']
########## zone 2 ###################
self.image_r[0:1856,3392:6784] += self.information['bias2_r']
########## zone 3
self.image_r[1856:3712,3392:6784] += self.information['bias3_r']
##########################################################################
self.log.info('Bias counts were added to the red image' )
###############################################################################
###############################################################################
def applyBleeding_yan(self):
"""
Apply bleeding along the CCD columns if the number of electrons in a pixel exceeds the full-well capacity.
Bleeding is modelled in the parallel direction only, because the CCD273s are assumed not to bleed in
serial direction.
:return: None
"""
if self.image_b.max()>self.information['fullwellcapacity']:
self.log.info('Applying column bleeding to blue CCD image...')
#loop over each column, as bleeding is modelled column-wise
for i, column in enumerate(self.image_b.T):
sum = 0.
for j, value in enumerate(column):
#first round - from bottom to top (need to half the bleeding)
overload = value - self.information['fullwellcapacity']
if overload > 0.:
overload /= 2.
#self.image[j, i] -= overload
self.image_b[j, i] -= overload
sum += overload
elif sum > 0.:
if -overload > sum:
overload = -sum
self.image_b[j, i] -= overload
sum += overload
for i, column in enumerate(self.image_b.T):
sum = 0.
for j, value in enumerate(column[::-1]):
#second round - from top to bottom (bleeding was half'd already, so now full)
overload = value - self.information['fullwellcapacity']
if overload > 0.:
#self.image[-j-1, i] -= overload
self.image_b[-j-1, i] -= overload
sum += overload
elif sum > 0.:
if -overload > sum:
overload = -sum
#self.image[-j-1, i] -= overload
self.image_b[-j-1, i] -= overload
sum += overload
print('Applying column bleeding to blue image finished.......')
######################################################################
if self.image_r.max()>self.information['fullwellcapacity']:
self.log.info('Applying column bleeding to red CCD image...')
for i, column in enumerate(self.image_r.T):
sum = 0.
for j, value in enumerate(column):
#first round - from bottom to top (need to half the bleeding)
overload = value - self.information['fullwellcapacity']
if overload > 0.:
overload /= 2.
#self.image[j, i] -= overload
self.image_r[j, i] -= overload
sum += overload
elif sum > 0.:
if -overload > sum:
overload = -sum
#self.image[j, i] -= overload
self.image_r[j, i] -= overload
sum += overload
for i, column in enumerate(self.image_r.T):
sum = 0.
for j, value in enumerate(column[::-1]):
#second round - from top to bottom (bleeding was half'd already, so now full)
overload = value - self.information['fullwellcapacity']
if overload > 0.:
#self.image[-j-1, i] -= overload
self.image_r[-j-1, i] -= overload
sum += overload
elif sum > 0.:
if -overload > sum:
overload = -sum
#self.image[-j-1, i] -= overload
self.image_r[-j-1, i] -= overload
sum += overload
print('Applying column bleeding to red image finished.......')
############################################################################
############################################################################
def discretise(self, max=2**16-1):
"""
Converts a floating point image array (self.image) to an integer array with max values
defined by the argument max.
:param max: maximum value the the integer array may contain [default 65k]
:type max: float
:return: None
"""
#avoid negative numbers in case bias level was not added
self.image_b[self.image_b < 0.0] = 0.
#cut of the values larger than max
self.image_b[self.image_b > max] = max
self.image_b = np.rint(self.image_b).astype(int)
self.log.info('Maximum and total values of the image are %i and %i, respectively' % (np.max(self.image_b),
np.sum(self.image_b)))
#avoid negative numbers in case bias level was not added
self.image_r[self.image_r < 0.0] = 0.
#cut of the values larger than max
self.image_r[self.image_r > max] = max
self.image_r = np.rint(self.image_r).astype(int)
self.log.info('Maximum and total values of the image are %i and %i, respectively' % (np.max(self.image_r),
np.sum(self.image_r)))
##################################################################################################
def applyImageShift(self):
np.random.seed()
ud= np.random.random() # Choose a random rotation
dx= 2* (ud-0.5) * self.information['shiftmax']
np.random.seed()
ud= np.random.random() # Choose a random rotation
dy= 2* (ud-0.5) * self.information['shiftmax']
self.image_b= ndimage.shift(self.image_b.copy(), [dy+self.information['shift_b_y'] , dx+self.information['shift_b_x']], order=0, mode='nearest')
self.image_r= ndimage.shift(self.image_r.copy(), [dy+self.information['shift_r_y'] , dx+self.information['shift_r_x']], order=0, mode='nearest')
self.log.info('Applied image shifting to g r i channels.')
self.information['ra'] = dx*self.information['pixel_size']
self.information['dec'] = dy*self.information['pixel_size']
######################################################################################################33
def applyImageRotate(self ):
np.random.seed()
ud = np.random.random() # Choose a random rotation
angle = 2* (ud-0.5) * self.information['tel_rotmax']
inputimg=self.image_b.copy()
rotimg=ndimage.rotate(inputimg,angle+self.information['rotate_b'],order=1,reshape=False) # here we choose reshape=False, the rotated image will
self.image_b=rotimg
inputimg=self.image_r.copy()
rotimg=ndimage.rotate(inputimg,angle+self.information['rotate_r'],order=1,reshape=False) # here we choose reshape=False, the rotated image will
self.image_r=rotimg
self.information['Tel_rot']=angle
self.log.info('Applied telescope rotation with angle (in degree)= %f.', angle)
###############################################################################
def CCDreadout(self):
imgb=self.image_b.copy()
temp=np.zeros((3712,6784))
########## zone 1
x1=0
x2=x1+1024
y1=0
y2=y1+2048
temp[x1:x2,y1:y2]=imgb[0:1024,0:2048]
########## zone 4 #################
x1=2688
x2=x1+1024
y1=0
y2=y1+2048
temp[x1:x2,y1:y2]=imgb[1024:2048,0:2048]
########## zone 2 ###################
x1=0
x2=x1+1024
y1=6784-2048
y2=y1+2048
temp[x1:x2,y1:y2]=imgb[0:1024,2048:4096]
########## zone 3
x1=2688
x2=x1+1024
y1=6784-2048
y2=y1+2048
temp[x1:x2,y1:y2]=imgb[1024:2048,2048:4096]
self.image_b=temp
##############################################################################
imgr=self.image_r.copy()
temp=np.zeros((3712,6784))
########## zone 1
x1=0
x2=x1+1536
y1=0
y2=y1+3072
temp[x1:x2,y1:y2]=imgr[0:1536, 0:3072]
########## zone 4 #################
x1=2176
x2=x1+1536
y1=0
y2=y1+3072
temp[x1:x2,y1:y2]=imgr[1536:3712,0:3072]
########## zone 2 ###################
x1=0
x2=x1+1536
y1=6784-3072
y2=y1+3072
temp[x1:x2,y1:y2]=imgr[0:1536,3072:6144]
########## zone 3
x1=2176
x2=x1+1536
y1=6784-3072
y2=y1+3072
temp[x1:x2,y1:y2]=imgr[1536:3072,3072:6144]
self.image_r=temp
return
##############################################################################
def writeOutputs(self):
"""
Writes out a FITS file using PyFITS and converts the image array to 16bit unsigned integer as
appropriate for VIS.
Updates header with the input values and flags used during simulation.
"""
## Readout information
self.source='sci'
now=datetime.now()
data_time=now.strftime("%Y-%m-%d %H:%M:%S")
exp_endtime=now.strftime("%Y%m%d%H%M%S")
start=now-timedelta(seconds=self.information['exptime'])
exp_starttime=start.strftime("%Y%m%d%H%M%S")
#write the actual file
obsid=300000000+1
data_time=self.dt.strftime("%Y-%m-%d %H:%M:%S")
exp_starttime=self.dt.strftime("%Y%m%d%H%M%S")
### exposure end time is t2 ;
t2=self.dt+timedelta(seconds=self.information['exptime'])
exp_endtime=t2.strftime("%Y%m%d%H%M%S")
t3=self.dt+timedelta(seconds=self.information['exptime'])+timedelta(seconds=self.information['readouttime'])
filename_b='CSST_IFS_B_'+self.source+'_'+exp_starttime+'_'+exp_endtime+'_'+str(obsid)+'_X_L0_VER_'+self.information['img_ver']+'.fits'
file_b=self.result_path+'/sky_Data/'+filename_b
filename_r='CSST_IFS_R_'+self.source+'_'+exp_starttime+'_'+exp_endtime+'_'+str(obsid)+'_X_L0_VER_'+self.information['img_ver']+'.fits'
file_r=self.result_path+'/sky_Data/'+filename_r
#create a new FITS file, using HDUList instance
ofd_b = fits.PrimaryHDU()
ofd_b.header['GROUPS']=( bool(False), 'always F')
ofd_b.header['DATE'] =( data_time, 'date this file was written' )
ofd_b.header['FILENAME']=(filename_b, ' file name C48 ')
ofd_b.header['OBSTYPE'] =( self.source, 'observation type raw,flt, mask, bias, dark, sci')
ofd_b.header['TELESCOP']=('CSST', 'always CSST')
ofd_b.header['INSTRUME']=( 'IFS', ' ')
ofd_b.header['RADECSYS']=('ICRS', ' always ICRS ')
ofd_b.header['EQUINOX'] =( float(2000.0), 'always 2000.0')
ofd_b.header['FITSCREA']=( '4.2.1', 'FITS create software version')
######### Object information #############
ofd_b.header['OBJECT']=( self.information['name_obj'], 'object name')
ofd_b.header['TARGET']=( (self.information['target']), 'target name, hhmmss+ddmmss')
ofd_b.header['OBJ_RA'] =(np.float64(self.information['ra_obj']) , 'RA of the object in deg')
ofd_b.header['OBJ_DEC']=(np.float64(self.information['dec_obj']) , 'DEC of the object in deg')
ofd_b.header['RA_PNT0']=( np.float64(self.information['ra_pnt0']) , 'RA of the pointing (degrees) at EXPSTART')
ofd_b.header['DEC_PNT0']=(np.float64(self.information['dec_pnt0']) , 'DEC of the pointing (degrees) at EXPSTART')
##############
ofd_b.header['OBSID'] =(str(obsid) , 'observation ID, 3+8bit')
######## Telescope information ###############
# ofd_b.header['COMMENT'] ='=========================================================='
# ofd_b.header['COMMENT'] ='Telescope information'
# ofd_b.header['COMMENT'] ='=========================================================='
ofd_b.header['REFFRAME']=('CSSTGSC-1.0' , 'guide star catalog version')
ofd_b.header['DATE-OBS']=(data_time , 'date of the observation (yyyy-mm-dd hh:mm:ss)')
ofd_b.header['EXPSTART']=(np.float64(time2jd(self.dt)), 'exposure start time')
ofd_b.header['SUNANGL0']=(np.float32(0.0) , 'angle between sun and optical axis at EXPSTART')
ofd_b.header['MOONANG0']=(np.float32(0.0) , 'angle between moon and optical axis at EXPSTART')
ofd_b.header['POS_ANG0']=(np.float64(0.0), 'angle between optical axis and the North Pole at EXPSTART in arcsec')
ofd_b.header['TEL_ALT0']=(np.float64(0.0), 'angle between optical axis and the ground- piston at EXPSTART in deg')
ofd_b.header['HOODSTA0']=(np.float32(0.0) , 'lens hood altitude at EXPSTART')
ofd_b.header['HOODANG0']=(np.float32(0.0), 'lens hood azimuth at EXPSTART')
ofd_b.header['POSI0_X'] =(np.float64(self.information['POSI0_X']) , 'the orbital position of CSST in X direction at EXPSTART')
ofd_b.header['POSI0_Y'] =(np.float64(self.information['POSI0_Y']) , 'the orbital position of CSST in Y direction at EXPSTART')
ofd_b.header['POSI0_Z'] =(np.float64(self.information['POSI0_Z']) , 'the orbital position of CSST in Z direction at EXPSTART')
ofd_b.header['VELO0_X']=( np.float64(self.information['VELO0_X']) , 'the orbital velocity of CSST in X direction at EXPSTART')
ofd_b.header['VELO0_Y']=( np.float64(self.information['VELO0_Y']) , 'the orbital velocity of CSST in Y direction at EXPSTART')
ofd_b.header['VELO0_Z']=( np.float64(self.information['VELO0_Z']) , 'the orbital velocity of CSST in Z direction at EXPSTART')
ofd_b.header['Euler0_1']=( np.float64(0.0), 'Euler angle 1 at EXPSTART')
ofd_b.header['Euler0_2']=( np.float64(0.0), 'Euler angle 2 at EXPSTART')
ofd_b.header['Euler0_3']=( np.float64(0.0), 'Euler angle 3 at EXPSTART')
ofd_b.header['EXPEND'] =(np.float64(time2jd(t2)) , 'exposure end time')
ofd_b.header['SUNANGL1']=(np.float32(0.0), 'angle between sun and optical axis at EXPEND')
ofd_b.header['MOONANG1']=(np.float32(0.0) , 'angle between moon and optical axis at EXPEND ')
ofd_b.header['POS_ANG1']=(np.float64(0.0) , 'angle between optical axis and the North Pole at EXPEND in arcsec')
ofd_b.header['TEL_ALT1']=(np.float64(0.0) , 'angle between optical axis and the ground- piston at EXPEND in deg ')
ofd_b.header['HOODSTA1']=(np.float32(0.0), 'lens hood altitude at EXPEND ')
ofd_b.header['HOODANG1']=(np.float32(0.0), 'lens hood azimuth at EXPEND ')
ofd_b.header['POSI1_X'] =(np.float64(self.information['POSI1_X']) , 'the orbital position of CSST in X direction at EXPEND')
ofd_b.header['POSI1_Y'] =(np.float64(self.information['POSI1_Y']) , 'the orbital position of CSST in Y direction at EXPEND')
ofd_b.header['POSI1_Z'] =(np.float64(self.information['POSI1_Z']) , 'the orbital position of CSST in Z direction at EXPEND')
ofd_b.header['VELO1_X']=(np.float64(self.information['VELO1_X']) , 'the orbital velocity of CSST in X direction at EXPEND')
ofd_b.header['VELO1_Y']=(np.float64(self.information['VELO1_Y']) , 'the orbital velocity of CSST in Y direction at EXPEND')
ofd_b.header['VELO1_Z']=(np.float64(self.information['VELO1_Z']) , 'the orbital velocity of CSST in Z direction at EXPEND')
ofd_b.header['Euler1_1']=( np.float64(0.0), 'Euler angle 1 at EXPEND')
ofd_b.header['Euler1_2']=( np.float64(0.0), 'Euler angle 2 at EXPEND')
ofd_b.header['Euler1_3']=( np.float64(0.0), 'Euler angle 3 at EXPEND')
ofd_b.header['RA_PNT1']=(np.float64(ofd_b.header['RA_PNT0']), 'RA of the pointing (degrees) at EXPEND in deg')
ofd_b.header['DEC_PNT1']=(np.float64(ofd_b.header['DEC_PNT0']), 'DEC of the pointing (degrees) at EXPEND in deg')
ofd_b.header['EXPTIME']=(self.information['exptime'], 'exposure duration')
ofd_b.header['EPOCH'] =(np.float32(0.0), 'coordinate epoch')
ofd_b.header['CHECKSUM']=( 0 , 'hdu-checksum')
########## finish header for 0 layer
#############################################################################3
##### header
b1= self.image_b[1856:3712,0:3392] #b4
b2= self.image_b[1856:3712,3392:6784] #b3
b3= self.image_b[0:1856,0:3392] #b1
b4= self.image_b[0:1856,3392:6784] #b2
####### do Flip the b2 b2 and b4 array in the up/down or left/right direction.
b2=np.fliplr(b2) ## left to right
b3=np.flipud(b3) ## down to up
b4=np.fliplr(b4) ## left to right and down to up
b4=np.flipud(b4)
bb=np.hstack((b1,b2,b3,b4))
#new image HDU, blue channel, layer 1
hdu_b =fits.ImageHDU(data=np.uint16(bb))
######### instrument information ######
#####
hdu_b.header['PMIRRPOS']=(bool(False), 'FSM pointing,T: to MCI, F: not to MCI')
if self.source =='sci':
hdu_b.header['CMIRRPOS']=(bool(False), 'position of calibration switch mirror,T: for calibration, F: not')
else:
hdu_b.header['CMIRRPOS']=(bool(True), 'position of calibration switch mirror,T: for calibration, F: not')
if self.source=='flat':
hdu_b.header['FLAMP'] =(int(1), 'status of flat lamp,0: off, 1: , 2: ')
else:
hdu_b.header['FLAMP'] =(int(0), 'status of flat lamp,0: off, 1: , 2: ')
if self.source=='lamp':
hdu_b.header['ALAMP'] =(int(1),'status of atomic emission line lamp,0: off, 1: , 2: ')
else:
hdu_b.header['ALAMP'] =(int(0),'status of atomic emission line lamp,0: off, 1: , 2: ')
#############
hdu_b.header['IFSMODE'] =(int(0), 'IFS working mode')
hdu_b.header['IFSTEMP'] =(float(0.0), 'IFS components temperature in degC')
hdu_b.header['IFSSTAT'] =(int(0), 'IFS components status parameter')
##############################################################################
################### detector information#############################
# hdu_b.header['COMMENT'] ='=========================================================='
# hdu_b.header['COMMENT'] ='Detector information'
# hdu_b.header['COMMENT'] ='=========================================================='
hdu_b.header['CAMERA'] =('Blue','camera of IFS')
hdu_b.header['DETNAM'] =('CCD231-c4','detector name')
hdu_b.header['DETSIZE'] =('', 'detector size')
hdu_b.header['DATASEC'] =('', 'data section')
hdu_b.header['PIXSCAL1']=(1856, 'pixel scale for axis 1')
hdu_b.header['PIXSCAL2']=(3392*4, 'pixel scale for axis 2')
hdu_b.header['PIXSIZE1']=(15, 'pixel size in um')
hdu_b.header['PIXSIZE2']=(15, 'pixel size in um')
hdu_b.header['NCHAN'] =(4, 'number of readout channels')
hdu_b.header['NCHAN1'] =(2, 'number of horizontal channels')
hdu_b.header['NCHAN2'] =(2, 'number of verticalchannels')
hdu_b.header['PSCAN1'] =(0, 'horizontal prescan width, per readout channel')
hdu_b.header['PSCAN2'] =(0, 'vertical prescan width, per readout channel')
hdu_b.header['OSCAN1'] =(0,' horizontal overscan width, per readout channel')
hdu_b.header['OSCAN2'] =(0, 'vertical overscan width, per readout channel')
## Readout information
# hdu_b.header['COMMENT'] ='============================================================='
# hdu_b.header['COMMENT'] ='Readout information'
# hdu_b.header['COMMENT'] ='============================================================='
hdu_b.header['READT0'] =(np.float64(time2jd(t2)),'read start time (UTC)')
hdu_b.header['READT1'] =(np.float64(time2jd(t3)), 'read end time (UTC)')
hdu_b.header['DETTEMP0']=(np.float32(0.0), 'detector temperature at READT0')
hdu_b.header['DETTEMP1']=(np.float32(0.0), 'detector temperature at READT1')
hdu_b.header['BIN_X'] =(0, 'bin number in X (wavelength)')
hdu_b.header['BIN_Y'] =(0, 'bin number in Y (spatial)')
hdu_b.header['GAIN1'] =(self.information['gain4_b'],'CCD gain (channel 1)')
hdu_b.header['GAIN2'] =(self.information['gain3_b'],'CCD gain (channel 2)')
hdu_b.header['GAIN3'] =(self.information['gain1_b'],'CCD gain (channel 3)')
hdu_b.header['GAIN4'] =(self.information['gain2_b'],'CCD gain (channel 4)')
hdu_b.header['DARK1'] =(self.information['dark4_b'],'CCD dark (channel 1)')
hdu_b.header['DARK2'] =(self.information['dark3_b'],'CCD dark (channel 2)')
hdu_b.header['DARK3'] =(self.information['dark1_b'],'CCD dark (channel 3)')
hdu_b.header['DARK4'] =(self.information['dark2_b'],'CCD dark (channel 4)')
hdu_b.header['RDNOIS1'] =(self.information['rn4_b'],'read noise (channel 1')
hdu_b.header['RDNOIS2'] =(self.information['rn3_b'],'read noise (channel 2')
hdu_b.header['RDNOIS3'] =(self.information['rn1_b'],'read noise (channel 3')
hdu_b.header['RDNOIS4'] =(self.information['rn2_b'],'read noise (channel 4')
hdu_b.header['DETBIA1'] =(self.information['bias4_b'],'amplifier bias voltage (channel1)')
hdu_b.header['DETBIA2'] =(self.information['bias3_b'],'amplifier bias voltage (channel2)')
hdu_b.header['DETBIA3'] =(self.information['bias1_b'],'amplifier bias voltage (channel3)')
hdu_b.header['DETBIA4'] =(self.information['bias2_b'],'amplifier bias voltage (channel4)')
hdu_b.header['RDSPEED'] =(100,'read speed (in MHz)')
hdu_b.header['EXPTIME'] =(self.information['exptime'],'exposure time in seconds')
hdu_b.header['Img_Ver'] =(self.information['img_ver'], 'IFS CCD image Version')
hdu_b.header['sky_obj'] =(self.skyfilepath,'input sky fits filepath')
##########################################################
#################### red camera ######################
#create a new FITS file, using HDUList instance
ofd_r = fits.PrimaryHDU()
ofd_r.header['GROUPS']=( bool(False), 'always F')
ofd_r.header['DATE'] =( data_time, 'date this file was written' )
ofd_r.header['FILENAME']=(filename_r, ' file name C48 ')
ofd_r.header['OBSTYPE'] =( self.source, 'observation type raw,flt, mask, bias, dark, sci')
ofd_r.header['TELESCOP']=('CSST' , 'always CSST')
ofd_r.header['INSTRUME']=( 'IFS', ' ')
ofd_r.header['RADECSYS']=('ICRS', ' always ICRS ')
ofd_r.header['EQUINOX'] =( float(2000.0), 'always 2000.0')
ofd_r.header['FITSCREA']=( '4.2.1' , 'FITS create software version')
######### Object information #############
######### Object information #############
# ofd_r.header['COMMENT']='======================================================================='
# ofd_r.header['COMMENT']='Object information'
# ofd_r.header['COMMENT']='======================================================================='
ofd_r.header['OBJECT']=( self.information['name_obj'], 'object name')
ofd_r.header['TARGET']=( (self.information['target']), 'target name, hhmmss+ddmmss')
ofd_r.header['OBJ_RA'] =(np.float64(self.information['ra_obj']) , 'RA of the object in deg')
ofd_r.header['OBJ_DEC']=(np.float64(self.information['dec_obj']) , 'DEC of the object in deg')
ofd_r.header['RA_PNT0']=( np.float64(self.information['ra_pnt0']) , 'RA of the pointing (degrees) at EXPSTART')
ofd_r.header['DEC_PNT0']=(np.float64(self.information['dec_pnt0']) , 'DEC of the pointing (degrees) at EXPSTART')
ofd_r.header['OBSID'] =(str(obsid) , 'observation ID, 3+8bit ')
######## Telescope information ###############
# ofd_r.header['COMMENT']='======================================================================='
# ofd_r.header['COMMENT']='Telescope information'
# ofd_r.header['COMMENT']='======================================================================='
ofd_r.header['REFFRAME']=('CSSTGSC-1.0' , 'guide star catalog version ')
ofd_r.header['DATE-OBS']=(data_time , 'date of the observation (yyyy-mm-dd hh:mm:ss)')
ofd_r.header['EXPSTART']=(np.float64(exp_starttime), 'exposure start time ')
ofd_r.header['SUNANGL0']=(np.float32(0.0) , 'angle between sun and optical axis at EXPSTART ')
ofd_r.header['MOONANG0']=(np.float32(0.0) , 'angle between moon and optical axis at EXPSTART ')
ofd_r.header['POS_ANG0']=(np.float64(0.0), 'angle between optical axis and the North Pole at EXPSTART in arcsec ')
ofd_r.header['TEL_ALT0']=(np.float64(0.0), 'angle between optical axis and the ground- piston at EXPSTART in deg ')
ofd_r.header['HOODSTA0']=(np.float32(0.0) , 'lens hood altitude at EXPSTART ')
ofd_r.header['HOODANG0']=(np.float32(0.0), 'lens hood azimuth at EXPSTART ')
ofd_r.header['POSI0_X'] =(np.float64(self.information['POSI0_X']) , 'the orbital position of CSST in X direction at EXPSTART')
ofd_r.header['POSI0_Y'] =(np.float64(self.information['POSI0_Y']) , 'the orbital position of CSST in Y direction at EXPSTART')
ofd_r.header['POSI0_Z'] =(np.float64(self.information['POSI0_Z']) , 'the orbital position of CSST in Z direction at EXPSTART')
ofd_r.header['VELO0_X']=( np.float64(self.information['VELO0_X']) , 'the orbital velocity of CSST in X direction at EXPSTART')
ofd_r.header['VELO0_Y']=( np.float64(self.information['VELO0_Y']) , 'the orbital velocity of CSST in Y direction at EXPSTART')
ofd_r.header['VELO0_Z']=( np.float64(self.information['VELO0_Z']) , 'the orbital velocity of CSST in Z direction at EXPSTART')
ofd_r.header['Euler0_1']=( np.float64(0.0), 'Euler angle 1 at EXPSTART')
ofd_r.header['Euler0_2']=( np.float64(0.0), 'Euler angle 2 at EXPSTART')
ofd_r.header['Euler0_3']=( np.float64(0.0), 'Euler angle 3 at EXPSTART')
ofd_r.header['EXPEND'] =(np.float64(exp_endtime) , 'exposure end time')
ofd_r.header['SUNANGL1']=(np.float32(0.0), 'angle between sun and optical axis at EXPEND')
ofd_r.header['MOONANG1']=(np.float32(0.0) , 'angle between moon and optical axis at EXPEND ')
ofd_r.header['POS_ANG1']=(np.float64(0.0) , 'angle between optical axis and the North Pole at EXPEND in arcsec')
ofd_r.header['TEL_ALT1']=(np.float64(0.0) , 'angle between optical axis and the ground- piston at EXPEND in deg ')
ofd_r.header['HOODSTA1']=(np.float32(0.0), 'lens hood altitude at EXPEND ')
ofd_r.header['HOODANG1']=(np.float32(0.0), 'lens hood azimuth at EXPEND ')
ofd_r.header['POSI1_X'] =(np.float64(self.information['POSI1_X']) , 'the orbital position of CSST in X direction at EXPEND')
ofd_r.header['POSI1_Y'] =(np.float64(self.information['POSI1_Y']) , 'the orbital position of CSST in Y direction at EXPEND')
ofd_r.header['POSI1_Z'] =(np.float64(self.information['POSI1_Z']) , 'the orbital position of CSST in Z direction at EXPEND')
ofd_r.header['VELO1_X']=(np.float64(self.information['VELO1_X']) , 'the orbital velocity of CSST in X direction at EXPEND')
ofd_r.header['VELO1_Y']=(np.float64(self.information['VELO1_Y']) , 'the orbital velocity of CSST in Y direction at EXPEND')
ofd_r.header['VELO1_Z']=(np.float64(self.information['VELO1_Z']) , 'the orbital velocity of CSST in Z direction at EXPEND')
ofd_r.header['Euler1_1']=( np.float64(0.0), 'Euler angle 1 at EXPEND')
ofd_r.header['Euler1_2']=( np.float64(0.0), 'Euler angle 2 at EXPEND')
ofd_r.header['Euler1_3']=( np.float64(0.0), 'Euler angle 3 at EXPEND')
ofd_r.header['RA_PNT1']=(np.float64(ofd_r.header['RA_PNT0']), 'RA of the pointing (degrees) at EXPEND in deg')
ofd_r.header['DEC_PNT1']=(np.float64(ofd_r.header['DEC_PNT0']), 'DEC of the pointing (degrees) at EXPEND in deg')
ofd_r.header['EXPTIME']=(self.information['exptime'], 'exposure duration')
ofd_r.header['EPOCH'] =(np.float32(0.0), 'coordinate epoch')
ofd_r.header['CHECKSUM']=( 0 , 'hdu-checksum')
### finish 0 layer header
########## finish header for 0 layer
# ########## blue zone 1--to--3
# self.image_r[0:1856,0:3392] += self.information['bias1_r']
# ########## zone 4 --to---1 #################
# self.image_r[1856:3712,0:3392] += self.information['bias4_r']
# ########## zone 2 ----to----4 ###################
# self.image_r[0:1856,3392:6784] += self.information['bias2_r']
# ########## zone 3 ---to------2
# self.image_r[1856:3712,3392:6784] += self.information['bias3_r']
#############################################################################3
##### header
b1= self.image_r[1856:3712,0:3392]
b2= self.image_r[1856:3712,3392:6784]
b3= self.image_r[0:1856,0:3392]
b4= self.image_r[0:1856,3392:6784]
####### do Flip the b2 b2 and b4 array in the up/down or left/right direction.
b2=np.fliplr(b2) ## left to right
b3=np.flipud(b3) ## down to up
b4=np.fliplr(b4) ## left to right and down to up
b4=np.flipud(b4)
rr=np.hstack((b1,b2,b3,b4))
#new image HDU, blue channel, layer 1
hdu_r =fits.ImageHDU(data=np.uint16(rr))
#########################################
######### instrument information ######
hdu_r.header['PMIRRPOS']=(bool(False), 'FSM pointing,T: to MCI, F: not to MCI')
hdu_r.header['CMIRRPOS']=(bool(False), 'position of calibration switch mirror,T: for calibration, F: not')
hdu_r.header['FLAMP'] =(int(0), 'status of flat lamp,0: off, 1: , 2: ')
hdu_r.header['ALAMP'] =(int(0),'status of atomic emission line lamp,0: off, 1: , 2: ')
hdu_r.header['IFSMODE'] =(int(0), 'IFS working mode')
hdu_r.header['IFSTEMP'] =(float(0.0), 'IFS components temperature in degC')
hdu_r.header['IFSSTAT'] =(int(0), 'IFS components status parameter')
################### detector information#############################
# hdu_r.header['COMMENT']='======================================================================='
# hdu_r.header['COMMENT']='Detector information'
# hdu_r.header['COMMENT']='======================================================================='
hdu_r.header['CAMERA'] =('Red','camera of IFS')
hdu_r.header['DETNAM'] =('CCD231-c4','detector name')
hdu_r.header['DETSIZE'] =('', 'detector size')
hdu_r.header['DATASEC'] =('', 'data section')
hdu_r.header['PIXSCAL1']=(1856, 'pixel scale for axis 1')
hdu_r.header['PIXSCAL2']=(3392*4, 'pixel scale for axis 2')
hdu_r.header['PIXSIZE1']=(15, 'pixel size in um')
hdu_r.header['PIXSIZE2']=(15, 'pixel size in um')
hdu_r.header['NCHAN'] =(4, 'number of readout channels')
hdu_r.header['NCHAN1'] =(2, 'number of horizontal channels')
hdu_r.header['NCHAN2'] =(2, 'number of verticalchannels')
hdu_r.header['PSCAN1'] =(0, 'horizontal prescan width, per readout channel')
hdu_r.header['PSCAN2'] =(0, 'vertical prescan width, per readout channel')
hdu_r.header['OSCAN1'] =(0,' horizontal overscan width, per readout channel')
hdu_r.header['OSCAN2'] =(0, 'vertical overscan width, per readout channel')
#####################################################################################################
## Readout information
# hdu_r.header['COMMENT']='======================================================================='
# hdu_r.header['COMMENT']='Readout information'
# hdu_r.header['COMMENT']='======================================================================='
hdu_r.header['READT0'] =(np.float64(time2jd(t2)),'read start time (UTC)')
hdu_r.header['READT1'] =(np.float64(time2jd(t3)), 'read end time (UTC)')
hdu_r.header['DETTEMP0']=(np.float32(0.0), 'detector temperature at READT0')
hdu_r.header['DETTEMP1']=(np.float32(0.0), 'detector temperature at READT1')
hdu_r.header['BIN_X'] =(0, 'bin number in X (wavelength)')
hdu_r.header['BIN_Y'] =(0, 'bin number in Y (spatial)')
hdu_r.header['GAIN1'] =(self.information['gain4_r'],'CCD gain (channel 1)')
hdu_r.header['GAIN2'] =(self.information['gain3_r'],'CCD gain (channel 2)')
hdu_r.header['GAIN3'] =(self.information['gain1_r'],'CCD gain (channel 3)')
hdu_r.header['GAIN4'] =(self.information['gain2_r'],'CCD gain (channel 4)')
hdu_r.header['DARK1'] =(self.information['dark4_r'],'CCD dark (channel 1)')
hdu_r.header['DARK2'] =(self.information['dark3_r'],'CCD dark (channel 2)')
hdu_r.header['DARK3'] =(self.information['dark1_r'],'CCD dark (channel 3)')
hdu_r.header['DARK4'] =(self.information['dark2_r'],'CCD dark (channel 4)')
hdu_r.header['RDNOIS1'] =(self.information['rn4_r'],'read noise (channel 1')
hdu_r.header['RDNOIS2'] =(self.information['rn3_r'],'read noise (channel 2')
hdu_r.header['RDNOIS3'] =(self.information['rn1_r'],'read noise (channel 3')
hdu_r.header['RDNOIS4'] =(self.information['rn2_r'],'read noise (channel 4')
hdu_r.header['DETBIA1'] =(self.information['bias4_r'],'amplifier bias voltage (channel1)')
hdu_r.header['DETBIA2'] =(self.information['bias3_r'],'amplifier bias voltage (channel2)')
hdu_r.header['DETBIA3'] =(self.information['bias1_r'],'amplifier bias voltage (channel3)')
hdu_r.header['DETBIA4'] =(self.information['bias2_r'],'amplifier bias voltage (channel4)')
hdu_r.header['RDSPEED'] =(100,'read speed (in MHz)')
hdu_r.header['EXPTIME'] =(self.information['exptime'],'exposure time in seconds')
hdu_r.header['Img_Ver'] =(self.information['img_ver'], 'IFS CCD image Version')
hdu_r.header['sky_obj'] =(self.skyfilepath,'input sky fits filename')
hdulist_b=fits.HDUList([ofd_b, hdu_b])
hdulist_b.writeto(file_b, overwrite=True)
#print('IFS_b.fits is created ')
hdulist_r=fits.HDUList([ofd_r, hdu_r])
hdulist_r.writeto(file_r, overwrite=True)
#print('IFS_r.fits is created ')
##################################################################################
def earthshine(self,theta):
"""
For given theta angle, return the earth-shine spectrum.
:param theta: angle (in degree) from the target to earth limb.
:return: the scaled solar spectrum
template_wave: unit in A
template_flux: unit in erg/s/cm^2/A/arcsec^2
"""
# read solar template
solar_template = pd.read_csv(self.information['indata_path']+'/refs/solar_spec.dat', sep='\s+',
header=None, comment='#')
template_wave = solar_template[0].values
template_flux = solar_template[1].values
# read earth shine surface brightness
earthshine_curve = pd.read_csv(self.information['indata_path']+'/refs/earthshine.dat',
header=None, comment='#')
angle = earthshine_curve[0].values
surface_brightness = earthshine_curve[1].values
# read V-band throughtput
cat_filter_V = pd.read_csv(self.information['indata_path']+'/refs/filter_Bessell_V.dat', sep='\s+',
header=None, comment='#')
filter_wave = cat_filter_V[0].values
filter_response = cat_filter_V[1].values
# interplate to the target wavelength in V-band
ind_filter = (template_wave >= np.min(filter_wave)) & (template_wave <= np.max(filter_wave))
filter_wave_interp = template_wave[ind_filter]
filter_response_interp = np.interp(filter_wave_interp, filter_wave, filter_response)
filter_constant = simps(filter_response_interp * filter_wave_interp, filter_wave_interp)
template_constant = simps(filter_response_interp * template_wave[ind_filter] * template_flux[ind_filter],
template_wave[ind_filter])
dwave = filter_wave_interp[1:] - filter_wave_interp[:-1]
wave_eff = np.nansum(dwave * filter_wave_interp[1:] * filter_response_interp[1:]) / \
np.nansum(dwave * filter_response_interp[1:])
# get the normalized value at theta.
u0 = np.interp(theta, angle, surface_brightness) # mag/arcsec^2
u0 = 10**((u0 + 48.6)/(-2.5)) # target flux in erg/s/cm^2/Hz unit
u0 = u0 * 3e18 / wave_eff**2 # erg/s/cm^2/A/arcsec^2
factor = u0 * filter_constant / template_constant
norm_flux = template_flux * factor # erg/s/cm^2/A/arcsec^2
self.earthshine_wave=template_wave # A
self.earthshine_flux=norm_flux
return
########################################################################################################################################################################################################################################################
def zodiacal(self, ra, dec, time):
"""
For given RA, DEC and TIME, return the interpolated zodical spectrum in Leinert-1998.
:param ra: RA in unit of degree, ICRS frame
:param dec: DEC in unit of degree, ICRS frame
:param time: the specified string that in ISO format i.e., yyyy-mm-dd.
:return:
wave_A: wavelength of the zodical spectrum
spec_mjy: flux of the zodical spectrum, in unit of MJy/sr
spec_erg: flux of the zodical spectrum, in unit of erg/s/cm^2/A/sr
"""
# get solar position
dt = datetime.fromisoformat(time)
jd = julian.to_jd(dt, fmt='jd')
t = Time(jd, format='jd', scale='utc')
astro_sun = get_sun(t)
ra_sun, dec_sun = astro_sun.gcrs.ra.deg, astro_sun.gcrs.dec.deg
radec_sun = SkyCoord(ra=ra_sun*u.degree, dec=dec_sun*u.degree, frame='gcrs')
lb_sun = radec_sun.transform_to('geocentrictrueecliptic')
# get offsets between the target and sun.
radec_obj = SkyCoord(ra=ra*u.degree, dec=dec*u.degree, frame='icrs')
lb_obj = radec_obj.transform_to('geocentrictrueecliptic')
beta = abs(lb_obj.lat.degree)
lamda = abs(lb_obj.lon.degree - lb_sun.lon.degree)
# interpolated zodical surface brightness at 0.5 um
zodi = pd.read_csv(self.information['indata_path']+'/refs/zodi_map.dat', sep='\s+', header=None, comment='#')
beta_angle = np.array([0, 5, 10, 15, 20, 25, 30, 45, 60, 75])
lamda_angle = np.array([0, 5, 10, 15, 20, 25, 30, 35, 40, 45,
60, 75, 90, 105, 120, 135, 150, 165, 180])
xx, yy = np.meshgrid(beta_angle, lamda_angle)
f = interpolate.interp2d(xx, yy, zodi, kind='linear')
zodi_obj = f(beta, lamda) # 10^−8 W m−2 sr−1 um−1
# read the zodical spectrum in the ecliptic
cat_spec = pd.read_csv(self.information['indata_path']+'/refs/solar_spec.dat', sep='\s+', header=None, comment='#')
wave = cat_spec[0].values # A
spec0 = cat_spec[1].values # 10^-8 W m^−2 sr^−1 μm^−1
zodi_norm = 252 # 10^-8 W m^−2 sr^−1 μm^−1
spec = spec0 * (zodi_obj / zodi_norm) * 1e-8 # W m^−2 sr^−1 μm^−1
# convert to the commonly used unit of MJy/sr, erg/s/cm^2/A/sr
wave_A = wave # A
spec_mjy = spec * 0.1 * wave_A**2 / 3e18 * 1e23 * 1e-6 # MJy/sr
spec_erg = spec * 0.1 # erg/s/cm^2/A/sr
spec_erg2 = spec_erg / 4.25452e10 # erg/s/cm^2/A/arcsec^2
self.zodiacal_wave=wave_A # in A
self.zodiacal_flux=spec_erg2
return
###################################################################################
##################################################################################
def CalskyNoise(self, lam):
# calculate sky noise;
planckh= 6.62620*10**-27 # erg s;
cc=2.99792458*10**17 # in nm/s
fov2=0.01 # arcsec^2
# lam is input wavelength in nm
##########################################
self.earthshine_wave # A
self.earthshine_flux # erg/s/cm^2/A/arcsec^2
earthshine_flux=np.interp(lam*10.0, self.earthshine_wave,self.earthshine_flux) # flux from zodiacal
self.zodiacal_wave # in A
self.zodiacal_flux # erg/s/cm^2/A/arcsec^2
zodiacal_flux=np.interp(lam*10.0, self.zodiacal_wave,self.zodiacal_flux) # flux from zodiacal
fluxlam_sky=(earthshine_flux+zodiacal_flux)*fov2 # erg/s/cm2/A
###############
ephoton=planckh*cc/lam # in erg/photon, cc与lambda单位需要一致；
Ns_skynoise=fluxlam_sky/ephoton # in unit of photons/cm2/s/A
return Ns_skynoise
#################################################################################
def sim_sky_img(self,skyfitsfilename, skyRa_shift, skyDec_shift, sky_rot, exposuretime):
############################################################################
### load fits file
indatafile=skyfitsfilename
a=fits.open(indatafile)
#######################################
self.information['name_obj']=a[0].header['OBJECT']
self.information['ra_obj'] =a[0].header['RA'] ### in degree
self.information['dec_obj'] =a[0].header['DEC'] ### in degree
disRa =(skyRa_shift)/3600.0 ##convert unit of degree to arcsec
disDec=(skyDec_shift)/3600.0 ##convert unit of degree to arcsec
self.information['ra_pnt0'] =a[0].header['RA'] + disRa/np.cos(a[0].header['DEC']/180.0*np.pi)
self.information['dec_pnt0']=a[0].header['DEC']+ disDec
self.earthshine(self.earthshine_theta)
self.zodiacal(self.information['ra_obj'], self.information['dec_obj'], self.zodiacal_time)
self.information['target']=deg2HMS(self.information['ra_obj'], self.information['dec_obj'])
### main input data
SpecCube=a[1].data ## spectrum data cube;
Wave=0.1*a[2].data # the relatived wavelength which is converted from Unit A to nm
#print('Wave data header', hdr)
######################################################################################
exptime=self.information['exptime'] #exposure time
dis_dx=disRa # image shift Ra in arcsec
dis_dy=disDec # image shift Dec in arcsec
sizeout= len(SpecCube[:,0,0])
blue_img=galsim.Image(np.zeros((2048,4096)),copy=True)
blue_img.scale=self.pixelscale
blue_img.setOrigin=(0,0)
red_img=galsim.Image(np.zeros((3072,6144)),copy=True)
red_img.scale=self.pixelscale
red_img.setOrigin=(0,0)
blue_sensor=galsim.Sensor()
red_sensor=galsim.Sensor()
deltalam=np.mean(np.diff(Wave))
energy=0.0
energy_blue=0.0
energy_red=0.0
width_blue=0
################################
############## doppler effect to photons.wavelength #############
#self.orbit_pars
x_sat=float(self.orbit_pars[self.orbit_exp_num,1])
y_sat=float(self.orbit_pars[self.orbit_exp_num,2])
z_sat=float(self.orbit_pars[self.orbit_exp_num,3])
vx_sat=float(self.orbit_pars[self.orbit_exp_num,4])
vy_sat=float(self.orbit_pars[self.orbit_exp_num,5])
vz_sat=float(self.orbit_pars[self.orbit_exp_num,6])
self.information['POSI0_X']=x_sat
self.information['POSI0_Y']=y_sat
self.information['POSI0_Z']=z_sat
self.information['VELO0_X']=vx_sat
self.information['VELO0_Y']=vy_sat
self.information['VELO0_Z']=vz_sat
theta1=beta_angle( x_sat,y_sat,z_sat,vx_sat,vy_sat,vz_sat, self.information['ra_obj'], self.information['dec_obj'])
v1=np.sqrt(vx_sat**2+vy_sat**2+vz_sat**2)*np.cos(theta1/180.0*np.pi) # velocity at stat exposure time
vv1=LSR_velocity(self.information['ra_obj'],self.information['dec_obj'],v1, self.TianCe_day)
#################################################
### exposure end time is t2 ;
t2=self.dt+timedelta(seconds=self.information['exptime'])
t2jd=time2jd(t2)
if self.orbit_pars[-1,0]<t2jd: ## orbit parameters are not in currenct txt file
self.orbit_file_num=self.orbit_file_num+1
fn=self.information['indata_path']+'/refs/orbit20160925/'+str(self.orbit_file_num)+'.txt'
self.orbit_pars=np.loadtxt(fn)
self.orbit_exp_num=0
for k in range(self.orbit_exp_num, len(self.orbit_pars),1):
if t2jd-self.orbit_pars[k,0]<0:
break
if k==0:
deltaT=jd2time(self.orbit_pars[k,0])-t2
p1x=self.orbit_pars[k,1]-(self.orbit_pars[k+1,1]-self.orbit_pars[k,1])*deltaT.seconds/120
p1y=self.orbit_pars[k,2]-(self.orbit_pars[k+1,2]-self.orbit_pars[k,2])*deltaT.seconds/120
p1z=self.orbit_pars[k,3]-(self.orbit_pars[k+1,3]-self.orbit_pars[k,3])*deltaT.seconds/120
p1vx=self.orbit_pars[k,4]-(self.orbit_pars[k+1,4]-self.orbit_pars[k,4])*deltaT.seconds/120
p1vx=self.orbit_pars[k,5]-(self.orbit_pars[k+1,5]-self.orbit_pars[k,5])*deltaT.seconds/120
p1vx=self.orbit_pars[k,6]-(self.orbit_pars[k+1,6]-self.orbit_pars[k,6])*deltaT.seconds/120
else:
deltaT=jd2time(self.orbit_pars[k,0])-t2
p1x=self.orbit_pars[k-1,1]+(self.orbit_pars[k,1]-self.orbit_pars[k-1,1])*deltaT.seconds/120
p1y=self.orbit_pars[k-1,2]+(self.orbit_pars[k,2]-self.orbit_pars[k-1,2])*deltaT.seconds/120
p1z=self.orbit_pars[k-1,3]+(self.orbit_pars[k,3]-self.orbit_pars[k-1,3])*deltaT.seconds/120
p1vx=self.orbit_pars[k-1,4]+(self.orbit_pars[k,4]-self.orbit_pars[k-1,4])*deltaT.seconds/120
p1vy=self.orbit_pars[k-1,5]+(self.orbit_pars[k,5]-self.orbit_pars[k-1,5])*deltaT.seconds/120
p1vz=self.orbit_pars[k-1,6]+(self.orbit_pars[k,6]-self.orbit_pars[k-1,6])*deltaT.seconds/120
#######
self.information['POSI1_X']=p1x
self.information['POSI1_Y']=p1y
self.information['POSI1_Z']=p1z
self.information['VELO1_X']=p1vx
self.information['VELO1_Y']=p1vy
self.information['VELO1_Z']=p1vz
theta2=beta_angle(p1x, p1y, p1z, p1vx, p1vy, p1vz,self.information['ra_obj'], self.information['dec_obj'])
v2=np.sqrt(p1vx**2+p1vy**2+p1vz**2)*np.cos(theta2/180.0*np.pi) # velocity at end exposure time
vv2=LSR_velocity(self.information['ra_obj'],self.information['dec_obj'], v2, self.TianCe_day)
#### get slice and slit mask
########################################################################################
for ilam in range(len(Wave)):
#print('ilam=', ilam)
if ilam%500==0:
self.log.info('ilam = %i' % ilam)
lam=Wave[ilam] # the wavelength of the i-th frame data
if lam<350:
continue
###############################################
Specimg=SpecCube[:,:,ilam] # get the i-th frame of the input SpecCube data;
Nspecimg=SpecCube2photon(Specimg,lam) # convert to photons/cm2/s/A
Nskynoise=self.CalskyNoise(lam)
### add sky noise
if self.skyback:
Nimg=Nspecimg+Nskynoise
else:
Nimg=Nspecimg
# multipe the tel area, exposure time, and bandwidh;
Nimg=Nimg*self.telarea*exptime*deltalam*10.0 # photons/cm2/s/A to photons, here Nimg size is 100*100
##########################################################################
### shift image with photons position and rotate them round the image true center
img=galsim.Image(Nimg,copy=True)
img.scale = self.pixelscale
img.setOrigin(0,0)
photons=galsim.PhotonArray.makeFromImage(img, max_flux=max(Nimg.max()/1000.0, 1.0))
### now shift and rotated photons
self.information['shift_dx']=dis_dx
self.information['shift_dy']=dis_dy
photons.x=photons.x-dis_dx ### apply shift to photons position in ra direction
photons.y=photons.y-dis_dy ### apply shift to photons position in dec direction
x2,y2=rotation_yan(img.true_center.x, img.true_center.y, photons.x/img.scale, photons.y/img.scale, sky_rot)
#print('rotation time=', t2-t1)
photons2=galsim.PhotonArray(N=len(x2), x=x2*img.scale, y=y2*img.scale, flux = photons.flux)
photons=photons2
## convert the photon image to galsim photons;
###################################################################
image0 = galsim.Image(sizeout,sizeout)
image0.scale = self.pixelscale
image0.setOrigin(0,0)
rotphotons=galsim.PhotonArray(N=len(photons.x), x=x2, y=y2, flux = photons2.flux)
sensor = galsim.Sensor()
sensor.accumulate(rotphotons, image0)
#####################################
#####################################################################
### do convolve image0 with PSF0 from primay CSST ###
### calculate the PSF0 at this wavelength
Q=lam*1e-3*self.Fnum/self.pixelsize
wavefront0=self.opd0/lam
wavefront=ndimage.rotate(wavefront0, -sky_rot, order=1, reshape=False) # here we choose reshape=False, rotate the wavefront
psf0=anySampledPSF(wavefront, self.pupil , Q , 64)
conv = fftconvolve(image0.array, psf0, mode='same')
conv[conv < 0.0] = 0.0
img0=conv
################ opd1 rms =0.075 @632.8nm
wavefront=self.opd1/lam
wavefront=ndimage.rotate(wavefront, -sky_rot, order=1, reshape=False) # here we choose reshape=False, rotate the wavefront
psf1=anySampledPSF(wavefront, self.pupil , Q , 32)
##### do convolve with psf1
conv = fftconvolve(img0, psf1, mode='same')
conv[conv < 0.0] = 0.0
img0=conv
#######################################################
energy=energy+img0[50-32:50+32,50-32:50+32].sum() ## calculate the slice image energy;
CCD_Qe_lam=np.interp(lam,self.CCD_Qe[:,0], self.CCD_Qe[:,1]) ## CCD quantum efficiency
Qe_blue = np.interp(lam, self.optical_blue_Q[:,0], self.optical_blue_Q[:,1]) # optical efficiency , convert the wavelength to A
Qe_red = np.interp(lam, self.optical_red_Q[:,0], self.optical_red_Q[:,1]) # optical efficiency, convert the wavelength to A
Qe_blue=Qe_blue*CCD_Qe_lam
Qe_red =Qe_red *CCD_Qe_lam
lam1=lam*(1+vv1/(3.0*1e5))
lam2=lam*(1+vv2/(3.0*1e5))
######## consider the slice optical efficiency in different slicer channel
da=fits.open(self.information['indata_path']+'/slicer_QE.fits')
slicer_Qe=da[0].data
img0=img0*slicer_Qe
########## do the slice effect ###################
for k in range(32):
#### do slice effect to get slice image
img1=img0*self.maskSlice[str(k)]
############ get opd2 and PSF2 ######################
wavefront2=self.opd2/lam
psf2=anySampledPSF(wavefront2, self.pupil , Q , 64)
psf2=ndimage.rotate(psf2, -sky_rot, order=1, reshape=False) # here we choose reshape=False, the rotated image will
##### do convolve
##################
conv = fftconvolve(img1, psf2, mode='same')
#suppress negative numbers
conv[conv < 0.0] = 0.0
img2=conv
##############do Slit Mask ###########################
img2=img2*self.maskSlit[str(k)]
######### get opd3 and PSF3 ##########################
wavefront3=self.opd3/lam
psf3=anySampledPSF(wavefront3, self.pupil , Q , 64)
psf3=ndimage.rotate(psf3, -sky_rot, order=1, reshape=False) # here we choose reshape=False, the rotated image will
##### do convolve
########################################
conv = fftconvolve(img2, psf3, mode='same')
#suppress negative numbers
conv[conv < 0.0] = 0.0
img3=conv
######################## get subimage #####################
subimage = galsim.Image(80,80)
subimage.array[:,:]=img3[int(sizeout/2)-40:int(sizeout/2)+40,int(sizeout/2)-40:int(sizeout/2)+40]
######################## get photons from sub-image #####################
subimage.scale = self.pixelscale
subimage.setOrigin(0,0)
photons=galsim.PhotonArray.makeFromImage(subimage, max_flux=max(img3.max()/1000.0,1.0))
#############################################################################
###### do something for each photons;
#######
idx0=np.where( photons.flux>1e-3)
energy=energy+sum(photons.flux[idx0]) ### totla energy for slice image
p_num=len(idx0[0])
###############################################################################################
############### find photons for blue channel, and make the flux multiple the optical and CCD efficiency
np.random.seed()
wavesample=lam1+(lam2-lam1)*np.random.rand(p_num)
if (lam>=350.0 and lam<=650.0):
## bulue channel
photons_blue=galsim.PhotonArray(N=p_num, x=photons.x[idx0], y=photons.y[idx0], flux = Qe_blue*photons.flux[idx0],wavelength=wavesample)
dx_blue, dy_blue = get_dx_dy_blue(wavesample)
photons_blue.x=photons_blue.x/self.pixelscale+dx_blue+self.slice_blue['px'][k]
photons_blue.y=photons_blue.y/self.pixelscale+dy_blue+self.slice_blue['py'][k]
blue_sensor.accumulate(photons_blue, blue_img)
energy_blue=energy_blue+sum(photons_blue.flux)
width_blue=width_blue+deltalam/32.0
if (lam>=560.0) & (lam<=1000.0 ):
## red channel
photons_red=galsim.PhotonArray(N=p_num,x=photons.x[idx0],y=photons.y[idx0], flux = Qe_red*photons.flux[idx0],wavelength=wavesample)
dx_red, dy_red = get_dx_dy_red(wavesample)
photons_red.x=photons_red.x/self.pixelscale+dx_red+self.slice_red['px'][k]
photons_red.y=photons_red.y/self.pixelscale+dy_red+self.slice_red['py'][k]
red_sensor.accumulate(photons_red, red_img)
energy_red=energy_red+sum(photons_red.flux)
####################################################################################
## stray light will cover 2% of input total light;
blue_img.array[:,:]=blue_img.array[:,:]+0.01*energy/2048/4096
red_img.array[:,:] =red_img.array[:,:]+ 0.01*energy/3072/6144
self.image_b=blue_img.array
self.image_r=red_img.array
return
#################################################################################################
#################################################################################################
def simulate(self, skyfitsin, skyRa_shift, skyDec_shift,sky_rot, exptime):
"""
Create a single simulated image of a quadrant defined by the configuration file.
"""
#self.configure() #print the configfile name and path;
self.dt=datetime.now()
self.information['exptime']=exptime
self.skyfilepath=skyfitsin
np.random.seed()
ud = np.random.random() # Choose a random
self.earthshine_theta=ud * 60 # in degree
##################################################################
#### load orbit parameters #####
flag=0
for k in range(1,50,1):
fn=self.information['indata_path']+'/refs/orbit20160925/'+str(k)+'.txt';
d=np.loadtxt(fn);
self.dt_num=int((self.information['exptime']+self.information['readouttime']+125)/120)
now_dt=datetime.utcnow()
now_jd=time2jd(now_dt)
for kk in range(len(d[:,0])):
if now_jd-d[kk,0]<=0:
flag=1
break
if flag==1:
break
#####################end for
self.orbit_pars=d
self.orbit_file_num=k
self.orbit_exp_num =kk
exptime_start_jd=d[kk,0] #### jd time, utc format
self.dt=julian.from_jd(exptime_start_jd, fmt='jd')
self.TianCe_day=self.dt.strftime("%Y-%m-%d") ###str(self.dt.year)+'-'+str(self.dt.month)+'-'+str(self.dt.day)
self.TianCe_exp_start=dt2hmd(self.dt)
self.zodiacal_time=self.TianCe_day
#######################################################################
self.sim_sky_img(skyfitsin, skyRa_shift, skyDec_shift, sky_rot, self.information['exptime'])
self.information['sky_rot']=sky_rot
self.information['skyRa_shift']=skyRa_shift
self.information['skyRa_shift']=skyDec_shift
###############################################################################
############ add some effect to images #######################################
if self.flatfieldM:
self.applyflatfield()
#print('Applying flatfieldM finished.......')
if self.darknoise:
self.applyDarkCurrent()
if self.cosmicRays:
self.addCosmicRays()
#print('Applying cosmicRays finished.......')
if self.bleeding:
self.applyBleeding_yan()
#print('Applying bleeding finished.......')
self.applyPoissonNoise()
if self.nonlinearity:
self.applyNonlinearity()
#print('Applying nonlinearity finished.......')
if self.radiationDamage:
self.applyRadiationDamage()
#print('Applying radiationDamage finished.......')
##### cut original CCD image to four parts by four read out channels and zones
self.CCDreadout()
if self.readoutNoise:
self.applyReadoutNoise()
#print('Applying readoutNoise finished.......')
self.electrons2ADU()
self.applyBias()
if self.cosmetics:
self.applyCosmetics()
#print('Applying cosmetics finished.......')
self.discretise()
self.writeOutputs()
self.log.info('Using the following input values:')
for key, value in self.information.items():
self.log.info('%s = %s' % (key, value))
self.log.info('Using the following booleans:')
for key, value in self.booleans.items():
self.log.info('%s = %s' % (key, value))
self.log.info('Finished the simulation.')
##############################################################################################
##############################################################################################
def runIFSsim(configfile):
opts, args = processArgs()
opts.configfile=configfile
simulate= IFSsimulator(opts)
skyfitsin=simulate.information['skyfitsin']
exptime=simulate.information['exptime']
sky_ra_dis=simulate.information['sky_ra_dis']
sky_dec_dis=simulate.information['sky_dec_dis']
sky_angle_dis=simulate.information['sky_angle_dis']
simulate.simulate(skyfitsin, sky_ra_dis, sky_dec_dis, sky_angle_dis, exptime)
########################## begin main fucntion #######################################
##############################################################################################
##############################################################################################
if __name__ == "__main__":
if len(sys.argv[:]) <2:
configfile='./ifs_data/ifs_sim_example.config'
###########################################################################################
if len(sys.argv[:]) >=2:
configfile=sys.argv[1]
if not os.path.exists(configfile):
print('The given input configfile path is wrong......')
sys.exit(1)
################################################
runIFSsim(configfile)
print('---The CSST-IFS simulation is successful!---')
#'/home/yan/IFS_FabuCode/InputData/IFS_sim_Fabu.config'
############################################################
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