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 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',datFn = '', filter = 'i', grating = "GI"): super(TestStraylight,self).__init__(methodName) self.pointingData = np.loadtxt(datFn, dtype=np.double) self.filter = filter self.grating = grating 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__': suit = unittest.TestSuite() # case1 = TestStraylight('test_EarthShineFilter',datFn = 'Straylight_test.dat', filter = 'i') # suit.addTest(case1) # case2 = TestStraylight('test_ZodiacalFilter',datFn = 'Straylight_test.dat',filter = 'i') # suit.addTest(case2) # case3 = TestStraylight('test_StarFilter', datFn='Straylight_test.dat', filter='i') # suit.addTest(case3) case4 = TestStraylight('test_GratingStraylight', datFn='Straylight_test.dat', grating = 'GI') suit.addTest(case4) unittest.TextTestRunner(verbosity=2).run(suit)