Skip to content
GitLab
Commit 9928a062 authored by Yan Zhaojun's avatar Yan Zhaojun
Browse files

Delete earthshine.py

parent 580aee74
Show whitespace changes
Inline Side-by-side
csst_mci_sim/earthshine.py deleted 100644 → 0
View file @ 580aee74
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)
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