combine_sources(int n_sel_sources, struct timespec *ref_time, double *offset,
                double *offset_sd, double *frequency, double *frequency_sd, double *skew)
{
  struct timespec src_ref_time;
  double src_offset, src_offset_sd, src_frequency, src_frequency_sd, src_skew;
  double src_root_delay, src_root_dispersion, sel_src_distance, elapsed;
  double offset_weight, sum_offset_weight, sum_offset, sum2_offset_sd;
  double frequency_weight, sum_frequency_weight, sum_frequency;
  double inv_sum2_frequency_sd, inv_sum2_skew;
  int i, index, combined;

  if (n_sel_sources == 1)
    return 1;

  sum_offset_weight = sum_offset = sum2_offset_sd = 0.0;
  sum_frequency_weight = sum_frequency = inv_sum2_frequency_sd = inv_sum2_skew = 0.0;

  sel_src_distance = sources[selected_source_index]->sel_info.root_distance;
  if (sources[selected_source_index]->type == SRC_NTP)
    sel_src_distance += reselect_distance;

  for (i = combined = 0; i < n_sel_sources; i++) {
    index = sel_sources[i];
    SST_GetTrackingData(sources[index]->stats, &src_ref_time,
                        &src_offset, &src_offset_sd,
                        &src_frequency, &src_frequency_sd, &src_skew,
                        &src_root_delay, &src_root_dispersion);

    /* Don't include this source if its distance is longer than the distance of
       the selected source multiplied by the limit, their estimated frequencies
       are not close, or it was recently marked as distant */

    if (index != selected_source_index &&
        (sources[index]->sel_info.root_distance > combine_limit * sel_src_distance ||
         fabs(*frequency - src_frequency) >
           combine_limit * (*skew + src_skew + LCL_GetMaxClockError()))) {
      /* Use a smaller penalty in first few updates */
      sources[index]->distant = sources[index]->reachability_size >= SOURCE_REACH_BITS ?
                                DISTANT_PENALTY : 1;
    } else if (sources[index]->distant) {
      sources[index]->distant--;
    }

    if (sources[index]->distant) {
      sources[index]->status = SRC_DISTANT;
      continue;
    }

    if (sources[index]->status == SRC_OK)
      sources[index]->status = SRC_UNSELECTED;

    elapsed = UTI_DiffTimespecsToDouble(ref_time, &src_ref_time);
    src_offset += elapsed * src_frequency;
    src_offset_sd += elapsed * src_frequency_sd;
    offset_weight = 1.0 / sources[index]->sel_info.root_distance;
    frequency_weight = 1.0 / SQUARE(src_frequency_sd);

    DEBUG_LOG("combining index=%d oweight=%e offset=%e osd=%e fweight=%e freq=%e fsd=%e skew=%e",
              index, offset_weight, src_offset, src_offset_sd,
              frequency_weight, src_frequency, src_frequency_sd, src_skew);

    sum_offset_weight += offset_weight;
    sum_offset += offset_weight * src_offset;
    sum2_offset_sd += offset_weight * (SQUARE(src_offset_sd) +
                                       SQUARE(src_offset - *offset));

    sum_frequency_weight += frequency_weight;
    sum_frequency += frequency_weight * src_frequency;
    inv_sum2_frequency_sd += 1.0 / SQUARE(src_frequency_sd);
    inv_sum2_skew += 1.0 / SQUARE(src_skew);

    combined++;
  }

  assert(combined);
  *offset = sum_offset / sum_offset_weight;
  *offset_sd = sqrt(sum2_offset_sd / sum_offset_weight);
  *frequency = sum_frequency / sum_frequency_weight;
  *frequency_sd = 1.0 / sqrt(inv_sum2_frequency_sd);
  *skew = 1.0 / sqrt(inv_sum2_skew);

  DEBUG_LOG("combined result offset=%e osd=%e freq=%e fsd=%e skew=%e",
            *offset, *offset_sd, *frequency, *frequency_sd, *skew);

  return combined;
}