Commit 22a98299 authored by Yan Zhaojun's avatar Yan Zhaojun
Browse files

Upload New File

parent 3b2de49d
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 stray_light(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]))
star_e = sl.caculateStarLightFilter(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(star_e, path)
# sample as mci wavelength
wave_mci = np.linspace(2500, 11000, 8501) # np.arange(2500, 11000, 1)
f2 = interp1d(earthshine_wave0, earthshine_flux0)
star_stray = f2(wave_mci)
return star_stray
#################################################################
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=stray_light(path,time_jd, x_sat, y_sat, z_sat, ra, dec)
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment