author | kergoth <kergoth> | 2002-01-25 22:14:26 (UTC) |
---|---|---|
committer | kergoth <kergoth> | 2002-01-25 22:14:26 (UTC) |
commit | 15318cad33835e4e2dc620d033e43cd930676cdd (patch) (unidiff) | |
tree | c2fa0399a2c47fda8e2cd0092c73a809d17f68eb /core/settings/citytime/sun.c | |
download | opie-15318cad33835e4e2dc620d033e43cd930676cdd.zip opie-15318cad33835e4e2dc620d033e43cd930676cdd.tar.gz opie-15318cad33835e4e2dc620d033e43cd930676cdd.tar.bz2 |
Initial revision
Diffstat (limited to 'core/settings/citytime/sun.c') (more/less context) (ignore whitespace changes)
-rw-r--r-- | core/settings/citytime/sun.c | 323 |
1 files changed, 323 insertions, 0 deletions
diff --git a/core/settings/citytime/sun.c b/core/settings/citytime/sun.c new file mode 100644 index 0000000..d3f3731 --- a/dev/null +++ b/core/settings/citytime/sun.c | |||
@@ -0,0 +1,323 @@ | |||
1 | /* | ||
2 | * Sun clock. X11 version by John Mackin. | ||
3 | * | ||
4 | * This program was derived from, and is still in part identical with, the | ||
5 | * Suntools Sun clock program whose author's comment appears immediately | ||
6 | * below. Please preserve both notices. | ||
7 | * | ||
8 | * The X11R3/4 version of this program was written by John Mackin, at the | ||
9 | * Basser Department of Computer Science, University of Sydney, Sydney, | ||
10 | * New South Wales, Australia; <john@cs.su.oz.AU>. This program, like | ||
11 | * the one it was derived from, is in the public domain: `Love is the | ||
12 | * law, love under will.' | ||
13 | */ | ||
14 | |||
15 | /* | ||
16 | |||
17 | Sun clock | ||
18 | |||
19 | Designed and implemented by John Walker in November of 1988. | ||
20 | |||
21 | Version for the Sun Workstation. | ||
22 | |||
23 | The algorithm used to calculate the position of the Sun is given in | ||
24 | Chapter 18 of: | ||
25 | |||
26 | "Astronomical Formulae for Calculators" by Jean Meeus, Third Edition, | ||
27 | Richmond: Willmann-Bell, 1985. This book can be obtained from: | ||
28 | |||
29 | Willmann-Bell | ||
30 | P.O. Box 35025 | ||
31 | Richmond, VA 23235 | ||
32 | USA | ||
33 | Phone: (804) 320-7016 | ||
34 | |||
35 | This program was written by: | ||
36 | |||
37 | John Walker | ||
38 | Autodesk, Inc. | ||
39 | 2320 Marinship Way | ||
40 | Sausalito, CA 94965 | ||
41 | USA | ||
42 | Fax: (415) 389-9418 | ||
43 | Voice: (415) 332-2344 Ext. 2829 | ||
44 | Usenet: {sun,well,uunet}!acad!kelvin | ||
45 | or: kelvin@acad.uu.net | ||
46 | |||
47 | modified for interactive maps by | ||
48 | |||
49 | Stephen Martin | ||
50 | Fujitsu Systems Business of Canada | ||
51 | smartin@fujitsu.ca | ||
52 | |||
53 | This program is in the public domain: "Do what thou wilt shall be the | ||
54 | whole of the law". I'd appreciate receiving any bug fixes and/or | ||
55 | enhancements, which I'll incorporate in future versions of the | ||
56 | program. Please leave the original attribution information intactso | ||
57 | that credit and blame may be properly apportioned. | ||
58 | |||
59 | Revision history: | ||
60 | |||
61 | 1.0 12/21/89 Initial version. | ||
62 | 8/24/89 Finally got around to submitting. | ||
63 | |||
64 | 1.1 8/31/94 Version with interactive map. | ||
65 | 1.2 10/12/94 Fixes for HP and Solaris, new icon bitmap | ||
66 | 1.3 11/01/94 Timezone now shown in icon | ||
67 | 1.4 03/29/98 Fixed city drawing, added icon animation | ||
68 | |||
69 | */ | ||
70 | |||
71 | #include "sun.h" | ||
72 | |||
73 | #include <qpe/qmath.h> | ||
74 | |||
75 | /* PROJILLUM -- Project illuminated area on the map. */ | ||
76 | |||
77 | void | ||
78 | projillum(wtab, xdots, ydots, dec) | ||
79 | short *wtab; | ||
80 | int xdots, ydots; | ||
81 | double dec; | ||
82 | { | ||
83 | int i, ftf = 1, ilon, ilat, lilon = 0, lilat = 0, xt; | ||
84 | double m, x, y, z, th, lon, lat, s, c; | ||
85 | |||
86 | /* Clear unoccupied cells in width table */ | ||
87 | |||
88 | for (i = 0; i < ydots; i++) | ||
89 | wtab[i] = -1; | ||
90 | |||
91 | /* Build transformation for declination */ | ||
92 | |||
93 | s = qSin(-dtr(dec)); | ||
94 | c = qCos(-dtr(dec)); | ||
95 | |||
96 | /* Increment over a semicircle of illumination */ | ||
97 | |||
98 | for (th = -(PI / 2); th <= PI / 2 + 0.001; | ||
99 | th += PI / TERMINC) { | ||
100 | |||
101 | /* Transform the point through the declination rotation. */ | ||
102 | |||
103 | x = -s * qSin(th); | ||
104 | y = qCos(th); | ||
105 | z = c * qSin(th); | ||
106 | |||
107 | /* Transform the resulting co-ordinate through the | ||
108 | map projection to obtain screen co-ordinates. */ | ||
109 | |||
110 | lon = (y == 0 && x == 0) ? 0.0 : rtd(qATan2(y, x)); | ||
111 | lat = rtd(qASin(z)); | ||
112 | |||
113 | ilat = ydots - (lat + 90) * (ydots / 180.0); | ||
114 | ilon = lon * (xdots / 360.0); | ||
115 | |||
116 | if (ftf) { | ||
117 | |||
118 | /* First time. Just save start co-ordinate. */ | ||
119 | |||
120 | lilon = ilon; | ||
121 | lilat = ilat; | ||
122 | ftf = 0; | ||
123 | } else { | ||
124 | |||
125 | /* Trace out the line and set the width table. */ | ||
126 | |||
127 | if (lilat == ilat) { | ||
128 | wtab[(ydots - 1) - ilat] = ilon == 0 ? 1 : ilon; | ||
129 | } else { | ||
130 | m = ((double) (ilon - lilon)) / (ilat - lilat); | ||
131 | for (i = lilat; i != ilat; i += sgn(ilat - lilat)) { | ||
132 | xt = lilon + qFloor((m * (i - lilat)) + 0.5); | ||
133 | wtab[(ydots - 1) - i] = xt == 0 ? 1 : xt; | ||
134 | } | ||
135 | } | ||
136 | lilon = ilon; | ||
137 | lilat = ilat; | ||
138 | } | ||
139 | } | ||
140 | |||
141 | /* Now tweak the widths to generate full illumination for | ||
142 | the correct pole. */ | ||
143 | |||
144 | if (dec < 0.0) { | ||
145 | ilat = ydots - 1; | ||
146 | lilat = -1; | ||
147 | } else { | ||
148 | ilat = 0; | ||
149 | lilat = 1; | ||
150 | } | ||
151 | |||
152 | for (i = ilat; i != ydots / 2; i += lilat) { | ||
153 | if (wtab[i] != -1) { | ||
154 | while (1) { | ||
155 | wtab[i] = xdots / 2; | ||
156 | if (i == ilat) | ||
157 | break; | ||
158 | i -= lilat; | ||
159 | } | ||
160 | break; | ||
161 | } | ||
162 | } | ||
163 | } | ||
164 | |||
165 | /* | ||
166 | * Sun clock - astronomical routines. | ||
167 | */ | ||
168 | |||
169 | /* JDATE -- Convert internal GMT date and time to Julian day | ||
170 | and fraction. */ | ||
171 | |||
172 | long | ||
173 | jdate(t) | ||
174 | struct tm *t; | ||
175 | { | ||
176 | long c, m, y; | ||
177 | |||
178 | y = t->tm_year + 1900; | ||
179 | m = t->tm_mon + 1; | ||
180 | if (m > 2) | ||
181 | m = m - 3; | ||
182 | else { | ||
183 | m = m + 9; | ||
184 | y--; | ||
185 | } | ||
186 | c = y / 100L; /* Compute century */ | ||
187 | y -= 100L * c; | ||
188 | return t->tm_mday + (c * 146097L) / 4 + (y * 1461L) / 4 + | ||
189 | (m * 153L + 2) / 5 + 1721119L; | ||
190 | } | ||
191 | |||
192 | /* JTIME -- Convert internal GMT date andtime to astronomical | ||
193 | Julian time (i.e. Julian date plus day fraction, | ||
194 | expressed as a double).*/ | ||
195 | |||
196 | double | ||
197 | jtime(t) | ||
198 | struct tm *t; | ||
199 | { | ||
200 | return (jdate(t) - 0.5) + | ||
201 | (((long) t->tm_sec) + | ||
202 | 60L * (t->tm_min + 60L * t->tm_hour)) / 86400.0; | ||
203 | } | ||
204 | |||
205 | /* KEPLER --Solve the equation of Kepler. */ | ||
206 | |||
207 | double | ||
208 | kepler(m, ecc) | ||
209 | double m, ecc; | ||
210 | { | ||
211 | double e, delta; | ||
212 | #define EPSILON 1E-6 | ||
213 | |||
214 | e = m = dtr(m); | ||
215 | do { | ||
216 | delta = e - ecc * qSin(e) - m; | ||
217 | e -= delta / (1 - ecc * qCos(e)); | ||
218 | } while (qFabs(delta) > EPSILON); | ||
219 | return e; | ||
220 | } | ||
221 | |||
222 | /* SUNPOS -- Calculate position of the Sun.JD is the Julian date | ||
223 | of the instant for which the position is desired and | ||
224 | APPARENT should be nonzero if the apparent position | ||
225 | (corrected for nutation and aberration) is desired. | ||
226 | The Sun's co-ordinates are returned in RA and DEC, | ||
227 | both specified in degrees (divide RA by 15 to obtain | ||
228 | hours). The radius vector to the Sun in astronomical | ||
229 | units is returned in RV and the Sun's longitude (true | ||
230 | or apparent, as desired) is returned as degrees in | ||
231 | SLONG.*/ | ||
232 | |||
233 | |||
234 | void | ||
235 | sunpos(jd, apparent, ra, dec, rv, slong) | ||
236 | double jd; | ||
237 | int apparent; | ||
238 | double *ra, *dec, *rv, *slong; | ||
239 | { | ||
240 | double t, t2, t3, l, m, e, ea, v, theta, omega, | ||
241 | eps; | ||
242 | |||
243 | /* Time, in Julian centuries of 36525 ephemeris days, | ||
244 | measured from the epoch 1900 January 0.5 ET. */ | ||
245 | |||
246 | t = (jd - 2415020.0) / 36525.0; | ||
247 | t2 = t * t; | ||
248 | t3 = t2 * t; | ||
249 | |||
250 | /* Geometric mean longitude of the Sun, referred to the | ||
251 | mean equinox of the date. */ | ||
252 | |||
253 | l = fixangle(279.69668 + 36000.76892 * t + 0.0003025 * t2); | ||
254 | |||
255 | /* Sun's mean anomaly. */ | ||
256 | |||
257 | m = fixangle(358.47583 + 35999.04975*t - 0.000150*t2 - 0.0000033*t3); | ||
258 | |||
259 | /* Eccentricity of the Earth's orbit. */ | ||
260 | |||
261 | e = 0.01675104 - 0.0000418 * t - 0.000000126 * t2; | ||
262 | |||
263 | /* Eccentric anomaly. */ | ||
264 | |||
265 | ea = kepler(m, e); | ||
266 | |||
267 | /* True anomaly */ | ||
268 | |||
269 | v = fixangle(2 * rtd(qATan(qSqrt((1 + e) / (1 - e)) * qTan(ea / 2)))); | ||
270 | |||
271 | /* Sun's true longitude. */ | ||
272 | |||
273 | theta = l + v - m; | ||
274 | |||
275 | /* Obliquity of the ecliptic. */ | ||
276 | |||
277 | eps = 23.452294 - 0.0130125 * t - 0.00000164 * t2 + 0.000000503 * t3; | ||
278 | |||
279 | /* Corrections for Sun's apparent longitude, if desired. */ | ||
280 | |||
281 | if (apparent) { | ||
282 | omega = fixangle(259.18 - 1934.142 * t); | ||
283 | theta = theta - 0.00569 - 0.00479 * qSin(dtr(omega)); | ||
284 | eps += 0.00256 * qCos(dtr(omega)); | ||
285 | } | ||
286 | |||
287 | /* Return Sun's longitude and radius vector */ | ||
288 | |||
289 | *slong = theta; | ||
290 | *rv = (1.0000002 * (1 - e * e)) / (1 + e * qCos(dtr(v))); | ||
291 | |||
292 | /* Determine solar co-ordinates. */ | ||
293 | |||
294 | *ra = | ||
295 | fixangle(rtd(qATan2(qCos(dtr(eps)) * qSin(dtr(theta)), qCos(dtr(theta))))); | ||
296 | *dec = rtd(qASin(sin(dtr(eps)) * qSin(dtr(theta)))); | ||
297 | } | ||
298 | |||
299 | /* GMST -- Calculate Greenwich Mean Siderial Time for a given | ||
300 | instant expressed as a Julian date and fraction.*/ | ||
301 | |||
302 | double | ||
303 | gmst(jd) | ||
304 | double jd; | ||
305 | { | ||
306 | double t, theta0; | ||
307 | |||
308 | |||
309 | /* Time, in Julian centuries of 36525 ephemeris days, | ||
310 | measured from the epoch 1900 January 0.5 ET. */ | ||
311 | |||
312 | t = ((qFloor(jd + 0.5) - 0.5) - 2415020.0) / 36525.0; | ||
313 | |||
314 | theta0 = 6.6460656 + 2400.051262 * t + 0.00002581 * t * t; | ||
315 | |||
316 | t = (jd + 0.5) - (qFloor(jd + 0.5)); | ||
317 | |||
318 | theta0 += (t * 24.0) * 1.002737908; | ||
319 | |||
320 | theta0 = (theta0 - 24.0 * (qFloor(theta0 / 24.0))); | ||
321 | |||
322 | return theta0; | ||
323 | } | ||