1 /* WMGlobe 0.5 - All the Earth on a WMaker Icon
2 * copyright (C) 1998,99 Jerome Dumonteil <jerome.dumonteil@capway.com>
3 * sunpos.c is taken from Xearth :
10 * code for calculating the position on the earth's surface for which
11 * the sun is directly overhead (adapted from _practical astronomy
12 * with your calculator, third edition_, peter duffett-smith,
13 * cambridge university press, 1988.)
16 * Copyright (C) 1989, 1990, 1993, 1994, 1995 Kirk Lauritz Johnson
18 * Parts of the source code (as marked) are:
19 * Copyright (C) 1989, 1990, 1991 by Jim Frost
20 * Copyright (C) 1992 by Jamie Zawinski <jwz@lucid.com>
22 * Permission to use, copy, modify and freely distribute xearth for
23 * non-commercial and not-for-profit purposes is hereby granted
24 * without fee, provided that both the above copyright notice and this
25 * permission notice appear in all copies and in supporting
28 * Unisys Corporation holds worldwide patent rights on the Lempel Zev
29 * Welch (LZW) compression technique employed in the CompuServe GIF
30 * image file format as well as in other formats. Unisys has made it
31 * clear, however, that it does not require licensing or fees to be
32 * paid for freely distributed, non-commercial applications (such as
33 * xearth) that employ LZW/GIF technology. Those wishing further
34 * information about licensing the LZW patent should contact Unisys
35 * directly at (lzw_info@unisys.com) or by writing to
38 * Welch Licensing Department
43 * The author makes no representations about the suitability of this
44 * software for any purpose. It is provided "as is" without express or
47 * THE AUTHOR DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
48 * INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS,
49 * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY SPECIAL, INDIRECT
50 * OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM
51 * LOSS OF USE, DATA OR PROFITS, WHETHER IN AN ACTION OF CONTRACT,
52 * NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF OR IN
53 * CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
56 /*************************************************************************/
62 #define PI 3.141592653
67 * the epoch upon which these astronomical calculations are based is
68 * 1990 january 0.0, 631065600 seconds since the beginning of the
69 * "unix epoch" (00:00:00 GMT, Jan. 1, 1970)
71 * given a number of seconds since the start of the unix epoch,
72 * DaysSinceEpoch() computes the number of days since the start of the
73 * astronomical epoch (1990 january 0.0)
76 #define EpochStart (631065600)
77 #define DaysSinceEpoch(secs) (((secs)-EpochStart)*(1.0/(24*3600)))
80 * assuming the apparent orbit of the sun about the earth is circular,
81 * the rate at which the orbit progresses is given by RadsPerDay --
82 * TWOPI radians per orbit divided by 365.242191 days per year:
85 #define RadsPerDay (TWOPI/365.242191)
88 * details of sun's apparent orbit at epoch 1990.0 (after
89 * duffett-smith, table 6, section 46)
91 * Epsilon_g (ecliptic longitude at epoch 1990.0) 279.403303 degrees
92 * OmegaBar_g (ecliptic longitude of perigee) 282.768422 degrees
93 * Eccentricity (eccentricity of orbit) 0.016713
96 #define Epsilon_g (279.403303*(TWOPI/360))
97 #define OmegaBar_g (282.768422*(TWOPI/360))
98 #define Eccentricity (0.016713)
101 * MeanObliquity gives the mean obliquity of the earth's axis at epoch
102 * 1990.0 (computed as 23.440592 degrees according to the method given
103 * in duffett-smith, section 27)
105 #define MeanObliquity (23.440592*(TWOPI/360))
108 * solve Kepler's equation via Newton's method
109 * (after duffett-smith, section 47)
111 static double solve_keplers_equation(double M
)
118 delta
= E
- Eccentricity
* sin(E
) - M
;
119 if (fabs(delta
) <= 1e-10)
121 E
-= delta
/ (1 - Eccentricity
* cos(E
));
129 * compute ecliptic longitude of sun (in radians)
130 * (after duffett-smith, section 47)
132 static double sun_ecliptic_longitude(time_t ssue
)
133 /* seconds since unix epoch */
139 D
= DaysSinceEpoch(ssue
);
146 M_sun
= N
+ Epsilon_g
- OmegaBar_g
;
150 E
= solve_keplers_equation(M_sun
);
151 v
= 2 * atan(sqrt((1 + Eccentricity
) / (1 - Eccentricity
)) * tan(E
/ 2));
153 return (v
+ OmegaBar_g
);
158 * convert from ecliptic to equatorial coordinates
159 * (after duffett-smith, section 27)
161 static void ecliptic_to_equatorial(double lambda
, double beta
, double *alpha
, double *delta
)
163 * double lambda; ecliptic longitude
164 * double beta; ecliptic latitude
165 * double *alpha; (return) right ascension
166 * double *delta; (return) declination
171 sin_e
= sin(MeanObliquity
);
172 cos_e
= cos(MeanObliquity
);
174 *alpha
= atan2(sin(lambda
) * cos_e
- tan(beta
) * sin_e
, cos(lambda
));
175 *delta
= asin(sin(beta
) * cos_e
+ cos(beta
) * sin_e
* sin(lambda
));
180 * computing julian dates (assuming gregorian calendar, thus this is
181 * only valid for dates of 1582 oct 15 or later)
182 * (after duffett-smith, section 4)
184 static double julian_date(int y
, int m
, int d
)
186 * int y; year (e.g. 19xx)
187 * int m; month (jan=1, feb=2, ...)
188 * int d; day of month
194 /* lazy test to ensure gregorian calendar */
199 if ((m
== 1) || (m
== 2)) {
205 C
= (int) (365.25 * y
);
206 D
= (int) (30.6001 * (m
+ 1));
208 JD
= B
+ C
+ D
+ d
+ 1720994.5;
215 * compute greenwich mean sidereal time (GST) corresponding to a given
216 * number of seconds since the unix epoch
217 * (after duffett-smith, section 12)
219 static double GST(time_t ssue
)
220 /*time_t ssue; seconds since unix epoch */
229 JD
= julian_date(tm
->tm_year
+ 1900, tm
->tm_mon
+ 1, tm
->tm_mday
);
230 T
= (JD
- 2451545) / 36525;
232 T0
= ((T
+ 2.5862e-5) * T
+ 2400.051336) * T
+ 6.697374558;
238 UT
= tm
->tm_hour
+ (tm
->tm_min
+ tm
->tm_sec
/ 60.0) / 60.0;
240 T0
+= UT
* 1.002737909;
250 * given a particular time (expressed in seconds since the unix
251 * epoch), compute position on the earth (lat, lon) such that sun is
254 void GetSunPos(time_t ssue
, double *lat
, double *lon
)
255 /* time_t ssue; seconds since unix epoch */
256 /* double *lat; (return) latitude */
257 /* double *lon; (return) longitude */
263 lambda
= sun_ecliptic_longitude(ssue
);
264 ecliptic_to_equatorial(lambda
, 0.0, &alpha
, &delta
);
266 tmp
= alpha
- (TWOPI
/ 24) * GST(ssue
);
271 } else if (tmp
> PI
) {