From 22a98299c5828f8e9925090bdf9fb4742f896fb3 Mon Sep 17 00:00:00 2001 From: Yan Zhaojun Date: Fri, 25 Oct 2024 06:20:20 +0000 Subject: [PATCH] Upload New File --- csst_mci_sim/straylight.py | 216 +++++++++++++++++++++++++++++++++++++ 1 file changed, 216 insertions(+) create mode 100644 csst_mci_sim/straylight.py diff --git a/csst_mci_sim/straylight.py b/csst_mci_sim/straylight.py new file mode 100644 index 0000000..cb37588 --- /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) -- GitLab