straylight.py

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']

filterMirrorEff = {'nuv':0.54,'u':0.68,'g':0.8,'r':0.8, 'i':0.8, 'z':0.8,'y':0.8}

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), ctypes.c_char_p]

        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), ctypes.c_char_p]

        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.argtypes=[ctypes.c_char_p,ctypes.c_char_p,ctypes.c_char_p,ctypes.c_char_p]
        self.deFn = "DE405"
        self.PSTFn = "PST"
        self.RFn = "R"
        self.ZolFn = "Zodiacal"
        self.brightStarTabFn = "BrightGaia_with_csst_mag"

        self.slcdll.Init(str.encode(self.deFn),str.encode(self.PSTFn),str.encode(self.RFn),str.encode(self.ZolFn))

    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, str.encode(self.brightStarTabFn))

        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 caculateEarthshineGrating(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 = band_earth_e2
        py = py2[:]
        if band_earth_e1<band_earth_e2:
            band_earth_e = band_earth_e1
            py = py1[:]
        # 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, py

    def caculateStarLightGrating(self, grating = 'GU', 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, str.encode(self.brightStarTabFn))

        filterPivotWaveList = np.zeros(7)
        bandRangeList = np.zeros(7)
        filterMirrorEffList = np.zeros(7)
        filterNameList = list(filterPivotWave.keys())
        for i in np.arange(7):
            filterPivotWaveList[i] = filterPivotWave[filterNameList[i]]
            filterMirrorEffList[i] = filterMirrorEff[filterNameList[i]]
            brange = bandRange[filterNameList[i]]
            bandRangeList[i] = brange[1] - brange[0]
        filterFlux_lamb = star_e1[:]/bandRangeList/filterMirrorEffList/(h*c/(filterPivotWaveList*1e-10))
        filterFlux_lambi = interpolate.interp1d(filterPivotWaveList,filterFlux_lamb,fill_value="extrapolate")

        lamb_g = np.arange(bandRange[grating][0], bandRange[grating][1],1)
        flux_g = filterFlux_lambi(lamb_g)
        # flux_total_g = np.trapz(flux_g,lamb_g)
        zoldialcal_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_g)
            beamsZol = np.trapz(flux_g*thp,lamb_g)*pixel_size_phy*1e-6*pixel_size_phy*1e-6
            zoldialcal_grating = zoldialcal_grating + beamsZol
            # print(beamsZol)



        # 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 zoldialcal_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 caculateStrayLightFilter(self, filter = 'i', pixel_size_phy = 10, pixel_scale = 0.074, sun_pos = np.array([0,0,0])):
        e1,py = self.caculateEarthShineFilter(filter = filter, pixel_size_phy = pixel_size_phy)
        e2, _ = self.calculateZaodiacalFilter2(filter = 'i', sun_pos=sun_pos, pixelsize = pixel_scale)
        e3 = self.caculateStarLightFilter(filter = 'i',pointYaxis = py, pixel_size_phy = pixel_size_phy)
        
        return e1+e2+e3
    
    def caculateStrayLightGrating(self, grating = 'GI', pixel_size_phy = 10, normFilter_es = 'g', sun_pos = np.array([0,0,0])):
        e1,py = self.caculateEarthshineGrating(grating = grating, pixel_size_phy = pixel_size_phy, normFilter = normFilter_es)
        e2 = self.caculateStarLightGrating(grating = grating, pointYaxis = py)
        spec, _ = self.calculateZaodicalSpec(longitude = self.ecliptic.lon.value, latitude = self.ecliptic.lat.value, sun_pos = sun_pos)
        
        return e1+e2, spec
 

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 

    e1,py = sl.caculateEarthshineGrating(grating = 'GI', pixel_size_phy = 10, normFilter = 'g')

    # e2 = sl.caculateStarLightGrating(grating = 'GV', pointYaxis = py)
    e2 = sl.caculateStarLightGrating(grating = 'GI', pointYaxis = py)
    e4 = sl.caculateStarLightFilter(filter = 'i',pointYaxis = py)

    e5=sl.caculateStrayLightFilter(filter = 'i', pixel_size_phy = 10, pixel_scale = 0.074, sun_pos = surveylist[i,6:9])
    e6,_=sl.caculateStrayLightGrating(grating = 'GI', normFilter_es = 'g', sun_pos = surveylist[i,6:9])

    sky_pix[i,0] = e1
    sky_pix[i,1] = e2
    sky_pix[i,2] = e3
    sky_pix[i,3] = e4
    print(e1+e2,e1_+e3+e4,e5,e6)

    # print(e1,e2,e3,e4)