From 91837361ca28b1511397ee31f9449a75860a0ea5 Mon Sep 17 00:00:00 2001 From: Cao Zihuang Date: Thu, 28 Dec 2023 07:16:03 +0000 Subject: [PATCH] Upload New File --- csst_common/epoch.py | 274 +++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 274 insertions(+) create mode 100644 csst_common/epoch.py diff --git a/csst_common/epoch.py b/csst_common/epoch.py new file mode 100644 index 0000000..f24d5f0 --- /dev/null +++ b/csst_common/epoch.py @@ -0,0 +1,274 @@ +""" +Identifier: csst-l1/msc/csst_msc_common/csst_msc_common/catalog/epoch.py +Name: epoch.py +Description: Coordinates system transformation due to motion. +Author: Tianmeng Zhang(zhangtm@nao.cas.cn), Jundan Nie(jdnie@nao.cas.cn), + Zihuang Cao (zhcao@nao.cas.cn) +Created: 2023-11-16 +Modified-History: + 2023-11-16, Tianmeng Zhang, created. + 2023-11-22, Li Shao, change the module structure, improve code format. + 2023-12-12, Zihuang Cao, Add "get_refcat_onOrbit". + 2023-12-19, Zihuang Cao, Change "get_refcat_onOrbit" code, followed Jundan Nie's suggestions. + 2023-12-21, Jundan Nie, impove code. +""" + +from typing import OrderedDict +import astropy.units as units +from astropy.coordinates import get_sun +from astropy.coordinates import SkyCoord +from astropy.coordinates import Distance +from astropy.time import Time +from astropy.table import Table +import numpy as np +from ctypes import c_double, c_int, POINTER, byref, cdll +import sys +import os +from astropy import units as u + +__all__ = ["get_refcat", "get_refcat_onOrbit"] + + +def get_refcat(catalog_table: np.ndarray, + pm_cor: bool = True, + epoch_ref: float = 2016.0, + epoch_now: float = 61605.18976, + ) -> Table: + """ + Function to transform coordination between different epochs. + + Parameters + ---------- + catalog_table : astropy.table.Table + The table contained ra, dec, pmra, pmdec, rv etc. + pm_cor : bool, optional + Proper motion correction. + epoch_now : float, optional + The mjd of exposure time from image.header["EXPSTART"], format='mjd'. + epoch_ref : float, optional + The epoch of reference catalog, format='jyear'. + + Returns + ------- + catalog_table : astropy.table.Table + The corrected catalog contained ra, dec, pmra, pmdec, rv etc. after transfer + and parallax correction. + """ + + ra = catalog_table["ra"] # For the catalog from DFS api + dec = catalog_table["dec"] + parallax = catalog_table["parallax"] + pmra = catalog_table["pmra"] + pmdec = catalog_table["pmdec"] + if pm_cor: + exptime = Time(epoch_now, format='mjd', scale='tcb') + else: + exptime = Time(epoch_ref, format='jyear', scale='tcb') + mask = (parallax > 0.) + parallax_cor = mask * 1 * parallax + (1 - mask * 1) * 1E-10 # mask and change the negative parallax to 1E-10 + catalog_tran = SkyCoord(ra=ra * units.degree, dec=dec * units.degree, frame='icrs', + pm_ra_cosdec=pmra * units.mas / units.yr, + pm_dec=pmdec * units.mas / units.yr, + obstime=Time(epoch_ref, format='jyear', scale='tcb'), + distance=Distance(parallax=parallax_cor * units.mas)) + coo_new = catalog_tran.apply_space_motion(Time(exptime, format='mjd', scale='tcb')) # correct the proper motion + sun_pos = get_sun(Time(exptime, format='mjd', scale='tcb')) + c_sun = SkyCoord(ra=sun_pos.ra, dec=sun_pos.dec, distance=sun_pos.distance) + x_new = coo_new.cartesian.x + c_sun.cartesian.x + y_new = coo_new.cartesian.y + c_sun.cartesian.y + z_new = coo_new.cartesian.z + c_sun.cartesian.z + coord_cor = SkyCoord(x_new, y_new, z_new, unit='pc', representation_type='cartesian') + coord_cor.representation_type = 'spherical' # correct the Annual triangular parallax + return coord_cor + + +def get_refcat_onOrbit( + catalog_table: np.ndarray, + epoch_ref: float = 2016.0, + obs_year: int = 2027, + obs_month: int = 7, + obs_day: int = 19, + obs_hour: int = 4, + obs_minute: int = 33, + obs_second: float = 15.2640001, + orbit_x: float = 0.0, + orbit_y: float = 0.0, + orbit_z: float = 0.0, + orbit_vx: float = 0.0, + orbit_vy: float = 0.0, + orbit_vz: float = 0.0, +) -> Table: + """ + Function to transform coordination between different epochs. + + Parameters + ---------- + catalog_table : astropy.table.Table + The table contained ra, dec, pmra, pmdec, rv etc. + epoch_ref : float, optional + The epoch of reference catalog, format='jyear'. + obs_year : int, optional + The year of the exposure time, derived from image.header["EXPSTART"], format='mjd'. + obs_month : int, optional + The month of the exposure time, derived from image.header["EXPSTART"], format='mjd'. + obs_day : int, optional + The day of the exposure time, derived from image.header["EXPSTART"], format='mjd'. + obs_hour : int, optional + The hour of the exposure time, derived from image.header["EXPSTART"], format='mjd'. + obs_minute : int, optional + The minute of the exposure time, derived from image.header["EXPSTART"], format='mjd'. + obs_second : float, optional + The second of the exposure time, derived from image.header["EXPSTART"], format='mjd'. + orbit_x : flaot + CSST position at X direction on orbit, Unit: km + orbit_y : float + CSST position at Y direction on orbit, Unit: km + orbit_z : float + CSST position at Z direction on orbit, Unit: km + orbit_vx : float + CSST speed at X direction on orbit, Unit: km + orbit_vy : float + CSST speed at Y direction on orbit, Unit: km + orbit_vz : float + CSST speed at Z direction on orbit, Unit: km + + Returns + ------- + catalog_table: astropy.table.Table + The catalog from DFS (gaia dr3) contained ra, dec, pmra, pmdec, rv etc. after transfer and parallax correction. + """ + + nstars = len(catalog_table) + ra = catalog_table["ra"] # For the catalog from DFS api + dec = catalog_table["dec"] + pmra = catalog_table["pmra"] + pmdec = catalog_table["pmdec"] + parallax = catalog_table["parallax"] + mask1 = (parallax > 0.) + parallax_cor = mask1 * 1 * parallax + (1 - mask1 * 1) * 1E-10 # mask and change the negative parallax to 1E-10 + rv = [0.0 for _ in range(nstars)] + rv_cor = rv + # Check input date and time + input_ok = True + if obs_year < 1900: + input_ok = False + if obs_year > 2100: + input_ok = False + if obs_month < 1: + input_ok = False + if obs_month > 12: + input_ok = False + if obs_day < 1: + input_ok = False + if obs_day > 31: + input_ok = False + if obs_hour < 0: + input_ok = False + if obs_hour > 23: + input_ok = False + if obs_minute < 0: + input_ok = False + if obs_minute > 59: + input_ok = False + if obs_second < 0: + input_ok = False + if obs_second > 60: + input_ok = False + # Convert km -> m + input_x = orbit_x * 1000.0 + input_y = orbit_y * 1000.0 + input_z = orbit_z * 1000.0 + input_vx = orbit_vx * 1000.0 + input_vy = orbit_vy * 1000.0 + input_vz = orbit_vz * 1000.0 + # Get path + path = os.path.join(os.path.dirname(__file__), "../lib") + # Inital dynamic lib + if (sys.platform == "darwin"): + print("Running in OSX Darwin ... ") + onOrbitObsLib = cdll.LoadLibrary(os.path.join(path, 'libOnOrbitObs.dylib')) + elif (sys.platform == "linux"): + print("Running in Linux ... ") + onOrbitObsLib = cdll.LoadLibrary(os.path.join(path, 'libOnOrbitObs.so')) + else: + print("Can not run in this system, exit ... ") + sys.exit(1) + onOrbitObsLib.onOrbitObs.restype = c_int + d3 = c_double * 3 + onOrbitObsLib.onOrbitObs.argtypes = [ + c_double, c_double, c_double, c_double, c_double, c_double, + c_int, c_int, c_int, c_int, c_int, c_double, c_double, + POINTER(d3), POINTER(d3), c_int, c_int, c_int, c_int, c_int, + c_double, c_double, POINTER(c_double), POINTER(c_double)] + output_ra_list = list() + output_dec_list = list() + for i in range(nstars): + # check input R.A. and Dec + if ra[i] < 0.0: + output_ra_list.append(-9999.0) + output_dec_list.append(-9999.0) + input_ok = False + continue + if ra[i] > 360.0: + output_ra_list.append(-9999.0) + output_dec_list.append(-9999.0) + input_ok = False + continue + if np.isnan(ra[i]): + output_ra_list.append(-9999.0) + output_dec_list.append(-9999.0) + input_ok = False + continue + if dec[i] < -90.0: + output_ra_list.append(-9999.0) + output_dec_list.append(-9999.0) + input_ok = False + continue + if dec[i] > 90.0: + output_ra_list.append(-9999.0) + output_dec_list.append(-9999.0) + input_ok = False + continue + if np.isnan(dec[i]): + output_ra_list.append(-9999.0) + output_dec_list.append(-9999.0) + input_ok = False + continue + print("Processing : "+str(round((i/(nstars-1)*100.0), 2))+"%") + print(str(i)+"/"+str(nstars-1)) + input_ra = c_double(ra[i]) + input_dec = c_double(dec[i]) + input_pmra = c_double(pmra[i]) + input_pmdec = c_double(pmdec[i]) + input_rv = c_double(rv_cor[i]) + input_parallax = c_double(parallax_cor[i]) + p3 = d3(input_x, input_y, input_z) + v3 = d3(input_vx, input_vy, input_vz) + input_year_c = c_int(obs_year) + input_month_c = c_int(obs_month) + input_day_c = c_int(obs_day) + input_hour_c = c_int(obs_hour) + input_minute_c = c_int(obs_minute) + input_second_c = c_double(obs_second) + DAT = c_double(37.0) + input_epoch = c_double(epoch_ref) + output_ra = c_double(0.0) + output_dec = c_double(0.0) + rs = onOrbitObsLib.onOrbitObs( + input_ra, input_dec, input_parallax, input_pmra, input_pmdec, input_rv, + input_year_c, input_month_c, input_day_c, input_hour_c, input_minute_c, + input_second_c, DAT, byref(p3), byref(v3), input_year_c, input_month_c, + input_day_c, input_hour_c, input_minute_c, input_second_c, input_epoch, + byref(output_ra), byref(output_dec)) + if input_ok and rs == 0: + output_ra_list.append(output_ra.value) + output_dec_list.append(output_dec.value) + print(output_ra.value) + print(output_dec.value) + else: + output_ra_list.append(-9999.0) + output_dec_list.append(-9999.0) + print(-9999.0) + print(-9999.0) + cur_obs_skyCoord_list = SkyCoord(ra=output_ra_list*u.degree, dec=output_dec_list*u.degree, frame='icrs', unit='deg') + return cur_obs_skyCoord_list -- GitLab