1*2d4e511cSCy Schubert /* 2*2d4e511cSCy Schubert * ntp_calgps.c - calendar for GPS/GNSS based clocks 3*2d4e511cSCy Schubert * 4*2d4e511cSCy Schubert * Written by Juergen Perlinger (perlinger@ntp.org) for the NTP project. 5*2d4e511cSCy Schubert * The contents of 'html/copyright.html' apply. 6*2d4e511cSCy Schubert * 7*2d4e511cSCy Schubert * -------------------------------------------------------------------- 8*2d4e511cSCy Schubert * 9*2d4e511cSCy Schubert * This module implements stuff often used with GPS/GNSS receivers 10*2d4e511cSCy Schubert */ 11*2d4e511cSCy Schubert 12*2d4e511cSCy Schubert #include <config.h> 13*2d4e511cSCy Schubert #include <sys/types.h> 14*2d4e511cSCy Schubert 15*2d4e511cSCy Schubert #include "ntp_types.h" 16*2d4e511cSCy Schubert #include "ntp_calendar.h" 17*2d4e511cSCy Schubert #include "ntp_calgps.h" 18*2d4e511cSCy Schubert #include "ntp_stdlib.h" 19*2d4e511cSCy Schubert #include "ntp_unixtime.h" 20*2d4e511cSCy Schubert 21*2d4e511cSCy Schubert #include "ntp_fp.h" 22*2d4e511cSCy Schubert #include "ntpd.h" 23*2d4e511cSCy Schubert #include "vint64ops.h" 24*2d4e511cSCy Schubert 25*2d4e511cSCy Schubert /* ==================================================================== 26*2d4e511cSCy Schubert * misc. helpers -- might go elsewhere sometime? 27*2d4e511cSCy Schubert * ==================================================================== 28*2d4e511cSCy Schubert */ 29*2d4e511cSCy Schubert 30*2d4e511cSCy Schubert l_fp 31*2d4e511cSCy Schubert ntpfp_with_fudge( 32*2d4e511cSCy Schubert l_fp lfp, 33*2d4e511cSCy Schubert double ofs 34*2d4e511cSCy Schubert ) 35*2d4e511cSCy Schubert { 36*2d4e511cSCy Schubert l_fp fpo; 37*2d4e511cSCy Schubert /* calculate 'lfp - ofs' as '(l_fp)(-ofs) + lfp': negating a 38*2d4e511cSCy Schubert * double is cheap, as it only flips one bit... 39*2d4e511cSCy Schubert */ 40*2d4e511cSCy Schubert ofs = -ofs; 41*2d4e511cSCy Schubert DTOLFP(ofs, &fpo); 42*2d4e511cSCy Schubert L_ADD(&fpo, &lfp); 43*2d4e511cSCy Schubert return fpo; 44*2d4e511cSCy Schubert } 45*2d4e511cSCy Schubert 46*2d4e511cSCy Schubert 47*2d4e511cSCy Schubert /* ==================================================================== 48*2d4e511cSCy Schubert * GPS calendar functions 49*2d4e511cSCy Schubert * ==================================================================== 50*2d4e511cSCy Schubert */ 51*2d4e511cSCy Schubert 52*2d4e511cSCy Schubert /* -------------------------------------------------------------------- 53*2d4e511cSCy Schubert * normalization functions for day/time and week/time representations. 54*2d4e511cSCy Schubert * Since we only use moderate offsets (leap second corrections and 55*2d4e511cSCy Schubert * alike) it does not really pay off to do a floor-corrected division 56*2d4e511cSCy Schubert * here. We use compare/decrement/increment loops instead. 57*2d4e511cSCy Schubert * -------------------------------------------------------------------- 58*2d4e511cSCy Schubert */ 59*2d4e511cSCy Schubert static void 60*2d4e511cSCy Schubert _norm_ntp_datum( 61*2d4e511cSCy Schubert TNtpDatum * datum 62*2d4e511cSCy Schubert ) 63*2d4e511cSCy Schubert { 64*2d4e511cSCy Schubert static const int32_t limit = SECSPERDAY; 65*2d4e511cSCy Schubert 66*2d4e511cSCy Schubert if (datum->secs >= limit) { 67*2d4e511cSCy Schubert do 68*2d4e511cSCy Schubert ++datum->days; 69*2d4e511cSCy Schubert while ((datum->secs -= limit) >= limit); 70*2d4e511cSCy Schubert } else if (datum->secs < 0) { 71*2d4e511cSCy Schubert do 72*2d4e511cSCy Schubert --datum->days; 73*2d4e511cSCy Schubert while ((datum->secs += limit) < 0); 74*2d4e511cSCy Schubert } 75*2d4e511cSCy Schubert } 76*2d4e511cSCy Schubert 77*2d4e511cSCy Schubert static void 78*2d4e511cSCy Schubert _norm_gps_datum( 79*2d4e511cSCy Schubert TGpsDatum * datum 80*2d4e511cSCy Schubert ) 81*2d4e511cSCy Schubert { 82*2d4e511cSCy Schubert static const int32_t limit = 7 * SECSPERDAY; 83*2d4e511cSCy Schubert 84*2d4e511cSCy Schubert if (datum->wsecs >= limit) { 85*2d4e511cSCy Schubert do 86*2d4e511cSCy Schubert ++datum->weeks; 87*2d4e511cSCy Schubert while ((datum->wsecs -= limit) >= limit); 88*2d4e511cSCy Schubert } else if (datum->wsecs < 0) { 89*2d4e511cSCy Schubert do 90*2d4e511cSCy Schubert --datum->weeks; 91*2d4e511cSCy Schubert while ((datum->wsecs += limit) < 0); 92*2d4e511cSCy Schubert } 93*2d4e511cSCy Schubert } 94*2d4e511cSCy Schubert 95*2d4e511cSCy Schubert /* -------------------------------------------------------------------- 96*2d4e511cSCy Schubert * Add an offset to a day/time and week/time representation. 97*2d4e511cSCy Schubert * 98*2d4e511cSCy Schubert * !!Attention!! the offset should be small, compared to the time period 99*2d4e511cSCy Schubert * (either a day or a week). 100*2d4e511cSCy Schubert * -------------------------------------------------------------------- 101*2d4e511cSCy Schubert */ 102*2d4e511cSCy Schubert void 103*2d4e511cSCy Schubert gpsntp_add_offset( 104*2d4e511cSCy Schubert TNtpDatum * datum, 105*2d4e511cSCy Schubert l_fp offset 106*2d4e511cSCy Schubert ) 107*2d4e511cSCy Schubert { 108*2d4e511cSCy Schubert /* fraction can be added easily */ 109*2d4e511cSCy Schubert datum->frac += offset.l_uf; 110*2d4e511cSCy Schubert datum->secs += (datum->frac < offset.l_uf); 111*2d4e511cSCy Schubert 112*2d4e511cSCy Schubert /* avoid integer overflow on the seconds */ 113*2d4e511cSCy Schubert if (offset.l_ui >= INT32_MAX) 114*2d4e511cSCy Schubert datum->secs -= (int32_t)~offset.l_ui + 1; 115*2d4e511cSCy Schubert else 116*2d4e511cSCy Schubert datum->secs += (int32_t)offset.l_ui; 117*2d4e511cSCy Schubert _norm_ntp_datum(datum); 118*2d4e511cSCy Schubert } 119*2d4e511cSCy Schubert 120*2d4e511cSCy Schubert void 121*2d4e511cSCy Schubert gpscal_add_offset( 122*2d4e511cSCy Schubert TGpsDatum * datum, 123*2d4e511cSCy Schubert l_fp offset 124*2d4e511cSCy Schubert ) 125*2d4e511cSCy Schubert { 126*2d4e511cSCy Schubert /* fraction can be added easily */ 127*2d4e511cSCy Schubert datum->frac += offset.l_uf; 128*2d4e511cSCy Schubert datum->wsecs += (datum->frac < offset.l_uf); 129*2d4e511cSCy Schubert 130*2d4e511cSCy Schubert 131*2d4e511cSCy Schubert /* avoid integer overflow on the seconds */ 132*2d4e511cSCy Schubert if (offset.l_ui >= INT32_MAX) 133*2d4e511cSCy Schubert datum->wsecs -= (int32_t)~offset.l_ui + 1; 134*2d4e511cSCy Schubert else 135*2d4e511cSCy Schubert datum->wsecs += (int32_t)offset.l_ui; 136*2d4e511cSCy Schubert _norm_gps_datum(datum); 137*2d4e511cSCy Schubert } 138*2d4e511cSCy Schubert 139*2d4e511cSCy Schubert /* ------------------------------------------------------------------- 140*2d4e511cSCy Schubert * API functions civil calendar and NTP datum 141*2d4e511cSCy Schubert * ------------------------------------------------------------------- 142*2d4e511cSCy Schubert */ 143*2d4e511cSCy Schubert 144*2d4e511cSCy Schubert static TNtpDatum 145*2d4e511cSCy Schubert _gpsntp_fix_gps_era( 146*2d4e511cSCy Schubert TcNtpDatum * in 147*2d4e511cSCy Schubert ) 148*2d4e511cSCy Schubert { 149*2d4e511cSCy Schubert /* force result in basedate era 150*2d4e511cSCy Schubert * 151*2d4e511cSCy Schubert * When calculating this directly in days, we have to execute a 152*2d4e511cSCy Schubert * real modulus calculation, since we're obviously not doing a 153*2d4e511cSCy Schubert * modulus by a power of 2. Executing this as true floor mod 154*2d4e511cSCy Schubert * needs some care and is done under explicit usage of one's 155*2d4e511cSCy Schubert * complement and masking to get mostly branchless code. 156*2d4e511cSCy Schubert */ 157*2d4e511cSCy Schubert static uint32_t const clen = 7*1024; 158*2d4e511cSCy Schubert 159*2d4e511cSCy Schubert uint32_t base, days, sign; 160*2d4e511cSCy Schubert TNtpDatum out = *in; 161*2d4e511cSCy Schubert 162*2d4e511cSCy Schubert /* Get base in NTP day scale. No overflows here. */ 163*2d4e511cSCy Schubert base = (basedate_get_gpsweek() + GPSNTP_WSHIFT) * 7 164*2d4e511cSCy Schubert - GPSNTP_DSHIFT; 165*2d4e511cSCy Schubert days = out.days; 166*2d4e511cSCy Schubert 167*2d4e511cSCy Schubert sign = (uint32_t)-(days < base); 168*2d4e511cSCy Schubert days = sign ^ (days - base); 169*2d4e511cSCy Schubert days %= clen; 170*2d4e511cSCy Schubert days = base + (sign & clen) + (sign ^ days); 171*2d4e511cSCy Schubert 172*2d4e511cSCy Schubert out.days = days; 173*2d4e511cSCy Schubert return out; 174*2d4e511cSCy Schubert } 175*2d4e511cSCy Schubert 176*2d4e511cSCy Schubert TNtpDatum 177*2d4e511cSCy Schubert gpsntp_fix_gps_era( 178*2d4e511cSCy Schubert TcNtpDatum * in 179*2d4e511cSCy Schubert ) 180*2d4e511cSCy Schubert { 181*2d4e511cSCy Schubert TNtpDatum out = *in; 182*2d4e511cSCy Schubert _norm_ntp_datum(&out); 183*2d4e511cSCy Schubert return _gpsntp_fix_gps_era(&out); 184*2d4e511cSCy Schubert } 185*2d4e511cSCy Schubert 186*2d4e511cSCy Schubert /* ----------------------------------------------------------------- */ 187*2d4e511cSCy Schubert static TNtpDatum 188*2d4e511cSCy Schubert _gpsntp_from_daytime( 189*2d4e511cSCy Schubert TcCivilDate * jd, 190*2d4e511cSCy Schubert l_fp fofs, 191*2d4e511cSCy Schubert TcNtpDatum * pivot, 192*2d4e511cSCy Schubert int warp 193*2d4e511cSCy Schubert ) 194*2d4e511cSCy Schubert { 195*2d4e511cSCy Schubert static const int32_t shift = SECSPERDAY / 2; 196*2d4e511cSCy Schubert 197*2d4e511cSCy Schubert TNtpDatum retv; 198*2d4e511cSCy Schubert 199*2d4e511cSCy Schubert /* set result based on pivot -- ops order is important here */ 200*2d4e511cSCy Schubert ZERO(retv); 201*2d4e511cSCy Schubert retv.secs = ntpcal_date_to_daysec(jd); 202*2d4e511cSCy Schubert gpsntp_add_offset(&retv, fofs); /* result is normalized */ 203*2d4e511cSCy Schubert retv.days = pivot->days; 204*2d4e511cSCy Schubert 205*2d4e511cSCy Schubert /* Manual periodic extension without division: */ 206*2d4e511cSCy Schubert if (pivot->secs < shift) { 207*2d4e511cSCy Schubert int32_t lim = pivot->secs + shift; 208*2d4e511cSCy Schubert retv.days -= (retv.secs > lim || 209*2d4e511cSCy Schubert (retv.secs == lim && retv.frac >= pivot->frac)); 210*2d4e511cSCy Schubert } else { 211*2d4e511cSCy Schubert int32_t lim = pivot->secs - shift; 212*2d4e511cSCy Schubert retv.days += (retv.secs < lim || 213*2d4e511cSCy Schubert (retv.secs == lim && retv.frac < pivot->frac)); 214*2d4e511cSCy Schubert } 215*2d4e511cSCy Schubert return warp ? _gpsntp_fix_gps_era(&retv) : retv; 216*2d4e511cSCy Schubert } 217*2d4e511cSCy Schubert 218*2d4e511cSCy Schubert /* ----------------------------------------------------------------- 219*2d4e511cSCy Schubert * Given the time-of-day part of a civil datum and an additional 220*2d4e511cSCy Schubert * (fractional) offset, calculate a full time stamp around a given pivot 221*2d4e511cSCy Schubert * time so that the difference between the pivot and the resulting time 222*2d4e511cSCy Schubert * stamp is less or equal to 12 hours absolute. 223*2d4e511cSCy Schubert */ 224*2d4e511cSCy Schubert TNtpDatum 225*2d4e511cSCy Schubert gpsntp_from_daytime2_ex( 226*2d4e511cSCy Schubert TcCivilDate * jd, 227*2d4e511cSCy Schubert l_fp fofs, 228*2d4e511cSCy Schubert TcNtpDatum * pivot, 229*2d4e511cSCy Schubert int/*BOOL*/ warp 230*2d4e511cSCy Schubert ) 231*2d4e511cSCy Schubert { 232*2d4e511cSCy Schubert TNtpDatum dpiv = *pivot; 233*2d4e511cSCy Schubert _norm_ntp_datum(&dpiv); 234*2d4e511cSCy Schubert return _gpsntp_from_daytime(jd, fofs, &dpiv, warp); 235*2d4e511cSCy Schubert } 236*2d4e511cSCy Schubert 237*2d4e511cSCy Schubert /* ----------------------------------------------------------------- 238*2d4e511cSCy Schubert * This works similar to 'gpsntp_from_daytime1()' and actually even uses 239*2d4e511cSCy Schubert * it, but the pivot is calculated from the pivot given as 'l_fp' in NTP 240*2d4e511cSCy Schubert * time scale. This is in turn expanded around the current system time, 241*2d4e511cSCy Schubert * and the resulting absolute pivot is then used to calculate the full 242*2d4e511cSCy Schubert * NTP time stamp. 243*2d4e511cSCy Schubert */ 244*2d4e511cSCy Schubert TNtpDatum 245*2d4e511cSCy Schubert gpsntp_from_daytime1_ex( 246*2d4e511cSCy Schubert TcCivilDate * jd, 247*2d4e511cSCy Schubert l_fp fofs, 248*2d4e511cSCy Schubert l_fp pivot, 249*2d4e511cSCy Schubert int/*BOOL*/ warp 250*2d4e511cSCy Schubert ) 251*2d4e511cSCy Schubert { 252*2d4e511cSCy Schubert vint64 pvi64; 253*2d4e511cSCy Schubert TNtpDatum dpiv; 254*2d4e511cSCy Schubert ntpcal_split split; 255*2d4e511cSCy Schubert 256*2d4e511cSCy Schubert pvi64 = ntpcal_ntp_to_ntp(pivot.l_ui, NULL); 257*2d4e511cSCy Schubert split = ntpcal_daysplit(&pvi64); 258*2d4e511cSCy Schubert dpiv.days = split.hi; 259*2d4e511cSCy Schubert dpiv.secs = split.lo; 260*2d4e511cSCy Schubert dpiv.frac = pivot.l_uf; 261*2d4e511cSCy Schubert return _gpsntp_from_daytime(jd, fofs, &dpiv, warp); 262*2d4e511cSCy Schubert } 263*2d4e511cSCy Schubert 264*2d4e511cSCy Schubert /* ----------------------------------------------------------------- 265*2d4e511cSCy Schubert * Given a calendar date, zap it into a GPS time format and then convert 266*2d4e511cSCy Schubert * that one into the NTP time scale. 267*2d4e511cSCy Schubert */ 268*2d4e511cSCy Schubert TNtpDatum 269*2d4e511cSCy Schubert gpsntp_from_calendar_ex( 270*2d4e511cSCy Schubert TcCivilDate * jd, 271*2d4e511cSCy Schubert l_fp fofs, 272*2d4e511cSCy Schubert int/*BOOL*/ warp 273*2d4e511cSCy Schubert ) 274*2d4e511cSCy Schubert { 275*2d4e511cSCy Schubert TGpsDatum gps; 276*2d4e511cSCy Schubert gps = gpscal_from_calendar_ex(jd, fofs, warp); 277*2d4e511cSCy Schubert return gpsntp_from_gpscal_ex(&gps, FALSE); 278*2d4e511cSCy Schubert } 279*2d4e511cSCy Schubert 280*2d4e511cSCy Schubert /* ----------------------------------------------------------------- 281*2d4e511cSCy Schubert * create a civil calendar datum from a NTP date representation 282*2d4e511cSCy Schubert */ 283*2d4e511cSCy Schubert void 284*2d4e511cSCy Schubert gpsntp_to_calendar( 285*2d4e511cSCy Schubert TCivilDate * cd, 286*2d4e511cSCy Schubert TcNtpDatum * nd 287*2d4e511cSCy Schubert ) 288*2d4e511cSCy Schubert { 289*2d4e511cSCy Schubert memset(cd, 0, sizeof(*cd)); 290*2d4e511cSCy Schubert ntpcal_rd_to_date( 291*2d4e511cSCy Schubert cd, 292*2d4e511cSCy Schubert nd->days + DAY_NTP_STARTS + ntpcal_daysec_to_date( 293*2d4e511cSCy Schubert cd, nd->secs)); 294*2d4e511cSCy Schubert } 295*2d4e511cSCy Schubert 296*2d4e511cSCy Schubert /* ----------------------------------------------------------------- 297*2d4e511cSCy Schubert * get day/tod representation from week/tow datum 298*2d4e511cSCy Schubert */ 299*2d4e511cSCy Schubert TNtpDatum 300*2d4e511cSCy Schubert gpsntp_from_gpscal_ex( 301*2d4e511cSCy Schubert TcGpsDatum * gd, 302*2d4e511cSCy Schubert int/*BOOL*/ warp 303*2d4e511cSCy Schubert ) 304*2d4e511cSCy Schubert { 305*2d4e511cSCy Schubert TNtpDatum retv; 306*2d4e511cSCy Schubert vint64 ts64; 307*2d4e511cSCy Schubert ntpcal_split split; 308*2d4e511cSCy Schubert TGpsDatum date = *gd; 309*2d4e511cSCy Schubert 310*2d4e511cSCy Schubert if (warp) { 311*2d4e511cSCy Schubert uint32_t base = basedate_get_gpsweek() + GPSNTP_WSHIFT; 312*2d4e511cSCy Schubert _norm_gps_datum(&date); 313*2d4e511cSCy Schubert date.weeks = ((date.weeks - base) & 1023u) + base; 314*2d4e511cSCy Schubert } 315*2d4e511cSCy Schubert 316*2d4e511cSCy Schubert ts64 = ntpcal_weekjoin(date.weeks, date.wsecs); 317*2d4e511cSCy Schubert ts64 = subv64u32(&ts64, (GPSNTP_DSHIFT * SECSPERDAY)); 318*2d4e511cSCy Schubert split = ntpcal_daysplit(&ts64); 319*2d4e511cSCy Schubert 320*2d4e511cSCy Schubert retv.frac = gd->frac; 321*2d4e511cSCy Schubert retv.secs = split.lo; 322*2d4e511cSCy Schubert retv.days = split.hi; 323*2d4e511cSCy Schubert return retv; 324*2d4e511cSCy Schubert } 325*2d4e511cSCy Schubert 326*2d4e511cSCy Schubert /* ----------------------------------------------------------------- 327*2d4e511cSCy Schubert * get LFP from ntp datum 328*2d4e511cSCy Schubert */ 329*2d4e511cSCy Schubert l_fp 330*2d4e511cSCy Schubert ntpfp_from_ntpdatum( 331*2d4e511cSCy Schubert TcNtpDatum * nd 332*2d4e511cSCy Schubert ) 333*2d4e511cSCy Schubert { 334*2d4e511cSCy Schubert l_fp retv; 335*2d4e511cSCy Schubert 336*2d4e511cSCy Schubert retv.l_uf = nd->frac; 337*2d4e511cSCy Schubert retv.l_ui = nd->days * (uint32_t)SECSPERDAY 338*2d4e511cSCy Schubert + nd->secs; 339*2d4e511cSCy Schubert return retv; 340*2d4e511cSCy Schubert } 341*2d4e511cSCy Schubert 342*2d4e511cSCy Schubert /* ------------------------------------------------------------------- 343*2d4e511cSCy Schubert * API functions GPS week calendar 344*2d4e511cSCy Schubert * 345*2d4e511cSCy Schubert * Here we use a calendar base of 1899-12-31, so the NTP epoch has 346*2d4e511cSCy Schubert * { 0, 86400.0 } in this representation. 347*2d4e511cSCy Schubert * ------------------------------------------------------------------- 348*2d4e511cSCy Schubert */ 349*2d4e511cSCy Schubert 350*2d4e511cSCy Schubert static TGpsDatum 351*2d4e511cSCy Schubert _gpscal_fix_gps_era( 352*2d4e511cSCy Schubert TcGpsDatum * in 353*2d4e511cSCy Schubert ) 354*2d4e511cSCy Schubert { 355*2d4e511cSCy Schubert /* force result in basedate era 356*2d4e511cSCy Schubert * 357*2d4e511cSCy Schubert * This is based on calculating the modulus to a power of two, 358*2d4e511cSCy Schubert * so signed integer overflow does not affect the result. Which 359*2d4e511cSCy Schubert * in turn makes for a very compact calculation... 360*2d4e511cSCy Schubert */ 361*2d4e511cSCy Schubert uint32_t base, week; 362*2d4e511cSCy Schubert TGpsDatum out = *in; 363*2d4e511cSCy Schubert 364*2d4e511cSCy Schubert week = out.weeks; 365*2d4e511cSCy Schubert base = basedate_get_gpsweek() + GPSNTP_WSHIFT; 366*2d4e511cSCy Schubert week = base + ((week - base) & (GPSWEEKS - 1)); 367*2d4e511cSCy Schubert out.weeks = week; 368*2d4e511cSCy Schubert return out; 369*2d4e511cSCy Schubert } 370*2d4e511cSCy Schubert 371*2d4e511cSCy Schubert TGpsDatum 372*2d4e511cSCy Schubert gpscal_fix_gps_era( 373*2d4e511cSCy Schubert TcGpsDatum * in 374*2d4e511cSCy Schubert ) 375*2d4e511cSCy Schubert { 376*2d4e511cSCy Schubert TGpsDatum out = *in; 377*2d4e511cSCy Schubert _norm_gps_datum(&out); 378*2d4e511cSCy Schubert return _gpscal_fix_gps_era(&out); 379*2d4e511cSCy Schubert } 380*2d4e511cSCy Schubert 381*2d4e511cSCy Schubert /* ----------------------------------------------------------------- 382*2d4e511cSCy Schubert * Given a calendar date, zap it into a GPS time format and the do a 383*2d4e511cSCy Schubert * proper era mapping in the GPS time scale, based on the GPS base date, 384*2d4e511cSCy Schubert * if so requested. 385*2d4e511cSCy Schubert * 386*2d4e511cSCy Schubert * This function also augments the century if just a 2-digit year 387*2d4e511cSCy Schubert * (0..99) is provided on input. 388*2d4e511cSCy Schubert * 389*2d4e511cSCy Schubert * This is a fail-safe against GPS receivers with an unknown starting 390*2d4e511cSCy Schubert * point for their internal calendar calculation and therefore 391*2d4e511cSCy Schubert * unpredictable (but reproducible!) rollover behavior. While there 392*2d4e511cSCy Schubert * *are* receivers that create a full date in the proper way, many 393*2d4e511cSCy Schubert * others just don't. The overall damage is minimized by simply not 394*2d4e511cSCy Schubert * trusting the era mapping of the receiver and doing the era assignment 395*2d4e511cSCy Schubert * with a configurable base date *inside* ntpd. 396*2d4e511cSCy Schubert */ 397*2d4e511cSCy Schubert TGpsDatum 398*2d4e511cSCy Schubert gpscal_from_calendar_ex( 399*2d4e511cSCy Schubert TcCivilDate * jd, 400*2d4e511cSCy Schubert l_fp fofs, 401*2d4e511cSCy Schubert int/*BOOL*/ warp 402*2d4e511cSCy Schubert ) 403*2d4e511cSCy Schubert { 404*2d4e511cSCy Schubert /* (-DAY_GPS_STARTS) (mod 7*1024) -- complement of cycle shift */ 405*2d4e511cSCy Schubert static const uint32_t s_compl_shift = 406*2d4e511cSCy Schubert (7 * 1024) - DAY_GPS_STARTS % (7 * 1024); 407*2d4e511cSCy Schubert 408*2d4e511cSCy Schubert TGpsDatum gps; 409*2d4e511cSCy Schubert TCivilDate cal; 410*2d4e511cSCy Schubert int32_t days, week; 411*2d4e511cSCy Schubert 412*2d4e511cSCy Schubert /* if needed, convert from 2-digit year to full year 413*2d4e511cSCy Schubert * !!NOTE!! works only between 1980 and 2079! 414*2d4e511cSCy Schubert */ 415*2d4e511cSCy Schubert cal = *jd; 416*2d4e511cSCy Schubert if (cal.year < 80) 417*2d4e511cSCy Schubert cal.year += 2000; 418*2d4e511cSCy Schubert else if (cal.year < 100) 419*2d4e511cSCy Schubert cal.year += 1900; 420*2d4e511cSCy Schubert 421*2d4e511cSCy Schubert /* get RDN from date, possibly adjusting the century */ 422*2d4e511cSCy Schubert again: if (cal.month && cal.monthday) { /* use Y/M/D civil date */ 423*2d4e511cSCy Schubert days = ntpcal_date_to_rd(&cal); 424*2d4e511cSCy Schubert } else { /* using Y/DoY date */ 425*2d4e511cSCy Schubert days = ntpcal_year_to_ystart(cal.year) 426*2d4e511cSCy Schubert + (int32_t)cal.yearday 427*2d4e511cSCy Schubert - 1; /* both RDN and yearday start with '1'. */ 428*2d4e511cSCy Schubert } 429*2d4e511cSCy Schubert 430*2d4e511cSCy Schubert /* Rebase to days after the GPS epoch. 'days' is positive here, 431*2d4e511cSCy Schubert * but it might be less than the GPS epoch start. Depending on 432*2d4e511cSCy Schubert * the input, we have to do different things to get the desired 433*2d4e511cSCy Schubert * result. (Since we want to remap the era anyway, we only have 434*2d4e511cSCy Schubert * to retain congruential identities....) 435*2d4e511cSCy Schubert */ 436*2d4e511cSCy Schubert 437*2d4e511cSCy Schubert if (days >= DAY_GPS_STARTS) { 438*2d4e511cSCy Schubert /* simply shift to days since GPS epoch */ 439*2d4e511cSCy Schubert days -= DAY_GPS_STARTS; 440*2d4e511cSCy Schubert } else if (jd->year < 100) { 441*2d4e511cSCy Schubert /* Two-digit year on input: add another century and 442*2d4e511cSCy Schubert * retry. This can happen only if the century expansion 443*2d4e511cSCy Schubert * yielded a date between 1980-01-01 and 1980-01-05, 444*2d4e511cSCy Schubert * both inclusive. We have at most one retry here. 445*2d4e511cSCy Schubert */ 446*2d4e511cSCy Schubert cal.year += 100; 447*2d4e511cSCy Schubert goto again; 448*2d4e511cSCy Schubert } else { 449*2d4e511cSCy Schubert /* A very bad date before the GPS epoch. There's not 450*2d4e511cSCy Schubert * much we can do, except to add the complement of 451*2d4e511cSCy Schubert * DAY_GPS_STARTS % (7 * 1024) here, that is, use a 452*2d4e511cSCy Schubert * congruential identity: Add the complement instead of 453*2d4e511cSCy Schubert * subtracting the value gives a value with the same 454*2d4e511cSCy Schubert * modulus. But of course, now we MUST to go through a 455*2d4e511cSCy Schubert * cycle fix... because the date was obviously wrong! 456*2d4e511cSCy Schubert */ 457*2d4e511cSCy Schubert warp = TRUE; 458*2d4e511cSCy Schubert days += s_compl_shift; 459*2d4e511cSCy Schubert } 460*2d4e511cSCy Schubert 461*2d4e511cSCy Schubert /* Splitting to weeks is simple now: */ 462*2d4e511cSCy Schubert week = days / 7; 463*2d4e511cSCy Schubert days -= week * 7; 464*2d4e511cSCy Schubert 465*2d4e511cSCy Schubert /* re-base on start of NTP with weeks mapped to 1024 weeks 466*2d4e511cSCy Schubert * starting with the GPS base day set in the calendar. 467*2d4e511cSCy Schubert */ 468*2d4e511cSCy Schubert gps.weeks = week + GPSNTP_WSHIFT; 469*2d4e511cSCy Schubert gps.wsecs = days * SECSPERDAY + ntpcal_date_to_daysec(&cal); 470*2d4e511cSCy Schubert gps.frac = 0; 471*2d4e511cSCy Schubert gpscal_add_offset(&gps, fofs); 472*2d4e511cSCy Schubert return warp ? _gpscal_fix_gps_era(&gps) : gps; 473*2d4e511cSCy Schubert } 474*2d4e511cSCy Schubert 475*2d4e511cSCy Schubert /* ----------------------------------------------------------------- 476*2d4e511cSCy Schubert * get civil date from week/tow representation 477*2d4e511cSCy Schubert */ 478*2d4e511cSCy Schubert void 479*2d4e511cSCy Schubert gpscal_to_calendar( 480*2d4e511cSCy Schubert TCivilDate * cd, 481*2d4e511cSCy Schubert TcGpsDatum * wd 482*2d4e511cSCy Schubert ) 483*2d4e511cSCy Schubert { 484*2d4e511cSCy Schubert TNtpDatum nd; 485*2d4e511cSCy Schubert 486*2d4e511cSCy Schubert memset(cd, 0, sizeof(*cd)); 487*2d4e511cSCy Schubert nd = gpsntp_from_gpscal_ex(wd, FALSE); 488*2d4e511cSCy Schubert gpsntp_to_calendar(cd, &nd); 489*2d4e511cSCy Schubert } 490*2d4e511cSCy Schubert 491*2d4e511cSCy Schubert /* ----------------------------------------------------------------- 492*2d4e511cSCy Schubert * Given the week and seconds in week, as well as the fraction/offset 493*2d4e511cSCy Schubert * (which should/could include the leap seconds offset), unfold the 494*2d4e511cSCy Schubert * weeks (which are assumed to have just 10 bits) into expanded weeks 495*2d4e511cSCy Schubert * based on the GPS base date derived from the build date (default) or 496*2d4e511cSCy Schubert * set by the configuration. 497*2d4e511cSCy Schubert * 498*2d4e511cSCy Schubert * !NOTE! This function takes RAW GPS weeks, aligned to the GPS start 499*2d4e511cSCy Schubert * (1980-01-06) on input. The output weeks will be aligned to NTPD's 500*2d4e511cSCy Schubert * week calendar start (1899-12-31)! 501*2d4e511cSCy Schubert */ 502*2d4e511cSCy Schubert TGpsDatum 503*2d4e511cSCy Schubert gpscal_from_gpsweek( 504*2d4e511cSCy Schubert uint16_t week, 505*2d4e511cSCy Schubert int32_t secs, 506*2d4e511cSCy Schubert l_fp fofs 507*2d4e511cSCy Schubert ) 508*2d4e511cSCy Schubert { 509*2d4e511cSCy Schubert TGpsDatum retv; 510*2d4e511cSCy Schubert 511*2d4e511cSCy Schubert retv.frac = 0; 512*2d4e511cSCy Schubert retv.wsecs = secs; 513*2d4e511cSCy Schubert retv.weeks = week + GPSNTP_WSHIFT; 514*2d4e511cSCy Schubert gpscal_add_offset(&retv, fofs); 515*2d4e511cSCy Schubert return _gpscal_fix_gps_era(&retv); 516*2d4e511cSCy Schubert } 517*2d4e511cSCy Schubert 518*2d4e511cSCy Schubert /* ----------------------------------------------------------------- 519*2d4e511cSCy Schubert * internal work horse for time-of-week expansion 520*2d4e511cSCy Schubert */ 521*2d4e511cSCy Schubert static TGpsDatum 522*2d4e511cSCy Schubert _gpscal_from_weektime( 523*2d4e511cSCy Schubert int32_t wsecs, 524*2d4e511cSCy Schubert l_fp fofs, 525*2d4e511cSCy Schubert TcGpsDatum * pivot 526*2d4e511cSCy Schubert ) 527*2d4e511cSCy Schubert { 528*2d4e511cSCy Schubert static const int32_t shift = SECSPERWEEK / 2; 529*2d4e511cSCy Schubert 530*2d4e511cSCy Schubert TGpsDatum retv; 531*2d4e511cSCy Schubert 532*2d4e511cSCy Schubert /* set result based on pivot -- ops order is important here */ 533*2d4e511cSCy Schubert ZERO(retv); 534*2d4e511cSCy Schubert retv.wsecs = wsecs; 535*2d4e511cSCy Schubert gpscal_add_offset(&retv, fofs); /* result is normalized */ 536*2d4e511cSCy Schubert retv.weeks = pivot->weeks; 537*2d4e511cSCy Schubert 538*2d4e511cSCy Schubert /* Manual periodic extension without division: */ 539*2d4e511cSCy Schubert if (pivot->wsecs < shift) { 540*2d4e511cSCy Schubert int32_t lim = pivot->wsecs + shift; 541*2d4e511cSCy Schubert retv.weeks -= (retv.wsecs > lim || 542*2d4e511cSCy Schubert (retv.wsecs == lim && retv.frac >= pivot->frac)); 543*2d4e511cSCy Schubert } else { 544*2d4e511cSCy Schubert int32_t lim = pivot->wsecs - shift; 545*2d4e511cSCy Schubert retv.weeks += (retv.wsecs < lim || 546*2d4e511cSCy Schubert (retv.wsecs == lim && retv.frac < pivot->frac)); 547*2d4e511cSCy Schubert } 548*2d4e511cSCy Schubert return _gpscal_fix_gps_era(&retv); 549*2d4e511cSCy Schubert } 550*2d4e511cSCy Schubert 551*2d4e511cSCy Schubert /* ----------------------------------------------------------------- 552*2d4e511cSCy Schubert * expand a time-of-week around a pivot given as week datum 553*2d4e511cSCy Schubert */ 554*2d4e511cSCy Schubert TGpsDatum 555*2d4e511cSCy Schubert gpscal_from_weektime2( 556*2d4e511cSCy Schubert int32_t wsecs, 557*2d4e511cSCy Schubert l_fp fofs, 558*2d4e511cSCy Schubert TcGpsDatum * pivot 559*2d4e511cSCy Schubert ) 560*2d4e511cSCy Schubert { 561*2d4e511cSCy Schubert TGpsDatum wpiv = * pivot; 562*2d4e511cSCy Schubert _norm_gps_datum(&wpiv); 563*2d4e511cSCy Schubert return _gpscal_from_weektime(wsecs, fofs, &wpiv); 564*2d4e511cSCy Schubert } 565*2d4e511cSCy Schubert 566*2d4e511cSCy Schubert /* ----------------------------------------------------------------- 567*2d4e511cSCy Schubert * epand a time-of-week around an pivot given as LFP, which in turn 568*2d4e511cSCy Schubert * is expanded around the current system time and then converted 569*2d4e511cSCy Schubert * into a week datum. 570*2d4e511cSCy Schubert */ 571*2d4e511cSCy Schubert TGpsDatum 572*2d4e511cSCy Schubert gpscal_from_weektime1( 573*2d4e511cSCy Schubert int32_t wsecs, 574*2d4e511cSCy Schubert l_fp fofs, 575*2d4e511cSCy Schubert l_fp pivot 576*2d4e511cSCy Schubert ) 577*2d4e511cSCy Schubert { 578*2d4e511cSCy Schubert vint64 pvi64; 579*2d4e511cSCy Schubert TGpsDatum wpiv; 580*2d4e511cSCy Schubert ntpcal_split split; 581*2d4e511cSCy Schubert 582*2d4e511cSCy Schubert /* get 64-bit pivot in NTP epoch */ 583*2d4e511cSCy Schubert pvi64 = ntpcal_ntp_to_ntp(pivot.l_ui, NULL); 584*2d4e511cSCy Schubert 585*2d4e511cSCy Schubert /* convert to weeks since 1899-12-31 and seconds in week */ 586*2d4e511cSCy Schubert pvi64 = addv64u32(&pvi64, (GPSNTP_DSHIFT * SECSPERDAY)); 587*2d4e511cSCy Schubert split = ntpcal_weeksplit(&pvi64); 588*2d4e511cSCy Schubert 589*2d4e511cSCy Schubert wpiv.weeks = split.hi; 590*2d4e511cSCy Schubert wpiv.wsecs = split.lo; 591*2d4e511cSCy Schubert wpiv.frac = pivot.l_uf; 592*2d4e511cSCy Schubert return _gpscal_from_weektime(wsecs, fofs, &wpiv); 593*2d4e511cSCy Schubert } 594*2d4e511cSCy Schubert 595*2d4e511cSCy Schubert /* ----------------------------------------------------------------- 596*2d4e511cSCy Schubert * get week/tow representation from day/tod datum 597*2d4e511cSCy Schubert */ 598*2d4e511cSCy Schubert TGpsDatum 599*2d4e511cSCy Schubert gpscal_from_gpsntp( 600*2d4e511cSCy Schubert TcNtpDatum * gd 601*2d4e511cSCy Schubert ) 602*2d4e511cSCy Schubert { 603*2d4e511cSCy Schubert TGpsDatum retv; 604*2d4e511cSCy Schubert vint64 ts64; 605*2d4e511cSCy Schubert ntpcal_split split; 606*2d4e511cSCy Schubert 607*2d4e511cSCy Schubert ts64 = ntpcal_dayjoin(gd->days, gd->secs); 608*2d4e511cSCy Schubert ts64 = addv64u32(&ts64, (GPSNTP_DSHIFT * SECSPERDAY)); 609*2d4e511cSCy Schubert split = ntpcal_weeksplit(&ts64); 610*2d4e511cSCy Schubert 611*2d4e511cSCy Schubert retv.frac = gd->frac; 612*2d4e511cSCy Schubert retv.wsecs = split.lo; 613*2d4e511cSCy Schubert retv.weeks = split.hi; 614*2d4e511cSCy Schubert return retv; 615*2d4e511cSCy Schubert } 616*2d4e511cSCy Schubert 617*2d4e511cSCy Schubert /* ----------------------------------------------------------------- 618*2d4e511cSCy Schubert * convert week/tow to LFP stamp 619*2d4e511cSCy Schubert */ 620*2d4e511cSCy Schubert l_fp 621*2d4e511cSCy Schubert ntpfp_from_gpsdatum( 622*2d4e511cSCy Schubert TcGpsDatum * gd 623*2d4e511cSCy Schubert ) 624*2d4e511cSCy Schubert { 625*2d4e511cSCy Schubert l_fp retv; 626*2d4e511cSCy Schubert 627*2d4e511cSCy Schubert retv.l_uf = gd->frac; 628*2d4e511cSCy Schubert retv.l_ui = gd->weeks * (uint32_t)SECSPERWEEK 629*2d4e511cSCy Schubert + (uint32_t)gd->wsecs 630*2d4e511cSCy Schubert - (uint32_t)SECSPERDAY * GPSNTP_DSHIFT; 631*2d4e511cSCy Schubert return retv; 632*2d4e511cSCy Schubert } 633*2d4e511cSCy Schubert 634*2d4e511cSCy Schubert /* -*-EOF-*- */ 635