Python numpy.radians() Examples

def HaversineDistance(lon1,lat1,lon2,lat2):
    """
    Function to calculate the great circle distance between two points
    using the Haversine formula
    """
    R = 6371. #Mean radius of the Earth

    # convert decimal degrees to radians
    lon1, lat1, lon2, lat2 = map(np.radians, [lon1, lat1, lon2, lat2])

    # haversine formula
    dlon = lon2 - lon1
    dlat = lat2 - lat1
    a = np.sin(dlat/2.)**2. + np.cos(lat1) * np.cos(lat2) * np.sin(dlon/2.)**2.
    c = 2.*np.arcsin(np.sqrt(a))
    distance = R * c

    return distance 

def calculate_dayside_reconnection(inst):
    """ Calculate the dayside reconnection rate (Milan et al. 2014)

    Parameters
    -----------
    inst : pysat.Instrument
        Instrument with OMNI HRO data, requires BYZ_GSM and clock_angle

    Notes
    --------
    recon_day = 3.8 Re (Vx / 4e5 m/s)^1/3 Vx B_yz (sin(theta/2))^9/2
    """
    rearth = 6371008.8
    sin_htheta = np.power(np.sin(np.radians(0.5 * inst['clock_angle'])), 4.5)
    byz = inst['BYZ_GSM'] * 1.0e-9
    vx = inst['flow_speed'] * 1000.0

    recon_day = 3.8 * rearth * vx * byz * sin_htheta * np.power((vx / 4.0e5),
                                                                1.0/3.0)
    inst['recon_day'] = pds.Series(recon_day, index=inst.data.index)
    return 

def coords_to_vec(lon, lat):
    """ Converts longitute and latitude coordinates to a unit 3-vector

    return array(3,n) with v_x[i],v_y[i],v_z[i] = directional cosines
    """
    phi = np.radians(lon)
    theta = (np.pi / 2) - np.radians(lat)
    sin_t = np.sin(theta)
    cos_t = np.cos(theta)

    xVals = sin_t * np.cos(phi)
    yVals = sin_t * np.sin(phi)
    zVals = cos_t

    # Stack them into the output array
    out = np.vstack((xVals, yVals, zVals)).swapaxes(0, 1)
    return out 

def rotipp(acceleration_x, time_step_x, acceleration_y, time_step_y, periods,
        percentile, damping=0.05, units="cm/s/s", method="Nigam-Jennings"):
    """
    Returns the rotationally independent spectrum RotIpp as defined by
    Boore (2010)
    """
    if np.fabs(time_step_x - time_step_y) > 1E-10:
        raise ValueError("Record pair must have the same time-step!")
    acceleration_x, acceleration_y = equalise_series(acceleration_x,
                                                     acceleration_y)
    target, rota, rotv, rotd, angles = rotdpp(acceleration_x, time_step_x,
                                              acceleration_y, time_step_y,
                                              periods, percentile, damping,
                                              units, method)
    locn, penalty = _get_gmrotd_penalty(
        np.hstack([target["PGA"],target["Pseudo-Acceleration"]]),
        rota)
    target_theta = np.radians(angles[locn])
    arotpp = acceleration_x * np.cos(target_theta) +\
        acceleration_y * np.sin(target_theta)
    spec = get_response_spectrum(arotpp, time_step_x, periods, damping, units,
        method)[0]
    spec["GMRot{:2.0f}".format(percentile)] = target
    return spec 

def _latlon_to_cart(lat, lon, R = typhon.constants.earth_radius):
    """
    Simple conversion of latitude and longitude to Cartesian coordinates.
    Approximates the Earth as sphere with radius :code:`R` and computes
    cartesian x, y, z coordinates with the center of the Earth as origin.

    Args:
        lat: Array of latitude coordinates.
        lon: Array of longitude coordinates.
        R: The radius to assume.
    Returns:
        Tuple :code:`(x, y, z)` of arrays :code:`x, y, z` containing the
        resulting x-, y- and z-coordinates.
    """
    lat = np.radians(lat)
    lon = np.radians(lon)
    x = R * np.cos(lat) * np.cos(lon)
    y = R * np.cos(lat) * np.sin(lon)
    z = R * np.sin(lat)
    return x, y, z 

def orthogonalization_matrix(lengths, angles):
    """Return orthogonalization matrix for crystallographic cell coordinates.

    Angles are expected in degrees.

    The de-orthogonalization matrix is the inverse.

    >>> O = orthogonalization_matrix([10, 10, 10], [90, 90, 90])
    >>> numpy.allclose(O[:3, :3], numpy.identity(3, float) * 10)
    True
    >>> O = orthogonalization_matrix([9.8, 12.0, 15.5], [87.2, 80.7, 69.7])
    >>> numpy.allclose(numpy.sum(O), 43.063229)
    True

    """
    a, b, c = lengths
    angles = numpy.radians(angles)
    sina, sinb, _ = numpy.sin(angles)
    cosa, cosb, cosg = numpy.cos(angles)
    co = (cosa * cosb - cosg) / (sina * sinb)
    return numpy.array([
        [ a*sinb*math.sqrt(1.0-co*co),  0.0,    0.0, 0.0],
        [-a*sinb*co,                    b*sina, 0.0, 0.0],
        [ a*cosb,                       b*cosa, c,   0.0],
        [ 0.0,                          0.0,    0.0, 1.0]]) 

def get_displacement_km(n1, x1, y1, n2, x2, y2):
    ''' Find displacement in kilometers using Haversine
        http://www.movable-type.co.uk/scripts/latlong.html
    Parameters
    ----------
        n1 : First Nansat object
        x1 : 1D vector - X coordinates of keypoints on image 1
        y1 : 1D vector - Y coordinates of keypoints on image 1
        n2 : Second Nansat object
        x1 : 1D vector - X coordinates of keypoints on image 2
        y1 : 1D vector - Y coordinates of keypoints on image 2
    Returns
    -------
        h : 1D vector - total displacement, km
    '''
    lon1, lat1 = n1.transform_points(x1, y1)
    lon2, lat2 = n2.transform_points(x2, y2)

    lt1, ln1, lt2, ln2 = map(np.radians, (lat1, lon1, lat2, lon2))
    dlat = lt2 - lt1
    dlon = ln2 - ln1
    d = (np.sin(dlat * 0.5) ** 2 +
         np.cos(lt1) * np.cos(lt2) * np.sin(dlon * 0.5) ** 2)
    return 2 * AVG_EARTH_RADIUS * np.arcsin(np.sqrt(d)) 

def __init__(self, rotation: Vector, translation: Vector, angle_unit: str, notation: str='XYZ'):
        self.rotation    = rotation
        self.translation = translation

        self.angle_unit = angle_unit
        if self.angle_unit == 'degrees':
            self.rotation = Vector(*[radians(alpha) for alpha in rotation])

        self.R_x = None
        self.R_y = None
        self.R_z = None
        self.T   = None

        self.matrix = None
        self.notation = notation

        self._update_matrix() 

def orthogonalization_matrix(lengths, angles):
    """Return orthogonalization matrix for crystallographic cell coordinates.

    Angles are expected in degrees.

    The de-orthogonalization matrix is the inverse.

    >>> O = orthogonalization_matrix([10, 10, 10], [90, 90, 90])
    >>> numpy.allclose(O[:3, :3], numpy.identity(3, float) * 10)
    True
    >>> O = orthogonalization_matrix([9.8, 12.0, 15.5], [87.2, 80.7, 69.7])
    >>> numpy.allclose(numpy.sum(O), 43.063229)
    True

    """
    a, b, c = lengths
    angles = numpy.radians(angles)
    sina, sinb, _ = numpy.sin(angles)
    cosa, cosb, cosg = numpy.cos(angles)
    co = (cosa * cosb - cosg) / (sina * sinb)
    return numpy.array([
        [ a*sinb*math.sqrt(1.0-co*co),  0.0,    0.0, 0.0],
        [-a*sinb*co,                    b*sina, 0.0, 0.0],
        [ a*cosb,                       b*cosa, c,   0.0],
        [ 0.0,                          0.0,    0.0, 1.0]]) 

def test_true_to_mean(self):
        # Data from Schlesinger & Udick, 1912
        data = [
            # ecc, M (deg), ta (deg)
            (0.0, 0.0, 0.0),
            (0.05, 10.0, 11.06),
            (0.06, 30.0, 33.67),
            (0.04, 120.0, 123.87),
            (0.14, 65.0, 80.50),
            (0.19, 21.0, 30.94),
            (0.35, 65.0, 105.71),
            (0.48, 180.0, 180.0),
            (0.75, 125.0, 167.57)
        ]
        for row in data:
            ecc, expected_M, ta = row

            M = angles.ta_to_M(radians(ta), ecc)

            self.assertAlmostEqual(degrees(M), expected_M, places=1) 

def test_mean_to_true(self):
        # Data from Schlesinger & Udick, 1912
        data = [
            # ecc, M (deg), ta (deg)
            (0.0, 0.0, 0.0),
            (0.05, 10.0, 11.06),
            (0.06, 30.0, 33.67),
            (0.04, 120.0, 123.87),
            (0.14, 65.0, 80.50),
            (0.19, 21.0, 30.94),
            (0.35, 65.0, 105.71),
            (0.48, 180.0, 180.0),
            (0.75, 125.0, 167.57)
        ]
        for row in data:
            ecc, M, expected_ta = row

            ta = angles.M_to_ta(radians(M), ecc)

            self.assertAlmostEqual(degrees(ta), expected_ta, places=2) 

def orthogonalization_matrix(lengths, angles):
    """Return orthogonalization matrix for crystallographic cell coordinates.
    Angles are expected in degrees.
    The de-orthogonalization matrix is the inverse.
    >>> O = orthogonalization_matrix((10., 10., 10.), (90., 90., 90.))
    >>> numpy.allclose(O[:3, :3], numpy.identity(3, float) * 10)
    True
    >>> O = orthogonalization_matrix([9.8, 12.0, 15.5], [87.2, 80.7, 69.7])
    >>> numpy.allclose(numpy.sum(O), 43.063229)
    True
    """
    a, b, c = lengths
    angles = numpy.radians(angles)
    sina, sinb, _ = numpy.sin(angles)
    cosa, cosb, cosg = numpy.cos(angles)
    co = (cosa * cosb - cosg) / (sina * sinb)
    return numpy.array((
        ( a*sinb*math.sqrt(1.0-co*co),  0.0,    0.0, 0.0),
        (-a*sinb*co,                    b*sina, 0.0, 0.0),
        ( a*cosb,                       b*cosa, c,   0.0),
        ( 0.0,                          0.0,    0.0, 1.0)),
        dtype=numpy.float64) 

def test_true_to_eccentric(self):
        # Data from NASA-TR-R-158
        data = [
            # ecc, E (deg), ta(deg)
            (0.0, 0.0, 0.0),
            (0.05, 10.52321, 11.05994),
            (0.10, 54.67466, 59.49810),
            (0.35, 142.27123, 153.32411),
            (0.61, 161.87359, 171.02189)
        ]
        for row in data:
            ecc, expected_E, ta = row

            E = angles.ta_to_E(radians(ta), ecc)

            self.assertAlmostEqual(degrees(E), expected_E, places=4) 

def orthogonalization_matrix(lengths, angles):
    """Return orthogonalization matrix for crystallographic cell coordinates.

    Angles are expected in degrees.

    The de-orthogonalization matrix is the inverse.

    >>> O = orthogonalization_matrix([10, 10, 10], [90, 90, 90])
    >>> numpy.allclose(O[:3, :3], numpy.identity(3, float) * 10)
    True
    >>> O = orthogonalization_matrix([9.8, 12.0, 15.5], [87.2, 80.7, 69.7])
    >>> numpy.allclose(numpy.sum(O), 43.063229)
    True

    """
    a, b, c = lengths
    angles = numpy.radians(angles)
    sina, sinb, _ = numpy.sin(angles)
    cosa, cosb, cosg = numpy.cos(angles)
    co = (cosa * cosb - cosg) / (sina * sinb)
    return numpy.array([
        [ a*sinb*math.sqrt(1.0-co*co),  0.0,    0.0, 0.0],
        [-a*sinb*co,                    b*sina, 0.0, 0.0],
        [ a*cosb,                       b*cosa, c,   0.0],
        [ 0.0,                          0.0,    0.0, 1.0]]) 

def orthogonalization_matrix(lengths, angles):
    """Return orthogonalization matrix for crystallographic cell coordinates.

    Angles are expected in degrees.

    The de-orthogonalization matrix is the inverse.

    >>> O = orthogonalization_matrix([10, 10, 10], [90, 90, 90])
    >>> numpy.allclose(O[:3, :3], numpy.identity(3, float) * 10)
    True
    >>> O = orthogonalization_matrix([9.8, 12.0, 15.5], [87.2, 80.7, 69.7])
    >>> numpy.allclose(numpy.sum(O), 43.063229)
    True

    """
    a, b, c = lengths
    angles = numpy.radians(angles)
    sina, sinb, _ = numpy.sin(angles)
    cosa, cosb, cosg = numpy.cos(angles)
    co = (cosa * cosb - cosg) / (sina * sinb)
    return numpy.array([
        [ a*sinb*math.sqrt(1.0-co*co),  0.0,    0.0, 0.0],
        [-a*sinb*co,                    b*sina, 0.0, 0.0],
        [ a*cosb,                       b*cosa, c,   0.0],
        [ 0.0,                          0.0,    0.0, 1.0]]) 

def orthogonalization_matrix(lengths, angles):
    """Return orthogonalization matrix for crystallographic cell coordinates.

    Angles are expected in degrees.

    The de-orthogonalization matrix is the inverse.

    >>> O = orthogonalization_matrix([10, 10, 10], [90, 90, 90])
    >>> numpy.allclose(O[:3, :3], numpy.identity(3, float) * 10)
    True
    >>> O = orthogonalization_matrix([9.8, 12.0, 15.5], [87.2, 80.7, 69.7])
    >>> numpy.allclose(numpy.sum(O), 43.063229)
    True

    """
    a, b, c = lengths
    angles = numpy.radians(angles)
    sina, sinb, _ = numpy.sin(angles)
    cosa, cosb, cosg = numpy.cos(angles)
    co = (cosa * cosb - cosg) / (sina * sinb)
    return numpy.array([
        [ a*sinb*math.sqrt(1.0-co*co),  0.0,    0.0, 0.0],
        [-a*sinb*co,                    b*sina, 0.0, 0.0],
        [ a*cosb,                       b*cosa, c,   0.0],
        [ 0.0,                          0.0,    0.0, 1.0]]) 

def orthogonalization_matrix(lengths, angles):
    """Return orthogonalization matrix for crystallographic cell coordinates.

    Angles are expected in degrees.

    The de-orthogonalization matrix is the inverse.

    >>> O = orthogonalization_matrix([10, 10, 10], [90, 90, 90])
    >>> numpy.allclose(O[:3, :3], numpy.identity(3, float) * 10)
    True
    >>> O = orthogonalization_matrix([9.8, 12.0, 15.5], [87.2, 80.7, 69.7])
    >>> numpy.allclose(numpy.sum(O), 43.063229)
    True

    """
    a, b, c = lengths
    angles = numpy.radians(angles)
    sina, sinb, _ = numpy.sin(angles)
    cosa, cosb, cosg = numpy.cos(angles)
    co = (cosa * cosb - cosg) / (sina * sinb)
    return numpy.array([
        [ a*sinb*math.sqrt(1.0-co*co),  0.0,    0.0, 0.0],
        [-a*sinb*co,                    b*sina, 0.0, 0.0],
        [ a*cosb,                       b*cosa, c,   0.0],
        [ 0.0,                          0.0,    0.0, 1.0]]) 

def orthogonalization_matrix(lengths, angles):
  """Return orthogonalization matrix for crystallographic cell coordinates.

  Angles are expected in degrees.

  The de-orthogonalization matrix is the inverse.

  >>> O = orthogonalization_matrix([10, 10, 10], [90, 90, 90])
  >>> numpy.allclose(O[:3, :3], numpy.identity(3, float) * 10)
  True
  >>> O = orthogonalization_matrix([9.8, 12.0, 15.5], [87.2, 80.7, 69.7])
  >>> numpy.allclose(numpy.sum(O), 43.063229)
  True

  """
  a, b, c = lengths
  angles = numpy.radians(angles)
  sina, sinb, _ = numpy.sin(angles)
  cosa, cosb, cosg = numpy.cos(angles)
  co = (cosa * cosb - cosg) / (sina * sinb)
  return numpy.array([
    [a * sinb * math.sqrt(1.0 - co * co), 0.0, 0.0, 0.0],
    [-a * sinb * co, b * sina, 0.0, 0.0],
    [a * cosb, b * cosa, c, 0.0],
    [0.0, 0.0, 0.0, 1.0]]) 

def test_rotate_simple(self):
        vec = np.array([1, 0, 0])

        assert_allclose(rotate(vec, 0, np.radians(90)), np.array([1, 0, 0]), atol=1e-16)
        assert_allclose(rotate(vec, 1, np.radians(90)), np.array([0, 0, -1]), atol=1e-16)
        assert_allclose(rotate(vec, 2, np.radians(90)), np.array([0, 1, 0]), atol=1e-16) 

def _view_rectspherebivariate(cls, data, data_view, target_view, coords_in_radians=False,
                                  kx=3, ky=3, s=0, x_tolerance=1e-3, y_tolerance=1e-3):
        t_xmin, t_ymin, t_xmax, t_ymax = target_view.bbox
        d_xmin, d_ymin, d_xmax, d_ymax = data_view.bbox
        nodata = target_view.nodata
        yx_tolerance = np.sqrt(x_tolerance**2 + y_tolerance**2)
        target_dx, target_dy = target_view.affine.a, target_view.affine.e
        data_dx, data_dy = data_view.affine.a, data_view.affine.e
        viewrows, viewcols = target_view.grid_indices(col_ascending=True,
                                                      row_ascending=True)
        rows, cols = data_view.grid_indices(col_ascending=True,
                                            row_ascending=True)
        viewrows += target_dy
        viewcols += target_dx
        rows += data_dy
        cols += data_dx
        row_bool = (rows <= t_ymax + y_tolerance) & (rows >= t_ymin - y_tolerance)
        col_bool = (cols <= t_xmax + x_tolerance) & (cols >= t_xmin - x_tolerance)
        if not coords_in_radians:
            rows = np.radians(rows) + np.pi/2
            cols = np.radians(cols) + np.pi
            viewrows = np.radians(viewrows) + np.pi/2
            viewcols = np.radians(viewcols) + np.pi
        rsbs_interpolator = (interpolate.
                            RectBivariateSpline(rows[row_bool],
                                                cols[col_bool],
                                                data[np.ix_(row_bool[::-1], col_bool)],
                                                kx=kx, ky=ky, s=s))
        xy_query = np.vstack(np.dstack(np.meshgrid(viewrows, viewcols, indexing='ij')))
        view = rsbs_interpolator.ev(xy_query[:,0], xy_query[:,1]).reshape(target_view.shape)
        return view 

def test_interp_healpix_nside(cst_efield_2freq_cut, cst_efield_2freq_cut_healpix):
    efield_beam = cst_efield_2freq_cut

    efield_beam.interpolation_function = "az_za_simple"

    # test calling interp with healpix parameters directly gives same result
    min_res = np.min(
        np.array(
            [np.diff(efield_beam.axis1_array)[0], np.diff(efield_beam.axis2_array)[0]]
        )
    )
    nside_min_res = np.sqrt(3 / np.pi) * np.radians(60.0) / min_res
    nside = int(2 ** np.ceil(np.log2(nside_min_res)))

    new_efield_beam = cst_efield_2freq_cut_healpix
    assert new_efield_beam.nside == nside

    new_efield_beam.interpolation_function = "healpix_simple"

    # check error with cut sky
    with pytest.raises(
        ValueError, match="simple healpix interpolation requires full sky healpix maps."
    ):
        new_efield_beam.interp(
            az_array=efield_beam.axis1_array,
            za_array=efield_beam.axis2_array,
            az_za_grid=True,
            new_object=True,
        ) 

def distance(self, x1, y1, x2, y2):
        if PAR.UNITS in ['lonlat']:
            dlat = np.radians(y2-y1)
            dlon = np.radians(x2-x1)
            a = np.sin(dlat/2) * np.sin(dlat/2) + np.cos(np.radians(y1)) \
                * np.cos(np.radians(y2)) * np.sin(dlon/2) * np.sin(dlon/2)
            D = 2 * np.arctan2(np.sqrt(a), np.sqrt(1-a))
            D *= 180/np.pi
            return D

        else:
            return ((x1-x2)**2 + (y1-y2)**2)**0.5 

def haversine(latlon1, latlon2):
    r"""
    Calculate the great circle distance between two points
    on the earth (specified in decimal degrees)

    Args:
        latlon1 (tuple): (lat, lon)
        latlon2 (tuple): (lat, lon)

    Returns:
        float : distance in kilometers

    References:
        en.wikipedia.org/wiki/Haversine_formula
        gis.stackexchange.com/questions/81551/matching-gps-tracks
        stackoverflow.com/questions/4913349/haversine-distance-gps-points

    Doctest:
        >>> from ibeis.algo.preproc.occurrence_blackbox import *  # NOQA
        >>> import scipy.spatial.distance as spdist
        >>> import functools
        >>> latlon1 = [-80.21895315, -158.81099213]
        >>> latlon2 = [  9.77816711,  -17.27471498]
        >>> kilometers = haversine(latlon1, latlon2)
        >>> result = ('kilometers = %s' % (kilometers,))
        >>> print(result)
        kilometers = 11930.909364189827
    """
    # convert decimal degrees to radians
    lat1, lon1 = np.radians(latlon1)
    lat2, lon2 = np.radians(latlon2)
    # haversine formula
    dlon = lon2 - lon1
    dlat = lat2 - lat1
    a = (np.sin(dlat / 2) ** 2) + np.cos(lat1) * np.cos(lat2) * (np.sin(dlon / 2) ** 2)
    c = 2 * np.arcsin(np.sqrt(a))
    # convert to kilometers
    EARTH_RADIUS_KM = 6367
    kilometers = EARTH_RADIUS_KM * c
    return kilometers 

def distance_between2pts(self, coord0, coord1, in_xy=False):
        """
        Calculate distance between two sets of coordinates.

        Parameters
        ----------
        coord0 : tuple of floats
            Coordinate pair in degrees
        coord1 : tuple of floats
            Coordinate pair in degrees
        in_xy : Boolean
            Default False.  If True, the pair is in (lon, lat)
            rather than (lat, lon)

        Returns
        -------
        distance : float
            Great circle distance in meters.

        """
        coord0 = np.radians(coord0)
        coord1 = np.radians(coord1)

        coord_order = {False: [0, 1],
                       True: [1, 0]}

        a, b = coord_order[in_xy]

        distance = (np.arccos(np.sin(coord1[a]) * np.sin(coord0[a]) +
                              np.cos(coord1[a]) * np.cos(coord0[a]) *
                              np.cos(coord0[b] - coord1[b])) * 6371) * 1000

        return distance 

def cmtest(*args, **kwargs):
    """
    Non parametric multi-sample test for equal medians. Similar to a
    Kruskal-Wallis test for linear data.

    H0: the s populations have equal medians
    HA: the s populations have unequal medians

    :param alpha1: angles in radians
    :param alpha2: angles in radians
    :returns: p-value and test statistic of the common median test


    References: [Fisher1995]_

    """
    axis = kwargs.get('axis', None)
    if axis is None:
        alpha = list(map(np.ravel, args))
    else:
        alpha = args

    s = len(alpha)
    n = [(0 * a + 1).sum(axis=axis) for a in alpha]
    N = sum(n)

    med = descriptive.median(np.concatenate(alpha, axis=axis), axis=axis)
    if axis is not None:
        med = np.expand_dims(med, axis=axis)

    m = [np.sum(descriptive.cdiff(a, med) < 0, axis=axis) for a in alpha]
    if np.any([nn < 10 for nn in n]):
        warnings.warn('Test not applicable. Sample size in at least one group to small.')
    M = sum(m)
    P = (N ** 2. / (M * (N - M))) * sum([mm ** 2. / nn for mm, nn in zip(m, n)]) - N * M / (N - M)
    pval = stats.chi2.sf(P, df=s - 1)
    return pval, P 

def _sample_cdf(alpha, resolution=100., axis=None):
    """

    Helper function for circ_kuipertest.
    Evaluates CDF of sample in thetas.

    :param alpha: sample (in radians)
    :param resolution: resolution at which the cdf is evaluated (default 100)
    :param axis: axis along which the cdf is computed
    :returns: points at which cdf is evaluated, cdf values

    """

    if axis is None:
        alpha = alpha.ravel()
        axis = 0
    bins = np.linspace(0, 2 * np.pi, resolution + 1)
    old_shape = alpha.shape
    alpha = alpha % (2 * np.pi)

    alpha = alpha.reshape((alpha.shape[0], int(np.prod(alpha.shape[1:])))).T
    cdf = np.array([np.histogram(a, bins=bins)[0]
                    for a in alpha]).cumsum(axis=1) / float(alpha.shape[1])
    cdf = cdf.T.reshape((len(bins) - 1,) + old_shape[1:])

    return bins[:-1], cdf 

def heading(p1, p2):
    p1 = np.radians(p1)
    p2 = np.radians(p2)
    lat1, lon1 = p1[1], p1[0]
    lat2, lon2 = p2[1], p2[0]
    aa = np.sin(lon2 - lon1) * np.cos(lat2)
    bb = np.cos(lat1) * np.sin(lat2) - np.sin(lat1) * np.cos(lat2) * np.cos(lon2 - lon1)
    return np.arctan2(aa, bb) + np.pi 

def rotate_coordinates(center: Tuple[float, float],
                       position: Union[Tuple[float, float], np.ndarray],
                       angle: float) -> Tuple[float, float]:
    """
    Function to rotate coordinates around the image center.

    Parameters
    ----------
    center : tuple(float, float)
        Image center (y, x).
    position : tuple(float, float)
        Position (y, x) in the image, or a 2D numpy array of positions.
    angle : float
        Angle (deg) to rotate in counterclockwise direction.

    Returns
    -------
    tuple(float, float)
        New position (y, x).
    """

    pos_y = (position[1] - center[1]) * math.sin(np.radians(angle)) + \
            (position[0] - center[0]) * math.cos(np.radians(angle))

    pos_x = (position[1] - center[1]) * math.cos(np.radians(angle)) - \
            (position[0] - center[0]) * math.sin(np.radians(angle))

    return center[0] + pos_y, center[1] + pos_x 

def polar_to_cartesian(image: np.ndarray,
                       sep: float,
                       ang: float) -> Tuple[float, float]:
    """
    Function to convert polar coordinates to pixel coordinates.

    Parameters
    ----------
    image : numpy.ndarray
        Input image (2D or 3D).
    sep : float
        Separation (pix).
    ang : float
        Position angle (deg), measured counterclockwise with respect to the positive y-axis.

    Returns
    -------
    tuple(float, float)
        Cartesian coordinates (y, x). The bottom left corner of the image is (-0.5, -0.5).
    """

    center = center_subpixel(image)  # (y, x)

    x_pos = center[1] + sep * math.cos(math.radians(ang + 90))
    y_pos = center[0] + sep * math.sin(math.radians(ang + 90))

    return y_pos, x_pos 

def rotate_into(self, direction, rotation):
        """ slerp the given rotation values into the direction given
        by the brush_state
        @param direction MVector: the target direction
        @param rotation MVector: current euler rotation """

        vector_weight = self.brush_state.settings['strength']
        up_vector = om.MVector(0, 1, 0)
        local_up = up_vector.rotateBy(om.MEulerRotation(math.radians(rotation.x),
                                                        math.radians(rotation.y),
                                                        math.radians(rotation.z)))

        target_rotation = om.MQuaternion(local_up, direction, vector_weight)

        util = om.MScriptUtil()
        x_rot = np.radians(rotation.x)
        y_rot = np.radians(rotation.y)
        z_rot = np.radians(rotation.z)
        util.createFromDouble(x_rot, y_rot, z_rot)
        rotation_ptr = util.asDoublePtr()
        mat = om.MTransformationMatrix()
        mat.setRotation(rotation_ptr, om.MTransformationMatrix.kXYZ)

        mat = mat.asMatrix() * target_rotation.asMatrix()
        rotation = om.MTransformationMatrix(mat).rotation()

        return om.MVector(math.degrees(rotation.asEulerRotation().x),
                        math.degrees(rotation.asEulerRotation().y),
                        math.degrees(rotation.asEulerRotation().z))