def great_circle(**kwargs):
    """
        Named arguments:
        distance  = distance to travel, or numpy array of distances
        azimuth   = angle, in DEGREES of HEADING from NORTH, or numpy array of azimuths
        latitude  = latitude, in DECIMAL DEGREES, or numpy array of latitudes
        longitude = longitude, in DECIMAL DEGREES, or numpy array of longitudes
        rmajor    = radius of earth's major axis. default=6378137.0 (WGS84)
        rminor    = radius of earth's minor axis. default=6356752.3142 (WGS84)

        Returns a dictionary with:
        'latitude' in decimal degrees
        'longitude' in decimal degrees
        'reverse_azimuth' in decimal degrees

    """

    distance  = kwargs.pop('distance')
    azimuth   = np.radians(kwargs.pop('azimuth'))
    latitude  = np.radians(kwargs.pop('latitude'))
    longitude = np.radians(kwargs.pop('longitude'))
    rmajor    = kwargs.pop('rmajor', 6378137.0)
    rminor    = kwargs.pop('rminor', 6356752.3142)
    f         = (rmajor - rminor) / rmajor

    vector_pt = np.vectorize(vinc_pt)
    lat_result, lon_result, angle_result = vector_pt(f, rmajor,
                                                     latitude,
                                                     longitude,
                                                     azimuth,
                                                     distance)
    return {'latitude': np.degrees(lat_result),
            'longitude': np.degrees(lon_result),
            'reverse_azimuth': np.degrees(angle_result)}

    def from_parameters(a, b, c, alpha, beta, gamma):
        """
        Create a Lattice using unit cell lengths and angles (in degrees).

        Args:
            a (float): *a* lattice parameter.
            b (float): *b* lattice parameter.
            c (float): *c* lattice parameter.
            alpha (float): *alpha* angle in degrees.
            beta (float): *beta* angle in degrees.
            gamma (float): *gamma* angle in degrees.

        Returns:
            Lattice with the specified lattice parameters.
        """

        alpha_r = radians(alpha)
        beta_r = radians(beta)
        gamma_r = radians(gamma)
        val = (np.cos(alpha_r) * np.cos(beta_r) - np.cos(gamma_r))\
            / (np.sin(alpha_r) * np.sin(beta_r))
        #Sometimes rounding errors result in values slightly > 1.
        val = val if abs(val) <= 1 else val / abs(val)
        gamma_star = np.arccos(val)
        vector_a = [a * np.sin(beta_r), 0.0, a * np.cos(beta_r)]
        vector_b = [-b * np.sin(alpha_r) * np.cos(gamma_star),
                    b * np.sin(alpha_r) * np.sin(gamma_star),
                    b * np.cos(alpha_r)]
        vector_c = [0.0, 0.0, float(c)]
        return Lattice([vector_a, vector_b, vector_c])

    def testNativeLonLatVector(self):
        """
        Test that nativeLonLatFromRaDec works in a vectorized way; we do this
        by performing a bunch of tansformations passing in ra and dec as numpy arrays
        and then comparing them to results computed in an element-wise way
        """

        obs = ObservationMetaData(pointingRA=123.0, pointingDec=43.0, mjd=53467.2)

        raPoint = 145.0
        decPoint = -35.0

        nSamples = 100
        np.random.seed(42)
        raList = np.random.random_sample(nSamples)*360.0
        decList = np.random.random_sample(nSamples)*180.0 - 90.0

        lonList, latList = nativeLonLatFromRaDec(raList, decList, obs)

        for rr, dd, lon, lat in zip(raList, decList, lonList, latList):
            lonControl, latControl = nativeLonLatFromRaDec(rr, dd, obs)
            distance = arcsecFromRadians(haversine(np.radians(lon), np.radians(lat),
                                                   np.radians(lonControl), np.radians(latControl)))

            self.assertLess(distance, 0.0001)

def b1950toj2000(ra, dec):
    """
    Convert B1950 to J2000 coordinates.

    This routine is based on the technique described at
    http://www.stargazing.net/kepler/b1950.html
    """

    # Convert to radians
    ra = np.radians(ra)
    dec = np.radians(dec)

    # Convert RA, Dec to rectangular coordinates
    x = np.cos(ra) * np.cos(dec)
    y = np.sin(ra) * np.cos(dec)
    z = np.sin(dec)

    # Apply the precession matrix
    x2 = P1[0, 0] * x + P1[1, 0] * y + P1[2, 0] * z
    y2 = P1[0, 1] * x + P1[1, 1] * y + P1[2, 1] * z
    z2 = P1[0, 2] * x + P1[1, 2] * y + P1[2, 2] * z

    # Convert the new rectangular coordinates back to RA, Dec
    ra = np.arctan2(y2, x2)
    dec = np.arcsin(z2)

    # Convert to degrees
    ra = np.degrees(ra)
    dec = np.degrees(dec)

    # Make sure ra is between 0. and 360.
    ra = np.mod(ra, 360.0)
    dec = np.mod(dec + 90.0, 180.0) - 90.0

    return ra, dec

    def __init__(self, end_effector_length):
        # create joint limit dicts
        self.joint_lim_dict = {}
        self.joint_lim_dict["right_arm"] = {
            "max": np.radians([120.00, 122.15, 77.5, 144.0, 122.0, 45.0, 45.0]),
            "min": np.radians([-47.61, -20.0, -77.5, 0.0, -80.0, -45.0, -45.0]),
        }

        self.joint_lim_dict["left_arm"] = {
            "max": np.radians([120.00, 20.0, 77.5, 144.0, 80.0, 45.0, 45.0]),
            "min": np.radians([-47.61, -122.15, -77.5, 0.0, -122.0, -45.0, -45.0]),
        }

        end_effector_length += 0.0135 + 0.04318  # add wrist linkange and FT sensor lengths
        self.setup_kdl_mekabot(end_effector_length)
        q_guess_pkl_l = (
            os.environ["HOME"]
            + "/svn/gt-ros-pkg/hrl/hrl_arm_control/hrl_cody_arms/src/hrl_cody_arms/q_guess_left_dict.pkl"
        )
        q_guess_pkl_r = (
            os.environ["HOME"]
            + "/svn/gt-ros-pkg/hrl/hrl_arm_control/hrl_cody_arms/src/hrl_cody_arms/q_guess_right_dict.pkl"
        )

        self.q_guess_dict_left = ut.load_pickle(q_guess_pkl_l)
        self.q_guess_dict_right = ut.load_pickle(q_guess_pkl_r)

 def draw_scater(self, sw_list):
     """
     Descript. : draws data collection item on the scatter
                 subwedge is represented as list:
                 collection_id, sw_id, first_image, num_images, 
                 osc_start, osc_full_range
     """
     self.axes.clear()
     col_count = 0
     for sw_index, sw in enumerate(sw_list):
         bars = self.axes.bar(
             np.radians(sw[4]), 1, width=np.radians(sw[5]), bottom=sw[0], color=Qt4_widget_colors.TASK_GROUP[sw[0]]
         )
         x_mid = bars[0].get_bbox().xmin + (bars[0].get_bbox().xmax - bars[0].get_bbox().xmin) / 2.0
         y_mid = bars[0].get_bbox().ymin + (bars[0].get_bbox().ymax - bars[0].get_bbox().ymin) / 2.0
         self.axes.text(
             x_mid,
             y_mid,
             "%d (%d:%d)" % (sw_index + 1, sw[0] + 1, sw[1] + 1),
             horizontalalignment="center",
             verticalalignment="center",
             weight="bold",
         )
         if sw[0] > col_count:
             col_count = sw[0]
     self.axes.set_yticks(np.arange(1, col_count + 2))
     self.fig.canvas.draw_idle()

    def xy2lonlat(self, x, y):
        """Calculate x,y in own projection from given lon,lat (scalars/arrays).
        """
        if self.projected is True:
            if self.proj.is_latlong():
                return x, y
            else:
                if 'ob_tran' in self.proj4:
                    logging.info('NB: Converting degrees to radians ' +
                                 'due to ob_tran srs')
                    x = np.radians(np.array(x))
                    y = np.radians(np.array(y))
                return self.proj(x, y, inverse=True)
        else:
            np.seterr(invalid='ignore')  # Disable warnings for nan-values
            y = np.atleast_1d(np.array(y))
            x = np.atleast_1d(np.array(x))

            # NB: mask coordinates outside domain
            x[x < self.xmin] = np.nan
            x[x > self.xmax] = np.nan
            y[y < self.ymin] = np.nan
            y[y < self.ymin] = np.nan

            lon = map_coordinates(self.lon, [y, x], order=1,
                                  cval=np.nan, mode='nearest')
            lat = map_coordinates(self.lat, [y, x], order=1,
                                  cval=np.nan, mode='nearest')
            return (lon, lat)

    def testGetRotTelPos(self):
        rotSkyList = np.random.random_sample(len(self.raList))*2.0*np.pi
        mjd=56789.3
        obsTemp = ObservationMetaData(mjd=mjd, site=Site(longitude=self.lon, latitude=self.lat, name='LSST'))

        rotTelRad = utils._getRotTelPos(self.raList, self.decList,
                                        obsTemp, rotSkyList)

        rotTelDeg = utils.getRotTelPos(np.degrees(self.raList),
                                       np.degrees(self.decList),
                                       obsTemp, np.degrees(rotSkyList))

        np.testing.assert_array_almost_equal(rotTelRad, np.radians(rotTelDeg), 10)

        rotTelRad = utils._getRotTelPos(self.raList, self.decList,
                                        obsTemp, rotSkyList[0])

        rotTelDeg = utils.getRotTelPos(np.degrees(self.raList),
                                       np.degrees(self.decList),
                                       obsTemp, np.degrees(rotSkyList[0]))

        np.testing.assert_array_almost_equal(rotTelRad, np.radians(rotTelDeg), 10)


        for ra, dec, rotSky in \
        zip(self.raList, self.decList, rotSkyList):

            obsTemp = ObservationMetaData(mjd=mjd, site=Site(longitude=self.lon, latitude=self.lat, name='LSST'))

            rotTelRad = utils._getRotTelPos(ra, dec, obsTemp, rotSky)

            rotTelDeg = utils.getRotTelPos(np.degrees(ra), np.degrees(dec),
                                           obsTemp, np.degrees(rotSky))

            self.assertAlmostEqual(rotTelRad, np.radians(rotTelDeg), 10)

    def testObservedFromICRS(self):
        obs = ObservationMetaData(pointingRA=35.0, pointingDec=-45.0,
                                  mjd=43572.0)
        for pmRaList in [self.pm_raList, None]:
            for pmDecList in [self.pm_decList, None]:
                for pxList in [self.pxList, None]:
                    for vRadList in [self.v_radList, None]:
                        for includeRefraction in [True, False]:


                            raRad, decRad = utils._observedFromICRS(self.raList, self.decList,
                                                                         pm_ra=pmRaList, pm_dec=pmDecList,
                                                                         parallax=pxList, v_rad=vRadList,
                                                                         obs_metadata=obs, epoch=2000.0,
                                                                         includeRefraction=includeRefraction)

                            raDeg, decDeg = utils.observedFromICRS(np.degrees(self.raList), np.degrees(self.decList),
                                                                         pm_ra=utils.arcsecFromRadians(pmRaList),
                                                                         pm_dec=utils.arcsecFromRadians(pmDecList),
                                                                         parallax=utils.arcsecFromRadians(pxList),
                                                                         v_rad=vRadList,
                                                                         obs_metadata=obs, epoch=2000.0,
                                                                     includeRefraction=includeRefraction)



                            dRa = utils.arcsecFromRadians(raRad-np.radians(raDeg))
                            np.testing.assert_array_almost_equal(dRa, np.zeros(self.nStars), 9)

                            dDec = utils.arcsecFromRadians(decRad-np.radians(decDeg))
                            np.testing.assert_array_almost_equal(dDec, np.zeros(self.nStars), 9)

def fk4_e_terms(equinox):
    """
    Return the e-terms of aberation vector

    Parameters
    ----------
    equinox : Time object
        The equinox for which to compute the e-terms
    """

    from . import earth_orientation as earth

    # Constant of aberration at J2000
    k = 0.0056932

    # Eccentricity of the Earth's orbit
    e = earth.eccentricity(equinox.jd)
    e = np.radians(e)

    # Mean longitude of perigee of the solar orbit
    g = earth.mean_lon_of_perigee(equinox.jd)
    g = np.radians(g)

    # Obliquity of the ecliptic
    o = earth.obliquity(equinox.jd, algorithm=1980)
    o = np.radians(o)

    return e * k * np.sin(g), \
           -e * k * np.cos(g) * np.cos(o), \
           -e * k * np.cos(g) * np.sin(o)

    def __init__(self, irf, alpha, min_fsig=0.0, redge="0.0/1.0"):

        self._irf = irf
        #        self._th68 = self._irf._psf.counts
        #        self._bkg_rate = (self._det.proton_wcounts_density +
        #                          self._det.electron_wcounts_density)/(50.0*Units.hr*
        #                                                               Units.deg2)

        self._loge_axis = Axis.create(np.log10(Units.gev) + 1.4, np.log10(Units.gev) + 3.6, 22)
        self._loge = self._loge_axis.center
        self._dloge = self._loge_axis.width

        rmin, rmax = [float(t) for t in redge.split("/")]

        self._psi_axis = Axis(np.linspace(np.radians(0.0), np.radians(1.0), 101))
        self._psi = self._psi_axis.center
        self._domega = np.pi * (np.power(self._psi_axis.edges[1:], 2) - np.power(self._psi_axis.edges[:-1], 2))

        self._aeff = self._irf.aeff(self._loge - Units.log10_mev)
        self._bkg_rate = irf.bkg(self._loge - Units.log10_mev) * self._dloge

        self._redge = [np.radians(rmin), np.radians(rmax)]
        self._msk = (self._psi > self._redge[0]) & (self._psi < self._redge[1])

        self._domega_sig = np.sum(self._domega[self._msk])
        self._domega_bkg = self._domega_sig / alpha

        self._min_fsig = min_fsig
        self._alpha = alpha
        self._alpha_bin = self._domega / self._domega_bkg

 def convert_arc(en):
     angles = [np.radians(en.data['50']), np.radians(en.data['51'])]
     C = [en.data['10'], en.data['20']] #, en.data['30']]
     R = en.data['40']
     points = angles_to_threepoint(angles, C[0:2], R).tolist()   
     entities.append(Arc(len(vertices) + np.arange(3), closed=False))
     vertices.extend(points)

def get_fringe(qubic,horns, SOURCE_POWER=1., SOURCE_THETA=np.radians(0), SOURCE_PHI=np.radians(45),
               NU=150e9,DELTANU_NU=0.25,npts=512,display=True,background=True):
    horn_i=horns[0]
    horn_j=horns[1]
    all_open=~qubic.horn.removed
    all_open_i=all_open.copy()
    all_open_i[horn_i[0],horn_i[1]] = False
    all_open_j=all_open.copy()
    all_open_j[horn_j[0],horn_j[1]] = False
    all_open_ij=all_open.copy()
    all_open_ij[horn_i[0],horn_i[1]] = False
    all_open_ij[horn_j[0],horn_j[1]] = False

    Sall,Nall=get_fpimage(qubic,SOURCE_POWER=SOURCE_POWER,SOURCE_THETA=SOURCE_THETA, SOURCE_PHI=SOURCE_PHI,
                          NU=NU,DELTANU_NU=DELTANU_NU,npts=npts,display=display,HORN_OPEN=all_open,
                          background=background)
    S_i,N_i=get_fpimage(qubic,SOURCE_POWER=SOURCE_POWER,SOURCE_THETA=SOURCE_THETA, SOURCE_PHI=SOURCE_PHI,
                          NU=NU,DELTANU_NU=DELTANU_NU,npts=npts,display=display,HORN_OPEN=all_open_i,
                          background=background)
    S_j,N_j=get_fpimage(qubic,SOURCE_POWER=SOURCE_POWER,SOURCE_THETA=SOURCE_THETA, SOURCE_PHI=SOURCE_PHI,
                          NU=NU,DELTANU_NU=DELTANU_NU,npts=npts,display=display,HORN_OPEN=all_open_j,
                          background=background)
    S_ij,N_ij=get_fpimage(qubic,SOURCE_POWER=SOURCE_POWER,SOURCE_THETA=SOURCE_THETA, SOURCE_PHI=SOURCE_PHI,
                          NU=NU,DELTANU_NU=DELTANU_NU,npts=npts,display=display,HORN_OPEN=all_open_ij,
                          background=background)

    saturated = (Sall == 0) | (S_i == 0) | (S_j == 0) | (S_ij == 0)
    Sout=Sall+S_ij-S_i-S_j
    Nout=np.sqrt(Nall**2+N_ij**2+N_i**2+N_j**2)
    Sout[saturated]=0
    Nout[saturated]=0
    return Sout,Nout,[Sall,S_i,S_j,S_ij]

def clock_image(time):

    # Create new image
    display = Image.new('RGBA', display_size, color=0)

    # Draw graphics
    draw = ImageDraw.Draw(display)

    # Add hour markers
    for angle in range(0, 360, 30):
        theta = np.radians(angle)
        points = radial_line(centre, theta, r1, r2, w)
        draw.polygon(points, fill=colors[0])

    # Add minute and hour hands
    h1_angle = np.radians(time.minute*360/60)
    h2_angle = np.radians((time.hour*360 + time.minute*6)/12)

    points = radial_line(centre, h1_angle, -w, r3, w)
    draw.polygon(points, fill=colors[1])

    points = radial_line(centre, h2_angle, -w, r4, w+2)
    draw.polygon(points, fill=colors[1])

    return display

    def _get_delta_tdb_tt(self, jd1=None, jd2=None):
        if not hasattr(self, '_delta_tdb_tt'):
            # If jd1 and jd2 are not provided (which is the case for property
            # attribute access) then require that the time scale is TT or TDB.
            # Otherwise the computations here are not correct.
            if jd1 is None or jd2 is None:
                if self.scale not in ('tt', 'tdb'):
                    raise ValueError('Accessing the delta_tdb_tt attribute '
                                     'is only possible for TT or TDB time '
                                     'scales')
                else:
                    jd1 = self._time.jd1
                    jd2 = self._time.jd2

            # First go from the current input time (which is either
            # TDB or TT) to an approximate UTC.  Since TT and TDB are
            # pretty close (few msec?), assume TT.
            njd1, njd2 = sofa_time.tt_tai(jd1, jd2)
            njd1, njd2 = sofa_time.tai_utc(njd1, njd2)
            # TODO: actually need to go to UT1 which needs DUT.
            ut = njd1 + njd2

            # Compute geodetic params needed for d_tdb_tt()
            phi = np.radians(self.lat)
            elon = np.radians(self.lon)
            xyz = sofa_time.iau_gd2gc(1, elon, phi, 0.0)
            u = np.sqrt(xyz[0] ** 2 + xyz[1] ** 2)
            v = xyz[2]

            self._delta_tdb_tt = sofa_time.d_tdb_tt(jd1, jd2, ut, elon, u, v)

        return self._delta_tdb_tt

def lb_to_azalt(gal_l, gal_b, lat, lon, unixtime=None):
    """
    Converts galactic coordiantes to az/alt coordinates. The input
    angles are expected to be in degrees.
    """
    # Convert degrees to radians
    gal_l = _np.radians(gal_l)
    gal_b = _np.radians(gal_b)

    # Location on unit sphere
    theta = _np.pi / 2 - gal_b
    phi   = gal_l
    x = _np.sin(theta) * _np.cos(phi)
    y = _np.sin(theta) * _np.sin(phi)
    z = _np.cos(theta)
    cartesian = _np.array([x, y, z])

    # Get the of-date ra/dec to convert to az/alt
    radec_cart = _np.dot(matrix_lb_radec(), cartesian)
    x, y, z = radec_cart
    r = _np.sqrt(x**2 + y**2 + z**2) # should be 1
    theta = _np.arccos(z / r)
    phi = _np.arctan2(y, x)
    ra = _np.degrees(phi)
    dec = _np.degrees(_np.pi / 2 - theta)
    if ra < 0:
        ra += 360.0
    ra, dec = get_radec_ofdate(ra, dec, unixtime)

    # Convert ra/dec to az/alt (in degrees)
    return radec_to_azalt(ra, dec, lat, lon, unixtime)

def fnc(ra, *args):
    """
    Function used in the minimisation problem.

    :param ra:
        Semi-axis of the ellipses used in the Yu et al.
    :returns:
        The absolute difference between the epicentral distance and the
        adjusted distance
    """
    #
    # epicentral distance
    repi = args[0]
    #
    # azimuth
    theta = args[1]
    #
    # magnitude
    mag = args[2]
    #
    # coefficients
    coeff = args[3]
    #
    # compute the difference between epicentral distances
    rb = rbf(ra, coeff, mag)
    t1 = ra**2 * (np.sin(np.radians(theta)))**2
    t2 = rb**2 * (np.cos(np.radians(theta)))**2
    xx = ra * rb / (t1+t2)**0.5
    return xx-repi

def azalt_to_lb(az, alt, lat, lon, unixtime=None):
    """
    Converts az/alt coordiantes to galactic coordinates. All inputs are
    expected to be in degrees.
    """
    # Get the of-date ra/dec
    ra, dec = azalt_to_radec(az, alt, lat, lon, unixtime)

    # Convert of-date ra/dec to J2000
    ra, dec = get_radec_j2000(ra, dec, unixtime)

    # Convert degrees to radians
    ra  = _np.radians(ra)
    dec = _np.radians(dec)

    # Location on unit sphere
    theta = _np.pi / 2 - dec
    phi   = ra
    x = _np.sin(theta) * _np.cos(phi)
    y = _np.sin(theta) * _np.sin(phi)
    z = _np.cos(theta)
    cartesian = _np.array([x, y, z])

    # Perform the final matrix multilication to get l/b
    lb_cart = _np.dot(matrix_radec_lb(), cartesian)
    x, y, z = lb_cart
    r = _np.sqrt(x**2 + y**2 + z**2) # should be 1
    theta = _np.arccos(z / r)
    phi = _np.arctan2(y, x)
    gal_l = _np.degrees(phi)
    gal_b = _np.degrees(_np.pi / 2 - theta)
    if gal_l < 0:
        gal_l += 360.0

    return (gal_l, gal_b)

    def update(self):

        if not self.targetFrame:
            return

        r = self.properties.getProperty("Distance (m)")
        theta = np.radians(90 - self.properties.getProperty("Elevation (deg)"))
        phi = np.radians(180 - self.properties.getProperty("Azimuth (deg)"))

        x = r * np.cos(phi) * np.sin(theta)
        y = r * np.sin(phi) * np.sin(theta)
        z = r * np.cos(theta)

        c = self.camera
        targetToWorld = self.targetFrame.transform

        currentPosition = np.array(c.GetPosition())
        desiredPosition = np.array(targetToWorld.TransformPoint([x, y, z]))

        smoothTime = self.properties.getProperty("Smooth Time (s)")

        newPosition, self.currentVelocity = smoothDamp(
            currentPosition, desiredPosition, self.currentVelocity, smoothTime, maxSpeed=100, deltaTime=self.dt
        )

        trackerToWorld = transformUtils.getLookAtTransform(targetToWorld.GetPosition(), newPosition)

        c.SetFocalPoint(targetToWorld.GetPosition())
        c.SetPosition(trackerToWorld.GetPosition())
        self.view.render()

def readGC(file,standard_cosmology=True):

    ra,dec,zs,zp =[],[],[],[]
    dl = []
    with open(file,'r') as f:
        if standard_cosmology: omega = lal.CreateCosmologicalParameters(0.7,0.3,0.7,-1.0,0.0,0.0)
        for line in f:
            fields = line.split(None)
            if 0.0 < np.float(fields[40]) > 0.0 or np.float(fields[41]) > 0.0:
                if not(standard_cosmology):
                    h = np.random.uniform(0.5,1.0)
                    om = np.random.uniform(0.0,1.0)
                    ol = 1.0-om
                    omega = lal.CreateCosmologicalParameters(h,om,ol,-1.0,0.0,0.0)
                
                ra.append(np.float(fields[0]))
                dec.append(np.float(fields[1]))
                zs.append(np.float(fields[40]))
                zp.append(np.float(fields[41]))
                if not(np.isnan(zs[-1])):
                    dl.append(lal.LuminosityDistance(omega,zs[-1]))
                elif not(np.isnan(zp[-1])):
                    dl.append(lal.LuminosityDistance(omega,zp[-1]))
                else:
                    dl.append(-1)
        f.close()
    return np.column_stack((np.radians(np.array(ra)),np.radians(np.array(dec)),np.array(dl)))