# #need add environment parameter UNIT_TEST_DATA_ROOT, link to "testData/" #linx and mac can run as follow, need modify the name of file directory #export UNIT_TEST_DATA_ROOT=/Users/zhangxin/Work/SlitlessSim/CSST_SIM/CSST_develop/csst-simulation/tests/testData # import unittest from ObservationSim.Straylight import Straylight import numpy as np import math import astropy.constants as cons import galsim from astropy.table import Table from scipy import interpolate import matplotlib.pyplot as plt import os hubbleAverZodiacal = {'nuv':0.0035,'u':0.0163,'g':0.1109,'r':0.1471,'i':0.1568,'z':0.0953,'y':0.0283} hubbleAverEarthShine = {'nuv':0.00024,'u':0.0051,'g':0.0506,'r':0.0591,'i':0.0568,'z':0.0315,'y':0.0090} 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 class TestStraylight(unittest.TestCase): def __init__(self,methodName='runTest', filter = 'i', grating = "GI"): super(TestStraylight,self).__init__(methodName) # print(file_name) # fn = os.path.join(os.getenv('UNIT_TEST_DATA_ROOT'), file_name) # self.pointingData = np.loadtxt(os.path.join(fn, 'Straylight_test.dat'), dtype=np.double) self.filePath('csst_fz_gc0') self.filter = filter self.grating = grating def filePath(self, file_name): fn = os.path.join(os.getenv('UNIT_TEST_DATA_ROOT'), file_name) self.pointingData = np.loadtxt(os.path.join(fn, 'Straylight_test.dat'), dtype=np.double) def test_EarthShineFilter(self): d_sh = self.pointingData.shape sl_e_pix = np.zeros([d_sh[0],3],dtype=np.double) for i in np.arange(d_sh[0]): # if i > 50: # continue ju = self.pointingData[i, 5] # pointing = transRaDec2D(self.pointingData[i, 0], self.pointingData[i, 1]) # print(ju, pointing, surveylist[i,3:9]) sl = Straylight(jtime=ju, sat_pos=self.pointingData[i, 6:9], pointing_radec=np.array([self.pointingData[i, 0], self.pointingData[i, 1]]),sun_pos=self.pointingData[i,9:12]) e1, py = sl.calculateEarthShineFilter(filter=self.filter) earthZenithAngle, isInSunSide = calculateAnglePwithEarth(sat=self.pointingData[i, 6:9], pointing= sl.pointing, sun=self.pointingData[i,9:12]) # e2, _ = sl.calculateZodiacalFilter2(filter='i', sun_pos=sl.sun_pos) # e3 = sl.calculateStarLightFilter(filter='i', pointYaxis=py) # e_all = sl.calculateStrayLightFilter(filter='i') # s_pix, spec = sl.calculateStrayLightGrating(grating='GI') sl_e_pix[i,0] = e1 sl_e_pix[i, 1] = earthZenithAngle sl_e_pix[i, 2] = isInSunSide median = np.median(sl_e_pix[:,0]) print(' average Earthshine %s: %e' % (self.filter, median)) self.assertTrue(median-hubbleAverEarthShine[self.filter] < 0.1) plt.figure() ids1 = sl_e_pix[:, 2] == 1 ids2 = sl_e_pix[:, 2] != 1 plt.plot(sl_e_pix[ids1, 0], sl_e_pix[ids1, 1], 'r.') plt.plot(sl_e_pix[ids2, 0], sl_e_pix[ids2, 1], 'b.') plt.legend(['In Sun Side', 'In Earths shadow']) plt.xlabel('straylight-earthshine(e-/pixel/s)') plt.ylabel('Angle with local zenith(degree)') plt.show() def test_ZodiacalFilter(self): d_sh = self.pointingData.shape sl_e_pix = np.zeros([d_sh[0],2],dtype=np.double) for i in np.arange(d_sh[0]): ju = self.pointingData[i, 5] sl = Straylight(jtime=ju, sat_pos=self.pointingData[i, 6:9], pointing_radec=np.array([self.pointingData[i, 0], self.pointingData[i, 1]]),sun_pos=self.pointingData[i,9:12]) e1, _ = sl.calculateZodiacalFilter2(filter=self.filter, sun_pos=sl.sun_pos) sl_e_pix[i,0] = e1 sl_e_pix[i,1] = getAngle132(x1=self.pointingData[i,9], y1=self.pointingData[i,10], z1=self.pointingData[i,11], x2=sl.pointing[0], y2=sl.pointing[1], z2=sl.pointing[2], x3=0, y3=0, z3=0) plt.figure() plt.plot(sl_e_pix[:, 0], sl_e_pix[:, 1], 'r.') plt.xlabel('straylight-zodiacal(e-/pixel/s)') plt.ylabel('Angle between pointing and sun(degree)') plt.show() median = np.median(sl_e_pix[:,0]) print(' average Zodiacal %s: %f' % (self.filter, median)) self.assertTrue(median-hubbleAverZodiacal[self.filter] < 0.1) def test_StarFilter(self): d_sh = self.pointingData.shape sl_e_pix = np.zeros(d_sh[0],dtype=np.double) tnum = 10 for i in np.arange(tnum): # if i > 50: # continue ju = self.pointingData[i, 5] # pointing = transRaDec2D(self.pointingData[i, 0], self.pointingData[i, 1]) # print(ju, pointing, surveylist[i,3:9]) sl = Straylight(jtime=ju, sat_pos=self.pointingData[i, 6:9], pointing_radec=np.array([self.pointingData[i, 0], self.pointingData[i, 1]]),sun_pos=self.pointingData[i,9:12]) e1, py = sl.calculateEarthShineFilter(filter=self.filter) # e2, _ = sl.calculateZodiacalFilter2(filter='i', sun_pos=sl.sun_pos) e3 = sl.calculateStarLightFilter(filter=self.filter, pointYaxis=py) # e_all = sl.calculateStrayLightFilter(filter='i') # s_pix, spec = sl.calculateStrayLightGrating(grating='GI') sl_e_pix[i] = e3 median = np.median(sl_e_pix[0:tnum]) print(' average Earthshine %s: %e' % (self.filter, median)) self.assertTrue(median-hubbleAverEarthShine[self.filter] < 0.2) def test_GratingStraylight(self): d_sh = self.pointingData.shape sl_e_pix = np.zeros(d_sh[0],dtype=np.double) tnum = 10 for i in np.arange(tnum): # if i > 50: # continue ju = self.pointingData[i, 5] # pointing = transRaDec2D(self.pointingData[i, 0], self.pointingData[i, 1]) # print(ju, pointing, surveylist[i,3:9]) sl = Straylight(jtime=ju, sat_pos=self.pointingData[i, 6:9], pointing_radec=np.array([self.pointingData[i, 0], self.pointingData[i, 1]]),sun_pos=self.pointingData[i,9:12]) # e1, py = sl.calculateEarthShineFilter(filter=self.filter) # e2, _ = sl.calculateZodiacalFilter2(filter='i', sun_pos=sl.sun_pos) # e3 = sl.calculateStarLightFilter(filter=self.filter, pointYaxis=py) # e_all = sl.calculateStrayLightFilter(filter='i') s_pix, spec = sl.calculateStrayLightGrating(grating=self.grating) sl_e_pix[i] = s_pix plt.figure() plt.plot(spec['WAVELENGTH'], spec['FLUX'], 'r') plt.xlabel('WAVELENGTH') plt.ylabel('F$\lambda$(erg/s/cm2/A/arcsec2)') plt.xlim(2000,10000) plt.show() median = np.median(sl_e_pix[0:tnum]) print(' average Earthshine %s: %e' % (self.grating, median)) self.assertTrue(median < 0.8) if __name__ == '__main__': os.environ['UNIT_TEST_DATA_ROOT']="/Users/zhangxin/Work/SlitlessSim/CSST_SIM/CSST_develop/csst-simulation/tests/testData" suit = unittest.TestSuite() case1 = TestStraylight('test_EarthShineFilter', filter = 'i') suit.addTest(case1) case2 = TestStraylight('test_ZodiacalFilter',filter = 'i') suit.addTest(case2) case3 = TestStraylight('test_StarFilter', filter='i') suit.addTest(case3) case4 = TestStraylight('test_GratingStraylight', grating = 'GI') suit.addTest(case4) unittest.TextTestRunner(verbosity=2).run(suit)