Skip to content

GitLab

Commit 5510ed4d authored by Niu Liu's avatar Niu Liu
Browse files

fix : fix minor bugs

parent 4ec120c5
Hide whitespace changes
Inline Side-by-side
aux_files/crf-format.txt 0 → 100644
View file @ 5510ed4d
------------
IERS ICRS-PC
------------
FORMAT FOR SUBMITTING SETS OF RADIOSOURCE COORDINATES
Field Contents all fields separated by one blank (1)
1 IAU name
2 Alternate name
3 hours of right ascension (2)
4 minutes of right ascension (2)
5 seconds of right ascension (2)
6 degrees of declination (3)
7 minutes of declination (3)
8 seconds of declination (3)
9 uncertainty on right ascension (s)
10 uncertainty on declination (")
11 correlation : right ascension, declination
12 weighted mean observation date (MJD) (4)
13 date of first observation (MJD)
14 date of last observation (MJD)
15 number of sessions
16 number of delays
17 number of delay rates
(5)
Notes
1 - All the 17 fields should be filled (with 0 when the information is
lacking).
2 - These three fields can be replaced by one field, giving the right
ascensions in degrees and decimal fraction (from 0. to 360.)
3 - These three fields can be replaced by one field, giving the declination
in degrees and decimal fraction (from -90. to +90.).
4 - The epoch of the catalog should be given in the general description of
the analysis.
5 - Any additional parameter(s) you feel appropriate to the description of
the result.
comp_cat.py deleted 100644 → 0
View file @ 4ec120c5
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
# File name: comp_cat.py
"""
Created on Fri Mar 22 15:58:02 2019
@author: Neo(liuniu@smail.nju.edu.cn)
"""
from astropy.table import Table, join, Column
import astropy.units as u
from astropy.coordinates import Angle
import numpy as np
from numpy import cos, sqrt
from functools import reduce
import time
__all__ = {"calc_cat_offset"}
# ----------------------------- FUNCTIONS -----------------------------
def root_sum_square(x, y):
"""Calculate the root-sum-square."""
return np.sqrt(x**2 + y**2)
def nor_sep_calc(dRA, dRA_err, dDC, dDC_err, C):
'''Calculate the normalized seperation.
Parameters
----------
dRA/dDC : Right Ascension / Declination differences in micro-as
dRA_err/dDC_err : formal uncertainty of dRA*cos(Dec)/dDC in micro-as
C : correlation coeffient between dRA*cos(Dec) and dDC.
Returns
----------
ang_sep : angular seperation, in micro-as
X_a / X_d : normalized coordinate differences in RA / DC, unit-less
X : Normalized separations, unit-less.
'''
# Angular seperations
ang_sep = sqrt(dRA**2 + dDC**2)
# Normalised coordinate difference
X_a = dRA / dRA_err
X_d = dDC / dDC_err
# Normalised separation - Mignard's statistics (considering covariance)
X = np.zeros_like(X_a)
for i, (X_ai, X_di, Ci) in enumerate(zip(X_a, X_d, C)):
if Ci == -1.:
Ci = -0.99
if Ci == 1.0:
Ci = 0.99
# print(Ci)
wgt = np.linalg.inv(np.mat([[1, Ci],
[Ci, 1]]))
Xmat = np.mat([X_ai, X_di])
X[i] = sqrt(reduce(np.dot, (Xmat, wgt, Xmat.T)))
return ang_sep, X_a, X_d, X
def calc_cat_offset(t_cat1, t_cat2, souname="iers_name"):
"""Calculate the radio source position differences between two catalogs.
Parameters
----------
t_cat1, t_cat2 : astropy.table object
radio source positions from two catalogs
Return
------
t_cat_oft : astropy.table object
radio source position differences
"""
# Copy the original tables and keep only the source position information.
t_cat3 = Table(t_cat1)
t_cat3.keep_columns([souname, "ra",
"dec", "ra_err", "dec_err", "ra_dec_corr",
"used_sess", "used_obs"])
t_cat4 = Table(t_cat2)
t_cat4.keep_columns([souname, "ra", "dec",
"ra_err", "dec_err", "ra_dec_corr"])
# Cross-match between two tables
soucom = join(t_cat3, t_cat4, keys=souname)
print("There are %d and %d sources in two catalogs, respectively,"
"between which %d are common."
% (len(t_cat1), len(t_cat2), len(soucom)))
# Calculate the offset and the uncertainties
arcfac = cos(Angle(soucom["dec_1"]).radian)
dra = (soucom["ra_1"] - soucom["ra_2"]) * arcfac
ddec = soucom["dec_1"] - soucom["dec_2"]
dra_err = root_sum_square(soucom["ra_err_1"], soucom["ra_err_2"])
ddec_err = root_sum_square(soucom["dec_err_1"], soucom["dec_err_2"])
dra_ddec_cov = soucom["ra_err_1"] * soucom["dec_err_1"] * soucom["ra_dec_corr_1"] + \
soucom["ra_err_2"] * soucom["dec_err_2"] * soucom["ra_dec_corr_2"]
dra_ddec_corr = dra_ddec_cov / \
root_sum_square(soucom["ra_err_1"], soucom["dec_err_1"]) / \
root_sum_square(soucom["ra_err_2"], soucom["dec_err_2"])
# Convert the unit
dra.convert_unit_to(u.uas)
dra_err.convert_unit_to(u.uas)
ddec.convert_unit_to(u.uas)
ddec_err.convert_unit_to(u.uas)
dra_ddec_corr.unit = None
dra_ddec_cov.unit = u.mas * u.mas
dra_ddec_cov.convert_unit_to(u.uas * u.uas)
# Calculate the angular seperation
ang_sep, X_a, X_d, X = nor_sep_calc(
dra, dra_err, ddec, ddec_err, dra_ddec_corr)
# sou_nb = len(soucom)
ang_sep = Column(ang_sep, unit=u.uas)
X_a = Column(X_a)
X_d = Column(X_d)
X = Column(X)
# Add these columns to the combined table.
t_cat_oft = Table([soucom[souname],
soucom["ra_1"], soucom["dec_1"],
dra, dra_err, ddec, ddec_err,
dra_ddec_cov, dra_ddec_corr,
ang_sep, X, X_a, X_d],
names=["ivs_name", "ra", "dec",
"dra", "dra_err", "ddec", "ddec_err",
"dra_ddec_cov", "dra_ddec_corr",
"ang_sep", "nor_sep",
"nor_sep_ra", "nor_sep_dec"])
return t_cat_oft
# --------------------------------- END --------------------------------
if __name__ == '__main__':
print("Nothing to do!")
pass
# --------------------------------- END --------------------------------
comp_eop.py deleted 100644 → 0
View file @ 4ec120c5
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
# File name: eob_comp.py
import time
"""
Created on Wed Feb 27 14:59:41 2019
@author: Neo(liuniu@smail.nju.edu.cn)
Compare the EOP series and calculate the offsets (or differences).
"""
from astropy.table import Table, join, Column
import astropy.units as u
from astropy.units import cds
import numpy as np
import sys
import os
import time
__all__ = {"calc_eop_offset", "save_eop_offset"}
# ----------------------------- FUNCTIONS -----------------------------
def root_sum_square(x, y):
"""Calculate the root-sum-square."""
return np.sqrt(x**2 + y**2)
def save_eop_offset(eopoft, oftfile):
"""Save the eop offset series into a text file.
Parameters
----------
eopoft : Table object
eop offset series
oftfile : string
name of file to store the data
"""
# Header
eopoft.meta["comments"] = [
" EOP Offset series",
" Columns Units Meaning",
" 1 day Time Tag for polar motion and UT1 (MJD)",
" 2 day Time Tag for Nutation (MJD)",
" 3 uas offset of X pole coordinate",
" 4 uas offset of Y pole coordinate",
" 5 musec offset of UT1",
" 6 musec offset of LOD",
" 7 uas offset of dX of Nutation offsets",
" 8 uas offset of dY of Nutation offsets",
" 9 uas formal uncertainty for offset of X pole coordinate",
" 10 uas formal uncertainty for offset of Y pole coordinate",
" 11 musec formal uncertainty for offset of UT1",
" 12 musec formal uncertainty for offset of LOD",
" 13 uas formal uncertainty for offset of Nutation dX",
" 14 uas formal uncertainty for offset of Nutation dY",
" Created date: %s." % time.strftime("%d/%m/%Y", time.localtime())]
eopoft.write(oftfile, format="ascii.fixed_width_no_header",
exclude_names=["db_name"],
formats={"epoch_pmr": "%14.6f", "epoch_nut": "%14.6f",
"dxp": "%+8.1f", "dyp": "%+8.1f",
"dut": "%+8.1f", "dlod": "%+8.3f",
"ddX": "%+8.1f", "ddY": "%+8.1f",
"dxp_err": "%8.1f", "dyp_err": "%8.1f",
"dut_err": "%8.1f", "dlod_err": "%8.3f",
"ddX_err": "%8.1f", "ddY_err": "%8.1f"},
delimiter="", overwrite=True)
def read_eop_offset(oftfile):
"""Read EOP offset data.
Parameters
----------
oftfile : string
EOP offset file
Returns
----------
eopoft : astropy.table object
"""
if not os.path.isfile(oftfile):
print("Couldn't find the file", oftfile)
sys.exit()
eopoft = Table.read(oftfile, format="ascii",
names=["epoch_pmr", "epoch_nut",
"dxp", "dyp", "dut", "dlod", "ddX", "ddY",
"dxp_err", "dyp_err", "dut_err", "dlod_err",
"ddX_err", "ddY_err"])
# Add unit information
eopoft["epoch_pmr"].unit = cds.MJD
eopoft["epoch_nut"].unit = cds.MJD
eopoft["dxp"].unit = u.uas
eopoft["dyp"].unit = u.uas
eopoft["dxp_err"].unit = u.uas
eopoft["dyp_err"].unit = u.uas
eopoft["dut"].unit = u.second / 1e6
eopoft["dut_err"].unit = u.second / 1e6
eopoft["dut"].unit = u.second / 1e6
eopoft["dut_err"].unit = u.second / 1e6
eopoft["ddX"].unit = u.uas
eopoft["ddY"].unit = u.uas
eopoft["ddX_err"].unit = u.uas
eopoft["ddY_err"].unit = u.uas
return eopoft
def calc_eop_offset(eop1, eop2, oftfile=None):
"""Calculate the EOP difference series between two solutions.
Parameters
----------
eop1, eop2 : astropy.table object
EOP series from two solutions
Return
------
eopoft : astropy.table object
EOP difference series
"""
# Copy the original tables and keep only the EOP information
eop3 = Table(eop1)
eop3.keep_columns(["db_name", "epoch_pmr", "epoch_nut",
"xp", "yp", "ut1_tai", "dX", "dY", "lod",
"xp_err", "yp_err", "ut1_err", "dX_err", "dY_err",
"lod_err"])
eop4 = Table(eop2)
eop4.keep_columns(["db_name",
"xp", "yp", "ut1_tai", "dX", "dY", "lod",
"xp_err", "yp_err", "ut1_err", "dX_err", "dY_err",
"lod_err"])
# Cross-match between two tables
eopcom = join(eop3, eop4, keys="db_name")
print("There are %d and %d points in series 1 and series 2, respectively,"
"between which %d are common."
% (len(eop1), len(eop2), len(eopcom)))
# Calculate the offset and the uncertainties
dxp = eopcom["xp_1"] - eopcom["xp_2"]
dyp = eopcom["yp_1"] - eopcom["yp_2"]
dut = eopcom["ut1_tai_1"] - eopcom["ut1_tai_2"]
ddX = eopcom["dX_1"] - eopcom["dX_2"]
ddY = eopcom["dY_1"] - eopcom["dY_2"]
dlod = eopcom["lod_1"] - eopcom["lod_2"]
dxperr = root_sum_square(eopcom["xp_err_1"], eopcom["xp_err_2"])
dyperr = root_sum_square(eopcom["yp_err_1"], eopcom["yp_err_2"])
duterr = root_sum_square(eopcom["ut1_err_1"], eopcom["ut1_err_2"])
ddXerr = root_sum_square(eopcom["dX_err_1"], eopcom["dX_err_2"])
ddYerr = root_sum_square(eopcom["dY_err_1"], eopcom["dY_err_2"])
dloderr = root_sum_square(eopcom["lod_err_1"], eopcom["lod_err_2"])
# Convert the unit
# Time tag
# from astropy.time import Time
# t_pmr_mjd = Time(eopcom["epoch_pmr"], format="mjd")
# t_pmr = Column(t_pmr_mjd.jyear, unit=u.year)
#
# t_nut_mjd = Time(eopcom["epoch_nut"], format="mjd")
# t_nut = Column(t_nut_mjd.jyear, unit=u.year)
# Polar motion (as -> uas)
dxp.convert_unit_to(u.uas)
dxperr.convert_unit_to(u.uas)
dyp.convert_unit_to(u.uas)
dyperr.convert_unit_to(u.uas)
# UT1-UTC and lod (s -> us)
dut.convert_unit_to(u.second / 1e6)
duterr.convert_unit_to(u.second / 1e6)
dlod.convert_unit_to(u.second / 1e6)
dloderr.convert_unit_to(u.second / 1e6)
# Nutation offset (as -> uas)
ddX.convert_unit_to(u.uas)
ddXerr.convert_unit_to(u.uas)
ddY.convert_unit_to(u.uas)
ddYerr.convert_unit_to(u.uas)
# Add these columns to the combined table.
eopoft = Table([eopcom["db_name"], eopcom["epoch_pmr"],
eopcom["epoch_nut"],
dxp, dyp, dut, dlod, ddX, ddY,
dxperr, dyperr, duterr, dloderr, ddXerr, ddYerr],
names=["db_name", "epoch_pmr", "epoch_nut",
"dxp", "dyp", "dut", "dlod", "ddX", "ddY",
"dxp_err", "dyp_err", "dut_err", "dlod_err",
"ddX_err", "ddY_err"])
# Save the EOP offset series
if oftfile is not None:
print("Save the EOP offset series in", oftfile)
save_eop_offset(eopoft, oftfile)
return eopoft
if __name__ == '__main__':
print('Nothing to do!')
# --------------------------------- END --------------------------------
comp_trf.py deleted 100644 → 0
View file @ 4ec120c5
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
# File name: comp_crf.py
"""
Created on Sun Mar 24 16:14:45 2019
@author: Neo(liuniu@smail.nju.edu.cn)
"""
from astropy.table import Table, join, Column
import astropy.units as u
import numpy as np
from numpy import cos, sqrt
from functools import reduce
import time
__all__ = {"calc_sta_offset", "calc_vel_offset"}
# ----------------------------- FUNCTIONS -----------------------------
def root_sum_square(x, y):
"""Calculate the root-sum-square."""
return np.sqrt(x**2 + y**2)
def calc_sta_offset(t_sta1, t_sta2):
"""Calculate station position differences.
Parameters
----------
t_sta1, t_sta2 : astropy.table object
station positions from two solutions
Return
------
t_sta_oft : astropy.table object
station position differences
"""
# Copy the original tables and keep only the source position information.
t_sta3 = Table(t_sta1)
t_sta3.keep_columns(["station", "xp", "xp_err",
"yp", "yp_err", "zp", "zp_err",
"xp_yp_corr", "xp_zp_corr", "yp_zp_corr"])
t_sta4 = Table(t_sta2)
t_sta4.keep_columns(["station", "xp", "xp_err",
"yp", "yp_err", "zp", "zp_err",
"xp_yp_corr", "xp_zp_corr", "yp_zp_corr"])
# Cross-match between two tables
t_sta_com = join(t_sta3, t_sta4, keys="station")
print("There are %d and %d stations in two sets, respectively, "
"between which %d are common."
% (len(t_sta1), len(t_sta2), len(t_sta_com)))
# Calculate the offset and the uncertainties
dxp = t_sta_com["xp_1"] - t_sta_com["xp_2"]
dyp = t_sta_com["yp_1"] - t_sta_com["yp_2"]
dzp = t_sta_com["zp_1"] - t_sta_com["zp_2"]
dxp_err = root_sum_square(t_sta_com["xp_err_1"], t_sta_com["xp_err_2"])
dyp_err = root_sum_square(t_sta_com["yp_err_1"], t_sta_com["yp_err_2"])
dzp_err = root_sum_square(t_sta_com["zp_err_1"], t_sta_com["zp_err_2"])
# Covariance
dxp_dyp_cov = t_sta_com["xp_err_1"] * t_sta_com["yp_err_1"] \
* t_sta_com["xp_yp_corr_1"] + t_sta_com["xp_err_2"] \
* t_sta_com["yp_err_2"] * t_sta_com["xp_yp_corr_2"]
dxp_dzp_cov = t_sta_com["xp_err_1"] * t_sta_com["zp_err_1"] \
* t_sta_com["xp_zp_corr_1"] + t_sta_com["xp_err_2"] \
* t_sta_com["zp_err_2"] * t_sta_com["xp_zp_corr_2"]
dyp_dzp_cov = t_sta_com["yp_err_1"] * t_sta_com["zp_err_1"] \
* t_sta_com["yp_zp_corr_1"] + t_sta_com["yp_err_2"] \
* t_sta_com["zp_err_2"] * t_sta_com["yp_zp_corr_2"]
# Add these columns to the combined table.
t_sta_oft = Table([t_sta_com["station"], t_sta_com["xp_1"],
t_sta_com["yp_1"], t_sta_com["zp_1"],
dxp, dxp_err, dyp, dyp_err, dzp, dzp_err,
dxp_dyp_cov, dxp_dzp_cov, dyp_dzp_cov],
names=["station", "xp", "yp", "zp",
"dxp", "dxp_err", "dyp", "dyp_err",
"dzp", "dzp_err",
"dxp_dyp_cov", "dxp_dzp_cov", "dyp_dzp_cov"])
return t_sta_oft
def calc_vel_offset(t_vel1, t_vel2):
"""Calculate station position differences.
Parameters
----------
t_vel1, t_vel2 : astropy.table object
station positions from two solutions
Return
------
t_vel_oft : astropy.table object
station position differences
"""
# Copy the original tables and keep only the source position information.
t_vel3 = Table(t_vel1)
t_vel3.keep_columns(["station", "xv", "xv_err",
"yv", "yv_err", "zv", "zv_err",
"xv_yv_corr", "xv_zv_corr", "yv_zv_corr"])
t_vel4 = Table(t_vel2)
t_vel4.keep_columns(["station", "xv", "xv_err",
"yv", "yv_err", "zv", "zv_err",
"xv_yv_corr", "xv_zv_corr", "yv_zv_corr"])
# Cross-match between two tables
t_vel_com = join(t_vel3, t_vel4, keys="station")
print("There are %d and %d stations in two sets, respectively, "
"between which %d are common."
% (len(t_vel1), len(t_vel2), len(t_vel_com)))
# Calculate the offset and the uncertainties
dxv = t_vel_com["xv_1"] - t_vel_com["xv_2"]
dyv = t_vel_com["yv_1"] - t_vel_com["yv_2"]
dzv = t_vel_com["zv_1"] - t_vel_com["zv_2"]
dxv_err = root_sum_square(t_vel_com["xv_err_1"], t_vel_com["xv_err_2"])
dyv_err = root_sum_square(t_vel_com["yv_err_1"], t_vel_com["yv_err_2"])
dzv_err = root_sum_square(t_vel_com["zv_err_1"], t_vel_com["zv_err_2"])
# Covariance
dxv_dyv_cov = t_vel_com["xv_err_1"] * t_vel_com["yv_err_1"] \
* t_vel_com["xv_yv_corr_1"] + t_vel_com["xv_err_2"] \
* t_vel_com["yv_err_2"] * t_vel_com["xv_yv_corr_2"]
dxv_dzv_cov = t_vel_com["xv_err_1"] * t_vel_com["zv_err_1"] \
* t_vel_com["xv_zv_corr_1"] + t_vel_com["xv_err_2"] \
* t_vel_com["zv_err_2"] * t_vel_com["xv_zv_corr_2"]
dyv_dzv_cov = t_vel_com["yv_err_1"] * t_vel_com["zv_err_1"] \
* t_vel_com["yv_zv_corr_1"] + t_vel_com["yv_err_2"] \
* t_vel_com["zv_err_2"] * t_vel_com["yv_zv_corr_2"]
# Add these columns to the combined table.
t_vel_oft = Table([t_vel_com["station"], t_vel_com["xv_1"],
t_vel_com["yv_1"], t_vel_com["zv_1"],
dxv, dxv_err, dyv, dyv_err, dzv, dzv_err,
dxv_dyv_cov, dxv_dzv_cov, dyv_dzv_cov],
names=["station", "xv", "yv", "zv",
"dxv", "dxv_err", "dyv", "dyv_err",
"dzv", "dzv_err",
"dxv_dyv_cov", "dxv_dzv_cov", "dyv_dzv_cov"])
return t_vel_oft
def calc_trf_offset(t_trf1, t_trf2):
"""Calculate station position and velocity differences.
Parameters
----------
t_trf1, t_trf2 : astropy.table object
station positions from two solutions
Return
------
t_trf_oft : astropy.table object
station position and velocity differences
"""
# Copy the original tables and keep only the source position information.
t_trf3 = Table(t_trf1)
t_trf3.keep_columns(["station", "xp", "xp_err",
"yp", "yp_err", "zp", "zp_err",
"xp_yp_corr", "xp_zp_corr", "yp_zp_corr",
"xv", "xv_err", "yv", "yv_err", "zv", "zv_err",
"xv_yv_corr", "xv_zv_corr", "yv_zv_corr"])
t_trf4 = Table(t_trf2)
t_trf4.keep_columns(["station", "xp", "xp_err",
"yp", "yp_err", "zp", "zp_err",
"xp_yp_corr", "xp_zp_corr", "yp_zp_corr",
"xv", "xv_err", "yv", "yv_err", "zv", "zv_err",
"xv_yv_corr", "xv_zv_corr", "yv_zv_corr"])
# Cross-match between two tables
t_trf_com = join(t_trf3, t_trf4, keys="station")
print("There are %d and %d stations in two sets, respectively, "
"between which %d are common."
% (len(t_trf1), len(t_trf2), len(t_trf_com)))
# Calculate the positional offset and the uncertainties
dxp = t_trf_com["xp_2"] - t_trf_com["xp_1"]
dyp = t_trf_com["yp_2"] - t_trf_com["yp_1"]
dzp = t_trf_com["zp_2"] - t_trf_com["zp_1"]
dxp_err = root_sum_square(t_trf_com["xp_err_1"], t_trf_com["xp_err_2"])
dyp_err = root_sum_square(t_trf_com["yp_err_1"], t_trf_com["yp_err_2"])
dzp_err = root_sum_square(t_trf_com["zp_err_1"], t_trf_com["zp_err_2"])