diff --git a/csst_mci_sim/CTI/__pycache__/CTI.cpython-311.pyc b/csst_mci_sim/CTI/__pycache__/CTI.cpython-311.pyc new file mode 100644 index 0000000000000000000000000000000000000000..a983316360d3e6a00ecc939a20238bbe7c6c991a Binary files /dev/null and b/csst_mci_sim/CTI/__pycache__/CTI.cpython-311.pyc differ diff --git a/csst_mci_sim/CTI/__pycache__/__init__.cpython-311.pyc b/csst_mci_sim/CTI/__pycache__/__init__.cpython-311.pyc new file mode 100644 index 0000000000000000000000000000000000000000..caed98cd22b974dc7d35c35afc730601a59b2f07 Binary files /dev/null and b/csst_mci_sim/CTI/__pycache__/__init__.cpython-311.pyc differ diff --git a/csst_mci_sim/csst_mci_sim.py b/csst_mci_sim/csst_mci_sim.py index de3a915f03a44b27a375dabefffc3a5565c31680..5b1c269c40177cc29917589c0d3adde1f56a4a9a 100644 --- a/csst_mci_sim/csst_mci_sim.py +++ b/csst_mci_sim/csst_mci_sim.py @@ -61,8 +61,8 @@ from datetime import datetime, timedelta from astropy.time import Time from astropy.coordinates import get_sun -import scipy.io as sio -import cmath +# import scipy.io as sio +# import cmath import numpy as np @@ -70,7 +70,7 @@ from tqdm import tqdm from astropy.table import Table -from astropy.wcs import WCS as WCS +# from astropy.wcs import WCS as WCS from astropy.io import fits @@ -80,9 +80,9 @@ import os, sys, math import configparser as ConfigParser -from optparse import OptionParser +# from optparse import OptionParser -from matplotlib import pyplot as plt +# from matplotlib import pyplot as plt from scipy import ndimage @@ -115,20 +115,6 @@ import pandas # sys.path.append('../MCI_inputData/SED_Code') - - - - -# import jax -#import jax.numpy as jnp - -# from jax import config -# config.update("jax_enable_x64", True) - -# os.environ['CUDA_VISIBLE_DEVICES'] = '0' -# devices = jax.local_devices() - - # set the folder FOLDER ='../' @@ -141,19 +127,6 @@ import astropy.coordinates as coord from scipy.interpolate import interp1d -#from sky_bkg import sky_bkg - - -filterPivotWave = {'nuv':2875.5,'u':3629.6,'g':4808.4,'r':6178.2, 'i':7609.0, 'z':9012.9,'y':9627.9} -filterIndex = {'nuv':0,'u':1,'g':2,'r':3, 'i':4, 'z':5,'y':6} -# filterCCD = {'nuv':'UV0','u':'UV0','g':'Astro_MB','r':'Astro_MB', 'i':'Basic_NIR', 'z':'Basic_NIR','y':'Basic_NIR'} -# bandRange = {'nuv':[2504.0,3230.0],'u':[3190.0,4039.0],'g':[3989.0,5498.0],'r':[5438.0,6956.0], 'i':[6886.0,8469.0], -# 'z':[8379.0,10855.0],'y':[9217.0, 10900.0], 'GU':[2550, 4000],'GV':[4000, 6200],'GI':[6200,10000]} -# Instrument_dir = '/Users/linlin/Data/csst/straylightsim-master/Instrument/' -# SpecOrder = ['-2','-1','0','1','2'] -# -# filterMirrorEff = {'nuv':0.54,'u':0.68,'g':0.8,'r':0.8, 'i':0.8, 'z':0.8,'y':0.8} - def transRaDec2D(ra, dec): # radec转为竞天程序里的ob, 赤道坐标系下的笛卡尔三维坐标xyz. @@ -163,53 +136,6 @@ def transRaDec2D(ra, dec): return np.array([x1, y1, z1]) ############################################################################### -# def flux2ill(wave, flux): - -# # erg/s/cm^2/A/arcsec^2 to W/m^2 - -# # 1 W/m^2/sr/μm = 0.10 erg/cm^2/s/sr/A -# # 1 sr = 1 rad^2 = 4.25452e10 arcsec^2 -# # 1 J/s = 1 W -# # 1 J = 10^7 erg - -# # convert erg/s/cm^2/A/arcsec^2 to erg/s/cm^2/A/sr -# flux1 = flux / (1/4.25452e10) -# # convert erg/s/cm^2/A/sr to W/m^2/sr/um -# flux2 = flux1 * 10 - -# # 对接收面积积分，输出单位 W/m^2/nm -# D = 2 # meter -# f = 28 # meter -# flux3 = flux2 * np.pi * D**2 / 4 / f**2 / 10**3 - -# # # # u_lambda: W/m^2/nm -# # # obj0: -# # V = 73 * 1e-8 # W/m^2/nm -# # Es = V * np.pi * D**2 / 4 / f**2 / 10**3 -# # c1 = 3.7418e-16 -# # c2 = 1.44e-2 -# # t = 5700 -# # # wave需要代入meter. -# # wave0 = np.arange(380, 780) # nm -# # delta_lamba0 = 1 # nm -# # u_lambda = c1 / (wave0*1e-9)**5 / (np.exp(c2 / (wave0*1e-9) / t) - 1) -# # f_lambda = (u_lambda / u_lambda[wave0 == 500]) * Es -# # E0 = np.sum(f_lambda * 1) -# # # plt.plot(wave, u_lambda) -# # # plt.show() - -# # 对波长积分 -# f = interp1d(wave, flux3) -# wave_interp = np.arange(3800, 7800) -# flux3_interp = f(wave_interp) -# # 输出单位 W/m^2 -# delta_lamba = 0.1 # nm -# E = np.sum(flux3_interp * delta_lamba) - -# # pdb.set_trace() - -# return E - ################################################################ def ill2flux(E,path): @@ -355,6 +281,8 @@ class StrayLight(object): def caculateStarLightFilter(self, filter='i'): + + filterIndex = {'nuv':0,'u':1,'g':2,'r':3, 'i':4, 'z':5,'y':6} sat = (ctypes.c_double*3)() sat[:] = self.sat @@ -376,6 +304,8 @@ class StrayLight(object): return max(band_star_e1, band_star_e2) def caculateEarthShineFilter(self, filter='i'): + + filterIndex = {'nuv':0,'u':1,'g':2,'r':3, 'i':4, 'z':5,'y':6} sat = (ctypes.c_double*3)() sat[:] = self.sat ob = (ctypes.c_double*3)() @@ -1893,9 +1823,6 @@ class MCIsimulator(): # EarthShine from straylight sl = StrayLight(self.information['dir_path'],jtime=time_jd, sat=np.array([x_sat, y_sat, z_sat]), radec=np.array([(ra*u.degree).value, (dec*u.degree).value])) - - - earth_e = sl.caculateEarthShineFilter(filter='r') star_e = sl.caculateStarLightFilter(filter='r') angle_earth = earth_angle(time_jd, x_sat, y_sat, z_sat, ra, dec) @@ -1905,8 +1832,6 @@ class MCIsimulator(): earthshine_wave0, earthshine_flux0 = ill2flux(earth_e+star_e, self.information['dir_path']) - - # sample as mci wavelength wave_mci = np.linspace(2500, 11000, 8501) #np.arange(2500, 11000, 1) f1 = interp1d(wave0, zodi0) diff --git a/csst_mci_sim/earthshine.py b/csst_mci_sim/earthshine.py new file mode 100644 index 0000000000000000000000000000000000000000..35f9cab0d53cada7acd13654242796dea5f5c229 --- /dev/null +++ b/csst_mci_sim/earthshine.py @@ -0,0 +1,216 @@ + +import numpy as np +import julian +from datetime import datetime +from astropy.time import Time +from astropy.coordinates import get_sun +from astropy.coordinates import SkyCoord +import astropy.coordinates as coord +import pandas as pd +from astropy import units as u + +from scipy.interpolate import interp1d + +from scipy import interpolate + +import ctypes + + +def transRaDec2D(ra, dec): + # radec转为竞天程序里的ob, 赤道坐标系下的笛卡尔三维坐标xyz. + x1 = np.cos(dec / 57.2957795) * np.cos(ra / 57.2957795) + y1 = np.cos(dec / 57.2957795) * np.sin(ra / 57.2957795) + z1 = np.sin(dec / 57.2957795) + return np.array([x1, y1, z1]) + + +def earth_angle(time_jd, x_sat, y_sat, z_sat, ra_obj, dec_obj): + + ra_sat = np.arctan2(y_sat, x_sat) / np.pi * 180 + dec_sat = np.arctan2(z_sat, np.sqrt(x_sat**2+y_sat**2)) / np.pi * 180 + radec_sat = SkyCoord(ra=ra_sat*u.degree, + dec=dec_sat*u.degree, frame='gcrs') + lb_sat = radec_sat.transform_to('geocentrictrueecliptic') + + # get the obj location + radec_obj = SkyCoord(ra=ra_obj*u.degree, + dec=dec_obj*u.degree, frame='gcrs') + lb_obj = radec_obj.transform_to('geocentrictrueecliptic') + + # calculate the angle between sub-satellite point and the earth side + earth_radius = 6371 # km + sat_height = np.sqrt(x_sat**2 + y_sat**2 + z_sat**2) + angle_a = np.arcsin(earth_radius/sat_height) / np.pi * 180 + + # calculate the angle between satellite position and the target position + angle_b = lb_sat.separation(lb_obj) + + # calculat the earth angle + angle = 180 - angle_a - angle_b.degree + + return angle + +############################################################################### + + +def ill2flux(E, path): + + # use template from sky_bkg (background_spec_hst.dat) + filename = path+'MCI_inputData/refs/background_spec_hst.dat' + cat_spec = pd.read_csv(filename, sep='\s+', header=None, comment='#') + wave0 = cat_spec[0].values # A + spec0 = cat_spec[2].values # erg/s/cm^2/A/arcsec^2 + + # convert erg/s/cm^2/A/arcsec^2 to erg/s/cm^2/A/sr + flux1 = spec0 / (1/4.25452e10) + # convert erg/s/cm^2/A/sr to W/m^2/sr/um + flux2 = flux1 * 10 + + # 对接收面积积分，输出单位 W/m^2/nm + D = 2 # meter + f = 28 # meter, 焦距，转换关系来源于王维notes. + flux3 = flux2 * np.pi * D**2 / 4 / f**2 / 10**3 + + f = interp1d(wave0, flux3) + wave_range = np.arange(3800, 7800) + flux3_mean = f(wave_range) + delta_lamba = 0.1 # nm + E0 = np.sum(flux3_mean * delta_lamba) + + factor = E / E0 + spec_scaled = factor * spec0 + + return wave0, spec_scaled + +############################################################## + + +def earthshine(path, time_jd, x_sat, y_sat, z_sat, ra, dec): + # EarthShine from straylight + sl = StrayLight(path, jtime=time_jd, sat=np.array( + [x_sat, y_sat, z_sat]), radec=np.array([(ra*u.degree).value, (dec*u.degree).value])) + + earth_e = sl.caculateEarthShineFilter(filter='r') + + angle_earth = earth_angle(time_jd, x_sat, y_sat, z_sat, ra, dec) + + if angle_earth < 0: + earth_e = 0 + + earthshine_wave0, earthshine_flux0 = ill2flux(earth_e, path) + + # sample as mci wavelength + wave_mci = np.linspace(2500, 11000, 8501) # np.arange(2500, 11000, 1) + + f2 = interp1d(earthshine_wave0, earthshine_flux0) + earthshine = f2(wave_mci) + + return earthshine +################################################################# + + +class StrayLight(object): + + def __init__(self, path, jtime=2460843., sat=np.array([0, 0, 0]), radec=np.array([0, 0])): + + self.path = path + self.jtime = jtime + self.sat = sat + self.equator = coord.SkyCoord( + radec[0]*u.degree, radec[1]*u.degree, frame='icrs') + self.ecliptic = self.equator.transform_to('barycentrictrueecliptic') + self.pointing = transRaDec2D(radec[0], radec[1]) + self.slcdll = ctypes.CDLL( + self.path+'MCI_inputData/refs/libstraylight.so') # dylib + + self.slcdll.Calculate.argtypes = [ctypes.c_double, ctypes.POINTER(ctypes.c_double), + ctypes.POINTER(ctypes.c_double), ctypes.POINTER( + ctypes.c_double), + ctypes.POINTER(ctypes.c_double), ctypes.c_char_p] + + self.slcdll.PointSource.argtypes = [ctypes.c_double, ctypes.POINTER(ctypes.c_double), + ctypes.POINTER(ctypes.c_double), ctypes.POINTER( + ctypes.c_double), + ctypes.POINTER(ctypes.c_double), ctypes.c_char_p] + + self.slcdll.EarthShine.argtypes = [ctypes.c_double, ctypes.POINTER(ctypes.c_double), + ctypes.POINTER(ctypes.c_double), ctypes.POINTER( + ctypes.c_double), + ctypes.POINTER(ctypes.c_double)] + + self.slcdll.Zodiacal.argtypes = [ctypes.c_double, ctypes.POINTER(ctypes.c_double), + ctypes.POINTER(ctypes.c_double)] + self.slcdll.ComposeY.argtypes = [ctypes.POINTER(ctypes.c_double), ctypes.POINTER(ctypes.c_double), + ctypes.POINTER(ctypes.c_double)] + self.slcdll.Init.argtypes = [ + ctypes.c_char_p, ctypes.c_char_p, ctypes.c_char_p, ctypes.c_char_p] + self.deFn = self.path+"MCI_inputData/refs/DE405" + self.PSTFn = self.path+"MCI_inputData/refs/PST" + self.RFn = self.path+"MCI_inputData/refs/R" + self.ZolFn = self.path+"MCI_inputData/refs/Zodiacal" + self.brightStarTabFn = self.path+"MCI_inputData/refs/BrightGaia_with_csst_mag" + + self.slcdll.Init(str.encode(self.deFn), str.encode( + self.PSTFn), str.encode(self.RFn), str.encode(self.ZolFn)) + + def caculateStarLightFilter(self, filter='i'): + filterIndex = {'nuv': 0, 'u': 1, 'g': 2, + 'r': 3, 'i': 4, 'z': 5, 'y': 6} + + sat = (ctypes.c_double*3)() + sat[:] = self.sat + ob = (ctypes.c_double*3)() + ob[:] = self.pointing + + py1 = (ctypes.c_double*3)() + py2 = (ctypes.c_double*3)() + self.slcdll.ComposeY(ob, py1, py2) + + star_e1 = (ctypes.c_double*7)() + self.slcdll.PointSource(self.jtime, sat, ob, py1, + star_e1, str.encode(self.brightStarTabFn)) + star_e2 = (ctypes.c_double*7)() + self.slcdll.PointSource(self.jtime, sat, ob, py2, + star_e2, str.encode(self.brightStarTabFn)) + + band_star_e1 = star_e1[:][filterIndex[filter]] + band_star_e2 = star_e2[:][filterIndex[filter]] + + return max(band_star_e1, band_star_e2) + + def caculateEarthShineFilter(self, filter='i'): + + filterIndex = {'nuv': 0, 'u': 1, 'g': 2, + 'r': 3, 'i': 4, 'z': 5, 'y': 6} + sat = (ctypes.c_double*3)() + sat[:] = self.sat + ob = (ctypes.c_double*3)() + ob[:] = self.pointing + + py1 = (ctypes.c_double*3)() + py2 = (ctypes.c_double*3)() + self.slcdll.ComposeY(ob, py1, py2) + + earth_e1 = (ctypes.c_double*7)() + self.slcdll.EarthShine(self.jtime, sat, ob, py1, + earth_e1) # e[7]代表7个波段的照度 + + earth_e2 = (ctypes.c_double*7)() + self.slcdll.EarthShine(self.jtime, sat, ob, py2, earth_e2) + + band_earth_e1 = earth_e1[:][filterIndex[filter]] + band_earth_e2 = earth_e2[:][filterIndex[filter]] + + return max(band_earth_e1, band_earth_e2) + + +############################################################################### +### test +# path='/home/yan/MCI/' +# time_jd = 2460417.59979167 +# x_sat = -4722.543136 +# y_sat = -1478.219213 +# z_sat = 4595.402769 +# ra = 116.18081536720157 +# dec= 39.42316681016602 +# earthshine0=earthshine(path,time_jd, x_sat, y_sat, z_sat, ra, dec) diff --git a/csst_mci_sim/shraylight.py b/csst_mci_sim/shraylight.py new file mode 100644 index 0000000000000000000000000000000000000000..d362715e24f429eec70764fe775d1b883b65ecb0 --- /dev/null +++ b/csst_mci_sim/shraylight.py @@ -0,0 +1,216 @@ + +import numpy as np +import julian +from datetime import datetime +from astropy.time import Time +from astropy.coordinates import get_sun +from astropy.coordinates import SkyCoord +import astropy.coordinates as coord +import pandas as pd +from astropy import units as u + +from scipy.interpolate import interp1d + +from scipy import interpolate + +import ctypes + + +def transRaDec2D(ra, dec): + # radec转为竞天程序里的ob, 赤道坐标系下的笛卡尔三维坐标xyz. + x1 = np.cos(dec / 57.2957795) * np.cos(ra / 57.2957795) + y1 = np.cos(dec / 57.2957795) * np.sin(ra / 57.2957795) + z1 = np.sin(dec / 57.2957795) + return np.array([x1, y1, z1]) + + +def earth_angle(time_jd, x_sat, y_sat, z_sat, ra_obj, dec_obj): + + ra_sat = np.arctan2(y_sat, x_sat) / np.pi * 180 + dec_sat = np.arctan2(z_sat, np.sqrt(x_sat**2+y_sat**2)) / np.pi * 180 + radec_sat = SkyCoord(ra=ra_sat*u.degree, + dec=dec_sat*u.degree, frame='gcrs') + lb_sat = radec_sat.transform_to('geocentrictrueecliptic') + + # get the obj location + radec_obj = SkyCoord(ra=ra_obj*u.degree, + dec=dec_obj*u.degree, frame='gcrs') + lb_obj = radec_obj.transform_to('geocentrictrueecliptic') + + # calculate the angle between sub-satellite point and the earth side + earth_radius = 6371 # km + sat_height = np.sqrt(x_sat**2 + y_sat**2 + z_sat**2) + angle_a = np.arcsin(earth_radius/sat_height) / np.pi * 180 + + # calculate the angle between satellite position and the target position + angle_b = lb_sat.separation(lb_obj) + + # calculat the earth angle + angle = 180 - angle_a - angle_b.degree + + return angle + +############################################################################### + + +def ill2flux(E, path): + + # use template from sky_bkg (background_spec_hst.dat) + filename = path+'MCI_inputData/refs/background_spec_hst.dat' + cat_spec = pd.read_csv(filename, sep='\s+', header=None, comment='#') + wave0 = cat_spec[0].values # A + spec0 = cat_spec[2].values # erg/s/cm^2/A/arcsec^2 + + # convert erg/s/cm^2/A/arcsec^2 to erg/s/cm^2/A/sr + flux1 = spec0 / (1/4.25452e10) + # convert erg/s/cm^2/A/sr to W/m^2/sr/um + flux2 = flux1 * 10 + + # 对接收面积积分，输出单位 W/m^2/nm + D = 2 # meter + f = 28 # meter, 焦距，转换关系来源于王维notes. + flux3 = flux2 * np.pi * D**2 / 4 / f**2 / 10**3 + + f = interp1d(wave0, flux3) + wave_range = np.arange(3800, 7800) + flux3_mean = f(wave_range) + delta_lamba = 0.1 # nm + E0 = np.sum(flux3_mean * delta_lamba) + + factor = E / E0 + spec_scaled = factor * spec0 + + return wave0, spec_scaled + +############################################################## + + +def stray_light(path, time_jd, x_sat, y_sat, z_sat, ra, dec): + # EarthShine from straylight + sl = StrayLight(path, jtime=time_jd, sat=np.array( + [x_sat, y_sat, z_sat]), radec=np.array([(ra*u.degree).value, (dec*u.degree).value])) + + star_e = sl.caculateStarLightFilter(filter='r') + + # angle_earth = earth_angle(time_jd, x_sat, y_sat, z_sat, ra, dec) + + # if angle_earth < 0: + # earth_e = 0 + + earthshine_wave0, earthshine_flux0 = ill2flux(star_e, path) + + # sample as mci wavelength + wave_mci = np.linspace(2500, 11000, 8501) # np.arange(2500, 11000, 1) + + f2 = interp1d(earthshine_wave0, earthshine_flux0) + star_stray = f2(wave_mci) + + return star_stray +################################################################# + + +class StrayLight(object): + + def __init__(self, path, jtime=2460843., sat=np.array([0, 0, 0]), radec=np.array([0, 0])): + + self.path = path + self.jtime = jtime + self.sat = sat + self.equator = coord.SkyCoord( + radec[0]*u.degree, radec[1]*u.degree, frame='icrs') + self.ecliptic = self.equator.transform_to('barycentrictrueecliptic') + self.pointing = transRaDec2D(radec[0], radec[1]) + self.slcdll = ctypes.CDLL( + self.path+'MCI_inputData/refs/libstraylight.so') # dylib + + self.slcdll.Calculate.argtypes = [ctypes.c_double, ctypes.POINTER(ctypes.c_double), + ctypes.POINTER(ctypes.c_double), ctypes.POINTER( + ctypes.c_double), + ctypes.POINTER(ctypes.c_double), ctypes.c_char_p] + + self.slcdll.PointSource.argtypes = [ctypes.c_double, ctypes.POINTER(ctypes.c_double), + ctypes.POINTER(ctypes.c_double), ctypes.POINTER( + ctypes.c_double), + ctypes.POINTER(ctypes.c_double), ctypes.c_char_p] + + self.slcdll.EarthShine.argtypes = [ctypes.c_double, ctypes.POINTER(ctypes.c_double), + ctypes.POINTER(ctypes.c_double), ctypes.POINTER( + ctypes.c_double), + ctypes.POINTER(ctypes.c_double)] + + self.slcdll.Zodiacal.argtypes = [ctypes.c_double, ctypes.POINTER(ctypes.c_double), + ctypes.POINTER(ctypes.c_double)] + self.slcdll.ComposeY.argtypes = [ctypes.POINTER(ctypes.c_double), ctypes.POINTER(ctypes.c_double), + ctypes.POINTER(ctypes.c_double)] + self.slcdll.Init.argtypes = [ + ctypes.c_char_p, ctypes.c_char_p, ctypes.c_char_p, ctypes.c_char_p] + self.deFn = self.path+"MCI_inputData/refs/DE405" + self.PSTFn = self.path+"MCI_inputData/refs/PST" + self.RFn = self.path+"MCI_inputData/refs/R" + self.ZolFn = self.path+"MCI_inputData/refs/Zodiacal" + self.brightStarTabFn = self.path+"MCI_inputData/refs/BrightGaia_with_csst_mag" + + self.slcdll.Init(str.encode(self.deFn), str.encode( + self.PSTFn), str.encode(self.RFn), str.encode(self.ZolFn)) + + def caculateStarLightFilter(self, filter='i'): + filterIndex = {'nuv': 0, 'u': 1, 'g': 2, + 'r': 3, 'i': 4, 'z': 5, 'y': 6} + + sat = (ctypes.c_double*3)() + sat[:] = self.sat + ob = (ctypes.c_double*3)() + ob[:] = self.pointing + + py1 = (ctypes.c_double*3)() + py2 = (ctypes.c_double*3)() + self.slcdll.ComposeY(ob, py1, py2) + + star_e1 = (ctypes.c_double*7)() + self.slcdll.PointSource(self.jtime, sat, ob, py1, + star_e1, str.encode(self.brightStarTabFn)) + star_e2 = (ctypes.c_double*7)() + self.slcdll.PointSource(self.jtime, sat, ob, py2, + star_e2, str.encode(self.brightStarTabFn)) + + band_star_e1 = star_e1[:][filterIndex[filter]] + band_star_e2 = star_e2[:][filterIndex[filter]] + + return max(band_star_e1, band_star_e2) + + def caculateEarthShineFilter(self, filter='i'): + + filterIndex = {'nuv': 0, 'u': 1, 'g': 2, + 'r': 3, 'i': 4, 'z': 5, 'y': 6} + sat = (ctypes.c_double*3)() + sat[:] = self.sat + ob = (ctypes.c_double*3)() + ob[:] = self.pointing + + py1 = (ctypes.c_double*3)() + py2 = (ctypes.c_double*3)() + self.slcdll.ComposeY(ob, py1, py2) + + earth_e1 = (ctypes.c_double*7)() + self.slcdll.EarthShine(self.jtime, sat, ob, py1, + earth_e1) # e[7]代表7个波段的照度 + + earth_e2 = (ctypes.c_double*7)() + self.slcdll.EarthShine(self.jtime, sat, ob, py2, earth_e2) + + band_earth_e1 = earth_e1[:][filterIndex[filter]] + band_earth_e2 = earth_e2[:][filterIndex[filter]] + + return max(band_earth_e1, band_earth_e2) + + +############################################################################### +## test +# path='/home/yan/MCI/' +# time_jd = 2460417.59979167 +# x_sat = -4722.543136 +# y_sat = -1478.219213 +# z_sat = 4595.402769 +# ra = 116.18081536720157 +# dec= 39.42316681016602 +# earthshine0=stray_light(path,time_jd, x_sat, y_sat, z_sat, ra, dec) diff --git a/csst_mci_sim/support/__pycache__/__init__.cpython-311.pyc b/csst_mci_sim/support/__pycache__/__init__.cpython-311.pyc new file mode 100644 index 0000000000000000000000000000000000000000..f365a255b2a496a8b7fd9fd6d637b20b85e235de Binary files /dev/null and b/csst_mci_sim/support/__pycache__/__init__.cpython-311.pyc differ diff --git a/csst_mci_sim/support/__pycache__/cosmicrays.cpython-311.pyc b/csst_mci_sim/support/__pycache__/cosmicrays.cpython-311.pyc new file mode 100644 index 0000000000000000000000000000000000000000..640fd538fa04cf3973401899b8f82dcbc53e512e Binary files /dev/null and b/csst_mci_sim/support/__pycache__/cosmicrays.cpython-311.pyc differ diff --git a/csst_mci_sim/support/__pycache__/logger.cpython-311.pyc b/csst_mci_sim/support/__pycache__/logger.cpython-311.pyc new file mode 100644 index 0000000000000000000000000000000000000000..667ec88207c5cbff63444d0e2650a66f96b87462 Binary files /dev/null and b/csst_mci_sim/support/__pycache__/logger.cpython-311.pyc differ diff --git a/csst_mci_sim/zodiacal.py b/csst_mci_sim/zodiacal.py new file mode 100644 index 0000000000000000000000000000000000000000..350340c951b0ae5b5bff971fb3fc481dfad7eaa2 --- /dev/null +++ b/csst_mci_sim/zodiacal.py @@ -0,0 +1,75 @@ +import numpy as np +import julian +from datetime import datetime +from astropy.time import Time +from astropy.coordinates import get_sun +from astropy.coordinates import SkyCoord +import pandas as pd +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. + + :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. + :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 + + """ + + # 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') + 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(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]) + 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='#') + 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^�? + + spec = spec0 * (zodi_obj / zodi_norm) * 1e-8 # W m^�? sr^�? μm^�? + + # 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 + + return wave_A, spec_erg2 + + +############################### +# path='/home/yan/MCI/' +# wave0, zodi0=zodiacal(10.0, 20.0, '2024-04-04', path) diff --git a/tests/test_earthshine.py b/tests/test_earthshine.py new file mode 100644 index 0000000000000000000000000000000000000000..8a5705e4a904364e4c01f42ffa661ed199840482 --- /dev/null +++ b/tests/test_earthshine.py @@ -0,0 +1,58 @@ +""" +Identifier: mci_sim/tests/test_earthshine.py +Name: test_earthshine.py +Description: Test earthshine sim. +Author: Zhaojun Yan +Created: 2024-04-09 +Modified-History: + 2024-04-09, Zhaojun Yan, created + +""" + +import unittest +import os +import sys +import faulthandler +from csst_mci_sim import earthshine as eth + + +class TestDemoFunction(unittest.TestCase): + def test_earthshine(self): + """ + Aim + --- + Test earthshine sim function: . + + 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/') + print(dir_path) + + # 获取当前工作目录 + # current_path = os.getcwd() + # print("当前路径:", current_path) + + time_jd = 2460417.59979167 + x_sat = -4722.543136 + y_sat = -1478.219213 + z_sat = 4595.402769 + ra = 116.18081536720157 + dec = 39.42316681016602 + earthshine0 = eth.earthshine( + dir_path, time_jd, x_sat, y_sat, z_sat, ra, dec) + + self.assertEqual( + 1, 1, + "case 1: EXDF sim passes.", + ) + ############################################# diff --git a/tests/test_straylight.py b/tests/test_straylight.py new file mode 100644 index 0000000000000000000000000000000000000000..ce98abcab3b87a1089f7f50893ab257c361a1cce --- /dev/null +++ b/tests/test_straylight.py @@ -0,0 +1,58 @@ +""" +Identifier: mci_sim/tests/test_straylight.py +Name: test_straylight.py +Description: Test straylight sim. +Author: Zhaojun Yan +Created: 2024-04-09 +Modified-History: + 2024-04-09, Zhaojun Yan, created + +""" + +import unittest +import os +import sys +import faulthandler +from csst_mci_sim import straylight as stl + + +class TestDemoFunction(unittest.TestCase): + def test_straylight(self): + """ + Aim + --- + Test straylight sim function: . + + 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/') + print(dir_path) + + # 获取当前工作目录 + # current_path = os.getcwd() + # print("当前路径:", current_path) + + time_jd = 2460417.59979167 + x_sat = -4722.543136 + y_sat = -1478.219213 + z_sat = 4595.402769 + ra = 116.18081536720157 + dec = 39.42316681016602 + straylight0 = stl.stray_light( + dir_path, time_jd, x_sat, y_sat, z_sat, ra, dec) + + self.assertEqual( + 1, 1, + "case 1: EXDF sim passes.", + ) + ############################################# diff --git a/tests/test_zodiacal.py b/tests/test_zodiacal.py new file mode 100644 index 0000000000000000000000000000000000000000..78efdf83985bd5a73a4e640b0e523960e9a58cc7 --- /dev/null +++ b/tests/test_zodiacal.py @@ -0,0 +1,51 @@ +""" +Identifier: mci_sim/tests/test_zodiacal.py +Name: test_zodiacal.py +Description: Test zodiacal sim. +Author: Zhaojun Yan +Created: 2024-04-09 +Modified-History: + 2024-04-09, Zhaojun Yan, created + +""" + +import unittest +import os +import sys +import faulthandler +from csst_mci_sim import zodiacal + + +class TestDemoFunction(unittest.TestCase): + def test_zodiacal(self): + """ + Aim + --- + Test zodiacal sim function: . + + 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/') + print(dir_path) + + # 获取当前工作目录 + # current_path = os.getcwd() + # print("当前路径:", current_path) + + zodiacal(10.0, 20.0, '2024-04-04', dir_path): + + self.assertEqual( + 1, 1, + "case 1: EXDF sim passes.", + ) + #############################################