def fit_angular_distribution(angles, rates, rate_errors, shape='pexp'):
    """Fits angular distribution of rates.

    Parameters
    ----------
    rates: numpy array
      with rates for all PMT combinations
    angles: numpy array
      with angles for all PMT combinations
    shape:
      which function to fit; exp for exponential or pexp for
      exponential_polinomial

    Returns
    -------
    fitted_rates: numpy array of fitted rates (fit_function(angles, popt...))

    """
    if shape == 'exp':
        fit_function = exponential
        # p0 = [-0.91871169,  2.72224241, -1.19065965,  1.48054122]
    if shape == 'pexp':
        fit_function = exponential_polinomial
        # p0 = [0.34921202, 2.8629577]

    cos_angles = np.cos(angles)
    popt, pcov = optimize.curve_fit(fit_function, cos_angles, rates)
    fitted_rates = fit_function(cos_angles, *popt)
    return fitted_rates, popt, pcov