home *** CD-ROM | disk | FTP | other *** search
/ nisttime.carsoncity.k12.mi.us / nisttime.carsoncity.k12.mi.us.tar / nisttime.carsoncity.k12.mi.us / pub / lockclock / lockclock.c < prev    next >
C/C++ Source or Header  |  1996-11-18  |  34KB  |  946 lines
  1. main(argc,argv)
  2. int argc;
  3. char *argv[];
  4. {
  5. #include <stdio.h>
  6. #include <math.h>
  7. #include <signal.h>
  8. #include <stdlib.h>
  9. #include <sys/types.h>
  10. #include <sys/stat.h>
  11. #include <fcntl.h>
  12. #include <sys/ioctl.h>
  13. #include <sys/param.h>
  14. #include <sys/wait.h>
  15. #include <sys/time.h>
  16. #include "timex.h"
  17. extern int errno;
  18. /*
  19.     this program locks the local clock using calibration
  20.     data received from some external source.  the actual
  21.     calibration process is implemented in a separate
  22.     subroutine and the details are not known to this
  23.     program.  the subroutine must return the current
  24.     time as a day number and time in seconds, followed
  25.     by the measured time difference, where a positive value
  26.     means that the local clock is fast.
  27. */
  28. FILE *lop;            /*used to open file lock.dat*/
  29. FILE *pop;            /*used to open file param.dat*/
  30. FILE *fopen();
  31. #include "lock.h"        /*get data structure definitions*/
  32. #include "sizint.h"
  33. struct lockst lock[NUMSTATES];    /*push-down stack of the system state*/
  34. struct lockpm  param;        /*lock parameter structure */
  35. char *rootdir="/local/etc";    /*directory for lockclock files*/
  36. char *sname="lock.dat";        /*name of state file */
  37. char *pname="param.dat";    /*name of parameter file*/
  38. void getlock();            /*reads state file into structure*/
  39. void pushlock();        /*push down lock structure by 1*/
  40. void lockerr();            /*prints message to file and exits*/
  41. void putleap();            /*writes leap second file for other tasks*/
  42. int adjt();            /*adjusts local time, returns 0 if ok*/
  43. int sett();            /*steps local time, returns 0 if ok*/
  44. LONG gettck();            /*returns usec when tick received */
  45. LONG jadj;            /*value returned by gettck*/
  46. LONG ladj = 0;            /*last value of jadj*/
  47. int i;
  48. LONG daym,fdaym;        /*time of calibration, MJD and seconds*/
  49. double delta;            /*time difference local-remote in s*/
  50. int lsflag;            /*leap second flag from getdif*/
  51. int lsdue;            /*leap second is due this month*/
  52. double sigmam;            /*sigma of differences returned by getdif*/
  53. int getdif();            /*returns time difference local-remote*/
  54. int j;
  55. unsigned int iwait;        /*nominal delay until next cycle*/
  56. unsigned int jwait;        /*nominal delay between time adjustments*/
  57. unsigned int kwait;        /*sleep interval for big adj. in mode 7*/
  58. double tcor,lptcor;        /*time adjustments in seconds */
  59. int jloop;            /*number of adjustments for each cycle*/
  60. double xavgt;            /*y-time constant normalized by delay*/
  61. char msg[80];            /*message to environment showing status*/
  62. char *s;            /*message pointer */
  63. double ddum;
  64. int conserr = 0;        /*number of consecutive errors from getdif*/
  65. int conterr;            /*number of consecutive errors from gettck*/
  66. int filefull = NUMSTATES;    /*counter to see when file is getting full*/
  67. int nloops = 0;            /*number of loops since dead-start */
  68. char filenam[50];        /*holds new file name copy command*/
  69. int pid;            /*return value from  wait call*/
  70. union wait *status;        /*return status from wait call*/
  71. int fstep;            /*reset counter in freq-step detector*/
  72. int fsteplim;            /*max value of fstep */
  73. struct timex tmx;        /*used for David Mills kernel calls */
  74. /*
  75.     parameter slew gives the time in seconds necessary to
  76.     correct the local clock by 1 second when subroutine
  77.     adjtime is used.  For ULTRIX/RISC, for example, adjtime
  78.     moves the clock 15 usec/tick.  Since there are 256
  79.     ticks/second in this system, the clock moves 256X15 =
  80.     3840 usec/sec, so that it takes about 260 seconds to
  81.     move 1 second.  This parameter is needed to know how
  82.     long it will take to adjust the clock by a big amount,
  83.     which is usually necessary only at startup.
  84. */
  85. int slew;
  86. /*
  87.     the first job is to convert ourselves into a daemon.
  88.     We enter the phone booth as mild-mannered Clark Kent,
  89.     we make a copy of ourselves, the original disappears,
  90.     the child emerges and then ....
  91. */
  92.     signal(SIGTTOU,SIG_IGN);    /*the keyboard can't stop us*/
  93.     signal(SIGTTIN,SIG_IGN);
  94.     signal(SIGTSTP,SIG_IGN);
  95.     if( (i=fork()) < 0)  abort();   /*fork of child failed*/
  96.     if(i > 0) 
  97.        {
  98.        printf("\n pid of process is %d",i);
  99.        exit(0);        /*the parent exits */
  100.     }
  101. /*
  102.     make ourselves into a unique process group so that we
  103.     can't re-acquire a control terminal that might be able
  104.     to stop us.
  105.     Then close our control terminal.
  106.     Nothing can stop us now ...
  107. */
  108.     if( setpgrp(0,getpid() ) == -1) abort();
  109.     if( (i= open("/dev/tty",O_RDWR) ) >= 0)
  110.        {
  111.        ioctl(i,TIOCNOTTY,(char *) NULL);    /*bye bye terminal*/
  112.        close(i);
  113.     }
  114.     for(i=0; i< NOFILE; i++) close(i);    /*close all files*/
  115.     errno=0;            /*reset errors, if any*/
  116.     chdir(rootdir);            /*set default directory*/
  117. /*
  118.     open parameter file and get system parameters
  119. */
  120.     if( (pop=fopen(pname,"r")) == NULL)    /*opened readonly*/
  121.        lockerr("Cannot open param.dat in lockclock.");
  122.     if(fscanf(pop," %f %d %d %d %f %d",¶m.xsigma,¶m.flags,
  123.        ¶m.wait0,¶m.iavgt,¶m.xoffst,&slew) !=6)
  124.        {
  125.        fclose(pop);
  126.        lockerr(" error reading parm.dat in lockclock.");
  127.     }
  128.     fclose(pop);
  129. /*
  130.     initialize leap second coming flag to off
  131. */
  132.     lsdue=0;
  133. /*
  134.     If David Mills' kernel routines are present and will be
  135.     used to set leap seconds then make sure that they are
  136.     fully disabled now.
  137. */
  138.     if(  (param.flags & 16) != 0)
  139.        {
  140.        tmx.modes=MOD_STATUS;        /*change status variable*/
  141.        tmx.status=STA_UNSYNC;        /*we start out un synched*/
  142.        syscall(SYS_ntp_adjtime,&tmx);    /*tell it to the kernel */
  143.     }
  144. /*
  145.     open system state file and get previous state
  146. */
  147.     if( (lop=fopen(sname,"r")) == NULL)    /*open read only */
  148.        lockerr("Cannot open file lock.dat in lockclock.");
  149.     for(i=1; i< NUMSTATES; i++) 
  150.        {
  151.        getlock(&lop,&lock[i]);  /*read structure from file via lop*/
  152.        if(lock[i].day == 0) break;   /*EOF on input */
  153.     }
  154.     if(fclose(lop) == EOF)
  155.         lockerr("Error closing lock.dat at startup");
  156.     while(i < NUMSTATES) lock[i++].day = 0;  /*zero remaining records*/
  157. /*
  158.     if the previous mode read from the file is 4, 5, or 6, then
  159.     set the current mode to 4, which will adjust the time if
  160.     necessary and then switch back to the locked modes after
  161.     the time is adjusted.  if the previous mode was anything
  162.     else, then start over again with mode 0
  163. */
  164.     if(lock[1].mode >= 4)
  165.        {
  166.        lock[0].mode=4;
  167.        lock[0].y=lock[1].y;
  168.        lock[0].d=lock[1].d;
  169.        lock[0].flags=lock[1].flags;
  170.        lock[0].waitn=lock[1].waitn=param.wait0;
  171.        lock[0].waitx=lock[1].waitx;
  172.        lock[0].errx=lock[1].errx;
  173.        lock[0].ybar=lock[1].ybar;
  174.        lock[0].sumerr=lock[1].sumerr;
  175.     }
  176.     else lock[0].mode=0;
  177. loop:            /*come back here when we wake up each time*/
  178. /*
  179.     now get the current measurement of the time difference.
  180.     on return, daym and fdaym are the day and time of the
  181.     measurement and delta is the time difference.
  182.     daym is the MJD as a long int;
  183.     fdaym is the time since midnight in seconds as a long int;
  184.     delta is the time difference in s as a double;
  185.       if delta > 0, the local clock is fast.
  186.  
  187.     the next parameter is the rejection threshold for the 
  188.     time difference data.  The subroutine computes three
  189.     consecutive time differences which must differ among
  190.     themselves by less than this value.  See getdif
  191.     for what happens if they don't.
  192.     
  193.     This rejection criterion is set more or less arbitrarily
  194.     at 3 X the measurement noise estimate read from the
  195.     parameter file.  If the measurement noise is white and
  196.     if param.xsigma is a good estimate of its std. deviation,
  197.     then the rejection criterion should operate about 1% of
  198.     the time or less.
  199.  
  200.     The next parameter (lsflag) is the leap second flag as read from
  201.     the ACTS transmissions
  202.  
  203.     the next paramater (sigmam) is an estimate of the measurement
  204.     noise found by getdif.  It is related to the variation among
  205.     the successive time differences found by getdif.
  206. */
  207.     j=getdif(&daym,&fdaym,&delta,(double)(3*param.xsigma),&lsflag,
  208.         &sigmam);
  209. /*
  210.     if j == 1 then the three time difference measurements agreed
  211.     among themselves to within 3 sigma.  Use the returned value
  212.     of sigmam to adjust the value of sigma from the value first
  213.     read from the file.
  214.     if j > 1 then the time difference data were too noisy based
  215.     on the current value of sigma which either means that one
  216.     of the points had a glitch or that sigma is a bit too small.
  217.     The returned value of sigmam will be larger than sigma in
  218.     this case so that sigma will grow.
  219.  
  220. */
  221.     if(j >= 1) param.xsigma= 0.8*param.xsigma + 0.2*sigmam;
  222. /*
  223.     if j<=0 then getdif had a problem connecting to acts
  224.     and the return values are not valid. wait a bit and
  225.     try again.
  226.  
  227.     if j ==3 then getdif connected to ACTS but the three
  228.     time-difference estimates disagree with each other
  229.     and there is something very wrong with the ACTS system
  230.     or with the local clock.
  231. */
  232.     if( (j <= 0)  || (j == 3) )
  233.        {
  234.        conserr++;    /*increment consecutive error counter */
  235.        if(conserr < 5)    /*wait 15 sec and retry up to 5 times*/
  236.           {
  237.           sleep(15);
  238.           goto loop;
  239.        }
  240.        else        /*retries don't fix the problem */
  241.           {        /*copy previous state to this one and wait*/
  242.           lock[0].day=lock[1].day;
  243.           lock[0].fday=lock[1].fday;
  244.           lock[0].mode=lock[1].mode;
  245.           lock[0].x=0;
  246.           lock[0].y=lock[1].y;
  247.           lock[0].d=lock[1].d;
  248.           lock[0].flags=j;
  249.           iwait=lock[0].waitn=lock[1].waitn;
  250.           lock[0].waitx=lock[1].waitx;
  251.           lock[0].errx=lock[1].errx;
  252.           lock[0].ybar=lock[1].ybar;
  253.           lock[0].sumerr=lock[1].sumerr;
  254.           goto finis;        /*hope for better times */
  255.        }
  256.     }        /*end of negative return status from getdif*/
  257. /*
  258.     test the returned leap-second status and compare it
  259.     to the value stored in lsdue to see if the status
  260.     this time is different from the last value that we
  261.     knew about.
  262.     if the status just received is > than the stored status
  263.     then either this is the first notification for an event
  264.     that is still a month away or the leap second is today.
  265.  
  266.     if the status just received is < than the stored status
  267.     then either an event was scheduled and is now over or
  268.     lsflag < 0 which means that the consecutive values of
  269.     this flag did not agree as received from ACTS. Either
  270.     situation should not happen normally since it implies 
  271.     either that we slept through the leap second,  or the
  272.     consecutive values read from ACTS did not agree.  Both
  273.     of these situations are errors. 
  274.  
  275.     if the current and stored status values are the same then
  276.     nothing happens here.
  277. */
  278.     if(lsflag > lsdue)
  279.        {
  280.        if(lsdue == 0) putleap(lsflag);    /*first notification */
  281.        lsdue=lsflag;
  282. /*
  283.        if the new status is >= 10 then this is a notice that
  284.        the leap second will happen today.  if David Mills'
  285.        kernel routines have been installed and enabled, then
  286.        leap seconds are handled by letting the kernel routine 
  287.        do it.  This leap second mode will be enabled by setting
  288.        the 5th bit of the flag word
  289. */
  290.        if( (lsdue >= 10)         /*leap second today */
  291.         && ( (param.flags & 16) != 0) )    /*use kernel routines*/
  292.           {
  293.           tmx.modes= MOD_STATUS;    /*set leap second flag*/
  294.           if(lsdue == 10) tmx.status=STA_INS;  /*insert leap second */
  295.           else            tmx.status=STA_DEL;  /*delete leap second */
  296.           syscall(SYS_ntp_adjtime,&tmx);
  297.        }
  298.     }
  299.     else if( (lsflag >= 0) &&         /*no error in lsflag*/
  300.              (lsflag < lsdue) ) lsdue=0;    /*event time is over*/
  301. /*
  302.     continue if return status from getdif is +1 or +2:
  303.     store this time difference in the state structure.
  304.  
  305.     The time difference value returned from getdif is
  306.     adjusted by subtracting a possible time offset read
  307.     from the parameter file.  This offset could be used
  308.     to correct for some systematic bias but is used on
  309.     the Sun version to move the time away from 0 to avoid
  310.     a possible race condition in the driver where the driver
  311.     latency varies as the time moves by small amounts near
  312.     a tick interrupt.  This race condition results in the
  313.     time jumping back and forth by 1 tick as the actual
  314.     difference either comes just before or just after a 
  315.     tick.
  316. */
  317.     lock[0].day=daym;    /*day in MJD */
  318.     lock[0].fday=fdaym;    /*time in seconds since 00:00*/
  319.     lock[0].x=delta - param.xoffst;    /* corrected time difference*/
  320.     lock[0].flags=j;    /*status of comparison*/
  321. /*
  322.     mode 0 is used for a totally cold start.  store
  323.     the time difference in the state structure and
  324.     advance to mode 1 for the next cycle.
  325.     the delay between cycles is set to the initial
  326.     value read from the parameter structure
  327. */
  328.     if(lock[1].mode == 0)    /*mode 0 is cold start */
  329.        {
  330.        lock[0].y=lock[0].d=lock[0].errx=lock[0].ybar=
  331.            lock[0].sumerr=0;
  332.        lock[0].mode=1;
  333.        iwait=lock[0].waitn=param.wait0;
  334.        goto finis;
  335.     }        /*end of mode 0*/
  336. /*
  337.     if this is not the first time of a cold start then
  338.     compute the actual time difference since the last
  339.     calibration in seconds.
  340.     xavgt is the time constant for frequency averaging
  341.     normalized by the current delay.  It is used in
  342.     the exponential filter for the frequency estimator.
  343.  
  344.     if the specified averaging time is longer than
  345.     3 samples times, then it is limited to 3 sample
  346.     times following a cold start to speed up the
  347.     time needed to get to equilibrium following a
  348.     re-start of the program.
  349. */
  350.     lock[0].waitx=86400*(lock[0].day - lock[1].day) +
  351.         (lock[0].fday - lock[1].fday);
  352.     xavgt=lock[0].waitx/( (double) param.iavgt);
  353.     if( (nloops < 10) && (xavgt < 0.33)) xavgt=0.33;
  354. /*
  355.     mode 1 is used to make an initial estimate for the
  356.     frequency during a cold start.  it is the mode used
  357.     for the second cycle during the cold start process
  358.  
  359.     if the time difference since the last loop is not 
  360.     much bigger than the estimate of the measurement noise
  361.     xsigma, then we did not wait long enough given the 
  362.     rate of the local clock. increase the wait time in
  363.     this case and do not advance to mode 2.  if this 
  364.     happens a second time then the rate is really very
  365.     small and it may take forever to lock up.  the true
  366.     rate is almost 0 in this case.
  367. */
  368.     if(lock[1].mode == 1)    /*make initial frequency estimate*/
  369.        {
  370.        lock[0].y=(lock[0].x - lock[1].x)/lock[0].waitx;
  371.        lock[0].d=lock[0].errx=lock[0].ybar=lock[0].sumerr=0;
  372.        if(fabs(lock[0].x - lock[1].x) < 2*param.xsigma)
  373.           {
  374.           if(lock[2].mode == 1)    /*if this is second time*/
  375.         {
  376.         lock[0].y=(lock[0].y + lock[1].y)/2;
  377.         lock[0].mode=2;
  378.         iwait=lock[0].waitn=lock[1].waitn;
  379.           }
  380.           else        /*this is only the first time here*/
  381.         {
  382.             lock[0].mode=1;
  383.             iwait=lock[0].waitn=1.5*lock[1].waitn;
  384.           }
  385.        }        /*end of first difference small */
  386.        else
  387.          {
  388.          lock[0].mode=2;
  389.          iwait=lock[0].waitn=lock[1].waitn;
  390.        }        /*end of first difference not small*/
  391.        goto finis;
  392.     }    /*end of mode 1 */
  393. /*
  394.     mode 2 is used for the third cycle of a cold start.
  395.     we make a second frequency estimate and see if it
  396.     is consistent with the previous one by checking 
  397.     how the first one would have worked as a predictor
  398.     for this round.
  399.     if the prediction error is large then either the
  400.     clock has lots of frequency noise, the channel has
  401.     lots of time noise or the drift is appreciable.
  402.     these questions should be addressed by measuring
  403.     the Allan Variance of the free running system.
  404.     the program tries to fix this problem by shortening
  405.     the interval between measurements.  Depending on what
  406.     is going on, this may not be the optimum strategy.
  407.     The initial estimate of ybar is simply the average
  408.     of the two frequency estimates that we have made
  409.     so far.
  410. */
  411.     if(lock[1].mode == 2)    /*check frequency estimate*/
  412.        {
  413.        lock[0].y=(lock[0].x - lock[1].x)/lock[0].waitx;
  414.        lock[0].d=lock[0].sumerr=0;
  415.        lock[0].ybar=(lock[1].y + lock[0].y)/2;
  416.        lock[0].errx=(lock[1].y-lock[0].y)*lock[0].waitx;
  417.        if(fabs(lock[0].errx) > 3*param.xsigma)
  418.          {
  419.          iwait=lock[0].waitn=0.9*(float)lock[1].waitn;
  420.          lock[0].mode=2;
  421.          goto finis;
  422.       }        /*end of large prediction error */
  423.       else
  424.        {
  425.        iwait=lock[0].waitn=lock[1].waitn;
  426.        lock[0].mode=3;
  427.        goto finis;
  428.       }        /*end of small prediction error */
  429.     }        /* end of mode 2*/
  430. /*
  431.     in mode 3, the exponential frequency estimate ybar
  432.     is used to predict the time difference between this
  433.     cycle and the last one.  The averaging time for this
  434.     filter is set by xavgt, which is determined from the
  435.     free-running Allan variance measurements.  The prediction 
  436.     error is stored in errx and integrated in sumerr.
  437.  
  438.     if the least significant bit of param.flags is set
  439.     the advance to mode 4 and set the local time after
  440.     4 cycles in mode 3.  If the least significant bit of
  441.     param.flags is not set, then stay in mode 3 forever.
  442.     This is the last mode that can run non-privileged and
  443.     is used for characterizing the local clock.
  444. */
  445.     if(lock[1].mode == 3)    /*begin using avg. freq for prediction*/
  446.        {
  447.        lock[0].y=(lock[0].x - lock[1].x)/lock[0].waitx;
  448.        lock[0].d=0;
  449.        lock[0].ybar=(lock[1].ybar + xavgt*lock[0].y)/(1+xavgt);
  450.        lock[0].errx=(lock[0].x - lock[1].x) 
  451.         - lock[0].waitx*lock[1].ybar;
  452.        lock[0].sumerr=lock[1].sumerr + lock[0].errx;
  453.        iwait=lock[0].waitn=lock[1].waitn;
  454.        lock[0].mode=3;
  455.        if( (param.flags & 1) == 0)
  456.                   goto finis;    /*never advance to mode 4*/
  457.        if( (lock[2].mode != 3) ||
  458.            (lock[3].mode != 3) ) goto finis;    /*wait longer*/
  459.        lock[0].mode=4;    /* advance immediately to set mode */
  460.        lock[0].sumerr=0;    /* start with clean slate */
  461.        lock[0].errx=fabs(lock[0].errx); /*value must be pos. def.*/
  462.     }        /*end of mode 3*/
  463. /*
  464.     mode 4 is used to set the clock.  If the third bit of
  465.     param.flags is set, then the clock is set to within 1
  466.     second in a single step using subroutine sett.  If the
  467.     third bit of param.flags is not set, then the clock is
  468.     adjusted gradually using adjtime. Since subroutine 
  469.     adjtime makes time adjustments slowly it may take
  470.     an appreciable time to set the clock correctly if it is 
  471.     far off.  
  472.     the integer parameter slew (defined in the header)
  473.     is the number of seconds that are needed to adjust
  474.     the local clock by 1 second.  The absolute value of
  475.     iwait/slew is therefore the maximum adjustment that
  476.     is possible in the nominal delay of iwait seconds.
  477.  
  478.     note that all of the parameters are frozen during
  479.     this time adjustment.
  480.  
  481.     when the time offset has been reduced so that it
  482.     can be done in one step, then perform this adjustment
  483.     and calculate how long it will take.  Set the delay
  484.     for this interval and switch to mode 5, which is not
  485.     a real mode, but just indicates that the final time
  486.     adjustment is in progress.  This mode is automatically
  487.     converted to mode 6 on the next round.
  488. */
  489.     if( (lock[1].mode == 4) || (lock[0].mode == 4) )
  490.        {
  491.        if( (param.flags & 4) != 0)        /*then set time in 1 step*/
  492.           {
  493.           if(fabs(lock[0].x) >= 1)        /*step time if big */
  494.         {
  495.         tcor= -rint(lock[0].x);        /* set to nearest second*/
  496.         sett(tcor);
  497.         lock[0].waitn=lock[1].waitn;
  498.         iwait=120;
  499.         goto finis;            /*wait for a bit*/
  500.           }
  501.           else goto fnladj;        /* error is less than 1 s */
  502.        }        /* end of third bit not zero */
  503. /*
  504.        if we are still here then the time will be set by
  505.        gradual adjustments rather than in a single step
  506. */
  507.        iwait=lock[0].waitn=lock[1].waitn;
  508.        if(lock[0].x <= -lock[0].waitn/slew)   /*error is large and neg.*/
  509.         {
  510.         tcor=lock[0].waitn/slew;
  511.         adjt(tcor);
  512.         goto finis;
  513.        }
  514.        if(lock[0].x >= lock[0].waitn/slew)   /*error is large and pos.*/
  515.         {
  516.         tcor= -lock[0].waitn/slew;
  517.         adjt(tcor);
  518.         goto finis;
  519.        }
  520. /*
  521.        if we get here, then the residual error can be 
  522.        removed in one more step.  perform this adjustment
  523.        and wait until it is finished, but not less than 10
  524.        minutes.  switch to the lock mode on the next round.
  525.        this intermediate step is called mode 5, but it is
  526.        not a true mode, since it will be converted to mode
  527.        6 on the next round.
  528.        note that in addition to the time adjustment remaining
  529.        to be done, we must add the amount by which the clock
  530.        will drift due to its frequency offset during the
  531.        time we will be waiting.
  532. */
  533. fnladj:                /*final adjustment is always a slew*/
  534.        tcor= -lock[0].x;
  535.        iwait=fabs(tcor)*slew;
  536.        if(iwait < 600) iwait = 600;
  537.        tcor -= iwait*lock[0].ybar;
  538.        adjt(tcor);
  539.        fstep=0;        /* reset freq. step detector*/
  540.        lock[0].mode=5;        /*see note above*/
  541. /*
  542.        if David Mills kernel stuff is enabled, tell it that
  543.        the clock is now okay.
  544. */
  545.        if( (param.flags & 16) != 0)
  546.           {
  547.           tmx.modes=MOD_STATUS | MOD_MAXERROR | MOD_ESTERROR;
  548.           tmx.status=0;        /*everything okay, pll off*/
  549.           tmx.maxerror=tmx.esterror=0;
  550.           syscall(SYS_ntp_adjtime,&tmx);
  551.        }
  552.        goto finis;
  553.     }        /* end of mode 4*/
  554. /*
  555.     mode 6 is used to keep the local clock on time.  the
  556.     average time offset should be 0 and any residual time
  557.     therefore means that the rate is not quite right.
  558.     this mode assumes that any time error last time was
  559.     removed then and that we have been making periodic
  560.     adjustments to compensate for the average frequency
  561.     offset.  Any time difference this time is therefore an
  562.     indication that the average frequency is wrong.
  563.  
  564.     There are two exceptions to this.  if the time error is 
  565.     very small (currently if it is less than xsigma/2 in 
  566.     absolute value), then it is treated as noise that is
  567.     not significantly different from 0 and is not used to 
  568.     update the frequency estimator here and is not removed 
  569.     below.  This will produce a step in what happens when
  570.     the time error is near xsigma/2 -- values slightly 
  571.     smaller are ignored and those slightly larger are used
  572.     at full strength.  
  573.  
  574.     The second exception is if the time error is very large
  575.     (currently if it is greater than 4 errx).  This is 
  576.     treated as a time step error and the response is a reset --
  577.     the frequency estimator is not changed since the time
  578.     is treated as a glitch, but the time error will be removed 
  579.     below. if the problem really is a time step, then this was
  580.     the right thing to do, but it won't work if there has been
  581.     a frequency step.  from observation, these are quite rare
  582.     although they sometimes happen on baldwin.  the parameter
  583.     fstep is set to the approximate value of the time constant
  584.     for frequency updates, and the reset algorithm is suspended
  585.     for that number of following cycles so that the update of
  586.     the frequency will be allowed to proceed.  this will allow
  587.     the frequency to adjust in case of a frequency step.
  588.  
  589.     note that the test for mode 5 will only be satisfied
  590.     once since it is the transition mode from 4.
  591. */
  592.     if( (lock[1].mode == 5) || (lock[1].mode == 6) )
  593.        {
  594.           lock[0].y=lock[1].ybar +lock[0].x/lock[0].waitx;
  595.        lock[0].d=0;
  596.        if(fabs(lock[0].x) <= 4.0*lock[1].errx)
  597.         fstep=0;    /*error is small -- no reset*/
  598.        else
  599.         {
  600.         if(fstep <= 0)    /*first time it is large*/
  601. /*
  602.                  initialize fstep to time constant for 
  603.             frequency updates and set fsteplim to this 
  604.             same value.  This is used below to show that 
  605.             this is the first round of a reset
  606. */
  607.           fstep = fsteplim= 1/xavgt + 1;
  608.         else
  609.           fstep--;    /*count number of consec. resets*/
  610.       }        /* end of error is large */
  611. /*
  612.       update frequency if error is not too small to be noise
  613.       and not too large to be reset or if this is not first
  614.       reset so that we are going to assume it is a frequency
  615.       step rather than a time step
  616. */
  617.       if( ( (fabs(lock[0].x) >= 0.5*param.xsigma) &&
  618.             (fabs(lock[0].x) <= 4.0*lock[1].errx)    ) ||
  619.           ( (fstep > 0) && (fstep < fsteplim)        )   )
  620.            lock[0].ybar=(lock[1].ybar + xavgt*lock[0].y)/(1+xavgt);
  621.        else  lock[0].ybar=lock[1].ybar;
  622. /*
  623.        errx is the average prediction error over the last
  624.        10 cycles and sumerr is simply the integral of the error.
  625. */
  626.        lock[0].errx=0.9*lock[1].errx + 0.1*fabs(lock[0].x);
  627.        lock[0].sumerr=lock[1].sumerr + lock[0].x;
  628.        iwait=lock[0].waitn=lock[1].waitn;
  629. /*
  630.        if the 4th bit of param.flags is set then switch
  631.        to mode 7 after 10 rounds.  this will synchronize
  632.        the system using the hardware tick via the CTS line
  633.        rather than the data from ACTS.  revert back to
  634.        mode 6 if this system doesn't work.
  635. */
  636.        if( ( (param.flags & 8) != 0) && (nloops >= 10) ) 
  637.           {
  638.           lock[0].mode=7;
  639.           conterr= 0;
  640.        }
  641.        else  lock[0].mode=6;
  642.        goto finis;
  643.     }    /* end of mode 5/6 */
  644. /*
  645.     mode 7 is used if the local machine has a direct 1 pps
  646.     input via the CTS line on a serial port.  It is like mode
  647.     6 in that the local clock is continuously adjusted, but
  648.     the adjustments use the external tick and are set to whatever
  649.     it takes to bring the local clock in coincidence with the
  650.     hardware tick.
  651.  
  652.     the periodic acts connections in this situation are used
  653.     only to be sure that we have not lost a whole second since
  654.     the hardware system cannot detect this.
  655.  
  656.     the frequency and average frequency are simply propagated
  657.     across without change in this mode.  Parameters errx and
  658.     sumerr are handled the same as in mode 6.
  659.  
  660.     if the error is less than 50 msec, then everything is
  661.     probably okay.  if it is larger than this value, then
  662.     something has gone wrong. switch back to mode 6 to try
  663.     and see what is happening. Also reset nloops to 0 so
  664.     that the program will be forced to stay in mode 6
  665.     for at least 10 loops to reacquire the frequency
  666.     if it has changed.
  667. */
  668.     if(lock[1].mode == 7)
  669.        {
  670.        lock[0].y=lock[1].y;
  671.        lock[0].ybar=lock[1].ybar;
  672.        lock[0].d=0;
  673.        lock[0].errx= 0.9*lock[1].errx + 0.1*fabs(lock[0].x);
  674.        lock[0].sumerr=lock[1].sumerr + lock[0].x;
  675.        iwait=lock[0].waitn= lock[1].waitn;
  676.        if(fabs(lock[0].x) <  0.05) lock[0].mode=7;
  677.        else 
  678.          {
  679.          nloops=0;
  680.          lock[0].mode=6;
  681.        }
  682.        goto  finis;        /*finished with mode 7*/
  683.     }        /* end of lock mode 7*/
  684. finis:        /*finished with this round*/
  685.     conserr = 0;        /*reset consecutive error counter*/
  686. /*
  687.     if the nominal delay is longer than 3000 seconds, then
  688.     limit the delay to 3000 seconds during the first 10 loops
  689.     following a re-start so that the program settles down
  690.     pretty quickly.  if iwait is already < 3000 then it is
  691.     short enough and need not be reduced.
  692. */
  693.     if( (nloops++ < 10) && (iwait > 3000) ) iwait = 3000;
  694.     if( (lop=fopen(sname,"r+")) == NULL) 
  695.        lockerr("Cannot open lock.dat for update in lockclock.");
  696.     for(i=0; i< NUMSTATES; i++) putlock(&lop,&lock[i]);
  697.     if(fclose(lop) == EOF)
  698.         lockerr("Error close lock.dat after update.");
  699.     pushlock(lock,0,NUMSTATES);    /*push down lock stack by 1*/
  700. /*
  701.     if generate new files names has been selected and if this
  702.     file is getting full, then copy it to anoter name.
  703.     the new name will be the same as the root name as stored
  704.     in string sname with the latest MJD added on to the end
  705.     to make it unique.
  706. */
  707.     if( (param.flags & 2) != 0)    /*generate new files enabled*/
  708.        {
  709.        if(filefull-- < 8)        /*file is almost full*/
  710.           {
  711. #if ALPHA == 1
  712.           sprintf(filenam,"cp %s %s_%d",sname,sname,lock[0].day);
  713. #else
  714.           sprintf(filenam,"cp %s %s_%ld",sname,sname,lock[0].day);
  715. #endif
  716.           system(filenam);
  717.           pid=wait(status);        /*wait for copy to finish*/
  718.           filefull = NUMSTATES;    /*reset for next round */
  719.        }
  720.     }
  721.     if(argc > 1)        /*store time of next update*/
  722.        {
  723.        i=(fdaym+iwait)/3600;    /*compute hr and min. of next cycle*/
  724.        j=((fdaym+iwait)%3600)/60;
  725.        sprintf(msg,"-> %d:%.2d >%.1f",i,j,1000*param.xsigma);
  726. /*
  727.        if a leap second is due today then announce it in the
  728.        argument list
  729. */
  730.        if(lsdue >= 10) msg[1]='!';
  731.        s=argv[1];
  732.        for(i=0; msg[i] != 0; i++) *s++ = msg[i];
  733.        *s='\0';            /*add trailing 0 byte */
  734.     }
  735.     if(lock[0].mode == 6)        /*this is continuous adjust mode*/
  736.        {
  737. /*
  738.        first remove the static time error if it is >= than
  739.        the measurement noise estimate.  if it is less than
  740.        the measurement noise then don't do the correction since
  741.        that is just as likely to make things worse.
  742. */
  743.        if(fabs(lock[0].x) >= 0.5*param.xsigma)
  744.           {
  745.           tcor= -lock[0].x;
  746.           jwait= (int)(fabs(tcor)*slew) + 1;
  747.           adjt(tcor);
  748.           sleep(jwait);
  749.           fdaym += jwait;
  750.        }
  751.     }        /* end of static time correction in mode 6*/
  752. /*
  753.     make a static adjustment in mode 7 if we have slipped an
  754.     entire second.
  755. */
  756.     if(  (lock[0].mode    == 7) &&
  757.          (fabs(lock[0].x) >= 1) )
  758.         {
  759.         tcor= -lock[0].x;
  760.         jwait= (int)(fabs(tcor)*slew) +1;
  761.         adjt(tcor);
  762.         sleep(jwait);
  763.         fdaym += jwait;
  764.     }
  765.     if( (lock[0].mode == 6) || (lock[0].mode == 7) )
  766.        {
  767. /*
  768.        if a leap second adjustment is due today and if it is
  769.        near the end of the day, then do it now.  Since it 
  770.        takes slew seconds to make the adjustment, if we are
  771.        within about slew seconds of the end of the day, then
  772.        we might just as well do it now.
  773.  
  774.        if bit 4 of the flag word is set then leap second
  775.        adjustments are made by calling David Mills' kernel
  776.        routines and no leap second adjustments are ever
  777.        made here.
  778. */
  779.        if( ( (param.flags & 16) == 0) && (lsdue >= 10) ) /*if due today */
  780.           {
  781.           if( (fdaym + slew) > 86200l)    /*do it now*/
  782.         {
  783.          if(lsdue == 10) lptcor= -1;    /* ls = 1 */
  784.         else            lptcor=  1;    /* ls = 2 */
  785.         adjt(lptcor);
  786.         sleep(slew);
  787.         fdaym += slew;
  788.         lsdue= 0;
  789.           }        /*end of do correction now*/
  790.        }        /* end of currection due today*/
  791.     }        /* end of mode is 6 or 7 */
  792.     if(lock[0].mode == 6)    /*begin periodic corrections in mode 6*/
  793.        {
  794. /*
  795.        find out how long we must wait before a correction
  796.        of 1 ms will be needed.  Then adjust this interval
  797.        so that there will be an integral number of them in
  798.        the calibration delay and then compute the exact
  799.        correction needed each time.
  800.  
  801.        if the offset frequency is very small, then the 
  802.        correction needed in the full iwait seconds may
  803.        be less than 1 ms.  Compute the correction needed
  804.        in this case and apply it in two steps.
  805.        
  806.        if the offset frequency is very large, then the
  807.        delay between 1 ms corrections may be less than
  808.        1 sec.  Since the maximum slew rate is on the order
  809.        of milliseconds/second, the best we can do is to
  810.        delay for 1 second and then adjust at the maximum
  811.        slew rate.  
  812. */
  813.        jwait=(LONG) fabs(0.001/lock[0].ybar);
  814.        if(jwait < 1) jwait = 1;
  815.        if(jwait > iwait) 
  816.           {
  817.           jwait=iwait/2;
  818.           jloop=2;
  819.        }            /*end of jwait too big */
  820.        else            /*jwait is not too big */
  821.          {
  822.          jloop=iwait/jwait;    /*number of adjustments per cycle*/
  823.          jwait=iwait/jloop;    /*actual delay between adjustments*/
  824.        }
  825.        tcor= -((double)jwait)*lock[0].ybar;
  826.        for(i=0; i<jloop; i++)
  827.           {
  828. /*
  829.         if a leap second is due today then do it when
  830.         we get within slew seconds of the end of the
  831.         day, since it takes slew seconds to do it. we
  832.         will make the correction on the last sleep cycle
  833.         before the end of the day.
  834. */
  835.           if( ( (param.flags & 16) == 0) /*leap sec via slew*/
  836.             && (lsdue >= 10) &&     /*leap second today*/
  837.           (fdaym + slew + jwait >=86400l) ) /*next cycle is too late*/
  838.         {
  839.         if(lsdue == 10) lptcor = -1;
  840.         else            lptcor =  1;
  841.         adjt(lptcor);            /*leap second adjustment*/
  842.         sleep(slew);
  843.         lsdue = 0;
  844.           }
  845.           adjt(tcor);
  846.           sleep(jwait);
  847.           fdaym += jwait;
  848.        }        /*end of periodic adjustment loop */
  849.        goto loop;    /*go around again */
  850.     }        /*end of incremental corrections in mode 6*/
  851.     if(lock[0].mode == 7)    /*incremental adjustments in mode 7*/
  852.        {
  853. /*
  854.        proceed as with mode 6.  find out how long between 
  855.        corrections assuming that the error is to be no more
  856.        than 0.2 msec -- about the jitter in consecutive ACTS
  857.        measurements.  then schedule a loop with this interval
  858.  
  859.        don't let the loop interval be less than 10 seconds 
  860.        and don't let there be less than 2 tests per acts call.
  861.        The acts calls are really just to be sure that we have
  862.        not jumped an entire second.
  863. */
  864.        jwait= (LONG) fabs(0.0002/lock[0].ybar);
  865.        if(jwait < 10) jwait = 10;
  866.        if(jwait > iwait/2) jwait=iwait/2;
  867.        jloop= iwait/jwait;    /*number of adjustments per ACTS cycle*/
  868.        jwait=iwait/jloop;    /*actual delay between adjustments */
  869.        for(i=0; i<jloop; i++)
  870.           {
  871.           if( ( (param.flags & 16) == 0) &&    /*leap sec via slew*/
  872.             (lsdue >= 10)     &&    /*if leap second today */
  873.           (fdaym + slew +jwait >= 86400l) )
  874.             {
  875.             if(lsdue == 10) lptcor = -1;
  876.             else        lptcor =  1;
  877.             adjt(lptcor);
  878.             sleep(slew);
  879.             lsdue = 0;
  880.           }
  881. /*
  882.           gettck returns microseconds portion of time when
  883.           tick is received.  if returned value is negative 
  884.           then system could not catch the tick.  Use the
  885.           previous correction in this case and increment
  886.           the error flag.  if we miss several in a row, then
  887.           something may be wrong with the tick system. try
  888.           to fix that by going back to mode 6.
  889. */
  890.           for(j=0; j<5; j++)
  891.         if( (jadj = gettck() ) > 0) break;    /*got it*/
  892.           if(jadj < 0)     /*5 errors in a row */
  893.         {
  894.         jadj= ladj;
  895.         if(conterr++ > 5)
  896.            {
  897.            lock[0].mode = 6;    /*go back to mode 6*/
  898.            nloops= 0;        /*may have to remeasure freq*/
  899.            goto loop;        /*get an acts calibration*/
  900.         }
  901.          }
  902.           else        /*no error from gettck */
  903.            {
  904.            conterr = 0;
  905.            ladj=jadj;    /*save this one for next time */
  906.         }
  907. /*
  908.           adjust time towards nearest second.  if microsecond
  909.           value is < 0.5 second, retard clock, otherwise advance
  910.           it to the next second.  This algorithm will fail if
  911.           the interval between calibrations is long enough that
  912.           the clock has drifted by 0.5 second or more, in which
  913.           case we will convert this error to a full second.
  914.  
  915.           if the time slew is not too large, then it will
  916.           complete in less than jwait seconds, so then wait
  917.           jwait seconds and loop again.
  918.           if the measured time error is very large so that it
  919.           will not complete in jwait seconds, then wait as long
  920.           as it will take.  This is almost certainly an error of 
  921.           some kind.  Nevertheless, it is important not to clobber 
  922.           the adjustment loop by going around too quickly.
  923. */
  924.           if(jadj < 500000l) tcor = -jadj;
  925.           else tcor= 1000000l - jadj;
  926.           tcor *= 1e-6;        /*convert usec to sec */
  927.           adjt(tcor);
  928. /*
  929.           if default loop time is greater than the time needed to
  930.           make the adjustment then all is well; set the sleep time
  931.           to the default loop time.  otherwise use a sleep time of
  932.           however long it will take based on the slew parameter
  933. */
  934.           if(jwait > 
  935.            (kwait=(int)(fabs(tcor)*slew + 1)) ) kwait=jwait;
  936.           sleep(kwait);    /* wait for next one */
  937.           fdaym += kwait;
  938.        }
  939.     }        /* end of correction in mode 7 */
  940.     else            /*mode is anything other than 6 or 7*/
  941.        {
  942.        sleep(iwait);    /*wait for next round*/
  943.     }
  944.     goto loop;        /*go around again*/
  945. }
  946.