def get_center(self):
     """
     Return the RA, Dec of the center of the circle bounding
     this trixel (RA, Dec both in degrees)
     """
     ra, dec = sphericalFromCartesian(self.bounding_circle[0])
     return np.degrees(ra), np.degrees(dec)

def star(a,b,c,alpha,beta,gamma):
    "Calculate unit cell volume, reciprocal cell volume, reciprocal lattice parameters"
    alpha=np.radians(alpha)
    beta=np.radians(beta)
    gamma=np.radians(gamma)
    V=2*a*b*c*\
        np.sqrt(np.sin((alpha+beta+gamma)/2)*\
               np.sin((-alpha+beta+gamma)/2)*\
               np.sin((alpha-beta+gamma)/2)*\
               np.sin((alpha+beta-gamma)/2))
    Vstar=(2*np.pi)**3/V;
    astar=2*np.pi*b*c*np.sin(alpha)/V
    bstar=2*np.pi*a*c*np.sin(beta)/V
    cstar=2*np.pi*b*a*np.sin(gamma)/V
    alphastar=np.arccos((np.cos(beta)*np.cos(gamma)-\
                        np.cos(alpha))/ \
                       (np.sin(beta)*np.sin(gamma)))
    betastar= np.arccos((np.cos(alpha)*np.cos(gamma)-\
                        np.cos(beta))/ \
                       (np.sin(alpha)*np.sin(gamma)))
    gammastar=np.arccos((np.cos(alpha)*np.cos(beta)-\
                        np.cos(gamma))/ \
                       (np.sin(alpha)*np.sin(beta)))
    V=V
    alphastar=np.degrees(alphastar)
    betastar=np.degrees(betastar)
    gammastar=np.degrees(gammastar)
    return astar,bstar,cstar,alphastar,betastar,gammastar

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 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 _lambet(self):
     """Lower overhead ecliptic longitude and latitude. [deg]"""
     from astropy.coordinates import Angle
     lam = np.arctan2(self._rot.T[1], self._rot.T[0])
     bet = np.arctan2(self._rot.T[2],
                      np.sqrt(self._rot.T[0]**2 + self._rot.T[1]**2))
     return np.degrees(lam), np.degrees(bet)

    def test_degrees(self):

        # the following doesn't make much sense in terms of the name of the
        # routine, but we check it gives the correct result.
        q1 = np.rad2deg(60. * u.degree)
        assert_allclose(q1.value, 60.)
        assert q1.unit == u.degree

        q2 = np.degrees(60. * u.degree)
        assert_allclose(q2.value, 60.)
        assert q2.unit == u.degree

        q3 = np.rad2deg(np.pi * u.radian)
        assert_allclose(q3.value, 180.)
        assert q3.unit == u.degree

        q4 = np.degrees(np.pi * u.radian)
        assert_allclose(q4.value, 180.)
        assert q4.unit == u.degree

        with pytest.raises(TypeError):
            np.rad2deg(3. * u.m)

        with pytest.raises(TypeError):
            np.degrees(3. * u.m)

def waypoints_to_commands(coords):
	#cmd = [[vx,az,time],etc]
	#Convert waypoints to value in stage
	lin_vel = 0.2
	ang_vel = math.radians(45)    #45 deg/s in rad/s
	init_ang = 0;
	move_ang = [0]
	move_dist = [0]
	for i in range(len(coords)-1):
		p1 = coords[i]
		p2 = coords[i+1]
		move_ang.append(math.atan2(p2[1]-p1[1],p2[0]-p1[0]))
		move_dist.append(math.sqrt((p2[1]-p1[1])**2+(p2[0]-p1[0])**2))

	print np.degrees(move_ang)
	print len(move_dist)
	move_cmd = []

	for i in range(len(move_ang)-1):
		ang_cmd = (move_ang[i+1]-move_ang[i])
		ang_time = ang_cmd/ang_vel
		dist_cmd =move_dist[i+1]-move_dist[i]
		dist_time = dist_cmd/lin_vel
		move_cmd.append([0,np.sign(ang_cmd),math.fabs(ang_time)])
		move_cmd.append([np.sign(dist_cmd),0,math.fabs(dist_time)])

	print move_cmd
	print len(move_cmd)
	return move_cmd

def sendSearchRequest(matches):
    ra=list(np.degrees(matches[matches['isinobs']==False]['can_ra']))
    dec=list(np.degrees(matches[matches['isinobs']==False]['can_dec']))
    bands='[g,r,i,z]'
    req='http://desdev3.cosmology.illinois.edu:8000/api?username=lzullo&password=lzu70chips&ra=%s&dec=%s&bands=%s' % (ra,dec,bands)
    submit = requests.get(req)
    return submit.json()['job']

 def draw(self):
     gL.glPushMatrix()
     if self.b:
         self.draw_rod(self.b)
         gL.glTranslatef(0, 0, self.b)
     gL.glRotatef(degrees(self.gamma), 0, 0, 1)
     if self.d:
         gL.glPushMatrix()
         gL.glRotatef(90, 0, 1, 0)
         self.draw_rod(self.d)
         gL.glPopMatrix()
         gL.glTranslatef(self.d, 0, 0)
     gL.glRotatef(degrees(self.alpha), 1, 0, 0)
     if self.r:
         self.draw_rod(self.r)
         gL.glTranslatef(0, 0, self.r)
     gL.glRotatef(degrees(self.theta), 0, 0, 1)
     if self.shift:
         gL.glPushMatrix()
         shift = self.shift * self.length
         self.draw_rod(shift)
         gL.glTranslatef(0, 0, shift)
         self.draw_joint()
         gL.glPopMatrix()
     else:
         self.draw_joint()
     for child in self.children:
         child.draw()
     gL.glPopMatrix()

def GetHealPixRectangles(nside, dbrange, nest):
    hpindex = np.arange(hp.nside2npix(nside))

    vec_corners = hp.boundaries(nside, hpindex, nest=nest)
    vec_corners = np.transpose(vec_corners, (0,2,1))
    vec_corners = np.reshape(vec_corners, (vec_corners.shape[0]*vec_corners.shape[1], vec_corners.shape[2]))
   
    theta_corners, phi_corners = hp.vec2ang(vec_corners)
    theta_corners = np.reshape(theta_corners, (theta_corners.shape[0]/4, 4))
    phi_corners = np.reshape(phi_corners, (phi_corners.shape[0]/4, 4))

    ra_corners = np.degrees(phi_corners)
    dec_corners = 90.0 - np.degrees(theta_corners)

    rainside = ( (ra_corners > dbrange[0]) & (ra_corners < dbrange[1]) )
    rakeep = np.sum(rainside, axis=-1)
    decinside = ( (dec_corners > dbrange[2]) & (dec_corners < dbrange[3]) )
    deckeep = np.sum(decinside, axis=-1)
    keep = ( (rakeep > 0) & (deckeep > 0) )
    ra_corners, dec_corners, hpindex = Cut(ra_corners, dec_corners, hpindex, cut=keep)

    ramin = np.amin(ra_corners, axis=-1)
    ramax = np.amax(ra_corners, axis=-1)
    decmin = np.amin(dec_corners, axis=-1)
    decmax = np.amax(dec_corners, axis=-1)

    return ramin, ramax, decmin, decmax, hpindex

def Jacobsen(h1, Xm_1, h2, Xm_2):
    alp = np.degrees(np.arccos(np.dot(h1, h2) /
                               (np.linalg.norm(h1) * np.linalg.norm(h2))))
    bet = np.degrees(np.arccos(np.dot(Xm_1, Xm_2) /
                               (np.linalg.norm(Xm_1) * np.linalg.norm(Xm_2))))
    if ((alp - bet)**2) > 1:
        print('check your indexing!')

    a = 3.567  # diamond lattice parameter
    # recip lattice par(note this is the mantid convention: no 2 pi)
    ast = (2 * np.pi) / a
    B = np.array([[ast, 0, 0], [0, ast, 0], [0, 0, ast]])
    Xm_g = np.cross(Xm_1, Xm_2)
    Xm = np.column_stack([Xm_1, Xm_2, Xm_g])

    # Vector Q1 is described in reciprocal space by its coordinate matrix h1
    Xa_1 = B.dot(h1)
    Xa_2 = B.dot(h2)
    Xa_g = np.cross(Xa_1, Xa_2)
    Xa = np.column_stack((Xa_1, Xa_2, Xa_g))

    R = Xa.dot(np.linalg.inv(Xm))
    U = np.linalg.inv(R)

    UB = U.dot(B)

    return UB

def cone(plunge, bearing, angle, segments=100):
    """
    Calculates the longitude and latitude of the small circle (i.e. a cone)
    centered at the given *plunge* and *bearing* with an apical angle of
    *angle*, all in degrees.

    Parameters
    ----------
    plunge : number or sequence of numbers
        The plunge of the center of the cone(s) in degrees. The plunge is
        measured in degrees downward from the end of the feature specified by
        the bearing.
    bearing : number or sequence of numbers
        The bearing (azimuth) of the center of the cone(s) in degrees.
    angle : number or sequence of numbers
        The apical angle (i.e. radius) of the cone(s) in degrees.
    segments : int, optional
        The number of vertices in the small circle.

    Returns
    -------
    lon, lat : arrays
        `num_measurements` x `num_segments` arrays of longitude and latitude in
        radians.
    """
    plunges, bearings, angles = np.atleast_1d(plunge, bearing, angle)
    lons, lats = [], []
    for plunge, bearing, angle in zip(plunges, bearings, angles):
        lat = (90 - angle) * np.ones(segments, dtype=float)
        lon = np.linspace(-180, 180, segments)
        lon, lat = _rotate(lon, lat, -plunge, axis='y')
        lon, lat = _rotate(np.degrees(lon), np.degrees(lat), bearing, axis='x')
        lons.append(lon)
        lats.append(lat)
    return np.vstack(lons), np.vstack(lats)

    def run(self, stars, visits, **kwargs):
        # XXX-Double check extinction is close to the Opsim transparency
        extinc_mags = visits['transparency']
        if extinc_mags != 0.:
            # need to decide on how to get extinc_mags from Opsim
            # Maybe push some of these params up to be setable?
            SFtheta, SFsf = self.SF.CloudSf(500., 300., 5., extinc_mags, .55)
            # Call the Clouds
            self.cloud.makeCloudImage(SFtheta,SFsf,extinc_mags, fov=self.fov)
            # Interpolate clouds to correct position.  Nearest neighbor for speed?
            nim = self.cloud.cloudimage[0,:].size
            # calc position in cloud image of each star
            starx_interp = (np.degrees(stars['x']) + self.fov/2.)*3600./ self.cloud.pixscale
            stary_interp = (np.degrees(stars['y']) + self.fov/2.)*3600./ self.cloud.pixscale

            # Round off position and make it an int
            starx_interp = np.round(starx_interp).astype(int)
            stary_interp = np.round(stary_interp).astype(int)

            # Handle any stars that are out of the field for some reason
            starx_interp[np.where(starx_interp < 0)] = 0
            starx_interp[np.where(starx_interp > nim-1)] = nim-1
            stary_interp[np.where(stary_interp < 0)] = 0
            stary_interp[np.where(stary_interp > nim-1)] = nim-1

            dmag = self.cloud.cloudimage[starx_interp,stary_interp]
        else:
            dmag = np.zeros(stars.size)
        return dmag

def plot_lm(d, snrs, l1s, m1s, outroot):
    """ Plot the lm coordinates (relative to phase center) for all candidates.
    """

    outname = os.path.join(d["workdir"], "plot_" + outroot + "_impeak.png")

    snrmin = 0.8 * min(d["sigma_image1"], d["sigma_image2"])
    fig4 = plt.Figure(figsize=(10, 10))
    ax4 = fig4.add_subplot(111)

    # plot positive
    good = n.where(snrs > 0)
    sizes = (snrs[good] - snrmin) ** 5  # set scaling to give nice visual sense of SNR
    xarr = 60 * n.degrees(l1s[good])
    yarr = 60 * n.degrees(m1s[good])
    ax4.scatter(xarr, yarr, s=sizes, facecolor="none", alpha=0.5, clip_on=False)
    # plot negative
    good = n.where(snrs < 0)
    sizes = (n.abs(snrs[good]) - snrmin) ** 5  # set scaling to give nice visual sense of SNR
    xarr = 60 * n.degrees(l1s[good])
    yarr = 60 * n.degrees(m1s[good])
    ax4.scatter(xarr, yarr, s=sizes, marker="x", edgecolors="k", alpha=0.5, clip_on=False)

    ax4.set_xlabel("Dec Offset (amin)")
    ax4.set_ylabel("RA Offset (amin)")
    fov = n.degrees(1.0 / d["uvres"]) * 60.0
    ax4.set_xlim(fov / 2, -fov / 2)
    ax4.set_ylim(-fov / 2, fov / 2)
    canvas4 = FigureCanvasAgg(fig4)
    canvas4.print_figure(outname)

  def pitch_roll(self, px, pz):
    """works out the pitch (rx) and roll (rz) to apply to an object
    on the surface of the map at this point

    * returns a tuple (pitch, roll) in degrees

    Arguments:
      *px*
        x location
      *pz*
        z location
    """
    px -= self.unif[0]
    pz -= self.unif[2]
    halfw = self.width/2.0
    halfd = self.depth/2.0
    dx = self.width/self.ix
    dz = self.depth/self.iy
    x0 = int(math.floor((halfw + px)/dx + 0.5))
    if x0 < 0: x0 = 0
    if x0 > self.ix-1: x0 = self.ix-1
    z0 = int(math.floor((halfd + pz)/dz + 0.5))
    if z0 < 0: z0 = 0
    if z0 > self.iy-1: z0 = self.iy-1
    normp = array(self.buf[0].normals[z0*self.ix + x0])
    # slight simplification to working out cross products as dirctn always 0,0,1
    #sidev = cross(normp, dirctn)
    sidev = array([normp[1], -normp[0], 0.0])
    sidev = sidev / sqrt(sidev.dot(sidev))
    #forwd = cross(sidev, normp)
    forwd = array([-normp[2]*normp[0], -normp[2]*normp[1],
                  normp[0]*normp[0] + normp[1]*normp[1]])
    forwd = forwd / sqrt(forwd.dot(forwd))
    return (degrees(arcsin(-forwd[1])), degrees(arctan2(sidev[1], normp[1])))

    def setUp(self):
        self.metadata={}

        #below are metadata values that need to be set in order for
        #get_getFocalPlaneCoordinates to work.  If we had been querying the database,
        #these would be set to meaningful values.  Because we are generating
        #an artificial set of inputs that must comport to the baseline SLALIB
        #inputs, these are set arbitrarily by hand
        self.metadata['pointingRA'] = (numpy.radians(200.0), float)
        self.metadata['pointingDec'] = (numpy.radians(-30.0), float)
        self.metadata['Opsim_rotskypos'] = (1.0, float)

        # these were the LSST site parameters as coded when this unit test was written
        self.test_site=Site(longitude=numpy.degrees(-1.2320792),
                            latitude=numpy.degrees(-0.517781017),
                            height=2650.0,
                            temperature=11.505,
                            pressure=749.3,
                            lapseRate=0.0065,
                            humidity=0.4)

        self.obs_metadata=ObservationMetaData(mjd=50984.371741,
                                     boundType='circle',
                                     boundLength=0.05,
                                     phoSimMetaData=self.metadata,
                                     site=self.test_site)

        self.tol=1.0e-5

 def proj(lons, lats, reverse=False):
     if not reverse:
         lambdas, phis = numpy.radians(lons), numpy.radians(lats)
         cos_phis = numpy.cos(phis)
         lambdas -= lambda0
         # calculate the sine of the distance between projection center
         # and each of the points to project
         sin_dist = numpy.sqrt(
             numpy.sin((phi0 - phis) / 2.0) ** 2.0
             + cos_phi0 * cos_phis * numpy.sin(lambdas / 2.0) ** 2.0
         )
         if (sin_dist > sin_pi_over_4).any():
             raise ValueError('some points are too far from the projection '
                              'center lon=%s lat=%s' %
                              (numpy.degrees(lambda0), numpy.degrees(phi0)))
         xx = numpy.cos(phis) * numpy.sin(lambdas)
         yy = cos_phi0 * numpy.sin(phis) \
              - sin_phi0 * cos_phis * numpy.cos(lambdas)
         return xx * EARTH_RADIUS, yy * EARTH_RADIUS
     else:
         # "reverse" mode, arguments are actually abscissae
         # and ordinates in 2d space
         xx, yy = lons / EARTH_RADIUS, lats / EARTH_RADIUS
         cos_c = numpy.sqrt(1 - (xx ** 2 + yy ** 2))
         phis = numpy.arcsin(cos_c * sin_phi0 + yy * cos_phi0)
         lambdas = numpy.arctan2(xx, cos_phi0 * cos_c - yy * sin_phi0)
         return numpy.degrees(lambda0 + lambdas), numpy.degrees(phis)

def makeObservationMetaData():
    #create a cartoon ObservationMetaData object
    mjd = 52000.0
    alt = numpy.pi/2.0
    az = 0.0
    band = 'r'
    testSite = Site(latitude=numpy.degrees(0.5), longitude=numpy.degrees(1.1), height=3000,
                    temperature=260.0, pressure=725.0, lapseRate=0.005, humidity=0.4)
    obsTemp = ObservationMetaData(site=testSite, mjd=mjd)
    centerRA, centerDec = _raDecFromAltAz(alt, az, obsTemp)
    rotTel = _getRotTelPos(centerRA, centerDec, obsTemp, 0.0)

    obsDict = calcObsDefaults(centerRA, centerDec, alt, az, rotTel, mjd, band,
                 testSite.longitude_rad, testSite.latitude_rad)

    obsDict['Opsim_expmjd'] = mjd
    radius = 0.1
    phoSimMetaData = OrderedDict([
                      (k, (obsDict[k],numpy.dtype(type(obsDict[k])))) for k in obsDict])

    obs_metadata = ObservationMetaData(boundType='circle', boundLength=2.0*radius,
                                       phoSimMetaData=phoSimMetaData, site=testSite)



    return obs_metadata

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 raDecFromVec(v):
    """
    Taken from
    http://www.math.montana.edu/frankw/ccp/multiworld/multipleIVP/spherical/learn.htm
    Search for "convert from Cartestion to spherical coordinates"

    Adapted because I'm dealing with declination which is defined
    with 90degrees at zenith
    """

    # Ensure v is a normal vector
    v /= np.linalg.norm(v)

    ra_deg = 0  # otherwise not in namespace0
    dec_rad = np.arcsin(v[2])
    s = np.hypot(v[0], v[1])
    if s == 0:
        ra_rad = 0
    else:
        ra_rad = np.arcsin(v[1] / s)
        ra_deg = np.degrees(ra_rad)
        if v[0] >= 0:
            if v[1] >= 0:
                pass
            else:
                ra_deg = 360 + ra_deg
        else:
            if v[1] > 0:
                ra_deg = 180 - ra_deg
            else:
                ra_deg = 180 - ra_deg

    raDec = ra_deg, np.degrees(dec_rad)
    return np.array(raDec)