test_Straylight.py

#
# 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 observation_sim.sky_background 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_msc_sim/test_sls_and_straylight')
        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(r'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)