docs: 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