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()