[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[avrlibccommit] [2369] Augmented the comments on several functions.
From: 
Mike Rice 
Subject: 
[avrlibccommit] [2369] Augmented the comments on several functions. 
Date: 
Sun, 28 Apr 2013 14:19:39 +0000 
Revision: 2369
http://svn.sv.gnu.org/viewvc/?view=rev&root=avrlibc&revision=2369
Author: swfltek
Date: 20130428 14:19:35 +0000 (Sun, 28 Apr 2013)
Log Message:

Augmented the comments on several functions. Hopefully the algorithms employed
will be clear.
Modified Paths:

trunk/avrlibc/libc/time/daylight_seconds.c
trunk/avrlibc/libc/time/equation_of_time.c
trunk/avrlibc/libc/time/gm_sidereal.c
trunk/avrlibc/libc/time/gmtime_r.c
trunk/avrlibc/libc/time/iso_weeknum.c
trunk/avrlibc/libc/time/mk_gmtime.c
trunk/avrlibc/libc/time/solar_declination.c
Modified: trunk/avrlibc/libc/time/daylight_seconds.c
===================================================================
 trunk/avrlibc/libc/time/daylight_seconds.c 20130427 15:49:22 UTC (rev
2368)
+++ trunk/avrlibc/libc/time/daylight_seconds.c 20130428 14:19:35 UTC (rev
2369)
@@ 29,8 +29,9 @@
/* $Id$ */
/*
 Determine the amount of time the sun is above the horizon. At high
latitudes, this can be zero,
 or >= ONE_DAY.
+ Determine the amount of time the sun is above the horizon. At high
latitudes, around the
+ solstices, this can be zero or greater than ONE_DAY.
+
*/
#include <time.h>
@@ 49,17 +50,20 @@
d = solar_declination(timer);
+ /* partial 'Sunrise Equation' */
d = tan(l) * tan(d);
 /* d may exceed a magnitude of 1.0 at high latitudes */
+ /* magnitude of d may exceed 1.0 at near solstices */
if (d > 1.0)
d = 1.0;
if (d < 1.0)
d = 1.0;
+ /* derive hour angle */
d = acos(d);
+ /* but for atmospheric refraction, this would be d /= M_PI */
d /= 3.112505;
n = ONE_DAY * d;
Modified: trunk/avrlibc/libc/time/equation_of_time.c
===================================================================
 trunk/avrlibc/libc/time/equation_of_time.c 20130427 15:49:22 UTC (rev
2368)
+++ trunk/avrlibc/libc/time/equation_of_time.c 20130428 14:19:35 UTC (rev
2369)
@@ 31,6 +31,19 @@
/*
The so called Equation of Time.
+ The eccentricity of Earths orbit contributes about 7.7 minutes of
variation to the result. It
+ has a period of 1 anomalous year, with zeroes at perihelion and aphelion.
+
+ The tilt of Earths rotational axis (obliquity) contributes about 9.9
minutes of variation. It
+ has a period of 1/2 tropical year, with zeroes at solstices and equinoxes.
The time of Earths
+ arrival at these events is influenced by the eccentricity, which causes it
to progress along its
+ orbital path faster as it approaches perihelion, imposing a 'modulation'
on the tropical phase.
+
+ The algorithm employed computes the orbital position with respect to
perihelion, deriving
+ from that a 'velocity correction factor'. The orbital position with
respect to the winter solstice
+ is then computed, as modulated by that factor. The individual
contributions of the obliquity and the
+ eccentricity components are then summed, and returned as an integer value
in seconds.
+
*/
#include <time.h>
@@ 58,14 +71,15 @@
s += SOLSTICE;
s *= 2;
sf = s;
+ sf /= TROP_CYCLE;
 sf /= TROP_CYCLE;
+ /* modulate to derive actual position */
sf += dV;
+ sf = sin(sf);
 sf = sin(sf);
+ /* compute contributions */
+ sf *= 592.2;
pf *= 459.6;
 sf *= 592.2;

s = pf + sf;
return s;
Modified: trunk/avrlibc/libc/time/gm_sidereal.c
===================================================================
 trunk/avrlibc/libc/time/gm_sidereal.c 20130427 15:49:22 UTC (rev
2368)
+++ trunk/avrlibc/libc/time/gm_sidereal.c 20130428 14:19:35 UTC (rev
2369)
@@ 29,9 +29,17 @@
/* $Id$ */
/*
 Greenwich Mean Sidereal Time. The ratio of sidereal to civil seconds is
 1.00273790935ss / 1.0cs , which when shifted left to fill 32 bits becomes
 0x8059B740.
+ Greenwich Mean Sidereal Time. A sidereal second is somewhat shorter than a
standard second,
+ about 1.002737909350795 sidereal seconds per standard second.
+
+ We resort to fixed point math due to the insufficient resolution of a
'double', using...
+
+ timestamp * ( 1.002737909350795 << 31 )
+  + Te
+ 1 << 31
+
+ Where Te is the sidereal time at the epoch.
+
*/
#include <time.h>
@@ 40,20 +48,13 @@
unsigned long
gm_sidereal(const time_t * timer)
{
 uint64_t tmp;
+ uint64_t tmp;
tmp = *timer;

 /* multiply by 1.00273790935 sidereal seconds */
tmp *= 0x8059B740;
+ tmp /= 0x80000000;
+ tmp += (uint64_t) 23991;
 /* divide by 1.0 civil second */
 tmp >>= 31;

 /* add sidereal time at epoch */
 tmp += (uint64_t)23991;

tmp %= ONE_DAY;
return tmp;
}

Modified: trunk/avrlibc/libc/time/gmtime_r.c
===================================================================
 trunk/avrlibc/libc/time/gmtime_r.c 20130427 15:49:22 UTC (rev 2368)
+++ trunk/avrlibc/libc/time/gmtime_r.c 20130428 14:19:35 UTC (rev 2369)
@@ 28,9 +28,7 @@
/* $Id$ */
/*
 Re entrant version of gmtime(), this function 'breaks down' a Y2K time
stamp .
*/
+/* Re entrant version of gmtime(). */
#include <time.h>
#include <stdlib.h>
@@ 39,88 +37,108 @@
void
gmtime_r(const time_t * timer, struct tm * timeptr)
{
 int32_t fract;
 ldiv_t lresult;
 div_t result;
 uint16_t days, n, leapyear, years;
+ int32_t fract;
+ ldiv_t lresult;
+ div_t result;
+ uint16_t days, n, leapyear, years;
 /* break down timer into whole and fractional parts of 1 day */
 days = *timer / 86400UL;
 fract = *timer % 86400UL;
+ /* break down timer into whole and fractional parts of 1 day */
+ days = *timer / 86400UL;
+ fract = *timer % 86400UL;
 /* Determine day of week ( the epoch was a Saturday ) */
 n = days + SATURDAY;
 n %= 7;
 timeptr>tm_wday = n;
+ /*
+ Extract hour, minute, and second from the fractional day
+ */
+ lresult = ldiv(fract, 60L);
+ timeptr>tm_sec = lresult.rem;
+ result = div(lresult.quot, 60);
+ timeptr>tm_min = result.rem;
+ timeptr>tm_hour = result.quot;
 /* map our place into the 100 year and 4 year leap cycles. */
 lresult = ldiv((long) days, 36525L);
 years = 100 * lresult.quot;
 lresult = ldiv(lresult.rem, 1461L);
 years += 4 * lresult.quot;
 days = lresult.rem;
 if (years > 100)
 days++;
+ /* Determine day of week ( the epoch was a Saturday ) */
+ n = days + SATURDAY;
+ n %= 7;
+ timeptr>tm_wday = n;
 /*
 * 'years' is now at a 4 year boundary
 * 'days' is an index into the current 1461 day leap cycle
 * If 'years' != 100, this cycle begins with a 366 day year.
+ /*
+ * Our epoch year has the property of being at the conjunction of all
three 'leap cycles',
+ * 4, 100, and 400 years ( though we can ignore the 400 year cycle in
this library).
+ *
+ * Using this property, we can easily 'map' the time stamp into the
leap cycles, quickly
+ * deriving the year and day of year, along with the fact of whether it
is a leap year.
+ */
+
+ /* map into a 100 year cycle */
+ lresult = ldiv((long) days, 36525L);
+ years = 100 * lresult.quot;
+
+ /* map into a 4 year cycle */
+ lresult = ldiv(lresult.rem, 1461L);
+ years += 4 * lresult.quot;
+ days = lresult.rem;
+ if (years > 100)
+ days++;
+
+ /*
+ * 'years' is now at the first year of a 4 year leap cycle, which will
always be a leap year,
+ * unless it is 100. 'days' is now an index into that cycle.
*/
 leapyear = 1;
 if (years == 100)
 leapyear = 0;
+ leapyear = 1;
+ if (years == 100)
+ leapyear = 0;
 n = 364 + leapyear;
+ /* compute length, in days, of first year of this cycle */
+ n = 364 + leapyear;
 /*
 Given the length of the first year of this cycle,
 we can proceed to break down year and day of year.
 */
 if (days > n) {
 days = leapyear;
 leapyear = 0;
 result = div(days, 365);
 years += result.quot;
 days = result.rem;
 }
 timeptr>tm_year = 100 + years;
 timeptr>tm_yday = days;
+ /*
+ * if the number of days remaining is greater than the length of the
+ * first year, we make one more division.
+ */
+ if (days > n) {
+ days = leapyear;
+ leapyear = 0;
+ result = div(days, 365);
+ years += result.quot;
+ days = result.rem;
+ }
+ timeptr>tm_year = 100 + years;
+ timeptr>tm_yday = days;
 /*
+ /*
Given the year, day of year, and leap year indicator, we can break
down the
 month and day of month.
+ month and day of month. If the day of year is less than 59 (or 60
if a leap year), then
+ we handle the Jan/Feb month pair as an exception.
*/
 n = 59 + leapyear;
+ n = 59 + leapyear;
+ if (days < n) {
+ /* special case: Jan/Feb month pair */
+ result = div(days, 31);
+ timeptr>tm_mon = result.quot;
+ timeptr>tm_mday = result.rem;
+ } else {
+ /*
+ The remaining 10 months form a regular pattern of 31 day months
alternating with 30 day
+ months, with a 'phase change' between July and August (153 days
after March 1).
+ We proceed by mapping our position into either MarchJuly or
AugustDecember.
+ */
+ days = n;
+ result = div(days, 153);
+ timeptr>tm_mon = 2 + result.quot * 5;
 if (days < n) {
 result = div(days, 31);
 timeptr>tm_mon = result.quot;
 timeptr>tm_mday = result.rem;
 } else {
 days = n;
 result = div(days, 153);
 timeptr>tm_mon = 2 + result.quot * 5;
 result = div(result.rem, 61);
 timeptr>tm_mon += result.quot * 2;
 result = div(result.rem, 31);
 timeptr>tm_mon += result.quot;
 timeptr>tm_mday = result.rem;
 }
+ /* map into a 61 day pair of months */
+ result = div(result.rem, 61);
+ timeptr>tm_mon += result.quot * 2;
 /*
 Break down hour, minute, and second from the fractional day
 */
 lresult = ldiv(fract, 60L);
 timeptr>tm_sec = lresult.rem;
 result = div(lresult.quot, 60);
 timeptr>tm_min = result.rem;
 timeptr>tm_hour = result.quot;
+ /* map into a month */
+ result = div(result.rem, 31);
+ timeptr>tm_mon += result.quot;
+ timeptr>tm_mday = result.rem;
+ }
 /*
+ /*
Cleanup and return
*/
 timeptr>tm_isdst = 0; /* gmt is never in DST */
 timeptr>tm_mday++; /* tm_mday is 1 based */
+ timeptr>tm_isdst = 0; /* gmt is never in DST */
+ timeptr>tm_mday++; /* tm_mday is 1 based */
}
Modified: trunk/avrlibc/libc/time/iso_weeknum.c
===================================================================
 trunk/avrlibc/libc/time/iso_weeknum.c 20130427 15:49:22 UTC (rev
2368)
+++ trunk/avrlibc/libc/time/iso_weeknum.c 20130428 14:19:35 UTC (rev
2369)
@@ 28,27 +28,40 @@
/* $Id$ */
/* Compute the ISO 8601 week number of the year. */
+/*
+ Compute the ISO 8601 week number of the year.
+ We return 0 if the week is the final one of the previous year.
+ We return 54 if the week is the first week of the following year.
+ Otherwise we return week numbers 1 to 53
+*/
#include <time.h>
uint8_t
iso_weeknum(const struct tm * timestruct)
{
 int d, w;
+ int d, w;
 d = timestruct>tm_wday;
 if (d == 0)
 d = 7;
 w = timestruct>tm_yday + 11  d;
 w /= 7;
 if (w == 53) {
 /* week 53 must include its thursday in the same year */
 d = timestruct>tm_mday  1;
 d = timestruct>tm_wday;
 d += THURSDAY;
 if (d > 30)
 w++; /* signal first week of the following year */
 }
 return w;
+ /* convert to a MONDAY based week */
+ d = timestruct>tm_wday;
+ if (d == 0)
+ d = 7;
+
+ /* compute tentative ISO 8601 week number */
+ w = timestruct>tm_yday + 11  d;
+ w /= 7;
+
+ /*
+ * handle the special case where week 53 may actually be week 1 of
+ * the following year
+ */
+ if (w == 53) {
+ /* week 53 must include its thursday in the same year */
+ d = timestruct>tm_mday  1;
+ d = timestruct>tm_wday;
+ d += THURSDAY;
+ if (d > 30)
+ w++; /* signal first week of the following year */
+ }
+ return w;
}
Modified: trunk/avrlibc/libc/time/mk_gmtime.c
===================================================================
 trunk/avrlibc/libc/time/mk_gmtime.c 20130427 15:49:22 UTC (rev
2368)
+++ trunk/avrlibc/libc/time/mk_gmtime.c 20130428 14:19:35 UTC (rev
2369)
@@ 31,6 +31,7 @@
/*
'Break down' a y2k time stamp into the elements of struct tm.
Unlike mktime(), this function does not 'normalize' the elements of
timeptr.
+
*/
#include <time.h>
@@ 44,7 +45,8 @@
int n, m, d, leaps;
/*
 Determine elapsed whole days, since Y2K, to the beginning of the
year
+ Determine elapsed whole days since the epoch to the beginning of this
year. Since our epoch is
+ at a conjunction of the leap cycles, we can do this rather quickly.
*/
n = timeptr>tm_year  100;
leaps = 0;
@@ 57,10 +59,13 @@
tmp = 365UL * n + leaps;
/*
 Derive the day of year from month and day of month
 */
+ Derive the day of year from month and day of month. We use the
pattern of 31 day months
+ followed by 30 day months to our advantage, but we must
'special case' Jan/Feb, and
+ account for a 'phase change' between July and August (153 days
after March 1).
+ */
d = timeptr>tm_mday  1; /* tm_mday is one based */
+ /* handle Jan/Feb as a special case */
if (timeptr>tm_mon < 2) {
if (timeptr>tm_mon)
d += 31;
@@ 68,19 +73,26 @@
} else {
n = 59 + is_leap_year(timeptr>tm_year + 1900);
d += n;
+ n = timeptr>tm_mon  MARCH;
 /* Other months exhibit a regular pattern */
 n = timeptr>tm_mon  2;
 if (n > 4)
+ /* account for phase change */
+ if (n > (JULY  MARCH))
d += 153;
+ n %= 5;
 n %= 5;
+ /*
+ * n is now an index into a group of alternating 31 and 30
+ * day months... 61 day pairs.
+ */
m = n / 2;
m *= 61;
d += m;
 n &= 1;
 if (n)
+ /*
+ * if n is odd, we are in the second half of the
+ * month pair
+ */
+ if (n & 1)
d += 31;
}
@@ 94,6 +106,7 @@
tmp *= ONE_HOUR;
tmp += timeptr>tm_min * 60UL;
tmp += timeptr>tm_sec;
+
ret += tmp;
return ret;
Modified: trunk/avrlibc/libc/time/solar_declination.c
===================================================================
 trunk/avrlibc/libc/time/solar_declination.c 20130427 15:49:22 UTC
(rev 2368)
+++ trunk/avrlibc/libc/time/solar_declination.c 20130428 14:19:35 UTC
(rev 2369)
@@ 28,6 +28,18 @@
/* $Id$ */
+/*
+ Were it not for the eccentricity of Earths orbit, this would be a trivial
function.
+
+ We compute the Earths orbital position with respect to perihelion, from
which we derive a
+ 'velocity correction factor'. We then compute the orbital angle with
respect to the
+ December solstice, as 'modulated' by that correction factor.
+
+ Due to the accumulation of rounding errors, the computed December solstice
of 2135 will lag
+ the actual solstice by many hours. A fudge factor, 'LAG', distributes the
error across
+ the 136 year range of this library.
+*/
+
#include <time.h>
#include <math.h>
#include "ephemera_common.h"
@@ 41,28 +53,25 @@
uint32_t fT, oV;
double dV, dT;
 /*
 * Determine approximate orbital angle, relative to the winter
 * solstice
 */
 fT = *timer % TROP_YEAR;
 fT += SOLSTICE;
 fT += LAG;
 dT = fT;
 dT /= TROP_CYCLE;

 /* Determine approximate orbital angle, relative to perihelion */
+ /* Determine orbital angle relative to perihelion of January 1999 */
oV = *timer % ANOM_YEAR;
oV += PERIHELION;
dV = oV;
dV /= ANOM_CYCLE;
 /* Derive a velocity correction factor from the perihelion angle */
+ /* Derive velocity correction factor from the perihelion angle */
dV = sin(dV);
dV *= DELTA_V;
 /* compute declination */
 dT = cos(dT + dV) * INCLINATION;
+ /* Determine orbital angle relative to solstice of December 1999 */
+ fT = *timer % TROP_YEAR;
+ fT += SOLSTICE + LAG;
+ dT = fT;
+ dT /= TROP_CYCLE;
+ dT += dV;
+ /* Finally having the solstice angle, we can compute the declination */
+ dT = cos(dT) * INCLINATION;
+
return dT;
}
[Prev in Thread] 
Current Thread 
[Next in Thread] 
 [avrlibccommit] [2369] Augmented the comments on several functions.,
Mike Rice <=