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

update

parent 9dc4a371
Expand all Hide whitespace changes
Inline Side-by-side
csst_mci_sim/CTI/CTI.py
View file @ 3d2380f0
......@@ -15,7 +15,7 @@ parameters in parallel and serial direction.
import numpy as np
#CDM03bidir
# CDM03bidir
class CDM03bidir():
"""
Class to run CDM03 CTI model, class Fortran routine to perform the actual CDM03 calculations.
......@@ -27,6 +27,7 @@ class CDM03bidir():
:param log: instance to Python logging
:type log: logging instance
"""
def __init__(self, settings, data, log=None):
"""
Class constructor.
......@@ -39,48 +40,50 @@ class CDM03bidir():
:type log: logging instance
"""
self.data = data
self.values = dict(quads=(0,1,2,3), xsize=2048, ysize=2066, dob=0.0, rdose=8.0e9)
self.values = dict(quads=(0, 1, 2, 3), xsize=2048,
ysize=2066, dob=0.0, rdose=8.0e9)
self.values.update(settings)
self.log = log
self._setupLogger()
#default CDM03 settings
# default CDM03 settings
self.params = dict(beta_p=0.6, beta_s=0.6, fwc=200000., vth=1.168e7, vg=6.e-11, t=20.48e-3,
sfwc=730000., svg=1.0e-10, st=5.0e-6, parallel=1., serial=0.0)
#update with inputs
# update with inputs
self.params.update(self.values)
#read in trap information
trapdata = np.loadtxt(self.values['dir_path']+self.values['paralleltrapfile'])
# read in trap information
trapdata = np.loadtxt(
self.values['dir_path']+self.values['paralleltrapfile'])
if trapdata.ndim > 1:
self.nt_p = trapdata[:, 0]
self.sigma_p = trapdata[:, 1]
self.taur_p = trapdata[:, 2]
else:
#only one trap species
# only one trap species
self.nt_p = [trapdata[0],]
self.sigma_p = [trapdata[1],]
self.taur_p = [trapdata[2],]
trapdata = np.loadtxt(self.values['dir_path']+self.values['serialtrapfile'])
trapdata = np.loadtxt(
self.values['dir_path']+self.values['serialtrapfile'])
if trapdata.ndim > 1:
self.nt_s = trapdata[:, 0]
self.sigma_s = trapdata[:, 1]
self.taur_s = trapdata[:, 2]
else:
#only one trap species
# only one trap species
self.nt_s = [trapdata[0],]
self.sigma_s = [trapdata[1],]
self.taur_s = [trapdata[2],]
#scale thibaut's values
# scale thibaut's values
if 'thibaut' in self.values['parallelTrapfile']:
self.nt_p /= 0.576 #thibaut's values traps / pixel
self.sigma_p *= 1.e4 #thibaut's values in m**2
self.nt_p /= 0.576 # thibaut's values traps / pixel
self.sigma_p *= 1.e4 # thibaut's values in m**2
if 'thibaut' in self.values['serialTrapfile']:
self.nt_s *= 0.576 #thibaut's values traps / pixel #should be division?
self.sigma_s *= 1.e4 #thibaut's values in m**2
self.nt_s *= 0.576 # thibaut's values traps / pixel #should be division?
self.sigma_s *= 1.e4 # thibaut's values in m**2
def _setupLogger(self):
"""
......@@ -90,7 +93,6 @@ class CDM03bidir():
# if self.log is None:
# self.logger = False
def applyRadiationDamage(self, data, iquadrant=0):
"""
Apply radian damage based on FORTRAN CDM03 model. The method assumes that
......@@ -127,10 +129,10 @@ class CDM03bidir():
array will be laid out in memory in C-style (row-major order).
:return: image that has been run through the CDM03 model
:rtype: ndarray
"""""
#return data
:rtype: ndarray """""
# return data
iflip = iquadrant / 2
jflip = iquadrant % 2
......@@ -154,24 +156,22 @@ class CDM03bidir():
self.log.info('jflip=%i' % jflip)
#################################################################################
###modify
#sys.path.append('../so')
# modify
# sys.path.append('../so')
# from ifs_so.cdm03.cpython-38-x86_64-linux-gnu import cdm03bidir
# import cdm03bidir
from .mci_so import cdm03bidir
CTIed = cdm03bidir.cdm03(np.asfortranarray(data),
jflip, iflip,
self.values['dob'], self.values['rdose'],
self.nt_p, self.sigma_p, self.taur_p,
self.nt_s, self.sigma_s, self.taur_s,
params,
[data.shape[0], data.shape[1], len(self.nt_p), len(self.nt_s), len(self.params)])
jflip, iflip,
self.values['dob'], self.values['rdose'],
self.nt_p, self.sigma_p, self.taur_p,
self.nt_s, self.sigma_s, self.taur_s,
params,
[data.shape[0], data.shape[1], len(self.nt_p), len(self.nt_s), len(self.params)])
return np.asanyarray(CTIed)
#################################################################################################################
#################################################################################################################
csst_mci_sim/csst_mci_sim.py
View file @ 3d2380f0
This diff is collapsed.
csst_mci_sim/earthshine.py
View file @ 3d2380f0
......@@ -177,8 +177,8 @@ class StrayLight(object):
# band_star_e2 = star_e2[:][filterIndex[filter]]
# return max(band_star_e1, band_star_e2)
###############################################################################
###############################################################################
def caculateEarthShineFilter(self, filter='i'):
......@@ -207,7 +207,7 @@ class StrayLight(object):
###############################################################################
### test
# test
# path='/home/yan/MCI/'
# time_jd = 2460417.59979167
# x_sat = -4722.543136
......
csst_mci_sim/straylight.py
View file @ 3d2380f0
......@@ -205,7 +205,7 @@ class StrayLight(object):
###############################################################################
## test
# test
# path='/home/yan/MCI/'
# time_jd = 2460417.59979167
# x_sat = -4722.543136
......
csst_mci_sim/support/MCIinstrumentModel.py
View file @ 3d2380f0
......@@ -21,12 +21,12 @@ def MCIinformation():
"""
Returns a dictionary describing MCI CCD. The following information is provided (id: value - reference)::
dob: 0 - CDM03 (Short et al. 2010)
fwc: 90000 - CCD spec EUCL-EST-RS-6-002 (for CDM03)
fwc: 90000 - CCD spec EUCL-EST-RS-6-002 (for CDM03)
rdose: 30000000000.0 - derived (above the PLM requirement)
sfwc: 730000.0 - CDM03 (Short et al. 2010), see also the CCD spec EUCL-EST-RS-6-002
st: 5e-06 - CDM03 (Short et al. 2010)
sfwc: 730000.0 - CDM03 (Short et al. 2010), see also the CCD spec EUCL-EST-RS-6-002
st: 5e-06 - CDM03 (Short et al. 2010)
svg: 1e-10 - CDM03 (Short et al. 2010)
t: 0.01024 - CDM03 (Short et al. 2010)
trapfile: cdm_euclid.dat - CDM03 (derived, refitted to CCD204 data)
......@@ -40,19 +40,19 @@ def MCIinformation():
"""
#########################################################################################################
out=dict()
out.update({'dob' : 0, 'rdose' : 8.0e9,
'parallelTrapfile' : 'cdm_euclid_parallel.dat', 'serialTrapfile' : 'cdm_euclid_serial.dat',
'beta_s' : 0.6, 'beta_p': 0.6, 'fwc' : 90000, 'vth' : 1.168e7, 't' : 20.48e-3, 'vg' : 6.e-11,
'st' : 5.0e-6, 'sfwc' : 730000., 'svg' : 1.0e-10})
out = dict()
out.update({'dob': 0, 'rdose': 8.0e9,
'parallelTrapfile': 'cdm_euclid_parallel.dat', 'serialTrapfile': 'cdm_euclid_serial.dat',
'beta_s': 0.6, 'beta_p': 0.6, 'fwc': 90000, 'vth': 1.168e7, 't': 20.48e-3, 'vg': 6.e-11,
'st': 5.0e-6, 'sfwc': 730000., 'svg': 1.0e-10})
return out
# def CCDnonLinearityModel(data, beta=6e-7):
# """
# """
# The non-linearity is modelled based on the results presented.
# The non-linearity is modelled based on the results presented.
# :param data: data to which the non-linearity model is being applied to
# :type data: ndarray
......@@ -60,7 +60,7 @@ def MCIinformation():
# :rtype: float or ndarray
# """
# out = data-beta*data**2
# return out
###################################################################
......
csst_mci_sim/support/cosmicrays.py
View file @ 3d2380f0
......@@ -31,6 +31,7 @@ class cosmicrays():
:param information: cosmic ray track information (file containing track length and energy information) and
exposure time.
"""
def __init__(self, log, image, exptime, crInfo=None, information=None):
"""
Cosmic ray generation class. Can either draw events from distributions or
......@@ -42,9 +43,9 @@ class cosmicrays():
:param information: cosmic ray track information (file containing track length and energy information) and
exposure time.
"""
self.exptime=exptime
self.exptime = exptime
self.log = log
#image and size
......@@ -56,11 +57,11 @@ class cosmicrays():
else:
self._readCosmicrayInformation()
##################################################
def _cosmicRayIntercepts(self, lum, x0, y0, l, phi):
##################################################
############
############
def _cosmicRayIntercepts(self, lum, x0, y0, dl, phi):
"""
Derive cosmic ray streak intercept points.
......@@ -73,18 +74,18 @@ class cosmicrays():
:return: cosmic ray map (image)
:rtype: nd-array
"""
#create empty array
# create empty array
crImage = np.zeros((self.ysize, self.xsize), dtype=np.float64)
#x and y shifts
dx = l * np.cos(phi) / 2.
dy = l * np.sin(phi) / 2.
# x and y shifts
dx = dl * np.cos(phi) / 2.
dy = dl * np.sin(phi) / 2.
mskdx = np.abs(dx) < 1e-8
mskdy = np.abs(dy) < 1e-8
dx[mskdx] = 0.
dy[mskdy] = 0.
#pixels in x-direction
# pixels in x-direction
ilo = np.round(x0.copy() - dx)
msk = ilo < 0.
ilo[msk] = 0
......@@ -95,7 +96,7 @@ class cosmicrays():
ihi[msk] = self.xsize
ihi = ihi.astype(int)
#pixels in y-directions
# pixels in y-directions
jlo = np.round(y0.copy() - dy)
msk = jlo < 0.
jlo[msk] = 0
......@@ -106,7 +107,7 @@ class cosmicrays():
jhi[msk] = self.ysize
jhi = jhi.astype(int)
#loop over the individual events
# loop over the individual events
for i, luminosity in enumerate(lum):
n = 0 # count the intercepts
......@@ -114,7 +115,7 @@ class cosmicrays():
x = []
y = []
#Compute X intercepts on the pixel grid
# Compute X intercepts on the pixel grid
if ilo[i] < ihi[i]:
for xcoord in range(ilo[i], ihi[i]):
ok = (xcoord - x0[i]) / dx[i]
......@@ -132,7 +133,7 @@ class cosmicrays():
x.append(xcoord)
y.append(y0[i] + ok * dy[i])
#Compute Y intercepts on the pixel grid
# Compute Y intercepts on the pixel grid
if jlo[i] < jhi[i]:
for ycoord in range(jlo[i], jhi[i]):
ok = (ycoord - y0[i]) / dy[i]
......@@ -150,13 +151,13 @@ class cosmicrays():
x.append(x0[i] + ok * dx[i])
y.append(ycoord)
#check if no intercepts were found
# check if no intercepts were found
if n < 1:
xc = int(np.floor(x0[i]))
yc = int(np.floor(y0[i]))
crImage[yc, xc] += luminosity
#Find the arguments that sort the intersections along the track
# Find the arguments that sort the intersections along the track
u = np.asarray(u)
x = np.asarray(x)
y = np.asarray(y)
......@@ -167,7 +168,7 @@ class cosmicrays():
x = x[args]
y = y[args]
#Decide which cell each interval traverses, and the path length
# Decide which cell each interval traverses, and the path length
for i in range(1, n - 1):
w = u[i + 1] - u[i]
cx = int(1 + np.floor((x[i + 1] + x[i]) / 2.))
......@@ -177,11 +178,11 @@ class cosmicrays():
crImage[cy, cx] += (w * luminosity)
return crImage
##################################################
############################################
############################################################################
#####################################
########################################
def _drawEventsToCoveringFactor(self, coveringFraction=3.0, limit=1000, verbose=False):
"""
Generate cosmic ray events up to a covering fraction and include it to a cosmic ray map (self.cosmicrayMap).
......@@ -198,50 +199,54 @@ class cosmicrays():
"""
self.cosmicrayMap = np.zeros((self.ysize, self.xsize))
#how many events to draw at once, too large number leads to exceeding the covering fraction
####cr_n = int(295 * self.exptime / 565. * coveringFraction / 1.4)
# how many events to draw at once, too large number leads to exceeding the covering fraction
# cr_n = int(295 * self.exptime / 565. * coveringFraction / 1.4)
cr_n = int(5000 * self.exptime / 565. * coveringFraction)
covering = 0.0
while covering < coveringFraction:
#pseudo-random numbers taken from a uniform distribution between 0 and 1
# pseudo-random numbers taken from a uniform distribution between 0 and 1
np.random.seed()
luck = np.random.rand(cr_n)
#draw the length of the tracks
ius = InterpolatedUnivariateSpline(self.cr['cr_cdf'], self.cr['cr_u'])
# draw the length of the tracks
ius = InterpolatedUnivariateSpline(
self.cr['cr_cdf'], self.cr['cr_u'])
self.cr['cr_l'] = ius(luck)
if limit is None:
ius = InterpolatedUnivariateSpline(self.cr['cr_cde'], self.cr['cr_v'])
ius = InterpolatedUnivariateSpline(
self.cr['cr_cde'], self.cr['cr_v'])
self.cr['cr_e'] = ius(luck)
else:
#set the energy directly to the limit
# set the energy directly to the limit
self.cr['cr_e'] = np.asarray([limit,])
#Choose the properties such as positions and an angle from a random Uniform dist
# Choose the properties such as positions and an angle from a random Uniform dist
np.random.seed()
cr_x = self.xsize * np.random.rand(int(np.floor(cr_n)))
np.random.seed()
cr_y = self.ysize * np.random.rand(int(np.floor(cr_n)))
np.random.seed()
cr_phi = np.pi * np.random.rand(int(np.floor(cr_n)))
#find the intercepts
self.cosmicrayMap += self._cosmicRayIntercepts(self.cr['cr_e'], cr_x, cr_y, self.cr['cr_l'], cr_phi)
# find the intercepts
self.cosmicrayMap += self._cosmicRayIntercepts(
self.cr['cr_e'], cr_x, cr_y, self.cr['cr_l'], cr_phi)
#count the covering factor
# count the covering factor
area_cr = np.count_nonzero(self.cosmicrayMap)
covering = 100.*area_cr / (self.xsize*self.ysize)
text = 'The cosmic ray covering factor is %i pixels i.e. %.3f per cent' % (area_cr, covering)
text = 'The cosmic ray covering factor is %i pixels i.e. %.3f per cent' % (
area_cr, covering)
self.log.info(text)
###################################################33
# 33
def addUpToFraction(self, coveringFraction, limit=None, verbose=False):
"""
......@@ -257,13 +262,10 @@ class cosmicrays():
:return: image with cosmic rays
:rtype: ndarray
"""
self._drawEventsToCoveringFactor(coveringFraction, limit=limit, verbose=verbose)
self._drawEventsToCoveringFactor(
coveringFraction, limit=limit, verbose=verbose)
#paste cosmic rays
# paste cosmic rays
self.image += self.cosmicrayMap
return self.image
csst_mci_sim/support/logger.py
View file @ 3d2380f0
......@@ -23,16 +23,17 @@ def setUpLogger(log_filename, loggername='logger'):
logger.setLevel(logging.DEBUG)
# Add the log message handler to the logger
handler = logging.handlers.RotatingFileHandler(log_filename)
#maxBytes=20, backupCount=5)
# maxBytes=20, backupCount=5)
# create formatter
formatter = logging.Formatter('%(asctime)s - %(module)s - %(funcName)s - %(levelname)s - %(message)s')
formatter = logging.Formatter(
'%(asctime)s - %(module)s - %(funcName)s - %(levelname)s - %(message)s')
# add formatter to ch
handler.setFormatter(formatter)
# add handler to logger
# add handler to logger
if (logger.hasHandlers()):
logger.handlers.clear()
logger.addHandler(handler)
return logger
......
csst_mci_sim/support/sed.py
View file @ 3d2380f0
......@@ -15,7 +15,8 @@ from astropy import units as u
from scipy.interpolate import interp1d
def Calzetti_Law(wave, Rv = 4.05):
def Calzetti_Law(wave, Rv=4.05):
"""Dust Extinction Curve by Calzetti et al. (2000)
Args:
......@@ -25,19 +26,20 @@ def Calzetti_Law(wave, Rv = 4.05):
Returns:
float: Extinction value E(B-V)
"""
wave_number = 1./(wave * 1e-4)
reddening_curve = np.zeros(len(wave))
idx = np.logical_and(wave >= 1200, wave <= 6300)
reddening_curve[idx] = 2.659 * ( -2.156 + 1.509 * wave_number[idx] - 0.198 * \
(wave_number[idx] ** 2)) + 0.011 * (wave_number[idx] **3 ) + Rv
reddening_curve[idx] = 2.659 * (-2.156 + 1.509 * wave_number[idx] - 0.198 *
(wave_number[idx] ** 2)) + 0.011 * (wave_number[idx] ** 3) + Rv
idx = np.logical_and(wave >= 6300, wave <= 22000)
reddening_curve[idx] = 2.659 * ( -1.857 + 1.040 * wave_number[idx]) + Rv
reddening_curve[idx] = 2.659 * (-1.857 + 1.040 * wave_number[idx]) + Rv
return reddening_curve
def reddening(wave, flux, ebv = 0.0, law = 'calzetti', Rv = 4.05):
def reddening(wave, flux, ebv=0.0, law='calzetti', Rv=4.05):
"""
Reddening an input spectra through a given reddening curve.
......@@ -52,10 +54,11 @@ def reddening(wave, flux, ebv = 0.0, law = 'calzetti', Rv = 4.05):
float: Flux of spectra after reddening.
"""
if law == 'calzetti':
curve = Calzetti_Law(wave, Rv = Rv)
curve = Calzetti_Law(wave, Rv=Rv)
fluxNew = flux / (10. ** (0.4 * ebv * curve))
return fluxNew
def flux_to_mag(wave, flux, path, band='GAIA_bp'):
"""Convert flux of given spectra to magnitude
......@@ -68,14 +71,14 @@ def flux_to_mag(wave, flux, path, band='GAIA_bp'):
float: value of magnitude
"""
# /home/yan/MCI_sim/MCI_input/SED_Code/data
##import os
###parent = os.path.dirname(os.path.realpath(__file__))
#
# parent = os.path.dirname(os.path.realpath(__file__))
band = ascii.read(path+'MCI_inputData/SED_Code/seddata/' + band + '.dat')
wave0= band['col1']
curv0= band['col2']
wave0 = band['col1']
curv0 = band['col2']
# Setting the response
func = interp1d(wave0, curv0)
response = np.copy(wave)
......@@ -86,10 +89,11 @@ def flux_to_mag(wave, flux, path, band='GAIA_bp'):
# Total Flux
Tflux = np.trapz(flux * response, wave) / np.trapz(response, wave)
return -2.5 * np.log10(Tflux)
def calibrate(wave, flux, mag, path,band='GAIA_bp'):
def calibrate(wave, flux, mag, path, band='GAIA_bp'):
"""
Calibrate the spectra according to the magnitude.
......@@ -102,78 +106,83 @@ def calibrate(wave, flux, mag, path,band='GAIA_bp'):
Returns:
float: Flux of calibrated spectra. Units: 1e-17 erg/s/A/cm^2
"""
inst_mag = flux_to_mag(wave, flux, path,band = band)
inst_mag = flux_to_mag(wave, flux, path, band=band)
instflux = 10 ** (-0.4 * inst_mag)
realflux = (mag * u.STmag).to(u.erg/u.s/u.cm**2/u.AA).value
# Normalization
flux_ratio = realflux / instflux
flux_calibrate = flux * flux_ratio * 1e17 # Units: 10^-17 erg/s/A/cm^2
# Units: 10^-17 erg/s/A/cm^2
flux_calibrate = flux * flux_ratio * 1e17
return flux_calibrate
# ------------
# SED Template
class Gal_Temp():
"""
Template of Galaxy SED
"""
def __init__(self,path):
###import os
###parent = os.path.dirname(os.path.realpath(__file__))
self.path=path
hdulist = fits.open(self.path+'MCI_inputData/SED_Code/seddata/galaxy_temp.fits')
def __init__(self, path):
#
# parent = os.path.dirname(os.path.realpath(__file__))
self.path = path
hdulist = fits.open(
self.path+'MCI_inputData/SED_Code/seddata/galaxy_temp.fits')
self.wave = hdulist[1].data['wave']
self.flux = hdulist[2].data
self.age_grid = hdulist[3].data['logAge']
self.feh_grid = hdulist[3].data['FeH']
def toMag(self, redshift = 0):
def toMag(self, redshift=0):
"""Calculating magnitude
Args:
redshift (float, optional): redshift of spectra. Defaults to 0.
"""
wave = self.wave * (1 + redshift)
self.umag = flux_to_mag(wave, self.flux, self.path,band='SDSS_u')
self.gmag = flux_to_mag(wave, self.flux, self.path,band='SDSS_g')
self.rmag = flux_to_mag(wave, self.flux, self.path,band='SDSS_r')
self.imag = flux_to_mag(wave, self.flux, self.path,band='SDSS_i')
self.zmag = flux_to_mag(wave, self.flux, self.path,band='SDSS_z')
self.umag = flux_to_mag(wave, self.flux, self.path, band='SDSS_u')
self.gmag = flux_to_mag(wave, self.flux, self.path, band='SDSS_g')
self.rmag = flux_to_mag(wave, self.flux, self.path, band='SDSS_r')
self.imag = flux_to_mag(wave, self.flux, self.path, band='SDSS_i')
self.zmag = flux_to_mag(wave, self.flux, self.path, band='SDSS_z')
class Star_Temp():
"""
Template of Stellar SED
"""
def __init__(self,path):
##import os
self.path=path
####parent = os.path.dirname(os.path.realpath(__file__))
###print("获取其父目录——" + parent) # 从当前文件路径中获取目录
hdulist = fits.open(path+'MCI_inputData/SED_Code/seddata/stellar_temp.fits')
def __init__(self, path):
#
self.path = path
# parent = os.path.dirname(os.path.realpath(__file__))
# print("获取其父目录——" + parent) # 从当前文件路径中获取目录
hdulist = fits.open(
path+'MCI_inputData/SED_Code/seddata/stellar_temp.fits')
self.wave = hdulist[1].data['wave']
self.flux = hdulist[2].data
self.Teff_grid = hdulist[3].data['Teff']
self.FeH_grid = hdulist[3].data['FeH']
self.bpmag = flux_to_mag(self.wave, self.flux, path,band='GAIA_bp')
self.rpmag = flux_to_mag(self.wave, self.flux, path,band='GAIA_rp')
self.bpmag = flux_to_mag(self.wave, self.flux, path, band='GAIA_bp')
self.rpmag = flux_to_mag(self.wave, self.flux, path, band='GAIA_rp')
def toMag(self):
wave = self.wave
self.bpmag = flux_to_mag(wave, self.flux, self.path,band='GAIA_bp')
self.rpmag = flux_to_mag(wave, self.flux, self.path,band='GAIA_rp')
self.bpmag = flux_to_mag(wave, self.flux, self.path, band='GAIA_bp')
self.rpmag = flux_to_mag(wave, self.flux, self.path, band='GAIA_rp')
# -------------
# SED Modelling
def Model_Stellar_SED(wave, bp, rp, temp):
"""Modelling stellar SED based on bp, rp magnitude
......@@ -184,19 +193,20 @@ def Model_Stellar_SED(wave, bp, rp, temp):
temp (class): Class of stellar template
Returns:
float array: Spectral energy distribution of stellar SED,
float array: Spectral energy distribution of stellar SED,
which have the same length to the input wave
"""
"""
color0 = bp - rp
colors = temp.bpmag - temp.rpmag
idx = np.argmin(np.abs(colors - color0))
flux0 = temp.flux[idx]
flux1 = np.interp(wave, temp.wave, flux0)
flux = calibrate(wave, flux1, rp, band = 'GAIA_rp')
flux = calibrate(wave, flux1, rp, band='GAIA_rp')
return flux
def Model_Galaxy_SED(wave, ugriz, z, temp, path):
"""Modelling galaxy SED based on u,g,r,i,z magnitude
......@@ -213,30 +223,30 @@ def Model_Galaxy_SED(wave, ugriz, z, temp, path):
sed = 10. ** (-0.4 * ugriz)
sed = sed / sed[2]
ntemp = len(temp.rmag)
dmag = np.zeros(ntemp)
dmag = np.zeros(ntemp)
for j in range(ntemp):
ugriz0 = np.array([temp.umag[j], temp.gmag[j], temp.rmag[j],
ugriz0 = np.array([temp.umag[j], temp.gmag[j], temp.rmag[j],
temp.imag[j], temp.zmag[j]])
sed0 = 10. ** (-0.4 * ugriz0)
sed0 = sed0 / sed0[2]
dmag[j] = np.sum(np.abs(sed - sed0))
idx = np.argmin(dmag)
flux0 = temp.flux[idx]
# Effect of E(B-V)
ri0 = ugriz[2] - ugriz[3]
ri = temp.rmag - temp.imag
ri = temp.rmag - temp.imag
dri = ri0 - ri[idx]
Alambda = Calzetti_Law(np.array([6213 / (1 + z), 7625 / (1 + z)]))
eri0 = (Alambda[0] - Alambda[1])
ebv = dri/eri0
if ebv<0:
ebv=0
if ebv>0.5:
ebv=0.5
flux1 = reddening(temp.wave, flux0, ebv = ebv)
if ebv < 0:
ebv = 0
if ebv > 0.5:
ebv = 0.5
flux1 = reddening(temp.wave, flux0, ebv=ebv)
flux2 = np.interp(wave, temp.wave * (1 + z), flux1)
flux = calibrate(wave, flux2, ugriz[2], path, band = 'SDSS_r')
return flux
\ No newline at end of file
flux = calibrate(wave, flux2, ugriz[2], path, band='SDSS_r')
return flux
csst_mci_sim/support/shao.py
View file @ 3d2380f0
from ctypes import *
def checkInputList(input_list, n):
if type(input_list) != type([1, 2, 3]):
if not isinstance(type(input_list), list):
# if type(input_list) != type([1, 2, 3]):
raise TypeError("Input type is not list!", input_list)
for i in input_list:
if type(i) != type(1.1):
if type(i) != type(1):
raise TypeError("Input list's element is not float or int!", input_list)
if not isinstance(type(i), float): # type(i) != type(1.1):
if not isinstance(type(i), int): # type(i) != type(1):
raise TypeError(
"Input list's element is not float or int!", input_list)
if len(input_list) != n:
raise RuntimeError("Length of input list is not equal to stars' number!", input_list)
def onOrbitObsPosition(path, input_ra_list, input_dec_list, input_pmra_list, input_pmdec_list, input_rv_list, \
input_parallax_list, input_nstars, input_x, input_y, input_z, input_vx, input_vy, \
raise RuntimeError(
"Length of input list is not equal to stars' number!", input_list)
def onOrbitObsPosition(path, input_ra_list, input_dec_list, input_pmra_list, input_pmdec_list, input_rv_list,
input_parallax_list, input_nstars, input_x, input_y, input_z, input_vx, input_vy,
input_vz, input_epoch, input_date_str, input_time_str):
#Check input parameters
if type(input_nstars) != type(1):
# Check input parameters
if not isinstance(type(input_nstars), int): # type(input_nstars) != type(1):
raise TypeError("Parameter 7 is not int!", input_nstars)
checkInputList(input_ra_list, input_nstars)
checkInputList(input_dec_list, input_nstars)
checkInputList(input_pmra_list, input_nstars)
checkInputList(input_pmdec_list, input_nstars)
checkInputList(input_rv_list, input_nstars)
checkInputList(input_parallax_list, input_nstars)
if type(input_x) != type(1.1):
if not isinstance(type(input_x), float): # type(input_x) != type(1.1):
raise TypeError("Parameter 8 is not double!", input_x)
if type(input_y) != type(1.1):
if not isinstance(type(input_y), float): # type(input_y) != type(1.1):
raise TypeError("Parameter 9 is not double!", input_y)
if type(input_z) != type(1.1):
if not isinstance(type(input_z), float): # type(input_z) != type(1.1):
raise TypeError("Parameter 10 is not double!", input_z)
if type(input_vx) != type(1.1):
if not isinstance(type(input_vx), float): # type(input_vx) != type(1.1):
raise TypeError("Parameter 11 is not double!", input_vx)
if type(input_vy) != type(1.1):
if not isinstance(type(input_vy), float): # type(input_vy) != type(1.1):
raise TypeError("Parameter 12 is not double!", input_vy)
if type(input_vz) != type(1.1):
if not isinstance(type(input_vz), float): # type(input_vz) != type(1.1):
raise TypeError("Parameter 13 is not double!", input_vz)
#Convert km -> m
input_x = input_x*1000.0
input_y = input_y*1000.0
# Convert km -> m
input_x = input_x*1000.0
input_y = input_y*1000.0
input_z = input_z*1000.0
input_vx = input_vx*1000.0
input_vy = input_vy*1000.0
input_vz = input_vz*1000.0
if type(input_date_str) != type("2025-03-05"):
input_vz = input_vz*1000.0
if not isinstance(type(input_date_str), str): # type(input_date_str) != type("2025-03-05"):
raise TypeError("Parameter 15 is not string!", input_date_str)
else:
input_date_str = input_date_str.strip()
if not (input_date_str[4]=="-" and input_date_str[7]=="-"):
if not (input_date_str[4] == "-" and input_date_str[7] == "-"):
raise TypeError("Parameter 15 format error (1)!", input_date_str)
else:
tmp = input_date_str.split("-")
if len(tmp) != 3:
raise TypeError("Parameter 15 format error (2)!", input_date_str)
raise TypeError(
"Parameter 15 format error (2)!", input_date_str)
input_year = int(tmp[0])
input_month = int(tmp[1])
input_day = int(tmp[2])
if not (input_year>=1900 and input_year<=2100):
raise TypeError("Parameter 15 year range error [1900 ~ 2100]!", input_year)
if not (input_month>=1 and input_month<=12):
raise TypeError("Parameter 15 month range error [1 ~ 12]!", input_month)
if not (input_day>=1 and input_day<=31):
raise TypeError("Parameter 15 day range error [1 ~ 31]!", input_day)
if type(input_time_str) != type("20:15:15.15"):
if not (input_year >= 1900 and input_year <= 2100):
raise TypeError(
"Parameter 15 year range error [1900 ~ 2100]!", input_year)
if not (input_month >= 1 and input_month <= 12):
raise TypeError(
"Parameter 15 month range error [1 ~ 12]!", input_month)
if not (input_day >= 1 and input_day <= 31):
raise TypeError(
"Parameter 15 day range error [1 ~ 31]!", input_day)
if not isinstance(type(input_time_str), str): # type(input_time_str) != type("20:15:15.15"):
raise TypeError("Parameter 16 is not string!", input_time_str)
else:
input_time_str = input_time_str.strip()
if not (input_time_str[2]==":" and input_time_str[5]==":"):
if not (input_time_str[2] == ":" and input_time_str[5] == ":"):
raise TypeError("Parameter 16 format error (1)!", input_time_str)
else:
tmp = input_time_str.split(":")
if len(tmp) != 3:
raise TypeError("Parameter 16 format error (2)!", input_time_str)
raise TypeError(
"Parameter 16 format error (2)!", input_time_str)
input_hour = int(tmp[0])
input_minute = int(tmp[1])
input_second = float(tmp[2])
if not (input_hour>=0 and input_hour<=23):
raise TypeError("Parameter 16 hour range error [0 ~ 23]!", input_hour)
if not (input_minute>=0 and input_minute<=59):
raise TypeError("Parameter 16 minute range error [0 ~ 59]!", input_minute)
if not (input_second>=0 and input_second<60.0):
raise TypeError("Parameter 16 second range error [0 ~ 60)!", input_second)
#Inital dynamic lib
if not (input_hour >= 0 and input_hour <= 23):
raise TypeError(
"Parameter 16 hour range error [0 ~ 23]!", input_hour)
if not (input_minute >= 0 and input_minute <= 59):
raise TypeError(
"Parameter 16 minute range error [0 ~ 59]!", input_minute)
if not (input_second >= 0 and input_second < 60.0):
raise TypeError(
"Parameter 16 second range error [0 ~ 60)!", input_second)
# Inital dynamic lib
import os
currfile=os.getcwd()
currfile = os.getcwd()
print(currfile)
shao = cdll.LoadLibrary(path+'MCI_inputData/TianCe/libshao.so')
shao.onOrbitObs.restype = c_int
d3 = c_double * 3
shao.onOrbitObs.argtypes = [c_double, c_double, c_double, c_double, c_double, c_double, \
c_int, c_int, c_int, c_int, c_int, c_double, \
c_double, POINTER(d3), POINTER(d3), \
c_int, c_int, c_int, c_int, c_int, c_double, \
POINTER(c_double), POINTER(c_double) ]
shao.onOrbitObs.argtypes = [c_double, c_double, c_double, c_double, c_double, c_double,
c_int, c_int, c_int, c_int, c_int, c_double,
c_double, POINTER(d3), POINTER(d3),
c_int, c_int, c_int, c_int, c_int, c_double,
POINTER(c_double), POINTER(c_double)]
output_ra_list = list()
output_dec_list = list()
for i in range(input_nstars):
......@@ -113,22 +125,21 @@ def onOrbitObsPosition(path, input_ra_list, input_dec_list, input_pmra_list, inp
v3 = d3(input_vx, input_vy, input_vz)
input_year_c = c_int(input_year)
input_month_c = c_int(input_month)
input_day_c = c_int(input_day)
input_day_c = c_int(input_day)
input_hour_c = c_int(input_hour)
input_minute_c = c_int(input_minute)
input_second_c = c_double(input_second)
input_second_c = c_double(input_second)
DAT = c_double(37.0)
output_ra = c_double(0.0)
output_dec = c_double(0.0)
rs = shao.onOrbitObs(input_ra, input_dec, input_parallax, input_pmra, input_pmdec, input_rv, \
input_year_c, input_month_c, input_day_c, input_hour_c, input_minute_c, input_second_c, \
DAT, byref(p3), byref(v3), \
input_year_c, input_month_c, input_day_c, input_hour_c, input_minute_c, input_second_c, \
rs = shao.onOrbitObs(input_ra, input_dec, input_parallax, input_pmra, input_pmdec, input_rv,
input_year_c, input_month_c, input_day_c, input_hour_c, input_minute_c, input_second_c,
DAT, byref(p3), byref(v3),
input_year_c, input_month_c, input_day_c, input_hour_c, input_minute_c, input_second_c,
byref(output_ra), byref(output_dec))
if rs != 0:
raise RuntimeError("Calculate error!")
output_ra_list.append(output_ra.value)
output_dec_list.append(output_dec.value)
return output_ra_list, output_dec_list
\ No newline at end of file
csst_mci_sim/zodiacal.py
View file @ 3d2380f0
......@@ -9,6 +9,7 @@ from astropy import units as u
from scipy import interpolate
def zodiacal(ra, dec, time, path):
"""
For given RA, DEC and TIME, return the interpolated zodical spectrum in Leinert-1998.
......@@ -19,7 +20,7 @@ def zodiacal(ra, dec, time, path):
:param path: the relative file path
:return:
wave_A: wavelength of the zodical spectrum
spec_erg: flux of the zodical spectrum, in unit of erg/s/cm^2/A/sr
"""
......@@ -27,14 +28,15 @@ def zodiacal(ra, dec, time, path):
# get solar position
dt = datetime.fromisoformat(time)
jd = julian.to_jd(dt, fmt='jd')
###jd = time2jd(dt)
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')
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.
......@@ -45,16 +47,18 @@ def zodiacal(ra, dec, time, path):
lamda = abs(lb_obj.lon.degree - lb_sun.lon.degree)
# interpolated zodical surface brightness at 0.5 um
zodi = pd.read_csv(path+'MCI_inputData/refs/zodi_map.dat', sep='\s+', header=None, comment='#')
zodi = pd.read_csv(path+'MCI_inputData/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])
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^�? W m�? sr�? um�?
# read the zodical spectrum in the ecliptic
cat_spec = pd.read_csv(path+'MCI_inputData/refs/solar_spec.dat', sep='\s+', header=None, comment='#')
cat_spec = pd.read_csv(
path+'MCI_inputData/refs/solar_spec.dat', sep='\s+', header=None, comment='#')
wave = cat_spec[0].values # A
spec0 = cat_spec[1].values # 10^-8 W m^�? sr^�? μm^�?
zodi_norm = 252 # 10^-8 W m^�? sr^�? μm^�?
......@@ -63,10 +67,11 @@ def zodiacal(ra, dec, time, path):
# 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_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
# erg/s/cm^2/A/arcsec^2
spec_erg2 = spec_erg / 4.25452e10
return wave_A, spec_erg2
......
tests/test_mci_sim.py
View file @ 3d2380f0
......@@ -15,9 +15,9 @@ import sys
import faulthandler
from csst_mci_sim import csst_mci_sim
class TestDemoFunction(unittest.TestCase):
def test_mci_sim_1(self):
"""
Aim
---
......@@ -33,76 +33,77 @@ class TestDemoFunction(unittest.TestCase):
This case aims to test whether the demo function returns `1` if input is `None`.
"""
faulthandler.enable()
# demo function test
dir_path = os.path.join(os.environ['UNIT_TEST_DATA_ROOT'],'mci_sim/')
dir_path = os.path.join(os.environ['UNIT_TEST_DATA_ROOT'], 'mci_sim/')
print(dir_path)
# 获取当前工作目录
# current_path = os.getcwd()
# current_path = os.getcwd()
# print("当前路径:", current_path)
configfile = './csst_mci_sim/mci_data/mci_all_9K.config'
sourcein = 'EXDF'
print(configfile)
debug=True
result_path = dir_path +'mci_sim_result/'
csst_mci_sim.runMCIsim(sourcein, configfile, dir_path, result_path, debug, 1)
debug = True
result_path = dir_path + 'mci_sim_result/'
csst_mci_sim.runMCIsim(sourcein, configfile,
dir_path, result_path, debug, 1)
self.assertEqual(
1 , 1,
1, 1,
"case 1: EXDF sim passes.",
)
############################################
def test_mci_sim_2(self):
"""
Aim
---
Test mci sim function: STAR case.
Criteria
--------
Pass if the demo function returns `1`.
Details
-------
The demo function returns the length of the input argument list.
This case aims to test whether the demo function returns `1` if input is `None`.
"""
faulthandler.enable()
# demo function test
dir_path = os.path.join(os.environ['UNIT_TEST_DATA_ROOT'],'mci_sim/')
dir_path = os.path.join(os.environ['UNIT_TEST_DATA_ROOT'], 'mci_sim/')
print(dir_path)
# 获取当前工作目录
# current_path = os.getcwd()
# current_path = os.getcwd()
# print("当前路径:", current_path)
configfile = './csst_mci_sim/mci_data/mci_all_9K.config'
sourcein = 'STAR'
print(configfile)
debug=True
result_path = dir_path +'mci_sim_result/'
csst_mci_sim.runMCIsim(sourcein, configfile, dir_path, result_path, debug, 1)
debug = True
result_path = dir_path + 'mci_sim_result/'
csst_mci_sim.runMCIsim(sourcein, configfile,
dir_path, result_path, debug, 1)
self.assertEqual(
1 , 1,
1, 1,
"case 2: STAR sim passes.",
)
)
#########################################################
def test_mci_sim_3(self):
"""
Aim
---
......@@ -118,33 +119,34 @@ class TestDemoFunction(unittest.TestCase):
This case aims to test whether the demo function returns `1` if input is `None`.
"""
faulthandler.enable()
# demo function test
dir_path = os.path.join(os.environ['UNIT_TEST_DATA_ROOT'],'mci_sim/')
dir_path = os.path.join(os.environ['UNIT_TEST_DATA_ROOT'], 'mci_sim/')
print(dir_path)
# 获取当前工作目录
# current_path = os.getcwd()
# current_path = os.getcwd()
# print("当前路径:", current_path)
configfile = './csst_mci_sim/mci_data/mci_all_9K.config'
sourcein = 'BIAS'
print(configfile)
debug=True
result_path = dir_path +'mci_sim_result/'
csst_mci_sim.runMCIsim(sourcein, configfile, dir_path, result_path, debug, 1)
debug = True
result_path = dir_path + 'mci_sim_result/'
csst_mci_sim.runMCIsim(sourcein, configfile,
dir_path, result_path, debug, 1)
self.assertEqual(
1 , 1,
1, 1,
"case 3: BIAS sim passes.",
)
#########################################################
def test_mci_sim_4(self):
"""
Aim
---
......@@ -160,34 +162,34 @@ class TestDemoFunction(unittest.TestCase):
This case aims to test whether the demo function returns `1` if input is `None`.
"""
faulthandler.enable()
# demo function test
dir_path = os.path.join(os.environ['UNIT_TEST_DATA_ROOT'],'mci_sim/')
dir_path = os.path.join(os.environ['UNIT_TEST_DATA_ROOT'], 'mci_sim/')
print(dir_path)
# 获取当前工作目录
# current_path = os.getcwd()
# current_path = os.getcwd()
# print("当前路径:", current_path)
configfile = './csst_mci_sim/mci_data/mci_all_9K.config'
sourcein = 'DARK'
print(configfile)
debug=True
result_path = dir_path +'mci_sim_result/'
csst_mci_sim.runMCIsim(sourcein, configfile, dir_path, result_path, debug, 1)
debug = True
result_path = dir_path + 'mci_sim_result/'
csst_mci_sim.runMCIsim(sourcein, configfile,
dir_path, result_path, debug, 1)
self.assertEqual(
1 , 1,
1, 1,
"case 4: DARK sim passes.",
)
)
#########################################################
def test_mci_sim_5(self):
"""
Aim
---
......@@ -203,32 +205,30 @@ class TestDemoFunction(unittest.TestCase):
This case aims to test whether the demo function returns `1` if input is `None`.
"""
faulthandler.enable()
# demo function test
dir_path = os.path.join(os.environ['UNIT_TEST_DATA_ROOT'],'mci_sim/')
dir_path = os.path.join(os.environ['UNIT_TEST_DATA_ROOT'], 'mci_sim/')
print(dir_path)
# 获取当前工作目录
# current_path = os.getcwd()
# current_path = os.getcwd()
# print("当前路径:", current_path)
configfile = './csst_mci_sim/mci_data/mci_all_9K.config'
configfile = './csst_mci_sim/mci_data/mci_all_9K.config'
sourcein = 'FLAT'
print(configfile)
debug=True
result_path = dir_path +'mci_sim_result/'
csst_mci_sim.runMCIsim(sourcein, configfile, dir_path, result_path, debug, 1)
debug = True
result_path = dir_path + 'mci_sim_result/'
csst_mci_sim.runMCIsim(sourcein, configfile,
dir_path, result_path, debug, 1)
self.assertEqual(
1 , 1,
1, 1,
"case 5: FLAT sim passes.",
)
# ############################################################################
\ No newline at end of file
# ############################################################################
tests/test_zodiacal.py
View file @ 3d2380f0
......@@ -42,7 +42,7 @@ class TestDemoFunction(unittest.TestCase):
# current_path = os.getcwd()
# print("当前路径:", current_path)
wave0, zodi0=zodiacal.zodiacal(10.0, 20.0, '2024-04-04', dir_path)
wave0, zodi0 = zodiacal.zodiacal(10.0, 20.0, '2024-04-04', dir_path)
self.assertEqual(
1, 1,
......
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