Module ephemeris_library.rotation_libraries.lib_j2k_geo_rotation
Expand source code
import os
import sys
import datetime
import numpy
from ephemeris_library.rotation_libraries import lib_quaternion as lib_quat
####
# Global Pathing:
Script_Directory = os.path.dirname(os.path.realpath(__file__))
IERS_File_Path = os.path.join(Script_Directory, "iers_eop.txt")
# Global Shared Vars:
IERS_Date, IERS_MET, IERS_X_Pole, IERS_Y_Pole, IERS_UT1_Minus_UTC, IERS_Delta_Psi, IERS_Delta_Eplison = [],[],[],[],[],[],[]
# DT, MET, Arcseconds, Arcseconds, Seconds, Arcseconds, Arcseconds.
##############################
# Global Convenience Functions:
def DT_to_MET( DT ): # Doesn't -have- to be UTC, but official MET follows UTC, like Unix-time.
return numpy.array([ ( dt - datetime.datetime(1968,5,24)).total_seconds() for dt in DT ])
def MET_to_DT( MET ):
return [ datetime.datetime(1968,5,24) + datetime.timedelta(seconds=met) for met in MET]
def MJD_to_DT( MDJ ):
# Epoch of 1858-11-17 00:00:00 UTC
return [ datetime.datetime(1858,11,17) + datetime.timedelta(days=mjd) for mjd in MDJ ]
def DT_to_MJD( DT ):
return numpy.array([ ( dt - datetime.datetime(1858,11,17)).total_seconds() // 86400 for dt in DT ])
##############################
# Coefficient-File Functions.
def Update_IERS_Data():
'''
Updates the IERS file found in the
lib_j2k_geo_rotation.IERS_File_Path variable.
(Default: <path>/rotation_libraries/iers_eop.txt
Documentation for the file found in this library's checkout:
<path>/rotation_libraries/iers_eop_documentation.txt file
Raises ConnectionRefusedError on inability to connect/update the file.
'''
# Standard Rapid EOP Data since 01. January 1992 (IAU1980)
# From https://www.iers.org/IERS/EN/DataProducts/EarthOrientationData/eop.html - finals.data (IAU1980)
# Metadata: https://datacenter.iers.org/versionMetadata.php?filename=latestVersionMeta/8_FINALS.DATA_IAU1980_V2013_018.txt
IERS_Default_URL = "https://datacenter.iers.org/data/8/finals.data"
try:
import urllib.request, ssl
# The IERS uses a self-signed certificate, which the URL-lib will raise exceptions about.
# As such, need to suppress the warning by overriding the default SSL-configuration context.
ssl._create_default_https_context = ssl._create_unverified_context
# Make the request:
URL_Request = urllib.request.urlopen( IERS_Default_URL )
IERS_Data = URL_Request.read().decode('utf-8') # Read & Decode the data
# And overwrite the IERS Text File.
FH = open( IERS_File_Path,'w' )
FH.write( IERS_Data )
FH.close()
# And update the contents!
Read_IERS_Data()
return
except:
pass # Don't raise exception here, re-raise below.
print("!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!")
print("!!! Could not fetch Earth Orientation Data from the IERS servers")
print("!!! Please update the IERS file at {:s}".format(IERS_File_Path))
print("!!! Using the EOP 14 C04 (IAU1980) publication.")
print("!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!")
raise ConnectionRefusedError()
def Sanity_Check_IERS_Data( MET_List ):
'''
Given a MET list, load & read the IERS EOP data,
& make sure that all times are within the given range.
Raises a ValueError should the requested list fall outside the constants file range.
'''
if len(IERS_Date) == 0: Read_IERS_Data() # Load the globals, if needed.
# If the requested times are outside the range of the IERS dataset, update the file.
if numpy.min(MET_List) < numpy.min(IERS_MET) or numpy.max(MET_List) > numpy.max(IERS_MET):
raise ValueError("Requested MET outside the range of the IERS_EOP file. Please update with .Update_IERS_Data()")
def Read_IERS_Data( iers_file_path = IERS_File_Path ):
'''
Loads the IERS coefficient data into memory.
Args:
iers_file_path: An abolute path to the given IERS file.
Defaults to ./<Script_Path>/iers_eop.txt
'''
# The file format used by IERS is... odd. Read the 'IERS_EOP_Documentation' file next to the EOP file for more details.
global IERS_Date, IERS_MET, IERS_X_Pole, IERS_Y_Pole, IERS_UT1_Minus_UTC, IERS_Delta_Psi, IERS_Delta_Eplison
FD = open( iers_file_path ).readlines()
# Will raise a FileNotFoundError if could not read the file.
# Read the MJD as to filter the input solutions.
MJD = numpy.array( [Line[7:15].strip() for Line in FD], dtype=float )
# Filter for entries from ~1 week before launch, to 1 week after the current time. Just reduced memory load.
Launch_DT, Current_DT = datetime.datetime( 2013, 9, 29 ), datetime.datetime.now().replace( hour=0, minute=0, second=0, microsecond=0)
Launch_MJD, Current_MJD = DT_to_MJD( [Launch_DT, Current_DT] )
_Valid_Mask = numpy.logical_and( Launch_MJD -7 <= MJD, MJD <= Current_MJD+7 )
FD = [ FD[i] for i in range(len(FD)) if _Valid_Mask[i] ]
# Don't bother with year/month, just read MJD.
MJD = [ Line[ 7: 15].strip() for Line in FD ] # Guaranteed populted for each line.
# Read the Bulletin A values.
Bulletin_A_X_Pole = [ Line[ 18: 27].strip() for Line in FD ] # Published about 6-10 months ahead of current time.
Bulletin_A_Y_Pole = [ Line[ 37: 46].strip() for Line in FD ] # "
Bulletin_A_UTC1_UTC = [ Line[ 58: 68].strip() for Line in FD ] # "
Bulletin_A_dPsi = [ Line[ 97:106].strip() for Line in FD ] # Published about 3 months ahead of current time.
Bulletin_A_dEps = [ Line[116:125].strip() for Line in FD ] # "
# And then the Bulletin B solutions.
Bulletin_B_X_Pole = [ Line[134:144].strip() for Line in FD ] # Should be populated a most a month behind the current time.
Bulletin_B_Y_Pole = [ Line[144:154].strip() for Line in FD ] # "
Bulletin_B_UTC1_UTC = [ Line[154:165].strip() for Line in FD ] # "
Bulletin_B_dPsi = [ Line[165:175].strip() for Line in FD ] # "
Bulletin_B_dEps = [ Line[175:185].strip() for Line in FD ] # "
# And then combine them, prefering the more accurate Bulletin-B solutions.
X_Pole = [ Bulletin_B_X_Pole [i] if Bulletin_B_X_Pole [i] != "" else Bulletin_A_X_Pole [i] for i in range(len(FD)) ]
Y_Pole = [ Bulletin_B_Y_Pole [i] if Bulletin_B_Y_Pole [i] != "" else Bulletin_A_Y_Pole [i] for i in range(len(FD)) ]
UTC1_UTC = [ Bulletin_B_UTC1_UTC[i] if Bulletin_B_UTC1_UTC[i] != "" else Bulletin_A_UTC1_UTC[i] for i in range(len(FD)) ]
dPsi = [ Bulletin_B_dPsi [i] if Bulletin_B_dPsi [i] != "" else Bulletin_A_dPsi [i] for i in range(len(FD)) ]
dEps = [ Bulletin_B_dEps [i] if Bulletin_B_dEps [i] != "" else Bulletin_A_dEps [i] for i in range(len(FD)) ]
# And then cast them all to floats!
try:
IERS_Date = MJD_to_DT( [ float(mjd) for mjd in MJD] )
IERS_MET = DT_to_MET( IERS_Date )
IERS_X_Pole = numpy.array( X_Pole, dtype=float) # Arcseconds.
IERS_Y_Pole = numpy.array( Y_Pole, dtype=float) # Arcseconds.
IERS_UT1_Minus_UTC = numpy.array( UTC1_UTC, dtype=float) # Seconds. (TIME)
IERS_Delta_Psi = numpy.array( dPsi, dtype=float) / 1000 # Milliarcseconds -> Arcseconds.
IERS_Delta_Eplison = numpy.array( dEps, dtype=float) / 1000 # Milliarcseconds -> Arcseconds.
except Exception as E:
raise ValueError("Could not parse the IERS file given! {:s}".format( iers_file_path) ) from E
##############################
# Coefficient-Fetching functions:
def Fetch_Delta_Symbols_MET( MET_List ):
'''
Returns the Delta_Psi and Delta_Epsilon values, linearly interploated.
Values returned as two NP arrays, in RADIAN units.
Args:
MET_List: np array of floats
Returns:
delta_Psi : numpy.array, float. Radians.
delta_Epsilon: numpy.array, float. Radians.
'''
if len(IERS_Date) == 0: Read_IERS_Data() # Load the globals, if needed.
# Just apply a linear interpolation.
delta_Psi = numpy.interp( MET_List, IERS_MET, IERS_Delta_Psi )
delta_Epsilon = numpy.interp( MET_List, IERS_MET, IERS_Delta_Eplison )
# These are given by the IERS & stored as arcseconds, so translate to radians.
delta_Psi = delta_Psi / 3600 * numpy.pi / 180
delta_Epsilon = delta_Epsilon / 3600 * numpy.pi / 180
return delta_Psi, delta_Epsilon
def Fetch_Pole_Wander_MET( MET_List ):
'''
Returns the X- and Y-Wander values, linearly interploated.
Values returned as two NP arrays, in RADIAN units.
Args:
MET_List: np array of floats
Returns:
X_Wander_Rad: numpy.array, float. Radians.
Y_Wander_Rad: numpy.array, float. Radians.
'''
if len(IERS_Date) == 0: Read_IERS_Data() # Load the globals, if needed.
# For now, just apply a linear interpolation.
X_Wander_Arcsec = numpy.interp( MET_List, IERS_MET, IERS_X_Pole )
Y_Wander_Arcsec = numpy.interp( MET_List, IERS_MET, IERS_Y_Pole )
# These are given and stored as arcseconds, so translate to radians.
X_Wander_Rad = X_Wander_Arcsec / 3600 * numpy.pi / 180
Y_Wander_Rad = Y_Wander_Arcsec / 3600 * numpy.pi / 180
return X_Wander_Rad, Y_Wander_Rad
def Fetch_UTC_UT1_Correction_MET( MET_List):
'''
Returns the UTC-UT1 values, linearly interploated.
Values returned as an NP arrays, in seconds, expressed as floats.
Args:
MET_List: np array of floats
Returns:
UT1_Minus_UTC_Seconds: numpy.array, float. Seconds.
'''
if len(IERS_Date) == 0: Read_IERS_Data() # Load the globals, if needed.
# Just apply a linear interpolation.
UT1_Minus_UTC_Seconds = numpy.interp( MET_List, IERS_MET, IERS_UT1_Minus_UTC )
# Can just return. They're floating-point seconds.
return UT1_Minus_UTC_Seconds
##############################
# Time Calculations.
def Calculate_Julian_Day( MET_List ):
'''
Returns the Julian Day numbers of the given MET list.
The given timestamps do not have to be UTC-keyed, like MET.
(Eg: TDT, TDT, TAI)
TLDR: Returns solar days since Noon Jan 1 4713 BC
J2K Epoch = 2000-01-01 12:00 UTC = MET 997444800 = 2,451,545.0 JD
MET Epoch = 1968-05-24 00:00 UTC = MET 0 = 2,440,000.5 JD
Args:
MET_List: np array of floats
Returns:
julian_day_list: np array of floats
'''
# JD at MET 0 == 2440000.5
# Thus just add the number of days (MET/86400) to the MET(0) JD value.
julian_day_list = 2440000.5 + MET_List / 86400
return julian_day_list
def Calculate_TDB( MET_UTC_List ):
'''
Derives the Barycentric Dynamical Time from UTC-MET.
TDB is the most accurate Time system for astronomical systems, as it accounts for t
Returns a NP array of float timestamps, expressed as seconds since the MET epoch.
Args:
MET_UTC_List: np array of floats. Seconds since the MET epoch. Must be UTC-keyed.
Returns:
MET_TDB_List: np array of floats. Seconds since the MET epoch. TDB-keyed
'''
# TDB = TDT + 0.001658(sec) * sin(g+0.10167 * sin(g))
# g = 2*pi*( 357.528* + 35999.50* * T) / 360*
# T = ( JD(TDT) - 3451545 ) / 36525) [Can use JD(TDB) to invert without significant error]
# TDT = TAI + 32.184s
# TAI = UTC + Leap_Seconds + 10
def Get_Leap_Seconds( MET_UTC_List ):
# Only considers leap-seconds after launch.
if numpy.min(MET_UTC_List) <= 1431302400: raise ValueError("Using MET before Launch.") # Sanity check.
Leap_Seconds = numpy.full(MET_UTC_List.shape, 25) # At launch, UTC-TAI was 35 (25+10)
Leap_Seconds[ MET_UTC_List >= 1518048000 ] += 1 # 2016-07-01
Leap_Seconds[ MET_UTC_List >= 1533945600 ] += 1 # 2017-01-01
return Leap_Seconds
MET = MET_UTC_List
TAI = MET + Get_Leap_Seconds(MET) + 10
TDT = TAI + 32.184
T = ( Calculate_Julian_Day(TDT) - 2451545) / 36525
g = 2*numpy.pi * (357.528 + 35999.050 * T) / 360 # Rad * Deg/Deg = Rad
_sec = 0.001658 * numpy.sin( g + 0.10167 * numpy.sin(g))
MET_TDB_List = TDT + _sec
return MET_TDB_List
def Calculate_Julian_Centures_Since_J2K( MET ):
'''
Given MET keyed to UTC, return Julian Centuries since the J2K epoch in the TDB timekey.
Julian Century = 36525 Julian days. (100x365 days)
J2K epoch:
2000-01-01 12:00:00.000 TT = 2451545.0 JD TT (Terrestrial Time)
2000-01-01 11:59:27.816 TAI
2000-01-01 11:58:55.816 UTC
Terrestrial Time = TAI + 32.184
Barycentric Dynamical Time = TAI + 32.184 +/- a few milliseconds
Calculations below use TDB, as following documentation.
Args:
MET: np array of floats. Seconds since the MET epoch. Must be UTC-keyed.
Returns:
JC: np array of floats. Julian centuries since the J2k epoch. TDB-keyed
'''
TDB = Calculate_TDB ( MET )
JD = Calculate_Julian_Day( TDB )
return (JD - 2451545.0) / 36525
def Calculate_Seconds_Into_Day( MET_UT1 ):
'''
Given MET keyed to UT1, return the current rotation fraction of the day.
As MET bound on the midnight, any MET_UT1 % 86400 = 0, which is midnight.
Just returns MET_UT1 % 86400
Args:
MET_UT1: np array of floats. Seconds since the MET epoch. Must be UT1-keyed.
Returns:
Seconds_into_day: np array of floats.
'''
return MET_UT1 % 86400
##############################
# Quaternion Calculations.
# The nominal R1/R2/R3 Roll/Pitch/Yaw rotation matricies.
# Returns Quaternions for each Radian entry given.
# Note that while given Radians, for consistancy, must translate to Degrees for the Quaternion Library.
# Also, the R1/R2/R3 values are sign flipped vs the Quaternion Library's requirement.
# These are kept sign-flipped in input as to follow the documentation for the R1/2/3 rotation matricies.
# Thus the Rad->Deg rotation in the function also flips the sign.
def R1_Quat( Rad ): # Roll-Quaternions
_Roll = Rad * 180 / numpy.pi
_Zeros = numpy.full(Rad.shape, 0)
YPR = numpy.column_stack(( _Zeros, _Zeros, _Roll ))
return lib_quat.YPR_to_Quaternion(YPR)
def R2_Quat( Rad ): # Pitch-Quaternions
_Pitch = Rad * 180 / numpy.pi
_Zeros = numpy.full(Rad.shape, 0)
YPR = numpy.column_stack(( _Zeros, _Pitch, _Zeros ))
return lib_quat.YPR_to_Quaternion(YPR)
def R3_Quat( Rad ): # Yaw-Quaternions
_Yaw = Rad * 180 / numpy.pi
_Zeros = numpy.full(Rad.shape, 0)
YPR = numpy.column_stack(( _Yaw, _Zeros, _Zeros ))
return lib_quat.YPR_to_Quaternion(YPR)
#########################################
# Frame Tranlsation Calculations.
# Precession: P
# Returns a Quaternion List with the Processing correction for each MET timestamp.
# Not Depenant on any external values.
def Calculate_Precession_Correction( MET ):
# Calculate Julian Centuries since J2K. (Float)
T = Calculate_Julian_Centures_Since_J2K( MET )
z = (2306.2181/3600 * T + 1.09468/3600 * T**2 + 0.018203/3600 * T**3) * numpy.pi/180 # Radians.
v = (2004.3109/3600 * T - 0.42665/3600 * T**2 - 0.041833/3600 * T**3) * numpy.pi/180
h = (2306.2181/3600 * T + 0.30188/3600 * T**2 + 0.017998/3600 * T**3) * numpy.pi/180
# h = Zeta (Greek).
# Returns Arrays of Quaternion.
qRot_1 = R3_Quat( -z )
qRot_2 = R2_Quat( v )
qRot_3 = R3_Quat( -h )
return lib_quat._qChain_Multiply((qRot_1, qRot_2, qRot_3)) # Returns Quaternions
# Nutation: N
# Return a Quaternion List for each MET-Entry given.
# Dependant on the IERS data file for the Delta Psi and Epison values.
def Calculate_Nutation_Correction( MET ):
# Calculate Julian Centuries since J2K. (Float)
T = Calculate_Julian_Centures_Since_J2K( MET )
Epsilon = (23 + 26/60 + 21.488/3600) - 46.8150/3600 * T - 0.00059/3600 * T**2 + 0.001813/3600 * T**3
Epsilon = Epsilon * numpy.pi / 180 # Translate to Radians.
delta_Psi, delta_Epsilon = Fetch_Delta_Symbols_MET( MET )
qRot_1 = R1_Quat( -(Epsilon + delta_Epsilon))
qRot_2 = R3_Quat( -delta_Psi )
qRot_3 = R1_Quat( Epsilon )
return lib_quat._qChain_Multiply((qRot_1, qRot_2, qRot_3)) # Returns Quaternions
# Diurnal Matrix: R_S
# Returns a Single-Axis quaternion, along +Z.
def Calculate_Diurnal_Correction( MET , Nutation_Quaternions = None ):
# Derive the UT1 time. It's still expressed as MET, to keep that in it's name.
UT1_MET = MET + Fetch_UTC_UT1_Correction_MET( MET )
# Then calculate the Julian centures since the J2K epoch in UT1 (By definition 2000-01-01 12h UT1 (Noon))
# Do not use Calc_Jul_Cent, as it translates UTC to TDB first.
T_u = (Calculate_Julian_Day( UT1_MET ) - 2451545.0 ) / 36525.0
# The Alpha-Epsilon is extracted from the Nutation Matrix:
# If not given, calculate it.
if Nutation_Quaternions is None:
Nutation_Quaternions = Calculate_Nutation_Correction( MET )
# Extract the Rotation matrix, use to extract the hour-angle between the True Equinox of Date and the Mean Equinox of date.
_Nutation_Rotation_Matrix = lib_quat.Quaternion_to_Rotation_Matrix( Nutation_Quaternions )
Alpha_Epsilon = numpy.arctan( _Nutation_Rotation_Matrix[:,0,1] / _Nutation_Rotation_Matrix[:,0,0] ) # Radians
# Then calculate the theta_G_0 value:
# Expressed as Terrestrial Time, as 86400s=360*. So 1h = 15*, 1m = 1/4* = 15', 1s = 1/240* == 15"
# Express theta_G_0 in -SECONDS- (The Time unit, == 15 arcseconds)
theta_G_0 = 6*60*60 + 41*60 + 50.54841 + 8640184.812866 * T_u + 0.093104 * T_u**2 - 6.2e-6 * T_u**3 # in "Seconds"
theta_G_0 = theta_G_0 / (24*60*60) * 360 * numpy.pi / 180 # "Seconds" -> Rotations -> -> Degrees -> Radians.
# Then, apply the UT1 observation time (The TIME component of the UT1 timestamp)
_UT1_Seconds_Into_Day = Calculate_Seconds_Into_Day( UT1_MET ) # Bascially just returns the modulo of 86400.
UT1_Angle_Into_Day = _UT1_Seconds_Into_Day / (24*60*60) * 360 * numpy.pi / 180 # "Seconds" -> Degrees -> Radians
# Change angle range from {0:360} to {-180:180}
UT1_Angle_Into_Day[ UT1_Angle_Into_Day >= numpy.pi ] -= 2*numpy.pi
# The given multiplier-float below has caused nothing but headaches.
# The replacement 1.0 decreases error, and increases consistancy. No clue as to why.
#theta_G = 1.002737909350795 * UT1_Angle_Into_Day + theta_G_0
theta_G = 1.0 * UT1_Angle_Into_Day + theta_G_0
#theta_G = (2-1.002737909350795) * UT1_Angle_Into_Day + theta_G_0
Theta_G = theta_G + Alpha_Epsilon # Theta = Capital theta. Can't exactly use greek symbols in python.
return R3_Quat( Theta_G )
# Polar Motion: R_M
def Calculate_Polar_Motion_Correction( MET ):
X_Wander, Y_Wander = Fetch_Pole_Wander_MET( MET ) # Radians
qX_Wander = R2_Quat( -X_Wander )
qY_Wander = R1_Quat( -Y_Wander )
return lib_quat._qMultiply(qX_Wander, qY_Wander)
#######################################################
#######################################################
def Build_ICRF_to_ITRF_Quaternions( MET_Timestamps ):
'''
Given a list of timestamps in MET,
Return the J2K->GEO == ICRF->ITRF quaternions at those timestamps.
Note that these are RIGHT-HANDED quaternions.
s.t. vITRF = qICRF->ITRF* . vICRF . qICRF->ITRF
'''
# r_CEP = N . P . r_ICRF
# r_ITRF = R_M . R_S . r_CEP
# Thus
# r_ITRF = R_M . R_S . ( N . P . r_ICRF )
# r_ITRF = R_M . R_S . N . P . r_ICRF
# Thus
# rITRF<-ICRF = R_M . R_S . N . P # Using Left-Handed Rotation-matrix chaining.
# But as building and returning Quaternions, rotation is RIGHT-handed.
# s.t.
# qICRF->ITRF = qP . qN . qR_S . qR_M
# vITRF = qICRF->ITRF* . vICRF . qICRF->ITRF
MET_Timestamps = numpy.array(MET_Timestamps, dtype=float)
# Apply sanity checks first:
Sanity_Check_IERS_Data( MET_Timestamps )
P = Calculate_Precession_Correction ( MET_Timestamps ) # ICRF->P
N = Calculate_Nutation_Correction ( MET_Timestamps ) # P ->N
R_S = Calculate_Diurnal_Correction ( MET_Timestamps ) # N ->R_S
R_M = Calculate_Polar_Motion_Correction( MET_Timestamps ) # R_S ->ITRF
qJ2K_to_GEO = qICRF_to_ITRF = lib_quat._qChain_Multiply(( P, N, R_S, R_M ))
return qICRF_to_ITRF
# Sanity for those who don't want to deal with quaternions.
# Returns a Nx3x3 NP array, where N is the given MET array length.
def Build_ITRF_from_ICRF_Rotation_Matricies( MET_Timestamps ):
'''
Given a list of timestamps in MET,
Return the rotation matrix of GEO<-J2K == ITRF<-ICRF at those timestamps.
Note that these are LEFT-HANDED rotation matricies.
s.t. vITRF = rITRF<-ICRF . vICRF
'''
qICRF_to_ITRF = Build_ICRF_to_ITRF_Quaternions( MET_Timestamps )
rITRF_from_ICRF = lib_quat.Quaternion_to_Rotation_Matrix( qICRF_to_ITRF ) # Change 'to/from' notation to re-inforce handedness change.
return rITRF_from_ICRF
############################################################
############################################################
# Sane transforms
# To rotate a Nx3x3 matrix-list by a Nx3 vector list in Numpy is a nuciance.
# Here the A->B quaternion is fetched, and the lib_quaternion.Vector_Rotate function used in it's stead
def GEO_to_J2K( MET, geo_x, geo_y, geo_z ):
'''
Given a list of times & GEO coordinates, returns the equivalent J2K coordinates.
Note: KNOWN ERROR. <50" at max, usually <10".
Args:
MET: nd.array of floats. MET timestamps.
geo_x: nd.array of floats. X-coordinates in GEO. (Arbitrary units)
geo_y: nd.array of floats. Y-coordinates in GEO. (Arbitrary units)
geo_z: nd.array of floats. Z-coordinates in GEO. (Arbitrary units)
Returns:
j2k_x: nd.array of floats. X-coordinates in J2K. (Same units as given in args)
j2k_y: nd.array of floats. Y-coordinates in J2K. (Same units as given in args)
j2k_z: nd.array of floats. Z-coordinates in J2K. (Same units as given in args)
'''
qJ2K_to_GEO = Build_ICRF_to_ITRF_Quaternions( MET ) # The rotation of J2K->GEO, qGEO<-J2k
qGEO_to_J2K = lib_quat.Conjugate( qJ2K_to_GEO ) # Conjugate = Quaternion inversion/transpose.
# Vector_Rotate takes an Nx3 matrix, so stack the Nx1 position columns:
geo_vector = numpy.column_stack(( geo_x, geo_y, geo_z ))
j2k_vector = lib_quat.Vector_Rotate( qGEO_to_J2K, geo_vector )
# And then split to the individual columns.
j2k_x, j2k_y, j2k_z = j2k_vector.T
return j2k_x, j2k_y, j2k_z
def J2K_to_GEO( MET, j2k_x, j2k_y, j2k_z ):
'''
Given a list of times & J2K coordinates, returns the equivalent GEO coordinates.
Note: KNOWN ERROR. <50" at max, usually <10".
Args:
MET: nd.array of floats. MET timestamps.
j2k_x: nd.array of floats. X-coordinates in J2K. (Arbitrary units)
j2k_y: nd.array of floats. Y-coordinates in J2K. (Arbitrary units)
j2k_z: nd.array of floats. Z-coordinates in J2K. (Arbitrary units)
Returns:
geo_x: nd.array of floats. X-coordinates in GEO. (Same units as given in args)
geo_y: nd.array of floats. Y-coordinates in GEO. (Same units as given in args)
geo_z: nd.array of floats. Z-coordinates in GEO. (Same units as given in args)
'''
qJ2K_to_GEO = Build_ICRF_to_ITRF_Quaternions( MET )
# Vector_Rotate takes an Nx3 matrix, so stack the Nx1 position columns:
j2k_vector = numpy.column_stack(( j2k_x, j2k_y, j2k_z ))
geo_vector = lib_quat.Vector_Rotate( qJ2K_to_GEO, j2k_vector )
# And then split to the individual columns.
geo_x, geo_y, geo_z = geo_vector.T
return geo_x, geo_y, geo_z
############################################################
############################################################
def _Test_Known_Ephemeris(J2K_Ephemeris, GEO_Ephemeris):
'''
Test harness for the J2K<->GEO rotation matricies.
Given GEO & J2K Ephemeris data over an equal region,
Derive the rotation, check accuracy.
Args
J2K_Ephemeris - A structured array from ephemeris.get_j2k_ephemeris()
GEO_Ephemeris - A structured array from ephemeris.get_geo_ephemeris()
Returns
None - Prints relevant data.
'''
# Make sure data is equivocable:
if len(J2K_Ephemeris) != len(GEO_Ephemeris):
raise ValueError("Given unequal lengths of ephemeris data!")
if not numpy.all( J2K_Ephemeris['met'] == GEO_Ephemeris['met'] ):
raise ValueError("Given timestamps lists are not equal!")
# Extract only the MET & Positional data.
MET = J2K_Ephemeris['met']
J2K_Pos = numpy.column_stack(( J2K_Ephemeris['j2k_x'], J2K_Ephemeris['j2k_y'], J2K_Ephemeris['j2k_z'] ))
GEO_Pos = numpy.column_stack(( GEO_Ephemeris['geo_x'], GEO_Ephemeris['geo_y'], GEO_Ephemeris['geo_z'] ))
# Fetch the rotation matricies.
qGEO_from_J2K = Build_ICRF_to_ITRF_Quaternions( MET )
rGEO_from_J2K = Build_ITRF_from_ICRF_Rotation_Matricies( MET )
# Rotate the J2K into the GEO frame:
GEO_Pos_Derived = lib_quat.Vector_Rotate( qGEO_from_J2K, J2K_Pos )
# Calculate the total position error:
Position_Error = numpy.sqrt(numpy.sum((GEO_Pos - GEO_Pos_Derived)**2, axis=1))
# Derive the Normalized GEO vectors to simplify angle-calculation
Norm_GEO, Norm_GEO_Derived = lib_quat.Normalize(GEO_Pos), lib_quat.Normalize(GEO_Pos_Derived)
# Angle = arccos( v1.v2 )
Angle_Error = numpy.arccos( [numpy.dot(Norm_GEO[i], Norm_GEO_Derived[i]) for i in range(len(Norm_GEO))] ) * 180/numpy.pi
print("Stats")
print("Avg Divergence: {:.2f} meters, {:.2f} arcseconds".format( numpy.mean(Position_Error), numpy.mean( Angle_Error) * 3600 ))
print("Min Divergence: {:.2f} meters, {:.2f} arcseconds".format( numpy.min (Position_Error), numpy.min ( Angle_Error) * 3600 ))
print("Max Divergence: {:.2f} meters, {:.2f} arcseconds".format( numpy.max (Position_Error), numpy.max ( Angle_Error) * 3600 ))
if __name__ == "__main__":
print("! lib_j2k_geo_rotation has no CLI interface. Use by importing.")
print("! Primary interface is Build_ITRF_from_ICRF_Quaternions and Build_ITRF_from_ICRF_Rotation_Matricies")
print("! Exiting")
sys.exit(1)
Functions
def Build_ICRF_to_ITRF_Quaternions(MET_Timestamps)
-
Given a list of timestamps in MET, Return the J2K->GEO == ICRF->ITRF quaternions at those timestamps. Note that these are RIGHT-HANDED quaternions. s.t. vITRF = qICRF->ITRF* . vICRF . qICRF->ITRF
Expand source code
def Build_ICRF_to_ITRF_Quaternions( MET_Timestamps ): ''' Given a list of timestamps in MET, Return the J2K->GEO == ICRF->ITRF quaternions at those timestamps. Note that these are RIGHT-HANDED quaternions. s.t. vITRF = qICRF->ITRF* . vICRF . qICRF->ITRF ''' # r_CEP = N . P . r_ICRF # r_ITRF = R_M . R_S . r_CEP # Thus # r_ITRF = R_M . R_S . ( N . P . r_ICRF ) # r_ITRF = R_M . R_S . N . P . r_ICRF # Thus # rITRF<-ICRF = R_M . R_S . N . P # Using Left-Handed Rotation-matrix chaining. # But as building and returning Quaternions, rotation is RIGHT-handed. # s.t. # qICRF->ITRF = qP . qN . qR_S . qR_M # vITRF = qICRF->ITRF* . vICRF . qICRF->ITRF MET_Timestamps = numpy.array(MET_Timestamps, dtype=float) # Apply sanity checks first: Sanity_Check_IERS_Data( MET_Timestamps ) P = Calculate_Precession_Correction ( MET_Timestamps ) # ICRF->P N = Calculate_Nutation_Correction ( MET_Timestamps ) # P ->N R_S = Calculate_Diurnal_Correction ( MET_Timestamps ) # N ->R_S R_M = Calculate_Polar_Motion_Correction( MET_Timestamps ) # R_S ->ITRF qJ2K_to_GEO = qICRF_to_ITRF = lib_quat._qChain_Multiply(( P, N, R_S, R_M )) return qICRF_to_ITRF
def Build_ITRF_from_ICRF_Rotation_Matricies(MET_Timestamps)
-
Given a list of timestamps in MET, Return the rotation matrix of GEO<-J2K == ITRF<-ICRF at those timestamps. Note that these are LEFT-HANDED rotation matricies. s.t. vITRF = rITRF<-ICRF . vICRF
Expand source code
def Build_ITRF_from_ICRF_Rotation_Matricies( MET_Timestamps ): ''' Given a list of timestamps in MET, Return the rotation matrix of GEO<-J2K == ITRF<-ICRF at those timestamps. Note that these are LEFT-HANDED rotation matricies. s.t. vITRF = rITRF<-ICRF . vICRF ''' qICRF_to_ITRF = Build_ICRF_to_ITRF_Quaternions( MET_Timestamps ) rITRF_from_ICRF = lib_quat.Quaternion_to_Rotation_Matrix( qICRF_to_ITRF ) # Change 'to/from' notation to re-inforce handedness change. return rITRF_from_ICRF
def Calculate_Diurnal_Correction(MET, Nutation_Quaternions=None)
-
Expand source code
def Calculate_Diurnal_Correction( MET , Nutation_Quaternions = None ): # Derive the UT1 time. It's still expressed as MET, to keep that in it's name. UT1_MET = MET + Fetch_UTC_UT1_Correction_MET( MET ) # Then calculate the Julian centures since the J2K epoch in UT1 (By definition 2000-01-01 12h UT1 (Noon)) # Do not use Calc_Jul_Cent, as it translates UTC to TDB first. T_u = (Calculate_Julian_Day( UT1_MET ) - 2451545.0 ) / 36525.0 # The Alpha-Epsilon is extracted from the Nutation Matrix: # If not given, calculate it. if Nutation_Quaternions is None: Nutation_Quaternions = Calculate_Nutation_Correction( MET ) # Extract the Rotation matrix, use to extract the hour-angle between the True Equinox of Date and the Mean Equinox of date. _Nutation_Rotation_Matrix = lib_quat.Quaternion_to_Rotation_Matrix( Nutation_Quaternions ) Alpha_Epsilon = numpy.arctan( _Nutation_Rotation_Matrix[:,0,1] / _Nutation_Rotation_Matrix[:,0,0] ) # Radians # Then calculate the theta_G_0 value: # Expressed as Terrestrial Time, as 86400s=360*. So 1h = 15*, 1m = 1/4* = 15', 1s = 1/240* == 15" # Express theta_G_0 in -SECONDS- (The Time unit, == 15 arcseconds) theta_G_0 = 6*60*60 + 41*60 + 50.54841 + 8640184.812866 * T_u + 0.093104 * T_u**2 - 6.2e-6 * T_u**3 # in "Seconds" theta_G_0 = theta_G_0 / (24*60*60) * 360 * numpy.pi / 180 # "Seconds" -> Rotations -> -> Degrees -> Radians. # Then, apply the UT1 observation time (The TIME component of the UT1 timestamp) _UT1_Seconds_Into_Day = Calculate_Seconds_Into_Day( UT1_MET ) # Bascially just returns the modulo of 86400. UT1_Angle_Into_Day = _UT1_Seconds_Into_Day / (24*60*60) * 360 * numpy.pi / 180 # "Seconds" -> Degrees -> Radians # Change angle range from {0:360} to {-180:180} UT1_Angle_Into_Day[ UT1_Angle_Into_Day >= numpy.pi ] -= 2*numpy.pi # The given multiplier-float below has caused nothing but headaches. # The replacement 1.0 decreases error, and increases consistancy. No clue as to why. #theta_G = 1.002737909350795 * UT1_Angle_Into_Day + theta_G_0 theta_G = 1.0 * UT1_Angle_Into_Day + theta_G_0 #theta_G = (2-1.002737909350795) * UT1_Angle_Into_Day + theta_G_0 Theta_G = theta_G + Alpha_Epsilon # Theta = Capital theta. Can't exactly use greek symbols in python. return R3_Quat( Theta_G )
def Calculate_Julian_Centures_Since_J2K(MET)
-
Given MET keyed to UTC, return Julian Centuries since the J2K epoch in the TDB timekey.
Julian Century = 36525 Julian days. (100x365 days) J2K epoch: 2000-01-01 12:00:00.000 TT = 2451545.0 JD TT (Terrestrial Time) 2000-01-01 11:59:27.816 TAI 2000-01-01 11:58:55.816 UTC
Terrestrial Time = TAI + 32.184 Barycentric Dynamical Time = TAI + 32.184 +/- a few milliseconds
Calculations below use TDB, as following documentation.
Args
MET
- np array of floats. Seconds since the MET epoch. Must be UTC-keyed.
Returns
JC
- np array of floats. Julian centuries since the J2k epoch. TDB-keyed
Expand source code
def Calculate_Julian_Centures_Since_J2K( MET ): ''' Given MET keyed to UTC, return Julian Centuries since the J2K epoch in the TDB timekey. Julian Century = 36525 Julian days. (100x365 days) J2K epoch: 2000-01-01 12:00:00.000 TT = 2451545.0 JD TT (Terrestrial Time) 2000-01-01 11:59:27.816 TAI 2000-01-01 11:58:55.816 UTC Terrestrial Time = TAI + 32.184 Barycentric Dynamical Time = TAI + 32.184 +/- a few milliseconds Calculations below use TDB, as following documentation. Args: MET: np array of floats. Seconds since the MET epoch. Must be UTC-keyed. Returns: JC: np array of floats. Julian centuries since the J2k epoch. TDB-keyed ''' TDB = Calculate_TDB ( MET ) JD = Calculate_Julian_Day( TDB ) return (JD - 2451545.0) / 36525
def Calculate_Julian_Day(MET_List)
-
Returns the Julian Day numbers of the given MET list. The given timestamps do not have to be UTC-keyed, like MET. (Eg: TDT, TDT, TAI)
TLDR: Returns solar days since Noon Jan 1 4713 BC J2K Epoch = 2000-01-01 12:00 UTC = MET 997444800 = 2,451,545.0 JD MET Epoch = 1968-05-24 00:00 UTC = MET 0 = 2,440,000.5 JD
Args
MET_List
- np array of floats
Returns
julian_day_list
- np array of floats
Expand source code
def Calculate_Julian_Day( MET_List ): ''' Returns the Julian Day numbers of the given MET list. The given timestamps do not have to be UTC-keyed, like MET. (Eg: TDT, TDT, TAI) TLDR: Returns solar days since Noon Jan 1 4713 BC J2K Epoch = 2000-01-01 12:00 UTC = MET 997444800 = 2,451,545.0 JD MET Epoch = 1968-05-24 00:00 UTC = MET 0 = 2,440,000.5 JD Args: MET_List: np array of floats Returns: julian_day_list: np array of floats ''' # JD at MET 0 == 2440000.5 # Thus just add the number of days (MET/86400) to the MET(0) JD value. julian_day_list = 2440000.5 + MET_List / 86400 return julian_day_list
def Calculate_Nutation_Correction(MET)
-
Expand source code
def Calculate_Nutation_Correction( MET ): # Calculate Julian Centuries since J2K. (Float) T = Calculate_Julian_Centures_Since_J2K( MET ) Epsilon = (23 + 26/60 + 21.488/3600) - 46.8150/3600 * T - 0.00059/3600 * T**2 + 0.001813/3600 * T**3 Epsilon = Epsilon * numpy.pi / 180 # Translate to Radians. delta_Psi, delta_Epsilon = Fetch_Delta_Symbols_MET( MET ) qRot_1 = R1_Quat( -(Epsilon + delta_Epsilon)) qRot_2 = R3_Quat( -delta_Psi ) qRot_3 = R1_Quat( Epsilon ) return lib_quat._qChain_Multiply((qRot_1, qRot_2, qRot_3)) # Returns Quaternions
def Calculate_Polar_Motion_Correction(MET)
-
Expand source code
def Calculate_Polar_Motion_Correction( MET ): X_Wander, Y_Wander = Fetch_Pole_Wander_MET( MET ) # Radians qX_Wander = R2_Quat( -X_Wander ) qY_Wander = R1_Quat( -Y_Wander ) return lib_quat._qMultiply(qX_Wander, qY_Wander)
def Calculate_Precession_Correction(MET)
-
Expand source code
def Calculate_Precession_Correction( MET ): # Calculate Julian Centuries since J2K. (Float) T = Calculate_Julian_Centures_Since_J2K( MET ) z = (2306.2181/3600 * T + 1.09468/3600 * T**2 + 0.018203/3600 * T**3) * numpy.pi/180 # Radians. v = (2004.3109/3600 * T - 0.42665/3600 * T**2 - 0.041833/3600 * T**3) * numpy.pi/180 h = (2306.2181/3600 * T + 0.30188/3600 * T**2 + 0.017998/3600 * T**3) * numpy.pi/180 # h = Zeta (Greek). # Returns Arrays of Quaternion. qRot_1 = R3_Quat( -z ) qRot_2 = R2_Quat( v ) qRot_3 = R3_Quat( -h ) return lib_quat._qChain_Multiply((qRot_1, qRot_2, qRot_3)) # Returns Quaternions
def Calculate_Seconds_Into_Day(MET_UT1)
-
Given MET keyed to UT1, return the current rotation fraction of the day. As MET bound on the midnight, any MET_UT1 % 86400 = 0, which is midnight. Just returns MET_UT1 % 86400
Args
MET_UT1
- np array of floats. Seconds since the MET epoch. Must be UT1-keyed.
Returns
Seconds_into_day
- np array of floats.
Expand source code
def Calculate_Seconds_Into_Day( MET_UT1 ): ''' Given MET keyed to UT1, return the current rotation fraction of the day. As MET bound on the midnight, any MET_UT1 % 86400 = 0, which is midnight. Just returns MET_UT1 % 86400 Args: MET_UT1: np array of floats. Seconds since the MET epoch. Must be UT1-keyed. Returns: Seconds_into_day: np array of floats. ''' return MET_UT1 % 86400
def Calculate_TDB(MET_UTC_List)
-
Derives the Barycentric Dynamical Time from UTC-MET. TDB is the most accurate Time system for astronomical systems, as it accounts for t
Returns a NP array of float timestamps, expressed as seconds since the MET epoch.
Args
MET_UTC_List
- np array of floats. Seconds since the MET epoch. Must be UTC-keyed.
Returns
MET_TDB_List
- np array of floats. Seconds since the MET epoch. TDB-keyed
Expand source code
def Calculate_TDB( MET_UTC_List ): ''' Derives the Barycentric Dynamical Time from UTC-MET. TDB is the most accurate Time system for astronomical systems, as it accounts for t Returns a NP array of float timestamps, expressed as seconds since the MET epoch. Args: MET_UTC_List: np array of floats. Seconds since the MET epoch. Must be UTC-keyed. Returns: MET_TDB_List: np array of floats. Seconds since the MET epoch. TDB-keyed ''' # TDB = TDT + 0.001658(sec) * sin(g+0.10167 * sin(g)) # g = 2*pi*( 357.528* + 35999.50* * T) / 360* # T = ( JD(TDT) - 3451545 ) / 36525) [Can use JD(TDB) to invert without significant error] # TDT = TAI + 32.184s # TAI = UTC + Leap_Seconds + 10 def Get_Leap_Seconds( MET_UTC_List ): # Only considers leap-seconds after launch. if numpy.min(MET_UTC_List) <= 1431302400: raise ValueError("Using MET before Launch.") # Sanity check. Leap_Seconds = numpy.full(MET_UTC_List.shape, 25) # At launch, UTC-TAI was 35 (25+10) Leap_Seconds[ MET_UTC_List >= 1518048000 ] += 1 # 2016-07-01 Leap_Seconds[ MET_UTC_List >= 1533945600 ] += 1 # 2017-01-01 return Leap_Seconds MET = MET_UTC_List TAI = MET + Get_Leap_Seconds(MET) + 10 TDT = TAI + 32.184 T = ( Calculate_Julian_Day(TDT) - 2451545) / 36525 g = 2*numpy.pi * (357.528 + 35999.050 * T) / 360 # Rad * Deg/Deg = Rad _sec = 0.001658 * numpy.sin( g + 0.10167 * numpy.sin(g)) MET_TDB_List = TDT + _sec return MET_TDB_List
def DT_to_MET(DT)
-
Expand source code
def DT_to_MET( DT ): # Doesn't -have- to be UTC, but official MET follows UTC, like Unix-time. return numpy.array([ ( dt - datetime.datetime(1968,5,24)).total_seconds() for dt in DT ])
def DT_to_MJD(DT)
-
Expand source code
def DT_to_MJD( DT ): return numpy.array([ ( dt - datetime.datetime(1858,11,17)).total_seconds() // 86400 for dt in DT ])
def Fetch_Delta_Symbols_MET(MET_List)
-
Returns the Delta_Psi and Delta_Epsilon values, linearly interploated. Values returned as two NP arrays, in RADIAN units.
Args
MET_List
- np array of floats
Returns: delta_Psi : numpy.array, float. Radians. delta_Epsilon: numpy.array, float. Radians.
Expand source code
def Fetch_Delta_Symbols_MET( MET_List ): ''' Returns the Delta_Psi and Delta_Epsilon values, linearly interploated. Values returned as two NP arrays, in RADIAN units. Args: MET_List: np array of floats Returns: delta_Psi : numpy.array, float. Radians. delta_Epsilon: numpy.array, float. Radians. ''' if len(IERS_Date) == 0: Read_IERS_Data() # Load the globals, if needed. # Just apply a linear interpolation. delta_Psi = numpy.interp( MET_List, IERS_MET, IERS_Delta_Psi ) delta_Epsilon = numpy.interp( MET_List, IERS_MET, IERS_Delta_Eplison ) # These are given by the IERS & stored as arcseconds, so translate to radians. delta_Psi = delta_Psi / 3600 * numpy.pi / 180 delta_Epsilon = delta_Epsilon / 3600 * numpy.pi / 180 return delta_Psi, delta_Epsilon
def Fetch_Pole_Wander_MET(MET_List)
-
Returns the X- and Y-Wander values, linearly interploated. Values returned as two NP arrays, in RADIAN units.
Args
MET_List
- np array of floats
Returns: X_Wander_Rad: numpy.array, float. Radians. Y_Wander_Rad: numpy.array, float. Radians.
Expand source code
def Fetch_Pole_Wander_MET( MET_List ): ''' Returns the X- and Y-Wander values, linearly interploated. Values returned as two NP arrays, in RADIAN units. Args: MET_List: np array of floats Returns: X_Wander_Rad: numpy.array, float. Radians. Y_Wander_Rad: numpy.array, float. Radians. ''' if len(IERS_Date) == 0: Read_IERS_Data() # Load the globals, if needed. # For now, just apply a linear interpolation. X_Wander_Arcsec = numpy.interp( MET_List, IERS_MET, IERS_X_Pole ) Y_Wander_Arcsec = numpy.interp( MET_List, IERS_MET, IERS_Y_Pole ) # These are given and stored as arcseconds, so translate to radians. X_Wander_Rad = X_Wander_Arcsec / 3600 * numpy.pi / 180 Y_Wander_Rad = Y_Wander_Arcsec / 3600 * numpy.pi / 180 return X_Wander_Rad, Y_Wander_Rad
def Fetch_UTC_UT1_Correction_MET(MET_List)
-
Returns the UTC-UT1 values, linearly interploated. Values returned as an NP arrays, in seconds, expressed as floats.
Args
MET_List
- np array of floats
Returns: UT1_Minus_UTC_Seconds: numpy.array, float. Seconds.
Expand source code
def Fetch_UTC_UT1_Correction_MET( MET_List): ''' Returns the UTC-UT1 values, linearly interploated. Values returned as an NP arrays, in seconds, expressed as floats. Args: MET_List: np array of floats Returns: UT1_Minus_UTC_Seconds: numpy.array, float. Seconds. ''' if len(IERS_Date) == 0: Read_IERS_Data() # Load the globals, if needed. # Just apply a linear interpolation. UT1_Minus_UTC_Seconds = numpy.interp( MET_List, IERS_MET, IERS_UT1_Minus_UTC ) # Can just return. They're floating-point seconds. return UT1_Minus_UTC_Seconds
def GEO_to_J2K(MET, geo_x, geo_y, geo_z)
-
Given a list of times & GEO coordinates, returns the equivalent J2K coordinates. Note: KNOWN ERROR. <50" at max, usually <10".
Args
MET
- nd.array of floats. MET timestamps.
geo_x
- nd.array of floats. X-coordinates in GEO. (Arbitrary units)
geo_y
- nd.array of floats. Y-coordinates in GEO. (Arbitrary units)
geo_z
- nd.array of floats. Z-coordinates in GEO. (Arbitrary units)
Returns
j2k_x
- nd.array of floats. X-coordinates in J2K. (Same units as given in args)
j2k_y
- nd.array of floats. Y-coordinates in J2K. (Same units as given in args)
j2k_z
- nd.array of floats. Z-coordinates in J2K. (Same units as given in args)
Expand source code
def GEO_to_J2K( MET, geo_x, geo_y, geo_z ): ''' Given a list of times & GEO coordinates, returns the equivalent J2K coordinates. Note: KNOWN ERROR. <50" at max, usually <10". Args: MET: nd.array of floats. MET timestamps. geo_x: nd.array of floats. X-coordinates in GEO. (Arbitrary units) geo_y: nd.array of floats. Y-coordinates in GEO. (Arbitrary units) geo_z: nd.array of floats. Z-coordinates in GEO. (Arbitrary units) Returns: j2k_x: nd.array of floats. X-coordinates in J2K. (Same units as given in args) j2k_y: nd.array of floats. Y-coordinates in J2K. (Same units as given in args) j2k_z: nd.array of floats. Z-coordinates in J2K. (Same units as given in args) ''' qJ2K_to_GEO = Build_ICRF_to_ITRF_Quaternions( MET ) # The rotation of J2K->GEO, qGEO<-J2k qGEO_to_J2K = lib_quat.Conjugate( qJ2K_to_GEO ) # Conjugate = Quaternion inversion/transpose. # Vector_Rotate takes an Nx3 matrix, so stack the Nx1 position columns: geo_vector = numpy.column_stack(( geo_x, geo_y, geo_z )) j2k_vector = lib_quat.Vector_Rotate( qGEO_to_J2K, geo_vector ) # And then split to the individual columns. j2k_x, j2k_y, j2k_z = j2k_vector.T return j2k_x, j2k_y, j2k_z
def J2K_to_GEO(MET, j2k_x, j2k_y, j2k_z)
-
Given a list of times & J2K coordinates, returns the equivalent GEO coordinates. Note: KNOWN ERROR. <50" at max, usually <10".
Args
MET
- nd.array of floats. MET timestamps.
j2k_x
- nd.array of floats. X-coordinates in J2K. (Arbitrary units)
j2k_y
- nd.array of floats. Y-coordinates in J2K. (Arbitrary units)
j2k_z
- nd.array of floats. Z-coordinates in J2K. (Arbitrary units)
Returns
geo_x
- nd.array of floats. X-coordinates in GEO. (Same units as given in args)
geo_y
- nd.array of floats. Y-coordinates in GEO. (Same units as given in args)
geo_z
- nd.array of floats. Z-coordinates in GEO. (Same units as given in args)
Expand source code
def J2K_to_GEO( MET, j2k_x, j2k_y, j2k_z ): ''' Given a list of times & J2K coordinates, returns the equivalent GEO coordinates. Note: KNOWN ERROR. <50" at max, usually <10". Args: MET: nd.array of floats. MET timestamps. j2k_x: nd.array of floats. X-coordinates in J2K. (Arbitrary units) j2k_y: nd.array of floats. Y-coordinates in J2K. (Arbitrary units) j2k_z: nd.array of floats. Z-coordinates in J2K. (Arbitrary units) Returns: geo_x: nd.array of floats. X-coordinates in GEO. (Same units as given in args) geo_y: nd.array of floats. Y-coordinates in GEO. (Same units as given in args) geo_z: nd.array of floats. Z-coordinates in GEO. (Same units as given in args) ''' qJ2K_to_GEO = Build_ICRF_to_ITRF_Quaternions( MET ) # Vector_Rotate takes an Nx3 matrix, so stack the Nx1 position columns: j2k_vector = numpy.column_stack(( j2k_x, j2k_y, j2k_z )) geo_vector = lib_quat.Vector_Rotate( qJ2K_to_GEO, j2k_vector ) # And then split to the individual columns. geo_x, geo_y, geo_z = geo_vector.T return geo_x, geo_y, geo_z
def MET_to_DT(MET)
-
Expand source code
def MET_to_DT( MET ): return [ datetime.datetime(1968,5,24) + datetime.timedelta(seconds=met) for met in MET]
def MJD_to_DT(MDJ)
-
Expand source code
def MJD_to_DT( MDJ ): # Epoch of 1858-11-17 00:00:00 UTC return [ datetime.datetime(1858,11,17) + datetime.timedelta(days=mjd) for mjd in MDJ ]
def R1_Quat(Rad)
-
Expand source code
def R1_Quat( Rad ): # Roll-Quaternions _Roll = Rad * 180 / numpy.pi _Zeros = numpy.full(Rad.shape, 0) YPR = numpy.column_stack(( _Zeros, _Zeros, _Roll )) return lib_quat.YPR_to_Quaternion(YPR)
def R2_Quat(Rad)
-
Expand source code
def R2_Quat( Rad ): # Pitch-Quaternions _Pitch = Rad * 180 / numpy.pi _Zeros = numpy.full(Rad.shape, 0) YPR = numpy.column_stack(( _Zeros, _Pitch, _Zeros )) return lib_quat.YPR_to_Quaternion(YPR)
def R3_Quat(Rad)
-
Expand source code
def R3_Quat( Rad ): # Yaw-Quaternions _Yaw = Rad * 180 / numpy.pi _Zeros = numpy.full(Rad.shape, 0) YPR = numpy.column_stack(( _Yaw, _Zeros, _Zeros )) return lib_quat.YPR_to_Quaternion(YPR)
def Read_IERS_Data(iers_file_path='/mnt/c/Users/Warren Holley/Desktop/ephemeris_library_3/ephemeris_library/rotation_libraries/iers_eop.txt')
-
Loads the IERS coefficient data into memory.
Args
iers_file_path
- An abolute path to the given IERS file.
Defaults to ./
/iers_eop.txt
Expand source code
def Read_IERS_Data( iers_file_path = IERS_File_Path ): ''' Loads the IERS coefficient data into memory. Args: iers_file_path: An abolute path to the given IERS file. Defaults to ./<Script_Path>/iers_eop.txt ''' # The file format used by IERS is... odd. Read the 'IERS_EOP_Documentation' file next to the EOP file for more details. global IERS_Date, IERS_MET, IERS_X_Pole, IERS_Y_Pole, IERS_UT1_Minus_UTC, IERS_Delta_Psi, IERS_Delta_Eplison FD = open( iers_file_path ).readlines() # Will raise a FileNotFoundError if could not read the file. # Read the MJD as to filter the input solutions. MJD = numpy.array( [Line[7:15].strip() for Line in FD], dtype=float ) # Filter for entries from ~1 week before launch, to 1 week after the current time. Just reduced memory load. Launch_DT, Current_DT = datetime.datetime( 2013, 9, 29 ), datetime.datetime.now().replace( hour=0, minute=0, second=0, microsecond=0) Launch_MJD, Current_MJD = DT_to_MJD( [Launch_DT, Current_DT] ) _Valid_Mask = numpy.logical_and( Launch_MJD -7 <= MJD, MJD <= Current_MJD+7 ) FD = [ FD[i] for i in range(len(FD)) if _Valid_Mask[i] ] # Don't bother with year/month, just read MJD. MJD = [ Line[ 7: 15].strip() for Line in FD ] # Guaranteed populted for each line. # Read the Bulletin A values. Bulletin_A_X_Pole = [ Line[ 18: 27].strip() for Line in FD ] # Published about 6-10 months ahead of current time. Bulletin_A_Y_Pole = [ Line[ 37: 46].strip() for Line in FD ] # " Bulletin_A_UTC1_UTC = [ Line[ 58: 68].strip() for Line in FD ] # " Bulletin_A_dPsi = [ Line[ 97:106].strip() for Line in FD ] # Published about 3 months ahead of current time. Bulletin_A_dEps = [ Line[116:125].strip() for Line in FD ] # " # And then the Bulletin B solutions. Bulletin_B_X_Pole = [ Line[134:144].strip() for Line in FD ] # Should be populated a most a month behind the current time. Bulletin_B_Y_Pole = [ Line[144:154].strip() for Line in FD ] # " Bulletin_B_UTC1_UTC = [ Line[154:165].strip() for Line in FD ] # " Bulletin_B_dPsi = [ Line[165:175].strip() for Line in FD ] # " Bulletin_B_dEps = [ Line[175:185].strip() for Line in FD ] # " # And then combine them, prefering the more accurate Bulletin-B solutions. X_Pole = [ Bulletin_B_X_Pole [i] if Bulletin_B_X_Pole [i] != "" else Bulletin_A_X_Pole [i] for i in range(len(FD)) ] Y_Pole = [ Bulletin_B_Y_Pole [i] if Bulletin_B_Y_Pole [i] != "" else Bulletin_A_Y_Pole [i] for i in range(len(FD)) ] UTC1_UTC = [ Bulletin_B_UTC1_UTC[i] if Bulletin_B_UTC1_UTC[i] != "" else Bulletin_A_UTC1_UTC[i] for i in range(len(FD)) ] dPsi = [ Bulletin_B_dPsi [i] if Bulletin_B_dPsi [i] != "" else Bulletin_A_dPsi [i] for i in range(len(FD)) ] dEps = [ Bulletin_B_dEps [i] if Bulletin_B_dEps [i] != "" else Bulletin_A_dEps [i] for i in range(len(FD)) ] # And then cast them all to floats! try: IERS_Date = MJD_to_DT( [ float(mjd) for mjd in MJD] ) IERS_MET = DT_to_MET( IERS_Date ) IERS_X_Pole = numpy.array( X_Pole, dtype=float) # Arcseconds. IERS_Y_Pole = numpy.array( Y_Pole, dtype=float) # Arcseconds. IERS_UT1_Minus_UTC = numpy.array( UTC1_UTC, dtype=float) # Seconds. (TIME) IERS_Delta_Psi = numpy.array( dPsi, dtype=float) / 1000 # Milliarcseconds -> Arcseconds. IERS_Delta_Eplison = numpy.array( dEps, dtype=float) / 1000 # Milliarcseconds -> Arcseconds. except Exception as E: raise ValueError("Could not parse the IERS file given! {:s}".format( iers_file_path) ) from E
def Sanity_Check_IERS_Data(MET_List)
-
Given a MET list, load & read the IERS EOP data, & make sure that all times are within the given range.
Raises a ValueError should the requested list fall outside the constants file range.
Expand source code
def Sanity_Check_IERS_Data( MET_List ): ''' Given a MET list, load & read the IERS EOP data, & make sure that all times are within the given range. Raises a ValueError should the requested list fall outside the constants file range. ''' if len(IERS_Date) == 0: Read_IERS_Data() # Load the globals, if needed. # If the requested times are outside the range of the IERS dataset, update the file. if numpy.min(MET_List) < numpy.min(IERS_MET) or numpy.max(MET_List) > numpy.max(IERS_MET): raise ValueError("Requested MET outside the range of the IERS_EOP file. Please update with .Update_IERS_Data()")
def Update_IERS_Data()
-
Updates the IERS file found in the lib_j2k_geo_rotation.IERS_File_Path variable. (Default:
/rotation_libraries/iers_eop.txt Documentation for the file found in this library's checkout:
/rotation_libraries/iers_eop_documentation.txt file Raises ConnectionRefusedError on inability to connect/update the file.
Expand source code
def Update_IERS_Data(): ''' Updates the IERS file found in the lib_j2k_geo_rotation.IERS_File_Path variable. (Default: <path>/rotation_libraries/iers_eop.txt Documentation for the file found in this library's checkout: <path>/rotation_libraries/iers_eop_documentation.txt file Raises ConnectionRefusedError on inability to connect/update the file. ''' # Standard Rapid EOP Data since 01. January 1992 (IAU1980) # From https://www.iers.org/IERS/EN/DataProducts/EarthOrientationData/eop.html - finals.data (IAU1980) # Metadata: https://datacenter.iers.org/versionMetadata.php?filename=latestVersionMeta/8_FINALS.DATA_IAU1980_V2013_018.txt IERS_Default_URL = "https://datacenter.iers.org/data/8/finals.data" try: import urllib.request, ssl # The IERS uses a self-signed certificate, which the URL-lib will raise exceptions about. # As such, need to suppress the warning by overriding the default SSL-configuration context. ssl._create_default_https_context = ssl._create_unverified_context # Make the request: URL_Request = urllib.request.urlopen( IERS_Default_URL ) IERS_Data = URL_Request.read().decode('utf-8') # Read & Decode the data # And overwrite the IERS Text File. FH = open( IERS_File_Path,'w' ) FH.write( IERS_Data ) FH.close() # And update the contents! Read_IERS_Data() return except: pass # Don't raise exception here, re-raise below. print("!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!") print("!!! Could not fetch Earth Orientation Data from the IERS servers") print("!!! Please update the IERS file at {:s}".format(IERS_File_Path)) print("!!! Using the EOP 14 C04 (IAU1980) publication.") print("!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!") raise ConnectionRefusedError()