[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[avr-libc-commit] [2369] Augmented the comments on several functions.
|
From: |
Mike Rice |
|
Subject: |
[avr-libc-commit] [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=avr-libc&revision=2369
Author: swfltek
Date: 2013-04-28 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/avr-libc/libc/time/daylight_seconds.c
trunk/avr-libc/libc/time/equation_of_time.c
trunk/avr-libc/libc/time/gm_sidereal.c
trunk/avr-libc/libc/time/gmtime_r.c
trunk/avr-libc/libc/time/iso_weeknum.c
trunk/avr-libc/libc/time/mk_gmtime.c
trunk/avr-libc/libc/time/solar_declination.c
Modified: trunk/avr-libc/libc/time/daylight_seconds.c
===================================================================
--- trunk/avr-libc/libc/time/daylight_seconds.c 2013-04-27 15:49:22 UTC (rev
2368)
+++ trunk/avr-libc/libc/time/daylight_seconds.c 2013-04-28 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/avr-libc/libc/time/equation_of_time.c
===================================================================
--- trunk/avr-libc/libc/time/equation_of_time.c 2013-04-27 15:49:22 UTC (rev
2368)
+++ trunk/avr-libc/libc/time/equation_of_time.c 2013-04-28 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/avr-libc/libc/time/gm_sidereal.c
===================================================================
--- trunk/avr-libc/libc/time/gm_sidereal.c 2013-04-27 15:49:22 UTC (rev
2368)
+++ trunk/avr-libc/libc/time/gm_sidereal.c 2013-04-28 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/avr-libc/libc/time/gmtime_r.c
===================================================================
--- trunk/avr-libc/libc/time/gmtime_r.c 2013-04-27 15:49:22 UTC (rev 2368)
+++ trunk/avr-libc/libc/time/gmtime_r.c 2013-04-28 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 March-July or
August-December.
+ */
+ 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/avr-libc/libc/time/iso_weeknum.c
===================================================================
--- trunk/avr-libc/libc/time/iso_weeknum.c 2013-04-27 15:49:22 UTC (rev
2368)
+++ trunk/avr-libc/libc/time/iso_weeknum.c 2013-04-28 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/avr-libc/libc/time/mk_gmtime.c
===================================================================
--- trunk/avr-libc/libc/time/mk_gmtime.c 2013-04-27 15:49:22 UTC (rev
2368)
+++ trunk/avr-libc/libc/time/mk_gmtime.c 2013-04-28 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/avr-libc/libc/time/solar_declination.c
===================================================================
--- trunk/avr-libc/libc/time/solar_declination.c 2013-04-27 15:49:22 UTC
(rev 2368)
+++ trunk/avr-libc/libc/time/solar_declination.c 2013-04-28 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] |
- [avr-libc-commit] [2369] Augmented the comments on several functions.,
Mike Rice <=