Python astropy.coordinates.Longitude() Examples

def test_healpix_to_lonlat(self):
        lon, lat = self.pix.healpix_to_lonlat([1, 2, 3])

        assert isinstance(lon, Longitude)
        assert isinstance(lat, Latitude)

        index = self.pix.lonlat_to_healpix(lon, lat)

        assert_equal(index, [1, 2, 3])

        lon, lat = self.pix.healpix_to_lonlat([1, 2, 3],
                                              dx=[0.1, 0.2, 0.3],
                                              dy=[0.5, 0.4, 0.7])

        assert isinstance(lon, Longitude)
        assert isinstance(lat, Latitude)

        index, dx, dy = self.pix.lonlat_to_healpix(lon, lat, return_offsets=True)

        assert_equal(index, [1, 2, 3])
        assert_allclose(dx, [0.1, 0.2, 0.3])
        assert_allclose(dy, [0.5, 0.4, 0.7]) 

def test_healpix_to_lonlat(order):
    lon, lat = healpix_to_lonlat([1, 2, 3], 4, order=order)

    assert isinstance(lon, Longitude)
    assert isinstance(lat, Latitude)

    index = lonlat_to_healpix(lon, lat, 4, order=order)

    assert_equal(index, [1, 2, 3])

    lon, lat = healpix_to_lonlat([1, 2, 3], 4,
                                 dx=[0.1, 0.2, 0.3],
                                 dy=[0.5, 0.4, 0.7], order=order)

    assert isinstance(lon, Longitude)
    assert isinstance(lat, Latitude)

    index, dx, dy = lonlat_to_healpix(lon, lat, 4, order=order, return_offsets=True)

    assert_equal(index, [1, 2, 3])
    assert_allclose(dx, [0.1, 0.2, 0.3])
    assert_allclose(dy, [0.5, 0.4, 0.7]) 

def get_observer_pos_wrt_earth(sat_observatories_data, obs_radec, site_codes):
    """Compute position of observer at Earth's surface, with respect
    to the Earth, in equatorial frame, during 3 distinct instants.

       Args:
           sat_observatories_data (string): path to file containing COSPAR satellite tracking stations data.
           obs_radec (1x3 SkyCoord array): three rad/dec observations
           site_codes (1x3 int array): COSPAR codes of observation sites

       Returns:
           R (1x3 array): cartesian position vectors (observer wrt Earth)
    """
    R = np.array((np.zeros((3,)),np.zeros((3,)),np.zeros((3,))))
    # load MPC observatory data
    obsite1 = get_station_data(site_codes[0], sat_observatories_data)
    obsite2 = get_station_data(site_codes[1], sat_observatories_data)
    obsite3 = get_station_data(site_codes[2], sat_observatories_data)

    R[0] = observerpos_sat(obsite1['Latitude'], obsite1['Longitude'], obsite1['Elev'], obs_radec[0].obstime)
    R[1] = observerpos_sat(obsite2['Latitude'], obsite2['Longitude'], obsite2['Elev'], obs_radec[1].obstime)
    R[2] = observerpos_sat(obsite3['Latitude'], obsite3['Longitude'], obsite3['Elev'], obs_radec[2].obstime)

    return R 

def _erfa_sidereal_time(self, model):
        """Calculate a sidereal time using a IAU precession/nutation model."""

        from astropy.coordinates import Longitude

        erfa_function = model['function']
        erfa_parameters = [getattr(getattr(self, scale)._time, jd_part)
                           for scale in model['scales']
                           for jd_part in ('jd1', 'jd2_filled')]

        sidereal_time = erfa_function(*erfa_parameters)

        if self.masked:
            sidereal_time[self.mask] = np.nan

        return Longitude(sidereal_time, u.radian).to(u.hourangle) 

def test_EarthLocation_basic():
    greenwichel = EarthLocation.of_site('greenwich')
    lon, lat, el = greenwichel.to_geodetic()
    assert_quantity_allclose(lon, Longitude('0:0:0', unit=u.deg),
                             atol=10*u.arcsec)
    assert_quantity_allclose(lat, Latitude('51:28:40', unit=u.deg),
                             atol=1*u.arcsec)
    assert_quantity_allclose(el, 46*u.m, atol=1*u.m)

    names = EarthLocation.get_site_names()
    assert 'greenwich' in names
    assert 'example_site' in names

    with pytest.raises(KeyError) as exc:
        EarthLocation.of_site('nonexistent site')
    assert exc.value.args[0] == "Site 'nonexistent site' not in database. Use EarthLocation.get_site_names to see available sites." 

def test_online_sites():
    reg = get_downloaded_sites()

    keck = reg['keck']
    lon, lat, el = keck.to_geodetic()
    assert_quantity_allclose(lon, -Longitude('155:28.7', unit=u.deg),
                             atol=0.001*u.deg)
    assert_quantity_allclose(lat, Latitude('19:49.7', unit=u.deg),
                             atol=0.001*u.deg)
    assert_quantity_allclose(el, 4160*u.m, atol=1*u.m)

    names = reg.names
    assert 'keck' in names
    assert 'ctio' in names

    with pytest.raises(KeyError) as exc:
        reg['nonexistent site']
    assert exc.value.args[0] == "Site 'nonexistent site' not in database. Use the 'names' attribute to see available sites."

    with pytest.raises(KeyError) as exc:
        reg['kec']
    assert exc.value.args[0] == "Site 'kec' not in database. Use the 'names' attribute to see available sites. Did you mean one of: 'keck'?'" 

def test_builtin_sites():
    reg = get_builtin_sites()

    greenwich = reg['greenwich']
    lon, lat, el = greenwich.to_geodetic()
    assert_quantity_allclose(lon, Longitude('0:0:0', unit=u.deg),
                             atol=10*u.arcsec)
    assert_quantity_allclose(lat, Latitude('51:28:40', unit=u.deg),
                             atol=1*u.arcsec)
    assert_quantity_allclose(el, 46*u.m, atol=1*u.m)

    names = reg.names
    assert 'greenwich' in names
    assert 'example_site' in names

    with pytest.raises(KeyError) as exc:
        reg['nonexistent site']
    assert exc.value.args[0] == "Site 'nonexistent site' not in database. Use the 'names' attribute to see available sites." 

def test_distance_change():

    ra = Longitude("4:08:15.162342", unit=u.hour)
    dec = Latitude("-41:08:15.162342", unit=u.degree)
    c1 = ICRS(ra, dec, Distance(1, unit=u.kpc))

    oldx = c1.cartesian.x.value
    assert (oldx - 0.35284083171901953) < 1e-10

    # first make sure distances are immutible
    with pytest.raises(AttributeError):
        c1.distance = Distance(2, unit=u.kpc)

    # now x should increase with a bigger distance increases
    c2 = ICRS(ra, dec, Distance(2, unit=u.kpc))
    assert c2.cartesian.x.value == oldx * 2 

def from_tree(cls, node, ctx):
        wrap_angle = node['wrap_angle']
        return Longitude(super().from_tree(node, ctx), wrap_angle=wrap_angle) 

def _gal2idx(self, gal):
        """
        Converts from Galactic coordinates to pixel indices.

        Args:
            gal (:obj:`astropy.coordinates.SkyCoord`): Galactic coordinates. Must
                store an array of coordinates (i.e., not be scalar).

        Returns:
            ``j, k, mask`` - Pixel indices of the coordinates, as well as a mask
            of in-bounds coordinates. Outputs have the same shape as the input
            coordinates.
        """

        # Make sure that l is in domain [-180 deg, 180 deg)
        l = coordinates.Longitude(gal.l, wrap_angle=180.*units.deg)

        j = (self._inv_pix_scale * (l.deg - self._l_bounds[0])).astype('i4')
        k = (self._inv_pix_scale * (gal.b.deg - self._b_bounds[0])).astype('i4')

        idx = (j < 0) | (j >= self._shape[0]) | (k < 0) | (k >= self._shape[1])

        if np.any(idx):
            j[idx] = -1
            k[idx] = -1

        return j, k, ~idx 

def test_skycoord_ra_dec(tmpdir):

    ra = Longitude([1, 2, 3], unit=u.deg)  # Could also use Angle
    dec = np.array([4.5, 5.2, 6.3]) * u.deg  # Astropy Quantity
    c = SkyCoord(ra, dec, frame='icrs')
    tree = dict(coord=c)
    assert_roundtrip_tree(tree, tmpdir)

    c = SkyCoord(frame=ICRS, ra=ra, dec=dec, obstime='2001-01-02T12:34:56')
    tree = dict(coord=c)
    assert_roundtrip_tree(tree, tmpdir) 

def test_longitude(tmpdir):
    tree = {'angle': Longitude(-100, u.deg, wrap_angle=180*u.deg)}
    assert_roundtrip_tree(tree, tmpdir) 

def test_icrs_basic(tmpdir):
    wrap_angle = Angle(1.5, unit=units.rad)
    ra = Longitude(25, unit=units.deg, wrap_angle=wrap_angle)
    dec = Latitude(45, unit=units.deg)

    tree = {'coord': ICRS(ra=ra, dec=dec)}

    assert_roundtrip_tree(tree, tmpdir) 

def test_hcrs_basic(tmpdir):
    ra = Longitude(25, unit=units.deg)
    dec = Latitude(45, unit=units.deg)

    tree = {'coord': ICRS(ra=ra, dec=dec)}

    assert_roundtrip_tree(tree, tmpdir) 

def test_basic():
    lon1 = Longitude(1.23, "radian", wrap_angle='180d')
    s = pickle.dumps(lon1)
    lon2 = pickle.loads(s) 

def test_distance_in_coordinates():
    """
    test that distances can be created from quantities and that cartesian
    representations come out right
    """

    ra = Longitude("4:08:15.162342", unit=u.hour)
    dec = Latitude("-41:08:15.162342", unit=u.degree)
    coo = ICRS(ra, dec, distance=2*u.kpc)

    cart = coo.cartesian

    assert isinstance(cart.xyz, u.Quantity) 

def setup(self):
        lon = Longitude(np.arange(0, 24, 4), u.hourangle)
        lat = Latitude(np.arange(-90, 91, 30), u.deg)
        # With same-sized arrays, no attributes
        self.s0 = ICRS(lon[:, np.newaxis] * np.ones(lat.shape),
                       lat * np.ones(lon.shape)[:, np.newaxis])
        # Make an AltAz frame since that has many types of attributes.
        # Match one axis with times.
        self.obstime = (Time('2012-01-01') +
                        np.arange(len(lon))[:, np.newaxis] * u.s)
        # And another with location.
        self.location = EarthLocation(20.*u.deg, lat, 100*u.m)
        # Ensure we have a quantity scalar.
        self.pressure = 1000 * u.hPa
        # As well as an array.
        self.temperature = np.random.uniform(
            0., 20., size=(lon.size, lat.size)) * u.deg_C
        self.s1 = AltAz(az=lon[:, np.newaxis], alt=lat,
                        obstime=self.obstime,
                        location=self.location,
                        pressure=self.pressure,
                        temperature=self.temperature)
        # For some tests, also try a GCRS, since that has representation
        # attributes.  We match the second dimension (via the location)
        self.obsgeoloc, self.obsgeovel = self.location.get_gcrs_posvel(
            self.obstime[0, 0])
        self.s2 = GCRS(ra=lon[:, np.newaxis], dec=lat,
                       obstime=self.obstime,
                       obsgeoloc=self.obsgeoloc,
                       obsgeovel=self.obsgeovel)
        # For completeness, also some tests on an empty frame.
        self.s3 = GCRS(obstime=self.obstime,
                       obsgeoloc=self.obsgeoloc,
                       obsgeovel=self.obsgeovel)
        # And make a SkyCoord
        self.sc = SkyCoord(ra=lon[:, np.newaxis], dec=lat, frame=self.s3) 

def datetime2sidereal(time: datetime, lon_radians: float, *, use_astropy: bool = True) -> float:
    """
    Convert ``datetime`` to local sidereal time

    from D. Vallado "Fundamentals of Astrodynamics and Applications"


    time : datetime.datetime
        time to convert
    lon_radians : float
        longitude (radians)
    use_astropy : bool, optional
        use AstroPy for conversion (False is Vallado)

    Results
    -------

    tsr : float
        Local sidereal time
    """
    if isinstance(time, (tuple, list)):
        return [datetime2sidereal(t, lon_radians) for t in time]

    if use_astropy and Time is not None:
        tsr = Time(time).sidereal_time(kind="apparent", longitude=Longitude(lon_radians, unit=u.radian)).radian
    else:
        jd = juliandate(str2dt(time))
        # %% Greenwich Sidereal time RADIANS
        gst = greenwichsrt(jd)
        # %% Algorithm 15 p. 188 rotate GST to LOCAL SIDEREAL TIME
        tsr = gst + lon_radians

    return tsr 

def observerpos_mpc(long, parallax_s, parallax_c, t_utc):
    """Compute geocentric observer position at UTC instant t_utc, for Sun-centered orbits,
    at a given observation site defined by its longitude, and parallax constants S and C.
    Formula taken from top of page 266, chapter 5, Orbital Mechanics book (Curtis).
    The parallax constants S and C are defined by:
    S=rho cos phi' C=rho sin phi', where
    rho: slant range
    phi': geocentric latitude

       Args:
           long (float): longitude of observing site
           parallax_s (float): parallax constant S of observing site
           parallax_c (float): parallax constant C of observing site
           t_utc (astropy.time.Time): UTC time of observation

       Returns:
           1x3 numpy array: cartesian components of observer's geocentric position
    """
    # Earth's equatorial radius in kilometers
    # Re = cts.R_earth.to(uts.Unit('km')).value

    # define Longitude object for the observation site longitude
    long_site = Longitude(long, uts.degree, wrap_angle=360.0*uts.degree)
    # compute sidereal time of observation at site
    t_site_lmst = t_utc.sidereal_time('mean', longitude=long_site)
    lmst_rad = t_site_lmst.rad # np.deg2rad(lmst_hrs*15.0) # radians

    # compute cartesian components of geocentric (non rotating) observer position
    x_gc = Re*parallax_c*np.cos(lmst_rad)
    y_gc = Re*parallax_c*np.sin(lmst_rad)
    z_gc = Re*parallax_s

    return np.array((x_gc,y_gc,z_gc)) 

def _coords2vec(self, coords):
        """
        Converts from sky coordinates to unit vectors. Before conversion to unit
        vectors, the coordiantes are transformed to the coordinate system used
        internally by the :obj:`UnstructuredDustMap`, which can be set during
        initialization of the class.

        Args:
            coords (:obj:`astropy.coordinates.SkyCoord`): Input coordinates to
                convert to unit vectors.

        Returns:
            Cartesian unit vectors corresponding to the input coordinates, after
            transforming to the coordinate system used internally by the
            :obj:`UnstructuredDustMap`.
        """

        # c = coords.transform_to(self._frame)
        # vec = np.empty((c.shape[0], 2), dtype='f8')
        # vec[:,0] = coordinates.Longitude(coords.l, wrap_angle=360.*units.deg).deg[:]
        # vec[:,1] = coords.b.deg[:]
        # return np.radians(vec)

        c = coords.transform_to(self._frame).represent_as('cartesian')
        vec_norm = np.sqrt(c.x**2 + c.y**2 + c.z**2)

        vec = np.empty((c.shape[0], 3), dtype=c.x.dtype)
        vec[:,0] = (c.x / vec_norm).value[:]
        vec[:,1] = (c.y / vec_norm).value[:]
        vec[:,2] = (c.z / vec_norm).value[:]

        return vec 

def test_unit_representation_subclass():

    class Longitude180(Longitude):
        def __new__(cls, angle, unit=None, wrap_angle=180*u.deg, **kwargs):
            self = super().__new__(cls, angle, unit=unit, wrap_angle=wrap_angle,
                                   **kwargs)
            return self

    class UnitSphericalWrap180Representation(UnitSphericalRepresentation):
        attr_classes = OrderedDict([('lon', Longitude180),
                                    ('lat', Latitude)])

    class SphericalWrap180Representation(SphericalRepresentation):
        attr_classes = OrderedDict([('lon', Longitude180),
                                    ('lat', Latitude),
                                    ('distance', u.Quantity)])

        _unit_representation = UnitSphericalWrap180Representation

    class MyFrame(ICRS):
        default_representation = SphericalWrap180Representation
        frame_specific_representation_info = {
            'spherical': [
                RepresentationMapping('lon', 'ra'),
                RepresentationMapping('lat', 'dec')]
        }
        frame_specific_representation_info['unitsphericalwrap180'] = \
            frame_specific_representation_info['sphericalwrap180'] = \
            frame_specific_representation_info['spherical']

    @frame_transform_graph.transform(FunctionTransform,
                                     MyFrame, astropy.coordinates.ICRS)
    def myframe_to_icrs(myframe_coo, icrs):
        return icrs.realize_frame(myframe_coo._data)

    f = MyFrame(10*u.deg, 10*u.deg)
    assert isinstance(f._data, UnitSphericalWrap180Representation)
    assert isinstance(f.ra, Longitude180)

    g = f.transform_to(astropy.coordinates.ICRS)
    assert isinstance(g, astropy.coordinates.ICRS)
    assert isinstance(g._data, UnitSphericalWrap180Representation)

    frame_transform_graph.remove_transform(MyFrame,
                                           astropy.coordinates.ICRS,
                                           None) 

def sidereal_time(self, kind, longitude=None, model=None):
        """Calculate sidereal time.

        Parameters
        ---------------
        kind : str
            ``'mean'`` or ``'apparent'``, i.e., accounting for precession
            only, or also for nutation.
        longitude : `~astropy.units.Quantity`, `str`, or `None`; optional
            The longitude on the Earth at which to compute the sidereal time.
            Can be given as a `~astropy.units.Quantity` with angular units
            (or an `~astropy.coordinates.Angle` or
            `~astropy.coordinates.Longitude`), or as a name of an
            observatory (currently, only ``'greenwich'`` is supported,
            equivalent to 0 deg).  If `None` (default), the ``lon`` attribute of
            the Time object is used.
        model : str or `None`; optional
            Precession (and nutation) model to use.  The available ones are:
            - {0}: {1}
            - {2}: {3}
            If `None` (default), the last (most recent) one from the appropriate
            list above is used.

        Returns
        -------
        sidereal time : `~astropy.coordinates.Longitude`
            Sidereal time as a quantity with units of hourangle
        """  # docstring is formatted below

        from astropy.coordinates import Longitude

        if kind.lower() not in SIDEREAL_TIME_MODELS.keys():
            raise ValueError('The kind of sidereal time has to be {}'.format(
                ' or '.join(sorted(SIDEREAL_TIME_MODELS.keys()))))

        available_models = SIDEREAL_TIME_MODELS[kind.lower()]

        if model is None:
            model = sorted(available_models.keys())[-1]
        else:
            if model.upper() not in available_models:
                raise ValueError(
                    'Model {} not implemented for {} sidereal time; '
                    'available models are {}'
                    .format(model, kind, sorted(available_models.keys())))

        if longitude is None:
            if self.location is None:
                raise ValueError('No longitude is given but the location for '
                                 'the Time object is not set.')
            longitude = self.location.lon
        elif longitude == 'greenwich':
            longitude = Longitude(0., u.degree,
                                  wrap_angle=180.*u.degree)
        else:
            # sanity check on input
            longitude = Longitude(longitude, u.degree,
                                  wrap_angle=180.*u.degree)

        gst = self._erfa_sidereal_time(available_models[model.upper()])
        return Longitude(gst + longitude, u.hourangle) 

def t_radec_res_vec_sat(x, inds, iod_object_data, sat_observatories_data, rov):
    """Compute vector of observed minus computed (O-C) residuals for ra/dec Earth-orbiting satellite observations
    with pre-computed observed radec values vector. Assumes ra/dec observed values vector
    is contained in rov, and they are stored as rov = [ra1, dec1, ra2, dec2, ...].

       Args:
           x (1x6 float array): set of orbital elements (a, e, taup, omega, I, Omega)
           inds (int array): line numbers of data in file
           iod_object_data (ndarray): observation data
           sat_observatories_data (ndarray): satellite tracking stations data
           rov (1xlen(inds) float-like array): vector of observed ra/dec values

       Returns:
           rv (1xlen(inds) array): vector of ra/dec (O-C) residuals.
           tv (1xlen(inds) array): vector of observation times.
    """
    rv = np.zeros((2*len(inds)))
    tv = np.zeros((len(inds)))
    for i in range(0,len(inds)):
        indm1 = inds[i]-1
        # extract observations data
        td = timedelta(hours=1.0*iod_object_data['hr'][indm1],
                       minutes=1.0*iod_object_data['min'][indm1],
                       seconds=(iod_object_data['sec'][indm1]+iod_object_data['msec'][indm1]/1000.0))
        timeobs = Time( datetime(iod_object_data['yr'][indm1],
                                 iod_object_data['month'][indm1],
                                 iod_object_data['day'][indm1]) + td )
        t_jd = timeobs.jd
        site_code = iod_object_data['station'][indm1]
        obsite = get_station_data(site_code, sat_observatories_data)
        # object position wrt to Earth
        xyz_obj = orbel2xyz(t_jd, cts.GM_earth.to(uts.Unit('km3 / day2')).value, x[0], x[1], x[2], x[3], x[4], x[5])
        # observer position wrt to Earth
        xyz_oe = observerpos_sat(obsite['Latitude'], obsite['Longitude'], obsite['Elev'], timeobs)
        # object position wrt observer (unnormalized LOS vector)
        rho_vec = xyz_obj - xyz_oe
        # compute normalized LOS vector
        rho_vec_norm = np.linalg.norm(rho_vec, ord=2)
        rho_vec_unit = rho_vec/rho_vec_norm
        # compute RA, Dec
        cosd_cosa = rho_vec_unit[0]
        cosd_sina = rho_vec_unit[1]
        sind = rho_vec_unit[2]
        # make sure computed RA (ra_comp) is always within [0.0, 2.0*np.pi]
        ra_comp = np.mod(np.arctan2(cosd_sina, cosd_cosa), 2.0*np.pi)
        dec_comp = np.arcsin(sind)
        #compute angle difference, taking always the smallest difference
        diff_ra = angle_diff_rad(rov[2*i-2], ra_comp)
        diff_dec = angle_diff_rad(rov[2*i-1], dec_comp)
        # store O-C residual into vector (O-C = "Observed minus Computed")
        rv[2*i-2], rv[2*i-1] = diff_ra, diff_dec
        tv[i] = t_jd
    return tv, rv

# compute auxiliary vector of observed ra,dec values 

def radec_res_vec_rov_sat(x, inds, iod_object_data, sat_observatories_data, rov):
    """Compute vector of observed minus computed (O-C) residuals for ra/dec Earth-orbiting satellite observations
    with pre-computed observed radec values vector. Assumes ra/dec observed values vector
    is contained in rov, and they are stored as rov = [ra1, dec1, ra2, dec2, ...].

       Args:
           x (1x6 float array): set of orbital elements (a, e, taup, omega, I, Omega)
           inds (int array): line numbers of data in file
           iod_object_data (ndarray): observation data
           sat_observatories_data (ndarray): satellite tracking stations data
           rov (1xlen(inds) float-like array): vector of observed ra/dec values

       Returns:
           rv (1xlen(inds) array): vector of ra/dec (O-C) residuals.
    """
    rv = np.zeros((2*len(inds)))
    for i in range(0,len(inds)):
        indm1 = inds[i]-1
        # extract observations data
        td = timedelta(hours=1.0*iod_object_data['hr'][indm1],
                       minutes=1.0*iod_object_data['min'][indm1],
                       seconds=(iod_object_data['sec'][indm1]+iod_object_data['msec'][indm1]/1000.0))
        timeobs = Time( datetime(iod_object_data['yr'][indm1],
                                 iod_object_data['month'][indm1],
                                 iod_object_data['day'][indm1]) + td )
        site_code = iod_object_data['station'][indm1]
        obsite = get_station_data(site_code, sat_observatories_data)
        # object position wrt to Earth
        xyz_obj = orbel2xyz(timeobs.jd,
                            cts.GM_earth.to(uts.Unit('km3 / day2')).value,
                            x[0], x[1], x[2], x[3], x[4], x[5])
        # observer position wrt to Earth
        xyz_oe = observerpos_sat(obsite['Latitude'], obsite['Longitude'], obsite['Elev'], timeobs)
        # object position wrt observer (unnormalized LOS vector)
        rho_vec = xyz_obj - xyz_oe
        # compute normalized LOS vector
        rho_vec_norm = np.linalg.norm(rho_vec, ord=2)
        rho_vec_unit = rho_vec/rho_vec_norm
        # compute RA, Dec
        cosd_cosa = rho_vec_unit[0]
        cosd_sina = rho_vec_unit[1]
        sind = rho_vec_unit[2]
        # make sure computed RA (ra_comp) is always within [0.0, 2.0*np.pi]
        ra_comp = np.mod(np.arctan2(cosd_sina, cosd_cosa), 2.0*np.pi)
        dec_comp = np.arcsin(sind)
        #compute angle difference, taking always the smallest difference
        diff_ra = angle_diff_rad(rov[2*i-2], ra_comp)
        diff_dec = angle_diff_rad(rov[2*i-1], dec_comp)
        # store O-C residual into vector (O-C = "Observed minus Computed")
        rv[2*i-2], rv[2*i-1] = diff_ra, diff_dec
    return rv

# compute residuals vector for ra/dec observations; return observation times and residual vector
# inds = obs_arr 

def observerpos_sat(lat, long, elev, t_utc):
    """Compute geocentric observer position at UTC instant t_utc, for Earth-centered orbits,
    at a given observation site defined by its longitude, geodetic latitude and elevation above reference ellipsoid.
    Formula taken from bottom of page 265 (Eq. 5.56), chapter 5, Orbital Mechanics book (Curtis).

       Args:
           lat (float): geodetic latitude (deg)
           long (float): longitude (deg)
           elev (float): elevation above reference ellipsoid (m)
           t_utc (astropy.time.Time): UTC time of observation

       Returns:
           1x3 numpy array: cartesian components of observer's geocentric position
    """

    # Earth's equatorial radius in kilometers
    # Re = cts.R_earth.to(uts.Unit('km')).value

    # define Longitude object for the observation site longitude
    long_site = Longitude(long, uts.degree, wrap_angle=180.0*uts.degree)
    # compute sidereal time of observation at site
    t_site_lmst = t_utc.sidereal_time('mean', longitude=long_site)
    lmst_rad = t_site_lmst.rad # np.deg2rad(lmst_hrs*15.0) # radians

    # latitude
    phi_rad = np.deg2rad(lat)

    # convert ellipsoid coordinates to cartesian
    cos_phi = np.cos( phi_rad )
    cos_phi_cos_theta = cos_phi*np.cos( lmst_rad )
    cos_phi_sin_theta = cos_phi*np.sin( lmst_rad )
    sin_phi = np.sin( phi_rad )
    denum = np.sqrt(1.0-(2.0*earth_f-earth_f**2)*sin_phi**2)
    r_xy = Re/denum+elev/1000.0
    r_z = Re*((1.0-earth_f)**2)/denum+elev/1000.0

    # compute cartesian components of geocentric (non-rotating) observer position
    x_gc = r_xy*cos_phi_cos_theta
    y_gc = r_xy*cos_phi_sin_theta
    z_gc = r_z*sin_phi

    return np.array((x_gc,y_gc,z_gc)) 

def healpix_to_lonlat(healpix_index, nside, dx=None, dy=None, order='ring'):
    """
    Convert HEALPix indices (optionally with offsets) to longitudes/latitudes.

    If no offsets (``dx`` and ``dy``) are provided, the coordinates will default
    to those at the center of the HEALPix pixels.

    Parameters
    ----------
    healpix_index : int or `~numpy.ndarray`
        HEALPix indices (as a scalar or array)
    nside : int or `~numpy.ndarray`
        Number of pixels along the side of each of the 12 top-level HEALPix tiles
    dx, dy : float or `~numpy.ndarray`, optional
        Offsets inside the HEALPix pixel, which must be in the range [0:1],
        where 0.5 is the center of the HEALPix pixels (as scalars or arrays)
    order : { 'nested' | 'ring' }, optional
        Order of HEALPix pixels

    Returns
    -------
    lon : :class:`~astropy.coordinates.Longitude`
        The longitude values
    lat : :class:`~astropy.coordinates.Latitude`
        The latitude values
    """

    _validate_nside(nside)

    if _validate_order(order) == 'ring':
        func = _core.healpix_ring_to_lonlat
    else:  # _validate_order(order) == 'nested'
        func = _core.healpix_nested_to_lonlat

    if dx is None:
        dx = 0.5
    else:
        _validate_offset('x', dx)
    if dy is None:
        dy = 0.5
    else:
        _validate_offset('y', dy)

    nside = np.asarray(nside, dtype=np.intc)

    lon, lat = func(healpix_index, nside, dx, dy)

    lon = Longitude(lon, unit=u.rad, copy=False)
    lat = Latitude(lat, unit=u.rad, copy=False)

    return lon, lat