Sunrise/Sunset Algorithm
Neal Bridges
http://quartus.net/files/PalmOS/Forth/Examples/sun.zip. It's bundled
with documentation and the original algorithm.
It's written for Quartus Forth, but may well be useful to others. If
you have any comments, please let me know!
Thanks,
--
Neal Bridges
Quartus Handheld Software
http://www.quartus.net
Home of the Quartus Forth on-board compiler for the Palm OS!
Ron Aaron
> I've written a routine to calculate sunrise/sunset -- it's at
> http://quartus.net/files/PalmOS/Forth/Examples/sun.zip. It's bundled
> with documentation and the original algorithm.
>
> It's written for Quartus Forth, but may well be useful to others. If
> you have any comments, please let me know!
Beatiful, Neil - I was contemplating converting some of my own C code for this
very purpose.
But -- it uses "float" words. Do you have plans for an integer-only version?
Best regards,
Ron
--
Reva Forth - http://ronware.org/reva/
Marcel Hendrix
[..]
Nice. Here is a translation for Win32.
-marcel
-- -----------
(*
* LANGUAGE : ANS Forth with extensions
* PROJECT : Forth Environments
* DESCRIPTION : A Forth library for Calculating Sunrise and Sunset
* CATEGORY : Utility
* AUTHOR : Neal Bridges, August 2005
* LAST CHANGE : August 19, 2005, Marcel Hendrix
*)
NEEDS -miscutil
NEEDS -dynlink
NEEDS -dates
REVISION -sunrise "ÄÄÄ Sunrise/sunset Version 1.00 ÄÄÄ"
PRIVATES
DOC
(*
The sun library module provides words for calculating sunrise/sunset times
for a specified latitude/longitude. The times are available both in UTC and
in the local timezone selected in Preferences on each PalmOS (v4 and
greater) device.
The source algorithm for the time calculation is from:
Almanac for Computers, 1990
Published by the Nautical Almanac Office
United States Naval Observatory
Washington, DC 20392.
Online:
http://williams.best.vwh.net/sunrise_sunset_algorithm.htm
Local copy: sunrise_sunset_algorithm.htm
This algorithm is valid for the years 1980-2050, and is generally accurate
to within a minute. Accuracy degrades when the location is further north or
south than 60 degrees of latitude; above 80 degrees of latitude, small
inaccuracies in the calculation of the sun's ecliptic will result in greater
errors in the time of sunrise and sunset.
The sun library module uses Quartus Forth internal single-precision floats.
MathLib is not required for this module.
Provides:
set-location ( F: longitude latitude -- )
Sets the location used by the sunrise/sunset words below. West longitudes
and south latitudes are negative; east longitudes and north latitudes are
positive.
set-zenith ( F: zenith -- )
Sets the zenith used by the sunrise/sunset words below. Applicable named
zenith words are provided, as follows:
astronomical-zenith ( -- )
civil-zenith ( -- )
nautical-zenith ( -- )
official-zenith ( -- ) The default when sun is loaded.
Each of these sets a value for zenith.
rising ( -- rising-flag )
setting ( -- setting-flag )
These flags control whether UTC-suntime calculates sunrise, or sunset.
UTC-suntime ( day month year set? -- ) ( F: -- h.m )
Calculates the time of the rising or the setting of the sun for the current
location on the specified date, using the current zenith. The value
returned is a floating-point hour.minute in UTC (GMT timezone). set? must
be either rising or setting.
If the sun does not rise or set at the selected lat/long on the given day,
-11 (result out of range) is thrown during the time calculation. Use CATCH
to trap that as required.
sunrise ( day month year -- ) ( F: -- h.m )
sunset ( day month year -- ) ( F: -- h.m )
These return sunrise and sunset, respectively, for the current location on
the specified date, using the current zenith, corrected for local time using
the Palm preferences settings for timezone and daylight-savings. The value
returned is a floating-point hour.minute.
time>mh ( F: h.m -- ) ( -- m h )
Converts a floating point hour.minute value into minute and hour integers on
the stack.
day-of-year ( day month year -- day-of-year )
Computes the day number for the specified date. January 1=1, January 2=2,
etc. The year is only used to determine leap-year status.
Example of use:
-79.4e 43.6e set-location
( west longitudes and south latitudes are negative )
official-zenith
21 August 2005 sunrise time>mh . .
1 July 1990 sunset time>mh . .
typedef struct _TIME_ZONE_INFORMATION { // tzi
LONG Bias;
WCHAR StandardName[ 32 ];
SYSTEMTIME StandardDate;
LONG StandardBias;
WCHAR DaylightName[ 32 ];
SYSTEMTIME DaylightDate;
LONG DaylightBias;
} TIME_ZONE_INFORMATION;
kernel32.dll:
DWORD GetTimeZoneInformation( LPTIME_ZONE_INFORMATION lpTimeZoneInformation );
returns TIME_ZONE_ID_DAYLIGHT ( 2 ) if operating under daylight saving conditions
All translations between UTC time and local time are based on the following formula:
UTC = local time + bias
*)
ENDDOC
-- ------------------------------------------------------------------------------------------------------------------------
2 =: TIME_ZONE_ID_DAYLIGHT PRIVATE
2 CELLS #32 2* 2* CHARS + 4 CELLS 2* + =: /daylightbias PRIVATE
0 =: /bias PRIVATE
/daylightbias CELL+ =: /TZIDINFO PRIVATE
CREATE timezoneinfo PRIVATE /TZIDINFO ALLOT
: TIMEZONEINFO@ ( -- )
0 0 LOCALS| 'tz@ dll |
S" kernel32.dll" LIBRARY-OPEN ABORT" library-open failed" TO dll
S" GetTimeZoneInformation" dll LIBRARY-FIND ABORT" library-search failed" TO 'tz@
LIBRARY-ERROR ?DUP IF TYPE EXIT ENDIF
timezoneinfo 1 'tz@ FOREIGN 0= ABORT" GetTimeZoneInformation :: failed" ( -- returnflag)
TIME_ZONE_ID_DAYLIGHT <> IF timezoneinfo /daylightbias + OFF ENDIF
dll LIBRARY-CLOSE ABORT" library-close failed" ;
-- Local latitude and longitude
0e FVALUE latitude -- west and south are negative, east and north are positive
0e FVALUE longitude
0e FVALUE zenith -- Sun's zenith for sunrise/sunset
: set-location ( F: long lat -- ) TO latitude TO longitude ;
: set-zenith ( F: zenith -- ) TO zenith ;
-- Builds zenith-setting words.
: ZENITH: ( F: f -- )
CREATE HERE 1 FLOATS ALLOT F!
DOES> ( -- ) F@ set-zenith ; PRIVATE
90.83333e ZENITH: official-zenith
96e ZENITH: civil-zenith
102e ZENITH: nautical-zenith
108e ZENITH: astronomical-zenith
-- Calculate the day-of-year number of a given date (January 1 == month 1, day 1).
-- Uses the DATES library.
: day-of-year ( d m y -- day )
DUP >R JDAY
1 1 ( January ) R> JDAY - 1+ ;
-- ------------------------------------------------------------------------------------------------------------------------
-- Adjust so the range is [0,360).
: range360 ( F: f1 -- f2 ) 360e FPREM FDUP F0< IF 360e F+ ENDIF ; PRIVATE
-- Adjust so the range is [0,23).
: range24 ( F: f1 -- f2 ) 24e FPREM FDUP F0< IF 24e F+ ENDIF ; PRIVATE
-- Round down to the nearest multiple of 90.
: floor90 ( F: f1 -- f2 ) 90e FP/REM FDROP S>F 90e F* ; PRIVATE
-- Convert a floating-point h.m time into integer minutes and hours.
: time>mh ( F: h.m -- ) ( -- min hour )
FDUP FLOOR FOVER FSWAP F-
60e F* F>S FLOOR F>S ;
FALSE =: rising
TRUE =: setting
: SIN ( F: deg -- r ) FRAD FSIN ; PRIVATE
: COS ( F: deg -- r ) FRAD FCOS ; PRIVATE
: TAN ( F: deg -- r ) FRAD FTAN ; PRIVATE
: ATAN ( F: r -- deg ) FATAN FDEG ; PRIVATE
: ASIN ( F: r -- deg ) FASIN FDEG ; PRIVATE
: ACOS ( F: r -- deg ) FACOS FDEG ; PRIVATE
-- Calculate the UTC sunrise or sunset time for a given day of the year,
-- using the location set in the longitude and latitude fvalues.
: UTC-suntime ( d m y set? -- ) ( F: -- h.m )
LOCAL set?
0e 0e 0e 0e 0e 0e 0e 0e 0e FLOCALS| H L RA T M lngHour sinDec cosDec cosH |
\ convert the longitude to hour value and calculate an approximate time
\ lngHour = longitude / 15
\ if rising time is desired: t = N + ((6 - lngHour) / 24)
\ if setting time is desired: t = N + ((18 - lngHour) / 24)
day-of-year S>F TO T
longitude 15e F/ TO lngHour
set? rising = IF 6e ELSE 18e ENDIF lngHour F- 24e F/ +TO T
\ calculate the Sun's mean anomaly M = (0.9856 * t) - 3.289
0.9856e T F* 3.289e F- TO M
\ calculate the Sun's true longitude
\ L = M + 1.916 * sin(M) + 0.020 * sin(2 * M) + 282.634
\ NOTE: L potentially needs to be adjusted into the range [0,360) by adding/subtracting 360
M M SIN 1.916e F* F+ M F2* SIN 0.02e F* F+ 282.634e F+ range360 TO L
\ calculate the Sun's right ascension
\ RA = atan(0.91764 * tan(L))
\ NOTE: RA potentially needs to be adjusted into the range [0,360) by adding/subtracting 360
L TAN 0.91764e F* ATAN range360 TO RA
\ right ascension value needs to be in the same quadrant as L
\ Lquadrant = (floor( L/90)) * 90
\ RAquadrant = (floor(RA/90)) * 90
\ RA = RA + (Lquadrant - RAquadrant)
L floor90 RA floor90 F- RA F+
\ right ascension value needs to be converted into hours
\ RA = RA / 15
15e F/ TO RA
\ calculate the Sun's declination
\ sinDec = 0.39782 * sin(L)
\ cosDec = cos(asin(sinDec))
L SIN 0.39782e F* FDUP TO sinDec
ASIN COS TO cosDec
\ calculate the Sun's local hour angle
\ cosH = (cos(zenith) - (sinDec * sin(latitude))) / (cosDec * cos(latitude))
\ if (cosH > 1) the sun never rises on this location (on the specified date)
\ if (cosH < -1) the sun never sets on this location (on the specified date)
zenith COS latitude SIN sinDec F* F-
latitude COS cosDec F* F/ TO cosH
cosH 1e F> IF -11 THROW ENDIF \ "the sun never rises on this location"
cosH -1e F< IF -12 THROW ENDIF \ "the sun never sets on this location"
\ finish calculating H and convert into hours
\ if rising time is desired: H = 360 - acos(cosH)
\ if setting time is desired: H = acos(cosH)
\ H = H / 15
cosH ACOS set? rising = IF 360e FSWAP F- ENDIF
15e F/ TO H
\ calculate local mean time of rising/setting
\ T = H + RA - (0.06571 * t) - 6.622
H RA F+ T 0.06571e F* F- 6.622e F- TO T
\ adjust back to UTC
\ UT = T - lngHour
\ NOTE: UT potentially needs to be adjusted into the range [0,24) by adding/subtracting 24
T lngHour F- range24 ;
-- Toronto, Canada: 43.6N 79.4W
-79.4e 43.6e set-location
official-zenith
20 7 ( July) 1989 setting UTC-suntime
CR .( Canada 20 July 1989, expect sunset at 0:53 -> ) time>mh 0 .R ':' EMIT .
20 7 ( July) 1989 rising UTC-suntime
CR .( Canada 20 July 1989, expect sunrise at 9:54 -> ) time>mh 0 .R ':' EMIT .
\ Netherlands 52'30'' North 5'45'' East
5.45e 52.3e set-location
official-zenith
19 8 ( Aug) 2005 setting UTC-suntime
CR .( Netherlands 19 Aug 2005, expect sunset at 18:55 -> ) time>mh 0 .R ':' EMIT .
19 8 ( Aug) 2005 rising UTC-suntime
CR .( Netherlands 19 Aug 2005, expect sunrise at 4:27 -> ) time>mh 0 .R ':' EMIT .
-- Return the total offset in minutes of the timezone and DST settings.
: local-offset ( -- local-offset )
timezoneinfo@ timezoneinfo /bias + @ timezoneinfo /daylightbias + @ + ;
-- Calculate sunrise or sunset time for the specified date, adjusting for the local timezone & DST.
: local-suntime ( d m y set? -- ) ( F: -- h.m ) UTC-suntime local-offset S>F 60e F/ F- range24 ;
: sunrise ( d m y -- ) ( F: -- h.m ) rising local-suntime ;
: sunset ( d m y -- ) ( F: -- h.m ) setting local-suntime ;
:ABOUT CR ." Try: ( d m y -- ) ( F: -- h.m ) sunrise -- sunrise in minutes at date"
CR ." ( d m y -- ) ( F: -- h.m ) sunset -- sunset in minutes at date" ;
.ABOUT -sunrise CR
DEPRIVE
(* End of Source *)
Anton Ertl
>"Neal Bridges" <quar...@gmail.com> writes:
>
>[..]
>
>Nice. Here is a translation for Win32.
>-- -----------
>(*
...
>DOC
...
>ENDDOC
>
>-- ------------------------------------------------------------------------------------------------------------------------
...
>FALSE =: rising
What do these iForthisms have to do with Win32?
>\ convert the longitude to hour value and calculate an approximate time
Did the program smell iForthish enough already, so you could stand the
presence of standard comments?
>-- Return the total offset in minutes of the timezone and DST settings.
Apparently not.
- anton
--
M. Anton Ertl http://www.complang.tuwien.ac.at/anton/home.html
comp.lang.forth FAQs: http://www.complang.tuwien.ac.at/forth/faq/toc.html
New standard: http://www.complang.tuwien.ac.at/forth/ansforth/forth200x.html
EuroForth 2005: http://www.complang.tuwien.ac.at/anton/euroforth2005/
Neal Bridges
Ron, I haven't planned to convert it. I did a quick look-around for an
integer formula initially, but found none. I'm not sure how this can
be done without the trig functions, and I've never seen an integer set
with sufficient resolution.
Ron Aaron
> Ron, I haven't planned to convert it. I did a quick look-around for an
> integer formula initially, but found none. I'm not sure how this can
> be done without the trig functions, and I've never seen an integer set
> with sufficient resolution.
My problem exactly :-)
It seems to me that since the algorithm is in any case only accurate to within
2 minutes, it should be possible to work out an integer (32bit) version.
Anyway, you've given me fresh impetus to work on the floating-point library in
Reva. :-)
All the best,
Paul E. Bennett
> Ron Aaron wrote:
>> On 16 Aug 2005 19:51:45 -0700, Neal Bridges <qua...@gmail.com> wrote:
>> > I've written a routine to calculate sunrise/sunset -- it's at
>> > http://quartus.net/files/PalmOS/Forth/Examples/sun.zip. It's bundled
>> > with documentation and the original algorithm.
>> >
>> > It's written for Quartus Forth, but may well be useful to others. If
>> > you have any comments, please let me know!
>>
>> Beatiful, Neil - I was contemplating converting some of my own C code for
>> this very purpose.
>>
>> But -- it uses "float" words. Do you have plans for an integer-only
>> version? Best regards,
>> Ron
>
> Ron, I haven't planned to convert it. I did a quick look-around for an
> integer formula initially, but found none. I'm not sure how this can
> be done without the trig functions, and I've never seen an integer set
> with sufficient resolution.
Is there a link to the description of the algorithm.
I was also wondering how fine grained you needed the resolution for the
calculations. 16 bits is enough for about 1.4 seconds of granularity in
timing (24 hours being a full cycle of the 16 bit cell). If you go to 32
bits the granularity gets you into microseconds. I'll agree that the maths
does require a bit more thought but I have already published a routine that
converts the resultant numbers to the required display format.
--
********************************************************************
Paul E. Bennett ....................<email://p...@amleth.demon.co.uk>
Forth based HIDECS Consultancy .....<http://www.amleth.demon.co.uk/>
Mob: +44 (0)7811-639972
Tel: +44 (0)1235-811095
Going Forth Safely ....EBA. http://www.electric-boat-association.org.uk/
********************************************************************
Neal Bridges
Neal Bridges
[snip]
> Is there a link to the description of the algorithm.
>
> I was also wondering how fine grained you needed the resolution for the
> calculations. 16 bits is enough for about 1.4 seconds of granularity in
> timing (24 hours being a full cycle of the 16 bit cell). If you go to 32
> bits the granularity gets you into microseconds. I'll agree that the maths
> does require a bit more thought but I have already published a routine that
> converts the resultant numbers to the required display format.
Paul, the algorithm is described in the .htm contained in the .zip
file. You can also find it here:
http://williams.best.vwh.net/sunrise_sunset_algorithm.htm
Krishna Myneni
A version for kForth (sun.4th) is available at
ftp://ccreweb.org/software/kforth/examples/
Look for both sun.4th and jd.4th. All fp calculations are in standard
Forth. See comments at top of sun.4th.
Krishna Myneni
--
-79.4e 43.6e set-location official-zenith
ok
20 July 1989 setting UTC-suntime time>mh . .
0 53 ok
20 July 1989 rising UTC-suntime time>mh . .
9 54 ok
Doug Hoffman
calculation-intensive such that it would be a good candidate for
something like Julian Noble's FORmulaTRANslator utility. Perhaps
slightly modified to use FVALUEs instead of FVARIABLEs. The following
comments could then become the source code:
M = (0.9856 * t) - 3.289
L = M + 1.916 * sin(M) + 0.020 * sin(2 * M) + 282.634
RA = RA / 15
sinDec = 0.39782 * sin(L)
cosDec = cos(asin(sinDec))
cosH = (cos(zenith) - (sinDec * sin(latitude))) / (cosDec *
cos(latitude))
T = H + RA - (0.06571 * t) - 6.622
Regards,
-Doug
Neal Bridges
In fact, this is exactly what I did in the original source. Have a
look at http://quartus.net/files/PalmOS/Forth/Examples/sun.zip. I used
Wil Baden's OPG to do the formula translation.
Paul E. Bennett
> Paul E. Bennett wrote:
> [snip]
>> Is there a link to the description of the algorithm.
>>
>> I was also wondering how fine grained you needed the resolution for the
>> calculations. 16 bits is enough for about 1.4 seconds of granularity in
>> timing (24 hours being a full cycle of the 16 bit cell). If you go to 32
>> bits the granularity gets you into microseconds. I'll agree that the
>> maths does require a bit more thought but I have already published a
>> routine that converts the resultant numbers to the required display
>> format.
>
> Paul, the algorithm is described in the .htm contained in the .zip
> file. You can also find it here:
> http://williams.best.vwh.net/sunrise_sunset_algorithm.htm
Thanks Neal,
I shall give that a thorough read a bit later. On a quick glance there
didn't seem to be too much to it. I guess I shall be dragging out my copy
of Hastings again to look up appropriate polynomial calculations to deal
with the trig functions.
Ron Aaron
> There is, of course, http://tinyboot.com/float.txt -- a starting point.
I'm using the native x86 FPU; so the only thing really holding me back is the
input and output words, e.g. ">float" and "f."
Alex McDonald
> On 19 Aug 2005 17:55:44 -0700, Neal Bridges <qua...@gmail.com> wrote:
>
>>There is, of course, http://tinyboot.com/float.txt -- a starting point.
>
>
> I'm using the native x86 FPU; so the only thing really holding me back is the
> input and output words, e.g. ">float" and "f."
>
Then pinch them from Win32Forth! And once you've rewritten them, send
them back...
--
Regards
Alex McDonald
Ron Aaron
<ar.mc...@virgin.nospam.net> wrote:
> Then pinch them from Win32Forth! And once you've rewritten them, send
> them back...
Good idea, but you probably won't want them back after I've mangled them :-)
Neal Bridges
> On Mon, 22 Aug 2005 00:36:34 GMT, Alex McDonald
> <ar.mc...@virgin.nospam.net> wrote:
>
> > Then pinch them from Win32Forth! And once you've rewritten them, send
> > them back...
>
> Good idea, but you probably won't want them back after I've mangled them :-)
Some already-invented wheels:
ftp.taygeta.com/pub/Forth/Applications/ANS/fpout18.f
http://tinyboot.com/float.txt
You'll find there sources for F. and >FLOAT, both.
Ron Aaron
> Some already-invented wheels:
>
> ftp.taygeta.com/pub/Forth/Applications/ANS/fpout18.f
> http://tinyboot.com/float.txt
>
> You'll find there sources for F. and >FLOAT, both.
>
Excellent leads, thanks!
Anton Ertl
>Neal Bridges wrote:
>> I've written a routine to calculate sunrise/sunset -- it's at
>> http://quartus.net/files/PalmOS/Forth/Examples/sun.zip. It's bundled
>> with documentation and the original algorithm.
>A version for kForth (sun.4th) is available at
>
>ftp://ccreweb.org/software/kforth/examples/
>
>Look for both sun.4th and jd.4th. All fp calculations are in standard
>Forth. See comments at top of sun.4th.
Thanks. I combined the two files and eliminated non-ANS Forth words
(DEG>RAD RAD>DEF F>); the result runs on Gforth, iForth and probably
other ANS Forth systems.
BTW, SPACE is a core word, so I am surprised that it's not built into
kForth. Similarly, D>S is a double word.
I have put the result on
<http://www.complang.tuwien.ac.at/forth/programs/sun.4th>
The program uses the following words
from CORE :
: here swap allot ; [ ] Constant spaces 2drop drop >r - dup 0< IF + r>
1- THEN * / r@ */ > negate 1+ cr mod DO i over = ." ELSE 0= LOOP
Create DOES> s>d
from CORE-EXT :
pick 2r> .r false true
from BLOCK-EXT :
\
from EXCEPTION-EXT :
abort"
from FACILITY-EXT :
time&date
from FILE :
(
from FLOAT :
fconstant fswap f< f/ FLiteral f* fvariable f! f@ fover f>d fdup f0<
f+ f- floor d>f
from FLOAT-EXT :
dfloats fsin fcos ftan fasin facos fatan fabs
Krishna Myneni
> ...
> BTW, SPACE is a core word, so I am surprised that it's not built into
> kForth. Similarly, D>S is a double word.
Originally, kForth's built-in vocabulary was restricted to 256 words,
since it compiles to byte code. Therefore, I had to be selective in
choosing the built-in words. Words such as SPACE and D>S were simple
to define in source (e.g. D>S is simply DROP for 2's complement system).
Recent version(s) of kForth removed the restriction of 256 intrinsic
words, although it still compiles to byte code! I should probably go
back and add some words I originally omitted.
Krishna
The Beez'
e.g. "SPACE" was expanded to "BL EMIT", e.g.
/*
This function compiles a NIP
*/
#ifndef ARCHAIC
static void DoNip (void)
#else
static void DoNip ()
#endif
{
InlineWords (2, SWAP, DROP);
}
/*
This function compiles a TUCK
*/
#ifndef ARCHAIC
static void DoTuck (void)
#else
static void DoTuck ()
#endif
{
InlineWords (2, SWAP, OVER);
}
/*
This function compiles a -ROT
*/
#ifndef ARCHAIC
static void DoMinRot (void)
#else
static void DoMinRot ()
#endif
{
InlineWords (2, ROT, ROT);
}
So I still have 256 tokens, but 150 of them are still free..!
Hans Bezemer