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 >
C/C++ Source or Header  |  2003-12-09  |  14KB  |  423 lines
  1. #include <stdio.h>
  2. #include <sys/types.h>
  3. #include <fcntl.h>
  4. #include <sys/stat.h>
  5. #include <sys/socket.h>
  6. #include <sys/time.h>
  7. #include <netinet/in.h>
  8. #include <arpa/inet.h>
  9. #include <string.h>
  10. #include <ctype.h>
  11. #include <signal.h>
  12. #include <math.h>
  13. #ifndef LINUX
  14. #include <stropts.h>
  15. #endif
  16. #include "autolock.h"
  17. int timeout;        /*set when SIGALRM happens*/
  18. int getdif(char *cp,struct autolockst *lock,int *lsr)
  19. {
  20. /*
  21.     this subroutine estimates the difference between the
  22.     local clock and a server.  The comparison is made using
  23.     NTP-format messages.
  24.     the integer lsr is set to the leap second flag received
  25.     from the remote host.
  26.     the possible values are:
  27.         0 -- server ok, no leap second
  28.         1 -- server ok, insert second
  29.         2 -- server ok, drop second
  30.         3 -- server not ok.
  31.  
  32.     if the server returns a value of 3, this subroutine
  33.     returns immediately with a status of -4.  the other
  34.     3 possibilities are simply passed back to the caller.
  35.  
  36.     the subroutine attempts to estimate the uncertainty in the
  37.     NTP time exchange and tries to adjust its performance so that
  38.     the measured uncertainty is consistent with the desired performance
  39.     as specified in the global parameters described below.
  40.  
  41.         This software was developed with US Government support
  42.         and it may not be sold, restricted or licensed.  You
  43.         may duplicate this program provided that this notice
  44.         remains in all of the copies, and you may give it to
  45.         others provided they understand and agree to this
  46.         condition.
  47.  
  48.         This program and the time protocol it uses are under
  49.         development and the implementation may change without
  50.         notice.
  51.  
  52.         For questions or additional information, contact:
  53.  
  54.         Judah Levine
  55.     JILA
  56.     Campus Box 440
  57.     University of Colorado
  58.         Boulder, Colorado 80309
  59.         (303) 492 7785
  60.         jlevine@utcnist.colorado.edu
  61.  
  62. */
  63. #include "ntp-packet.h"
  64. struct ntppkt msg;        /* the message we send out*/
  65. struct ntppkt ans;        /* the response*/
  66. struct timeval tvv;        /* get local system time*/
  67. int address;            /* holds ip address */
  68. int pserv = 123;        /* ntp port number on server */
  69. unsigned int ts,fts;        /* transmission time of outbound packet*/
  70. unsigned int tr,ftr;        /* time and which response was received*/
  71. unsigned int tx,ftx;        /* remote transmission and received times*/
  72. unsigned int trr,ftrr;
  73. long int dly;            /* one-way delay in microseconds*/
  74. long int tdiff[NCALMAX];    /* time difference in microseconds*/
  75. /*
  76.     ntpdif computes the difference between two times
  77.     expressed in NTP format
  78. */
  79. long int ntpdif (unsigned int, unsigned int, unsigned int, unsigned int);
  80. long int median (long int[], int);
  81. void smooth(long int[], int, double *, double *);
  82. int oops();            /*receives time-out alarms*/
  83. extern int ncal;        /*number of pings to a server*/
  84. double val,rms;
  85. struct sockaddr_in sin;         /* socket address structure */
  86. int s;                          /* socket number */
  87. int length;                     /* size of message */
  88. int i,j;
  89. int kval;            /*number of packets actually received.*/
  90. int leapind;            /*received leap indicator*/
  91. int mjd;            /*cycle time converted to MJD */
  92. int secoday;            /*the time in seconds after 0000*/
  93. int msosec;            /*the fraction of a second in milliseconds*/
  94.     signal(SIGALRM,(void *) &oops);    /*set routine to handle time-outs*/
  95. /*
  96.         convert address to internal format
  97.     then store address and port number in socket structure
  98. */
  99.         if( (address=inet_addr(cp) ) == -1) return(-1);    /*bad address*/
  100.         bzero( (char *)&sin,sizeof(sin));
  101.         sin.sin_addr.s_addr=address;
  102.         sin.sin_family=AF_INET;
  103.         sin.sin_port=htons(pserv);
  104. /*
  105.         create the socket then connect to the server.
  106.     return error status of -2 if socket creation fails
  107.     and error status of -3 if connetion fails.
  108. */
  109.         if( (s=socket(AF_INET,SOCK_DGRAM,0) ) < 0) return(-2);
  110.         if(connect(s, (struct sockaddr *) &sin, sizeof(sin) ) < 0) return(-3);
  111.     kval=0;            /*index showing number of measurements*/
  112.     for(j=0; j<ncal; j++)     /*ping 'em again and again*/
  113.        {
  114. /*
  115.        build an NTP packet and send it off.
  116.        
  117. */
  118.        msg.li_vn_mode=033;        /* 00 011 011 = no leap, v3, client*/
  119.        msg.stratum= 0;        /* 0 means unspecified*/
  120.        msg.ppoll= 8;        /* poll interval = 256s*/
  121.        msg.precision= 0;        /* 0= 1s, no precision at all here*/
  122.        msg.rootdelay=0;        /* don't really know or care*/
  123.        msg.rootdispersion=0;    /* ditto*/
  124.        msg.refid = 0;        /* we don't have a reference*/
  125. /*
  126.        system request returns time as seconds and microseconds since
  127.        1970.0.  Must convert this to NTP time format, which uses
  128.        1900.0 as the origin and fractions of a second instead of
  129.        microseconds.  Thus must add seconds 1900.0 -> 1970.0 to
  130.        system time in seconds since 1970.
  131.        The fractional part divides 1 second into 2**32 parts, so
  132.        that one least count is 232.8306 ps.  Thus the fraction is
  133.        system time in microseconds * 10**6/232.8306 = 
  134.        system time * 4294.9673
  135. */
  136.        gettimeofday(&tvv,(struct timezone *)NULL);
  137.        ts= 2208988800u + (unsigned int) tvv.tv_sec;
  138.        fts= (unsigned int) ( 4294.9673*(double)tvv.tv_usec);
  139.        msg.reftime.sys[0]=htonl(ts-1);
  140.        msg.reftime.sys[1]=htonl(fts);
  141.        msg.org.sys[0]= htonl(ts);
  142.        msg.org.sys[1]= htonl(fts);
  143.        alarm(5);        /*5 seconds -- should be long enough*/
  144.        timeout=0;
  145.        length=write(s,&msg,sizeof(msg));
  146. /*
  147.        read response and get local time as soon as read completes.
  148. */
  149.        do
  150.          {
  151.            length=recvfrom(s,&ans,sizeof(ans),0,
  152.         (struct sockaddr *) NULL,(size_t *) NULL);
  153.          gettimeofday(&tvv,(struct timezone *)NULL);
  154.        } while ( (length < sizeof(ans)) && (timeout == 0) );
  155.        alarm(0);        /*turn off the timer*/
  156.        if(timeout != 0)  continue;    /* try again if time out*/
  157.        tr= 2208988800u + (unsigned int) tvv.tv_sec;
  158.        ftr= (unsigned int) (4294.9673*(double)tvv.tv_usec);
  159.        leapind=ans.li_vn_mode;    /*copy to an integer*/
  160.        leapind >>= 6;
  161.        leapind &= 3;        /*isolate leap indicator bits*/
  162.        if(leapind == 3) return(-4);    /*server is not healthy*/
  163.        *lsr= leapind;        /*send leap second bits to caller*/
  164.        trr=ntohl(ans.rec.sys[0]);
  165.        ftrr=ntohl(ans.rec.sys[1]);
  166.        tx=ntohl(ans.xmt.sys[0]);
  167.        ftx=ntohl(ans.xmt.sys[1]);
  168. /*
  169.        dly is total elapsed time as measured on our local clock
  170.        - time spent in remote server, which is the transmit
  171.        time of the remote server - the received time.
  172. */
  173.        dly= ntpdif(tr,ftr,ts,fts) - ntpdif(tx,ftx,trr,ftrr);
  174.        tdiff[kval++]=ntpdif(ts,fts,trr,ftrr) + dly/2; /*local- remote*/
  175.     }    /* end of loop over j */
  176.     close(s);        /*done with the network connection*/
  177.     if(kval == 0) return(-5);    /* nobody home ? */
  178.     if(kval == 1)
  179.        {
  180.        val=tdiff[0];    /*only one value*/
  181.        rms= 999;        /*so rms has no meaning*/
  182.     }
  183.     else if(kval == 2)
  184.        {
  185.        val= (tdiff[0]+tdiff[1])/2;
  186.        rms= fabs(tdiff[0]-tdiff[1])/1.414;
  187.     }
  188.     else         /*more than 2 responses*/
  189.        {
  190.        smooth(tdiff,kval,&val,&rms);
  191.     }
  192.     mjd=ts/86400 + 15020;        /*compute MJD*/
  193.     secoday= ts%86400;
  194.     msosec= ((double) fts)/4294968.1;
  195.     lock->day=mjd;
  196.     lock->fday= 1000*secoday + msosec;
  197.     lock->x= 1e-6*val;  /*convert from usec to sec*/
  198.     lock->rms= 1e-6*rms; /*ditto*/
  199.     lock->nval=kval;
  200.     return(1);
  201. }
  202. long int ntpdif(unsigned int t1,unsigned int ft1,
  203.           unsigned int t2,unsigned int ft2)
  204. {
  205. /*
  206.     this subroutine computes the difference between the first
  207.     and second pairs of NTP-format times and returns the
  208.     difference in microseconds.
  209.     The input parameter pairs are seconds and fraction of
  210.     a second, respectively, in standard NTP format.
  211.     Since the fraction divides 1 second into 2**32 parts,
  212.     a least-count of the fraction is 232.8306 ps.  The
  213.     value in microseconds is therefore (value)*232.83/1000000
  214.  
  215.     note that NTP-format numbers are unsigned so that 
  216.     subtractions will produce strange results if performed
  217.     directly.
  218. */
  219. long int temp;
  220. long int ta;
  221. double tb;
  222.     ta= (long int) t1 - (long int) t2;
  223.     tb= (double) ( (long int) ft1 - (long int) ft2 );
  224.     temp= 1000000*ta + 
  225.         (long int)(tb/4294.9673);
  226.     return (temp);
  227. }
  228. long int median(long int ivals[],int indx)
  229. {
  230. /*
  231.     this subroutine computes the median of the input
  232.     array ivals.  There are indx values in the array.
  233.     the median is returned and the input data are
  234.     not altered.
  235.  
  236.     the median is that number for which there are imid
  237.     entries both bigger and smaller than it.  there is a
  238.     single index value that will satisfy this requirement
  239.     if the number of entries is odd and two adjacent ones
  240.     will satisfy it if the number is even. the average of
  241.     these two adjacent values are returned in this case.
  242.  
  243.     if there are a significant number of repeated values,
  244.     then no entry may satisfy this requirement.  this case
  245.     can be dealt with approximately by finding an entry that
  246.     approximately satisfies the test, but the median may not
  247.     do the "right" thing in this case.  this subroutine is
  248.     complicated by the facts that there are likely to be a
  249.     significant number of duplicates, that the input array
  250.     must be left alone and that the subroutine is in the
  251.     main loop and must be as fast as possible.
  252. */
  253. int j;
  254. int imid;        /*index of mid point value */
  255. char iflag[NCALMAX];    /*static array for speed */
  256. int imin;        /*index of minimum value*/
  257. long int mmin;        /*current minimum value*/
  258. int ifound =0;        /*index number of min just found */
  259. int iodd;        /*1 if indx is even, 0 if odd */
  260. long int m1;        /*lower median if size is even*/
  261.     iodd= indx & 1;
  262.     imid=indx/2;     /*index of mid-point value for odd size*/
  263.     for(j=0; j<indx; j++) iflag[j]=0;
  264. /*
  265.     loop over all elements and find the minimum one.
  266.     set the flag to show that it has been found and
  267.     then loop again to find the minimum one in the
  268.     ones that are left.  This process will always
  269.     converge even if all of the values are the same,
  270.     but it may not do the "right" thing if there are
  271.     only a few distinct values that appear many times.
  272.  
  273.     first find first value that has not yet been tagged
  274.     as a minimum in some round.
  275. */
  276. loop:
  277.     for(imin=0; imin<indx; imin++)
  278.        if(iflag[imin] == 0) 
  279.           {
  280.           mmin=ivals[imin];
  281.           break;
  282.        }
  283. /*
  284.     now compare this value to all other untagged
  285.     ones to find the minimum of the lot.  the surviving 
  286.     minimum would go in the next place in an output array 
  287.     that was ordered by size.  its index in that array would 
  288.     be ifound.  if there are an odd number of points and 
  289.     ifound=imid then we have found the median in this 
  290.     minimum and we are finished.  if there are an even
  291.     number of points, then we must average the values
  292.     found before and after imid to compute the median.
  293.     if this is the one before then just save it.  if
  294.     this is the one after then compute average and
  295.     return.
  296.     if ifound is too small in either case, then
  297.     we have to go around again. tag the surviving minimum 
  298.     and go around to find the minimum of the remaining
  299.     untagged folks.
  300. */
  301.     for(j=imin+1; j<indx; j++) 
  302.        if( (iflag[j] == 0) && (ivals[j] < mmin) )
  303.           {
  304.           mmin=ivals[j];
  305.           imin=j;
  306.        }
  307.     if(iodd == 1) 
  308.        {
  309.        if (ifound == imid) return (mmin);
  310.     }
  311.     else            /*size is even*/
  312.        {
  313.        if(ifound == (imid - 1) ) m1=mmin;
  314.        if(ifound == (imid + 1) ) return( ( (m1+mmin) )/2 );
  315.     }
  316.     ifound++;
  317.     iflag[imin]=1;
  318.     goto loop;
  319. }
  320. void smooth(long int ivals[],int indx,double *sval,double *rms)
  321. {
  322. /*
  323.     this subroutine smooths the input array ivals which has
  324.     indx members.  the smoothed value is returned in sval
  325.     and the rms deviation is returned in rms.
  326.  
  327.     the smoothing consists of first throwing away all values that
  328.     are more than 10 ms = 10000 us away from the median of the data
  329.     array.  then compute the average and the rms deviation
  330.     of the data that remain.
  331. */
  332. long int median(long int[], int);        /*compute median*/
  333. long int mm;            /*computed median*/
  334. long int ml,mu;            /*lower and upper acceptance bounds*/
  335. int jm;                /*number of data in mean */
  336. int j;
  337. long int mt;            /*temporary */
  338. double xx;            /*another temp*/
  339. int iflag[1024];        /*static array for speed*/
  340. double sw,swmax;            /*largest value of outlier*/
  341. int jrej;            /*index of rejected value*/
  342.     for(j=0; j<indx; j++) iflag[j] = 0;    /*everybody ok to start*/
  343.     mm=median(ivals,indx);    /*compute median of data */
  344.     ml=mm - 10000;
  345.     mu=mm + 10000;
  346. /*
  347.     first compute mean of all data that differ from
  348.     median by less than rejection threshold specified
  349.     by ml and mu above.
  350. */
  351. loop:
  352.     *sval = 0;
  353.     jm=0;
  354.     for(j=0; j<indx; j++)
  355.        {
  356.        if(iflag[j] == 1) continue;    /* value already rejected*/
  357.        mt=ivals[j];
  358.        if( (mt >= ml) && (mt <= mu) )
  359.           {
  360.           jm++;
  361.           *sval += (double) mt;
  362.        }
  363.        else iflag[j] = 1;        /*this value bad*/
  364.     }
  365. /*
  366.     if no points satisfy the criterion, then return
  367.     9999 for rms.
  368. */
  369.     if(jm == 0)
  370.        {
  371.        *rms=9999;
  372.        return;
  373.     }
  374.     *sval /= jm;
  375. /*
  376.     now compute rms deviation of these folks
  377.     if there is only one point left then rms is
  378.     simply set to the absolute value of the mean
  379. */
  380.     if(jm == 1) 
  381.        {
  382.        *rms=fabs(*sval);
  383.        return;
  384.     }
  385.     *rms=0;
  386.     jm=0;
  387.     for(j=0; j<indx; j++)
  388.        {
  389.        if(iflag[j] == 1) continue;    /*this value is already bad*/
  390.        mt=ivals[j];
  391.        xx= (double) mt - *sval;
  392.        *rms += xx*xx;
  393.        jm++;
  394.     }
  395.     *rms= sqrt(*rms/(double) (jm-1));
  396. /*
  397.     if any value differs from the mean by more than 
  398.     3 sigma, then reject the largest value and go around 
  399.     again.
  400. */
  401.     swmax = 0;
  402.     for(j=0; j<indx; j++)
  403.        {
  404.        if(iflag[j] == 1) continue;    /*value already bad*/
  405.        sw=fabs( (ivals[j] - *sval)/ *rms);
  406.        if(sw > swmax) 
  407.           {
  408.           swmax=sw;        /*this one is larger*/
  409.           jrej=j;
  410.        }
  411.     }
  412.     if(swmax > 3)
  413.        {
  414.        iflag[jrej] = 1;    /*reject the largest*/
  415.        goto loop;        /*and go around again.  */
  416.     }
  417. }
  418. int oops()
  419. {
  420.     timeout++;
  421.     return (1);
  422. }
  423.