Loading csst_mci_sim/straylight.pydeleted 100644 → 0 +0 −216 Original line number Diff line number Diff line 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) Loading
csst_mci_sim/straylight.pydeleted 100644 → 0 +0 −216 Original line number Diff line number Diff line 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)
