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 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'] 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 = np.array([0,0,0]), radec = np.array([0,0])): 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('./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)] 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)] 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() def getFilterAndCCD_Q(self, filter = 'i'): ccd_fn = Instrument_dir + 'ccd/' + filterCCD[filter] + '.txt' filter_fn = Instrument_dir + 'filters/' + filter + '.txt' q_ccd_f = np.loadtxt(ccd_fn) q_fil_f = np.loadtxt(filter_fn) band_s = 2000 band_e = 11000 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 caculateEarthShineFilter(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 calculateZaodiacalFilter2(self,filter = 'i', aper = 2, pixelsize = 0.074, sun_pos = np.array([0,0,0])): spec, v_mag = self.calculateZaodicalSpec(longitude = self.ecliptic.lon.value, latitude = self.ecliptic.lat.value, sun_pos = sun_pos) # spec = self.calculateZaodicalSpec(longitude = lon, latitude = lat) throughputFn = Instrument_dir + 'throughputs/' + filter + '_throughput.txt' 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 caculateStarLightFilter(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) 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 caculateStrayLightFilter(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) zodical_e = (ctypes.c_double*7)() self.slcdll.Zodiacal(self.jtime,ob,zodical_e) band_earth_e1 = earth_e1[:][filterIndex[filter]] band_earth_e2 = earth_e2[:][filterIndex[filter]] # band_earth_e1 = 0 # band_earth_e2 = 0 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_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 pix_zodical_e = band_zodical_e/(h*c/(p_lambda*1e-10))*pixel_size_phy*1e-6*pixel_size_phy*1e-6*q # star_e1 = (ctypes.c_double*7)() # self.slcdll.PointSource(self.jtime,sat,ob,py1,star_e1) # # star_e2 = (ctypes.c_double*7)() # # self.slcdll.PointSource(self.jtime,sat,ob,py2,star_e2) # band_star_e1 = star_e1[:][filterIndex[filter]] # # band_star_e2 = star_e2[:][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 pix_star_e1 = 0 # pix_star_e2 = band_star_e2/(h*c/(p_lambda*1e-10))*pixel_size_phy*1e-6*pixel_size_phy*1e-6*q return pix_earth_e1+pix_zodical_e+pix_star_e1, pix_zodical_e, pix_earth_e1, pix_earth_e2 def caculateStrayLightGrating(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 = np.min([band_earth_e1, band_earth_e2]) # band_earth_e1 = 0 # band_earth_e2 = 0 # band_zodical_e = zodical_e[:][filterIndex[normFilter]] q=self.getFilterAndCCD_Q(filter=normFilter) p_lambda = filterPivotWave[normFilter] c = cons.c.value h = cons.h.value pix_earth_e = band_earth_e/(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 # pix_zodical_e = band_zodical_e/(h*c/(p_lambda*1e-10))*pixel_size_phy*1e-6*pixel_size_phy*1e-6*q # pix_earth_e = np.min([pix_earth_e1, pix_earth_e2]) # zodical_v, zodical_spec = self.calculatSkylightBySpec(specType = 'zodical', filter = 'g', aper = 2, pixelsize = 0.074) earthshine_v, earthshine_spec = self.calculatSkylightBySpec(specType = 'earthshine', filter = 'g', aper = aper, pixelsize = pixelsize) lamb_earth = earthshine_spec['WAVELENGTH'] flux_earth = earthshine_spec['FLUX']*pix_earth_e/earthshine_v print(pix_earth_e,earthshine_v) earth_v_grating = 0 for s_order in SpecOrder: thpFn = Instrument_dir + 'sls_conf/' + grating + '.Throughput.' + s_order + 'st.fits' thp_ = Table.read(thpFn) thpFn_i = interpolate.interp1d(thp_['WAVELENGTH'], thp_['SENSITIVITY']) thp = thpFn_i(lamb_earth) beamsEarth = np.trapz(flux_earth*thp,lamb_earth)* math.pi*aper*aper/4 * pixelsize * pixelsize earth_v_grating = earth_v_grating + beamsEarth print(beamsEarth) # print(earthshine_v, pix_earth_e, earth_v_grating) return earth_v_grating def calculatSkylightBySpec(self, specType = 'earthshine', filter = 'g', aper = 2, pixelsize = 0.074, s = 2000, e = 11000): specFn = '' if specType == 'zodical': specFn=Instrument_dir + 'sky/zodiacal.dat' elif specType == 'earthshine': specFn= Instrument_dir + 'sky/earthShine.dat' spec = np.loadtxt(specFn) throughputFn = Instrument_dir + 'throughputs/' + filter + '_throughput.txt' throughput = np.loadtxt(throughputFn) deltL = 0.5 lamb = np.arange(bandRange[filter][0], bandRange[filter][1], 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 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 calculateZaodicalSpec(self,longitude = 50, latitude = 60, sun_pos = np.array([0,0,0])): from scipy.interpolate import interp2d from scipy.interpolate import griddata z_map_fn = Instrument_dir + 'Zodiacal_map1.dat' 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') zodical_fn = Instrument_dir + 'sky/zodiacal.dat' 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 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.calculateZaodiacalFilter2(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.calculateZaodiacalFilter2(filter = 'i', sun_pos=surveylist[i,6:9]) e4 = sl.caculateStarLightFilter(filter = 'i',pointYaxis = py) # e4 = 0 # e4 = sl.caculateStrayLightGrating(grating = 'GI', pixel_size_phy = 10, normFilter = 'g') sky_pix[i,0] = e1 sky_pix[i,1] = e3 sky_pix[i,2] = e4 sky_pix[i,3] = v_mag print(e1,e3,e4) # print(e1,e2,e3,e4)