test_Straylight.py

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)