I need a function (eg IsDark(time_t,lat,long))that returns boolean
whether its dark or not, when input with time(GMT), lat and long.  And
it has to use only integer mathematics.  The alogrithm has to be
implementable on a small processor, in the C language, and run
reasonably swiftly.  Absolute accuracy is not required.

---

Request for Question Clarification by answerguru-ga on 05 Jul 2002 08:10 PDT

Hi there,

Sorry to let you know that you cannot find out what you want with just
those parameters. Don't forget also that the earth revolves around the
sun, so sunset times will be different for the same place depending on
what time of year it is.

answerguru-ga

---

Request for Question Clarification by jeanluis-ga on 05 Jul 2002 08:47 PDT

Actually according to the man page for the C function time_t
time(time_t *t):
"time  returns the time since the Epoch (00:00:00 UTC, January 1,
1970), measured in seconds."

So that should be sufficient input data, however the real question is,
how much memory can we use? Because this could be done quite fast with
a table lookup function, but it would use alot of memory... Using a
table would also satisify the other requirement, that this function be
all integer.

I am also just curious however, why must it be all integer? (Just
speed issues? or does the target platform not have floating point
functions?) and what is the target platform?
Let me know
--jld

---

Clarification of Question by kaptainkazoo-ga on 05 Jul 2002 11:10 PDT

This is for a small embedded system, and I dont have any requirement
(nor an available library) for floating point other than this one odd
requirement.  Its on an 80188 processor, which has no FP support.  I'd
like the target footprint (inc tables) to be < 5Kbyte.  Ive seen
several solutions (eg http://bfn.org/~aa089/suntimes.htm,
http://www.hamiltonlabs.com/userguid/samples/sunrise.htm), but these
appear to need a lot of supporting math functions.

---

Clarification of Question by kaptainkazoo-ga on 07 Jul 2002 04:26 PDT

I'd have though ten minutes of error would be accepable.

---

Clarification of Question by kaptainkazoo-ga on 08 Jul 2002 00:19 PDT

Does it becomer easier if there is no performance requirement?  Its
just occured I can run this as a background task at low priority one a
minute, rather than as an inline function, so any need to complete
swiftly goes away.

I see, well due to the fact that I don't have time to code up a
solution to this problem I am posting this as a comment. Here are my
thoughts:
I would assemble a boolean table for both equinoxes (and possibly days
between the equinoxes, depending on how accurate the results you want,
and how much space you have available. Note the following trade off:
the more data you have in the table, the more accurate the results) At
the very least I would start with data for just the 2 equinoxes.
Anyhow in this table you would have boolean data for each hour of the
day, at each long/lat coordinate. The bool would represent weather or
not it is dark at each point in the table. For example:

+---------------+---------------+
|VernalEquinox  |AutumnalEquinox|
+---------------+---------------+
|(x,y,t):true   |(x,y,t):true   |
|(x,y,t+1):true |(x,y,t+1):false|
| ...           | ...           |
+---------------+---------------+

I dont know if that clearly illustrates it or not, but that is
supposed to represent a multi-dimensional table with data points in it
for all long/lat for every hour of the day, for 2 days. If you run out
of space then you could use fewer data points, by maybe using only
even integer long/lat points or something like that. But if you use a
cleaver indexing scheme you wouldn't have to store all the (x,y,t)
data, each data point could be 1 bit (for example, you could setup the
table to be indexed by EquinoxOffset+360*x+180*y+t, where
EquinoxOffset is the offset to the start of the Vernal or Autumnal
coloumn in the table, and assuming you have 360 longitude points, and
180 latitude points, and t is 0-23, then the input data could act as
an index into the table). Anyhow what you could then do, is when a
person inputs time/long/lat you could do a table lookup in the
VernalEquinox col, and find the closest long/lat/hour data point, and
then do the same in the AutumnalEquinox coloumn, and do a weighted
average of the two. Meaning if they are both false, then its dark, if
they are both true, then it is light, if one is dark, and one is
light, then you should take the value of the one that is closest in
time to the input time. (Meaning if the input time is 2 days after
Vernal Equinox, and it is (true)light on Vernal, but (false)Dark on
Autumnal, then you return true.)

The hardest part of this problem will be collecting all the data
points, and putting them into the table correctly... But I suppose you
could use one of the programs you mentioned to get the data, but it
doesn't sound like fun.

Anyhow I hope you understand this, and if it works for you let me
know.

By the way it just occured to me, that if the function is to be
called: IsDark() then I have my bools backwards, I was using true to
mean light, and false to mean dark. Sorry about that...

Jeanluis, don't you mean summer and winter solstice?  I don't see how
having data on the equinoxes would allow for any kind of
interpolation, which I assume is what you're aiming for.

dscotton, Yes, you are correct, thats pretty embarrassing...
I don't know what I was thinking...
Thanks for pointing that out... :)
--jld

For a clarification, what precision is required?  +/- 1, 5, 15, 60
min?

References will come later, just working from memory for now.

I agree, a table is the way to go, but you probably want a sine or tan
table, judging from the algorithm in "suntimes".  You can also use
that to calculate the arctan.

There are many symmetries for the earth and the seasons.  The nothern
hemisphere is the opposite of the southern.  The time from winter to
summer is the mirror of summer to winter.  The time from winter to
spring is the "inverse" of spring to summer

There are some variances due to the earth's orbit not being circular. 
I think it's at most +/- 15 minutes.  That can also be corrected.

Longitude is not really a factor.  It is a simple calculation for an
offset for the time.  (I recommend the book "Longitude").  24hr = 360
degrees.  You can remove that as an index for the table.

An alternative table would be to store sunrise/sunset times and
interpolate that.  There are locations and times that the sun doesn't
set or rise.

Don't forget about leap years.

No new code here, but I'll try to summerize how to modify this code
from Geo. Anderson.  The following code except for lines beginning
with ; is from the following website.
http://www3.edgenet.net/lingling/x10-sun.html
; figure out which day of the year
N = Int(275 * MonthFor / 9) - 2 * Int((MonthFor + 9) / 12) + DayFor -
30
; use that to figure out the angle we are from the sun
; if we ignore longitude, we are at most one day off
; use a table based on N (or higher resolution to include longitude)
;   to precompute most of the following
LO = 4.8771 + 0.0172 * (N + 0.5 - Longitude / 360)
C = 0.03342 * Sin(LO + 1.345)
; ??? something looks odd taking the atn(tan(x)) = x
C2 = (1 / RadiansPerDegree) * (Atn(Tan(LO + C)) - Atn(0.9175 * Tan(LO
+C)) - C)
SD = 0.3978 * Sin(LO + C)
CD = Sqr(1 - SD * SD)
; SD (and CD) are precomputed for N
; Sin(lat) is also precomputed. Linear interpolate between values in
table.
; Cos(lat) = Sin(lat+90deg)
; you might consider SC = tan(D)*tan(lat) + (sun below
horizon/(cos(lat)*cos(D))
; as latitude increases toward the poles, the tan(lat) and the 
;    second term becomes more significant (no sun or midnight sun)
; the sun appears early because of atmospheric bending and the sun is
a disk
'SC = (SD * Sin(Latitude * RadiansPerDegree) + 0.0145) / (Cos(Latitude
* RadiansPerDegree) * CD) 'this one calculates true astronomical
sunrise and sunset
; you might want this to indicate whether or not it is dark
SC = (SD * Sin(Latitude * RadiansPerDegree) + 0.104539) /
(Cos(Latitude * RadiansPerDegree) * CD) 'this one calculates "civil
twilight"
; since you have a sine table, compute the arcsin(SC)
C3 = (1 / RadiansPerDegree) * Atn(SC / Sqr(1 - SC * SC))
; I'm not sure what happens for the midnight sun here.
; TimeZone will be zero since you are using GMT
; longitude takes care of the offset
Dawn = DateFor + (6 - TimeZone - (Longitude + C2 + C3) / 15) / 24
Dusk = DateFor + (18 - TimeZone - (Longitude + C2 - C3) / 15) / 24

You can apply some the techniques suggested above to the other
algorithms you have.

Website to test your algorithm.
http://aa.usno.navy.mil/data/docs/RS_OneYear.html

I seem to remember some suggestions in the "Graphics Gems" books for
trig and interger math.
http://www.acm.org/tog/GraphicsGems/

Fixed Point.
http://dec.bournemouth.ac.uk/staff/awatson/micro/articles/9508_034.pdf

Some discussion on trig functions.
http://www.karlscalculus.org/trig_id.html
http://www.ganssle.com/articles/atrig.htm
http://www.restena.lu/convict/Jeunes/Math/arctan.htm

More of the same algorithms.
http://www.prairienet.org/cuas/faqcalc.shtml
http://www.btinternet.com/~kburnett/kepler/
http://skyandtelescope.com/resources/software/programs/sunup.bas