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 hgs_to_hcc(heliogcoord, heliocframe):
    """
    Convert from Heliographic Stonyhurst to Heliograpic Carrington.
    """
    hglon = heliogcoord.lon
    hglat = heliogcoord.lat
    r = heliogcoord.radius.to(u.m)

    l0b0_pair = [heliocframe.L0, heliocframe.B0]

    l0_rad = l0b0_pair[0].to(u.rad)
    b0_deg = l0b0_pair[1]

    lon = np.deg2rad(hglon)
    lat = np.deg2rad(hglat)

    cosb = np.cos(b0_deg.to(u.rad))
    sinb = np.sin(b0_deg.to(u.rad))

    lon = lon - l0_rad

    cosx = np.cos(lon)
    sinx = np.sin(lon)
    cosy = np.cos(lat)
    siny = np.sin(lat)

    x = r * cosy * sinx
    y = r * (siny * cosb - cosy * cosx * sinb)
    zz = r * (siny * sinb + cosy * cosx * cosb)

    representation = CartesianRepresentation(x.to(u.km), y.to(u.km), zz.to(u.km))
    return heliocframe.realize_frame(representation)

    def to_SQL(self, RAname, DECname):

        if self.DECdeg != 90.0 and self.DECdeg != -90.0:
            RAmax = self.RAdeg + \
            360.0 * np.arcsin(np.sin(0.5*self.radius) / np.cos(self.DEC))/np.pi
            RAmin = self.RAdeg - \
            360.0 * np.arcsin(np.sin(0.5*self.radius) / np.cos(self.DEC))/np.pi
        else:
           #just in case, for some reason, we are looking at the poles
           RAmax = 360.0
           RAmin = 0.0

        DECmax = self.DECdeg + self.radiusdeg
        DECmin = self.DECdeg - self.radiusdeg

        #initially demand that all objects are within a box containing the circle
        #set from the DEC1=DEC2 and RA1=RA2 limits of the haversine function
        bound = ("%s between %f and %f and %s between %f and %f "
                     % (RAname, RAmin, RAmax, DECname, DECmin, DECmax))

        #then use the Haversine function to constrain the angular distance form boresite to be within
        #the desired radius.  See http://en.wikipedia.org/wiki/Haversine_formula
        bound = bound + ("and 2 * ASIN(SQRT( POWER(SIN(0.5*(%s - %s) * PI() / 180.0),2)" % (DECname,self.DECdeg))
        bound = bound + ("+ COS(%s * PI() / 180.0) * COS(%s * PI() / 180.0) " % (DECname, self.DECdeg))
        bound = bound + ("* POWER(SIN(0.5 * (%s - %s) * PI() / 180.0),2)))" % (RAname, self.RAdeg))
        bound = bound + (" < %s " % self.radius)

        return bound

    def get_gabors(self, rf):
        lams =  float(rf[0])/self.sfs # lambda = 1./sf  #1./np.array([.1,.25,.4])
        sigma = rf[0]/2./np.pi
        # rf = [100,100]
        gabors = np.zeros(( len(oris),len(phases),len(lams), rf[0], rf[1] ))

        i = np.arange(-rf[0]/2+1,rf[0]/2+1)
        #print i
        j = np.arange(-rf[1]/2+1,rf[1]/2+1)
        ii,jj = np.meshgrid(i,j)
        for o, theta in enumerate(self.oris):
            x = ii*np.cos(theta) + jj*np.sin(theta)
            y = -ii*np.sin(theta) + jj*np.cos(theta)

            for p, phase in enumerate(self.phases):
                for s, lam in enumerate(lams):
                    fxx = np.cos(2*np.pi*x/lam + phase) * np.exp(-(x**2+y**2)/(2*sigma**2))
                    fxx -= np.mean(fxx)
                    fxx /= np.linalg.norm(fxx)

                    #if p==0:
                        #plt.subplot(len(oris),len(lams),count+1)
                        #plt.imshow(fxx,cmap=mpl.cm.gray,interpolation='bicubic')
                        #count+=1

                    gabors[o,p,s,:,:] = fxx
        plt.show()
        return gabors

def delta(phase,inc, ecc = 0, omega=0):
    """
    Compute the distance center-to-center between planet and host star.
    ___

    INPUT:

    phase: orbital phase in radian
    inc: inclination of the system in radian

    OPTIONAL INPUT:

    ecc:
    omega:

    //
    OUTPUT:

    distance center-to-center, double-float number.
    ___


    """
    phase = 2*np.pi*phase
    if ecc == 0 and omega == 0:
        delta = np.sqrt(1-(np.cos(phase)**2)*(np.sin(inc)**2))
    else:
        delta = (1.-ecc**2.)/(1.-ecc*np.sin(phase-omega))* np.sqrt((1.-(np.cos(phase))**2.*(np.sin(inc))**2))

    return delta

def rv_pqw(k, p, ecc, nu):
    """Returns r and v vectors in perifocal frame.

    """
    r_pqw = (np.array([cos(nu), sin(nu), 0 * nu]) * p / (1 + ecc * cos(nu))).T
    v_pqw = (np.array([-sin(nu), (ecc + cos(nu)), 0]) * sqrt(k / p)).T
    return r_pqw, v_pqw

def lomb(t, y, freq):
    r"""Calculates Lomb periodogram."""
    # Sets constants.
    nfreq = len(freq)
    fmax, fmin = freq[-1], freq[0]
    power = np.zeros(nfreq)
    f4pi = freq * 4 * np.pi
    pi2 = np.pi * 2.
    n = len(y)
    cosarg = np.zeros(n)
    sinarg = np.zeros(n)
    argu = np.zeros(n)
    var = np.cov(y)  # Variance.
    yn = y - y.mean()
    # Do one Lomb loop.
    for fi in range(nfreq):
        sinsum = np.sum(np.sin(f4pi[fi]) * t)
        cossum = np.sum(np.cos(f4pi[fi]) * t)

        tau = np.arctan2(sinsum, cossum)
        argu = pi2 * freq[fi] * (t - tau)

        cosarg = np.cos(argu)
        cfi = np.sum(yn * cosarg)
        cosnorm = np.sum(cosarg ** 2)

        sinarg = np.sin(argu)
        sfi = np.sum(yn * sinarg)
        sinnorm = np.sum(sinarg ** 2)

        power[fi] = (cfi ** 2 / cosnorm + sfi ** 2 / sinnorm) / 2 * var

    return power

def rotate (x, y, angle) :
    angle = angle*2*np.pi/360.
    rotMatrix = np.array([[np.cos(angle), -np.sin(angle)], 
                       [np.sin(angle),  np.cos(angle)]], dtype=np.float64)
    newx = x*rotMatrix[0,0] + y*rotMatrix[0,1]
    newy = x*rotMatrix[1,0] + y*rotMatrix[1,1]
    return newx,newy

def calc_Tb(thetak=np.pi/3., phik=np.pi/8., thetan=np.pi/3., phin=np.pi/4., 
            delta=0., Ts=11.1, Tg=57.23508, z=20, verbose=False,
            xalpha=34.247221, xc=0.004176, xB=0.365092, x1s=1.):
    
    """ Calculates brightness-temperature fluctuation T[K] from eq 1 of Paper II. 
    NOTE: Magnetic-field direction is along the z-axis!  It takes x's (all unitless), temperatures in [K], and angles in [rad]."""
    
    k_dot_n = np.cos(thetan)*np.cos(thetak) + np.sin(thetan)*np.sin(thetak)*np.cos(phin)*np.cos(phik) + np.sin(thetan)*np.sin(thetak)*np.sin(phin)*np.sin(phik)

    summ = 0.
    for i,m in enumerate( np.array([-2,-1,0,1,2]) ):
        summand = Y2( m,thetak,phik ) * np.conjugate( Y2( m,thetan,phin ) ) / (1. + xalpha + xc - 1j*m*xB)
        summ += summand.real

    first_term = 1 + delta + delta*k_dot_n**2
    second_term = 1 + 2.*delta + 2.*delta*k_dot_n**2 - delta*4.*np.pi/75.*summ
    
    res = x1s * ( 1 - Tg/Ts ) * np.sqrt( (1 + z)/10. ) * ( 26.4 * first_term - 0.128 * x1s * ( Tg/Ts ) * np.sqrt( (1 + z)/10. ) * second_term)
    
    if verbose:
        print '\n'
        print 'xalpha = %f' % xalpha
        print 'xc = %f' % xc
        print 'xB = %f' % xB
        print 'k_dot_n=%f' % k_dot_n
        print 'summ=%f' % summ
        print 'first=%f' % 26.4*first_term
        print 'second=%f' % second_term
        
    return res/1000. #this is to make it to K from mK.

    def array2raster(newRasterfn,rasterOrigin,pixelWidth,pixelHeight, data, variables, rotate=0):
        """Convert data dictionary (of arrays) into a multiband GeoTiff

        :param newRasterfn: filename to save to
        :param rasterOrigin: location of top left corner
        :param pixelWidth: e-w pixel size
        :param pixelHeight: n-s pixel size
        :param data: dictionary containing the data arrays
        :param variables: list of which keys from the dictionary to output
        :param rotate: Optional rotation angle (in radians)
        """
        cols = len(data['longitude'])
        rows = len(data['latitude'])
        originX = rasterOrigin[0]
        originY = rasterOrigin[1]

        we_res = np.cos(rotate) * pixelWidth
        rotX = np.sin(rotate) * pixelWidth
        rotY = -np.sin(rotate) * pixelHeight
        ns_res = np.cos(rotate) * pixelHeight

        driver = gdal.GetDriverByName('GTiff')
        nbands = len(variables)
        outRaster = driver.Create(newRasterfn, cols, rows, nbands, gdal.GDT_Float32)
        outRaster.SetGeoTransform((originX, we_res, rotX, originY, rotY, ns_res))
        for band,key in enumerate(variables, 1):
            outband = outRaster.GetRasterBand(band)
            outband.SetNoDataValue(0)
            outband.WriteArray(data[key])
            outband.FlushCache()
        outRasterSRS = osr.SpatialReference()
        outRasterSRS.ImportFromEPSG(4326)
        outRaster.SetProjection(outRasterSRS.ExportToWkt())

def spherical_to_cartesian(lons, lats, depths):
    """
    Return the position vectors (in Cartesian coordinates) of list of spherical
    coordinates.

    For equations see: http://mathworld.wolfram.com/SphericalCoordinates.html.

    Parameters are components of spherical coordinates in a form of scalars,
    lists or numpy arrays. ``depths`` can be ``None`` in which case it's
    considered zero for all points.

    :returns:
        ``numpy.array`` of 3d vectors representing points' coordinates in
        Cartesian space. The array has the same shape as parameter arrays.
        In particular it means that if ``lons`` and ``lats`` are scalars,
        the result is a single 3d vector. Vector of length ``1`` represents
        distance of 1 km.

    See also :func:`cartesian_to_spherical`.
    """
    phi = numpy.radians(lons)
    theta = numpy.radians(lats)
    if depths is None:
        rr = EARTH_RADIUS
    else:
        rr = EARTH_RADIUS - numpy.array(depths)
    cos_theta_r = rr * numpy.cos(theta)
    xx = cos_theta_r * numpy.cos(phi)
    yy = cos_theta_r * numpy.sin(phi)
    zz = rr * numpy.sin(theta)
    vectors = numpy.array([xx.transpose(), yy.transpose(), zz.transpose()]) \
                   .transpose()
    return vectors

def kdtree_fast(latvar,lonvar,lat0,lon0):
    '''
    :param latvar:
    :param lonvar:
    :param lat0:
    :param lon0:
    :return:
    '''
    rad_factor = pi/180.0 # for trignometry, need angles in radians
    # Read latitude and longitude from file into numpy arrays
    latvals = latvar[:] * rad_factor
    lonvals = lonvar[:] * rad_factor
    ny,nx = latvals.shape
    clat,clon = cos(latvals),cos(lonvals)
    slat,slon = sin(latvals),sin(lonvals)
    # Build kd-tree from big arrays of 3D coordinates
    triples = list(zip(ravel(clat*clon), ravel(clat*slon), ravel(slat)))
    kdt = cKDTree(triples)
    lat0_rad = lat0 * rad_factor
    lon0_rad = lon0 * rad_factor
    clat0,clon0 = cos(lat0_rad),cos(lon0_rad)
    slat0,slon0 = sin(lat0_rad),sin(lon0_rad)
    dist_sq_min, minindex_1d = kdt.query([clat0*clon0, clat0*slon0, slat0])
    iy_min, ix_min = unravel_index(minindex_1d, latvals.shape)
    return iy_min,ix_min

def tunnel_fast(latvar,lonvar,lat0,lon0):
    '''
    Find closest point in a set of (lat,lon) points to specified point
    latvar - 2D latitude variable from an open netCDF dataset
    lonvar - 2D longitude variable from an open netCDF dataset
    lat0, lon0 - query point
    Returns iy,ix such that the square of the tunnel distance
    between (latval[iy,ix], lonval[iy,ix]) and (lat0, lon0)
    is minimum.

    :param latvar:
    :param lonvar:
    :param lat0:
    :param lon0:

    :return:

    '''
    rad_factor = pi/180.0 # for trignometry, need angles in radians
    # Read latitude and longitude from file into numpy arrays
    latvals = latvar[:] * rad_factor
    lonvals = lonvar[:] * rad_factor
    ny,nx = latvals.shape
    lat0_rad = lat0 * rad_factor
    lon0_rad = lon0 * rad_factor
    # Compute numpy arrays for all values, no loops
    clat,clon = cos(latvals),cos(lonvals)
    slat,slon = sin(latvals),sin(lonvals)
    delX = cos(lat0_rad)*cos(lon0_rad) - clat*clon
    delY = cos(lat0_rad)*sin(lon0_rad) - clat*slon
    delZ = sin(lat0_rad) - slat;
    dist_sq = delX**2 + delY**2 + delZ**2
    minindex_1d = dist_sq.argmin()  # 1D index of minimum element
    iy_min,ix_min = unravel_index(minindex_1d, latvals.shape)
    return iy_min,ix_min

def random_rotation(x, rg, row_axis=1, col_axis=2, channel_axis=0,
                    fill_mode='nearest', cval=0.):
    """Performs a random rotation of a Numpy image tensor.

    # Arguments
        x: Input tensor. Must be 3D.
        rg: Rotation range, in degrees.
        row_axis: Index of axis for rows in the input tensor.
        col_axis: Index of axis for columns in the input tensor.
        channel_axis: Index of axis for channels in the input tensor.
        fill_mode: Points outside the boundaries of the input
            are filled according to the given mode
            (one of `{'constant', 'nearest', 'reflect', 'wrap'}`).
        cval: Value used for points outside the boundaries
            of the input if `mode='constant'`.

    # Returns
        Rotated Numpy image tensor.
    """
    theta = np.pi / 180 * np.random.uniform(-rg, rg)
    rotation_matrix = np.array([[np.cos(theta), -np.sin(theta), 0],
                                [np.sin(theta), np.cos(theta), 0],
                                [0, 0, 1]])

    h, w = x.shape[row_axis], x.shape[col_axis]
    transform_matrix = transform_matrix_offset_center(rotation_matrix, h, w)
    x = apply_transform(x, transform_matrix, channel_axis, fill_mode, cval)
    return x

def rotateAboutZaxis (x, y, z, alpha, verbose = 0) :
    if verbose : print "\t z axis rotation of ", alpha, "given  ", x[0], y[0], z[0]
    alpha = alpha*2*np.pi/360.
    xp = x*np.cos(alpha) - y*np.sin(alpha)
    yp = x*np.sin(alpha) + y*np.cos(alpha)
    zp = z
    return xp,yp,zp

def sphToCartesian(ra0, ra, dec, r=1) :
    ra = (ra-(ra0-90))*2*np.pi/360.
    dec = dec*2*np.pi/360.
    x = r * np.cos(ra)*np.cos(dec)
    y = r * np.sin(ra)*np.cos(dec)
    z = r * np.sin(dec)
    return x,y,z

def signaltest_xyt1(coastlines=False):
	"""
	Generate
	"""
	nt = 100
	nx = 128
	ny = 128
	t1d = np.linspace(0, 20 * np.pi, nt)
	x1d = np.linspace(0, 2 * np.pi, nx)
	y1d = np.linspace(0, 2 * np.pi, ny)
	t, y, x = np.meshgrid(t1d, y1d, x1d, indexing='ij')
	# Create four times modulation with
	m1 = np.cos(1.5 * t)
	m2 = np.cos(2 * t)
	m3 = np.cos(0.5 * t)
	m4 = np.cos(t)
	# Create a spatio-temporal gaussian noise
	noise = 0.8 * np.random.normal(0, 0.2, (nt, ny, nx))
	# Create four spatial patterns
	z1 = np.sin(x) * np.sin(y) * m1
	z2 = np.sin(2.5 * x) * np.sin(y) * m2
	z3 = np.sin(x) * np.sin(2.5 * y) * m3
	z4 = np.sin(2.5 * x) * np.sin(2.5 * y) * m4
	z = z1 + z2 + z3 + z4 + noise
	if coastlines:
		z[:, 0:ny/4, 0:nx/4] = np.nan
	return xr.DataArray(z, coords=[t1d, y1d, x1d], dims=['time', 'y', 'x'], name='signal')

def rotateAboutXaxis (x, y, z, alpha, verbose = 0) :
    if verbose : print "\t x axis rotation of ", alpha, "given  ", x[0], y[0], z[0]
    alpha = alpha*2*np.pi/360.
    xp = x
    yp = y*np.cos(alpha) - z*np.sin(alpha)
    zp = y*np.sin(alpha) + z*np.cos(alpha)
    return xp,yp,zp

        def get_new_cell(self):
            """Returns new basis vectors"""
            a = np.sqrt(self.a)
            b = np.sqrt(self.b)
            c = np.sqrt(self.c)

            ad = self.atoms.cell[0] / np.linalg.norm(self.atoms.cell[0])

            Z = np.cross(self.atoms.cell[0], self.atoms.cell[1])
            Z /= np.linalg.norm(Z)
            X = ad - np.dot(ad, Z) * Z
            X /= np.linalg.norm(X)
            Y = np.cross(Z, X)

            alpha = np.arccos(self.x / (2 * b * c))
            beta = np.arccos(self.y / (2 * a * c))
            gamma = np.arccos(self.z / (2 * a * b))

            va = a * np.array([1, 0, 0])
            vb = b * np.array([np.cos(gamma), np.sin(gamma), 0])
            cx = np.cos(beta)
            cy = (np.cos(alpha) - np.cos(beta) * np.cos(gamma)) \
                / np.sin(gamma)
            cz = np.sqrt(1. - cx * cx - cy * cy)
            vc = c * np.array([cx, cy, cz])

            abc = np.vstack((va, vb, vc))
            T = np.vstack((X, Y, Z))
            return np.dot(abc, T)

def compute_slat_lift(slat_angle,sweep_angle):
    """ SUAVE.Methods.Aerodynamics.compute_slat_lift(vehicle):
        Computes the increase of lift due to slat deployment

        Inputs:
            slat_angle  - Slat deflection angle - [rad]
            sweep_angle - Wing sweep angle      - [rad]

        Outputs:
            dcl_slat    - Lift coefficient increase due to slat

        Assumptions:
            if needed

    """

    #unpack
    sa = slat_angle  / Units.deg
    sw = sweep_angle

    #---AA241 Method
    dcl_slat = (sa/23.)*(np.cos(sw))**1.4 * np.cos(sa * Units.deg)**2

    #returning dcl_slat
    return dcl_slat