docs: def fit_circular_gaussian(samples, high=np.pi, low=0):
      """Compute the circular mean for samples in a range
  
      Args:
          samples (ndarray): a one or two dimensional array. If one dimensional we calculate the fit using all
              values. If two dimensional, we fit the Gaussian for every set of samples over the first dimension.
          high (float): The maximum wrap point
          low (float): The minimum wrap point
      """
      cl_func = SimpleCLFunction.from_string('''
          void compute(global mot_float_type* samples,
                       global mot_float_type* means,
                       global mot_float_type* stds,
                       int nmr_samples,
                       int low,
                       int high){ 
          
              double cos_mean = 0;
              double sin_mean = 0;
              double ang;
  
              for(uint i = 0; i < nmr_samples; i++){
                  ang = (samples[i] - low)*2*M_PI / (high - low);
  
                  cos_mean += (cos(ang) - cos_mean) / (i + 1);
                  sin_mean += (sin(ang) - sin_mean) / (i + 1);
              }
  
              double R = hypot(cos_mean, sin_mean);
              if(R > 1){
                  R = 1;
              }
  
              double stds = 1/2. * sqrt(-2 * log(R));
  
              double res = atan2(sin_mean, cos_mean);
              if(res < 0){
                   res += 2 * M_PI;
              }
  
              *(means) = res*(high - low)/2.0/M_PI + low;
              *(stds) = ((high - low)/2.0/M_PI) * sqrt(-2*log(R));
          }
      ''')
  
      def run_cl(samples):
          data = {'samples': Array(samples, 'mot_float_type'),
                  'means': Zeros(samples.shape[0], 'mot_float_type'),
                  'stds': Zeros(samples.shape[0], 'mot_float_type'),
                  'nmr_samples': Scalar(samples.shape[1]),
                  'low': Scalar(low),
                  'high': Scalar(high),
                  }
  
          cl_func.evaluate(data, samples.shape[0])
          return data['means'].get_data(), data['stds'].get_data()
  
      if len(samples.shape) == 1:
          mean, std = run_cl(samples[None, :])
          return mean[0], std[0]
      return run_cl(samples)