1 /*!
2 * \file rc-astro.c
3 * \brief Astronomical functions necessary for managing the fixed dates.
4 */
5 /*
6 * Copyright (c) 1994, 95, 96, 1997, 2000 Thomas Esken
7 * Copyright (c) 2010, 2011, 2013 Free Software Foundation, Inc.
8 *
9 * This software doesn't claim completeness, correctness or usability.
10 * On principle I will not be liable for ANY damages or losses (implicit
11 * or explicit), which result from using or handling my software.
12 * If you use this software, you agree without any exception to this
13 * agreement, which binds you LEGALLY !!
14 *
15 * This program is free software; you can redistribute it and/or modify
16 * it under the terms of the `GNU General Public License' as published by
17 * the `Free Software Foundation'; either version 3, or (at your option)
18 * any later version.
19 *
20 * You should have received a copy of the `GNU General Public License'
21 * along with this program; if not, write to the:
22 *
23 */
24
25
26
27 /*
28 * Include definition header file to see whether USE_RC is defined there.
29 * Compile this module only if USE_RC is defined, otherwise skip it.
30 */
31 #include "tailor.h"
32
33
34
35 #if USE_RC
36
37
38 /*
39 * Include header files.
40 */
41 # if HAVE_MATH_H && HAVE_LIBM
42 # include <math.h>
43 # endif
44 # include "common.h"
45 # include "rc-defs.h"
46 # include "globals.h"
47 # include "hd-astro.h"
48 # include "rc-astro.h"
49 # include "utils.h"
50 # include "rc-utils.h"
51
52
53
54 /*
55 * static functions prototypes.
56 */
57 __BEGIN_DECLARATIONS
58 /*
59 ************************************************** Defined in `rc-astro.c'.
60 */
61 static double
62 atmospheric_refraction __P_ ((const double altitude,
63 double pressure, const double temperature));
64 static int moon_charpos __P_ ((const double x, const int lines));
65 static double
66 internal_moon_rise_set __P_ ((const Aevent_enum event,
67 int day,
68 int month,
69 int year, Coor_struct * coordinates));
70 __END_DECLARATIONS
71 /*
72 * Function implementations.
73 */
74 static double
atmospheric_refraction(altitude,pressure,temperature)75 atmospheric_refraction (altitude, pressure, temperature)
76 const double altitude;
77 double pressure;
78 const double temperature;
79 /*!
80 Returns the approximate atmospheric refraction for the given true radians
81 ALTITUDE in radians, the atmospheric PRESSURE in Newton per sqare meter
82 (Nm^-2), and the atmospheric TEMPERATURE in degrees Celsius, which results
83 the apparent altitude when added to the true altitude.
84 Atmospheric refraction is calculated for altitude values that are
85 greater equal -2.0 degrees upto +90.0 degrees. If PRESSURE is set
86 to a value less equal zero, no atmospheric refraction will be calculated,
87 i.e. this function always returns 0.0 (radians).
88 This function is taken from the `aa.arc' archive, which is a nice and
89 precise ephemeris calculation software written by Stephen L. Moshier
90 <moshier@mediaone.net>; adapted and pretty-printed according to Gcal's
91 needs.
92 References:
93 * U.S. Government Printing Office, "The Astronomical Almanac" (AA), 1986.
94 */
95 {
96 auto double degrees_altitude = TODEG (altitude);
97
98
99 /*
100 Convert PRESSURE given in Newton per square meter (Nm^-2) to Millibar.
101 */
102 pressure /= 100.0;
103 if (pressure <= 0.0
104 || degrees_altitude < -2.0 || degrees_altitude >= DEGS_PER_06_HOURS)
105 /*
106 No atmospheric refraction for such altitudes.
107 */
108 return (0.0);
109 if (degrees_altitude < DEGS_PER_HOUR)
110 {
111 auto double old_altitude = degrees_altitude;
112 auto double old_degrees = 0.0;
113 auto double new_altitude = degrees_altitude;
114 auto double new_degrees = 0.0;
115 auto double new_pressure = (pressure - 80.0) / 930.0;
116 auto double new_temperature = 0.0048 * (temperature - 10.0);
117 auto double n;
118 register int i;
119
120
121 /*
122 Manage altitudes below 15 degrees (== 1.0 hour) above horizon.
123 The formula for low altitudes is from the "Almanac for Computers".
124 It gives the correction for the observed altitude, so has to be
125 inverted numerically to get the observed from the true.
126 Accuracy about 0.2' for -20 C < TEMPERATURE < +40 C
127 and 970 mb < PRESSURE < 1050 mb.
128 */
129 for (i = 0; i < 4; i++)
130 {
131 n =
132 1.0 / tan (TORAD (new_altitude + (7.31 / (new_altitude + 4.4))));
133 new_degrees =
134 n * new_pressure / (60.0 + new_temperature * (n + 39.0));
135 n = new_altitude - old_altitude;
136 /*
137 Denominator of derivative.
138 */
139 old_altitude = new_degrees - old_degrees - n;
140 if ((n != 0.0) && (old_altitude != 0.0))
141 /*
142 Newton's iteration with numerically estimated derivative.
143 */
144 n =
145 new_altitude - n * (degrees_altitude + new_degrees -
146 new_altitude) / old_altitude;
147 else
148 /*
149 Can't do it on first pass.
150 */
151 n = degrees_altitude + new_degrees;
152 old_altitude = new_altitude;
153 old_degrees = new_degrees;
154 new_altitude = n;
155 }
156
157 return (TORAD (new_degrees));
158 }
159
160 /*
161 Manage altitudes above 15 degrees above horizon.
162 */
163 return (0.00007888888 * pressure /
164 ((273.0 + temperature) * tan (altitude)));
165 }
166
167
168
169 static int
moon_charpos(x,lines)170 moon_charpos (x, lines)
171 const double x;
172 const int lines;
173 /*!
174 Computes the position where to place the next character
175 of the Moon text graphics image.
176 */
177 {
178 register int i;
179
180
181 i = (int) (x * (double) lines - 0.5);
182 if (i + lines < 1)
183 i = 0;
184 else
185 i += lines;
186
187 return (i);
188 }
189
190
191
192 double
gd_latitude2gc_latitude(gd_latitude,meters_above_sea_level,gc_latitude)193 gd_latitude2gc_latitude (gd_latitude, meters_above_sea_level, gc_latitude)
194 const double gd_latitude;
195 const int meters_above_sea_level;
196 double *gc_latitude;
197 /*!
198 Reduces the geodetic latitude GD_LATITUDE of the Earth (which must be
199 given in radians) to a geocentric latitude respecting the altitude
200 of the location above/below sea level by using METERS_ABOVE_SEA_LEVEL.
201 Returns the calculated true Earth radius in meters for the location
202 and the calculated geocentric latitude in radians via the address of
203 GC_LATITUDE.
204 References:
205 * U.S. Government Printing Office, "The Astronomical Almanac" (AA), 1986.
206 */
207 {
208 auto double x1;
209 auto double x2;
210 auto double x3;
211 auto double cos_gd_latitude = cos (gd_latitude);
212 auto double cos2_gd_latitude;
213 auto double sin2_gd_latitude;
214
215
216 sin2_gd_latitude = sin (gd_latitude);
217 sin2_gd_latitude *= sin2_gd_latitude;
218 cos2_gd_latitude = cos_gd_latitude * cos_gd_latitude;
219 x1 =
220 (1.0 / sqrt (cos2_gd_latitude + FLATTENING_OF_EARTH * sin2_gd_latitude)) *
221 EQUATOR_EARTH_RADIUS;
222 x2 = x1 + meters_above_sea_level;
223 x3 = x1 * FLATTENING_OF_EARTH + meters_above_sea_level;
224 x1 = sqrt (x2 * x2 * cos2_gd_latitude + x3 * x3 * sin2_gd_latitude);
225 x3 = acos (x2 * cos_gd_latitude / x1);
226 if (SGN (gd_latitude) < 0)
227 x3 = -x3;
228 if (gc_latitude != (double *) NULL)
229 *gc_latitude = x3;
230
231 return (x1);
232 }
233
234
235
236 double
sun_rise_set(event,is_limited,day,month,year,coordinates)237 sun_rise_set (event, is_limited, day, month, year, coordinates)
238 const Aevent_enum event;
239 const Bool is_limited;
240 const int day;
241 const int month;
242 const int year;
243 Coor_struct *coordinates;
244 /*!
245 Returns the approximate timezone-based (local time) rise time of a definite
246 Sun altitude for the given date in hours plus fraction if EVENT is `RIse',
247 otherwise the approximate timezone-based (local time) set time of a
248 definite Sun altitude for the given date in hours plus fraction, resp.,
249 SPECIAL_VALUE if the calculated Sun rise/set time refers to a yesterdays
250 or tomorrows date, 2*SPECIAL_VALUE if the Sun is always below the queried
251 altitude on the given date, or 3*SPECIAL_VALUE if the Sun is always above
252 the queried altitude on the given date.
253 If the global `adjust_value' variable is set within a range of -90...+90
254 degree, its value is respected as reference altitude for rise/set-based
255 times and data, respectively, reference shadow length factor instead of
256 the pre-defined values otherwise used.
257 If data is computed that refers to rise/set times, and a rise/set time
258 does not occur on the given date, the ``error'' indicating return values
259 are the same as before but converted to seconds per hour for distinction
260 purposes.
261 If IS_LIMITED is TRUE and a calculated timezone-based rise/set time
262 does not occur on the given date, i.e. is prior or past the given date,
263 SPECIAL_VALUE is returned. If IS_LIMITED is FALSE and a calculated
264 timezone-based rise/set time is prior or past the given date (such can
265 always happen at high latitudes; locations that are within --or close to,
266 if you calculate the rise/set times for reference altitudes other than
267 the standard ones-- the solar Arctic Circle where the Sun occasionally
268 rises/sets more than once during a single day, respectively, never rises
269 or sets during a single day), the ``normalized'' --back-converted to the
270 given date-- timezone-based time in hours plus fraction is returned.
271 The longitude co-ordinates west of the zero meridian have a negative sign.
272 The longitude co-ordinates east of the zero meridian have a positive sign.
273 The latitude co-ordinates north of the equator have a positive sign.
274 The latitude co-ordinates south of the equator have a negative sign.
275 For negative numbers, all three of `*_deg', `*_min' and `*_sec' should
276 be negative. For example, the ISO 6709 co-ordinate `-202233+1100010'
277 (==+|-Latitude+|-Longitude) must be defined as `-20 -22 -33 +110 0 +10'.
278 Define which data of the Sun must be computed via the member variable
279 `the_mode'. Calculations for non-rise/set/transit/shadow length based
280 data are done for a line at a definite meridian expressed as a time value.
281 This time value is given by the member variable `time_zone_in_mins' for
282 the timezone that is respected at the given co-ordinate, and by the global
283 `time_offset' variable, which serves as the displacement value that is
284 added to midnight's time. If `time_offset' is set to zero, calculations
285 are made for Universal Time (UT/GMT). If `time_offset' has a positive
286 sign, UT/GMT calculations are made for meridians East of Greenwich,
287 otherwise for meridians West of Greenwich.
288 Here are all valid mode values and their description:
289 0 == This case is managed specially.
290 Hour plus fraction when the Sun's center passes the meridian opposite
291 the noon direction under or above the horizon near the lowermost
292 culmination point (meridian anti-transit, it is Sun's midnight).
293 1 == This case is managed specially.
294 Hour plus fraction when the Sun's center passes the meridian opposite
295 the midnight direction under or above the horizon near the uppermost
296 culmination point (meridian transit, it is Sun's noon).
297 2 == 0.0 degrees.
298 Hour plus fraction when the center of Sun's disk touches the horizon,
299 uncorrected for atmospheric refraction (mathematical horizon).
300 3 == approximately -16arc_minutes/60 => -0.26~ degrees.
301 Hour plus fraction when the Sun's upper limb touches the horizon,
302 uncorrected for atmospheric refraction (mathematical horizon).
303 4 == approximately -34arc_minutes/60 => -0.56~ degrees.
304 Hour plus fraction when the center of Sun's disk touches the horizon,
305 corrected for atmospheric refraction.
306 5 == approximately -50arc_minutes/60 => -0.83~ degrees.
307 Hour plus fraction when the Sun's upper limb touches the horizon,
308 corrected for atmospheric refraction (standard rise/set).
309 6 == -6.0 degrees.
310 Hour plus fraction of civil twilight (one can no longer read outside
311 without artificial illumination).
312 7 == -12.0 degrees.
313 Hour plus fraction of nautical twilight (navigation using a sea
314 horizon no longer possible).
315 8 == -15.0 degrees.
316 Hour plus fraction of amateur astronomical twilight (the sky is dark
317 enough for most astronomical observations).
318 9 == -18.0 degrees.
319 Hour plus fraction of astronomical twilight (the sky is completely
320 dark).
321 10 == This case is managed specially.
322 The topocentric apparent height of the Sun's center under or above
323 the sea-level horizon (elevation angle) in degrees plus fraction
324 that appears at midnight on the given date, corrected for atmospheric
325 refraction.
326 11 == This case is managed specially.
327 The Sun's topocentric apparent azimuth angle in degrees plus fraction
328 that appears at midnight on the given date.
329 12 == This case is managed specially.
330 The Sun's topocentric apparent declination angle in degrees plus
331 fraction that appears at midnight on the given date.
332 13 == This case is managed specially.
333 The Sun's topocentric apparent ecliptic longitude in degrees plus
334 fraction that appears at midnight on the given date.
335 14 == This case is managed specially.
336 The Sun's topocentric apparent right ascension angle in hours plus
337 fraction that appears at midnight on the given date.
338 15 == This case is managed specially.
339 The topocentric Sun/Earth distance in astronomical units plus fraction
340 that appears at midnight on the given date.
341 16 == This case is managed specially.
342 The Sun's topocentric apparent horizontal parallax in degrees plus
343 fraction that appears at midnight on the given date.
344 17 == This case is managed specially.
345 The Sun's topocentric apparent semidiameter in degrees plus fraction
346 that appears at midnight on the given date.
347 18 == The atmospheric refraction which is calculated for correcting
348 the Sun's topocentric true elevation angle to the Sun's topocentric
349 apparent elevation angle in degrees plus fraction that appears at
350 midnight on the given date.
351 19 == This case is managed specially.
352 The geocentric apparent height of the Sun's center under or above
353 the imaginary geocentric horizon (elevation angle) in degrees plus
354 fraction that appears at midnight on the given date.
355 20 == This case is managed specially.
356 The Sun's geocentric apparent azimuth angle in degrees plus fraction
357 that appears at midnight on the given date.
358 21 == This case is managed specially.
359 The Sun's geocentric apparent declination angle in degrees plus
360 fraction that appears at midnight on the given date.
361 22 == This case is managed specially.
362 The Sun's geocentric apparent ecliptic longitude in degrees plus
363 fraction that appears at midnight on the given date.
364 23 == This case is managed specially.
365 The Sun's geocentric apparent right ascension angle in hours plus
366 fraction that appears at midnight on the given date.
367 24 == This case is managed specially.
368 The geocentric Sun/Earth distance in astronomical units plus fraction
369 that appears at midnight on the given date.
370 25 == This case is managed specially.
371 The Sun's geocentric apparent horizontal parallax in degrees plus
372 fraction that appears at midnight on the given date.
373 26 == This case is managed specially.
374 The Sun's geocentric apparent semidiameter in degrees plus fraction
375 that appears at midnight on the given date.
376 27 == This case is managed specially.
377 The approximate ``Delta-T'' in seconds plus fraction, which is the
378 difference between Terrestrial Dynamical Time, TDT (formerly called
379 Ephemeris Time, ET), and Universal Time, UT, so `Delta_T = TDT - UT',
380 that appears at midnight on the given date.
381 28 == This case is managed specially.
382 The local sidereal time in hours plus fraction
383 that appears at midnight on the given date.
384 29 == This case is managed specially.
385 The globally given time offset value in hours plus fraction.
386 30 == This case is managed specially.
387 The local ``Julian Date'' in days plus fraction
388 that appears at midnight on the given date.
389 31 == This case is managed specially.
390 The approximate local ``Julian Ephemeris Date'' in days plus fraction
391 that appears at midnight on the given date.
392 32 == This case is managed specially.
393 The ``Equation of Time'' in hours plus fraction
394 that appears at midnight on the given date.
395 33 == This case is managed specially.
396 The difference between the Sun's and the Moon's topocentric
397 apparent elevation angles in degrees plus fraction that appears
398 at midnight on the given date.
399 34 == This case is managed specially.
400 The difference between the Sun's and the Moon's topocentric
401 apparent azimuth angles in degrees plus fraction that appears
402 at midnight on the given date.
403 35 == This case is managed specially.
404 The difference between the Sun's and the Moon's geocentric
405 apparent elevation angles in degrees plus fraction that appears
406 at midnight on the given date.
407 36 == This case is managed specially.
408 The difference between the Sun's and the Moon's geocentric
409 apparent azimuth angles in degrees plus fraction that appears
410 at midnight on the given date.
411 37 == This case is managed specially.
412 The difference between the Sun's and the Moon's topocentric
413 apparent elevation angles in degrees plus fraction that appears
414 at Sun's standard rise/set time.
415 38 == This case is managed specially.
416 The difference between the Sun's and the Moon's topocentric
417 apparent azimuth angles in degrees plus fraction that appears
418 at Sun's standard rise/set time.
419 39 == This case is managed specially.
420 The difference between the Sun's and the Moon's geocentric
421 apparent elevation angles in degrees plus fraction that appears
422 at Sun's standard rise/set time.
423 40 == This case is managed specially.
424 The difference between the Sun's and the Moon's geocentric
425 apparent azimuth angles in degrees plus fraction that appears
426 at Sun's standard rise/set time.
427 41 == This case is managed specially.
428 The difference between the Sun's and the Moon's meridian
429 anti-transit time in hours plus fraction that appears at
430 Sun's meridian anti-transit time.
431 42 == This case is managed specially.
432 The difference between the Sun's and the Moon's meridian
433 transit time in hours plus fraction that appears at Sun's
434 meridian transit time.
435 43 == This case is managed specially.
436 The difference between the Sun's and the Moon's standard rise/set time
437 in hours plus fraction that appears at Sun's standard rise/set time.
438 44 == This case is managed specially.
439 The Sun's topocentric apparent elevation angle (Sun's midnight height)
440 in degrees plus fraction that appears when the Sun's center
441 anti-transits the meridian under or above the sea-level horizon,
442 corrected for atmospheric refraction.
443 45 == This case is managed specially.
444 The Sun's topocentric apparent elevation angle (Sun's noon height)
445 in degrees plus fraction that appears when the Sun's center
446 transits the meridian under or above the sea-level horizon,
447 corrected for atmospheric refraction.
448 By the way, you can calculate the shadow length of an object
449 by using the formula:
450 Shadow_Length := Object_Height / tangent Sun_elevation_angle
451 and if the object length is 1, directly by:
452 Shadow_Length := cotangent Sun_elevation_angle
453 where the ``cotangent'' can be defined by:
454 cotangent X := cosine X / sine X
455 or by:
456 cotangent X := 1 / tangent X.
457 46 == This case is managed specially.
458 The Sun's topocentric apparent elevation angle under or above the
459 true horizon (the given height at the location) in degrees plus
460 fraction that appears at Sun's standard rise/set time,
461 corrected for atmospheric refraction.
462 47 == This case is managed specially.
463 The Sun's topocentric apparent azimuth angle in degrees plus fraction
464 that appears at Sun's standard rise/set time.
465 48 == This case is managed specially.
466 The Sun's geocentric apparent elevation angle (Sun's midnight height)
467 in degrees plus fraction that appears when the Sun's center
468 anti-transits the meridian under or above the sea-level horizon.
469 49 == This case is managed specially.
470 The Sun's geocentric apparent elevation angle (Sun's noon height)
471 in degrees plus fraction that appears when the Sun's center
472 transits the meridian under or above the sea-level horizon.
473 50 == This case is managed specially.
474 The Sun's geocentric apparent elevation angle under or above the
475 true horizon (the given height at the location) in degrees plus
476 fraction that appears at Sun's standard rise/set time.
477 51 == This case is managed specially.
478 The Sun's geocentric apparent azimuth angle in degrees plus fraction
479 that appears at Sun's standard rise/set time.
480 52 == This case is managed specially.
481 Hour plus fraction when the shadow length of an object at
482 true horizon (the given height at the location) is twice
483 the height of the object.
484 53 == This case is managed specially.
485 Hour plus fraction when the shadow length of an object at
486 true horizon (the given height at the location) is equal
487 the height of the object.
488 *** CAREFUL ***
489 Keep in mind that all values computed here are approximations only, based
490 on low-precision formula terms; and some values are computed recursively!
491 *** CAREFUL ***
492 State `meters_above_sea_level' if you know them for the co-ordinate,
493 otherwise 0. Specify `time_zone_in_mins' in +/- minutes MMM the
494 co-ordinate is before/behind UT/GMT, so that the calculated UT/GMT-based
495 rise/set times can be properly converted to local time times. For example,
496 if a co-ordinates timezone is `GMT-2' (==CEST), use `+120'.
497 References:
498 * "Practical Astronomy with your Calculator", Peter Duffet-Smith,
499 3rd edition. Cambridge University Press 1988, ISBN 0-521-35699-7.
500 * "Astronomical Formulae for Calculators", Jean Meeus, 4th ed,
501 Willmann-Bell 1988, ISBN 0-943396-22-0.
502 * "Astronomical Algorithms", Jean Meeus, 1st ed, Willmann-Bell 1991,
503 ISBN 0-943396-35-2.
504 * The WWW documents "How to compute planetary positions" and
505 "How to compute rise/set times and altitude above horizon"
506 by Paul Schlyter, <http://welcome.to/pausch> or
507 <http://hotel04.ausys.se/pausch>.
508 Btw., Keith Burnett publishes a lot of diverse --good quality--
509 astronomical sources and links at <http://www.xylem.demon.co.uk/kepler>.
510 */
511 {
512 auto double x1 = 0.0;
513 auto double x2;
514 auto double x3;
515 auto double x4;
516 auto double longitude;
517 auto double latitude;
518 auto double sin_latitude;
519 auto double cos_latitude;
520 auto double mjd;
521 auto double mjd_0ut;
522 auto double jc;
523 auto double local_sidereal_time;
524 auto double obliquity_of_ecliptic;
525 auto double sin_obliquity_of_ecliptic;
526 auto double cos_obliquity_of_ecliptic;
527 auto double argument_of_perihelion;
528 auto double mean_anomaly;
529 auto double eccentricity;
530 auto double anomaly_of_eccentric;
531 auto double divisor_term;
532 auto double true_anomaly;
533 auto double distance;
534 auto double ecliptic_longitude;
535 auto double xe;
536 auto double ye;
537 auto double ze;
538 auto double right_ascension;
539 auto double declination;
540 auto double meridian_transit_time = SECS_PER_DAY;
541 auto double the_time = SECS_PER_DAY;
542 auto double geodetic_height = 0.0;
543 register int n = 0;
544 register int tz_offset;
545 register int save_the_mode;
546 auto Bool iter_mt = FALSE;
547
548
549 /*
550 Reduce the given timezone value in minutes to a single day.
551 */
552 tz_offset = (int) (SGN (coordinates->time_zone_in_mins)
553 * (abs (coordinates->time_zone_in_mins) % MINS_PER_DAY));
554 x2 = DAY2HH (time_offset);
555 x4 = ((int) DAY2MM (time_offset)) % MINS_PER_HOUR;
556 switch (coordinates->the_mode)
557 {
558 case 27:
559 /*
560 Return ``Delta-T'' in seconds plus fraction.
561 */
562 return (delta_t (day, month, year, (int) x2, (int) x4));
563 case 29:
564 /*
565 Return the globally given time offset value in hours plus fraction.
566 */
567 return (DAY2HH (time_offset));
568 case 31:
569 /*
570 Return the local ``Julian Ephemeris Date'' in days plus fraction.
571 */
572 x1 = SS2DAY (delta_t (day, month, year, (int) x2, (int) x4));
573 /* Fallthrough. */
574 case 30:
575 /*
576 Return the local ``Julian Date'' in days plus fraction.
577 */
578 return (date2num (day, month, year) + (MIN_BCE_TO_1_CE - 1.5) + x1 +
579 time_offset);
580 case 32:
581 longitude = 0.0;
582 break;
583 case 33:
584 case 34:
585 case 35:
586 case 36:
587 save_the_mode = coordinates->the_mode;
588 switch (save_the_mode)
589 {
590 case 33:
591 /*
592 Calculate the Sun's topocentric apparent elevation angle
593 under or above the sea-level horizon in degrees plus
594 fraction, corrected for atmospheric refraction.
595 */
596 coordinates->the_mode = 10;
597 break;
598 case 34:
599 /*
600 Calculate the Sun's topocentric apparent azimuth angle
601 in degrees plus fraction.
602 */
603 coordinates->the_mode = 11;
604 break;
605 case 35:
606 /*
607 Calculate the Sun's geocentric apparent elevation angle
608 under or above the imaginary geocentric horizon
609 in degrees plus fraction.
610 */
611 coordinates->the_mode = 19;
612 break;
613 default:
614 /*
615 Calculate the Sun's geocentric apparent azimuth angle
616 in degrees plus fraction.
617 */
618 coordinates->the_mode = 20;
619 break;
620 }
621 the_time =
622 sun_rise_set (event, is_limited, day, month, year, coordinates);
623 switch (save_the_mode)
624 {
625 case 33:
626 case 35:
627 if (save_the_mode == 33)
628 /*
629 Calculate the difference between the Sun's and the Moon's
630 topocentric apparent elevation angles in degrees plus
631 fraction that appears at midnight on the given date.
632 */
633 coordinates->the_mode = 10;
634 else
635 /*
636 Calculate the difference between the Sun's and the Moon's
637 geocentric apparent elevation angles in degrees plus
638 fraction that appears at midnight on the given date.
639 */
640 coordinates->the_mode = 23;
641 the_time -=
642 internal_moon_rise_set (event, day, month, year, coordinates);
643 break;
644 default:
645 if (save_the_mode == 34)
646 /*
647 Calculate the difference between the Sun's and the Moon's
648 topocentric apparent azimuth angles in degrees plus
649 fraction that appears at midnight on the given date.
650 */
651 coordinates->the_mode = 11;
652 else
653 /*
654 Calculate the difference between the Sun's and the Moon's
655 geocentric apparent azimuth angles in degrees plus
656 fraction that appears at midnight on the given date.
657 */
658 coordinates->the_mode = 24;
659 x4 =
660 the_time - internal_moon_rise_set (event, day, month, year,
661 coordinates);
662 /*
663 Reduce the resulting angle to -180.0...+180.0 degrees.
664 */
665 if (abs (x4) <= DEGS_PER_12_HOURS)
666 the_time = x4;
667 else if (x4 > DEGS_PER_12_HOURS)
668 the_time = -DEGS_PER_24_HOURS + x4;
669 else
670 the_time = DEGS_PER_24_HOURS + x4;
671 }
672 coordinates->the_mode = save_the_mode;
673 return (the_time);
674 case 41:
675 case 42:
676 case 43:
677 case 44:
678 case 45:
679 case 48:
680 case 49:
681 x1 = time_offset;
682 save_the_mode = coordinates->the_mode;
683 switch (save_the_mode)
684 {
685 case 43:
686 /*
687 Perform calculations for Sun's standard rise/set time.
688 */
689 coordinates->the_mode = 5;
690 break;
691 case 41:
692 case 44:
693 case 48:
694 /*
695 Perform calculations for Sun's astronomical midnight.
696 */
697 coordinates->the_mode = 0;
698 break;
699 default:
700 /*
701 Perform calculations for Sun's astronomical noon.
702 */
703 coordinates->the_mode = 1;
704 }
705 time_offset =
706 sun_rise_set (event, is_limited, day, month, year, coordinates);
707 if (time_offset > SPECIAL_VALUE)
708 {
709 time_offset = HH2DAY (time_offset);
710 if (save_the_mode <= 42)
711 {
712 if (save_the_mode == 41)
713 /*
714 Calculate the difference between the Sun's and the Moon's
715 meridian anti-transit time in hours plus fraction that
716 appears at Sun's meridian anti-transit time.
717 */
718 coordinates->the_mode = 0;
719 else
720 /*
721 Calculate the difference between the Sun's and the Moon's
722 meridian transit time in hours plus fraction that appears
723 at Sun's meridian transit time.
724 */
725 coordinates->the_mode = 1;
726 the_time = moon_rise_set (event, day, month, year, coordinates);
727 if (the_time > SPECIAL_VALUE)
728 the_time = DAY2HH (time_offset) - the_time;
729 else
730 the_time = HH2SS (the_time);
731 }
732 else if (save_the_mode == 43)
733 {
734 /*
735 Calculate the difference between the Sun's and the Moon's
736 standard rise/set time in hours plus fraction that appears
737 at Sun's standard rise/set time.
738 */
739 coordinates->the_mode = 5;
740 the_time = moon_rise_set (event, day, month, year, coordinates);
741 if (the_time > SPECIAL_VALUE)
742 the_time = DAY2HH (time_offset) - the_time;
743 else
744 the_time = HH2SS (the_time);
745 }
746 else
747 {
748 switch (save_the_mode)
749 {
750 case 44:
751 case 45:
752 /*
753 Calculate the Sun's topocentric apparent elevation angle
754 under or above the sea-level horizon in degrees plus
755 fraction, corrected for atmospheric refraction.
756 */
757 coordinates->the_mode = 10;
758 break;
759 default:
760 /*
761 Calculate the Sun's geocentric apparent elevation angle
762 under or above the imaginary geocentric horizon in
763 degrees plus fraction.
764 */
765 coordinates->the_mode = 19;
766 }
767 the_time =
768 sun_rise_set (event, is_limited, day, month, year,
769 coordinates);
770 }
771 }
772 else
773 the_time = HH2SS (time_offset);
774 coordinates->the_mode = save_the_mode;
775 time_offset = x1;
776 return (the_time);
777 case 37:
778 case 38:
779 case 39:
780 case 40:
781 case 46:
782 case 47:
783 case 50:
784 case 51:
785 x1 = time_offset;
786 save_the_mode = coordinates->the_mode;
787 coordinates->the_mode = 5;
788 time_offset =
789 sun_rise_set (event, is_limited, day, month, year, coordinates);
790 if (time_offset > SPECIAL_VALUE)
791 {
792 time_offset = HH2DAY (time_offset);
793 switch (save_the_mode)
794 {
795 case 37:
796 case 46:
797 /*
798 Calculate the Sun's topocentric apparent elevation angle
799 under or above the true horizon in degrees plus fraction
800 that appears at Sun's standard rise/set time,
801 corrected for atmospheric refraction.
802 */
803 coordinates->the_mode = 10;
804 break;
805 case 38:
806 case 47:
807 /*
808 Calculate the Sun's topocentric apparent azimuth angle
809 in degrees plus fraction that appears at Sun's
810 standard rise/set time.
811 */
812 coordinates->the_mode = 11;
813 break;
814 case 39:
815 case 50:
816 /*
817 Calculate the Sun's geocentric apparent elevation angle
818 under or above the true horizon in degrees plus fraction
819 that appears at Sun's standard rise/set time.
820 */
821 coordinates->the_mode = 19;
822 break;
823 default:
824 /*
825 Calculate the Sun's geocentric apparent azimuth angle
826 in degrees plus fraction that appears at Sun's
827 standard rise/set time.
828 */
829 coordinates->the_mode = 20;
830 }
831 the_time =
832 sun_rise_set (event, is_limited, day, month, year, coordinates);
833 if ((save_the_mode >= 37) && (save_the_mode <= 40))
834 switch (save_the_mode)
835 {
836 case 37:
837 case 39:
838 if (save_the_mode == 37)
839 /*
840 Calculate the difference between the Sun's and the Moon's
841 topocentric apparent elevation angles in degrees plus
842 fraction that appears at Sun's standard rise/set time.
843 */
844 coordinates->the_mode = 10;
845 else
846 /*
847 Calculate the difference between the Sun's and the Moon's
848 geocentric apparent elevation angles in degrees plus
849 fraction that appears at Sun's standard rise/set time.
850 */
851 coordinates->the_mode = 23;
852 the_time -=
853 internal_moon_rise_set (event, day, month, year,
854 coordinates);
855 break;
856 default:
857 if (save_the_mode == 38)
858 /*
859 Calculate the difference between the Sun's and the Moon's
860 topocentric apparent azimuth angles in degrees plus
861 fraction that appears at Sun's standard rise/set time.
862 */
863 coordinates->the_mode = 11;
864 else
865 /*
866 Calculate the difference between the Sun's and the Moon's
867 geocentric apparent azimuth angles in degrees plus
868 fraction that appears at Sun's standard rise/set time.
869 */
870 coordinates->the_mode = 24;
871 x4 =
872 the_time - internal_moon_rise_set (event, day, month, year,
873 coordinates);
874 /*
875 Reduce the resulting angle to -180.0...+180.0 degrees.
876 */
877 if (abs (x4) <= DEGS_PER_12_HOURS)
878 the_time = x4;
879 else if (x4 > DEGS_PER_12_HOURS)
880 the_time = -DEGS_PER_24_HOURS + x4;
881 else
882 the_time = DEGS_PER_24_HOURS + x4;
883 }
884 }
885 else
886 the_time = HH2SS (time_offset);
887 coordinates->the_mode = save_the_mode;
888 time_offset = x1;
889 return (the_time);
890 case 0:
891 case 1:
892 iter_mt = TRUE;
893 /* Fallthrough. */
894 default:
895 longitude = coordinates->lon_deg
896 + MM2DEG (coordinates->lon_min) + SS2DEG (coordinates->lon_sec);
897 }
898 latitude = TORAD (coordinates->lat_deg
899 + MM2DEG (coordinates->lat_min)
900 + SS2DEG (coordinates->lat_sec));
901 if (coordinates->meters_above_sea_level > 0)
902 geodetic_height =
903 TORAD (0.0347 * sqrt (coordinates->meters_above_sea_level));
904 /*
905 Epoch starts 01-Jan-2000 12:00:00 UT, convert UT to TDT.
906 */
907 x1 = SS2DAY (delta_t (day, month, year, 0, 0));
908 mjd = mjd_0ut = date2num (day, month, year) + x1 - 730122.5;
909 /*
910 Respect timezone and time displacement value for
911 Sun's non-rise/set/transit/shadow length based modes.
912 */
913 if ((coordinates->the_mode >= 10) && (coordinates->the_mode <= 36))
914 mjd += (time_offset - MM2DAY (tz_offset)
915 + SS2DAY (delta_t (day, month, year, (int) x2, (int) x4)) - x1);
916 do
917 {
918 LABEL_sun_iter_mt:
919 local_sidereal_time = FIXANGLE (100.46061837
920 + 0.9856472275 * mjd + longitude);
921 jc = mjd / 36525.0;
922 obliquity_of_ecliptic = TORAD (FIXANGLE (23.43929111
923 - (0.013004167 +
924 (0.00000016389 -
925 0.0000005036 * jc) * jc) *
926 jc));
927 argument_of_perihelion =
928 TORAD (FIXANGLE (282.93735 + (0.71953 + (0.00046 * jc)) * jc));
929 mean_anomaly =
930 TORAD (FIXANGLE
931 (357.52910 +
932 (35999.05030 - (0.0001559 + 0.00000048 * jc) * jc) * jc));
933 eccentricity = 0.016708617 - (0.000042037 + (0.0000001236 * jc)) * jc;
934 x2 = mean_anomaly + eccentricity
935 * sin (mean_anomaly) * (1.0 + eccentricity * cos (mean_anomaly));
936 do
937 {
938 x1 = x2;
939 x2 = x1 - (x1 - eccentricity * sin (x1) - mean_anomaly)
940 / (1.0 - eccentricity * cos (x1));
941 }
942 while (abs (abs (x1) - abs (x2)) > 0.000001);
943 anomaly_of_eccentric = x2;
944 divisor_term = cos (anomaly_of_eccentric) - eccentricity;
945 true_anomaly =
946 sqrt (1.0 - eccentricity * eccentricity) * sin (anomaly_of_eccentric);
947 /*
948 Calculate the geocentric Sun/Earth distance.
949 */
950 distance =
951 sqrt (divisor_term * divisor_term + true_anomaly * true_anomaly);
952 switch (coordinates->the_mode)
953 {
954 case 24:
955 /*
956 Return the geocentric Sun/Earth distance
957 in astronomical units plus fraction.
958 */
959 return (distance);
960 case 25:
961 /*
962 Return the Sun's geocentric apparent horizontal parallax
963 in degrees plus fraction.
964 */
965 return (SS2DEG (8.794) * distance);
966 case 26:
967 /*
968 Return the Sun's geocentric apparent semidiameter
969 in degrees plus fraction.
970 */
971 return (SS2DEG (959.63) * distance);
972 default:
973 ; /* Void, nothing to do! */
974 }
975 true_anomaly = atan (true_anomaly / divisor_term);
976 if (divisor_term < 0.0)
977 true_anomaly += MY_PI;
978 ecliptic_longitude =
979 TORAD (FIXANGLE
980 (TODEG (true_anomaly + argument_of_perihelion) + jc - 0.00569 -
981 0.00479 *
982 sin (TORAD (124.90 - (1934.134 - (0.002063 * jc)) * jc))));
983 if (coordinates->the_mode == 22)
984 /*
985 Return the Sun's geocentric apparent ecliptic longitude
986 in degrees plus fraction.
987 */
988 return (TODEG (ecliptic_longitude));
989 /*
990 Calculate geocentric rectangular equatorial co-ordinates.
991 */
992 cos_obliquity_of_ecliptic = cos (obliquity_of_ecliptic);
993 x1 = sin (ecliptic_longitude) * distance;
994 xe = cos (ecliptic_longitude) * distance;
995 ye = cos_obliquity_of_ecliptic * x1;
996 /*
997 Calculate the Sun's geocentric apparent right ascension angle.
998 */
999 right_ascension = FIXANGLE (TODEG (my_atan2 (ye, xe)));
1000 if (coordinates->the_mode == 23)
1001 /*
1002 Return the Sun's geocentric apparent right ascension angle
1003 in hours plus fraction.
1004 */
1005 return (DEG2HH (right_ascension));
1006 x4 = meridian_transit_time;
1007 meridian_transit_time =
1008 FIXANGLE (right_ascension - local_sidereal_time);
1009 if (coordinates->the_mode == 32)
1010 /*
1011 Return the ``Equation of Time'' in hours plus fraction.
1012 */
1013 return (DEG2HH (DEGS_PER_12_HOURS - meridian_transit_time));
1014 if (iter_mt)
1015 {
1016 /*
1017 Calculate the meridian transit time (Sun's astronomical noon).
1018 */
1019 mjd = mjd_0ut + DEG2DAY (meridian_transit_time);
1020 /*
1021 Exit iteration loop if the actually computed value is less
1022 than 1/1000.0 hours (== 3.6 seconds [0.015 degrees]) smaller
1023 the previously computed value. A ``next'' iteration loop
1024 would not change the result significantly -- the only change
1025 would be by the order of fractions of seconds.
1026 Well, such a ``huge :-)'' epsilon value as well as a smaller
1027 epsilon value can be really critical if the iterated time is
1028 very near midnight, but what's perfect?
1029 */
1030 if (abs (abs (meridian_transit_time) - abs (x4)) > 0.015)
1031 goto LABEL_sun_iter_mt;
1032 the_time = meridian_transit_time;
1033 break;
1034 }
1035 sin_obliquity_of_ecliptic = sin (obliquity_of_ecliptic);
1036 ze = sin_obliquity_of_ecliptic * x1;
1037 /*
1038 Calculate the Sun's geocentric apparent declination angle.
1039 */
1040 declination = asin (ze / distance);
1041 if (coordinates->the_mode == 21)
1042 /*
1043 Return the Sun's geocentric apparent declination angle
1044 in degrees plus fraction.
1045 */
1046 return (TODEG (declination));
1047 sin_latitude = sin (latitude);
1048 cos_latitude = cos (latitude);
1049 if ((coordinates->the_mode >= 10) && (coordinates->the_mode <= 28))
1050 {
1051 /*
1052 Calculate the local sidereal time of moment.
1053 */
1054 local_sidereal_time = FIXANGLE (local_sidereal_time
1055 + DAY2DEG (time_offset -
1056 MM2DAY (tz_offset)));
1057 switch (coordinates->the_mode)
1058 {
1059 case 19:
1060 case 20:
1061 x1 = TORAD (FIXANGLE (local_sidereal_time - right_ascension));
1062 if (coordinates->the_mode == 19)
1063 /*
1064 Return the geocentric apparent height of the Sun's center
1065 under or above the imaginary geocentric horizon
1066 (elevation angle) in degrees plus fraction.
1067 */
1068 return (TODEG (asin (sin_latitude * sin (declination)
1069 +
1070 cos_latitude * cos (declination) *
1071 cos (x1))));
1072 /*
1073 Return the Sun's geocentric apparent azimuth angle
1074 in degrees plus fraction.
1075 */
1076 return (FIXANGLE (TODEG (my_atan2 (sin (x1),
1077 sin_latitude * cos (x1) -
1078 cos_latitude *
1079 tan (declination)) +
1080 MY_PI)));
1081 case 28:
1082 /*
1083 Return the local sidereal time of moment
1084 in hours plus fraction.
1085 */
1086 return (DEG2HH (local_sidereal_time));
1087 default:
1088 /*
1089 Calculate topocentric rectangular equatorial co-ordinates.
1090 */
1091 x4 = gd_latitude2gc_latitude (latitude,
1092 (coordinates->
1093 meters_above_sea_level >
1094 0) ? coordinates->
1095 meters_above_sea_level : 0,
1096 &x3) / EQUATOR_EARTH_RADIUS;
1097 x4 /= (distance * EARTH_RADII_PER_AU);
1098 x1 = cos (x3) * x4;
1099 x2 = TORAD (local_sidereal_time);
1100 xe -= (cos (x2) * x1);
1101 ye -= (sin (x2) * x1);
1102 ze -= (sin (x3) * x4);
1103 }
1104 /*
1105 Calculate the topocentric Sun/Earth distance in astronomical
1106 units plus fraction.
1107 */
1108 distance = sqrt (xe * xe + ye * ye + ze * ze);
1109 switch (coordinates->the_mode)
1110 {
1111 case 15:
1112 /*
1113 Return the topocentric Sun/Earth distance
1114 in astronomical units plus fraction.
1115 */
1116 return (distance);
1117 case 16:
1118 /*
1119 Return the Sun's topocentric apparent horizontal parallax
1120 in degrees plus fraction.
1121 */
1122 return (SS2DEG (8.794) * distance);
1123 case 17:
1124 /*
1125 Return the Sun's topocentric apparent semidiameter
1126 in degrees plus fraction.
1127 */
1128 return (SS2DEG (959.63) * distance);
1129 default:
1130 ; /* Void, nothing to do! */
1131 }
1132 /*
1133 Calculate the topocentric right ascension angle,
1134 uncorrected for atmospheric refraction.
1135 */
1136 right_ascension = FIXANGLE (TODEG (my_atan2 (ye, xe)));
1137 /*
1138 Calculate the topocentric declination angle,
1139 uncorrected for atmospheric refraction.
1140 */
1141 declination = asin (ze / distance);
1142 if (coordinates->the_mode == 13)
1143 /*
1144 Return the Sun's topocentric apparent ecliptic longitude
1145 in degrees plus fraction.
1146 */
1147 return (FIXANGLE
1148 (TODEG
1149 (my_atan2
1150 (cos_obliquity_of_ecliptic *
1151 sin (TORAD (right_ascension)) +
1152 sin_obliquity_of_ecliptic * tan (declination),
1153 cos (TORAD (right_ascension))))));
1154 /*
1155 Calculate the topocentric apparent height of the Sun's center
1156 under or above the sea-level horizon (elevation angle),
1157 corrected for atmospheric refraction.
1158 */
1159 x1 = TORAD (FIXANGLE (local_sidereal_time - right_ascension));
1160 x3 = cos (x1);
1161 x2 =
1162 asin (sin_latitude * sin (declination) +
1163 cos_latitude * cos (declination) * x3);
1164 /*
1165 Calculate the atmospheric refraction for correcting the Sun's
1166 topocentric true elevation angle to the Sun's topocentric
1167 apparent elevation angle.
1168 */
1169 x4 = atmospheric_refraction (x2, atm_pressure, atm_temperature);
1170 if (coordinates->the_mode == 18)
1171 /*
1172 Return the atmospheric refraction for correcting the Sun's
1173 topocentric true elevation angle to the Sun's topocentric
1174 apparent elevation angle in degrees plus fraction.
1175 */
1176 return (TODEG (x4));
1177 x2 += x4;
1178 if (x2 > MY_HALF_PI)
1179 x2 = MY_PI - x2;
1180 if (coordinates->the_mode == 10)
1181 /*
1182 Return the topocentric apparent height of the Sun's center
1183 under or above the sea-level horizon (elevation angle)
1184 in degrees plus fraction, corrected for atmospheric
1185 refraction.
1186 */
1187 return (TODEG (x2));
1188 /*
1189 Calculate the Sun's topocentric apparent azimuth angle
1190 in degrees plus fraction.
1191 */
1192 x3 = FIXANGLE (TODEG (my_atan2 (sin (x1),
1193 sin_latitude * x3 -
1194 cos_latitude * tan (declination)) +
1195 MY_PI));
1196 if (coordinates->the_mode == 11)
1197 return (x3);
1198 /*
1199 Reconvert the topocentric apparent azimuth angle and the
1200 topocentric apparent elevation angle to topocentric
1201 apparent right ascension angle and topocentric apparent
1202 declination angle.
1203 */
1204 x3 = TORAD (x3);
1205 xe = cos (x2);
1206 ye = sin (x2);
1207 ze = cos (x3);
1208 x1 = xe * ze;
1209 x2 = -xe * sin (x3);
1210 x3 = cos_latitude * ye - sin_latitude * x1;
1211 if (coordinates->the_mode == 12)
1212 /*
1213 Return the Sun's topocentric apparent declination angle
1214 in degrees plus fraction.
1215 */
1216 return (TODEG (asin (sin_latitude * ye + cos_latitude * x1)));
1217 /*
1218 Return the Sun's topocentric apparent right ascension angle
1219 in hours plus fraction.
1220 */
1221 return (DEG2HH
1222 (FIXANGLE
1223 (local_sidereal_time - TODEG (my_atan2 (x2, x3)))));
1224 }
1225 /*
1226 Select the proper reference altitude
1227 according to the requested rise/set mode of the Sun.
1228 */
1229 if (adjust_value != DEGS_PER_24_HOURS)
1230 switch (coordinates->the_mode)
1231 {
1232 case 52:
1233 case 53:
1234 x1 = adjust_value;
1235 break;
1236 case 3:
1237 case 5:
1238 /*
1239 Sun's upper limb based altitude.
1240 */
1241 x1 = TORAD (adjust_value - (SS2DEG (959.63) * distance));
1242 break;
1243 default:
1244 /*
1245 All other Sun's center based altitudes.
1246 */
1247 x1 = TORAD (adjust_value);
1248 }
1249 else
1250 switch (coordinates->the_mode)
1251 {
1252 case 2:
1253 case 4:
1254 /*
1255 Sun's center based altitude (0.0 degree == 0.0 radians).
1256 */
1257 x1 = 0.0;
1258 break;
1259 case 3:
1260 case 5:
1261 /*
1262 Sun's upper limb based altitude.
1263 */
1264 x1 = TORAD (SS2DEG (-959.63) * distance);
1265 break;
1266 case 6:
1267 /*
1268 Civil twilight altitude (-6.0 degree == -0.104719755119659771 radians).
1269 */
1270 x1 = -0.104719755119659771;
1271 break;
1272 case 7:
1273 /*
1274 Nautical twilight altitude (-12.0 degree == -0.209439510239319541 radians).
1275 */
1276 x1 = -0.209439510239319541;
1277 break;
1278 case 8:
1279 /*
1280 Amateur astronomical twilight altitude (-15.0 degree == -0.261799387799149426 radians).
1281 */
1282 x1 = -0.261799387799149426;
1283 break;
1284 case 9:
1285 /*
1286 Astronomical twilight altitude (-18.0 degree == -0.314159265358979312 radians).
1287 */
1288 x1 = -0.314159265358979312;
1289 break;
1290 case 52:
1291 /*
1292 Shadow length at true horizon is twice the height of object.
1293 */
1294 x1 = 2.0;
1295 break;
1296 case 53:
1297 /*
1298 Shadow length at true horizon is equal the height of object.
1299 */
1300 x1 = 1.0;
1301 }
1302 switch (coordinates->the_mode)
1303 {
1304 case 52:
1305 case 53:
1306 x2 = latitude - declination;
1307 if (SGN (latitude) < 0)
1308 {
1309 declination = -declination;
1310 if (SGN (x2) < 0)
1311 x2 = -x2;
1312 }
1313 x2 = (sin (my_acot (x1 + tan (x2)) - geodetic_height)
1314 - abs (sin_latitude) * sin (declination))
1315 / (abs (cos_latitude) * cos (declination));
1316 break;
1317 case 4:
1318 case 5:
1319 /*
1320 Correct the true altitude giving the apparent altitude.
1321 */
1322 if (adjust_value != DEGS_PER_24_HOURS
1323 || atm_pressure != DEFAULT_PRESSURE
1324 || atm_temperature != DEFAULT_TEMPERATURE)
1325 x2 = atmospheric_refraction (x1, atm_pressure, atm_temperature);
1326 else
1327 /*
1328 Standard refraction of 34arc_minutes/60 == 0.00989019909463453388 radians.
1329 */
1330 x2 = 0.00989019909463453388;
1331 x1 -= x2;
1332 /* Fallthrough. */
1333 default:
1334 x2 = (sin (x1 - geodetic_height)
1335 - sin_latitude * sin (declination))
1336 / (cos_latitude * cos (declination));
1337 }
1338 if (abs (x2) > 1.0)
1339 {
1340 if (n)
1341 {
1342 if (x2 > 0.0)
1343 /*
1344 No rise/set for the given date
1345 because the Sun is always below the queried altitude.
1346 */
1347 return (2 * SPECIAL_VALUE);
1348 /*
1349 No rise/set for the given date
1350 because the Sun is always above the queried altitude.
1351 */
1352 return (3 * SPECIAL_VALUE);
1353 }
1354 /*
1355 Try a next iteration step by using the Sun's meridian transit time.
1356 */
1357 mjd = mjd_0ut + DEG2DAY (meridian_transit_time);
1358 n++;
1359 continue;
1360 }
1361 /*
1362 Calculate the Sun's rise/set time.
1363 */
1364 x2 = TODEG (acos (x2));
1365 x1 = the_time;
1366 if (event == RIse)
1367 the_time = meridian_transit_time - x2;
1368 else
1369 the_time = meridian_transit_time + x2;
1370 mjd = mjd_0ut + DEG2DAY (the_time);
1371 /*
1372 Exit iteration loop if the actually computed value is less
1373 than 1/1000.0 hours (== 3.6 seconds [0.015 degrees]) smaller
1374 the previously computed value. A ``next'' iteration loop
1375 would not change the result significantly -- the only change
1376 would be by the order of fractions of seconds.
1377 Well, such a ``huge :-)'' epsilon value as well as a smaller
1378 epsilon value can be really critical if the iterated time is
1379 very near midnight, but what's perfect?
1380 */
1381 }
1382 while (abs (abs (the_time) - abs (x1)) > 0.015);
1383 /*
1384 Convert sun's UT/GMT-based rise/set/transit position to timezone-based position.
1385 */
1386 the_time += HH2DEG (MM2HH (tz_offset));
1387 if (!coordinates->the_mode)
1388 {
1389 /*
1390 Calculate the meridian anti-transit time (Sun's astronomical midnight).
1391 */
1392 x1 = time_offset;
1393 coordinates->the_mode = 32;
1394 time_offset = DEG2DAY (the_time);
1395 x2 = sun_rise_set (event, is_limited, day, month, year, coordinates);
1396 if (the_time >= DEGS_PER_12_HOURS)
1397 the_time -= DEGS_PER_12_HOURS;
1398 else
1399 the_time += DEGS_PER_12_HOURS;
1400 time_offset = DEG2DAY (the_time);
1401 x3 = sun_rise_set (event, is_limited, day, month, year, coordinates);
1402 /*
1403 Calculation by using the ``Equation of Time'' is precise
1404 enough here in this low-precision function!
1405 */
1406 the_time += HH2DEG (x2 - x3);
1407 coordinates->the_mode = 0;
1408 time_offset = x1;
1409 }
1410 if (is_limited && (the_time < 0.0 || the_time >= DEGS_PER_24_HOURS))
1411 /*
1412 No rise/set of the Sun for the given date due to limitation to the
1413 given date, but on the day prior or after. Gcal ignores such an
1414 event because it would place a ``yesterday's'' or ``tomorrow's''
1415 rise/set time to ``today's'' date, and not to a ``yesterday's''
1416 or ``tomorrow's'' date as it would be the astronomical *TRUTH*!
1417 Gcal must behave like this as long as it is still unable to put
1418 more than a single rise/set time into output for each given
1419 special text in a resource file that enables such a rise/set time.
1420 */
1421 return (SPECIAL_VALUE);
1422 /*
1423 Correct the calculated (position) data to circle range.
1424 */
1425 while (the_time >= DEGS_PER_24_HOURS)
1426 the_time -= DEGS_PER_24_HOURS;
1427 while (the_time < 0.0)
1428 the_time += DEGS_PER_24_HOURS;
1429
1430 /*
1431 Return the rise/set/transit times of the Sun in hours plus fraction.
1432 */
1433 return (DEG2HH (the_time));
1434 }
1435
1436
1437
1438 double
moon_rise_set(event,day,month,year,coordinates)1439 moon_rise_set (event, day, month, year, coordinates)
1440 const Aevent_enum event;
1441 int day;
1442 int month;
1443 int year;
1444 Coor_struct *coordinates;
1445 /*!
1446 Driver/wrapper function to calculate timezone-based (local time)
1447 rise/set times of the Moon correctly.
1448 See the `internal_moon_rise_set()' function for more details!
1449 */
1450 {
1451 auto double the_time;
1452 register int tz_offset;
1453
1454
1455 /*
1456 Reduce the given timezone value to a single day.
1457 */
1458 tz_offset = SGN (coordinates->time_zone_in_mins)
1459 * (abs (coordinates->time_zone_in_mins) % MINS_PER_DAY);
1460 /*
1461 Calculate a UTC/GMT-based rise/set time or other data.
1462 */
1463 the_time = internal_moon_rise_set (event, day, month, year, coordinates);
1464 if (tz_offset && (coordinates->the_mode < 6))
1465 {
1466 /*
1467 Convert the calculated UTC/GMT-based rise/set time to local time.
1468 */
1469 if (the_time > SPECIAL_VALUE)
1470 {
1471 the_time += MM2HH (tz_offset);
1472 if (the_time < 0.0 || the_time >= HOURS_PER_DAY)
1473 {
1474 if (coordinates->time_zone_in_mins < 0)
1475 (void) next_date (&day, &month, &year);
1476 else
1477 (void) prev_date (&day, &month, &year);
1478 the_time =
1479 internal_moon_rise_set (event, day, month, year, coordinates);
1480 if (the_time <= SPECIAL_VALUE)
1481 return (the_time);
1482 the_time += MM2HH (tz_offset);
1483 if ((the_time >= 0.0) && (the_time < HOURS_PER_DAY))
1484 return (SPECIAL_VALUE);
1485 /*
1486 Correct the calculated time to day range.
1487 */
1488 while (the_time >= HOURS_PER_DAY)
1489 the_time -= HOURS_PER_DAY;
1490 while (the_time < 0.0)
1491 the_time += HOURS_PER_DAY;
1492 }
1493 }
1494 else if (the_time == SPECIAL_VALUE)
1495 {
1496 if (coordinates->time_zone_in_mins < 0)
1497 (void) next_date (&day, &month, &year);
1498 else
1499 (void) prev_date (&day, &month, &year);
1500 the_time =
1501 internal_moon_rise_set (event, day, month, year, coordinates);
1502 if (the_time > SPECIAL_VALUE)
1503 {
1504 the_time += MM2HH (tz_offset);
1505 if ((the_time > 0.0) && (the_time <= HOURS_PER_DAY))
1506 return (SPECIAL_VALUE);
1507 /*
1508 Correct the calculated time to day range.
1509 */
1510 while (the_time >= HOURS_PER_DAY)
1511 the_time -= HOURS_PER_DAY;
1512 while (the_time < 0.0)
1513 the_time += HOURS_PER_DAY;
1514 }
1515 }
1516 }
1517
1518 return (the_time);
1519 }
1520
1521
1522
1523 static double
internal_moon_rise_set(event,day,month,year,coordinates)1524 internal_moon_rise_set (event, day, month, year, coordinates)
1525 const Aevent_enum event;
1526 int day;
1527 int month;
1528 int year;
1529 Coor_struct *coordinates;
1530 /*!
1531 Returns the approximate UT/GMT-based rise time of a definite Moon altitude
1532 for the given date in hours plus fraction if EVENT is `RIse', otherwise
1533 the approximate UT/GMT-based set time of a definite Moon altitude for the
1534 given date in hours plus fraction, resp., SPECIAL_VALUE if the Moon
1535 does not rise/set on the given date, 2*SPECIAL_VALUE if the Moon is
1536 always below the queried altitude on the given date, or 3*SPECIAL_VALUE
1537 if the Moon is always above the queried altitude on the given date.
1538 If the global `adjust_value' variable is set within a range of -90...+90
1539 degree, its value is respected as reference altitude for rise/set-based
1540 times and data instead of the pre-defined values otherwise used.
1541 If data is computed that refers to rise/set/transit times, and a
1542 rise/set/transit time does not occur on the given date, the ``error''
1543 indicating return values are the same as before but converted to
1544 seconds per hour for distinction purposes.
1545 The longitude co-ordinates west of the zero meridian have a negative sign.
1546 The longitude co-ordinates east of the zero meridian have a positive sign.
1547 The latitude co-ordinates north of the equator have a positive sign.
1548 The latitude co-ordinates south of the equator have a negative sign.
1549 For negative numbers, all three of `*_deg', `*_min' and `*_sec' should
1550 be negative. For example, the ISO 6709 co-ordinate `-202233+1100010'
1551 (==+|-Latitude+|-Longitude) must be defined as `-20 -22 -33 +110 0 +10'.
1552 Define which data of the Moon must be computed via the member variable
1553 `the_mode'. Calculations for non-rise/set/transit based data are done
1554 for a line at a definite meridian expressed as a time value. This time
1555 value is given by the member variable `time_zone_in_mins' for the
1556 timezone that is respected at the given co-ordinate, and by the global
1557 `time_offset' variable, which serves as the displacement value that is
1558 added to midnight's time. If `time_offset' is set to zero, calculations
1559 are made for Universal Time (UT/GMT). If `time_offset' has a positive
1560 sign, UT/GMT calculations are made for meridians East of Greenwich,
1561 otherwise for meridians West of Greenwich.
1562 Here are all valid mode values and their description:
1563 0 == This case is managed specially.
1564 Hour plus fraction when the Moon's center passes the meridian opposite
1565 the noon direction under or above the horizon near the lowermost
1566 culmination point (meridian anti-transit, it is Moon's midnight).
1567 1 == This case is managed specially.
1568 Hour plus fraction when the Moon's center passes the meridian opposite
1569 the midnight direction under or above the horizon near the uppermost
1570 culmination point (meridian transit, it is Moon's noon).
1571 2 == approximately 54...61arc_minutes/60 ~=> +0.9000...+1.0163~ degrees.
1572 Hour plus fraction when the center of Moon's disk touches the horizon,
1573 uncorrected for atmospheric refraction (mathematical horizon).
1574 3 == approximately 39...44arc_minutes/60 ~=> +0.6500...+0.7333~ degrees.
1575 Hour plus fraction when the Moon's upper limb touches the horizon,
1576 uncorrected for atmospheric refraction (mathematical horizon).
1577 4 == approximately 20...27arc_minutes/60 ~=> +0.333~...+0.45000 degrees.
1578 Hour plus fraction when the center of Moon's disk touches the horizon,
1579 corrected for atmospheric refraction.
1580 5 == approximately 5...10arc_minutes/60 ~=> +0.083~...+0.1666~ degrees.
1581 Hour plus fraction when the Moon's upper limb touches the horizon,
1582 corrected for atmospheric refraction (standard rise/set).
1583 6 == This case is managed specially.
1584 The Moon's topocentric apparent horizontal parallax in degrees plus
1585 fraction that appears at midnight on the given date.
1586 7 == This case is managed specially.
1587 The Moon's topocentric apparent semidiameter in degrees plus fraction
1588 that appears at midnight on the given date.
1589 8 == This case is managed specially.
1590 The Moon's topocentric apparent brightness in magnitude units
1591 plus fraction that appears at midnight on the given date.
1592 9 == This case is managed specially.
1593 The topocentric apparent phase angle of the Moon in range
1594 +1.0...0.0 that appears at midnight on the given date.
1595 10 == This case is managed specially.
1596 The topocentric apparent height of the Moon's center under or above
1597 the sea-level horizon (elevation angle) in degrees plus fraction that
1598 appears at midnight on the given date, corrected for atmospheric
1599 refraction.
1600 11 == This case is managed specially.
1601 The Moon's topocentric apparent azimuth angle in degrees plus fraction
1602 that appears at midnight on the given date.
1603 12 == This case is managed specially.
1604 The Moon's topocentric apparent declination angle in degrees
1605 plus fraction that appears at midnight on the given date,
1606 corrected for atmospheric refraction.
1607 13 == This case is managed specially.
1608 The Moon's topocentric apparent ecliptic longitude in degrees
1609 plus fraction that appears at midnight on the given date.
1610 14 == This case is managed specially.
1611 The Moon's topocentric apparent ecliptic latitude in degrees
1612 plus fraction that appears at midnight on the given date.
1613 15 == This case is managed specially.
1614 The Moon's topocentric apparent right ascension angle in hours
1615 plus fraction that appears at midnight on the given date,
1616 corrected for atmospheric refraction.
1617 16 == This case is managed specially.
1618 The topocentric Moon/Earth distance in Earth equator radii plus
1619 fraction that appears at midnight on the given date.
1620 17 == This case is managed specially.
1621 The topocentric Sun/Moon elongation in degrees plus fraction
1622 that appears at midnight on the given date.
1623 18 == This case is managed specially.
1624 The atmospheric refraction which is calculated for correcting
1625 the Moon's topocentric true elevation angle to the Moon's topocentric
1626 apparent elevation angle in degrees plus fraction that appears
1627 at midnight on the given date.
1628 19 == This case is managed specially.
1629 The Moon's geocentric apparent horizontal parallax in degrees plus
1630 fraction that appears at midnight on the given date.
1631 20 == This case is managed specially.
1632 The Moon's geocentric apparent semidiameter in degrees plus fraction
1633 that appears at midnight on the given date.
1634 21 == This case is managed specially.
1635 The Moon's geocentric apparent brightness in magnitude units
1636 plus fraction that appears at midnight on the given date.
1637 22 == This case is managed specially.
1638 The geocentric apparent phase angle of the Moon in range
1639 +1.0...0.0 under or above the imaginary geocentric horizon
1640 that appears at midnight on the given date.
1641 23 == This case is managed specially.
1642 The geocentric apparent height of the Moon's center under or above
1643 the imaginary geocentric horizon (elevation angle) in degrees plus
1644 fraction that appears at midnight on the given date.
1645 24 == This case is managed specially.
1646 The Moon's geocentric apparent azimuth angle in degrees plus fraction
1647 that appears at midnight on the given date.
1648 25 == This case is managed specially.
1649 The Moon's geocentric apparent declination angle in degrees plus
1650 fraction that appears at midnight on the given date.
1651 26 == This case is managed specially.
1652 The Moon's geocentric apparent ecliptic longitude in degrees plus
1653 fraction that appears at midnight on the given date.
1654 27 == This case is managed specially.
1655 The Moon's geocentric apparent ecliptic latitude in degrees plus
1656 fraction that appears at midnight on the given date.
1657 28 == This case is managed specially.
1658 The Moon's geocentric apparent right ascension angle in hours plus
1659 fraction that appears at midnight on the given date.
1660 29 == This case is managed specially.
1661 The geocentric Moon/Earth distance in Earth equator radii plus
1662 fraction that appears at midnight on the given date.
1663 30 == This case is managed specially.
1664 The geocentric Sun/Moon elongation in degrees plus fraction
1665 that appears at midnight on the given date.
1666 31 == This case is managed specially.
1667 The approximate ``Delta-T'' in seconds plus fraction, which is the
1668 difference between Terrestrial Dynamical Time, TDT (formerly called
1669 Ephemeris Time, ET), and Universal Time, UT, so `Delta_T = TDT - UT',
1670 that appears at midnight on the given date.
1671 32 == This case is managed specially.
1672 The local sidereal time in hours plus fraction
1673 that appears at midnight on the given date.
1674 33 == This case is managed specially.
1675 The globally given time offset value in hours plus fraction.
1676 34 == This case is managed specially.
1677 The local ``Julian Date'' in days plus fraction
1678 that appears at midnight on the given date.
1679 35 == This case is managed specially.
1680 The approximate local ``Julian Ephemeris Date'' in days plus fraction
1681 that appears at midnight on the given date.
1682 36 == This case is managed specially.
1683 The difference between the Moon's and the Sun's topocentric
1684 apparent elevation angles in degrees plus fraction that appears
1685 at midnight on the given date.
1686 37 == This case is managed specially.
1687 The difference between the Moon's and the Sun's topocentric
1688 apparent azimuth angles in degrees plus fraction that appears
1689 at midnight on the given date.
1690 38 == This case is managed specially.
1691 The difference between the Moon's and the Sun's geocentric
1692 apparent elevation angles in degrees plus fraction that appears
1693 at midnight on the given date.
1694 39 == This case is managed specially.
1695 The difference between the Moon's and the Sun's geocentric
1696 apparent azimuth angles in degrees plus fraction that appears
1697 at midnight on the given date.
1698 40 == This case is managed specially.
1699 The difference between the Moon's and the Sun's topocentric
1700 apparent elevation angles in degrees plus fraction that appears
1701 at Moon's standard rise/set time.
1702 41 == This case is managed specially.
1703 The difference between the Moon's and the Sun's topocentric
1704 apparent azimuth angles in degrees plus fraction that appears
1705 at Moon's standard rise/set time.
1706 42 == This case is managed specially.
1707 The difference between the Moon's and the Sun's geocentric
1708 apparent elevation angles in degrees plus fraction that appears
1709 at Moon's standard rise/set time.
1710 43 == This case is managed specially.
1711 The difference between the Moon's and the Sun's geocentric
1712 apparent azimuth angles in degrees plus fraction that appears
1713 at Moon's standard rise/set time.
1714 44 == This case is managed specially.
1715 The difference between the Moon's and the Sun's meridian
1716 anti-transit time in hours plus fraction that appears at
1717 Moon's meridian anti-transit time.
1718 45 == This case is managed specially.
1719 The difference between the Moon's and the Sun's meridian
1720 transit time in hours plus fraction that appears at
1721 Moon's meridian transit time.
1722 46 == This case is managed specially.
1723 The difference between the Moon's and the Sun's standard rise/set time
1724 in hours plus fraction that appears at Moon's standard rise/set time.
1725 47 == This case is managed specially.
1726 The Moon's topocentric apparent elevation angle (midnight height)
1727 in degrees plus fraction that appears when the Moon's center
1728 anti-transits the meridian under or above the sea-level horizon,
1729 corrected for atmospheric refraction.
1730 48 == This case is managed specially.
1731 The topocentric apparent phase angle (Moon's midnight phase)
1732 in range +1.0...0.0 that appears when the Moon's center
1733 anti-transits the meridian under or above the sea-level horizon.
1734 49 == This case is managed specially.
1735 The Moon's topocentric apparent elevation angle (noon height)
1736 in degrees plus fraction that appears when the Moon's center
1737 transits the meridian under or above the sea-level horizon,
1738 corrected for atmospheric refraction.
1739 50 == This case is managed specially.
1740 The topocentric apparent phase angle (Moon's noon phase)
1741 in range +1.0...0.0 that appears when the Moon's center
1742 transits the meridian under or above the sea-level horizon.
1743 51 == This case is managed specially.
1744 The Moon's topocentric apparent elevation angle under or above the
1745 true horizon (the given height at the location) in degrees plus
1746 fraction that appears at Moon's standard rise/set time,
1747 corrected for atmospheric refraction.
1748 52 == This case is managed specially.
1749 The Moon's topocentric apparent azimuth angle in degrees plus fraction
1750 that appears at Moon's standard rise/set time.
1751 53 == This case is managed specially.
1752 The Moon's topocentric apparent phase angle in range
1753 +1.0...0.0 that appears at Moon's standard rise/set time.
1754 54 == This case is managed specially.
1755 The Moon's geocentric apparent elevation angle (midnight height)
1756 in degrees plus fraction that appears when the Moon's center
1757 anti-transits the meridian under or above the imaginary geocentric
1758 horizon.
1759 55 == This case is managed specially.
1760 The geocentric apparent phase angle (Moon's midnight phase) in
1761 range +1.0...0.0 that appears when the Moon's center anti-transits
1762 the meridian under or above the imaginary geocentric horizon.
1763 56 == This case is managed specially.
1764 The Moon's geocentric apparent elevation angle (noon height)
1765 in degrees plus fraction that appears when the Moon's center
1766 transits the meridian under or above the imaginary geocentric horizon.
1767 57 == This case is managed specially.
1768 The geocentric apparent phase angle (Moon's noon phase) in
1769 range +1.0...0.0 that appears when the Moon's center transits
1770 the meridian under or above the imaginary geocentric horizon.
1771 58 == This case is managed specially.
1772 The Moon's geocentric apparent elevation angle under or above the
1773 true horizon (the given height at the location) in degrees plus
1774 fraction that appears at Moon's standard rise/set time.
1775 59 == This case is managed specially.
1776 The Moon's geocentric apparent azimuth angle in degrees plus fraction
1777 that appears at Moon's standard rise/set time.
1778 60 == This case is managed specially.
1779 The Moon's geocentric apparent phase angle in range
1780 +1.0...0.0 that appears at Moon's standard rise/set time.
1781 *** CAREFUL ***
1782 Keep in mind that all values computed here are approximations only, based
1783 on low-precision formula terms; and some values are computed recursively!
1784 *** CAREFUL ***
1785 It can happen that this function misses A SECOND rise/set time if the
1786 Moon rises/sets more than once during a single day at some dates at high
1787 latitudes; locations that are within --or close to, if you calculate the
1788 rise/set times for reference altitudes other than the standard ones--
1789 the lunar Arctic Circle.
1790 *** CAREFUL ***
1791 For latitudes that are within the lunar Arctic Circle or nearby
1792 (depending on the reference altitude that is used to calculate the
1793 rise/set times), it occasionally can happen that this function simply
1794 returns SPECIAL_VALUE instead of the --more the astronomically ``truth''
1795 reflecting-- [2|3]*SPECIAL_VALUE value when the Moon does not rise/set
1796 on the given date if there is no merdian transit of the Moon on the
1797 given date!
1798 *** CAREFUL ***
1799 State `meters_above_sea_level' if you know them for the co-ordinate,
1800 otherwise 0. Specify `time_zone_in_mins' in +/- minutes MMM the
1801 co-ordinate is before/behind UT/GMT, so that the calculated UT/GMT-based
1802 rise/set times can be properly converted to local time times. For example,
1803 if a co-ordinates timezone is `GMT-2' (==CEST), use `+120'.
1804 References:
1805 * "Practical Astronomy with your Calculator", Peter Duffet-Smith,
1806 3rd edition. Cambridge University Press 1988, ISBN 0-521-35699-7.
1807 * "Astronomical Formulae for Calculators", Jean Meeus, 4th ed,
1808 Willmann-Bell 1988, ISBN 0-943396-22-0.
1809 - Geocentric apparent ecliptic longitude claimed to 10".
1810 - Geocentric apparent ecliptic latitude claimed to 3".
1811 - Horizontal apparent parallax claimed to 0.2".
1812 - Positions referred to mean equinox of date
1813 but nutation correction for geocentric apparent ecliptic longitude.
1814 * "Astronomical Algorithms", Jean Meeus, 1st ed, Willmann-Bell 1991,
1815 ISBN 0-943396-35-2.
1816 * The WWW documents "How to compute planetary positions" and
1817 "How to compute rise/set times and altitude above horizon"
1818 by Paul Schlyter, <http://welcome.to/pausch> or
1819 <http://hotel04.ausys.se/pausch>.
1820 Btw., Keith Burnett publishes a lot of diverse --good quality--
1821 astronomical sources and links at <http://www.xylem.demon.co.uk/kepler>.
1822 */
1823 {
1824 auto double x1 = 0.0;
1825 auto double x2;
1826 auto double x3 = 0.0;
1827 auto double x4;
1828 auto double longitude;
1829 auto double latitude;
1830 auto double sin_latitude;
1831 auto double cos_latitude;
1832 auto double jc;
1833 auto double mjd1;
1834 auto double mjd1_0ut;
1835 auto double mjd2;
1836 auto double mjd2_0ut;
1837 auto double local_sidereal_time;
1838 auto double obliquity_of_ecliptic;
1839 auto double sin_obliquity_of_ecliptic;
1840 auto double cos_obliquity_of_ecliptic;
1841 auto double argument_of_perihelion;
1842 auto double sma_sun_mean_anomaly;
1843 auto double sma_2;
1844 auto double e_eccentricity;
1845 auto double e_e;
1846 auto double anomaly_of_eccentric;
1847 auto double divisor_term;
1848 auto double true_anomaly;
1849 auto double s_distance;
1850 auto double s_longitude;
1851 auto double m_longitude;
1852 auto double m_mean_longitude;
1853 auto double mma_moon_mean_anomaly;
1854 auto double mma_2;
1855 auto double mma_3;
1856 auto double mme_moon_mean_elongation;
1857 auto double mme_2;
1858 auto double mme_4;
1859 auto double m_latitude;
1860 auto double m_2_latitude;
1861 auto double m_3_latitude;
1862 auto double m_parallax;
1863 auto double m_semidiameter;
1864 auto double m_distance;
1865 auto double xg;
1866 auto double yg;
1867 auto double zg;
1868 auto double xe;
1869 auto double ye;
1870 auto double ze;
1871 auto double right_ascension;
1872 auto double declination;
1873 auto double meridian_transit_time = SECS_PER_DAY;
1874 auto double the_time = SECS_PER_DAY;
1875 auto double tt = SECS_PER_DAY;
1876 auto double geodetic_height = 0.0;
1877 register int i = 0;
1878 register int j = 0;
1879 register int n = 0;
1880 register int tz_offset;
1881 register int save_the_mode;
1882 auto Bool ok = FALSE;
1883 auto Bool is_adjusted = FALSE;
1884 auto Bool iter_mt = FALSE;
1885
1886
1887 /*
1888 Reduce the given timezone value in minutes to a single day.
1889 */
1890 tz_offset = (int) (SGN (coordinates->time_zone_in_mins)
1891 * (abs (coordinates->time_zone_in_mins) % MINS_PER_DAY));
1892 x2 = DAY2HH (time_offset);
1893 x4 = ((int) DAY2MM (time_offset)) % MINS_PER_HOUR;
1894 switch (coordinates->the_mode)
1895 {
1896 case 0:
1897 /*
1898 Calculate the meridian anti-transit time (Moon's astronomical midnight).
1899 */
1900 coordinates->the_mode = 1;
1901 x1 = internal_moon_rise_set (event, day, month, year, coordinates);
1902 (void) next_date (&day, &month, &year);
1903 x2 = internal_moon_rise_set (event, day, month, year, coordinates);
1904 the_time = x1 + ((x2 - x1) + HOURS_PER_DAY) * 0.5;
1905 (void) prev_date (&day, &month, &year);
1906 if (the_time >= HOURS_PER_DAY
1907 || x1 == SPECIAL_VALUE || x2 == SPECIAL_VALUE)
1908 {
1909 if (x1 == SPECIAL_VALUE)
1910 x1 = x2 + HOURS_PER_DAY;
1911 (void) prev_date (&day, &month, &year);
1912 x2 = internal_moon_rise_set (event, day, month, year, coordinates);
1913 the_time = x2 + ((x1 - x2) + HOURS_PER_DAY) * 0.5;
1914 if (the_time < HOURS_PER_DAY)
1915 the_time = SPECIAL_VALUE;
1916 else
1917 {
1918 the_time -= HOURS_PER_DAY;
1919 (void) next_date (&day, &month, &year);
1920 }
1921 }
1922 if (the_time != SPECIAL_VALUE)
1923 {
1924 x1 = time_offset;
1925 coordinates->the_mode = 24;
1926 n = 0;
1927 e_e = MM2DAY (tz_offset);
1928 /*
1929 Start binary-search with margins set to +/- 2 minutes
1930 around the estimated value previously calculated.
1931 */
1932 x3 = HH2DAY (the_time);
1933 x4 = HHMM2DAY (0, 2);
1934 x2 = x3 - x4;
1935 x3 += x4;
1936 x4 = (x3 + x2) * 0.5;
1937 LOOP
1938 {
1939 time_offset = x4 + e_e;
1940 the_time =
1941 internal_moon_rise_set (event, day, month, year, coordinates);
1942 the_time = FIXANGLE (the_time + DEGS_PER_06_HOURS);
1943 if (the_time > DEGS_PER_12_HOURS)
1944 the_time = -(the_time - DEGS_PER_24_HOURS);
1945 /*
1946 Exit interpolation loop if the actually computed value is less
1947 than 0.0001 degrees smaller the previously computed value.
1948 A ``next'' interpolation step would not change the result
1949 significantly -- the only change would be by the order of
1950 small fractions of seconds.
1951 */
1952 if (abs (DEGS_PER_06_HOURS - the_time) > 0.0001)
1953 {
1954 /*
1955 Yes, perform next interpolation step with adjusted margins.
1956 */
1957 if (the_time > DEGS_PER_06_HOURS)
1958 x3 = x4;
1959 else
1960 x2 = x4;
1961 x4 = (x3 + x2) * 0.5;
1962 }
1963 else
1964 /*
1965 No, the interpolation result is precise enough now.
1966 */
1967 break;
1968 /*
1969 Avoid neverending loop if interpolation divergence is detected.
1970 */
1971 if (++n >= 100)
1972 my_error (ERR_INTERNAL_C_FUNC_FAILURE, __FILE__,
1973 ((long) __LINE__) - 1L,
1974 "internal_moon_rise_set():interpolation<", n);
1975 }
1976 the_time = DAY2HH (x4);
1977 time_offset = x1;
1978 }
1979 coordinates->the_mode = 0;
1980 return (the_time);
1981 case 31:
1982 /*
1983 Return ``Delta-T'' in seconds plus fraction.
1984 */
1985 return (delta_t (day, month, year, (int) x2, (int) x4));
1986 case 33:
1987 /*
1988 Return the globally given time offset value in hours plus fraction.
1989 */
1990 return (DAY2HH (time_offset));
1991 case 35:
1992 /*
1993 Return the local ``Julian Ephemeris Date'' in days plus fraction.
1994 */
1995 x1 = SS2DAY (delta_t (day, month, year, (int) x2, (int) x4));
1996 /* Fallthrough. */
1997 case 34:
1998 /*
1999 Return the local ``Julian Date'' in days plus fraction.
2000 */
2001 return (date2num (day, month, year) + (MIN_BCE_TO_1_CE - 1.5) + x1 +
2002 time_offset);
2003 case 36:
2004 case 37:
2005 case 38:
2006 case 39:
2007 save_the_mode = coordinates->the_mode;
2008 switch (save_the_mode)
2009 {
2010 case 36:
2011 /*
2012 Calculate the Moon's topocentric apparent elevation angle
2013 under or above the sea-level horizon in degrees plus
2014 fraction, corrected for atmospheric refraction.
2015 */
2016 coordinates->the_mode = 10;
2017 break;
2018 case 37:
2019 /*
2020 Calculate the Moon's topocentric apparent azimuth angle
2021 in degrees plus fraction.
2022 */
2023 coordinates->the_mode = 11;
2024 break;
2025 case 38:
2026 /*
2027 Calculate the Moon's geocentric apparent elevation angle
2028 under or above the imaginary geocentric horizon
2029 in degrees plus fraction.
2030 */
2031 coordinates->the_mode = 23;
2032 break;
2033 default:
2034 /*
2035 Calculate the Moon's geocentric apparent azimuth angle
2036 in degrees plus fraction.
2037 */
2038 coordinates->the_mode = 24;
2039 break;
2040 }
2041 the_time =
2042 internal_moon_rise_set (event, day, month, year, coordinates);
2043 switch (save_the_mode)
2044 {
2045 case 36:
2046 case 38:
2047 if (save_the_mode == 36)
2048 /*
2049 Calculate the difference between the Moon's and the Sun's
2050 topocentric apparent elevation angles in degrees plus
2051 fraction that appears at midnight on the given date.
2052 */
2053 coordinates->the_mode = 10;
2054 else
2055 /*
2056 Calculate the difference between the Moon's and the Sun's
2057 geocentric apparent elevation angles in degrees plus
2058 fraction that appears at midnight on the given date.
2059 */
2060 coordinates->the_mode = 19;
2061 the_time -=
2062 sun_rise_set (event, FALSE, day, month, year, coordinates);
2063 break;
2064 default:
2065 if (save_the_mode == 37)
2066 /*
2067 Calculate the difference between the Moon's and the Sun's
2068 topocentric apparent azimuth angles in degrees plus
2069 fraction that appears at midnight on the given date.
2070 */
2071 coordinates->the_mode = 11;
2072 else
2073 /*
2074 Calculate the difference between the Moon's and the Sun's
2075 geocentric apparent azimuth angles in degrees plus
2076 fraction that appears at midnight on the given date.
2077 */
2078 coordinates->the_mode = 20;
2079 x4 =
2080 the_time - sun_rise_set (event, FALSE, day, month, year,
2081 coordinates);
2082 /*
2083 Reduce the resulting angle to -180.0...+180.0 degrees.
2084 */
2085 if (abs (x4) <= DEGS_PER_12_HOURS)
2086 the_time = x4;
2087 else if (x4 > DEGS_PER_12_HOURS)
2088 the_time = -DEGS_PER_24_HOURS + x4;
2089 else
2090 the_time = DEGS_PER_24_HOURS + x4;
2091 }
2092 coordinates->the_mode = save_the_mode;
2093 return (the_time);
2094 case 44:
2095 case 45:
2096 case 46:
2097 case 47:
2098 case 48:
2099 case 49:
2100 case 50:
2101 case 54:
2102 case 55:
2103 case 56:
2104 case 57:
2105 x1 = time_offset;
2106 save_the_mode = coordinates->the_mode;
2107 switch (save_the_mode)
2108 {
2109 case 46:
2110 /*
2111 Perform calculations for Moon's standard rise/set time.
2112 */
2113 coordinates->the_mode = 5;
2114 break;
2115 case 44:
2116 case 47:
2117 case 48:
2118 case 54:
2119 case 55:
2120 /*
2121 Perform calculations for Moon's astronomical midnight.
2122 */
2123 coordinates->the_mode = 0;
2124 break;
2125 default:
2126 /*
2127 Perform calculations for Moon's astronomical noon.
2128 */
2129 coordinates->the_mode = 1;
2130 }
2131 time_offset = moon_rise_set (event, day, month, year, coordinates);
2132 if (time_offset > SPECIAL_VALUE)
2133 {
2134 time_offset = HH2DAY (time_offset);
2135 if (save_the_mode <= 45)
2136 {
2137 if (save_the_mode == 44)
2138 /*
2139 Calculate the difference between the Moon's and the Sun's
2140 meridian anti-transit time in hours plus fraction that
2141 appears at Moon's meridian anti-transit time.
2142 */
2143 coordinates->the_mode = 0;
2144 else
2145 /*
2146 Calculate the difference between the Moon's and the Sun's
2147 meridian transit time in hours plus fraction that appears
2148 at Moon's meridian transit time.
2149 */
2150 coordinates->the_mode = 1;
2151 the_time =
2152 sun_rise_set (event, FALSE, day, month, year, coordinates);
2153 if (the_time > SPECIAL_VALUE)
2154 the_time = DAY2HH (time_offset) - the_time;
2155 else
2156 the_time = HH2SS (the_time);
2157 }
2158 else if (save_the_mode == 46)
2159 {
2160 /*
2161 Calculate the difference between the Moon's and the Sun's
2162 standard rise/set time in hours plus fraction that appears
2163 at Moon's standard rise/set time.
2164 */
2165 coordinates->the_mode = 5;
2166 the_time =
2167 sun_rise_set (event, FALSE, day, month, year, coordinates);
2168 if (the_time > SPECIAL_VALUE)
2169 the_time = DAY2HH (time_offset) - the_time;
2170 else
2171 the_time = HH2SS (the_time);
2172 }
2173 else
2174 {
2175 switch (save_the_mode)
2176 {
2177 case 47:
2178 case 49:
2179 /*
2180 Calculate the Moon's topocentric apparent elevation angle
2181 under or above the sea-level horizon in degrees plus
2182 fraction, corrected for atmospheric refraction.
2183 */
2184 coordinates->the_mode = 10;
2185 break;
2186 case 48:
2187 case 50:
2188 /*
2189 Calculate the topocentric apparent phase angle in range
2190 +1.0...0.0 under or above the sea-level horizon.
2191 */
2192 coordinates->the_mode = 9;
2193 break;
2194 case 54:
2195 case 56:
2196 /*
2197 Calculate the Moon's geocentric apparent elevation angle
2198 under or above the sea-level horizon in degrees plus
2199 fraction.
2200 */
2201 coordinates->the_mode = 23;
2202 break;
2203 default:
2204 /*
2205 Calculate the geocentric apparent phase angle in range
2206 +1.0...0.0 under or above the imaginary geocentric
2207 horizon.
2208 */
2209 coordinates->the_mode = 22;
2210 }
2211 the_time = moon_rise_set (event, day, month, year, coordinates);
2212 }
2213 }
2214 else
2215 the_time = HH2SS (time_offset);
2216 coordinates->the_mode = save_the_mode;
2217 time_offset = x1;
2218 return (the_time);
2219 case 40:
2220 case 41:
2221 case 42:
2222 case 43:
2223 case 51:
2224 case 52:
2225 case 53:
2226 case 58:
2227 case 59:
2228 case 60:
2229 x1 = time_offset;
2230 save_the_mode = coordinates->the_mode;
2231 coordinates->the_mode = 5;
2232 time_offset = moon_rise_set (event, day, month, year, coordinates);
2233 if (time_offset > SPECIAL_VALUE)
2234 {
2235 time_offset = HH2DAY (time_offset);
2236 switch (save_the_mode)
2237 {
2238 case 40:
2239 case 51:
2240 /*
2241 Calculate the Moon's topocentric apparent elevation angle
2242 under or above the true horizon in degrees plus fraction
2243 that appears at Moon's standard rise/set time,
2244 corrected for atmospheric refraction.
2245 */
2246 coordinates->the_mode = 10;
2247 break;
2248 case 41:
2249 case 52:
2250 /*
2251 Calculate the Moon's topocentric apparent azimuth angle
2252 in degrees plus fraction that appears at Moon's
2253 standard rise/set time.
2254 */
2255 coordinates->the_mode = 11;
2256 break;
2257 case 53:
2258 /*
2259 Calculate the topocentric apparent phase angle
2260 of the Moon in range +1.0...0.0 that appears
2261 at Moon's standard rise/set time.
2262 */
2263 coordinates->the_mode = 9;
2264 break;
2265 case 42:
2266 case 58:
2267 /*
2268 Calculate the Moon's geocentric apparent elevation angle
2269 under or above the true horizon in degrees plus fraction
2270 that appears at Moon's standard rise/set time.
2271 */
2272 coordinates->the_mode = 23;
2273 break;
2274 case 43:
2275 case 59:
2276 /*
2277 Calculate the Moon's geocentric apparent azimuth angle
2278 in degrees plus fraction that appears at Moon's
2279 standard rise/set time.
2280 */
2281 coordinates->the_mode = 24;
2282 break;
2283 default:
2284 /*
2285 Calculate the geocentric apparent phase angle
2286 of the Moon in range +1.0...0.0 that appears
2287 at Moon's standard rise/set time.
2288 */
2289 coordinates->the_mode = 22;
2290 }
2291 the_time = moon_rise_set (event, day, month, year, coordinates);
2292 if ((save_the_mode >= 40) && (save_the_mode <= 43))
2293 switch (save_the_mode)
2294 {
2295 case 40:
2296 case 42:
2297 if (save_the_mode == 40)
2298 /*
2299 Calculate the difference between the Moon's and the Sun's
2300 topocentric apparent elevation angles in degrees plus
2301 fraction that appears at Moon's standard rise/set time.
2302 */
2303 coordinates->the_mode = 10;
2304 else
2305 /*
2306 Calculate the difference between the Moon's and the Sun's
2307 geocentric apparent elevation angles in degrees plus
2308 fraction that appears at Moon's standard rise/set time.
2309 */
2310 coordinates->the_mode = 19;
2311 the_time -=
2312 sun_rise_set (event, FALSE, day, month, year, coordinates);
2313 break;
2314 default:
2315 if (save_the_mode == 41)
2316 /*
2317 Calculate the difference between the Moon's and the Sun's
2318 topocentric apparent azimuth angles in degrees plus
2319 fraction that appears at Moon's standard rise/set time.
2320 */
2321 coordinates->the_mode = 11;
2322 else
2323 /*
2324 Calculate the difference between the Moon's and the Sun's
2325 geocentric apparent azimuth angles in degrees plus
2326 fraction that appears at Moon's standard rise/set time.
2327 */
2328 coordinates->the_mode = 20;
2329 x4 =
2330 the_time - sun_rise_set (event, FALSE, day, month, year,
2331 coordinates);
2332 /*
2333 Reduce the resulting angle to -180.0...+180.0 degrees.
2334 */
2335 if (abs (x4) <= DEGS_PER_12_HOURS)
2336 the_time = x4;
2337 else if (x4 > DEGS_PER_12_HOURS)
2338 the_time = -DEGS_PER_24_HOURS + x4;
2339 else
2340 the_time = DEGS_PER_24_HOURS + x4;
2341 }
2342 }
2343 else
2344 the_time = HH2SS (time_offset);
2345 coordinates->the_mode = save_the_mode;
2346 time_offset = x1;
2347 return (the_time);
2348 case 1:
2349 iter_mt = TRUE;
2350 /* Fallthrough. */
2351 default:
2352 longitude = coordinates->lon_deg
2353 + MM2DEG (coordinates->lon_min) + SS2DEG (coordinates->lon_sec);
2354 }
2355 latitude = TORAD (coordinates->lat_deg
2356 + MM2DEG (coordinates->lat_min)
2357 + SS2DEG (coordinates->lat_sec));
2358 if (coordinates->meters_above_sea_level > 0)
2359 geodetic_height =
2360 TORAD (0.0347 * sqrt (coordinates->meters_above_sea_level));
2361 /*
2362 Convert UT date to TDT date.
2363 */
2364 e_e = SS2DAY (delta_t (day, month, year, 0, 0));
2365 x1 = date2num (day, month, year) + e_e;
2366 /*
2367 Epoch-1 starts 01-Jan-2000 12:00:00 UT.
2368 */
2369 mjd1 = mjd1_0ut = x1 - 730122.5;
2370 /*
2371 Epoch-2 starts 31-Dec-1899 12:00:00 UT.
2372 */
2373 mjd2 = mjd2_0ut = x1 - 693597.5;
2374 /*
2375 Respect timezone and time displacement value
2376 for Moon's non-rise/set/transit based modes.
2377 */
2378 if ((coordinates->the_mode >= 6) && (coordinates->the_mode <= 39))
2379 {
2380 x1 = time_offset - MM2DAY (tz_offset)
2381 + SS2DAY (delta_t (day, month, year, (int) x2, (int) x4)) - e_e;
2382 mjd1 += x1;
2383 mjd2 += x1;
2384 }
2385 while (!ok)
2386 {
2387 do
2388 {
2389 LABEL_moon_iter_mt:
2390 local_sidereal_time = FIXANGLE (100.46061837
2391 + 0.9856472275 * mjd1 + longitude);
2392 /*
2393 Compute Sun data.
2394 */
2395 jc = mjd1 / 36525.0;
2396 obliquity_of_ecliptic = TORAD (FIXANGLE (23.43929111
2397 - (0.013004167 +
2398 (0.00000016389 -
2399 0.0000005036 * jc) *
2400 jc) * jc));
2401 argument_of_perihelion =
2402 TORAD (FIXANGLE (282.93735 + (0.71953 + (0.00046 * jc)) * jc));
2403 sma_sun_mean_anomaly =
2404 TORAD (FIXANGLE
2405 (357.52910 +
2406 (35999.05030 - (0.0001559 + 0.00000048 * jc) * jc) * jc));
2407 e_eccentricity =
2408 0.016708617 - (0.000042037 + (0.0000001236 * jc)) * jc;
2409 x2 =
2410 sma_sun_mean_anomaly +
2411 e_eccentricity * sin (sma_sun_mean_anomaly) * (1.0 +
2412 e_eccentricity *
2413 cos
2414 (sma_sun_mean_anomaly));
2415 do
2416 {
2417 x1 = x2;
2418 x2 =
2419 x1 - (x1 - e_eccentricity * sin (x1) -
2420 sma_sun_mean_anomaly) / (1.0 -
2421 e_eccentricity * cos (x1));
2422 }
2423 while (abs (abs (x1) - abs (x2)) > 0.000001);
2424 anomaly_of_eccentric = x2;
2425 divisor_term = cos (anomaly_of_eccentric) - e_eccentricity;
2426 s_distance =
2427 sqrt (1.0 -
2428 e_eccentricity * e_eccentricity) *
2429 sin (anomaly_of_eccentric);
2430 true_anomaly = atan (s_distance / divisor_term);
2431 if (divisor_term < 0.0)
2432 true_anomaly += MY_PI;
2433 s_longitude =
2434 TORAD (FIXANGLE
2435 (TODEG (true_anomaly + argument_of_perihelion) + jc -
2436 0.00569 -
2437 0.00479 *
2438 sin (TORAD
2439 (124.90 - (1934.134 - (0.002063 * jc)) * jc))));
2440 /*
2441 Compute Moon data.
2442 */
2443 jc = mjd2 / 36525.0;
2444 /*
2445 Precompute some constant values first.
2446 */
2447 x1 = sin (TORAD (FIXANGLE (51.2 + 20.2 * jc)));
2448 x2 = 0.0039640 * sin (TORAD (FIXANGLE (346.560
2449 + (132.870 -
2450 (0.0091731 * jc)) * jc)));
2451 m_longitude =
2452 TORAD (FIXANGLE
2453 (259.183275 -
2454 (1934.142 - (0.002078 + 0.0000022 * jc) * jc) * jc));
2455 /*
2456 Precompute more constant values.
2457 */
2458 x3 = sin (m_longitude);
2459 x4 = TORAD (275.05 - 2.30 * jc);
2460 /*
2461 Compute orbit data.
2462 */
2463 m_mean_longitude = FIXANGLE (270.434164
2464 + (481267.8831 -
2465 (0.001133 -
2466 0.0000019 * jc) * jc) * jc +
2467 0.000233 * x1 + x2 + 0.001964 * x3);
2468 # if 0
2469 /*
2470 That's Meeus' original term, but here we use the
2471 above computed one for speed-up purposes because
2472 the results do not differ very much.
2473 */
2474 sma_sun_mean_anomaly = TORAD (FIXANGLE (358.475833
2475 + (35999.0498 -
2476 (0.00015 +
2477 0.0000033 * jc) * jc) *
2478 jc - 0.001778 * x1));
2479 # endif
2480 mma_moon_mean_anomaly = TORAD (FIXANGLE (296.104608
2481 + (477198.8491 +
2482 (0.009192 +
2483 0.0000144 * jc) * jc) *
2484 jc + 0.000817 * x1 + x2 +
2485 0.002541 * x3));
2486 mme_moon_mean_elongation =
2487 TORAD (FIXANGLE
2488 (350.737486 +
2489 (445267.1142 - (0.001436 - 0.0000019 * jc) * jc) * jc +
2490 0.002011 * x1 + x2 + 0.001964 * x3));
2491 m_latitude =
2492 TORAD (FIXANGLE
2493 (11.250889 +
2494 (483202.0251 - (0.003211 + 0.0000003 * jc) * jc) * jc +
2495 x2 - 0.024691 * x3 - 0.004328 * sin (m_longitude + x4)));
2496 e_eccentricity = 1.0 - (0.002495 + (0.00000752 * jc)) * jc;
2497 /*
2498 Precompute more constant values.
2499 */
2500 mme_2 = 2.0 * mme_moon_mean_elongation;
2501 mme_4 = 4.0 * mme_moon_mean_elongation;
2502 sma_2 = 2.0 * sma_sun_mean_anomaly;
2503 mma_2 = 2.0 * mma_moon_mean_anomaly;
2504 mma_3 = 3.0 * mma_moon_mean_anomaly;
2505 m_2_latitude = 2.0 * m_latitude;
2506 m_3_latitude = 3.0 * m_latitude;
2507 e_e = e_eccentricity * e_eccentricity;
2508 zg = m_latitude + mme_4;
2509 yg = mme_4 - m_latitude;
2510 xg = mme_2 - m_latitude;
2511 ze = mme_2 + m_latitude;
2512 ye = mme_4 - mma_moon_mean_anomaly;
2513 xe = mme_4 - sma_sun_mean_anomaly;
2514 x1 = mme_2 - sma_sun_mean_anomaly;
2515 x2 = mme_2 - mma_moon_mean_anomaly;
2516 /*
2517 Calculate the Moon's geocentric apparent horizontal parallax.
2518 */
2519 m_parallax = 0.950724
2520 + 0.051818 * cos (mma_moon_mean_anomaly)
2521 + 0.009531 * cos (x2)
2522 + 0.007843 * cos (mme_2)
2523 + 0.002824 * cos (mma_2)
2524 + 0.000857 * cos (mme_2 + mma_moon_mean_anomaly)
2525 + 0.000533 * cos (x1) * e_eccentricity
2526 + 0.000401 * cos (x1 - mma_moon_mean_anomaly) * e_eccentricity
2527 + 0.000320 * cos (mma_moon_mean_anomaly -
2528 sma_sun_mean_anomaly) * e_eccentricity -
2529 0.000271 * cos (mme_moon_mean_elongation) -
2530 0.000264 * cos (sma_sun_mean_anomaly +
2531 mma_moon_mean_anomaly) * e_eccentricity -
2532 0.000198 * cos (m_2_latitude - mma_moon_mean_anomaly) +
2533 0.000173 * cos (mma_3) + 0.000167 * cos (ye) -
2534 0.000111 * cos (sma_sun_mean_anomaly) * e_eccentricity +
2535 0.000103 * cos (mme_4 - mma_2) - 0.000084 * cos (mma_2 - mme_2) -
2536 0.000083 * cos (mme_2 + sma_sun_mean_anomaly) * e_eccentricity +
2537 0.000079 * cos (mme_2 + mma_2) + 0.000072 * cos (mme_4) +
2538 0.000064 * cos (x1 + mma_moon_mean_anomaly) * e_eccentricity -
2539 0.000063 * cos (mme_2 + sma_sun_mean_anomaly -
2540 mma_moon_mean_anomaly) * e_eccentricity +
2541 0.000041 * cos (sma_sun_mean_anomaly +
2542 mme_moon_mean_elongation) * e_eccentricity +
2543 0.000035 * cos (mma_2 - sma_sun_mean_anomaly) * e_eccentricity -
2544 0.000033 * cos (mma_3 - mme_2) -
2545 0.000030 * cos (mma_moon_mean_anomaly +
2546 mme_moon_mean_elongation) -
2547 0.000029 * cos (m_2_latitude - mme_2) - 0.000029 * cos (mma_2 +
2548 sma_sun_mean_anomaly)
2549 * e_eccentricity + 0.000026 * cos (mme_2 - sma_2) * e_e -
2550 0.000023 * cos (m_2_latitude - mme_2 + mma_moon_mean_anomaly) +
2551 0.000019 * cos (xe - mma_moon_mean_anomaly) * e_eccentricity;
2552 if (coordinates->the_mode == 19)
2553 /*
2554 Return the Moon's geocentric apparent horizontal parallax
2555 in degrees plus fraction.
2556 */
2557 return (m_parallax);
2558 m_parallax = TORAD (m_parallax);
2559 /*
2560 Calculate the geocentric Moon/Earth distance in Earth equator
2561 radii plus fraction.
2562 */
2563 m_distance = 1.0 / sin (m_parallax);
2564 if (coordinates->the_mode == 29)
2565 return (m_distance);
2566 /*
2567 Calculate the Moon's geocentric apparent semidiameter
2568 in degrees plus fraction.
2569 */
2570 m_semidiameter = SS2DEG (936.85 * 60.0 / m_distance);
2571 if (coordinates->the_mode == 20)
2572 return (m_semidiameter);
2573 /*
2574 Calculate the Moon's geocentric apparent ecliptic latitude.
2575 */
2576 m_latitude = 5.128189 * sin (m_latitude)
2577 + 0.280606 * sin (mma_moon_mean_anomaly + m_latitude)
2578 + 0.277693 * sin (mma_moon_mean_anomaly - m_latitude)
2579 + 0.173238 * sin (xg)
2580 + 0.055413 * sin (ze - mma_moon_mean_anomaly)
2581 + 0.046272 * sin (xg - mma_moon_mean_anomaly)
2582 + 0.032573 * sin (ze)
2583 + 0.017198 * sin (mma_2 + m_latitude)
2584 + 0.009267 * sin (mme_2 + mma_moon_mean_anomaly - m_latitude)
2585 + 0.008823 * sin (mma_2 - m_latitude)
2586 + 0.008247 * sin (x1 - m_latitude) * e_eccentricity
2587 + 0.004323 * sin (xg - mma_2)
2588 + 0.004200 * sin (ze + mma_moon_mean_anomaly)
2589 + 0.003372 * sin (m_latitude - sma_sun_mean_anomaly -
2590 mme_2) * e_eccentricity + 0.002472 * sin (ze -
2591 sma_sun_mean_anomaly
2592 -
2593 mma_moon_mean_anomaly)
2594 * e_eccentricity + 0.002222 * sin (ze -
2595 sma_sun_mean_anomaly) *
2596 e_eccentricity + 0.002072 * sin (xg - sma_sun_mean_anomaly -
2597 mma_moon_mean_anomaly) *
2598 e_eccentricity + 0.001877 * sin (m_latitude -
2599 sma_sun_mean_anomaly +
2600 mma_moon_mean_anomaly) *
2601 e_eccentricity + 0.001828 * sin (yg - mma_moon_mean_anomaly) -
2602 0.001803 * sin (m_latitude +
2603 sma_sun_mean_anomaly) * e_eccentricity -
2604 0.001750 * sin (m_3_latitude) +
2605 0.001570 * sin (mma_moon_mean_anomaly - sma_sun_mean_anomaly -
2606 m_latitude) * e_eccentricity -
2607 0.001487 * sin (m_latitude + mme_moon_mean_elongation) -
2608 0.001481 * sin (m_latitude + sma_sun_mean_anomaly +
2609 mma_moon_mean_anomaly) * e_eccentricity +
2610 0.001417 * sin (m_latitude - sma_sun_mean_anomaly -
2611 mma_moon_mean_anomaly) * e_eccentricity +
2612 0.001350 * sin (m_latitude -
2613 sma_sun_mean_anomaly) * e_eccentricity +
2614 0.001330 * sin (m_latitude - mme_moon_mean_elongation) +
2615 0.001106 * sin (m_latitude + mma_3) + 0.001020 * sin (yg) +
2616 0.000833 * sin (zg - mma_moon_mean_anomaly) +
2617 0.000781 * sin (mma_moon_mean_anomaly - m_3_latitude) +
2618 0.000670 * sin (zg - mma_2) + 0.000606 * sin (mme_2 -
2619 m_3_latitude) +
2620 0.000597 * sin (mme_2 + mma_2 - m_latitude) +
2621 0.000492 * sin (mme_2 + mma_moon_mean_anomaly -
2622 sma_sun_mean_anomaly -
2623 m_latitude) * e_eccentricity +
2624 0.000450 * sin (mma_2 - m_latitude - mme_2) +
2625 0.000439 * sin (mma_3 - m_latitude) + 0.000423 * sin (m_latitude +
2626 mme_2 +
2627 mma_2) +
2628 0.000422 * sin (xg - mma_3) -
2629 0.000367 * sin (sma_sun_mean_anomaly + m_latitude + mme_2 -
2630 mma_moon_mean_anomaly) * e_eccentricity -
2631 0.000353 * sin (sma_sun_mean_anomaly + m_latitude +
2632 mme_2) * e_eccentricity + 0.000331 * sin (zg) +
2633 0.000317 * sin (ze - sma_sun_mean_anomaly +
2634 mma_moon_mean_anomaly) * e_eccentricity +
2635 0.000306 * sin (mme_2 - sma_2 - m_latitude) * e_e -
2636 0.000283 * sin (mma_moon_mean_anomaly + m_3_latitude)
2637 /* Nutation term. */
2638 * (1.0 - (0.0004664 * cos (m_longitude)) -
2639 (0.0000754 * cos (m_longitude + x4)));
2640 m_latitude = TORAD (m_latitude);
2641 /*
2642 Correction term for far away dates past/prior the epoch.
2643 */
2644 if (m_latitude < 0.0)
2645 m_latitude -= (jc * 0.0000111);
2646 else
2647 m_latitude += (jc * 0.0000111);
2648 if (coordinates->the_mode == 27)
2649 /*
2650 Return the Moon's geocentric apparent ecliptic latitude
2651 in degrees plus fraction.
2652 */
2653 return (TODEG (m_latitude));
2654 /*
2655 Calculate the Moon's geocentric apparent ecliptic longitude.
2656 */
2657 m_longitude = m_mean_longitude
2658 + 6.288750 * sin (mma_moon_mean_anomaly)
2659 + 1.274018 * sin (x2)
2660 + 0.658309 * sin (mme_2)
2661 + 0.213616 * sin (mma_2)
2662 - 0.185596 * sin (sma_sun_mean_anomaly) * e_eccentricity
2663 - 0.114336 * sin (m_2_latitude)
2664 + 0.058793 * sin (mme_2 - mma_2)
2665 + 0.057212 * sin (x1 - mma_moon_mean_anomaly) * e_eccentricity
2666 + 0.053320 * sin (mme_2 + mma_moon_mean_anomaly)
2667 + 0.045874 * sin (x1) * e_eccentricity
2668 + 0.041024 * sin (mma_moon_mean_anomaly -
2669 sma_sun_mean_anomaly) * e_eccentricity -
2670 0.034718 * sin (mme_moon_mean_elongation) -
2671 0.030465 * sin (sma_sun_mean_anomaly +
2672 mma_moon_mean_anomaly) * e_eccentricity +
2673 0.015326 * sin (mme_2 - m_2_latitude) -
2674 0.012528 * sin (m_2_latitude + mma_moon_mean_anomaly) -
2675 0.010980 * sin (m_2_latitude - mma_moon_mean_anomaly) +
2676 0.010674 * sin (ye) + 0.010034 * sin (mma_3) +
2677 0.008548 * sin (mme_4 - mma_2) -
2678 0.007910 * sin (sma_sun_mean_anomaly - mma_moon_mean_anomaly +
2679 mme_2) * e_eccentricity - 0.006783 * sin (mme_2 +
2680 sma_sun_mean_anomaly)
2681 * e_eccentricity + 0.005162 * sin (mma_moon_mean_anomaly -
2682 mme_moon_mean_elongation) +
2683 0.005000 * sin (sma_sun_mean_anomaly +
2684 mme_moon_mean_elongation) * e_eccentricity +
2685 0.004049 * sin (mma_moon_mean_anomaly - sma_sun_mean_anomaly +
2686 mme_2) * e_eccentricity + 0.003996 * sin (mma_2 +
2687 mme_2) +
2688 0.003862 * sin (mme_4) + 0.003665 * sin (mme_2 - mma_3) +
2689 0.002695 * sin (mma_2 - sma_sun_mean_anomaly) * e_eccentricity +
2690 0.002602 * sin (mma_moon_mean_anomaly - m_2_latitude - mme_2) +
2691 0.002396 * sin (x1 - mma_2) * e_eccentricity -
2692 0.002349 * sin (mma_moon_mean_anomaly +
2693 mme_moon_mean_elongation) +
2694 0.002249 * sin (mme_2 - sma_2) * e_e - 0.002125 * sin (mma_2 +
2695 sma_sun_mean_anomaly)
2696 * e_eccentricity - 0.002079 * sin (sma_2) * e_e +
2697 0.002059 * sin (x2 - sma_2) * e_e -
2698 0.001773 * sin (mma_moon_mean_anomaly + mme_2 - m_2_latitude) -
2699 0.001595 * sin (m_2_latitude + mme_2) + 0.001220 * sin (xe -
2700 mma_moon_mean_anomaly)
2701 * e_eccentricity - 0.001110 * sin (mma_2 + m_2_latitude) +
2702 0.000892 * sin (mma_moon_mean_anomaly -
2703 3.0 * mme_moon_mean_elongation) -
2704 0.000811 * sin (sma_sun_mean_anomaly + mma_moon_mean_anomaly +
2705 mme_2) * e_eccentricity + 0.000761 * sin (xe -
2706 mma_2) *
2707 e_eccentricity + 0.000717 * sin (mma_moon_mean_anomaly -
2708 sma_2) * e_e +
2709 0.000704 * sin (mma_moon_mean_anomaly - sma_2 - mme_2) * e_e +
2710 0.000693 * sin (sma_sun_mean_anomaly - mma_2 +
2711 mme_2) * e_eccentricity + 0.000598 * sin (x1 -
2712 m_2_latitude)
2713 * e_eccentricity + 0.000550 * sin (mma_moon_mean_anomaly +
2714 mme_4) +
2715 0.000538 * sin (4.0 * mma_moon_mean_anomaly) +
2716 0.000521 * sin (xe) * e_eccentricity + 0.000486 * sin (mma_2 -
2717 mme_moon_mean_elongation)
2718 /* Nutation term. */
2719 + SS2DEG (-17.2 * x3 -
2720 1.3 * sin (2.0 *
2721 TORAD (FIXANGLE
2722 (279.6967 +
2723 (36000.7689 +
2724 (0.000303 * jc)) * jc))))
2725 /* Correction term for far away dates past/prior the epoch. */
2726 - 0.0005915 * jc * jc;
2727 m_longitude = FIXANGLE (m_longitude);
2728 if (coordinates->the_mode == 26)
2729 /*
2730 Return the Moon's geocentric apparent ecliptic longitude
2731 in degrees plus fraction.
2732 */
2733 return (m_longitude);
2734 m_longitude = TORAD (m_longitude);
2735 /*
2736 Calculate the geocentric Sun/Moon elongation.
2737 */
2738 x4 = acos (cos (s_longitude - m_longitude) * cos (m_latitude));
2739 if (coordinates->the_mode == 30)
2740 /*
2741 Return the geocentric Sun/Moon elongation
2742 in degrees plus fraction.
2743 */
2744 return (TODEG (x4));
2745 /*
2746 Calculate geocentric rectangular ecliptic co-ordinates.
2747 */
2748 x1 = cos (m_latitude) * m_distance;
2749 xg = cos (m_longitude) * x1;
2750 yg = sin (m_longitude) * x1;
2751 zg = sin (m_latitude) * m_distance;
2752 /*
2753 Calculate geocentric rectangular equatorial co-ordinates.
2754 */
2755 sin_obliquity_of_ecliptic = sin (obliquity_of_ecliptic);
2756 cos_obliquity_of_ecliptic = cos (obliquity_of_ecliptic);
2757 xe = xg;
2758 ye =
2759 cos_obliquity_of_ecliptic * yg - sin_obliquity_of_ecliptic * zg;
2760 ze =
2761 sin_obliquity_of_ecliptic * yg + cos_obliquity_of_ecliptic * zg;
2762 /*
2763 Calculate the Moon's geocentric apparent right ascension angle.
2764 */
2765 right_ascension = FIXANGLE (TODEG (my_atan2 (ye, xe)));
2766 if (coordinates->the_mode == 28)
2767 /*
2768 Return the Moon's geocentric apparent right ascension angle
2769 in hours plus fraction.
2770 */
2771 return (DEG2HH (right_ascension));
2772 /*
2773 Calculate the Moon's geocentric apparent declination angle.
2774 */
2775 declination = asin (ze / m_distance);
2776 if (coordinates->the_mode == 25)
2777 /*
2778 Return the Moon's geocentric apparent declination angle
2779 in degrees plus fraction.
2780 */
2781 return (TODEG (declination));
2782 sin_latitude = sin (latitude);
2783 cos_latitude = cos (latitude);
2784 if (coordinates->the_mode >= 6)
2785 {
2786 /*
2787 Calculate the local sidereal time of moment.
2788 */
2789 local_sidereal_time = FIXANGLE (local_sidereal_time
2790 + DAY2DEG (time_offset -
2791 MM2DAY (tz_offset)));
2792 switch (coordinates->the_mode)
2793 {
2794 case 21:
2795 case 22:
2796 x2 = MY_PI - x4;
2797 if (coordinates->the_mode == 22)
2798 /*
2799 Return the geocentric apparent phase angle of the Moon
2800 in range +1.0...0.0.
2801 */
2802 return ((1.0 + cos (x2)) * 0.5);
2803 /*
2804 Calculate the Moon's geocentric apparent brightness
2805 in magnitude units plus fraction.
2806 */
2807 x2 = TODEG (x2);
2808 s_distance =
2809 sqrt (divisor_term * divisor_term +
2810 s_distance * s_distance);
2811 return (-21.62 + 5.0 * log10 (s_distance * m_distance) +
2812 0.026 * x2 + 0.000000004 * x2 * x2 * x2 * x2);
2813 case 23:
2814 case 24:
2815 x1 =
2816 TORAD (FIXANGLE (local_sidereal_time - right_ascension));
2817 if (coordinates->the_mode == 23)
2818 /*
2819 Return the geocentric apparent height of the Moon's
2820 center under or above the imaginary geocentric horizon
2821 (elevation angle) in degrees plus fraction.
2822 */
2823 return (TODEG (asin (sin_latitude * sin (declination)
2824 +
2825 cos_latitude * cos (declination) *
2826 cos (x1))));
2827 /*
2828 Return the Moon's geocentric apparent azimuth angle
2829 in degrees plus fraction.
2830 */
2831 return (FIXANGLE (TODEG (my_atan2 (sin (x1),
2832 sin_latitude * cos (x1) -
2833 cos_latitude *
2834 tan (declination)) +
2835 MY_PI)));
2836 case 32:
2837 /*
2838 Return the local sidereal time of moment
2839 in hours plus fraction.
2840 */
2841 return (DEG2HH (local_sidereal_time));
2842 default:
2843 /*
2844 Calculate topocentric rectangular equatorial co-ordinates.
2845 */
2846 x4 = gd_latitude2gc_latitude (latitude,
2847 (coordinates->
2848 meters_above_sea_level >
2849 0) ? coordinates->
2850 meters_above_sea_level : 0,
2851 &x3) / EQUATOR_EARTH_RADIUS;
2852 x1 = cos (x3) * x4;
2853 x2 = TORAD (local_sidereal_time);
2854 xe -= (cos (x2) * x1);
2855 ye -= (sin (x2) * x1);
2856 ze -= (sin (x3) * x4);
2857 }
2858 /*
2859 Calculate the topocentric Moon/Earth distance in Earth equator
2860 radii plus fraction.
2861 */
2862 m_distance = sqrt (xe * xe + ye * ye + ze * ze);
2863 switch (coordinates->the_mode)
2864 {
2865 case 6:
2866 /*
2867 Return the Moon's topocentric apparent horizontal parallax
2868 in degrees plus fraction.
2869 */
2870 return (TODEG (asin (1.0 / m_distance)));
2871 case 7:
2872 /*
2873 Return the Moon's topocentric apparent semidiameter
2874 in degrees plus fraction.
2875 */
2876 return (SS2DEG (936.85 * 60.0 / m_distance));
2877 case 16:
2878 /*
2879 Return the topocentric Moon/Earth distance in Earth equator
2880 radii plus fraction.
2881 */
2882 return (m_distance);
2883 default:
2884 ; /* Void, nothing to do! */
2885 }
2886 /*
2887 Calculate the topocentric right ascension angle,
2888 uncorrected for atmospheric refraction.
2889 */
2890 right_ascension = FIXANGLE (TODEG (my_atan2 (ye, xe)));
2891 /*
2892 Calculate the topocentric declination angle,
2893 uncorrected for atmospheric refraction.
2894 */
2895 declination = asin (ze / m_distance);
2896 switch (coordinates->the_mode)
2897 {
2898 case 8:
2899 case 9:
2900 case 13:
2901 case 14:
2902 case 17:
2903 /*
2904 Calculate the Moon's topocentric apparent ecliptic longitude.
2905 */
2906 m_longitude =
2907 FIXANGLE (TODEG
2908 (my_atan2
2909 (cos_obliquity_of_ecliptic *
2910 sin (TORAD (right_ascension)) +
2911 sin_obliquity_of_ecliptic * tan (declination),
2912 cos (TORAD (right_ascension)))));
2913 if (coordinates->the_mode == 13)
2914 /*
2915 Return the Moon's topocentric apparent ecliptic longitude
2916 in degrees plus fraction.
2917 */
2918 return (m_longitude);
2919 /*
2920 Calculate the Moon's topocentric apparent ecliptic latitude.
2921 */
2922 m_latitude =
2923 asin (cos_obliquity_of_ecliptic * sin (declination) -
2924 sin_obliquity_of_ecliptic * cos (declination) *
2925 sin (TORAD (right_ascension)));
2926 if (coordinates->the_mode == 14)
2927 /*
2928 Return the Moon's topocentric apparent ecliptic latitude
2929 in degrees plus fraction.
2930 */
2931 return (TODEG (m_latitude));
2932 /*
2933 Calculate the topocentric Sun/Moon elongation.
2934 */
2935 save_the_mode = coordinates->the_mode;
2936 /*
2937 Get the Sun's topocentric apparent ecliptic longitude.
2938 */
2939 coordinates->the_mode = 13;
2940 s_longitude =
2941 sun_rise_set (event, FALSE, day, month, year,
2942 coordinates);
2943 coordinates->the_mode = save_the_mode;
2944 x4 =
2945 acos (cos (TORAD (s_longitude - m_longitude)) *
2946 cos (m_latitude));
2947 if (coordinates->the_mode == 17)
2948 /*
2949 Return the topocentric Sun/Moon elongation
2950 in degrees plus fraction.
2951 */
2952 return (TODEG (x4));
2953 x2 = MY_PI - x4;
2954 if (coordinates->the_mode == 9)
2955 /*
2956 Return the topocentric apparent phase angle of the Moon
2957 in range +1.0...0.0.
2958 */
2959 return ((1.0 + cos (x2)) * 0.5);
2960 /*
2961 Calculate the Moon's topocentric apparent brightness
2962 in magnitude units plus fraction.
2963 */
2964 x2 = TODEG (x2);
2965 save_the_mode = coordinates->the_mode;
2966 /*
2967 Get the Sun's topocentric distance.
2968 */
2969 coordinates->the_mode = 15;
2970 s_distance =
2971 sun_rise_set (event, FALSE, day, month, year,
2972 coordinates);
2973 coordinates->the_mode = save_the_mode;
2974 return (-21.62
2975 + 5.0 * log10 (s_distance * m_distance)
2976 + 0.026 * x2 + 0.000000004 * x2 * x2 * x2 * x2);
2977 default:
2978 ; /* Void, nothing to do! */
2979 }
2980 /*
2981 Calculate the topocentric apparent height of the Moon's center
2982 under or above the sea-level horizon (elevation angle),
2983 corrected for atmospheric refraction.
2984 */
2985 x1 = TORAD (FIXANGLE (local_sidereal_time - right_ascension));
2986 x3 = cos (x1);
2987 x2 =
2988 asin (sin_latitude * sin (declination) +
2989 cos_latitude * cos (declination) * x3);
2990 /*
2991 Calculate the atmospheric refraction for correcting the Moon's
2992 topocentric true elevation angle to the Moon's topocentric
2993 apparent elevation angle.
2994 */
2995 x4 = atmospheric_refraction (x2, atm_pressure, atm_temperature);
2996 if (coordinates->the_mode == 18)
2997 /*
2998 Return the atmospheric refraction for correcting the Moon's
2999 topocentric true elevation angle to the Moon's topocentric
3000 apparent elevation angle in degrees plus fraction.
3001 */
3002 return (TODEG (x4));
3003 x2 += x4;
3004 if (x2 > MY_HALF_PI)
3005 x2 = MY_PI - x2;
3006 if (coordinates->the_mode == 10)
3007 /*
3008 Return the topocentric apparent height of the Moon's center
3009 under or above the sea-level horizon (elevation angle)
3010 in degrees plus fraction, corrected for atmospheric
3011 refraction.
3012 */
3013 return (TODEG (x2));
3014 /*
3015 Calculate the Moon's topocentric apparent azimuth angle
3016 in degrees plus fraction.
3017 */
3018 x3 = FIXANGLE (TODEG (my_atan2 (sin (x1),
3019 sin_latitude * x3 -
3020 cos_latitude *
3021 tan (declination)) + MY_PI));
3022 if (coordinates->the_mode == 11)
3023 return (x3);
3024 /*
3025 Reconvert the topocentric apparent azimuth angle and the
3026 topocentric apparent elevation angle to topocentric
3027 apparent right ascension angle and topocentric apparent
3028 declination angle.
3029 */
3030 x3 = TORAD (x3);
3031 xe = cos (x2);
3032 ye = sin (x2);
3033 ze = cos (x3);
3034 x1 = xe * ze;
3035 x2 = -xe * sin (x3);
3036 x3 = cos_latitude * ye - sin_latitude * x1;
3037 if (coordinates->the_mode == 12)
3038 /*
3039 Return the Moon's topocentric apparent declination angle
3040 in degrees plus fraction.
3041 */
3042 return (TODEG (asin (sin_latitude * ye + cos_latitude * x1)));
3043 /*
3044 Return the Moon's topocentric apparent right ascension angle
3045 in hours plus fraction.
3046 */
3047 return (DEG2HH
3048 (FIXANGLE
3049 (local_sidereal_time - TODEG (my_atan2 (x2, x3)))));
3050 }
3051 x4 = meridian_transit_time;
3052 meridian_transit_time =
3053 FIXANGLE (right_ascension - local_sidereal_time);
3054 if (iter_mt)
3055 {
3056 /*
3057 Calculate the meridian transit time (Moon's astronomical noon).
3058 */
3059 if ((i > 3)
3060 && (abs (abs (meridian_transit_time) - abs (x4)) >
3061 DEGS_PER_12_HOURS))
3062 return (SPECIAL_VALUE);
3063 i++;
3064 mjd1 = mjd1_0ut + DEG2DAY (meridian_transit_time);
3065 mjd2 = mjd2_0ut + DEG2DAY (meridian_transit_time);
3066 /*
3067 Exit iteration loop if the actually computed value is less
3068 than 1/1000.0 hours (== 3.6 seconds [0.015 degrees]) smaller
3069 the previously computed value. A ``next'' iteration loop
3070 would not change the result significantly -- the only change
3071 would be by the order of small fractions of seconds.
3072 Well, such a ``huge :-)'' epsilon value as well as a smaller
3073 epsilon value can be really critical if the iterated time is
3074 very near midnight, but what's perfect?
3075 */
3076 if (abs (abs (meridian_transit_time) - abs (x4)) > 0.015)
3077 goto LABEL_moon_iter_mt;
3078 return (DEG2HH (FIXANGLE (meridian_transit_time)));
3079 }
3080 /*
3081 Select the proper reference altitude
3082 according to the requested rise/set mode of the Moon.
3083 */
3084 switch (coordinates->the_mode)
3085 {
3086 case 2:
3087 case 4:
3088 /*
3089 Moon's center based altitude (0.0 degree == 0.0 radians).
3090 */
3091 if (adjust_value != DEGS_PER_24_HOURS)
3092 x1 = TORAD (adjust_value);
3093 else
3094 x1 = 0.0;
3095 break;
3096 case 3:
3097 case 5:
3098 /*
3099 Moon's upper limb based altitude.
3100 */
3101 if (adjust_value != DEGS_PER_24_HOURS)
3102 x1 = TORAD (adjust_value - m_semidiameter);
3103 else
3104 x1 = TORAD (-m_semidiameter);
3105 }
3106 switch (coordinates->the_mode)
3107 {
3108 case 4:
3109 case 5:
3110 /*
3111 Correct the true altitude giving the apparent altitude.
3112 */
3113 if (adjust_value != DEGS_PER_24_HOURS
3114 || atm_pressure != DEFAULT_PRESSURE
3115 || atm_temperature != DEFAULT_TEMPERATURE)
3116 x2 =
3117 atmospheric_refraction (x1, atm_pressure, atm_temperature);
3118 else
3119 /*
3120 Standard refraction of 34arc_minutes/60 == 0.00989019909463453388 radians.
3121 */
3122 x2 = 0.00989019909463453388;
3123 x1 -= x2;
3124 }
3125 x3 = (sin (x1 + m_parallax - geodetic_height)
3126 - sin_latitude * sin (declination))
3127 / (cos_latitude * cos (declination));
3128 if (abs (x3) > 1.0)
3129 {
3130 if (!n && (is_adjusted || the_time == SECS_PER_DAY))
3131 tt = x3;
3132 if (n > 2)
3133 {
3134 if (the_time != SECS_PER_DAY)
3135 break;
3136 if (x3 > 0.0)
3137 /*
3138 No rise/set for the given date
3139 because the Moon is always below the queried altitude.
3140 */
3141 return (2 * SPECIAL_VALUE);
3142 /*
3143 No rise/set for the given date
3144 because the Moon is always above the queried altitude.
3145 */
3146 return (3 * SPECIAL_VALUE);
3147 }
3148 /*
3149 Try a next iteration step by using the Moon's meridian transit time.
3150 */
3151 n++;
3152 x2 = DEG2DAY (meridian_transit_time);
3153 mjd1 = mjd1_0ut + x2;
3154 mjd2 = mjd2_0ut + x2;
3155 continue;
3156 }
3157 x2 = TODEG (acos (x3));
3158 x1 = the_time;
3159 if (event == RIse)
3160 the_time = meridian_transit_time - x2;
3161 else
3162 the_time = meridian_transit_time + x2;
3163 if (x1 != SECS_PER_DAY)
3164 {
3165 if (event == RIse)
3166 {
3167 if ((x1 < 0.0) && (the_time >= 0.0))
3168 {
3169 if (i > 2)
3170 {
3171 the_time = tt;
3172 tt = SECS_PER_DAY;
3173 break;
3174 }
3175 tt = x1;
3176 i++;
3177 }
3178 else
3179 if (!i
3180 && abs (abs (meridian_transit_time) - abs (x4)) >
3181 DEGS_PER_12_HOURS)
3182 {
3183 if (j > 2)
3184 /*
3185 No Moon rise for the given date.
3186 */
3187 return (SPECIAL_VALUE);
3188 j++;
3189 }
3190 }
3191 else
3192 if ((x1 < DEGS_PER_24_HOURS)
3193 && (the_time >= DEGS_PER_24_HOURS))
3194 {
3195 if (i > 2)
3196 {
3197 the_time = tt;
3198 tt = SECS_PER_DAY;
3199 break;
3200 }
3201 tt = the_time;
3202 i++;
3203 }
3204 else
3205 if (!i
3206 && abs (abs (meridian_transit_time) - abs (x4)) >
3207 DEGS_PER_12_HOURS)
3208 {
3209 if (j > 2)
3210 /*
3211 No Moon set for the given date.
3212 */
3213 return (SPECIAL_VALUE);
3214 j++;
3215 }
3216 }
3217 mjd1 = mjd1_0ut + DEG2DAY (the_time);
3218 mjd2 = mjd2_0ut + DEG2DAY (the_time);
3219 /*
3220 Exit iteration loop if the actually computed value is less
3221 than 1/1000.0 hours (== 3.6 seconds [0.015 degrees]) smaller
3222 the previously computed value. A ``next'' iteration loop
3223 would not change the result significantly -- the only change
3224 would be by the order of small fractions of seconds.
3225 Well, such a ``huge :-)'' epsilon value as well as a smaller
3226 epsilon value can be really critical if the iterated time is
3227 very near midnight, but what's perfect?
3228 */
3229 }
3230 while (abs (abs (the_time) - abs (x1)) > 0.015);
3231 /*
3232 And post-process the iterated rise/set/transit data of the Moon.
3233 */
3234 if (event == RIse)
3235 {
3236 if (!is_adjusted && (the_time < 0.0))
3237 {
3238 /*
3239 Try the day after the given date.
3240 */
3241 mjd1_0ut++;
3242 mjd2_0ut++;
3243 if (n > 1)
3244 tt = SECS_PER_DAY;
3245 }
3246 else
3247 ok = TRUE;
3248 }
3249 else if (!is_adjusted && (the_time >= DEGS_PER_24_HOURS))
3250 {
3251 /*
3252 Try the day prior the given date.
3253 */
3254 mjd1_0ut--;
3255 mjd2_0ut--;
3256 if (n > 1)
3257 tt = SECS_PER_DAY;
3258 }
3259 else
3260 ok = TRUE;
3261 if (!ok)
3262 {
3263 if (n == 1 || tt == SECS_PER_DAY)
3264 {
3265 /*
3266 Restart the calculation for the day prior or after the given date.
3267 */
3268 is_adjusted = TRUE;
3269 tt = the_time = SECS_PER_DAY;
3270 i = j = n = 0;
3271 }
3272 else if (n > 1)
3273 {
3274 if (tt > 1.0)
3275 return (2 * SPECIAL_VALUE);
3276 else if (tt < -1.0)
3277 return (3 * SPECIAL_VALUE);
3278 return (SPECIAL_VALUE);
3279 }
3280 else if (!n)
3281 return (SPECIAL_VALUE);
3282 }
3283 }
3284 /*
3285 Ok, all iteration-stuff done, so perform a final check whether the
3286 calculated rise/set/transit data of the Moon is valid for the given date!
3287 */
3288 if (is_adjusted && (n != 1) && (tt != SECS_PER_DAY))
3289 {
3290 if (n)
3291 {
3292 if (tt > 1.0)
3293 return (2 * SPECIAL_VALUE);
3294 else if (tt < -1.0)
3295 return (3 * SPECIAL_VALUE);
3296 }
3297 return (SPECIAL_VALUE);
3298 }
3299 if (!is_adjusted && (n > 1))
3300 {
3301 if (x3 > 1.0)
3302 return (2 * SPECIAL_VALUE);
3303 else if (x3 < -1.0)
3304 return (3 * SPECIAL_VALUE);
3305 return (SPECIAL_VALUE);
3306 }
3307 /*
3308 Yes, the rise/set/transit data of the Moon is valid for the given date,
3309 so correct the calculated (position) data to circle range.
3310 */
3311 while (the_time >= DEGS_PER_24_HOURS)
3312 the_time -= DEGS_PER_24_HOURS;
3313 while (the_time < 0.0)
3314 the_time += DEGS_PER_24_HOURS;
3315
3316 return (DEG2HH (the_time));
3317 }
3318
3319
3320
3321 int
moondisk(is_full_new,day,month,year,hour,min)3322 moondisk (is_full_new, day, month, year, hour, min)
3323 Bool *is_full_new;
3324 int day;
3325 int month;
3326 int year;
3327 const int hour;
3328 const int min;
3329 /*!
3330 Calculates the approximate illuminated fraction of the Moon's disk for
3331 the given Julian/Gregorian date (range 00010101...99991231) and returns
3332 it represented by an integer value within the range of -100...0...+100,
3333 which has a negative sign in case the Moon wanes, otherwise the sign
3334 is positive. The New Moon phase is around the 0.0 value and the Full
3335 Moon phase is around the +/-100 value. The address of IS_FULL_NEW is
3336 set to TRUE in case the Moon is in the "Full" or in the "New" phase
3337 at that date, otherwise to FALSE. Calculations are made for a definite
3338 meridian expressed as a time value in HOUR and MIN. If HOUR and MIN are
3339 set to zero, calculations are made for Universal Time (UTC/GMT). If HOUR
3340 and MIN have a positive sign, UTC/GMT calculations are made for meridians
3341 East of Greenwich, otherwise for meridians West of Greenwich.
3342 For simplification and speed-up purposes, this function assumes the
3343 illumination change rate of the Moon's disk as the linear result as
3344 expressed by a quadratic sine curve, where the +50% value (the Moon's
3345 First Quarter) and the -50% value (the Moon's Last Quarter) are treated
3346 as the mid-points respective in between successive New and Full Moon
3347 phases. This modell does not match the true phenomenon very well, but
3348 it produces acceptable results. The maximum error is less than 15%
3349 absolute (about 7% relative) near curve regions where the absolute
3350 gradient of the curve is at maximum (the +|-50% value).
3351 */
3352 {
3353 auto double new_moon;
3354 auto double full_moon;
3355 auto double len_lunation_half;
3356 auto double illuminated_fraction;
3357 auto Ulint lunation;
3358 auto Ulint mjd = date2num (day, month, year);
3359
3360
3361 *is_full_new = FALSE;
3362 if (mjd <= 13L)
3363 {
3364 lunation = 1L;
3365 new_moon =
3366 moonphase (MPHASE_NEW, FALSE, NULL, &lunation, &day, &month, &year,
3367 hour, min);
3368 /*
3369 The last Full Moon phase prior the Christian Era was on
3370 29-Dec-0001 BCE, 14:55 UT, approximately.
3371 */
3372 full_moon = -2.3784722;
3373 }
3374 else
3375 {
3376 lunation = 0L;
3377 new_moon =
3378 moonphase (MPHASE_NEW, FALSE, NULL, &lunation, &day, &month, &year,
3379 -hour, -min);
3380 lunation--;
3381 full_moon =
3382 moonphase (MPHASE_FUL, FALSE, NULL, &lunation, &day, &month, &year,
3383 -hour, -min);
3384 if (date2num (day, month, year) >= mjd)
3385 new_moon =
3386 moonphase (MPHASE_NEW, FALSE, NULL, &lunation, &day, &month, &year,
3387 -hour, -min);
3388 }
3389 len_lunation_half = abs (full_moon - new_moon);
3390 if (mjd > (Ulint) new_moon)
3391 len_lunation_half = (mjd - new_moon) / len_lunation_half;
3392 else
3393 len_lunation_half =
3394 (-len_lunation_half + (mjd - full_moon)) / len_lunation_half;
3395 illuminated_fraction = sin (len_lunation_half * MY_HALF_PI);
3396 illuminated_fraction *= illuminated_fraction;
3397 if (len_lunation_half < 0.0)
3398 illuminated_fraction = -illuminated_fraction;
3399 if (mjd == (Ulint) new_moon || mjd == (Ulint) full_moon)
3400 *is_full_new = TRUE;
3401
3402 return ((int) ROUND (illuminated_fraction * 100.0));
3403 }
3404
3405
3406
3407 void
draw_moon(age,lines,string)3408 draw_moon (age, lines, string)
3409 const int age;
3410 const int lines;
3411 char **string;
3412 /*!
3413 Creates a text graphics image of the Moon according to its `age', which is
3414 expressed as the illuminated fraction of the Moon's disk as a value within
3415 the range of -100...0...+100, which has a negative sign in case the Moon
3416 wanes, otherwise the sign is positive. Uses the delivered `&string' for
3417 storing this image (each line with a leading newline character represented
3418 by Gcal's RC_NL_CHAR=='~') and returns it. The caller has to guarantee
3419 that enough `&string' space is allocated, which means there must be two
3420 Bytes allocated minimum in this case, because this function increases
3421 the `&string' properly by reallocating it internally when used within Gcal.
3422 */
3423 {
3424 register Uint slen = 0;
3425 register int i;
3426 register int j;
3427 register int k;
3428 register int end;
3429 auto double counter = (double) age;
3430 auto const double step = 2.0 / (double) lines;
3431 auto double squisher;
3432 auto double horizon;
3433 auto double terminator;
3434 static char buffer[80]; /* Proper if [MOONIMAGE_MIN]6<=lines<=30[MOONIMAGE_MAX] */
3435 auto char *ptr_buffer;
3436
3437
3438 **string = '\0';
3439 if (counter < 0.0)
3440 counter = -counter;
3441 if (age < 0)
3442 counter = (100.0 - (double) (((Uint) counter >> 1) + 50)) / 100.0;
3443 else
3444 counter /= 200.0;
3445 squisher = cos (counter * MY_TWO_PI);
3446 for (counter = 0.93; counter > -1.0; counter -= step)
3447 {
3448 sprintf (buffer, "%c", RC_NL_CHAR);
3449 strcat (*string, buffer);
3450 slen++;
3451 horizon = sqrt (1.0 - counter * counter);
3452 i = j = moon_charpos (horizon, lines + 6);
3453 for (ptr_buffer = buffer; i--;)
3454 *ptr_buffer++ = ' ';
3455 i = j;
3456 slen += (i + 1);
3457 buffer[i] = *MOONIMAGE_REDGE;
3458 buffer[i + 1] = '\0';
3459 j = moon_charpos (-horizon, lines + 6);
3460 buffer[j] = *MOONIMAGE_LEDGE;
3461 terminator = horizon * squisher;
3462 if (abs (age) > 6)
3463 {
3464 k = moon_charpos (terminator, lines + 6);
3465 if (age > 0)
3466 {
3467 end = i;
3468 i = k;
3469 }
3470 else
3471 {
3472 i = j;
3473 end = k;
3474 }
3475 while (i <= end)
3476 buffer[i++] = *MOONIMAGE_BRIGHT;
3477 }
3478 if (slen + 1 >= maxlen_max)
3479 resize_all_strings (maxlen_max << 1, FALSE, __FILE__,
3480 (long) __LINE__);
3481 strcat (*string, buffer);
3482 }
3483 }
3484 #endif /* USE_RC */
3485