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 / autolock.c < prev    next >
C/C++ Source or Header  |  2003-12-09  |  61KB  |  1,519 lines
  1. #include <stdio.h>
  2. #include <math.h>
  3. #include <signal.h>
  4. #include <stdlib.h>
  5. #include <sys/types.h>
  6. #include <sys/stat.h>
  7. #include <fcntl.h>
  8. #include <sys/ioctl.h>
  9. #include <sys/param.h>
  10. #include <sys/wait.h>
  11. #include <sys/time.h>
  12. #include <sys/timex.h>
  13. #include <string.h>
  14. #include <unistd.h>
  15. #include "sizint.h"
  16. #include "autolock.h"
  17. /*
  18.     the following global variables define the operation of
  19.     the server. they are initialized from the parameter file
  20.     during the cold-start procedure. If they are not defined
  21.     in the parameter file then the values specified here are
  22.     used by default.
  23. */
  24. int sigma_desired;    /*the desired time performance in ms. */
  25. int ncalmax = NCALMAX;    /*the maximum number of ntp messages each cycle*/
  26. int ncalmin = 3;    /*the minimum number of ntp messages each cycle*/
  27. int intvlmin = 1000;    /*the minimum interval between cycles*/
  28. int intvlmax = 80000;    /*the maximum interval between cycles*/
  29. int opmode;        /*the initial mode of the program */
  30. int ncal= 5;        /*initial number of ntp messages on each cycle*/
  31. int intvl= 2000;    /*initial interval between cycles*/
  32. int evalint= 259200;    /*default is 3 days in seconds*/
  33. int tconf= 12000;    /*time constant (in seconds) for ybar update*/
  34. int ttot= 0;        /*total running time in seconds*/
  35. #ifdef KERNEL        /*if jl version of kernel routines available*/
  36. int kernel= 1;        /*default is to use kernel routines if possible*/
  37. struct timex tmx;    /*structure for linkage to kernel routines*/
  38. #endif    /*KERNEL*/
  39. #ifdef LINUX
  40. struct timex tmx;    /*used for adjstimex routine in Linux*/
  41. #endif    /*LINUX*/
  42. /*
  43.     opmode = 1    means characterize the network link and
  44.             compare it to the noise in the local clock
  45.             itself. exit after running for evalint seconds.
  46.  
  47.     opmode =2    means change to lock mode after the evaluation
  48.             is completed or else do a warm re-start and
  49.             switch to lock mode immediately using the parameters
  50.             from a previous evaluation.
  51. */
  52. main(int argc,char *argv[])
  53. {
  54. /*
  55.     function prototypes
  56. */
  57. int getdif(char [],struct autolockst *,int *);
  58.                 /*gets time difference local-server*/
  59. int sw(int, char* [],char *, long int *);  /*parses command line*/
  60. void getautolock(FILE **, struct autolockst *);    /*reads data from file*/
  61. void putautolock(FILE **, struct autolockst[]);    /*writes data to file*/
  62. void putmsg(char *);        /*writes message to log file*/
  63. void pushautolock(struct autolockst[],int,int);    /*pushes down data stack*/
  64. int sett(double);        /*steps time in one step using settimeofday*/
  65. int adjt(double);        /*adjusts time by slewing clock*/
  66. double adev(double [],int ,int); /*computes Allan deviation */
  67.  
  68. char serv[NUMSRV][20];     /*addresses of primary and alternate servers*/
  69. int i,j;
  70. int j1,j2;
  71. int nhr,nmin;            /*time of next cycle*/
  72. long int v;            /*values returned from command line switch*/
  73. char s;                /*command-line switch letter */
  74. int debug = 0;            /*flag for debug mode, default is off*/
  75. extern int errno;        /*error number from the system*/
  76. FILE *pop;            /*used to read from parameter file*/
  77. FILE *lop;            /*used to read from state file*/
  78. struct autolockst autolock[NUMSTATES];        /*push-down stack of states*/
  79. int numsrv = 0;            /*number of servers defined, up to NUMSRV*/
  80. char errmsg[100];        /*holds error message text*/
  81. char keywoid[30];        /*holds keyword from parameter file*/
  82. char keyval[30];        /*holds value from parameter file*/
  83. double fract;            /*time constant parameter*/
  84. int ntry;            /*number of consecutive calls to getdif*/
  85. int ntry2;            /*number of consecutive jitter failures*/
  86. int ntry3;            /*number of consecutive prediction failures*/
  87. int fstsrv;            /*index of server to try first*/
  88. int altsrv;            /*number of alternate servers we have tried*/
  89. double alterr[NUMSRV];        /*time diff from alternate servers*/
  90. double altrms[NUMSRV];        /*rms of time differences from alternates*/
  91. int filefull= NUMSTATES;    /*counter to see when file is getting full*/
  92. char filenam[80];        /*used in the file rename department*/
  93. int pid;                        /*return value from  wait call*/
  94. union wait *status;             /*return status from wait call*/
  95. int eval_elap= EVALPERIOD;    /*interval since last performance evaluation*/
  96. int tloop2;            /*number of calibration cycles in mode 2*/
  97. int tloop3;            /*number of calibration cycles in mode 3*/
  98. int tloop4;            /*number of cycles in mode 4*/
  99. char *cp;            /*pointer to command line parameter*/
  100. double tcor;            /*correction in integer seconds for sett*/
  101. double terr;            /*absolute value of x*/
  102. double tsig;            /*standard deviation of multiple servers*/
  103. double ddum;            /*temporary double*/
  104. int linno;            /*line number of configuration file*/
  105. double y[NUMSTATES];        /*frequency estimates for Avar calculation*/
  106. int numby;            /*pointer to next free location in y array*/
  107. double dtemp;            /*temporary double*/
  108. int itemp;            /*temporary integer*/
  109. double avar[10];        /*Allan deviation as function of lag*/
  110. int nvar;            /*index into avar array*/
  111. int jwait;            /*time between periodic corrections*/
  112. int jloop;            /*number of periodic corrections per cycle*/
  113. /*
  114.     lsl and lsr give the local and remote notions of the leap
  115.     second flag (respectively).  lsr is normally set by getdif
  116.     based on a message from a server; this change is normally
  117.     propagated into lsl.  What happens next depends on whether
  118.     or not the kernel knows how to do a leap second by itself.
  119.     see the comments below.
  120. */
  121. int lsl= 0;            /*local notion of the leap second flag*/
  122. int lsr=0;            /*server notion of the leap second flag*/
  123.     while( sw(argc,argv,&s,&v) != 0)  /*decode command line switches*/
  124.            {
  125.            switch (s)
  126.                 {
  127.                 case 'd':        /*turn on debug mode*/
  128.            debug=1;
  129.            fprintf(stderr,"\n Entering debug mode.");
  130.            break;
  131.         default:
  132.            fprintf(stderr,"\n Switch %c not recognized.\n",s);
  133.            exit(-1);
  134.        }
  135.        argc--;    /*one less parameter*/
  136.        argv++;    /*advance to next parameter*/
  137.     }
  138. /*
  139.     if debug is not defined, then we will run in normal mode.
  140.         first convert ourselves into a daemon.
  141.         We enter the phone booth as mild-mannered Clark Kent,
  142.         we make a copy of ourselves, the original disappears,
  143.         the child emerges and then ....
  144. */
  145.     if(debug == 1) goto getparam;    /*remain a normal program in debug */
  146.         signal(SIGTTOU,SIG_IGN);        /*the keyboard can't stop us*/
  147.         signal(SIGTTIN,SIG_IGN);
  148.         signal(SIGTSTP,SIG_IGN);
  149.         if( (i=fork()) < 0)  abort();   /*fork of child failed*/
  150.         if(i > 0)
  151.            {
  152.            printf("\n pid of process is %d",i);
  153.            exit(0);             /*the parent exits */
  154.         }
  155. /*
  156.         make ourselves into a unique process group so that we
  157.         can't re-acquire a control terminal that might be able
  158.         to stop us.
  159.         Then close our control terminal.
  160.         Nothing can stop us now ...
  161. */
  162. #ifdef LINUX
  163.     if( setpgid(0,getpid() ) == -1) abort();
  164. #else    /*LINUX*/
  165.     if( setpgrp(0,getpid() ) == -1) abort();
  166. #endif    /*LINUX*/
  167.         if( (i= open("/dev/tty",O_RDWR) ) >= 0)
  168.            {
  169.            ioctl(i,TIOCNOTTY,(char *) NULL);    /*bye bye terminal*/
  170.            close(i);
  171.         }
  172.         for(i=0; i< NOFILE; i++) close(i);      /*close all files*/
  173.         errno=0;                        /*reset errors, if any*/
  174. /*
  175.     now get the default parameters for the process
  176. */
  177. getparam:    
  178. /*
  179.     change to default directory, open parameter file and
  180.     get parameters.
  181. */
  182.     if(chdir(BINDIR) < 0)
  183.        {
  184.        sprintf(errmsg,"Oops - change directory to %s failed.\n",BINDIR);
  185.        putmsg(errmsg);
  186.        exit(-1);
  187.     }
  188.     if( (pop=fopen(PNAME,"r")) == NULL) /*open read only*/
  189.       {
  190.       sprintf(errmsg,"Oops -- can not open parameter file %s\n",PNAME);
  191.       putmsg(errmsg);
  192.       exit(-1);
  193.     }
  194. /*
  195.     the parameter file is a free-form file of the form:
  196.  
  197.     keyword value <newline>
  198.  
  199.     if the keyword at the start of a line is # then the line is
  200.     treated as a comment.  if there is a keyword but it is not 
  201.     recognized then the line is treated as a comment and a
  202.     comment is written to the log file.
  203.     the keyword comparison is done in lower case, and upper
  204.     case keywords are converted before comparing them.
  205. */
  206.     linno=0;    /*reset line counter*/
  207. readparam:
  208.     linno++;    /*increment line counter*/
  209.     if(fscanf(pop," %s",keywoid) == EOF) goto finparam;
  210.     if(keywoid[0] == 043)    goto nxtlin;    /* # starts a comment line*/
  211.     if(fscanf(pop," %s",keyval) != 1) goto fmterr;
  212.     for(i=0; keywoid[i] != '\0'; i++) 
  213.         keywoid[i]=tolower( (int) keywoid[i]);
  214.     if(strncmp(keywoid,"sigma",5) == 0) 
  215.        {
  216.        if(sscanf(keyval," %d",&sigma_desired) != 1) goto fmterr;
  217.        goto nxtlin;
  218.     }
  219.     if(strncmp(keywoid,"nmax",4) == 0)
  220.        {
  221.        if(sscanf(keyval," %d",&ncalmax) != 1) goto fmterr;
  222.        goto nxtlin;
  223.     }
  224.     if(strncmp(keywoid,"nmin",4) == 0)
  225.        {
  226.        if(sscanf(keyval," %d",&ncalmin) != 1) goto fmterr;
  227.        goto nxtlin;
  228.     }
  229.     if(strncmp(keywoid,"tmin",4) == 0)
  230.        {
  231.        if(sscanf(keyval," %d",&intvlmin) != 1) goto fmterr;
  232.        goto nxtlin;
  233.     }
  234.     if(strncmp(keywoid,"tmax",4) == 0)
  235.        {
  236.        if(sscanf(keyval," %d",&intvlmax) != 1) goto fmterr;
  237.        goto nxtlin;
  238.     }
  239.     if(strncmp(keywoid,"tzero",5) == 0)
  240.        {
  241.        if(sscanf(keyval," %d",&intvl) != 1) goto fmterr;
  242.        goto nxtlin;
  243.     }
  244.     if(strncmp(keywoid,"mode",4) == 0)
  245.        {
  246.        if(sscanf(keyval," %d",&opmode) != 1) goto fmterr;
  247.        goto nxtlin;
  248.     }
  249.     if(strncmp(keywoid,"server",6) == 0)
  250.        {
  251.        if(numsrv >= NUMSRV) goto nxtlin;    /*server list is full*/
  252.        if(sscanf(keyval," %s",&serv[numsrv][0]) != 1) goto fmterr;
  253.        numsrv++;    /*total number of servers we know about*/
  254.        goto nxtlin;
  255.     }
  256.     if(strncmp(keywoid,"evalt",5) == 0)
  257.        {
  258.        if(sscanf(keyval," %d",&evalint) != 1) goto fmterr;
  259.        goto nxtlin;
  260.     }
  261.     if(strncmp(keywoid,"tconf",5) == 0)
  262.        {
  263.        if(sscanf(keyval," %d",&tconf) != 1) goto fmterr;
  264.        goto nxtlin;
  265.     }
  266. #ifdef KERNEL
  267. /*
  268.     the keywoid "kernel" will only be recognized if
  269.     KERNEL is defined.  otherwise, this code will not be
  270.     present and it will be treated as a format error.
  271. */
  272.     if(strncmp(keywoid,"kernel",6) == 0)
  273.        {
  274.        if(sscanf(keyval," %d",&kernel) != 1) goto fmterr;
  275.        goto nxtlin;
  276.     }
  277. #endif    /*KERNEL*/
  278. fmterr:            /*format error or did not recognize keyword*/
  279.     sprintf(errmsg,    
  280.         "Format error reading line %d starting with %s",
  281.                 linno,keywoid);
  282.     putmsg(errmsg);
  283. nxtlin:            /*skip remainder of line, go back and read next one*/
  284.     while( (s=getc(pop)) != '\n') if(s == EOF) goto finparam;
  285.     goto readparam;
  286. finparam:        /*done with the parameter file*/
  287.     fclose(pop);        /*finished reading parameters*/
  288. /*
  289.     now open state file and read the most recent states, which
  290.     are the first ones in the file. The read ends either when
  291.     NUMSTATES have been read in or when a state with a day number
  292.     of 0 is found.  This can also happen when the physical EOF
  293.     is reached, but normally there should be an explicit record
  294.     with a day number of 0.  If fewer than NUMSTATES are read
  295.     from the file, then the remaining entries are set to 0
  296.     by the following while loop.
  297.     note that the most recent state with index 0 is not read
  298.     by  this process.  It will be changed in response to the
  299.     getdif request.
  300. */
  301.     if( (lop=fopen(SNAME,"r")) == NULL) /*open read only*/
  302.        {
  303.        sprintf(errmsg,"Oops -- error opening state file %s.\n",SNAME);
  304.        putmsg(errmsg);
  305.        exit(-1);
  306.     }
  307.     for(i=1; i < NUMSTATES; i++)
  308.        {
  309.        getautolock(&lop,&autolock[i]);
  310.        if(autolock[i].day == 0) break;    /*logical end of data*/
  311.     }
  312.     fclose(lop);
  313.     while(i < NUMSTATES) autolock[i++].day = 0; /*set remaining ones to 0*/
  314. /*
  315.     if the most recent state read from the file does not have an MJD
  316.     value of 0, then transfer the static parameters from that state
  317.     to the current state, which is number 0
  318.     there is no point in transferring those parameters that will
  319.     be changed during the call to getdif below.
  320.  
  321.     note that the time constant for the frequency filter is
  322.     copied from the existing state vector, and that this value
  323.     may not be the same as the value read from the initial parameter
  324.     file or the default value in the program definition. This version
  325.     of the software may adjust fcon based on the results of an Allan
  326.     variance calculation -- see the call to adev below.
  327.  
  328.     what happens next depends on opmode.  if opmode = 1 then we
  329.     are just trying to evaluate the local clock, and we are going
  330.     to start over again as if this was a cold start even though
  331.     we might have read previous data from the file.
  332.     if opmode == 2 then we are going to proceed to a lock state.
  333.     advance to mode 4 if the previous mode is 3 -- i.e., lock the
  334.     clock if we have read a previous evaluation.  otherwise go
  335.     back and do a complete cold-start evaluation.
  336. */
  337.     if(autolock[1].day != 0)
  338.        {
  339.        autolock[0].mode=autolock[1].mode;
  340.        autolock[0].rms_avg=autolock[1].rms_avg;
  341.        autolock[0].d= autolock[1].d;    /*normally set to 0*/
  342.        intvl=autolock[0].waitn= autolock[1].waitn;
  343.        autolock[0].ybar= autolock[1].ybar;
  344.        autolock[0].fcon= autolock[1].fcon;
  345.        autolock[0].sumerr= autolock[1].sumerr;
  346.        if( (opmode == 1) ||     /*either we are just evaluating, or*/
  347.            (autolock[0].mode < 3) ) /*lock coming but eval not finished*/
  348.           {
  349.           autolock[0].mode= 0;    /*then go back and do a cold start*/
  350.           autolock[0].waitn= intvl;
  351.        }        /*end of opmode == 1 or mode < 3*/
  352.        else 
  353.           {
  354.           autolock[0].mode= 4;    /*opmode ==2 and mode >= 3*/
  355.           tloop4= 0;        /*initialize loop counter for mode 4*/
  356.        }    /*end of opmode != 1 and mode >= 3*/
  357.     }        /*end of read a day number != 0*/
  358.     else        /*day number == 0 means absolute cold start*/
  359.        {
  360.        autolock[0].waitn=intvl;
  361.        autolock[0].mode= 0;
  362.     }        /*end of absolute cold start*/
  363. #ifdef KERNEL
  364. /*
  365.     if the kernel routines are available and if they have been
  366.     selected, then reset all of the kernel parameters to 0,
  367.         and set the health to bad.
  368.         note that MOD_LOCKCLOCK must always be specified in the jl
  369.     version or else a request to change anything will be ignored.
  370.         the expected return status is TIME_ERROR because we have
  371.         just set the clock to unsynchronized.  if we don't get this,
  372.         then there is something wrong with the linkage to the
  373.         kernel.
  374. */
  375.     if(kernel == 1)        /*kernel routines selected*/
  376.          {
  377.            tmx.modes= MOD_LOCKCLOCK | MOD_STATUS | MOD_OFFSET | MOD_FREQUENCY |
  378.                         MOD_MAXERROR | MOD_ESTERROR | MOD_TIMECONST |
  379.                         MOD_DST | MOD_HEALTH | MOD_LS;
  380.            tmx.status=STA_UNSYNC;               /*we start out un synched*/
  381.            tmx.offset=tmx.freq=tmx.maxerror=tmx.esterror=tmx.constant=0;
  382.            tmx.dst=tmx.ls=0;
  383.            tmx.health=3;    /*and un healthy*/
  384.            j=syscall(SYS_ntp_adjtime,&tmx);     /*tell it to the kernel */
  385.            if(j != TIME_ERROR)
  386.               {
  387.               sprintf(errmsg,
  388.                  "103-return status from syscall is %d, expected=%d",
  389.                         j,TIME_ERROR);
  390.               putmsg(errmsg);
  391.            }        /*end of return status not as expected*/
  392.         }        /*end of kernel routines selected*/
  393. #endif         /*KERNEL*/
  394. #ifdef LINUX
  395.     tmx.modes= ADJ_OFFSET | ADJ_FREQUENCY | ADJ_MAXERROR |
  396.         ADJ_ESTERROR | ADJ_STATUS | ADJ_TIMECONST;
  397.     tmx.offset=tmx.freq=tmx.maxerror=tmx.esterror=0;
  398.     tmx.constant=tmx.precision=tmx.tolerance=0;
  399.     tmx.status= TIME_BAD;
  400.     j=adjtimex(&tmx);
  401.     if(j != TIME_BAD)
  402.        {
  403.        sprintf(errmsg,
  404.         "103-return status from adjtimex is %d, expected=%d",
  405.         j,TIME_BAD);
  406.        putmsg(errmsg);
  407.     }
  408. #endif    /*LINUX*/
  409. /*
  410.     now we are ready for serious stuff.
  411.     infinite loop of main code begins here.
  412.  
  413.     each time we make another estimate of the frequency,
  414.     we save the value in array y, and numby is the index
  415.     into this array.  the index is reset if the value
  416.     of waitn changes so that the array contains only
  417.     values at a constant value of waitn.
  418.     these data are then used to calculate an estimate of
  419.     Avar once a day or so.  this estimate is only approximate
  420.     since it will not be robust if there have been a lot
  421.     of errors during the day. There is no good way of dealing
  422.     with errors since at the very least they introduce gaps
  423.     into the array. When the clock is free-running during the
  424.     evaluation phase, the avar estimate is used to estimate a
  425.     value for fcon, the time-constant for the frequency estimator.
  426.     Since this is the only purpose in this version, there is not 
  427.     much point in calculating the avar if the clock is locked.
  428. */
  429.     numby= 0;        /*initialize the array index */
  430.     while(1)
  431.        {
  432. /*
  433.        query the primary server (the first one).  If the response is
  434.        +1 then the operation was okay and the subroutine has stored the 
  435.        time difference local-server in the x parameter of the 0th autolock 
  436.        record in units of seconds.
  437.        in addition the subroutine has filled in the MJD of the
  438.        comparison in autolock[0].day, the time in milliseconds since
  439.        0000 UTC in autolock[0].fday and the rms in units of seconds.
  440.        the nval parameter is set to the number of responses that were
  441.        actually received from the server.
  442.        if the return status is anything else then something is wrong with
  443.        that server.  Either there is something wrong with its address, it
  444.        is unreachable, it does not respond or it is not healthy.  See
  445.        getdif for details.  Try the next server (if there is one) in
  446.        any of these cases. If we run out of servers then the status
  447.        of the last one is the one that is used below.
  448. */
  449.        ntry=ntry2=ntry3=0;    /*initialize number of consecutive tries*/
  450.        fstsrv= 0;        /*point to 0th server first*/
  451.        altsrv= 0;        /*no alternate servers yet*/
  452. gogetem:
  453. /*
  454.        if the kernel routine support is available and if using
  455.        those routines is enabled, then esterr is set to 1 to
  456.        show that a calibration cycle is in progress and maxerr
  457.        starts a watch dog timer.
  458.        reset esterr as soon as we return from getdif to show
  459.        that we are not hung up there. 
  460. */
  461. #if KERNEL == 1
  462.     if(kernel == 1)
  463.        {
  464.            tmx.modes= MOD_LOCKCLOCK | MOD_MAXERROR | MOD_ESTERROR;
  465.            tmx.maxerror= -10240l;         /* -20*512 => 20 s*/
  466.            tmx.esterror=1;
  467.            syscall(SYS_ntp_adjtime,&tmx);  /*tell the kernel*/
  468.     }    /*end kernel == 1*/
  469. #endif          /*KERNEL*/
  470.  
  471. /*
  472.        start a loop over all of the time servers that we
  473.        know about.  the loop starts from fstsrv (which is
  474.        normally 0) and continues through all numsrv systems.
  475.        the operation ends with the data from the first time
  476.        server that advertises itself as healthy. if we are
  477.        currently in locked mode, then the lock loop may change
  478.        fstsrv if it looks like there is a consistency problem.
  479.        see below ...
  480. */
  481.        for(j=fstsrv; j<numsrv; j++) 
  482.           if( (i=getdif(&serv[j][0],&autolock[0],&lsr)) == 1) break;
  483.  
  484. #if KERNEL == 1
  485.     if(kernel == 1)
  486.        {
  487.        tmx.modes= MOD_LOCKCLOCK | MOD_ESTERROR;
  488.        tmx.esterror= 0;
  489.        syscall(SYS_ntp_adjtime,&tmx);
  490.     }    /*end kernel == 1*/
  491. #endif         /*KERNEL*/
  492.  
  493. /*
  494.         when we get here, we queried server number j and 
  495.         received status i.  If i is negative then we couldn't
  496.         get anything from any server.  wait a while and try
  497.         again.  If this doesn't work after 5 times (waiting
  498.         a bit longer each time), then something is really
  499.         broken.  So write out the current state and wait
  500.         for next time using a shorter interval than the
  501.         standard waitn seconds. 
  502. */
  503.         autolock[0].srv=j;      /*store server number */
  504.         autolock[0].flags=i;     /* and state of that request*/
  505.         if(i < 1)        /*error return for all servers*/
  506.            {
  507.            if(ntry++ > 5) 
  508.          {
  509.          intvl= 600;
  510.          goto finfin;    /*give up*/
  511.            }
  512.            sleep(15 + 15*ntry);    /*wait progressively longer*/
  513.            goto gogetem;        /*go try again*/
  514.         }
  515. /*
  516.         if the previous mode is 0 then this is a cold start
  517.         just store the time difference and go around again.
  518. */
  519.         if(autolock[0].mode == 0)
  520.         {
  521.         autolock[0].mode = 1;    /*advance to mode 1*/
  522.         autolock[0].rms_avg=autolock[0].rms;
  523.         autolock[0].d= 0;    /*frequency aging turned off*/
  524.         autolock[0].fcon= tconf;    /*time constant for ybar*/
  525.         intvl=autolock[0].waitn;    /*standard interval*/
  526.           goto finfin;        /*finished with this cycle*/
  527.         }      /*end of processing in mode 0*/
  528. /*
  529.         if we are still here then this is not the first cycle following
  530.         a complete cold start. compute the actual elapsed time since the 
  531.         last cycle and use this to compute the frequency of the local 
  532.         clock with respect to the calibration source. this estimate
  533.         is simple if the clock is free-running and has not been disturbed 
  534.         since the last cycle.
  535.         this estimate will not be correct if the local clock has 
  536.         been adjusted in some way since the last cycle -- if the
  537.         system has been re-booted, for example.  For the same
  538.         reason, this estimate will not be correct when we are in any 
  539.         of the lock modes, since the clock is being adjusted
  540.         all of time time there.
  541.         if we are in mode 5 then the clock is being adjusted
  542.         using the current average frequency ybar.  the effective
  543.         frequency during the last interval is then with respect to
  544.         ybar.
  545. */
  546.             autolock[0].waitx= 86400*(autolock[0].day-autolock[1].day)
  547.                 + 0.001*(autolock[0].fday - autolock[1].fday);
  548.         if(autolock[0].mode < 4)    /*clock is guaranteed free-running*/
  549.                autolock[0].y=(autolock[0].x - autolock[1].x)/autolock[0].waitx;
  550.         else if(autolock[0].mode == 4)    /*if adjustment in progress */
  551.         autolock[0].y=autolock[0].ybar;    /*then lock the frequency*/
  552.         else autolock[0].y= autolock[0].ybar +
  553.              autolock[0].x/autolock[0].waitx;
  554. /*
  555.             if this is mode 1 then this is the second cycle following
  556.             a cold start.  we have already estimated the frequency above,
  557.             so initialize the average frequency to the value we just
  558.             got. At this point, there is no way of knowing if the
  559.             data are consistent since we don't have any idea of the
  560.             rate offset of the local clock.
  561.         use the average of the rms on the two cycles for rms_avg.
  562. */
  563.             if(autolock[0].mode == 1)
  564.                 {
  565.                 autolock[0].mode=2;
  566.                 autolock[0].sumerr= 0;
  567.                 autolock[0].ybar= autolock[0].y;
  568.                 autolock[0].rms_avg= 
  569.             0.5*autolock[0].rms + 0.5*autolock[1].rms_avg;
  570.                 tloop2= 0;               /*initialize counter for mode 2*/
  571.         intvl=autolock[0].waitn;    /*standard interval*/
  572.         y[numby++]=autolock[0].y;    /*save this freq. estimate*/
  573.                 goto finfin;        /*that's all for this round*/
  574.             }     /*end of processing in mode 1*/
  575. /*
  576.         if we are still here then there have been at least two previous
  577.         cycles since the cold start, and we have used those values to
  578.         compute an initial estimate of the frequency offset of the local 
  579.         clock with respect to the calibration source.  
  580.         1.  look at the rms of these measurements and reject this
  581.             group if the rms is much larger than the last two times.
  582.         2.  use this average value to compute the prediction error
  583.             over this interval.
  584.         3.  continue to average the frequency value for the number of 
  585.             loops given by the time constant fcon.
  586.         4.  compute the average of the rms measurement noise using
  587.             a filter with a 12 hour time constant.
  588.         note that glitches at this point are not rejected -- they simply
  589.         are included in  the value for ybar and for rms.
  590.  
  591.         thus the primary purpose of mode 2 is to get an  estimate for
  592.         ybar and for the rms noise of the link. We should have a
  593.         reasonable estimate of ybar in fcon seconds, so switch
  594.         to mode 3 after that.
  595. */
  596.         if(autolock[0].mode == 2)
  597.         {
  598.         if(autolock[0].rms > 8*autolock[1].rms_avg) /*very noisy*/
  599.            {
  600.            if(ntry2++ > 5)
  601.               {
  602.               autolock[0].rms_avg= 2*autolock[1].rms_avg;
  603.               autolock[0].flags= 2;    /*show rms increase*/
  604.               intvl= 600;
  605.               goto finfin;                /*wait for next time*/
  606.            }                /*end of the fifth try*/
  607.            sleep(15+ 15*ntry2);        /*wait progressively longer*/
  608.            goto gogetem;        /*and then try again*/
  609.         }           /*end of current rms much larger than average*/
  610. /*
  611.         compute the time prediction error.  this is the 
  612.         differnce between the actual frequency measured over
  613.         the most recent interval in y and the low-passed
  614.         average value in ybar. this frequency difference is
  615.         converted to a time dispersion by multiplying by
  616.         waitx, which is the actual elapsed time over the
  617.         most recent time interval.
  618. */
  619.         autolock[0].errx=
  620.             fabs(autolock[0].y-autolock[1].ybar)*autolock[0].waitx;
  621. /*
  622.         now update ybar and sumerr using an exponential filter with 
  623.         a time constant of fcon seconds. the scaling value for this
  624.         filter is waitx/fcon -- the ratio of the actual elapsed
  625.         time over the last interval to the time constant.  both
  626.         parameters are in seconds.
  627.         since waitx is an adjustable parameter, it might 
  628.         eventually exceed fcon.  Do not allow the value of
  629.         fract to exceed 0.4 no matter what the delay is.
  630.         note that this strategy can push the effective value of
  631.         the time constant to values that may be larger than the 
  632.         optimum value given the actual Allan deviation of this 
  633.         clock. However, the alternative is allow fract to be 1
  634.         or greater, which makes the update more sensitive to 
  635.         frequency glitches.  You pays your money ...
  636. */
  637.         y[numby++]= autolock[0].y;    /*save this value*/
  638.         fract= ((double)autolock[0].waitx)/((double)autolock[0].fcon);
  639.         if(fract > 0.4) fract= 0.4;
  640.         autolock[0].ybar= (1-fract)*autolock[0].ybar +
  641.                         fract*autolock[0].y;
  642.                 autolock[0].sumerr= (1-fract)*autolock[0].sumerr +
  643.                         fract*autolock[0].errx;
  644.         fract= ((double)autolock[0].waitx)/432e+2; /*12 hr time const.*/
  645. /*
  646.         as above, do not alow the value of the time constant
  647.         to exceed 1 even if the delay becomes larger than
  648.         12 hours (which actually isn't very likely).
  649. */
  650.         if(fract > 1) fract= 1;
  651.         autolock[0].rms_avg=
  652.                         (1-fract)*autolock[1].rms_avg  +
  653.                         fract*autolock[0].rms;
  654.         tloop2 += autolock[0].waitx;    /*running time in mode 2*/
  655.         if(tloop2 > autolock[0].fcon)     /*then switch to mode 3*/
  656.            {
  657.            autolock[0].mode=3;
  658.            tloop3= 0;
  659.         }
  660.         intvl= autolock[0].waitn;
  661.         goto finfin;        /*finished for this round*/
  662.          }         /*end of processing in mode 2*/
  663. /*
  664.          the primary purpose of mode 3 is to continue to get
  665.          a better estimate of ybar and to start using it to
  666.          look at the prediction error in a serious way.
  667.          in mode 2 we calculate the prediction error but don't
  668.          do anything we these estimates beyond averaging them.
  669.          here we are going to get serious and reject a measurement
  670.          if it is not consistent with the average of the prediction 
  671.          error we have been estimating, which is being accumulated
  672.          in sumerr.
  673.          note that we have two noise parameters: rms, which measures
  674.          the scatter in a closely-spaced group of measurements and
  675.          errx, which measures the prediction error of the algorithm
  676.          based on the current value of ybar.  Both of these parameters
  677.          are averaged -- rms into rms_avg and errx into sumerr.
  678. */
  679.          if(autolock[0].mode == 3)
  680.         {
  681.         if(autolock[0].rms > 4*autolock[1].rms_avg)
  682.            {
  683.            if(ntry2++ > 5)
  684.               {
  685.               autolock[0].rms_avg= 2*autolock[1].rms_avg;
  686.               autolock[0].flags= 2;        /*show rms increase*/
  687.               intvl= 600;
  688.               goto finfin;        /*wait for next time*/
  689.            }        /*end of the fifth try*/
  690.            sleep(15+ 15*ntry2);    /*wait progressively longer*/
  691.            goto gogetem;        /*and then try again*/
  692.         } /*end of current rms larger than the average*/
  693.             fract= ((double)autolock[0].waitx)/432e+2;/*12 hr time const.*/
  694.         if(fract > 1) fract= 1;
  695.             autolock[0].rms_avg=
  696.             (1-fract)*autolock[1].rms_avg  +
  697.             fract*autolock[0].rms;
  698.         autolock[0].errx=
  699.                     fabs(autolock[0].y-autolock[1].ybar)*autolock[0].waitx;
  700.         if(autolock[0].errx > 5*autolock[0].sumerr)
  701.            {
  702.            if(ntry3++ > 5) 
  703.               {
  704.               autolock[0].sumerr = 2*autolock[1].sumerr;
  705.               autolock[0].flags=3;    /*show data noisy*/
  706.               intvl= 600;
  707.               goto finfin;      /*give up*/
  708.            }
  709.            sleep(15 + 15*ntry3);     /*wait progressively longer*/
  710.            goto gogetem;            /*go try again*/
  711.         }        /*end of prediction error very noisy*/
  712.         y[numby++]= autolock[0].y;    /*save this value*/
  713.                 fract= ((double)autolock[0].waitx)/((double)autolock[0].fcon);
  714.         if(fract > 0.4) fract= 0.4;
  715.                 autolock[0].ybar= (1-fract)*autolock[0].ybar +
  716.                         fract*autolock[0].y;
  717.         autolock[0].sumerr= (1-fract)*autolock[0].sumerr +
  718.             fract*autolock[0].errx;
  719.         eval_elap -= autolock[0].waitx;
  720.         tloop3 += autolock[0].waitx;    /*time in mode 3*/
  721. /*
  722.         if opmode is set to 2 then switch to mode 4 after
  723.         fcon seconds in mode 3.  if opmode is something else
  724.         then stay in mode 3 forever.
  725. */
  726.         if( (opmode == 2) && (tloop3 > autolock[0].fcon) )
  727.            {
  728.            autolock[0].mode =4;
  729.            tloop4= 0;      /*initialize loop counter for mode 4*/
  730.         }    /*end of time to switch to mode 4*/
  731.         intvl=autolock[0].waitn;
  732.         goto finfin;            /*that's all for this round*/
  733.          }    /*end of mode 3*/
  734.          if(autolock[0].mode == 4)        /*adjust clock to be on time*/
  735.         {
  736. /*
  737.         first check to see if the rms scatter in the current
  738.         group is reasonable compared to previous values.
  739.         if it is not then wait a bit and go around again
  740.         up to 5 times. if that still doesn't work then
  741.         increase the average value of the rms and wait
  742.         a while before trying again.
  743.         if the current scatter is reasoanble then set the
  744.         clock to be on time.  this is done in a single step
  745.         if the time offset is large and by a slew if the offset
  746.         is small.  If the offset is large, it will probably
  747.         take 2 rounds of this code to get it right -- the first
  748.         round will take out the big part of the offset and the
  749.         second round will take out the remainder.
  750.         note that the slew correction is always made using
  751.         a fixed delay of 10 minutes and that the offset of
  752.         the local clock is accounted for during this time.
  753. */
  754.             if(autolock[0].rms > 4*autolock[1].rms_avg)
  755.                    {
  756.                    if(ntry2++ > 5)
  757.                       {
  758.                       autolock[0].rms_avg= 2*autolock[1].rms_avg;
  759.                       autolock[0].flags= 2;             /*show rms increase*/
  760.                       intvl= 600;
  761.                       goto finfin;              /*wait for next time*/
  762.                    }            /*end of the fifth try*/
  763.                    sleep(15+ 15*ntry2); /*wait progressively longer*/
  764.                    goto gogetem;                /*and then try again*/
  765.                 } /*end of current rms larger than the average*/
  766. /*
  767.         if this is the first time in mode 4 then the clock
  768.         has been free-running up to this point and we can
  769.         check the data using the usual errx test.  we can't
  770.         do this after the first time because we are adjusting
  771.         the clock then and ybar is no longer a good estimator
  772.         of the average frequency over this interval.
  773. */
  774.         if(tloop4 == 0)        /*first time in mode 4*/
  775.            {
  776.            autolock[0].errx=
  777.                     fabs(autolock[0].y-autolock[1].ybar)*autolock[0].waitx;
  778.                    if(autolock[0].errx > 5*autolock[0].sumerr)
  779.                       {
  780.                       if(ntry3++ > 5)     /*give up, too many tries*/
  781.                          {
  782.                          autolock[0].sumerr = 2*autolock[1].sumerr;
  783.                          autolock[0].flags=3;      /*show data noisy*/
  784.                          intvl= 600;
  785.                          goto finfin;      /*give up*/
  786.                       }        /*end of give up, too many tries*/
  787.                       sleep(15 + 15*ntry3);    /*wait progressively longer*/
  788.                       goto gogetem;            /*and go try again*/
  789.                    }               /*end of prediction error very noisy*/
  790.         }        /*end of first time in mode 4*/
  791.         tloop4 += autolock[0].waitx;    /*actual time in mode 4*/
  792.         terr= fabs(autolock[0].x);    /*absolute value of offset*/
  793. /*
  794.                 if the kernel mods are supported and enabled then change 
  795.         the health to 2 at this point -- getting better.
  796. */
  797. #ifdef KERNEL
  798.         if(kernel == 1)
  799.                    {
  800.            tmx.modes= MOD_LOCKCLOCK | MOD_HEALTH;
  801.                    tmx.health=2;
  802.                    syscall(SYS_ntp_adjtime,&tmx);
  803.         }
  804. #endif /* KERNEL*/
  805. /*
  806.         if the time offset of the clock is large then set it
  807.         in one shot using a time step. This may not work
  808.         to better than 1 tick, and so we may be left with
  809.         a small time error.  go around again in mode 4 in
  810.         this case to check on this.
  811. */
  812.         if(terr >= 0.5)    /*step clock for large time offset*/
  813.            {
  814.            tcor= -autolock[0].x;
  815.            j=sett(tcor);
  816.            if(j != 0)    /*error trying to set clock*/
  817.              {
  818.              sprintf(errmsg,"Error %d setting clock by sett.",j);
  819.              putmsg(errmsg);
  820.              break;
  821.            } 
  822.            intvl=120;    /*wait just a short time*/
  823.            goto finfin;    /*and go around again*/
  824.         }
  825. /*
  826.         if the time error is significant in the sense that it is
  827.         larger than the average rms of the measurement noise but
  828.         at the same time is small enough to be removed by a clock 
  829.         slew rather than a step then make a final adjustment
  830.         using the clock slew system. use a fixed delay of 600s
  831.         for this, and include the fact that the clock will drift
  832.         because of its frequency offset during this adjustment
  833.         period.  
  834. */
  835.         if(terr >= 2*autolock[0].rms_avg)/*small but significant*/
  836.            {
  837.            tcor= -autolock[0].x -600*autolock[0].ybar;
  838.            j=adjt(tcor);/*slew clock for small but significant offset*/
  839.            if(j != 0)
  840.               {
  841.               sprintf(errmsg,"Error %d adjusting clock by adjt.",j);
  842.               putmsg(errmsg);
  843.               break;
  844.            }
  845.            intvl= 600;
  846. /*
  847.                    if the kernel mods are supported and enabled then change
  848.                    the health to 1 at this point -- getting better.
  849. */
  850. #ifdef KERNEL
  851.                    if(kernel == 1)
  852.                       {
  853.                       tmx.modes= MOD_LOCKCLOCK | MOD_HEALTH;
  854.                       tmx.health=1;
  855.                       syscall(SYS_ntp_adjtime,&tmx);
  856.                    }
  857. #endif /* KERNEL*/
  858. /*
  859.                     if we have been adjusting the clock for some
  860.                     time now and we are still here, then it is
  861.                     probable that the frequency that we are using
  862.                     is not quite right. This will never be fixed
  863.                     in this mode, since we don't update the frequency
  864.                     here. fall through in this case and switch to mode 5.
  865.  
  866.                     this is most likely to happen when the rms_avg is
  867.                     very small, since even a small time offset will
  868.                     be treated as "significant" by the test above.
  869. */
  870.                    if(tloop4 < 3*autolock[0].waitn) goto finfin;
  871.         }        /*end of small but significant*/
  872. /*
  873.         finally, switch to the mode 5 lock system if the clock time 
  874.         offset is small enough that it is no longer statistically 
  875.         significant.
  876.         although we have just set the clock, it is going to
  877.         start drifting off again at a rate of ybar. since we
  878.         are not yet controlling the frequency of the clock,
  879.         we are going to find an error on the next round of
  880.         intvl*ybar, where intvl is the length of time we
  881.         wait. choose intvl so that this dispersion will on the 
  882.         order of the average measurement error which is rms_avg,
  883.         but not more than 10 minutes.
  884. */
  885.         autolock[0].mode= 5;    /*switch to locked mode*/
  886.         eval_elap= EVALPERIOD;    /*reset evaluation timer*/
  887.           intvl= fabs(autolock[0].rms_avg/autolock[0].ybar);
  888.         if(intvl > 600) intvl= 600;
  889. /*
  890.         do not let the interval between queries become too
  891.         small even if the rms error goes to 0
  892. */
  893.         if(intvl <  intvlmin) intvl= intvlmin; 
  894.         goto finfin;            /*done with this round*/
  895.          }    /*end of processing in mode 4 -- set clock to be on time*/
  896. /*
  897.          when we get to mode 5, we are going to use the calibration
  898.          data to keep the local clock locked.  The first step is
  899.          to perform the usual sanity checks on the data.
  900. */
  901.          if(autolock[0].mode == 5)    /*maintain lock*/
  902.         {
  903. /*
  904.         if the rms of the current group of measurements
  905.         is much larger than the average rms then go around
  906.         and try again, waiting a little bit longer each time.
  907.  
  908.         if the current rms is consistent with the previous
  909.         values, then update the rms with a 12-hour time
  910.         constant as above.
  911. */
  912.               if(autolock[0].rms > 4*autolock[1].rms_avg)
  913.                    {
  914.                    if(ntry2++ > 5)
  915.                       {
  916.                       autolock[0].rms_avg= 2*autolock[1].rms_avg;
  917.                       autolock[0].flags= 2;             /*show rms increase*/
  918.                       intvl= 600;
  919.                       goto finfin;              /*wait for next time*/
  920.                    }            /*end of the fifth try*/
  921.                    sleep(15+ 15*ntry2); /*wait progressively longer*/
  922.                    goto gogetem;                /*and then try again*/
  923.                 } /*end of current rms larger than the average*/
  924.                 fract= ((double)autolock[0].waitx)/432e+2;/*12 hr time const.*/
  925.                 if(fract > 1) fract= 1;
  926.                 autolock[0].rms_avg=
  927.                         (1-fract)*autolock[1].rms_avg  +
  928.                         fract*autolock[0].rms;
  929. /*
  930.         as above, the prediction error is the difference between 
  931.         the effective freqency over the last interval and 
  932.         the average frequency, where this frequency difference is 
  933.         converted to a time by multiplying by the time interval
  934.         since the last calibration cycle.  The local clock has
  935.         been steered by applying ybar to it. since the observed
  936.         frequency over the last interval is ybar + x/waitx, the
  937.         difference between the two multipled by waitx really reduces 
  938.         to just x.  Another way of saying this is that the predicted 
  939.         value of x is 0 and the prediction error is therefore whatever 
  940.         non-zero value we find.  note that this will not be true on
  941.         the first round through this part, since the local clock
  942.         frequency is not controlled yet.
  943. */
  944.         autolock[0].errx= autolock[0].x;
  945. /*
  946.         the algorithm splits into two parts at this point.  If
  947.         the prediction error is larger than the desired performance
  948.         level as specified by sigma_desired, then the prediction
  949.         error is treated as a problem because the algorithm is
  950.         not working as well as was specified.  The first question
  951.         in this case is whether or not this big prediction error
  952.         is a one-time glitch or the real thing. To test for this
  953.         possibility, try the whole thing all over again, waiting 
  954.         a bit between each try.  if this doesn't fix the problem 
  955.         after a number of tries, then more serious analysis is
  956.         needed.  there are four possibilities:
  957.  
  958.         1.  a local time step
  959.         2.  a local frequency step
  960.         3.  the distant server is broken.
  961.         4.  a network glitch
  962.  
  963.         first check possibilities 3 and 4 by switching to
  964.         another server, if we know about one and if we
  965.         haven't tried it yet. the current data were received
  966.         from autolock[0].srv. this will normally point to
  967.         server 0, unless that one was unreachable or unhealthy.
  968.         set the fstsrv pointer to the one after that if we
  969.         know about it, reset the flag and go around again.
  970.         save the time difference from the current server in
  971.         array alterr and the rms scatter of the measurements
  972.            in altrms.
  973.  
  974.         if, on the other hand, the prediction error is not
  975.         larger than the desired performance level, then this
  976.         is not necessarily a glitch.  The prediction error
  977.         will always be larger than the rms_avg in this case
  978.         because we are pushing the averaging time to longer
  979.         values than is optimum based on white-noise considerations,
  980.         and the average prediction error will no longer be
  981.         a good expectation of the observed data.
  982.         note that the units of sigma_desired are in milliseconds,
  983.         whereas everything else is in seconds.
  984. */
  985.                 if( (fabs(autolock[0].errx) > 5*autolock[0].rms_avg) &&
  986.             (1000*fabs(autolock[0].errx) > 3*sigma_desired) )
  987.                    {
  988.                    if(ntry3++ > 5)    /*looks like it is not going away*/
  989.                       {
  990.               altrms[altsrv]  = autolock[0].rms;
  991.               alterr[altsrv++]= autolock[0].errx;
  992.               fstsrv= autolock[0].srv + 1;    /* point to next one*/
  993.               if(fstsrv < numsrv)        /* it is available*/
  994.             {
  995.             ntry3= 0;    /*reset the loop counter*/
  996.             goto gogetem;    /*try this one*/
  997.               }
  998. /*
  999.               when we get here we have gotten a prediction error
  1000.               using all ofthe servers that we know about. the
  1001.               array alterr has all of the prediction errors that
  1002.               we measured and altsrv is the number of them we have
  1003.               received. if there are several of these, see
  1004.               if they are consistent. if the data from the different
  1005.               servers seem to agree then assume we have a time
  1006.               step. If this is the second time step in a row then
  1007.               it is more likely a frequency step.
  1008. */
  1009.               if(altsrv > 1)        /*more than one server*/
  1010.             {
  1011.             terr= 0;    /*initialize the mean */
  1012.             for(j=0; j < altsrv; j++) terr += alterr[j];
  1013.             terr /= altsrv;        /*average difference*/
  1014.             tsig= 0;    /*intialize the std. dev.*/
  1015.             for(j=0; j < altsrv; j++) /*compute sum square*/
  1016.                {
  1017.                ddum= alterr[j]-terr;
  1018.                tsig += ddum*ddum;
  1019.             }    /*end of for loop to compute sum square*/
  1020.             tsig= sqrt(tsig/(altsrv-1)); /*normalize ss to sigma*/
  1021.             if(fabs(terr) > 2*tsig)        /*looks real*/
  1022.                {
  1023.                autolock[0].x= terr;        /*average value*/
  1024.                autolock[0].rms= tsig;
  1025.                autolock[0].flags= 4;    /*time step*/
  1026.                if(autolock[1].flags == 4) autolock[0].flags=6;
  1027.                goto update;            /*adjust local clock*/
  1028.             }        /* end of looks real*/
  1029.             else        /*multiple servers don't agree*/
  1030.                {
  1031. /*
  1032.                write out all of the time messages for
  1033.                later debugging.
  1034. */
  1035.                  for(j=0; j<altsrv; j++)
  1036.                                 {
  1037.                                 sprintf(errmsg,"server %d diff=%le",
  1038.                                         j,alterr[j]);
  1039.                                 putmsg(errmsg);
  1040.                }
  1041.                            sprintf(errmsg,"mean = %le var=%le",terr,tsig);
  1042.                            putmsg(errmsg);
  1043. /*
  1044.                if the values from the different servers
  1045.                don't agree then either one of the distant
  1046.                servers is broken or else one of the paths has
  1047.                become very asymmetric so that the delay
  1048.                estimate has been fooled.
  1049.                there is nothing we can do about this if
  1050.                we have only two other servers. set the flag
  1051.                to show the problem, and wait a bit to
  1052.                see if the problem gets fixed.
  1053. */
  1054.                if(altsrv == 2)   /*only 2 servers -- tough luck*/
  1055.                  {
  1056.                  autolock[0].flags= 5;    /*servers don't agree*/
  1057.                  intvl= 600;
  1058.                  goto finfin;    /*wait a while and try again*/
  1059.                }        /*end of only two servers*/
  1060. /*
  1061.                if we are here then we have more than 2 servers
  1062.                and they don't agree.  see if we can find an
  1063.                outlier. recompute the mean time difference
  1064.                dropping a different server each time.  in
  1065.                the following loops, we drop server j1 on each
  1066.                round and then loop j2 over the remaining ones
  1067.                to compute the mean.
  1068. */
  1069.                for(j1=0; j1< altsrv; j1++)
  1070.                 {
  1071.                 terr=0;    /*initialize the mean value*/
  1072.                 for(j2=0; j2< altsrv; j2++)
  1073.                    {
  1074.                    if(j2 == j1) continue; /*drop this server*/
  1075.                    terr += alterr[j2];
  1076.                 }    /*end of loop over remaining ones*/
  1077.                 terr /= altsrv - 1; /*we dropped one*/
  1078.                                tsig= 0;        /*intialize the std. dev.*/
  1079.                             for(j2=0; j2 < altsrv; j2++) /*compute sum sq.*/
  1080.                                   {
  1081.                    if(j2 == j1) continue;  /*drop this one*/
  1082.                                   ddum= alterr[j2]-terr;
  1083.                                   tsig += ddum*ddum;
  1084.                             }       /*end of for loop to compute sum sq*/
  1085.                             tsig= sqrt(tsig/(altsrv-2)); /*ss -> sigma*/
  1086.                 sprintf(errmsg,
  1087.                 "Drop server %d new mean, sigma= %le %le",
  1088.                     j1,terr,tsig);
  1089.                 putmsg(errmsg);
  1090.                                if(fabs(terr) > 2*tsig)         /*looks real*/
  1091.                                   {
  1092.                                   autolock[0].x= terr;        /*average value*/
  1093.                                   autolock[0].rms= tsig;
  1094.                                   autolock[0].flags= 4;        /*time step*/
  1095.                                   if(autolock[1].flags == 4) 
  1096.                         autolock[0].flags=6;
  1097.                    sprintf(errmsg,"drop server %d",j1);
  1098.                    putmsg(errmsg);
  1099.                    sprintf(errmsg," new mean= %le new sig= %le",
  1100.                     terr,tsig);
  1101.                    putmsg(errmsg);
  1102.                                   goto update;           /*adjust local clock*/
  1103.                 } /*end of looks real*/
  1104. /*
  1105.                 if we get to here then the servers don't
  1106.                 agree and dropping one of them didn't seem
  1107.                 to help. looks like that's all for now 
  1108.                 folks.
  1109. */
  1110.                               autolock[0].flags= 5;
  1111.                                 intvl= 600;
  1112.                                 goto finfin;     /*wait a while and try again*/
  1113.                }        /*end of loop to drop server j1*/
  1114.             }        /*end of multiple servers dont agree*/
  1115.               }            /*end of more than one server*/
  1116. /*
  1117.               if we get here then either we only know about one
  1118.               server or we can only find one server that is healthy.
  1119.               it is not clear who is right in this case, so
  1120.               do not use these data to update our frequency or our
  1121.               time estimates. increase the prediction error tolerance, 
  1122.               set the flag to show that the data are noisy and go 
  1123.               around again.
  1124. */
  1125.                       autolock[0].rms_avg= 2*autolock[1].rms_avg;
  1126.                       autolock[0].flags=3;      /*show data noisy*/
  1127.                       intvl= 600;
  1128.                       goto finfin;      /*give up*/
  1129.                    }        /*end of the error persists after 5th try*/
  1130.                    sleep(15 + 15*ntry3);     /*wait progressively longer*/
  1131.                    goto gogetem;            /*go try again*/
  1132.                 }               /*end of prediction error very noisy*/
  1133. /*
  1134.         when we get to update then either the data received on
  1135.         the current cycle are about what we expected or else
  1136.         we have confirmed time step because several servers
  1137.         agree on the time difference.
  1138. */
  1139. update:            
  1140.                 fract= ((double)autolock[0].waitx)/((double)autolock[0].fcon);
  1141.                 if(fract > 0.4) fract= 0.4;
  1142. /*
  1143.         compute the updated frequency, noting that in this
  1144.         mode, y= ybar + x/waitx. using the usual exponential filter,
  1145.         the update to the frequency is simply x/waitx scaled
  1146.         by the time constant.
  1147.         do not update the frequency if iflags == 4, since that
  1148.         means we have a time step. the frequency estimate is 
  1149.         updated if we think this is a frequency step.
  1150. */
  1151.         if(autolock[0].flags != 4) autolock[0].ybar +=
  1152.             fract*autolock[0].x/autolock[0].waitx;
  1153.                 autolock[0].sumerr= (1-fract)*autolock[0].sumerr +
  1154.                         fract*fabs(autolock[0].errx);
  1155.                 eval_elap -= autolock[0].waitx;
  1156.         intvl=autolock[0].waitn;    /*default delay*/
  1157. #ifdef KERNEL
  1158. /*
  1159.         if kernel mods are supported and enabled, then our estimated 
  1160.         frequency offset can be transmitted directly to the kernel 
  1161.         where it will be done automatically. the units of the kernel 
  1162.         frequency variable are ps/s or e-12, whereas ybar is
  1163.         dimensionless and thus has units of s/s.
  1164. */
  1165.         if(kernel == 1)
  1166.            {
  1167.            tmx.modes=MOD_LOCKCLOCK | MOD_FREQUENCY | MOD_HEALTH | 
  1168.                     MOD_STATUS;
  1169.            tmx.freq= -1e+12*autolock[0].ybar;
  1170.            tmx.health=0;
  1171.            tmx.status=TIME_OK;
  1172. /*
  1173.            if the leap second status from the remote server is
  1174.            different from our local notion of the leap second
  1175.            flag, then do the right thing.
  1176. */
  1177.            if(lsr != lsl)    /*if flags are different*/
  1178.               {
  1179.               lsl=lsr;    /*set local to value of remote*/
  1180.               if(lsl == 1) tmx.status |= STA_INS;    /*insert*/
  1181.               if(lsl == 2) tmx.status |= STA_DEL;    /*drop*/
  1182.            }        /*end of flags are different*/
  1183. /*
  1184.            remove a time error if it is significantly 
  1185.            larger than the uncertainty of the calibration
  1186.            data.
  1187.            the offset parameter is in units of ps. it is
  1188.            defined as a long int, which means that it is
  1189.            64 bits on an Alpha.  It is not likely that any
  1190.            reasonable time offset could cause a 64 bit time 
  1191.            difference to overflow -- even if it is in units of ps.
  1192.            nevertheless, limit it just to be sure.
  1193. */
  1194.            if(fabs(autolock[0].x) > autolock[0].rms)
  1195.             {
  1196.             tmx.modes |= MOD_OFFSET;
  1197.             if(autolock[0].x >  1e+7) autolock[0].x=  1e+7;
  1198.             if(autolock[0].x < -1e+7) autolock[0].x= -1e+7;
  1199.             tmx.offset= -1e+12*autolock[0].x; /*convert to ps*/
  1200.            }
  1201.            syscall(SYS_ntp_adjtime,&tmx);
  1202.         }    /*end of if kernel == 1*/
  1203. #endif        /*KERNEL*/
  1204. #ifdef LINUX
  1205. /*    
  1206.     if LINUX is defined then we have to do the same jobs as
  1207.     above, but we do them slightly differently.
  1208.  
  1209.     First, if the leap second status from the remote server is
  1210.         different from our local notion of the leap second
  1211.         flag, then do the right thing.
  1212.  
  1213.     then remove a time error as above. Note that the units
  1214.     for time offset in LINUX are microseconds.
  1215. */
  1216.         tmx.status= TIME_OK;
  1217.         tmx.modes=ADJ_STATUS;
  1218.         if(lsr != lsl)       /*if flags are different*/
  1219.                    {
  1220.                    lsl=lsr;  /*set local to value of remote*/
  1221.                    if(lsl == 1) tmx.status |= STA_INS;       /*insert*/
  1222.                    if(lsl == 2) tmx.status |= STA_DEL;       /*drop*/
  1223.                 }            /*end of flags are different*/
  1224.         if(fabs(autolock[0].x) > autolock[0].rms)
  1225.            {
  1226.            tmx.modes |= ADJ_OFFSET_SINGLESHOT;
  1227.            if(autolock[0].x > 100) tmx.offset=-1e+7;
  1228.            else if(autolock[0].x < -100) tmx.offset= 1e+7;
  1229.            else tmx.offset=-1e+6*autolock[0].x;
  1230.         }    /*end of time adjustment is significant*/
  1231.         adjtimex(&tmx);
  1232. #endif    /*LINUX*/
  1233.         goto finfin;            /*end of this round*/
  1234.          }
  1235. finfin:        /*write current states to file, ... */
  1236.             if( (lop=fopen(SNAME,"w")) == NULL) 
  1237.                {
  1238.                sprintf(errmsg,"Oops -- error opening state file %s.\n",SNAME);
  1239.            putmsg(errmsg);
  1240.                exit(-1);
  1241.         }
  1242.         putautolock(&lop,autolock);    /*save current state*/
  1243.         fclose(lop);
  1244. /*
  1245.         if we are running in opmode = 1 then the purpose is
  1246.         simply to evaluate the network and the clock.  This is
  1247.         done by running for a fixed time and then exiting.
  1248.         stop the music when this time expires.
  1249. */
  1250.         if( (opmode == 1) && (ttot >= evalint) ) break;
  1251. /*
  1252.         see if it is time to examine how well we have been
  1253.         doing for the last day.
  1254. */
  1255.         if(eval_elap <= 0)        /*yes, time for a retrospective look*/
  1256.         {
  1257. /*
  1258.         if the local clock is still free-running then
  1259.         estimate the Allan variance using the y array.
  1260.         the estimate is made for lags that are powers of 2.
  1261.         and the operation ends when the adev subroutine 
  1262.         returns 0 meaning that the lag is now too big
  1263.         for the number of frequency estimates currently
  1264.         in the y array.
  1265.         note that in the following, the first calculation
  1266.         is done with lag = 1 and is stored in the avar
  1267.         array at index 0.
  1268.         if we are in mode 4 or later then the local clock
  1269.         is being adjusted and the calculation doesn't mean
  1270.         much.  so skip it.
  1271. */
  1272.         if(autolock[0].mode >= 4) goto novar;    /*skip adev if yes*/
  1273.         nvar=0;        /*initialize avar index*/
  1274.         for(j=1; j < NUMSTATES; j *= 2)     /*j sets the lag*/
  1275.             {
  1276.             if( (dtemp=adev(y,numby,j)) == 0) break;    
  1277.             sprintf(errmsg," Adev at lag %d is %le",
  1278.             j*autolock[0].waitn,dtemp);
  1279.             putmsg(errmsg);
  1280.             avar[nvar++]= dtemp;    /*save this value*/
  1281.         }        /*end of adev calculation*/
  1282. /*
  1283.         the avar calculations are used to evaluate fcon 
  1284.         which is the effective time constant for averaging
  1285.         y to produce ybar. the value of nvar gives the
  1286.         number of adev calculations computed above.
  1287.         look for the break-point in the variation of adev
  1288.         with lag, but don't do anything if nvar <= 1, which
  1289.         means we have fewer than two estimates.
  1290.         we expect that avar will decrease at least as fast
  1291.         as the square root of j before we get into the flicker
  1292.         FM domain. look for the break in the slope that signals
  1293.         that this.
  1294.         if we can find this point see if it is roughly equal
  1295.         to fcon, and adjust fcon if this is not true.  this
  1296.         is only a rough calculation, and don't change things
  1297.         if fcon is right to within a factor of 2.
  1298.         note that the j-th estimate is calculated using an 
  1299.         averaging time of 2^(j)*waitn, where the first estimate 
  1300.         uses j=0 -- see above.
  1301.         if this routine changes the value then put a message
  1302.         in the log file giving the new value.
  1303. */
  1304.         for(j=1; j<nvar; j++) /*loop over all adev calculations*/
  1305.            {
  1306.            if(avar[j-1]/avar[j] < 1.4) break;    /*break is here*/
  1307.         }        /*end of loop over all adev calculations*/
  1308.         if( (j > 1) && 
  1309.             (j < nvar))        /*break found inside array*/
  1310.            {
  1311.            itemp= autolock[0].waitn << j; /*multiply by 2^j*/
  1312.            if(autolock[0].fcon > 2*itemp) autolock[0].fcon=itemp;
  1313.            else 
  1314.             if(autolock[0].fcon < itemp/4) autolock[0].fcon=itemp/2;
  1315.            if(autolock[0].fcon != autolock[1].fcon)
  1316.               {
  1317.               sprintf(errmsg,
  1318.                 "value for fcon is now %d",autolock[0].fcon);
  1319.               putmsg(errmsg);
  1320.            }    /*end of write out value if different*/
  1321.         }    /*end of break found inside array*/
  1322. novar:            /*mode >= 4, don't compute Allan variance*/
  1323. /*
  1324.         first question: is the algorithm providing an average
  1325.         prediction error that is better, worse or about right
  1326.         relative to the desired performance?
  1327.         if not root 2 better and not root 2 worse then leave
  1328.         the parameters alone -- they must be about right.
  1329.         note that the units of sumerr are seconds whereas
  1330.         sigma_desired is in milliseconds.
  1331.  
  1332.         if the performance is more than root 2 better then
  1333.         increase the time interval by a factor of 2 unless
  1334.         we are already at the maximum value specified in the
  1335.         parameter file, in which case don't increase past
  1336.         that point.
  1337.         note that it is quite likely that the prediction error
  1338.         on the next cycle will be worse than we think because
  1339.         we have just increased averaging time.  Increase
  1340.         the value of sumerr in this case as well. If the 
  1341.         interval is already at the maximum value then don't
  1342.         increase sumerr.
  1343. */
  1344.         if(1414*autolock[0].sumerr < sigma_desired) /*root 2 better*/
  1345.            {
  1346.            numby= 0;        /*reset y index for new waitn*/
  1347.            autolock[0].waitn *= 2; /*increase interval between cycles*/
  1348.            if(autolock[0].waitn > intvlmax) 
  1349.             autolock[0].waitn= intvlmax; /*not more than max*/
  1350.            else
  1351.             autolock[0].sumerr *= 1.4;
  1352.            intvl= autolock[0].waitn;    /*set new value for delay*/
  1353.         }     /*end of performance is better*/
  1354.         else if(707*autolock[0].sumerr > sigma_desired) /*root 2 worse*/
  1355.            {
  1356.            numby= 0;        /*reset y index for new waitn*/
  1357.            autolock[0].waitn /= 2;/*decrease interval between cycles*/
  1358.            if(autolock[0].waitn < intvlmin) 
  1359.             autolock[0].waitn= intvlmin; /*but not less than min*/
  1360.            intvl= autolock[0].waitn;    /*set new value for delay*/
  1361.         }    /*end of performance is worse*/
  1362. /*
  1363.         second question:
  1364.         is prediction error better, worse or about the same as
  1365.         average measurement error on each cycle. if the two are
  1366.         about the same then the process is dominated by measurement
  1367.         noise and the clock noise is not too important. This means
  1368.         that we might improve the algorithm by increasing the
  1369.         number of measurements in each group.
  1370.         tests:
  1371.         
  1372.         1.  if the measurement noise is much less than the prediction
  1373.             error then we are dominated by clock noise and the extra
  1374.             measurements are not worth the trouble. So decrease
  1375.             the number of measurements in each group until we get
  1376.             down to the minimum.
  1377.         2.  if the measurement noise is greater than the desired
  1378.             performance, then the loop is not doing as well as
  1379.             we specify.  increase the number of calibrations in
  1380.             each group until we get up to the maximum.
  1381. */
  1382.         if(4*autolock[0].rms_avg < autolock[0].sumerr)
  1383.            {
  1384.            ncal /= 2;
  1385.            if(ncal < ncalmin) ncal= ncalmin;
  1386.         }
  1387.         else if(autolock[0].rms_avg > 2*sigma_desired)
  1388.            {
  1389.            ncal *= 2;
  1390.            if(ncal > ncalmax) ncal= ncalmax;
  1391.         }
  1392.         eval_elap= EVALPERIOD;    /*reset the counter for the next round*/
  1393.         }        /*end of retrospective look at performance*/
  1394. /*
  1395.         display time of the next calibration cycle
  1396.         if not in debug mode. note that this value will
  1397.         be wrong if the calibration cycle ended in an error
  1398.         since fday was not updated and still contains the
  1399.         time of the last successful cycle.
  1400. */
  1401.         if( (argc > 1) && (debug == 0) )
  1402.         {
  1403.         cp=argv[1];
  1404.             while(*cp != '\0') *(cp++) = ' ';
  1405.         nmin=(autolock[0].fday/1000 + intvl)/60; /*convert to minutes*/
  1406.         nhr=nmin/60;
  1407.         nmin %= 60;
  1408.         sprintf(argv[1],"-> %2d:%2.2d",nhr,nmin);
  1409.         }
  1410. /*
  1411.         if the kernel routines are present and enabled then set
  1412.         a watch dog timer in the kernel for the number of seconds
  1413.         that we are going to sleep.  if that timer becomes positive
  1414.         then the time has expired. if we don't reset it then we
  1415.         have crashed. since our calculated frequency offset is
  1416.         being removed automatically in this mode, we can just
  1417.         go back to sleep without doing anything else.
  1418. */
  1419. #ifdef KERNEL
  1420.         if(kernel == 1)
  1421.           {
  1422.               tmx.modes= MOD_LOCKCLOCK | MOD_MAXERROR;
  1423.               tmx.maxerror= -512*(intvl);  /*reset time out counter*/
  1424.               syscall(SYS_ntp_adjtime,&tmx);  /*tell the kernel*/
  1425.           sleep( (unsigned int) intvl);    /*the sleep of the just*/
  1426.           ttot += intvl;        /* add this to the running time*/
  1427.           goto awake;        /* rise and shine! */
  1428.         }    /*end kernel == 1*/
  1429. #endif          /*KERNEL*/
  1430. /*
  1431.         if we get here then kernel mods are either not supported
  1432.         on this system or they have been disabled.
  1433.         if we are not in a lock mode, then we can just sleep the
  1434.         uninterrupted sleep of the just until next time.
  1435. */
  1436.         if(autolock[0].mode < 5)
  1437.         {
  1438.             sleep( (unsigned int) intvl); /*and sleep until next time*/
  1439.             ttot += intvl;        /*add this sleep interval */
  1440.         goto awake;
  1441.         }
  1442. /*
  1443.         if we get here then kernel mods are not supported and we 
  1444.         are in a lock mode. Instead of sleeping through the whole 
  1445.         time we have to schedule periodic corrections to the local 
  1446.         clock to keep it on time based on our estimate of ybar.
  1447.             find out how long we must wait before a correction of 0.1 ms 
  1448.         will be needed.  Then adjust this interval so that there will 
  1449.         be an integral number of them in the calibration delay and then
  1450.             compute the exact correction needed each time.
  1451.  
  1452.            if the offset frequency is very small, then the correction needed 
  1453.        in the full iwait seconds may be less than 1 ms.  Compute the 
  1454.        correction needed in this case and apply it in two steps.
  1455.  
  1456.            if the offset frequency is very large, then the delay between 
  1457.        0.1 ms corrections may be less than 1 sec.  Since the maximum 
  1458.        slew rate is on the order of milliseconds/second, the best we can 
  1459.        do is to delay for 1 second and then adjust at the maximum
  1460.            slew rate.
  1461.        in all cases, fall through to awake when the adjustments are
  1462.        finished.
  1463. */
  1464.        jwait= fabs(0.0001/autolock[0].ybar);    /*time for 0.1 ms*/
  1465.        if(jwait < 1) jwait= 1;    /*not less than 1 second*/
  1466.        if(jwait > intvl) jwait=intvl/2;    /*if longer than whole cycle*/
  1467.        jloop= intvl/jwait;        /*number of adjustments in whole cycle*/
  1468.        tcor= -(double)jwait*autolock[0].ybar;
  1469.        for(i=0; i< jloop; i++)    /*correction loop*/
  1470.         {
  1471.         adjt(tcor);
  1472.         sleep(jwait);
  1473.         }
  1474. awake:        /*come here when we wake up after our sleep*/
  1475. /*
  1476.         when we wake up, we have to push the autolock states down
  1477.         by 1; the data acquired on the next cycle is then put in
  1478.         autolock[0]. note that the values previously in the last
  1479.         state will be lost and the values previously in the 0th
  1480.         slot are duplicated into the first slot.
  1481. */
  1482.          pushautolock(autolock,0,NUMSTATES);
  1483. /*
  1484.          slide y array down 1 slot if it is full
  1485.          numby incremented after each transfer and thus points to
  1486.          next free slot. Thus array is full when numby == NUMSTATES
  1487.          test is on >= just in case we miss ==.
  1488.          once we switch to any mode after 3 then numby doesn't 
  1489.          change anymore since there is no point evaluating the
  1490.          frequency stability of the clock when it is being
  1491.          adjusted.
  1492. */
  1493.          if(numby >= NUMSTATES)    /*y array is full*/
  1494.         {
  1495.         for(i=1; i < NUMSTATES; i++) y[i-1]=y[i];
  1496.         numby= NUMSTATES - 1;     /*last element in array*/
  1497.          }        /*end of y array is full*/
  1498. /*
  1499.          decrement filefull each time we push the autolock down,
  1500.          and copy the file when filefull is decremented to 
  1501.          a small value. This will prevent data from getting
  1502.          lost by being pushed off the bottom of the file. 
  1503. */ 
  1504.       
  1505.              if(filefull-- < 20)           /*file is almost full*/
  1506.                 {
  1507. #if ALPHA == 1
  1508.                 sprintf(filenam,"cp %s %s_%d",SNAME,SNAME,autolock[0].day);
  1509. #else
  1510.                 sprintf(filenam,"cp %s %s_%ld",SNAME,SNAME,autolock[0].day);
  1511. #endif
  1512.                 system(filenam);
  1513.                 pid=wait(status);         /*wait for copy to finish*/
  1514.                 filefull = NUMSTATES;     /*reset for next round */
  1515.           }     /*end of file is almost full*/
  1516.     }        /* end of the main infinite while loop*/
  1517.     exit(0);    /* can get here only on a break out of while loop*/
  1518. }
  1519.