Upload New File

import numpy as np

import julian

from datetime import datetime

from astropy.time import Time

from astropy.coordinates import get_sun

from astropy.coordinates import SkyCoord

import astropy.coordinates as coord

import pandas as pd

from astropy import units as u

from scipy.interpolate import interp1d

from scipy import interpolate

import ctypes

def transRaDec2D(ra, dec):

    # radec转为竞天程序里的ob, 赤道坐标系下的笛卡尔三维坐标xyz.

    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 earth_angle(time_jd, x_sat, y_sat, z_sat, ra_obj, dec_obj):

    ra_sat = np.arctan2(y_sat, x_sat) / np.pi * 180

    dec_sat = np.arctan2(z_sat, np.sqrt(x_sat**2+y_sat**2)) / np.pi * 180

    radec_sat = SkyCoord(ra=ra_sat*u.degree,

                         dec=dec_sat*u.degree, frame='gcrs')

    lb_sat = radec_sat.transform_to('geocentrictrueecliptic')

    # get the obj location

    radec_obj = SkyCoord(ra=ra_obj*u.degree,

                         dec=dec_obj*u.degree, frame='gcrs')

    lb_obj = radec_obj.transform_to('geocentrictrueecliptic')

    # calculate the angle between sub-satellite point and the earth side

    earth_radius = 6371     # km

    sat_height = np.sqrt(x_sat**2 + y_sat**2 + z_sat**2)

    angle_a = np.arcsin(earth_radius/sat_height) / np.pi * 180

    # calculate the angle between satellite position and the target position

    angle_b = lb_sat.separation(lb_obj)

    # calculat the earth angle

    angle = 180 - angle_a - angle_b.degree

    return angle

###############################################################################

def ill2flux(E, path):

    # use template from sky_bkg (background_spec_hst.dat)

    filename = path+'MCI_inputData/refs/background_spec_hst.dat'

    cat_spec = pd.read_csv(filename, sep='\s+', header=None, comment='#')

    wave0 = cat_spec[0].values       # A

    spec0 = cat_spec[2].values      # erg/s/cm^2/A/arcsec^2

    # convert erg/s/cm^2/A/arcsec^2 to erg/s/cm^2/A/sr

    flux1 = spec0 / (1/4.25452e10)

    # convert erg/s/cm^2/A/sr to W/m^2/sr/um

    flux2 = flux1 * 10

    # 对接收面积积分，输出单位 W/m^2/nm

    D = 2   # meter

    f = 28  # meter, 焦距，转换关系来源于王维notes.

    flux3 = flux2 * np.pi * D**2 / 4 / f**2 / 10**3

    f = interp1d(wave0, flux3)

    wave_range = np.arange(3800, 7800)

    flux3_mean = f(wave_range)

    delta_lamba = 0.1   # nm

    E0 = np.sum(flux3_mean * delta_lamba)

    factor = E / E0

    spec_scaled = factor * spec0

    return wave0, spec_scaled

##############################################################

def earthshine(path, time_jd, x_sat, y_sat, z_sat, ra, dec):

    # EarthShine from straylight

    sl = StrayLight(path, jtime=time_jd, sat=np.array(

        [x_sat, y_sat, z_sat]), radec=np.array([(ra*u.degree).value, (dec*u.degree).value]))

    earth_e = sl.caculateEarthShineFilter(filter='r')

    angle_earth = earth_angle(time_jd, x_sat, y_sat, z_sat, ra, dec)

    if angle_earth < 0:

        earth_e = 0

    earthshine_wave0, earthshine_flux0 = ill2flux(earth_e, path)

    # sample as mci wavelength

    wave_mci = np.linspace(2500, 11000, 8501)  # np.arange(2500, 11000, 1)

    f2 = interp1d(earthshine_wave0, earthshine_flux0)

    earthshine = f2(wave_mci)

    return earthshine

#################################################################

class StrayLight(object):

    def __init__(self, path, jtime=2460843., sat=np.array([0, 0, 0]), radec=np.array([0, 0])):

        self.path = path

        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(

            self.path+'MCI_inputData/refs/libstraylight.so')  # 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 = self.path+"MCI_inputData/refs/DE405"

        self.PSTFn = self.path+"MCI_inputData/refs/PST"

        self.RFn = self.path+"MCI_inputData/refs/R"

        self.ZolFn = self.path+"MCI_inputData/refs/Zodiacal"

        self.brightStarTabFn = self.path+"MCI_inputData/refs/BrightGaia_with_csst_mag"

        self.slcdll.Init(str.encode(self.deFn), str.encode(

            self.PSTFn), str.encode(self.RFn), str.encode(self.ZolFn))

    # def caculateStarLightFilter(self, filter='i'):

    #     filterIndex = {'nuv': 0, 'u': 1, 'g': 2,

    #                    'r': 3, 'i': 4, 'z': 5, 'y': 6}

    #     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)

    #     star_e1 = (ctypes.c_double*7)()

    #     self.slcdll.PointSource(self.jtime, sat, ob, py1,

    #                             star_e1, str.encode(self.brightStarTabFn))

    #     star_e2 = (ctypes.c_double*7)()

    #     self.slcdll.PointSource(self.jtime, sat, ob, py2,

    #                             star_e2, str.encode(self.brightStarTabFn))

    #     band_star_e1 = star_e1[:][filterIndex[filter]]

    #     band_star_e2 = star_e2[:][filterIndex[filter]]

    #     return max(band_star_e1, band_star_e2)

###############################################################################

    def caculateEarthShineFilter(self, filter='i'):

        filterIndex = {'nuv': 0, 'u': 1, 'g': 2,

                       'r': 3, 'i': 4, 'z': 5, 'y': 6}

        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)          # e[7]代表7个波段的照度

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

        return max(band_earth_e1, band_earth_e2)

###############################################################################

# test

# path='/home/yan/MCI/'

# time_jd = 2460417.59979167

# x_sat = -4722.543136

# y_sat = -1478.219213

# z_sat = 4595.402769

# ra = 116.18081536720157

# dec= 39.42316681016602

# earthshine0=earthshine(path,time_jd, x_sat, y_sat, z_sat, ra, dec)