debug

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):
 
@@ -356,6 +282,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
         ob = (ctypes.c_double*3)()
@@ -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)

csst_mci_sim/earthshine.py

+
+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)

csst_mci_sim/shraylight.py

+
+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)

csst_mci_sim/zodiacal.py

+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)

tests/test_earthshine.py

+"""
+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.",
+        )
+    #############################################

tests/test_straylight.py

+"""
+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.",
+        )
+    #############################################

tests/test_zodiacal.py

+"""
+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.",
+        )
+    #############################################