Include offset correction error in dispersion

diff --git a/chrony.texi b/chrony.texi
index 51d59e06155c9e55124a038b14a4abd444b04dc0..60f2924b39dc9f1d59e77a545199072bd0dd38d1 100644 (file)
--- a/chrony.texi
+++ b/chrony.texi
@@ -2002,7 +2002,7 @@ An example line (which actually appears as a single line in the file)
 from the refclocks log file is shown below.
 
 @example
-2009-11-30 14:33:27.000000 PPS2    7 N 1  4.900000e-07 -6.741777e-07
+2009-11-30 14:33:27.000000 PPS2    7 N 1  4.900000e-07 -6.741777e-07  1.000e-06
 @end example
 
 The columns are as follows (the quantities in square brackets are the
@@ -2030,6 +2030,8 @@ Local clock error measured by refclock driver. [4.900000e-07]
 @item
 Local clock error with applied corrections.  Positive indicates
 that the local clock is slow. [-6.741777e-07]
+@item
+Assumed dispersion of the sample. [1.000e-06]
 @end enumerate
 
 A banner is periodically written to the log file to indicate the

diff --git a/ntp_core.c b/ntp_core.c
index 7b2d7174d47f656a9049b956bc64b066a844fb2f..c0eada3d53ab7477c97baf14e3a2ef1af98fb9a2 100644 (file)
--- a/ntp_core.c
+++ b/ntp_core.c
@@ -754,7 +754,7 @@ transmit_timeout(void *arg)
 /* ================================================== */
 
 static void
-receive_packet(NTP_Packet *message, struct timeval *now, NCR_Instance inst, int do_auth)
+receive_packet(NTP_Packet *message, struct timeval *now, double now_err, NCR_Instance inst, int do_auth)
 {
   int pkt_leap;
   int source_is_synchronized;
@@ -928,7 +928,7 @@ receive_packet(NTP_Packet *message, struct timeval *now, NCR_Instance inst, int
     skew = source_freq_hi - source_freq_lo;
     
     /* and then calculate peer dispersion */
-    epsilon = LCL_GetSysPrecisionAsQuantum() + skew * local_interval;
+    epsilon = LCL_GetSysPrecisionAsQuantum() + now_err + skew * local_interval;
     
   } else {
     /* If test3 failed, we probably can't calculate these quantities
@@ -1318,6 +1318,7 @@ static void
 process_known
 (NTP_Packet *message,           /* the received message */
  struct timeval *now,           /* timestamp at time of receipt */
+ double now_err,
  NCR_Instance inst,             /* the instance record for this peer/server */
  int do_auth                   /* whether the received packet allegedly contains
                                    authentication info*/
@@ -1406,7 +1407,7 @@ process_known
         case MODE_ACTIVE:
           /* Ordinary symmetric peering */
           CLG_LogNTPPeerAccess(&inst->remote_addr.ip_addr, (time_t) now->tv_sec);
-          receive_packet(message, now, inst, do_auth);
+          receive_packet(message, now, now_err, inst, do_auth);
           break;
         case MODE_PASSIVE:
           /* In this software this case should not arise, we don't
@@ -1416,7 +1417,7 @@ process_known
           /* This is where we have the remote configured as a server and he has
              us configured as a peer - fair enough. */
           CLG_LogNTPPeerAccess(&inst->remote_addr.ip_addr, (time_t) now->tv_sec);
-          receive_packet(message, now, inst, do_auth);
+          receive_packet(message, now, now_err, inst, do_auth);
           break;
         case MODE_SERVER:
           /* Nonsense - we can't have a preconfigured server */
@@ -1437,14 +1438,14 @@ process_known
         case MODE_ACTIVE:
           /* Slightly bizarre combination, but we can still process it */
           CLG_LogNTPPeerAccess(&inst->remote_addr.ip_addr, (time_t) now->tv_sec);
-          receive_packet(message, now, inst, do_auth);
+          receive_packet(message, now, now_err, inst, do_auth);
           break;
         case MODE_PASSIVE:
           /* We have no passive peers in this software */
           break;
         case MODE_CLIENT:
           /* Standard case where he's a server and we're the client */
-          receive_packet(message, now, inst, do_auth);
+          receive_packet(message, now, now_err, inst, do_auth);
           break;
         case MODE_SERVER:
           /* RFC1305 error condition. */
@@ -1465,7 +1466,7 @@ process_known
           /* This would arise if we have the remote configured as a peer and
              he does not have us configured */
           CLG_LogNTPPeerAccess(&inst->remote_addr.ip_addr, (time_t) now->tv_sec);
-          receive_packet(message, now, inst, do_auth);
+          receive_packet(message, now, now_err, inst, do_auth);
           break;
         case MODE_PASSIVE:
           /* Error condition in RFC1305.  Also, we can't have any
@@ -1474,7 +1475,7 @@ process_known
           break;
         case MODE_CLIENT:
           /* This is a wierd combination - how could it arise? */
-          receive_packet(message, now, inst, do_auth);
+          receive_packet(message, now, now_err, inst, do_auth);
           break;
         case MODE_SERVER:
           /* Error condition in RFC1305 */
@@ -1507,10 +1508,10 @@ process_known
    and it relates to a source we have an ongoing protocol exchange with */
 
 void
-NCR_ProcessNoauthKnown(NTP_Packet *message, struct timeval *now, NCR_Instance inst)
+NCR_ProcessNoauthKnown(NTP_Packet *message, struct timeval *now, double now_err, NCR_Instance inst)
 {
 
-  process_known(message, now, inst, 0);
+  process_known(message, now, now_err, inst, 0);
 
 }
 
@@ -1519,7 +1520,7 @@ NCR_ProcessNoauthKnown(NTP_Packet *message, struct timeval *now, NCR_Instance in
    and we do not recognize its source */
 
 void
-NCR_ProcessNoauthUnknown(NTP_Packet *message, struct timeval *now, NTP_Remote_Address *remote_addr)
+NCR_ProcessNoauthUnknown(NTP_Packet *message, struct timeval *now, double now_err, NTP_Remote_Address *remote_addr)
 {
 
   NTP_Mode his_mode;
@@ -1576,9 +1577,9 @@ NCR_ProcessNoauthUnknown(NTP_Packet *message, struct timeval *now, NTP_Remote_Ad
    exchange with */
 
 void
-NCR_ProcessAuthKnown(NTP_Packet *message, struct timeval *now, NCR_Instance data)
+NCR_ProcessAuthKnown(NTP_Packet *message, struct timeval *now, double now_err, NCR_Instance data)
 {
-  process_known(message, now, data, 1);
+  process_known(message, now, now_err, data, 1);
 
 }
 
@@ -1587,7 +1588,7 @@ NCR_ProcessAuthKnown(NTP_Packet *message, struct timeval *now, NCR_Instance data
    the network, and we do not recognize its source */
 
 void
-NCR_ProcessAuthUnknown(NTP_Packet *message, struct timeval *now, NTP_Remote_Address *remote_addr)
+NCR_ProcessAuthUnknown(NTP_Packet *message, struct timeval *now, double now_err, NTP_Remote_Address *remote_addr)
 {
 
   NTP_Mode his_mode;

diff --git a/ntp_core.h b/ntp_core.h
index fa127863eb8caab7f25ddbb4d85aea79d7e1ce37..ee22e466705d17d7fed66e5df71bd21cceaff77a 100644 (file)
--- a/ntp_core.h
+++ b/ntp_core.h
@@ -57,20 +57,20 @@ extern void NCR_DestroyInstance(NCR_Instance instance);
 
 /* This routine is called when a new packet arrives off the network,
    and it relates to a source we have an ongoing protocol exchange with */
-extern void NCR_ProcessNoauthKnown(NTP_Packet *message, struct timeval *now, NCR_Instance data);
+extern void NCR_ProcessNoauthKnown(NTP_Packet *message, struct timeval *now, double now_err, NCR_Instance data);
 
 /* This routine is called when a new packet arrives off the network,
    and we do not recognize its source */
-extern void NCR_ProcessNoauthUnknown(NTP_Packet *message, struct timeval *now, NTP_Remote_Address *remote_addr);
+extern void NCR_ProcessNoauthUnknown(NTP_Packet *message, struct timeval *now, double now_err, NTP_Remote_Address *remote_addr);
 
 /* This routine is called when a new authenticated packet arrives off
    the network, and it relates to a source we have an ongoing protocol
    exchange with */
-extern void NCR_ProcessAuthKnown(NTP_Packet *message, struct timeval *now, NCR_Instance data);
+extern void NCR_ProcessAuthKnown(NTP_Packet *message, struct timeval *now, double now_err, NCR_Instance data);
 
 /* This routine is called when a new authenticated packet arrives off
    the network, and we do not recognize its source */
-extern void NCR_ProcessAuthUnknown(NTP_Packet *message, struct timeval *now, NTP_Remote_Address *remote_addr);
+extern void NCR_ProcessAuthUnknown(NTP_Packet *message, struct timeval *now, double now_err, NTP_Remote_Address *remote_addr);
 
 /* Slew receive and transmit times in instance records */
 extern void NCR_SlewTimes(NCR_Instance inst, struct timeval *when, double dfreq, double doffset);

diff --git a/ntp_io.c b/ntp_io.c
index 1841b9c585831fecf4eb061149ccefcd1253890b..7ccc459049d3909e5422c9f1cf968bf5913f0318 100644 (file)
--- a/ntp_io.c
+++ b/ntp_io.c
@@ -356,11 +356,11 @@ read_from_socket(void *anything)
 
     if (status == NTP_NORMAL_PACKET_SIZE) {
 
-      NSR_ProcessReceive((NTP_Packet *) &message.ntp_pkt, &now, &remote_addr);
+      NSR_ProcessReceive((NTP_Packet *) &message.ntp_pkt, &now, now_err, &remote_addr);
 
     } else if (status == sizeof(NTP_Packet)) {
 
-      NSR_ProcessAuthenticatedReceive((NTP_Packet *) &message.ntp_pkt, &now, &remote_addr);
+      NSR_ProcessAuthenticatedReceive((NTP_Packet *) &message.ntp_pkt, &now, now_err, &remote_addr);
 
     } else {
 

diff --git a/ntp_sources.c b/ntp_sources.c
index caf27ac3af129ee774fa758b11f3e49bb863c100..9a058f3d0a0cac1fd731a139db44618b95bcc746 100644 (file)
--- a/ntp_sources.c
+++ b/ntp_sources.c
@@ -264,7 +264,7 @@ NSR_RemoveSource(NTP_Remote_Address *remote_addr)
 
 /* This routine is called by ntp_io when a new packet arrives off the network.*/
 void
-NSR_ProcessReceive(NTP_Packet *message, struct timeval *now, NTP_Remote_Address *remote_addr)
+NSR_ProcessReceive(NTP_Packet *message, struct timeval *now, double now_err, NTP_Remote_Address *remote_addr)
 {
   int slot, found;
 
@@ -278,9 +278,9 @@ NSR_ProcessReceive(NTP_Packet *message, struct timeval *now, NTP_Remote_Address
   
   find_slot(remote_addr, &slot, &found);
   if (found == 2) { /* Must match IP address AND port number */
-    NCR_ProcessNoauthKnown(message, now, records[slot].data);
+    NCR_ProcessNoauthKnown(message, now, now_err, records[slot].data);
   } else {
-    NCR_ProcessNoauthUnknown(message, now, remote_addr);
+    NCR_ProcessNoauthUnknown(message, now, now_err, remote_addr);
   }
 }
 
@@ -288,7 +288,7 @@ NSR_ProcessReceive(NTP_Packet *message, struct timeval *now, NTP_Remote_Address
 
 /* This routine is called by ntp_io when a new packet with an authentication tail arrives off the network */
 void
-NSR_ProcessAuthenticatedReceive(NTP_Packet *message, struct timeval *now, NTP_Remote_Address *remote_addr)
+NSR_ProcessAuthenticatedReceive(NTP_Packet *message, struct timeval *now, double now_err, NTP_Remote_Address *remote_addr)
 {
   int slot, found;
 
@@ -296,9 +296,9 @@ NSR_ProcessAuthenticatedReceive(NTP_Packet *message, struct timeval *now, NTP_Re
 
   find_slot(remote_addr, &slot, &found);
   if (found == 2) {
-    NCR_ProcessAuthKnown(message, now, records[slot].data);
+    NCR_ProcessAuthKnown(message, now, now_err, records[slot].data);
   } else {
-    NCR_ProcessAuthUnknown(message, now, remote_addr);
+    NCR_ProcessAuthUnknown(message, now, now_err, remote_addr);
   }
 }
 

diff --git a/ntp_sources.h b/ntp_sources.h
index 5cfcac6d55bdd89ebde3c1f7ccbaf689649914e9..d708949843d4f8b67be8ce955d5a8ab948a3598a 100644 (file)
--- a/ntp_sources.h
+++ b/ntp_sources.h
@@ -62,10 +62,10 @@ extern NSR_Status NSR_AddPeer(NTP_Remote_Address *remote_addr, SourceParameters
 extern NSR_Status NSR_RemoveSource(NTP_Remote_Address *remote_addr);
 
 /* This routine is called by ntp_io when a new packet arrives off the network */
-extern void NSR_ProcessReceive(NTP_Packet *message, struct timeval *now, NTP_Remote_Address *remote_addr);
+extern void NSR_ProcessReceive(NTP_Packet *message, struct timeval *now, double now_err, NTP_Remote_Address *remote_addr);
 
 /* This routine is called by ntp_io when a new packet with an authentication tail arrives off the network */
-extern void NSR_ProcessAuthenticatedReceive(NTP_Packet *message, struct timeval *now, NTP_Remote_Address *remote_addr);
+extern void NSR_ProcessAuthenticatedReceive(NTP_Packet *message, struct timeval *now, double now_err, NTP_Remote_Address *remote_addr);
 
 /* Initialisation function */
 extern void NSR_Initialise(void);

diff --git a/refclock.c b/refclock.c
index f4e42a1bb38832ed040d7ac39401b030813301b5..1a99f53e67c8d49baacb2106da73a6b733de0bb1 100644 (file)
--- a/refclock.c
+++ b/refclock.c
@@ -43,6 +43,7 @@ extern RefclockDriver RCL_PPS_driver;
 
 struct FilterSample {
   double offset;
+  double dispersion;
   struct timeval sample_time;
 };
 
@@ -52,6 +53,7 @@ struct MedianFilter {
   int used;
   int last;
   struct FilterSample *samples;
+  int *sort_array;
 };
 
 struct RCL_Instance_Record {
@@ -89,13 +91,13 @@ static int pps_stratum(RCL_Instance instance, struct timeval *tv);
 static void poll_timeout(void *arg);
 static void slew_samples(struct timeval *raw, struct timeval *cooked, double dfreq, double afreq,
              double doffset, int is_step_change, void *anything);
-static void log_sample(RCL_Instance instance, struct timeval *sample_time, int pulse, double raw_offset, double cooked_offset);
+static void log_sample(RCL_Instance instance, struct timeval *sample_time, int pulse, double raw_offset, double cooked_offset, double dispersion);
 
 static void filter_init(struct MedianFilter *filter, int length);
 static void filter_fini(struct MedianFilter *filter);
 static void filter_reset(struct MedianFilter *filter);
-static void filter_add_sample(struct MedianFilter *filter, struct timeval *sample_time, double offset);
-static int filter_get_last_sample(struct MedianFilter *filter, struct timeval *sample_time, double *offset);
+static void filter_add_sample(struct MedianFilter *filter, struct timeval *sample_time, double offset, double dispersion);
+static int filter_get_last_sample(struct MedianFilter *filter, struct timeval *sample_time, double *offset, double *dispersion);
 static int filter_get_sample(struct MedianFilter *filter, struct timeval *sample_time, double *offset, double *dispersion);
 static void filter_slew_samples(struct MedianFilter *filter, struct timeval *when, double dfreq, double doffset);
 
@@ -317,11 +319,12 @@ RCL_GetDriverOption(RCL_Instance instance, char *name)
 int
 RCL_AddSample(RCL_Instance instance, struct timeval *sample_time, double offset, NTP_Leap leap_status)
 {
-  double correction, err;
+  double correction, dispersion;
   struct timeval cooked_time;
 
-  LCL_GetOffsetCorrection(sample_time, &correction, &err);
+  LCL_GetOffsetCorrection(sample_time, &correction, &dispersion);
   UTI_AddDoubleToTimeval(sample_time, correction, &cooked_time);
+  dispersion += LCL_GetSysPrecisionAsQuantum();
 
   if (!valid_sample_time(instance, sample_time))
     return 0;
@@ -331,10 +334,10 @@ RCL_AddSample(RCL_Instance instance, struct timeval *sample_time, double offset,
       offset, offset - correction + instance->offset);
 #endif
 
-  filter_add_sample(&instance->filter, &cooked_time, offset - correction + instance->offset);
+  filter_add_sample(&instance->filter, &cooked_time, offset - correction + instance->offset, dispersion);
   instance->leap_status = leap_status;
 
-  log_sample(instance, &cooked_time, 0, offset, offset - correction + instance->offset);
+  log_sample(instance, &cooked_time, 0, offset, offset - correction + instance->offset, dispersion);
 
   return 1;
 }
@@ -342,18 +345,13 @@ RCL_AddSample(RCL_Instance instance, struct timeval *sample_time, double offset,
 int
 RCL_AddPulse(RCL_Instance instance, struct timeval *pulse_time, double second)
 {
-  double correction, err, offset;
+  double correction, dispersion, offset;
   struct timeval cooked_time;
   int rate;
 
-  struct timeval ref_time;
-  int is_synchronised, stratum;
-  double root_delay, root_dispersion, distance;
-  NTP_Leap leap;
-  unsigned long ref_id;
-
-  LCL_GetOffsetCorrection(pulse_time, &correction, &err);
+  LCL_GetOffsetCorrection(pulse_time, &correction, &dispersion);
   UTI_AddDoubleToTimeval(pulse_time, correction, &cooked_time);
+  dispersion += LCL_GetSysPrecisionAsQuantum();
 
   if (!valid_sample_time(instance, pulse_time))
     return 0;
@@ -372,10 +370,10 @@ RCL_AddPulse(RCL_Instance instance, struct timeval *pulse_time, double second)
 
   if (instance->lock_ref != -1) {
     struct timeval ref_sample_time;
-    double sample_diff, ref_offset, shift;
+    double sample_diff, ref_offset, ref_dispersion, shift;
 
     if (!filter_get_last_sample(&refclocks[instance->lock_ref].filter,
-          &ref_sample_time, &ref_offset))
+          &ref_sample_time, &ref_offset, &ref_dispersion))
       return 0;
 
     UTI_DiffTimevalsToDouble(&sample_diff, &cooked_time, &ref_sample_time);
@@ -390,7 +388,7 @@ RCL_AddPulse(RCL_Instance instance, struct timeval *pulse_time, double second)
 
     offset += shift;
 
-    if (fabs(ref_offset - offset) >= 0.2 / rate)
+    if (fabs(ref_offset - offset) + ref_dispersion + dispersion >= 0.2 / rate)
       return 0;
 
 #if 0
@@ -398,6 +396,12 @@ RCL_AddPulse(RCL_Instance instance, struct timeval *pulse_time, double second)
         second, offset, ref_offset - offset, sample_diff);
 #endif
   } else {
+    struct timeval ref_time;
+    int is_synchronised, stratum;
+    double root_delay, root_dispersion, distance;
+    NTP_Leap leap;
+    unsigned long ref_id;
+
     /* Ignore the pulse if we are not well synchronized */
 
     REF_GetReferenceParams(&cooked_time, &is_synchronised, &leap, &stratum,
@@ -420,10 +424,10 @@ RCL_AddPulse(RCL_Instance instance, struct timeval *pulse_time, double second)
       second, offset);
 #endif
 
-  filter_add_sample(&instance->filter, &cooked_time, offset);
+  filter_add_sample(&instance->filter, &cooked_time, offset, dispersion);
   instance->leap_status = LEAP_Normal;
 
-  log_sample(instance, &cooked_time, 1, second, offset);
+  log_sample(instance, &cooked_time, 1, second, offset, dispersion);
 
   return 1;
 }
@@ -545,7 +549,7 @@ slew_samples(struct timeval *raw, struct timeval *cooked, double dfreq, double a
 }
 
 static void
-log_sample(RCL_Instance instance, struct timeval *sample_time, int pulse, double raw_offset, double cooked_offset)
+log_sample(RCL_Instance instance, struct timeval *sample_time, int pulse, double raw_offset, double cooked_offset, double dispersion)
 {
   char sync_stats[4] = {'N', '+', '-', '?'};
 
@@ -564,11 +568,11 @@ log_sample(RCL_Instance instance, struct timeval *sample_time, int pulse, double
 
   if (((logwrites++) % 32) == 0) {
     fprintf(logfile,
-        "====================================================================\n"
-        "   Date (UTC) Time         Refid  DP L P  Raw offset   Cooked offset\n"
-        "====================================================================\n");
+        "===============================================================================\n"
+        "   Date (UTC) Time         Refid  DP L P  Raw offset   Cooked offset      Disp.\n"
+        "===============================================================================\n");
   }
-  fprintf(logfile, "%s.%06d %-5s %3d %1c %1d %13.6e %13.6e\n",
+  fprintf(logfile, "%s.%06d %-5s %3d %1c %1d %13.6e %13.6e %10.3e\n",
       UTI_TimeToLogForm(sample_time->tv_sec),
       (int)sample_time->tv_usec,
       UTI_RefidToString(instance->ref_id),
@@ -576,7 +580,8 @@ log_sample(RCL_Instance instance, struct timeval *sample_time, int pulse, double
       sync_stats[instance->leap_status],
       pulse,
       raw_offset,
-      cooked_offset);
+      cooked_offset,
+      dispersion);
   fflush(logfile);
 }
 
@@ -591,12 +596,14 @@ filter_init(struct MedianFilter *filter, int length)
   filter->used = 0;
   filter->last = -1;
   filter->samples = MallocArray(struct FilterSample, filter->length);
+  filter->sort_array = MallocArray(int, filter->length);
 }
 
 static void
 filter_fini(struct MedianFilter *filter)
 {
   Free(filter->samples);
+  Free(filter->sort_array);
 }
 
 static void
@@ -607,7 +614,7 @@ filter_reset(struct MedianFilter *filter)
 }
 
 static void
-filter_add_sample(struct MedianFilter *filter, struct timeval *sample_time, double offset)
+filter_add_sample(struct MedianFilter *filter, struct timeval *sample_time, double offset, double dispersion)
 {
   filter->index++;
   filter->index %= filter->length;
@@ -617,23 +624,30 @@ filter_add_sample(struct MedianFilter *filter, struct timeval *sample_time, doub
 
   filter->samples[filter->index].sample_time = *sample_time;
   filter->samples[filter->index].offset = offset;
+  filter->samples[filter->index].dispersion = dispersion;
 }
 
 static int
-filter_get_last_sample(struct MedianFilter *filter, struct timeval *sample_time, double *offset)
+filter_get_last_sample(struct MedianFilter *filter, struct timeval *sample_time, double *offset, double *dispersion)
 {
   if (filter->last < 0)
     return 0;
 
   *sample_time = filter->samples[filter->last].sample_time;
   *offset = filter->samples[filter->last].offset;
+  *dispersion = filter->samples[filter->last].dispersion;
   return 1;
 }
 
+static const struct FilterSample *tmp_sorted_array;
+
 static int
 sample_compare(const void *a, const void *b)
 {
-  const struct FilterSample *s1 = a, *s2 = b;
+  const struct FilterSample *s1, *s2;
+
+  s1 = &tmp_sorted_array[*(int *)a];
+  s2 = &tmp_sorted_array[*(int *)b];
 
   if (s1->offset < s2->offset)
     return -1;
@@ -651,44 +665,66 @@ filter_get_sample(struct MedianFilter *filter, struct timeval *sample_time, doub
   if (filter->used == 1) {
     *sample_time = filter->samples[filter->index].sample_time;
     *offset = filter->samples[filter->index].offset;
-    *dispersion = 0.0;
+    *dispersion = filter->samples[filter->index].dispersion;
   } else {
-    int i, from, to;
-    double x, x1, y, d;
+    struct FilterSample *s;
+    int i, j, from, to;
+    double x, x1, y, d, e, min_dispersion;
+
+    /* find minimum dispersion */
+    for (i = 1, min_dispersion = filter->samples[0].dispersion; i < filter->used; i++) {
+      if (min_dispersion > filter->samples[i].dispersion)
+        min_dispersion = filter->samples[i].dispersion;
+    }
 
-    /* sort samples by offset */
-    qsort(filter->samples, filter->used, sizeof (struct FilterSample), sample_compare);
+    /* select samples with dispersion better than 1.5 * minimum */
+    for (i = j = 0; i < filter->used; i++) {
+      if (filter->samples[i].dispersion <= 1.5 * min_dispersion)
+        filter->sort_array[j++] = i;
+    }
 
-    /* average the half of the samples closest to the median */ 
-    if (filter->used > 2) {
-      from = (filter->used + 2) / 4;
-      to = filter->used - from;
+    assert(j > 0);
+    
+    /* and sort their indexes by offset */
+    tmp_sorted_array = filter->samples;
+    qsort(filter->sort_array, j, sizeof (int), sample_compare);
+
+    /* select half of the samples closest to the median */ 
+    if (j > 2) {
+      from = (j + 2) / 4;
+      to = j - from;
     } else {
       from = 0;
-      to = filter->used;
+      to = j;
     }
 
+    /* average offset and sample time */ 
     for (i = from, x = y = 0.0; i < to; i++) {
+      s = &filter->samples[filter->sort_array[i]];
 #if 0
-      LOG(LOGS_INFO, LOGF_Refclock, "refclock averaging offset %.9f [%s]",
-          filter->samples[i].offset, UTI_TimevalToString(&filter->samples[i].sample_time));
+      LOG(LOGS_INFO, LOGF_Refclock, "refclock averaging sample: offset %.9f dispersion %.9f [%s]",
+          s->offset, s->dispersion, UTI_TimevalToString(&filter->samples[i].sample_time));
 #endif
-      UTI_DiffTimevalsToDouble(&x1, &filter->samples[i].sample_time, &filter->samples[0].sample_time);
+      UTI_DiffTimevalsToDouble(&x1, &s->sample_time, &filter->samples[0].sample_time);
       x += x1;
-      y += filter->samples[i].offset;
+      y += s->offset;
     }
 
     x /= to - from;
     y /= to - from;
 
-    for (i = from, d = 0.0; i < to; i++)
-      d += (filter->samples[i].offset - y) * (filter->samples[i].offset - y);
+    for (i = from, d = e = 0.0; i < to; i++) {
+      s = &filter->samples[filter->sort_array[i]];
+      d += (s->offset - y) * (s->offset - y);
+      e += s->dispersion;
+    }
 
     d = sqrt(d / (to - from));
+    e /= to - from;
 
     UTI_AddDoubleToTimeval(&filter->samples[0].sample_time, x, sample_time);
     *offset = y;
-    *dispersion = d;
+    *dispersion = d + e;
   }
 
   return 1;