text: # -*- coding: utf-8 -*-
  # Licensed under a 3-clause BSD style license - see LICENSE.rst
  """
  This module contains functions/values used repeatedly in different modules of
  the ``builtin_frames`` package.
  """
  
  import warnings
  
  import erfa
  import numpy as np
  
  from astropy import units as u
  from astropy.time import Time
  from astropy.utils import iers
  from astropy.utils.exceptions import AstropyWarning
  from ..representation import CartesianRepresentation, CartesianDifferential
  
  
  # We use tt as the time scale for this equinoxes, primarily because it is the
  # convention for J2000 (it is unclear if there is any "right answer" for B1950)
  # while #8600 makes this the default behavior, we show it here to ensure it's
  # clear which is used here
  EQUINOX_J2000 = Time('J2000', scale='tt')
  EQUINOX_B1950 = Time('B1950', scale='tt')
  
  # This is a time object that is the default "obstime" when such an attribute is
  # necessary.  Currently, we use J2000.
  DEFAULT_OBSTIME = Time('J2000', scale='tt')
  
  PIOVER2 = np.pi / 2.
  
  # comes from the mean of the 1962-2014 IERS B data
  _DEFAULT_PM = (0.035, 0.29)*u.arcsec
  
  
  def get_polar_motion(time):
      """
      gets the two polar motion components in radians for use with apio13
      """
      # Get the polar motion from the IERS table
      iers_table = iers.earth_orientation_table.get()
      xp, yp, status = iers_table.pm_xy(time, return_status=True)
  
      wmsg = (
          'Tried to get polar motions for times {} IERS data is '
          'valid. Defaulting to polar motion from the 50-yr mean for those. '
          'This may affect precision at the arcsec level'
      )
      if np.any(status == iers.TIME_BEFORE_IERS_RANGE):
          xp[status == iers.TIME_BEFORE_IERS_RANGE] = _DEFAULT_PM[0]
          yp[status == iers.TIME_BEFORE_IERS_RANGE] = _DEFAULT_PM[1]
  
          warnings.warn(wmsg.format('before'), AstropyWarning)
  
      if np.any(status == iers.TIME_BEYOND_IERS_RANGE):
  
          xp[status == iers.TIME_BEYOND_IERS_RANGE] = _DEFAULT_PM[0]
          yp[status == iers.TIME_BEYOND_IERS_RANGE] = _DEFAULT_PM[1]
  
          warnings.warn(wmsg.format('after'), AstropyWarning)
  
      return xp.to_value(u.radian), yp.to_value(u.radian)
  
  
  def _warn_iers(ierserr):
      """
      Generate a warning for an IERSRangeerror
  
      Parameters
      ----------
      ierserr : An `~astropy.utils.iers.IERSRangeError`
      """
      msg = '{0} Assuming UT1-UTC=0 for coordinate transformations.'
      warnings.warn(msg.format(ierserr.args[0]), AstropyWarning)
  
  
  def get_dut1utc(time):
      """
      This function is used to get UT1-UTC in coordinates because normally it
      gives an error outside the IERS range, but in coordinates we want to allow
      it to go through but with a warning.
      """
      try:
          return time.delta_ut1_utc
      except iers.IERSRangeError as e:
          _warn_iers(e)
          return np.zeros(time.shape)
  
  
  def get_jd12(time, scale):
      """
      Gets ``jd1`` and ``jd2`` from a time object in a particular scale.
  
      Parameters
      ----------
      time : `~astropy.time.Time`
          The time to get the jds for
      scale : str
          The time scale to get the jds for
  
      Returns
      -------
      jd1 : float
      jd2 : float
      """
      if time.scale == scale:
          newtime = time
      else:
          try:
              newtime = getattr(time, scale)
          except iers.IERSRangeError as e:
              _warn_iers(e)
              newtime = time
  
      return newtime.jd1, newtime.jd2
  
  
  def norm(p):
      """
      Normalise a p-vector.
      """
      return p / np.sqrt(np.einsum('...i,...i', p, p))[..., np.newaxis]
  
  
  def pav2pv(p, v):
      """
      Combine p- and v- vectors into a pv-vector.
      """
      pv = np.empty(np.broadcast(p, v).shape[:-1], erfa.dt_pv)
      pv['p'] = p
      pv['v'] = v
      return pv
  
  
  def get_cip(jd1, jd2):
      """
      Find the X, Y coordinates of the CIP and the CIO locator, s.
  
      Parameters
      ----------
      jd1 : float or `np.ndarray`
          First part of two part Julian date (TDB)
      jd2 : float or `np.ndarray`
          Second part of two part Julian date (TDB)
  
      Returns
      --------
      x : float or `np.ndarray`
          x coordinate of the CIP
      y : float or `np.ndarray`
          y coordinate of the CIP
      s : float or `np.ndarray`
          CIO locator, s
      """
      # classical NPB matrix, IAU 2006/2000A
      rpnb = erfa.pnm06a(jd1, jd2)
      # CIP X, Y coordinates from array
      x, y = erfa.bpn2xy(rpnb)
      # CIO locator, s
      s = erfa.s06(jd1, jd2, x, y)
      return x, y, s
  
  
  def aticq(srepr, astrom):
      """
      A slightly modified version of the ERFA function ``eraAticq``.
  
      ``eraAticq`` performs the transformations between two coordinate systems,
      with the details of the transformation being encoded into the ``astrom`` array.
  
      There are two issues with the version of aticq in ERFA. Both are associated
      with the handling of light deflection.
  
      The companion function ``eraAtciqz`` is meant to be its inverse. However, this
      is not true for directions close to the Solar centre, since the light deflection
      calculations are numerically unstable and therefore not reversible.
  
      This version sidesteps that problem by artificially reducing the light deflection
      for directions which are within 90 arcseconds of the Sun's position. This is the
      same approach used by the ERFA functions above, except that they use a threshold of
      9 arcseconds.
  
      In addition, ERFA's aticq assumes a distant source, so there is no difference between
      the object-Sun vector and the observer-Sun vector. This can lead to errors of up to a
      few arcseconds in the worst case (e.g a Venus transit).
  
      Parameters
      ----------
      srepr : `~astropy.coordinates.SphericalRepresentation`
          Astrometric GCRS or CIRS position of object from observer
      astrom : eraASTROM array
          ERFA astrometry context, as produced by, e.g. ``eraApci13`` or ``eraApcs13``
  
      Returns
      --------
      rc : float or `~numpy.ndarray`
          Right Ascension in radians
      dc : float or `~numpy.ndarray`
          Declination in radians
      """
      # ignore parallax effects if no distance, or far away
      srepr_distance = srepr.distance
      ignore_distance = srepr_distance.unit == u.one
  
      # RA, Dec to cartesian unit vectors
      pos = erfa.s2c(srepr.lon.radian, srepr.lat.radian)
  
      # Bias-precession-nutation, giving GCRS proper direction.
      ppr = erfa.trxp(astrom['bpn'], pos)
  
      # Aberration, giving GCRS natural direction
      d = np.zeros_like(ppr)
      for j in range(2):
          before = norm(ppr-d)
          after = erfa.ab(before, astrom['v'], astrom['em'], astrom['bm1'])
          d = after - before
      pnat = norm(ppr-d)
  
      # Light deflection by the Sun, giving BCRS coordinate direction
      d = np.zeros_like(pnat)
      for j in range(5):
          before = norm(pnat-d)
          if ignore_distance:
              # No distance to object, assume a long way away
              q = before
          else:
              # Find BCRS direction of Sun to object.
              # astrom['eh'] and astrom['em'] contain Sun to observer unit vector,
              # and distance, respectively.
              eh = astrom['em'][..., np.newaxis] * astrom['eh']
              # unit vector from Sun to object
              q = eh + srepr_distance[..., np.newaxis].to_value(u.au) * before
              sundist, q = erfa.pn(q)
              sundist = sundist[..., np.newaxis]
              # calculation above is extremely unstable very close to the sun
              # in these situations, default back to ldsun-style behaviour,
              # since this is reversible and drops to zero within stellar limb
              q = np.where(sundist > 1.0e-10, q, before)
  
          after = erfa.ld(1.0, before, q, astrom['eh'], astrom['em'], 1e-6)
          d = after - before
      pco = norm(pnat-d)
  
      # ICRS astrometric RA, Dec
      rc, dc = erfa.c2s(pco)
      return erfa.anp(rc), dc
  
  
  def atciqz(srepr, astrom):
      """
      A slightly modified version of the ERFA function ``eraAtciqz``.
  
      ``eraAtciqz`` performs the transformations between two coordinate systems,
      with the details of the transformation being encoded into the ``astrom`` array.
  
      There are two issues with the version of atciqz in ERFA. Both are associated
      with the handling of light deflection.
  
      The companion function ``eraAticq`` is meant to be its inverse. However, this
      is not true for directions close to the Solar centre, since the light deflection
      calculations are numerically unstable and therefore not reversible.
  
      This version sidesteps that problem by artificially reducing the light deflection
      for directions which are within 90 arcseconds of the Sun's position. This is the
      same approach used by the ERFA functions above, except that they use a threshold of
      9 arcseconds.
  
      In addition, ERFA's atciqz assumes a distant source, so there is no difference between
      the object-Sun vector and the observer-Sun vector. This can lead to errors of up to a
      few arcseconds in the worst case (e.g a Venus transit).
  
      Parameters
      ----------
      srepr : `~astropy.coordinates.SphericalRepresentation`
          Astrometric ICRS position of object from observer
      astrom : eraASTROM array
          ERFA astrometry context, as produced by, e.g. ``eraApci13`` or ``eraApcs13``
  
      Returns
      --------
      ri : float or `~numpy.ndarray`
          Right Ascension in radians
      di : float or `~numpy.ndarray`
          Declination in radians
      """
      # ignore parallax effects if no distance, or far away
      srepr_distance = srepr.distance
      ignore_distance = srepr_distance.unit == u.one
  
      # BCRS coordinate direction (unit vector).
      pco = erfa.s2c(srepr.lon.radian, srepr.lat.radian)
  
      # Find BCRS direction of Sun to object
      if ignore_distance:
          # No distance to object, assume a long way away
          q = pco
      else:
          # Find BCRS direction of Sun to object.
          # astrom['eh'] and astrom['em'] contain Sun to observer unit vector,
          # and distance, respectively.
          eh = astrom['em'][..., np.newaxis] * astrom['eh']
          # unit vector from Sun to object
          q = eh + srepr_distance[..., np.newaxis].to_value(u.au) * pco
          sundist, q = erfa.pn(q)
          sundist = sundist[..., np.newaxis]
          # calculation above is extremely unstable very close to the sun
          # in these situations, default back to ldsun-style behaviour,
          # since this is reversible and drops to zero within stellar limb
          q = np.where(sundist > 1.0e-10, q, pco)
  
      # Light deflection by the Sun, giving BCRS natural direction.
      pnat = erfa.ld(1.0, pco, q, astrom['eh'], astrom['em'], 1e-6)
  
      # Aberration, giving GCRS proper direction.
      ppr = erfa.ab(pnat, astrom['v'], astrom['em'], astrom['bm1'])
  
      # Bias-precession-nutation, giving CIRS proper direction.
      # Has no effect if matrix is identity matrix, in which case gives GCRS ppr.
      pi = erfa.rxp(astrom['bpn'], ppr)
  
      # CIRS (GCRS) RA, Dec
      ri, di = erfa.c2s(pi)
      return erfa.anp(ri), di
  
  
  def prepare_earth_position_vel(time):
      """
      Get barycentric position and velocity, and heliocentric position of Earth
  
      Parameters
      -----------
      time : `~astropy.time.Time`
          time at which to calculate position and velocity of Earth
  
      Returns
      --------
      earth_pv : `np.ndarray`
          Barycentric position and velocity of Earth, in au and au/day
      earth_helio : `np.ndarray`
          Heliocentric position of Earth in au
      """
      # this goes here to avoid circular import errors
      from astropy.coordinates.solar_system import (
          get_body_barycentric,
          get_body_barycentric_posvel,
          solar_system_ephemeris,
      )
      # get barycentric position and velocity of earth
  
      ephemeris = solar_system_ephemeris.get()
  
      # if we are using the builtin erfa based ephemeris,
      # we can use the fact that epv00 already provides all we need.
      # This avoids calling epv00 twice, once
      # in get_body_barycentric_posvel('earth') and once in
      # get_body_barycentric('sun')
      if ephemeris == 'builtin':
          jd1, jd2 = get_jd12(time, 'tdb')
          earth_pv_heliocentric, earth_pv = erfa.epv00(jd1, jd2)
          earth_heliocentric = earth_pv_heliocentric['p']
  
      # all other ephemeris providers probably don't have a shortcut like this
      else:
          earth_p, earth_v = get_body_barycentric_posvel('earth', time)
  
          # get heliocentric position of earth, preparing it for passing to erfa.
          sun = get_body_barycentric('sun', time)
          earth_heliocentric = (earth_p - sun).get_xyz(xyz_axis=-1).to_value(u.au)
  
          # Also prepare earth_pv for passing to erfa, which wants it as
          # a structured dtype.
          earth_pv = pav2pv(
              earth_p.get_xyz(xyz_axis=-1).to_value(u.au),
              earth_v.get_xyz(xyz_axis=-1).to_value(u.au / u.d)
          )
  
      return earth_pv, earth_heliocentric
  
  
  def get_offset_sun_from_barycenter(time, include_velocity=False, reverse=False):
      """
      Returns the offset of the Sun center from the solar-system barycenter (SSB).
  
      Parameters
      ----------
      time : `~astropy.time.Time`
          Time at which to calculate the offset
      include_velocity : `bool`
          If ``True``, attach the velocity as a differential.  Defaults to ``False``.
      reverse : `bool`
          If ``True``, return the offset of the barycenter from the Sun.  Defaults to ``False``.
  
      Returns
      -------
      `~astropy.coordinates.CartesianRepresentation`
          The offset
      """
      if include_velocity:
          # Import here to avoid a circular import
          from astropy.coordinates.solar_system import get_body_barycentric_posvel
          offset_pos, offset_vel = get_body_barycentric_posvel('sun', time)
          if reverse:
              offset_pos, offset_vel = -offset_pos, -offset_vel
          offset_vel = offset_vel.represent_as(CartesianDifferential)
          offset_pos = offset_pos.with_differentials(offset_vel)
  
      else:
          # Import here to avoid a circular import
          from astropy.coordinates.solar_system import get_body_barycentric
          offset_pos = get_body_barycentric('sun', time)
          if reverse:
              offset_pos = -offset_pos
  
      return offset_pos
  