def _dots_per_unit(self, units):
     """
     Return a scale factor for converting from units to pixels
     """
     ax = self.ax
     if units in ("x", "y", "xy"):
         if units == "x":
             dx0 = ax.viewLim.width
             dx1 = ax.bbox.width
         elif units == "y":
             dx0 = ax.viewLim.height
             dx1 = ax.bbox.height
         else:  # 'xy' is assumed
             dxx0 = ax.viewLim.width
             dxx1 = ax.bbox.width
             dyy0 = ax.viewLim.height
             dyy1 = ax.bbox.height
             dx1 = np.hypot(dxx1, dyy1)
             dx0 = np.hypot(dxx0, dyy0)
         dx = dx1 / dx0
     else:
         if units == "width":
             dx = ax.bbox.width
         elif units == "height":
             dx = ax.bbox.height
         elif units == "dots":
             dx = 1.0
         elif units == "inches":
             dx = ax.figure.dpi
         else:
             raise ValueError("unrecognized units")
     return dx

 def get_error(self, data, model):
     (center, diameter, angle) = model
     dis = npy.hypot(data[:, 0] - center[0], data[:, 1] - center[1])
     theta = npy.arctan2(data[:, 1] - center[1], data[:, 0] - center[0]) - angle
     ellipse_radius = npy.hypot(diameter[0] / 2 * npy.cos(theta), 
                diameter[1] / 2 * npy.sin(theta))
     return npy.abs(dis - ellipse_radius)

    def zoomToAll(self):
        if self.m_nImgs < 1:
            return

        posA=N.array(self.m_imgPosArr)
        sizA=N.array(self.m_imgSizeArr)
        a=N.array([N.minimum.reduce(posA),
                   N.maximum.reduce(posA+sizA),
                   ])
        from .all import U

        MC = N.array([0.5, 0.5]) # mosaic viewer's center (0.5, 0.5)
        a -= MC
        hypot = N.array((N.hypot(a[0][0], a[0][1]),
                         N.hypot(a[1][0], a[1][1])))
        theta = N.array((N.arctan2(a[0][1], a[0][0]),
                         N.arctan2(a[1][1], a[1][0]))) # radians
        phi = theta + U.deg2rad(self.m_rot)
        mimXY = N.array((hypot[0]*N.cos(phi[0]), hypot[0]*N.sin(phi[0])))
        maxXY = N.array((hypot[1]*N.cos(phi[1]), hypot[1]*N.sin(phi[1])))
        a = N.array((mimXY, maxXY))
        a.sort(0)
        if self.m_aspectRatio == -1:
            a = N.array(([a[0][0],-a[1][1]],[a[1][0],-a[0][1]]))

        self.zoomToRect(x0=a[0][0], y0=a[0][1],
                        x1=a[-1][0],y1=a[-1][1])

    def get_bezier_points_at(self, at, grid=256):
        at = np.asarray(at)

        # The Bezier curve is parameterized by a value t which ranges from 0
        # to 1. However, there is a nonlinear relationship between this value
        # and arclength. We want to parameterize by t', which measures
        # normalized arclength. To do this, we have to calculate the function
        # arclength(t), and then invert it.
        t = np.linspace(0, 1, grid)

        x, y = Bezier(list(zip(self._xp, self._yp)), t).T
        x_deltas = np.diff(x)
        y_deltas = np.diff(y)

        arclength_deltas = np.empty(t.shape)
        arclength_deltas[0] = 0

        np.hypot(x_deltas, y_deltas, out=arclength_deltas[1:])
        arclength = np.cumsum(arclength_deltas)
        arclength /= arclength[-1]

        # Now (t, arclength) is a LUT describing the t -> arclength mapping
        # Invert it to get at -> t
        at_t = np.interp(at, arclength, t)

        # And finally look up at the Bezier values at at_t
        # (Might be quicker to np.interp againts x and y, but eh, doesn't
        # really matter.)
        return Bezier(list(zip(self._xp, self._yp)), at_t).T

def plot_winds(m, lat,lon, image, u, v, name, step=16):
    m.pcolormesh(lon, lat, image, latlon=True, cmap='gray')
    
    h,w = lat.shape[:2]
    y,x = np.mgrid[step/2:h:step,step/2:w:step].reshape(2,-1)
    
    #extract data
    u = u[y,x]
    v = v[y,x]
    lon = lon[y,x]
    lat = lat[y,x]
    
    #calculate map positionf lat lons
    x,y = m(lon,lat)
    
    # Calculate the orientation of the vectors
    x1, y1 = m(lon+u, lat+v)
    u_map, v_map = x1-x, y1-y
    
    # Rescale the magnitudes of the vectors...
    mag_scale = np.hypot(u_map, v_map) / np.hypot(u, v)
    u_map /= mag_scale
    v_map /= mag_scale
    
    # Draw barbs
    #m.barbs(x,y,u_map,v_map, length=5, color='red')
    m.quiver(x,y,u_map,v_map,scale=200, color='red')

    # Draw some grid lines for reference
    m.drawparallels(np.arange(-90,90,2),labels=[1,1,0,0])
    m.drawmeridians(np.arange(0,360,2),labels=[0,0,0,1])
    plt.savefig(name,bbox_inches='tight')
    plt.close()

def get_impact_parameters_cpair(cpair,ra0,dec0,cosmo,npoints=1000):
    """Returns the impact parameter between the straight line given by a
    cluster-pair and (ra0,dec0), as well as the distance along the
    cpair axis to the closest cluster of the pair.

    """
    
    ra1  = cpair['ra1']
    dec1 = cpair['dec1']
    z1   = cpair['z1']
    ra2  = cpair['ra2']
    dec2 = cpair['dec2']
    z2   = cpair['z2']
    zmed = cpair['redshift']
    
    ra_mpc,dec_mpc,sep_mpc,dv = get_line_radec_sepdv(ra1,dec1,z1,ra2,dec2,z2,ra0,dec0,zmed,cosmo,npoints=npoints)
    
    #get hypotenuse
    sep = np.hypot(ra_mpc,dec_mpc)

    #get minimum distance
    sep_min = np.min(sep)
    
    #lets add the distance along the cpair axis to the closest cluster
    #find the closest cluster distance first
    sep_cl_mpc1 = np.hypot(ra_mpc[0],dec_mpc[0])
    sep_cl_mpc2 = np.hypot(ra_mpc[-1],dec_mpc[-1])
    sep_cl_mpc = np.min([sep_cl_mpc1.value,sep_cl_mpc2.value])
    #then, project along the cluster axis 
    sep_x_mpc = np.sqrt(sep_cl_mpc**2 - sep_min.value**2)
    
    return sep_min.value, sep_x_mpc

def go(start, stores, items, price_of_gas, perishable, solved):
  if (start, str(sorted(items)), perishable) in solved:
    return solved[(start, str(sorted(items)), perishable)]
  gas_home = price_of_gas * hypot(0 - start[0], 0 - start[1])
  if len(items) == 0:
    return gas_home 
  if start == (0,0):
    perishable = 0
  
  costs = list()
 
  if perishable == 1:
    costs.append(gas_home + go((0,0), stores, items, price_of_gas, 0, solved))
    filtered_stores = [x for x in stores if x[0] == start]
  else:
    filtered_stores = stores

  for i in range(len(items)):
    p = perishable
    if items[i][-1] == '!':
      p = 1
      item = items[i][:-1]
    else:
      item = items[i]

    for j in range(len(filtered_stores)):
      filtered_items = [x for x in filtered_stores[j][1] if item in x]
      if len(filtered_items) == 1:
        cost = int(filtered_items[0].split(":")[1])
        gas_here = price_of_gas * hypot(start[0] - filtered_stores[j][0][0], start[1] - filtered_stores[j][0][1])
        costs.append( cost + gas_here + go(filtered_stores[j][0], stores, [x for x in items if item not in x], price_of_gas, p, solved))
  solved[(start,str(sorted(items)),perishable)] = min(costs)
  return min(costs)

def get_plotsense_dict(file, NTOP=None, NANTS=None):
	antpos = np.loadtxt(file)
	if NANTS is None:
		NANTS = antpos.shape[0]

	bl_gps = {}
	for i in xrange(NANTS-1):
		a1pos = antpos[i]
		for j in xrange(i+1,NANTS):
			a2pos = antpos[j]
			blx, bly, blz = a2pos-a1pos
			has_entry = False
			for key in bl_gps.keys():
				if np.hypot(key[0]-blx, key[1]-bly) < 1:
					has_entry = True
					bl_gps[key] = bl_gps.get(key, []) + [(i,j)]
					break
			if not has_entry:
				bl_gps[(blx, bly)] = [(i,j)]

	n_unique = len(bl_gps.keys())
	n_total = NANTS*(NANTS-1)/2
	
	print "Found %d classes among %d total baselines" % (n_unique, n_total)
	print "Sorting dictionary"
	sorted_keys = sorted(bl_gps.keys(), key=lambda x: np.hypot(*x))
	if NTOP is None: 
		NTOP = n_unique
	key_pool = np.arange(NTOP)
	top_dict = {}
	for i, key in enumerate(sorted_keys[:NTOP]):
		mult = len(bl_gps[key])
		label = str(bl_gps[key][0][0])+'_'+str(bl_gps[key][0][1])
		top_dict[label] = ((key[0], key[1], 0.), mult) #dict[example_bl] = ((blx, bly, blz), multiplicity)
	return top_dict, bl_gps

def solve_for_nearest( px,py,rx,ry ):
  dpx = polyder(px)
  dpy = polyder(py)
  cp = polymul( dpx, px ) + polymul( dpy, py )
  cp = polyadd( cp, -rx*dpx )      
  cp = polyadd( cp, -ry*dpy )
  t = roots(cp)
  t = real(t[isreal(t)])
  t = t[ (t>=0) * (t<=1) ]

  ##tt = linspace(0,1,100)
  ##from pylab import plot
  ##plot( polyval(px,tt), polyval(py,tt), 'k', hold = 0 )
  ##plot( [rx],[ry], 'r.' )
  ##plot( polyval(px,t[isreal(t)*(real(t)>=0)*(real(t)<=1)]),
  ##      polyval(py,t[isreal(t)*(real(t)>=0)*(real(t)<=1)]), 'o' )
  ##pdb.set_trace()

  if len(t):
    if len(t) == 1:
      return t[0]
    else:
      ux = polyval( px, t )
      uy = polyval( py, t )
      d = hypot( ux - rx, uy - ry )
      return t[ d==d.min() ][0]
  else:
    t = array([0.0,1.0])
    ux = polyval( px, t )
    uy = polyval( py, t )
    d = hypot( ux - rx, uy - ry )
    if d[0] < d[1]:
      return 0.0
    else:
      return 1.0

def norm_dist(src1, src2):
    """
    Calculate the normalised distance between two sources.
    Sources are elliptical Gaussians.

    The normalised distance is calculated as the GCD distance between the centers,
    divided by quadrature sum of the radius of each ellipse along a line joining the two ellipses.

    For ellipses that touch at a single point, the normalized distance will be 1/sqrt(2).

    Parameters
    ----------
    src1, src2 : object
        The two positions to compare. Objects must have the following parameters: (ra, dec, a, b, pa).

    Returns
    -------
    dist: float
        The normalised distance.

    """
    if np.all(src1 == src2):
        return 0
    dist = gcd(src1.ra, src1.dec, src2.ra, src2.dec) # degrees

    # the angle between the ellipse centers
    phi = bear(src1.ra, src1.dec, src2.ra, src2.dec) # Degrees
    # Calculate the radius of each ellipse along a line that joins their centers.
    r1 = src1.a*src1.b / np.hypot(src1.a * np.sin(np.radians(phi - src1.pa)),
                                  src1.b * np.cos(np.radians(phi - src1.pa)))
    r2 = src2.a*src2.b / np.hypot(src2.a * np.sin(np.radians(180 + phi - src2.pa)),
                                  src2.b * np.cos(np.radians(180 + phi - src2.pa)))
    R = dist / (np.hypot(r1, r2) / 3600)
    return R

    def sky2pix_ellipse(self, pos, a, b, pa):
        """
        Convert an ellipse from sky to pixel corrds
        a/b vectors are calculated at an origin pos=(ra,dec)
        All input parameters are in degrees
        Output parameters are:
        x,y - the x,y pixels corresponding to the ra/dec position
        sx, sy - the major minor axes (FWHM) in pixels
        theta - the position angle in degrees

        :param pos: [ra,dec] of the ellipse center
        :param a: major axis
        :param b: minor axis
        :param pa: position angle
        :return: x, y, sx, sy, theta
        """
        ra, dec = pos
        x, y = self.sky2pix(pos)

        x_off, y_off = self.sky2pix(translate(ra, dec, a, pa))
        sx = np.hypot((x - x_off),(y - y_off))
        theta = np.arctan2((y_off - y), (x_off - x))

        x_off, y_off = self.sky2pix(translate(ra, dec, b, pa-90))
        sy = np.hypot((x - x_off), (y - y_off))
        theta2 = np.arctan2((y_off - y), (x_off - x)) - np.pi/2

        # The a/b vectors are perpendicular in sky space, but not always in pixel space
        # so we have to account for this by calculating the angle between the two vectors
        # and modifying the minor axis length
        defect = theta - theta2
        sy *= abs(np.cos(defect))

        return x, y, sx, sy, np.degrees(theta)

def plateASTundo(ast, p, q):
    """ Map gnomonic coordinates to theta, phi. """

    M = ast.M

    # Calculate beta and gamma
    bet = np.arcsin(np.hypot(p, q))
    gam = np.arctan2(q, p)

    if bet > np.pi:
        return False

    # Init vector v
    v = np.zeros(3)

    v[0] = np.sin(bet)*np.cos(gam)
    v[1] = np.sin(bet)*np.sin(gam)
    v[2] = np.cos(bet)

    # Calculate vector u
    u = np.zeros(3)        

    u[0] = M[0,0]*v[0] + M[0,1]*v[1] + M[0,2]*v[2]
    u[1] = M[1,0]*v[0] + M[1,1]*v[1] + M[1,2]*v[2]
    u[2] = M[2,0]*v[0] + M[2,1]*v[1] + M[2,2]*v[2]

    # Convert to theta, phi
    th  = np.arctan2(np.hypot(u[0], u[1]), u[2])
    phi = np.arctan2(u[1], u[0])

    return th, phi

def getDr(d, IO, type_):
    """Compute distortion with given radial distance."""
    # Define radial distance range
    xp = np.linspace(0, d, d*50)
    yp = 0

    # Compute distortion corrections
    x0 = IO["x0"]
    y0 = IO["y0"]
    xbar = xp - x0
    ybar = yp - y0
    r = np.hypot(xbar, ybar)

    if type_ == "symmetric":
        k1 = IO["k1"]
        k2 = IO["k2"]
        k3 = IO["k3"]
        dx = xbar * (r**2 * k1 + r**4 * k2 + r**6 * k3)
        dy = ybar * (r**2 * k1 + r**4 * k2 + r**6 * k3)
    elif type_ == "decentering":
        p1 = IO["p1"]
        p2 = IO["p2"]
        dx = (p1 * (r**2 + 2 * xbar**2) + 2 * p2 * xbar * ybar)
        dy = (2 * p1 * xbar * ybar + p2 * (r**2 + 2 * ybar**2))

    dr = np.hypot(dx, dy)

    return xp, dr

def get_dH2(lab1, lab2):
    """squared hue difference term occurring in deltaE_cmc and deltaE_ciede94

    Despite its name, "dH" is not a simple difference of hue values.  We avoid
    working directly with the hue value, since differencing angles is
    troublesome.  The hue term is usually written as:
        c1 = sqrt(a1**2 + b1**2)
        c2 = sqrt(a2**2 + b2**2)
        term = (a1-a2)**2 + (b1-b2)**2 - (c1-c2)**2
        dH = sqrt(term)

    However, this has poor roundoff properties when a or b is dominant.
    Instead, ab is a vector with elements a and b.  The same dH term can be
    re-written as:
        |ab1-ab2|**2 - (|ab1| - |ab2|)**2
    and then simplified to:
        2*|ab1|*|ab2| - 2*dot(ab1, ab2)
    """
    lab1 = np.asarray(lab1)
    lab2 = np.asarray(lab2)
    a1, b1 = np.rollaxis(lab1, -1)[1:3]
    a2, b2 = np.rollaxis(lab2, -1)[1:3]

    # magnitude of (a, b) is the chroma
    C1 = np.hypot(a1, b1)
    C2 = np.hypot(a2, b2)

    term = (C1 * C2) - (a1 * a2 + b1 * b2)
    return 2 * term

def cart2sph(x, y, z):
    """Converts cartesian coordinates x, y, z to spherical coordinates az, el, r."""
    hxy = _np.hypot(x, y)
    r = _np.hypot(hxy, z)
    el = _np.arctan2(z, hxy)
    az = _np.arctan2(y, x)
    return az, el, r

def bending_angle(xc,yc,x1,y1,x2,y2):

    # Compute the bending angle (in radians) between three points in pixel space

    '''
    Points are:

        - xc,yc:  x,y center of IR counterpart
        - x1,y1:  x,y center of 1st radio lobe
        - x2,y2:  x,y center of 2nd radio lobe
    '''

    r1 = np.array([x1,y1])
    r2 = np.array([x2,y2])
    center = np.array([xc,yc])

    r1diff = r1-center
    r2diff = r2-center

    r1len = np.hypot(r1diff[0],r1diff[1])
    r2len = np.hypot(r2diff[0],r2diff[1])

    alpha = np.arccos(np.dot(r1diff,r2diff) / (r1len*r2len))

    return alpha

def _cartesian_to_sphere(x, y, z):
    """Transform cartesian coordinates to spherical"""
    hypotxy = np.hypot(x, y)
    r = np.hypot(hypotxy, z)
    elev = np.arctan2(z, hypotxy)
    az = np.arctan2(y, x)
    return az, elev, r

def position_angle(xc,yc,x1,y1,x2,y2):

    # Compute the position angle (in radians, with respect to north) between three points in pixel space

    '''
    Points are:

        - xc,yc:    x,y center of bending angle
        - x1,y1:  x,y center of 1st component
        - x2,y2:  x,y center of 2nd component
    '''

    r1 = np.array([x1,y1])
    r2 = np.array([x2,y2])
    center = np.array([xc,yc])

    r12sum = (r1-center) + (r2-center)
    r12len = np.hypot(r12sum[0],r12sum[1])

    north = np.array([0,1])
    northlen = np.hypot(north[0],north[1])

    alpha = np.arccos(np.dot(r12sum,north) / (r12len*northlen))

    # Measure CCW from north

    if r12sum[0] > 0.:
        alpha = 2*np.pi - alpha

    return alpha

def averageNSLSdata (ee, mc, mb, dc, db, n) :
    i = ee.size / n # total number of energy points
    e = ee[0::n]

    # [0,:] is the first energy point
    # Reshape to put energy on one axis and readings on the other.
    mc.shape = ((i, n))
    dc.shape = ((i, n))
    mb.shape = ((i, n))
    db.shape = ((i, n))

    # Average over readings.
    mca = np.average (mc, axis=1) 
    dca = np.average (dc, axis=1) 
    mba = np.average (mb, axis=1) 
    dba = np.average (db, axis=1) 
    # Note that the average might not be best for cases where there 
    # are outliers. Perhaps a median of some kind? 

    # Find the standard deviation of detector current measurements.
    mcs = np.std (mc, axis=1)
    dcs = np.std (dc, axis=1)
    mbs = np.std (mb, axis=1)
    dbs = np.std (db, axis=1)
    ms = np.hypot (mcs, mbs)
    ds = np.hypot (dcs, dbs)
    return e, mca, mba, ms, dca, dba, ds

def edgedist(pixel, shape=(512,512)):
    x = pixel[0]
    y= pixel[1]
    circumference = 2.0*shape[0] + 2.0*shape[1]
    totaldist = np.sum([np.hypot(x-a, y-b) for a in range(shape[0]) for b in [0,shape[1]]])
    totaldist2 = np.sum([np.hypot(x-a, y-b) for a in [0,shape[0]] for b in range(shape[1])])
    return (totaldist + totaldist2) / circumference, circumference