diff --git a/csst_mci_sim/straylight.py b/csst_mci_sim/straylight.py
new file mode 100644
index 0000000000000000000000000000000000000000..cb37588ccabd777a8b991bb6a5ad11bc0a67992f
--- /dev/null
+++ b/csst_mci_sim/straylight.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)