home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
nisttime.carsoncity.k12.mi.us
/
nisttime.carsoncity.k12.mi.us.tar
/
nisttime.carsoncity.k12.mi.us
/
pub
/
autolock
/
getdif.c
< prev
next >
Wrap
C/C++ Source or Header
|
2003-12-09
|
14KB
|
423 lines
#include <stdio.h>
#include <sys/types.h>
#include <fcntl.h>
#include <sys/stat.h>
#include <sys/socket.h>
#include <sys/time.h>
#include <netinet/in.h>
#include <arpa/inet.h>
#include <string.h>
#include <ctype.h>
#include <signal.h>
#include <math.h>
#ifndef LINUX
#include <stropts.h>
#endif
#include "autolock.h"
int timeout; /*set when SIGALRM happens*/
int getdif(char *cp,struct autolockst *lock,int *lsr)
{
/*
this subroutine estimates the difference between the
local clock and a server. The comparison is made using
NTP-format messages.
the integer lsr is set to the leap second flag received
from the remote host.
the possible values are:
0 -- server ok, no leap second
1 -- server ok, insert second
2 -- server ok, drop second
3 -- server not ok.
if the server returns a value of 3, this subroutine
returns immediately with a status of -4. the other
3 possibilities are simply passed back to the caller.
the subroutine attempts to estimate the uncertainty in the
NTP time exchange and tries to adjust its performance so that
the measured uncertainty is consistent with the desired performance
as specified in the global parameters described below.
This software was developed with US Government support
and it may not be sold, restricted or licensed. You
may duplicate this program provided that this notice
remains in all of the copies, and you may give it to
others provided they understand and agree to this
condition.
This program and the time protocol it uses are under
development and the implementation may change without
notice.
For questions or additional information, contact:
Judah Levine
JILA
Campus Box 440
University of Colorado
Boulder, Colorado 80309
(303) 492 7785
jlevine@utcnist.colorado.edu
*/
#include "ntp-packet.h"
struct ntppkt msg; /* the message we send out*/
struct ntppkt ans; /* the response*/
struct timeval tvv; /* get local system time*/
int address; /* holds ip address */
int pserv = 123; /* ntp port number on server */
unsigned int ts,fts; /* transmission time of outbound packet*/
unsigned int tr,ftr; /* time and which response was received*/
unsigned int tx,ftx; /* remote transmission and received times*/
unsigned int trr,ftrr;
long int dly; /* one-way delay in microseconds*/
long int tdiff[NCALMAX]; /* time difference in microseconds*/
/*
ntpdif computes the difference between two times
expressed in NTP format
*/
long int ntpdif (unsigned int, unsigned int, unsigned int, unsigned int);
long int median (long int[], int);
void smooth(long int[], int, double *, double *);
int oops(); /*receives time-out alarms*/
extern int ncal; /*number of pings to a server*/
double val,rms;
struct sockaddr_in sin; /* socket address structure */
int s; /* socket number */
int length; /* size of message */
int i,j;
int kval; /*number of packets actually received.*/
int leapind; /*received leap indicator*/
int mjd; /*cycle time converted to MJD */
int secoday; /*the time in seconds after 0000*/
int msosec; /*the fraction of a second in milliseconds*/
signal(SIGALRM,(void *) &oops); /*set routine to handle time-outs*/
/*
convert address to internal format
then store address and port number in socket structure
*/
if( (address=inet_addr(cp) ) == -1) return(-1); /*bad address*/
bzero( (char *)&sin,sizeof(sin));
sin.sin_addr.s_addr=address;
sin.sin_family=AF_INET;
sin.sin_port=htons(pserv);
/*
create the socket then connect to the server.
return error status of -2 if socket creation fails
and error status of -3 if connetion fails.
*/
if( (s=socket(AF_INET,SOCK_DGRAM,0) ) < 0) return(-2);
if(connect(s, (struct sockaddr *) &sin, sizeof(sin) ) < 0) return(-3);
kval=0; /*index showing number of measurements*/
for(j=0; j<ncal; j++) /*ping 'em again and again*/
{
/*
build an NTP packet and send it off.
*/
msg.li_vn_mode=033; /* 00 011 011 = no leap, v3, client*/
msg.stratum= 0; /* 0 means unspecified*/
msg.ppoll= 8; /* poll interval = 256s*/
msg.precision= 0; /* 0= 1s, no precision at all here*/
msg.rootdelay=0; /* don't really know or care*/
msg.rootdispersion=0; /* ditto*/
msg.refid = 0; /* we don't have a reference*/
/*
system request returns time as seconds and microseconds since
1970.0. Must convert this to NTP time format, which uses
1900.0 as the origin and fractions of a second instead of
microseconds. Thus must add seconds 1900.0 -> 1970.0 to
system time in seconds since 1970.
The fractional part divides 1 second into 2**32 parts, so
that one least count is 232.8306 ps. Thus the fraction is
system time in microseconds * 10**6/232.8306 =
system time * 4294.9673
*/
gettimeofday(&tvv,(struct timezone *)NULL);
ts= 2208988800u + (unsigned int) tvv.tv_sec;
fts= (unsigned int) ( 4294.9673*(double)tvv.tv_usec);
msg.reftime.sys[0]=htonl(ts-1);
msg.reftime.sys[1]=htonl(fts);
msg.org.sys[0]= htonl(ts);
msg.org.sys[1]= htonl(fts);
alarm(5); /*5 seconds -- should be long enough*/
timeout=0;
length=write(s,&msg,sizeof(msg));
/*
read response and get local time as soon as read completes.
*/
do
{
length=recvfrom(s,&ans,sizeof(ans),0,
(struct sockaddr *) NULL,(size_t *) NULL);
gettimeofday(&tvv,(struct timezone *)NULL);
} while ( (length < sizeof(ans)) && (timeout == 0) );
alarm(0); /*turn off the timer*/
if(timeout != 0) continue; /* try again if time out*/
tr= 2208988800u + (unsigned int) tvv.tv_sec;
ftr= (unsigned int) (4294.9673*(double)tvv.tv_usec);
leapind=ans.li_vn_mode; /*copy to an integer*/
leapind >>= 6;
leapind &= 3; /*isolate leap indicator bits*/
if(leapind == 3) return(-4); /*server is not healthy*/
*lsr= leapind; /*send leap second bits to caller*/
trr=ntohl(ans.rec.sys[0]);
ftrr=ntohl(ans.rec.sys[1]);
tx=ntohl(ans.xmt.sys[0]);
ftx=ntohl(ans.xmt.sys[1]);
/*
dly is total elapsed time as measured on our local clock
- time spent in remote server, which is the transmit
time of the remote server - the received time.
*/
dly= ntpdif(tr,ftr,ts,fts) - ntpdif(tx,ftx,trr,ftrr);
tdiff[kval++]=ntpdif(ts,fts,trr,ftrr) + dly/2; /*local- remote*/
} /* end of loop over j */
close(s); /*done with the network connection*/
if(kval == 0) return(-5); /* nobody home ? */
if(kval == 1)
{
val=tdiff[0]; /*only one value*/
rms= 999; /*so rms has no meaning*/
}
else if(kval == 2)
{
val= (tdiff[0]+tdiff[1])/2;
rms= fabs(tdiff[0]-tdiff[1])/1.414;
}
else /*more than 2 responses*/
{
smooth(tdiff,kval,&val,&rms);
}
mjd=ts/86400 + 15020; /*compute MJD*/
secoday= ts%86400;
msosec= ((double) fts)/4294968.1;
lock->day=mjd;
lock->fday= 1000*secoday + msosec;
lock->x= 1e-6*val; /*convert from usec to sec*/
lock->rms= 1e-6*rms; /*ditto*/
lock->nval=kval;
return(1);
}
long int ntpdif(unsigned int t1,unsigned int ft1,
unsigned int t2,unsigned int ft2)
{
/*
this subroutine computes the difference between the first
and second pairs of NTP-format times and returns the
difference in microseconds.
The input parameter pairs are seconds and fraction of
a second, respectively, in standard NTP format.
Since the fraction divides 1 second into 2**32 parts,
a least-count of the fraction is 232.8306 ps. The
value in microseconds is therefore (value)*232.83/1000000
note that NTP-format numbers are unsigned so that
subtractions will produce strange results if performed
directly.
*/
long int temp;
long int ta;
double tb;
ta= (long int) t1 - (long int) t2;
tb= (double) ( (long int) ft1 - (long int) ft2 );
temp= 1000000*ta +
(long int)(tb/4294.9673);
return (temp);
}
long int median(long int ivals[],int indx)
{
/*
this subroutine computes the median of the input
array ivals. There are indx values in the array.
the median is returned and the input data are
not altered.
the median is that number for which there are imid
entries both bigger and smaller than it. there is a
single index value that will satisfy this requirement
if the number of entries is odd and two adjacent ones
will satisfy it if the number is even. the average of
these two adjacent values are returned in this case.
if there are a significant number of repeated values,
then no entry may satisfy this requirement. this case
can be dealt with approximately by finding an entry that
approximately satisfies the test, but the median may not
do the "right" thing in this case. this subroutine is
complicated by the facts that there are likely to be a
significant number of duplicates, that the input array
must be left alone and that the subroutine is in the
main loop and must be as fast as possible.
*/
int j;
int imid; /*index of mid point value */
char iflag[NCALMAX]; /*static array for speed */
int imin; /*index of minimum value*/
long int mmin; /*current minimum value*/
int ifound =0; /*index number of min just found */
int iodd; /*1 if indx is even, 0 if odd */
long int m1; /*lower median if size is even*/
iodd= indx & 1;
imid=indx/2; /*index of mid-point value for odd size*/
for(j=0; j<indx; j++) iflag[j]=0;
/*
loop over all elements and find the minimum one.
set the flag to show that it has been found and
then loop again to find the minimum one in the
ones that are left. This process will always
converge even if all of the values are the same,
but it may not do the "right" thing if there are
only a few distinct values that appear many times.
first find first value that has not yet been tagged
as a minimum in some round.
*/
loop:
for(imin=0; imin<indx; imin++)
if(iflag[imin] == 0)
{
mmin=ivals[imin];
break;
}
/*
now compare this value to all other untagged
ones to find the minimum of the lot. the surviving
minimum would go in the next place in an output array
that was ordered by size. its index in that array would
be ifound. if there are an odd number of points and
ifound=imid then we have found the median in this
minimum and we are finished. if there are an even
number of points, then we must average the values
found before and after imid to compute the median.
if this is the one before then just save it. if
this is the one after then compute average and
return.
if ifound is too small in either case, then
we have to go around again. tag the surviving minimum
and go around to find the minimum of the remaining
untagged folks.
*/
for(j=imin+1; j<indx; j++)
if( (iflag[j] == 0) && (ivals[j] < mmin) )
{
mmin=ivals[j];
imin=j;
}
if(iodd == 1)
{
if (ifound == imid) return (mmin);
}
else /*size is even*/
{
if(ifound == (imid - 1) ) m1=mmin;
if(ifound == (imid + 1) ) return( ( (m1+mmin) )/2 );
}
ifound++;
iflag[imin]=1;
goto loop;
}
void smooth(long int ivals[],int indx,double *sval,double *rms)
{
/*
this subroutine smooths the input array ivals which has
indx members. the smoothed value is returned in sval
and the rms deviation is returned in rms.
the smoothing consists of first throwing away all values that
are more than 10 ms = 10000 us away from the median of the data
array. then compute the average and the rms deviation
of the data that remain.
*/
long int median(long int[], int); /*compute median*/
long int mm; /*computed median*/
long int ml,mu; /*lower and upper acceptance bounds*/
int jm; /*number of data in mean */
int j;
long int mt; /*temporary */
double xx; /*another temp*/
int iflag[1024]; /*static array for speed*/
double sw,swmax; /*largest value of outlier*/
int jrej; /*index of rejected value*/
for(j=0; j<indx; j++) iflag[j] = 0; /*everybody ok to start*/
mm=median(ivals,indx); /*compute median of data */
ml=mm - 10000;
mu=mm + 10000;
/*
first compute mean of all data that differ from
median by less than rejection threshold specified
by ml and mu above.
*/
loop:
*sval = 0;
jm=0;
for(j=0; j<indx; j++)
{
if(iflag[j] == 1) continue; /* value already rejected*/
mt=ivals[j];
if( (mt >= ml) && (mt <= mu) )
{
jm++;
*sval += (double) mt;
}
else iflag[j] = 1; /*this value bad*/
}
/*
if no points satisfy the criterion, then return
9999 for rms.
*/
if(jm == 0)
{
*rms=9999;
return;
}
*sval /= jm;
/*
now compute rms deviation of these folks
if there is only one point left then rms is
simply set to the absolute value of the mean
*/
if(jm == 1)
{
*rms=fabs(*sval);
return;
}
*rms=0;
jm=0;
for(j=0; j<indx; j++)
{
if(iflag[j] == 1) continue; /*this value is already bad*/
mt=ivals[j];
xx= (double) mt - *sval;
*rms += xx*xx;
jm++;
}
*rms= sqrt(*rms/(double) (jm-1));
/*
if any value differs from the mean by more than
3 sigma, then reject the largest value and go around
again.
*/
swmax = 0;
for(j=0; j<indx; j++)
{
if(iflag[j] == 1) continue; /*value already bad*/
sw=fabs( (ivals[j] - *sval)/ *rms);
if(sw > swmax)
{
swmax=sw; /*this one is larger*/
jrej=j;
}
}
if(swmax > 3)
{
iflag[jrej] = 1; /*reject the largest*/
goto loop; /*and go around again. */
}
}
int oops()
{
timeout++;
return (1);
}