import ctypes import numpy as np import astropy.constants as cons from scipy import interpolate import math from astropy.table import Table import astropy.coordinates as coord from astropy import units as u import sys try: import importlib.resources as pkg_resources except ImportError: # Try backported to PY<37 'importlib_resources' import importlib_resources as pkg_resources 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/zhangxin/Work/SlitlessSim/CSST_SIM/CSST_C6/straylight/straylight/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): 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 getAngle132(x1 = 0, y1 = 0, z1 = 0, x2 = 0, y2 = 0, z2 = 0, x3 = 0, y3 = 0, z3 = 0): cosValue = 0; angle = 0; x11 = x1-x3; y11 = y1-y3; z11 = z1-z3; x22 = x2-x3; y22 = y2-y3; z22 = z2-z3; tt = np.sqrt((x11*x11 + y11*y11 + z11* z11) * (x22*x22 + y22*y22 + z22*z22)); if(tt==0): return 0; cosValue = (x11*x22+y11*y22+z11*z22)/tt; if (cosValue > 1): cosValue = 1; if (cosValue < -1): cosValue = -1; angle = math.acos(cosValue); return angle * 360 / (2 * math.pi); def calculateAnglePwithEarth(sat = np.array([0,0,0]), pointing = np.array([0,0,0]), sun = np.array([0,0,0])): modSat = np.sqrt(sat[0]*sat[0] + sat[1]*sat[1]+sat[2]*sat[2]) modPoint = np.sqrt(pointing[0]*pointing[0] + pointing[1]*pointing[1] + pointing[2]*pointing[2]) withLocalZenithAngle = (pointing[0] * sat[0] + pointing[1] * sat[1] + pointing[2] * sat[2]) / (modPoint*modSat) innerM_sat_sun = sat[0] * sun[0] + sat[1] * sun[1] + sat[2] * sun[2] cosAngle = innerM_sat_sun / (modSat * cons.au.value/1000) isInSunSide = 1 if (cosAngle < -0.3385737): #cos109.79 isInSunSide = -1; elif cosAngle >= -0.3385737 and cosAngle <= 0.3385737: isInSunSide = 0; return math.acos(withLocalZenithAngle)*180/math.pi,isInSunSide # /** # * *eCoor = ra, *eCoor+1 = dec # */ def Cartesian2Equatorial(carCoor = np.array([0,0,0])): eCoor = np.zeros(2) if (carCoor[0] > 0 and carCoor[1] >= 0): eCoor[0] = math.atan(carCoor[1] / carCoor[0]) * 360 / (2 * math.pi) elif (carCoor[0] < 0): eCoor[0] = (math.atan(carCoor[1] / carCoor[0]) + math.pi) * 360 / (2 * math.pi) elif (carCoor[0] > 0 and carCoor[1] < 0): eCoor[0] = (math.atan(carCoor[1] / carCoor[0]) + 2 * math.pi) * 360 / (2 * math.pi) elif (carCoor[0] == 0 and carCoor[1] < 0): eCoor[0] = 270 elif (carCoor[0] == 0 and carCoor[1] > 0): eCoor[0] = 90 eCoor[1] = math.atan(carCoor[2] / np.sqrt(carCoor[0] * carCoor[0] + carCoor[1] * carCoor[1])) * 360 / (2 * math.pi) return eCoor class Straylight(object): def __init__(self, jtime = 2460843., sat_pos = np.array([0,0,0]), pointing_radec = np.array([0,0]), sun_pos = np.array([0,0,0])): self.jtime = jtime self.sat = sat_pos self.sun_pos = sun_pos self.equator = coord.SkyCoord(pointing_radec[0]*u.degree, pointing_radec[1]*u.degree,frame='icrs') self.ecliptic = self.equator.transform_to('barycentrictrueecliptic') self.pointing = transRaDec2D(pointing_radec[0], pointing_radec[1]) platForm = sys.platform if platForm == 'darwin': try: with pkg_resources.files('ObservationSim.Straylight.lib').joinpath('libstraylight.dylib') as dllFn: self.slcdll = ctypes.CDLL(dllFn) except AttributeError: with pkg_resources.path('ObservationSim.Straylight.lib', 'libstraylight.dylib') as dllFn: self.slcdll = ctypes.CDLL(dllFn) elif platForm == 'linux': try: with pkg_resources.files('ObservationSim.Straylight.lib').joinpath('libstraylight.so') as dllFn: self.slcdll = ctypes.CDLL(dllFn) except AttributeError: with pkg_resources.path('ObservationSim.Straylight.lib', 'libstraylight.so') as dllFn: self.slcdll = ctypes.CDLL(dllFn) # self.slcdll=ctypes.CDLL('./libstraylight.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] try: with pkg_resources.files('ObservationSim.Straylight.lib').joinpath('DE405') as tFn: self.deFn = tFn.as_uri()[7:] except AttributeError: with pkg_resources.path('ObservationSim.Straylight.lib', 'DE405') as tFn: self.deFn = tFn.as_uri()[7:] try: with pkg_resources.files('ObservationSim.Straylight.lib').joinpath('PST') as tFn: self.PSTFn = tFn.as_uri()[7:] except AttributeError: with pkg_resources.path('ObservationSim.Straylight.lib', 'PST') as tFn: self.PSTFn = tFn.as_uri()[7:] try: with pkg_resources.files('ObservationSim.Straylight.lib').joinpath('R') as tFn: self.RFn = tFn.as_uri()[7:] except AttributeError: with pkg_resources.path('ObservationSim.Straylight.lib', 'R') as tFn: self.RFn = tFn.as_uri()[7:] try: with pkg_resources.files('ObservationSim.Straylight.lib').joinpath('Zodiacal') as tFn: self.ZolFn = tFn.as_uri()[7:] except AttributeError: with pkg_resources.path('ObservationSim.Straylight.lib', 'Zodiacal') as tFn: self.ZolFn = tFn.as_uri()[7:] try: with pkg_resources.files('ObservationSim.Straylight.lib').joinpath('BrightGaia_with_csst_mag') as tFn: self.brightStarTabFn = tFn.as_uri()[7:] except AttributeError: with pkg_resources.path('ObservationSim.Straylight.lib', 'BrightGaia_with_csst_mag') as tFn: self.brightStarTabFn = tFn.as_uri()[7:] print(self.deFn) self.slcdll.Init(str.encode(self.deFn), str.encode(self.PSTFn), str.encode(self.RFn), str.encode(self.ZolFn)) def getFilterAndCCD_Q(self, filter = 'i'): try: with pkg_resources.files('ObservationSim.Instrument.data.ccd').joinpath(filterCCD[filter] + '.txt') as ccd_fn: q_ccd_f = np.loadtxt(ccd_fn) except AttributeError: with pkg_resources.path('ObservationSim.Instrument.data.ccd', filterCCD[filter] + '.txt') as ccd_fn: q_ccd_f = np.loadtxt(ccd_fn) try: with pkg_resources.files('ObservationSim.Instrument.data.filters').joinpath(filter + '.txt') as filter_fn: q_fil_f = np.loadtxt(filter_fn) except AttributeError: with pkg_resources.path('ObservationSim.Instrument.data.filters', filter + '.txt') as filter_fn: q_fil_f = np.loadtxt(filter_fn) band_s = 2000 band_e = 11000 q_ccd_f[:,0] = q_ccd_f[:,0]*10 q_ccd = np.zeros([q_ccd_f.shape[0]+2,q_ccd_f.shape[1]]) q_ccd[1:-1,:] = q_ccd_f q_ccd[0] = [band_s,0] q_ccd[-1] = [band_e,0] q_fil = np.zeros([q_fil_f.shape[0]+2,q_fil_f.shape[1]]) q_fil[1:-1,:] = q_fil_f q_fil[0] = [band_s,0] q_fil[-1] = [band_e,0] q_fil_i = interpolate.interp1d(q_fil[:,0], q_fil[:,1]) q_ccd_i = interpolate.interp1d(q_ccd[:,0], q_ccd[:,1]) bands = np.arange(bandRange[filter][0], bandRange[filter][1],0.5) q_ccd_fil = q_fil_i(bands)*q_ccd_i(bands) return np.trapz(q_ccd_fil, bands)/(bandRange[filter][1]-bandRange[filter][0]) def calculateEarthShineFilter(self, filter = 'i', pixel_size_phy = 10 ): 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) 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]] q=self.getFilterAndCCD_Q(filter=filter) p_lambda = filterPivotWave[filter] c = cons.c.value h = cons.h.value pix_earth_e1 = band_earth_e1/(h*c/(p_lambda*1e-10))*pixel_size_phy*1e-6*pixel_size_phy*1e-6*q pix_earth_e2 = band_earth_e2/(h*c/(p_lambda*1e-10))*pixel_size_phy*1e-6*pixel_size_phy*1e-6*q if pix_earth_e1< pix_earth_e2: return pix_earth_e1, py1[:] else: return pix_earth_e2, py2[:] """ calculate zodiacal call c++ program, seems to have some problem """ def calculateZodiacalFilter1(self, filter = 'i', pixel_size_phy = 10 ): sat = (ctypes.c_double*3)() sat[:] = self.sat ob = (ctypes.c_double*3)() ob[:]=self.pointing zodical_e = (ctypes.c_double*7)() self.slcdll.Zodiacal(self.jtime,ob,zodical_e) ob1 = (ctypes.c_double*2)() ob1[:] = np.array([self.ecliptic.lon.value, self.ecliptic.lat.value]) zodical_e1 = (ctypes.c_double*7)() self.slcdll.Zodiacal1(ob1,zodical_e1) band_zodical_e = zodical_e[:][filterIndex[filter]] q=self.getFilterAndCCD_Q(filter=filter) p_lambda = filterPivotWave[filter] c = cons.c.value h = cons.h.value pix_zodical_e = band_zodical_e/(h*c/(p_lambda*1e-10))*pixel_size_phy*1e-6*pixel_size_phy*1e-6*q return pix_zodical_e, band_zodical_e """ calculate zodiacal use python """ def calculateZodiacalFilter2(self,filter = 'i', aper = 2, pixelsize = 0.074, sun_pos = np.array([0,0,0])): spec, v_mag = self.calculateZodicalSpec(longitude = self.ecliptic.lon.value, latitude = self.ecliptic.lat.value, sun_pos = sun_pos) # spec = self.calculateZodicalSpec(longitude = lon, latitude = lat) try: with pkg_resources.files('ObservationSim.Instrument.data.throughputs').joinpath(filter + '_throughput.txt') as throughputFn: throughput = np.loadtxt(throughputFn) except AttributeError: with pkg_resources.path('ObservationSim.Instrument.data.throughputs', filter + '_throughput.txt') as throughputFn: throughput = np.loadtxt(throughputFn) deltL = 0.5 lamb = np.arange(bandRange[filter][0], bandRange[filter][1], deltL) speci = interpolate.interp1d(spec['WAVELENGTH'], spec['FLUX']) y = speci(lamb) # erg/s/cm2/A --> photo/s/m2/A flux = y * lamb / (cons.h.value * cons.c.value) * 1e-13 throughput_i = interpolate.interp1d(throughput[:, 0], throughput[:, 1]) throughput_ = throughput_i(lamb) sky_pix = np.trapz(flux*throughput_, lamb) * math.pi * aper*aper/4 * pixelsize * pixelsize # sky_pix_e = np.trapz(y, lamb) * math.pi * aper*aper/4 * pixelsize * pixelsize/(10*10*1e-6*1e-6)*1e-7*1e4 return sky_pix, v_mag#, sky_pix_e def calculateStarLightFilter(self, filter = 'i', pointYaxis = np.array([1,1,1]), pixel_size_phy = 10 ): sat = (ctypes.c_double*3)() sat[:] = self.sat ob = (ctypes.c_double*3)() ob[:]=self.pointing py = (ctypes.c_double*3)() py[:] = pointYaxis q=self.getFilterAndCCD_Q(filter=filter) p_lambda = filterPivotWave[filter] c = cons.c.value h = cons.h.value star_e1 = (ctypes.c_double*7)() self.slcdll.PointSource(self.jtime,sat,ob,py,star_e1, str.encode(self.brightStarTabFn)) band_star_e1 = star_e1[:][filterIndex[filter]] pix_star_e1 = band_star_e1/(h*c/(p_lambda*1e-10))*pixel_size_phy*1e-6*pixel_size_phy*1e-6*q return pix_star_e1 def calculateEarthshineGrating(self, grating = 'GU', pixel_size_phy = 10, normFilter = 'g', aper = 2, pixelsize = 0.074): 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) earth_e2 = (ctypes.c_double*7)() self.slcdll.EarthShine(self.jtime,sat,ob,py2,earth_e2) # zodical_e = (ctypes.c_double*7)() # self.slcdll.Zodiacal(self.jtime,ob,zodical_e) band_earth_e1 = earth_e1[:][filterIndex[normFilter]] band_earth_e2 = earth_e2[:][filterIndex[normFilter]] band_earth_e = band_earth_e2 py = py2[:] if band_earth_e1 photo/s/m2/A flux = y * lamb / (cons.h.value * cons.c.value) * 1e-13 throughput_i = interpolate.interp1d(throughput[:, 0], throughput[:, 1]) throughput_ = throughput_i(lamb) sky_pix = np.trapz(flux*throughput_, lamb) * math.pi * aper*aper/4 * pixelsize * pixelsize lamb = np.arange(s, e, deltL) speci = interpolate.interp1d(spec[:, 0], spec[:, 1]) y = speci(lamb) # erg/s/cm2/A --> photo/s/m2/A flux = y * lamb / (cons.h.value * cons.c.value) * 1e-13 return sky_pix, Table(np.array([lamb, flux]).T,names=('WAVELENGTH', 'FLUX')) def calculateZodicalSpec(self,longitude = 50, latitude = 60, sun_pos = np.array([0,0,0])): from scipy.interpolate import interp2d from scipy.interpolate import griddata try: with pkg_resources.files('ObservationSim.Straylight.data').joinpath('Zodiacal_map1.dat') as z_map_fn: ZL = np.loadtxt(z_map_fn) except AttributeError: with pkg_resources.path('ObservationSim.Straylight.data', 'Zodiacal_map1.dat') as z_map_fn: ZL = np.loadtxt(z_map_fn) # zl_sh = ZL.shape # x = np.arange(0,zl_sh[1],1) # y = np.arange(0,zl_sh[0],1) x = ZL[0,1:] y = ZL[1:,0] X,Y = np.meshgrid(x,y) # f_sur = interp2d(X,Y,ZL,kind='linear') sun_radec = Cartesian2Equatorial(sun_pos) sun_eclip = coord.SkyCoord(sun_radec[0]*u.degree, sun_radec[1]*u.degree,frame='icrs') sun_equtor = sun_eclip.transform_to('barycentrictrueecliptic') longitude = longitude - (sun_equtor.lon*u.degree).value longitude = np.abs(longitude) # print((sun_equtor.lon*u.degree).value) if (longitude > 180): longitude = 360 - longitude latitude = np.abs(latitude) lo = longitude la = latitude zl = griddata((X.flatten(),Y.flatten()),ZL[1:,1:].flatten(),(la,lo), method='cubic').min() zl = zl*(math.pi*math.pi)/(180*180)/(3600*3600)*1e-4*1e7*1e-8*1e-4 # print(zl , '\n') try: with pkg_resources.files('ObservationSim.Straylight.data.sky').joinpath('zodiacal.dat') as zodical_fn: spec = np.loadtxt(zodical_fn) except AttributeError: with pkg_resources.path('ObservationSim.Straylight.data.sky', 'zodiacal.dat') as zodical_fn: spec = np.loadtxt(zodical_fn) speci = interpolate.interp1d(spec[:, 0], spec[:, 1]) flux5000 = speci(5000) f_ration = zl/flux5000 v_mag = np.log10(f_ration)*(-2.5)+22.1 # print("factor:", v_mag, lo, la) return Table(np.array([spec[:,0], spec[:,1]*f_ration]).T,names=('WAVELENGTH', 'FLUX')), v_mag def calculateStrayLightFilter(self, filter = 'i', pixel_size_phy = 10, pixel_scale = 0.074): e1,py = self.calculateEarthShineFilter(filter = filter, pixel_size_phy = pixel_size_phy) e2, _ = self.calculateZodiacalFilter2(filter = filter, sun_pos=self.sun_pos, pixelsize = pixel_scale) e3 = self.calculateStarLightFilter(filter = filter,pointYaxis = py, pixel_size_phy = pixel_size_phy) return e1+e2+e3 def calculateStrayLightGrating(self, grating = 'GI', pixel_size_phy = 10, normFilter_es = 'g'): e1,py = self.calculateEarthshineGrating(grating = grating, pixel_size_phy = pixel_size_phy, normFilter = normFilter_es) e2 = self.calculateStarLightGrating(grating = grating, pointYaxis = py) spec, _ = self.calculateZodicalSpec(longitude = self.ecliptic.lon.value, latitude = self.ecliptic.lat.value, sun_pos = self.sun_pos) return e1+e2, spec def testZodiacal(lon = 285.04312526255366, lat = 30.): c_eclip = coord.SkyCoord(lon*u.degree, lat*u.degree,frame='barycentrictrueecliptic') c_equtor = c_eclip.transform_to('icrs') sl = Straylight(jtime = 2459767.00354975, sat = np.array([]), radec = np.array([(c_equtor.ra*u.degree).value, (c_equtor.dec*u.degree).value])) e_zol, v_mag = sl.calculateZodiacalFilter2(filter = 'i', sun_pos=np.array([-3.70939436e+07, 1.35334903e+08, 5.86673104e+07])) print(e_zol) # ju=2.4608437604166665e+06 # sat = (ctypes.c_double*3)() # sat[:] = np.array([5873.752, -1642.066, 2896.744]) # ob = (ctypes.c_double*3)() # ob[:]=np.array([0.445256,0.76061,-0.47246]) # sl = StrayLight(jtime = ju, sat = np.array([5873.752, -1642.066, 2896.744]), pointing = np.array([-0.445256,-0.76061,0.47246])) # fn = '/Users/zhangxin/Work/SurveyPlan/point/csst_survey_sim_20211028/E17.5_b17.5_beta_11.6_opt_transtime_1_CMG_1_dp_2_0.25_da_10_Texp_1.5_DEC60_500_0.1_800_1000_+5deg.dat' # # surveylist = np.loadtxt(fn) # sky_pix = np.zeros([surveylist.shape[0],7]) # # # i = 693438 # c_eclip = coord.SkyCoord(surveylist[:,2]*u.degree, surveylist[:,1]*u.degree,frame='barycentrictrueecliptic') # c_equtor = c_eclip.transform_to('icrs') # # # # pointing = transRaDec2D((c_equtor[i].ra*u.degree).value, (c_equtor[i].dec*u.degree).value) # # # print(ju, pointing, surveylist[i,3:9]) # # ju = surveylist[i,0] # # sl = StrayLight(jtime = ju, sat = surveylist[i,3:6], pointing = pointing) # # sl.caculateStrayLightGrating(grating = 'GI', pixel_size_phy = 10, normFilter = 'g') # # for i in np.arange(surveylist.shape[0]): # print(i) # if i > 300: # break # # if i != 300: # # continue # # if i != 693438: # # continue # ju = surveylist[i,0] # pointing = transRaDec2D((c_equtor[i].ra*u.degree).value, (c_equtor[i].dec*u.degree).value) # # print(ju, pointing, surveylist[i,3:9]) # sl = StrayLight(jtime = ju, sat = surveylist[i,3:6], radec = np.array([(c_equtor[i].ra*u.degree).value, (c_equtor[i].dec*u.degree).value])) # # strayl_i,s_zoldical ,s_earth, s_earth1 = sl.caculateStrayLightFilter(filter = 'i') # # print(i,strayl_i,s_zoldical,s_earth, s_earth1) # p_cart= transRaDec2D((c_equtor[i].ra*u.degree).value, (c_equtor[i].dec*u.degree).value) # sky_pix[i,6] = getAngle132(x1 = surveylist[i,6], y1 = surveylist[i,7], z1 = surveylist[i,8], x2 = p_cart[0], y2 = p_cart[1], z2 = p_cart[2], x3 = 0, y3 = 0, z3 = 0) # # earthZenithAngle,isInSunSide = calculateAnglePwithEarth(sat = surveylist[i,3:6], pointing = pointing, sun = surveylist[i,6:9]) # sky_pix[i,4] = earthZenithAngle # sky_pix[i,5] = isInSunSide # # e1_,py = sl.caculateEarthShineFilter(filter = 'i') # # e2, e2_ = sl.calculateZodiacalFilter1(filter = 'i') # e3, v_mag = sl.calculateZodiacalFilter2(filter = 'i', sun_pos=surveylist[i,6:9]) # # e4 = sl.caculateStarLightFilter(filter = 'i',pointYaxis = py) # # e4 = 0 # # e1,py = sl.caculateEarthshineGrating(grating = 'GI', pixel_size_phy = 10, normFilter = 'g') # # # e2 = sl.caculateStarLightGrating(grating = 'GV', pointYaxis = py) # e2 = sl.caculateStarLightGrating(grating = 'GI', pointYaxis = py) # e4 = sl.caculateStarLightFilter(filter = 'i',pointYaxis = py) # # e5=sl.caculateStrayLightFilter(filter = 'i', pixel_size_phy = 10, pixel_scale = 0.074, sun_pos = surveylist[i,6:9]) # e6,_=sl.caculateStrayLightGrating(grating = 'GI', normFilter_es = 'g', sun_pos = surveylist[i,6:9]) # # sky_pix[i,0] = e1 # sky_pix[i,1] = e2 # sky_pix[i,2] = e3 # sky_pix[i,3] = e4 # print(e1+e2,e1_+e3+e4,e5,e6) # # # print(e1,e2,e3,e4)