The SSH host key has changed on 8 April, 2022 to this one: SHA256:573uTBSeh74kvOo0HJXi5ijdzRm8me27suzNEDlGyrQ
A GLib style wrapper library around the Swiss Ephemeris library, created by Astrodienst.
You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
 
 
 
 
swe-glib/swe/src/swecl.c

5410 lines
183 KiB

/* SWISSEPH
$Header: /home/dieter/sweph/RCS/swecl.c,v 1.75 2008/08/26 07:23:27 dieter Exp $
Ephemeris computations
Author: Dieter Koch
************************************************************/
/* Copyright (C) 1997 - 2008 Astrodienst AG, Switzerland. All rights reserved.
License conditions
------------------
This file is part of Swiss Ephemeris.
Swiss Ephemeris is distributed with NO WARRANTY OF ANY KIND. No author
or distributor accepts any responsibility for the consequences of using it,
or for whether it serves any particular purpose or works at all, unless he
or she says so in writing.
Swiss Ephemeris is made available by its authors under a dual licensing
system. The software developer, who uses any part of Swiss Ephemeris
in his or her software, must choose between one of the two license models,
which are
a) GNU public license version 2 or later
b) Swiss Ephemeris Professional License
The choice must be made before the software developer distributes software
containing parts of Swiss Ephemeris to others, and before any public
service using the developed software is activated.
If the developer choses the GNU GPL software license, he or she must fulfill
the conditions of that license, which includes the obligation to place his
or her whole software project under the GNU GPL or a compatible license.
See http://www.gnu.org/licenses/old-licenses/gpl-2.0.html
If the developer choses the Swiss Ephemeris Professional license,
he must follow the instructions as found in http://www.astro.com/swisseph/
and purchase the Swiss Ephemeris Professional Edition from Astrodienst
and sign the corresponding license contract.
The License grants you the right to use, copy, modify and redistribute
Swiss Ephemeris, but only under certain conditions described in the License.
Among other things, the License requires that the copyright notices and
this notice be preserved on all copies.
Authors of the Swiss Ephemeris: Dieter Koch and Alois Treindl
The authors of Swiss Ephemeris have no control or influence over any of
the derived works, i.e. over software or services created by other
programmers which use Swiss Ephemeris functions.
The names of the authors or of the copyright holder (Astrodienst) must not
be used for promoting any software, product or service which uses or contains
the Swiss Ephemeris. This copyright notice is the ONLY place where the
names of the authors can legally appear, except in cases where they have
given special permission in writing.
The trademarks 'Swiss Ephemeris' and 'Swiss Ephemeris inside' may be used
for promoting such software, products or services.
*/
#include "swejpl.h"
#include "swephexp.h"
#include "sweph.h"
#include "swephlib.h"
#include <time.h>
#define SEFLG_EPHMASK (SEFLG_JPLEPH|SEFLG_SWIEPH|SEFLG_MOSEPH)
static int find_maximum(double y00, double y11, double y2, double dx,
double *dxret, double *yret);
static int find_zero(double y00, double y11, double y2, double dx,
double *dxret, double *dxret2);
static double calc_dip(double geoalt, double atpress, double attemp, double lapse_rate);
static double calc_astronomical_refr(double geoalt,double atpress, double attemp);
static TLS double const_lapse_rate = SE_LAPSE_RATE; /* for refraction */
#if 0
#define DSUN (1391978489.9 / AUNIT) /* this value is consistent with
* 959.63 arcsec at AU distance (Astr. Alm.) */
#else
#define DSUN (1392000000.0 / AUNIT)
#endif
#define DMOON (3476300.0 / AUNIT)
#define DEARTH (6378140.0 * 2 / AUNIT)
#define RSUN (DSUN / 2)
#define RMOON (DMOON / 2)
#define REARTH (DEARTH / 2)
/*#define SEI_OCC_FAST (16 * 1024L)*/
static int32 eclipse_where( double tjd_ut, int32 ipl, char *starname, int32 ifl, double *geopos,
double *dcore, char *serr);
static int32 eclipse_how( double tjd_ut, int32 ipl, char *starname, int32 ifl,
double geolon, double geolat, double geohgt,
double *attr, char *serr);
static int32 eclipse_when_loc(double tjd_start, int32 ifl, double *geopos,
double *tret, double *attr, AS_BOOL backward, char *serr);
static int32 occult_when_loc(double tjd_start, int32 ipl, char *starname, int32 ifl,
double *geopos, double *tret, double *attr, AS_BOOL backward, char *serr);
static int32 lun_eclipse_how(double tjd_ut, int32 ifl, double *attr,
double *dcore, char *serr);
static int32 calc_mer_trans(
double tjd_ut, int32 ipl, int32 epheflag, int32 rsmi,
double *geopos,
char *starname,
double *tret,
char *serr);
static int32 calc_planet_star(double tjd_et, int32 ipl, char *starname, int32 iflag, double *x, char *serr);
struct saros_data {int series_no; double tstart;};
#define SAROS_CYCLE 6585.3213
#define NSAROS_SOLAR 181
struct saros_data saros_data_solar[NSAROS_SOLAR] = {
{0, 641886.5}, /* 23 May -2955 */
{1, 672214.5}, /* 04 Jun -2872 */
{2, 676200.5}, /* 04 May -2861 */
{3, 693357.5}, /* 24 Apr -2814 */
{4, 723685.5}, /* 06 May -2731 */
{5, 727671.5}, /* 04 Apr -2720 */
{6, 744829.5}, /* 27 Mar -2673 */
{7, 775157.5}, /* 08 Apr -2590 */
{8, 779143.5}, /* 07 Mar -2579 */
{9, 783131.5}, /* 06 Feb -2568 */
{10, 820044.5}, /* 28 Feb -2467 */
{11, 810859.5}, /* 06 Jan -2492 */
{12, 748993.5}, /* 20 Aug -2662 */
{13, 792492.5}, /* 23 Sep -2543 */
{14, 789892.5}, /* 11 Aug -2550 */
{15, 787294.5}, /* 01 Jul -2557 */
{16, 824207.5}, /* 23 Jul -2456 */
{17, 834779.5}, /* 03 Jul -2427 */
{18, 838766.5}, /* 02 Jun -2416 */
{19, 869094.5}, /* 15 Jun -2333 */
{20, 886251.5}, /* 05 Jun -2286 */
{21, 890238.5}, /* 05 May -2275 */
{22, 927151.5}, /* 28 May -2174 */
{23, 937722.5}, /* 07 May -2145 */
{24, 941709.5}, /* 06 Apr -2134 */
{25, 978623.5}, /* 30 Apr -2033 */
{26, 989194.5}, /* 08 Apr -2004 */
{27, 993181.5}, /* 09 Mar -1993 */
{28, 1023510.5}, /* 22 Mar -1910 */
{29, 1034081.5}, /* 01 Mar -1881 */
{30, 972214.5}, /* 12 Oct -2051 */
{31, 1061811.5}, /* 31 Jan -1805 */
{32, 1006529.5}, /* 24 Sep -1957 */
{33, 997345.5}, /* 02 Aug -1982 */
{34, 1021088.5}, /* 04 Aug -1917 */
{35, 1038245.5}, /* 25 Jul -1870 */
{36, 1042231.5}, /* 23 Jun -1859 */
{37, 1065974.5}, /* 25 Jun -1794 */
{38, 1089716.5}, /* 26 Jun -1729 */
{39, 1093703.5}, /* 26 May -1718 */
{40, 1117446.5}, /* 28 May -1653 */
{41, 1141188.5}, /* 28 May -1588 */
{42, 1145175.5}, /* 28 Apr -1577 */
{43, 1168918.5}, /* 29 Apr -1512 */
{44, 1192660.5}, /* 30 Apr -1447 */
{45, 1196647.5}, /* 30 Mar -1436 */
{46, 1220390.5}, /* 01 Apr -1371 */
{47, 1244132.5}, /* 02 Apr -1306 */
{48, 1234948.5}, /* 08 Feb -1331 */
{49, 1265277.5}, /* 22 Feb -1248 */
{50, 1282433.5}, /* 11 Feb -1201 */
{51, 1207395.5}, /* 02 Sep -1407 */
{52, 1217968.5}, /* 14 Aug -1378 */
{53, 1254881.5}, /* 06 Sep -1277 */
{54, 1252282.5}, /* 25 Jul -1284 */
{55, 1262855.5}, /* 06 Jul -1255 */
{56, 1293182.5}, /* 17 Jul -1172 */
{57, 1297169.5}, /* 17 Jun -1161 */
{58, 1314326.5}, /* 07 Jun -1114 */
{59, 1344654.5}, /* 19 Jun -1031 */
{60, 1348640.5}, /* 18 May -1020 */
{61, 1365798.5}, /* 10 May -0973 */
{62, 1396126.5}, /* 22 May -0890 */
{63, 1400112.5}, /* 20 Apr -0879 */
{64, 1417270.5}, /* 11 Apr -0832 */
{65, 1447598.5}, /* 24 Apr -0749 */
{66, 1444999.5}, /* 12 Mar -0756 */
{67, 1462157.5}, /* 04 Mar -0709 */
{68, 1492485.5}, /* 16 Mar -0626 */
{69, 1456959.5}, /* 09 Dec -0724 */
{70, 1421434.5}, /* 05 Sep -0821 */
{71, 1471518.5}, /* 19 Oct -0684 */
{72, 1455748.5}, /* 16 Aug -0727 */
{73, 1466320.5}, /* 27 Jul -0698 */
{74, 1496648.5}, /* 08 Aug -0615 */
{75, 1500634.5}, /* 07 Jul -0604 */
{76, 1511207.5}, /* 18 Jun -0575 */
{77, 1548120.5}, /* 11 Jul -0474 */
{78, 1552106.5}, /* 09 Jun -0463 */
{79, 1562679.5}, /* 21 May -0434 */
{80, 1599592.5}, /* 13 Jun -0333 */
{81, 1603578.5}, /* 12 May -0322 */
{82, 1614150.5}, /* 22 Apr -0293 */
{83, 1644479.5}, /* 05 May -0210 */
{84, 1655050.5}, /* 14 Apr -0181 */
{85, 1659037.5}, /* 14 Mar -0170 */
{86, 1695950.5}, /* 06 Apr -0069 */
{87, 1693351.5}, /* 23 Feb -0076 */
{88, 1631484.5}, /* 06 Oct -0246 */
{89, 1727666.5}, /* 04 Feb 0018 */
{90, 1672384.5}, /* 28 Sep -0134 */
{91, 1663200.5}, /* 06 Aug -0159 */
{92, 1693529.5}, /* 19 Aug -0076 */
{93, 1710685.5}, /* 09 Aug -0029 */
{94, 1714672.5}, /* 09 Jul -0018 */
{95, 1738415.5}, /* 11 Jul 0047 */
{96, 1755572.5}, /* 01 Jul 0094 */
{97, 1766144.5}, /* 11 Jun 0123 */
{98, 1789887.5}, /* 12 Jun 0188 */
{99, 1807044.5}, /* 03 Jun 0235 */
{100, 1817616.5}, /* 13 May 0264 */
{101, 1841359.5}, /* 15 May 0329 */
{102, 1858516.5}, /* 05 May 0376 */
{103, 1862502.5}, /* 04 Apr 0387 */
{104, 1892831.5}, /* 17 Apr 0470 */
{105, 1903402.5}, /* 27 Mar 0499 */
{106, 1887633.5}, /* 23 Jan 0456 */
{107, 1924547.5}, /* 15 Feb 0557 */
{108, 1921948.5}, /* 04 Jan 0550 */
{109, 1873251.5}, /* 07 Sep 0416 */
{110, 1890409.5}, /* 30 Aug 0463 */
{111, 1914151.5}, /* 30 Aug 0528 */
{112, 1918138.5}, /* 31 Jul 0539 */
{113, 1935296.5}, /* 22 Jul 0586 */
{114, 1959038.5}, /* 23 Jul 0651 */
{115, 1963024.5}, /* 21 Jun 0662 */
{116, 1986767.5}, /* 23 Jun 0727 */
{117, 2010510.5}, /* 24 Jun 0792 */
{118, 2014496.5}, /* 24 May 0803 */
{119, 2031654.5}, /* 15 May 0850 */
{120, 2061982.5}, /* 27 May 0933 */
{121, 2065968.5}, /* 25 Apr 0944 */
{122, 2083126.5}, /* 17 Apr 0991 */
{123, 2113454.5}, /* 29 Apr 1074 */
{124, 2104269.5}, /* 06 Mar 1049 */
{125, 2108256.5}, /* 04 Feb 1060 */
{126, 2151755.5}, /* 10 Mar 1179 */
{127, 2083302.5}, /* 10 Oct 0991 */
{128, 2080704.5}, /* 29 Aug 0984 */
{129, 2124203.5}, /* 03 Oct 1103 */
{130, 2121603.5}, /* 20 Aug 1096 */
{131, 2132176.5}, /* 01 Aug 1125 */
{132, 2162504.5}, /* 13 Aug 1208 */
{133, 2166490.5}, /* 13 Jul 1219 */
{134, 2177062.5}, /* 22 Jun 1248 */
{135, 2207390.5}, /* 05 Jul 1331 */
{136, 2217962.5}, /* 14 Jun 1360 */
{137, 2228534.5}, /* 25 May 1389 */
{138, 2258862.5}, /* 06 Jun 1472 */
{139, 2269434.5}, /* 17 May 1501 */
{140, 2273421.5}, /* 16 Apr 1512 */
{141, 2310334.5}, /* 19 May 1613 */
{142, 2314320.5}, /* 17 Apr 1624 */
{143, 2311722.5}, /* 07 Mar 1617 */
{144, 2355221.5}, /* 11 Apr 1736 */
{145, 2319695.5}, /* 04 Jan 1639 */
{146, 2284169.5}, /* 19 Sep 1541 */
{147, 2314498.5}, /* 12 Oct 1624 */
{148, 2325069.5}, /* 21 Sep 1653 */
{149, 2329056.5}, /* 21 Aug 1664 */
{150, 2352799.5}, /* 24 Aug 1729 */
{151, 2369956.5}, /* 14 Aug 1776 */
{152, 2380528.5}, /* 26 Jul 1805 */
{153, 2404271.5}, /* 28 Jul 1870 */
{154, 2421428.5}, /* 19 Jul 1917 */
{155, 2425414.5}, /* 17 Jun 1928 */
{156, 2455743.5}, /* 01 Jul 2011 */
{157, 2472900.5}, /* 21 Jun 2058 */
{158, 2476886.5}, /* 20 May 2069 */
{159, 2500629.5}, /* 23 May 2134 */
{160, 2517786.5}, /* 13 May 2181 */
{161, 2515187.5}, /* 01 Apr 2174 */
{162, 2545516.5}, /* 15 Apr 2257 */
{163, 2556087.5}, /* 25 Mar 2286 */
{164, 2487635.5}, /* 24 Oct 2098 */
{165, 2504793.5}, /* 16 Oct 2145 */
{166, 2535121.5}, /* 29 Oct 2228 */
{167, 2525936.5}, /* 06 Sep 2203 */
{168, 2543094.5}, /* 28 Aug 2250 */
{169, 2573422.5}, /* 10 Sep 2333 */
{170, 2577408.5}, /* 09 Aug 2344 */
{171, 2594566.5}, /* 01 Aug 2391 */
{172, 2624894.5}, /* 13 Aug 2474 */
{173, 2628880.5}, /* 12 Jul 2485 */
{174, 2646038.5}, /* 04 Jul 2532 */
{175, 2669780.5}, /* 05 Jul 2597 */
{176, 2673766.5}, /* 04 Jun 2608 */
{177, 2690924.5}, /* 27 May 2655 */
{178, 2721252.5}, /* 09 Jun 2738 */
{179, 2718653.5}, /* 28 Apr 2731 */
{180, 2729226.5}, /* 08 Apr 2760 */
};
#define NSAROS_LUNAR 180
struct saros_data saros_data_lunar[NSAROS_LUNAR] = {
{1, 782437.5}, /* 14 Mar -2570 */
{2, 799593.5}, /* 03 Mar -2523 */
{3, 783824.5}, /* 30 Dec -2567 */
{4, 754884.5}, /* 06 Oct -2646 */
{5, 824724.5}, /* 22 Dec -2455 */
{6, 762857.5}, /* 04 Aug -2624 */
{7, 773430.5}, /* 16 Jul -2595 */
{8, 810343.5}, /* 08 Aug -2494 */
{9, 807743.5}, /* 26 Jun -2501 */
{10, 824901.5}, /* 17 Jun -2454 */
{11, 855229.5}, /* 29 Jun -2371 */
{12, 859215.5}, /* 28 May -2360 */
{13, 876373.5}, /* 20 May -2313 */
{14, 906701.5}, /* 01 Jun -2230 */
{15, 910687.5}, /* 30 Apr -2219 */
{16, 927845.5}, /* 21 Apr -2172 */
{17, 958173.5}, /* 04 May -2089 */
{18, 962159.5}, /* 02 Apr -2078 */
{19, 979317.5}, /* 24 Mar -2031 */
{20, 1009645.5}, /* 05 Apr -1948 */
{21, 1007046.5}, /* 22 Feb -1955 */
{22, 1017618.5}, /* 02 Feb -1926 */
{23, 1054531.5}, /* 25 Feb -1825 */
{24, 979493.5}, /* 16 Sep -2031 */
{25, 976895.5}, /* 06 Aug -2038 */
{26, 1020394.5}, /* 09 Sep -1919 */
{27, 1017794.5}, /* 28 Jul -1926 */
{28, 1028367.5}, /* 09 Jul -1897 */
{29, 1058695.5}, /* 21 Jul -1814 */
{30, 1062681.5}, /* 19 Jun -1803 */
{31, 1073253.5}, /* 30 May -1774 */
{32, 1110167.5}, /* 23 Jun -1673 */
{33, 1114153.5}, /* 22 May -1662 */
{34, 1131311.5}, /* 13 May -1615 */
{35, 1161639.5}, /* 25 May -1532 */
{36, 1165625.5}, /* 24 Apr -1521 */
{37, 1176197.5}, /* 03 Apr -1492 */
{38, 1213111.5}, /* 27 Apr -1391 */
{39, 1217097.5}, /* 26 Mar -1380 */
{40, 1221084.5}, /* 24 Feb -1369 */
{41, 1257997.5}, /* 18 Mar -1268 */
{42, 1255398.5}, /* 04 Feb -1275 */
{43, 1186946.5}, /* 07 Sep -1463 */
{44, 1283128.5}, /* 06 Jan -1199 */
{45, 1227845.5}, /* 29 Aug -1351 */
{46, 1225247.5}, /* 19 Jul -1358 */
{47, 1255575.5}, /* 31 Jul -1275 */
{48, 1272732.5}, /* 21 Jul -1228 */
{49, 1276719.5}, /* 21 Jun -1217 */
{50, 1307047.5}, /* 03 Jul -1134 */
{51, 1317619.5}, /* 13 Jun -1105 */
{52, 1328191.5}, /* 23 May -1076 */
{53, 1358519.5}, /* 05 Jun -0993 */
{54, 1375676.5}, /* 26 May -0946 */
{55, 1379663.5}, /* 25 Apr -0935 */
{56, 1409991.5}, /* 07 May -0852 */
{57, 1420562.5}, /* 16 Apr -0823 */
{58, 1424549.5}, /* 16 Mar -0812 */
{59, 1461463.5}, /* 09 Apr -0711 */
{60, 1465449.5}, /* 08 Mar -0700 */
{61, 1436509.5}, /* 13 Dec -0780 */
{62, 1493179.5}, /* 08 Feb -0624 */
{63, 1457653.5}, /* 03 Nov -0722 */
{64, 1435298.5}, /* 20 Aug -0783 */
{65, 1452456.5}, /* 11 Aug -0736 */
{66, 1476198.5}, /* 12 Aug -0671 */
{67, 1480184.5}, /* 11 Jul -0660 */
{68, 1503928.5}, /* 14 Jul -0595 */
{69, 1527670.5}, /* 15 Jul -0530 */
{70, 1531656.5}, /* 13 Jun -0519 */
{71, 1548814.5}, /* 04 Jun -0472 */
{72, 1579142.5}, /* 17 Jun -0389 */
{73, 1583128.5}, /* 16 May -0378 */
{74, 1600286.5}, /* 07 May -0331 */
{75, 1624028.5}, /* 08 May -0266 */
{76, 1628015.5}, /* 07 Apr -0255 */
{77, 1651758.5}, /* 09 Apr -0190 */
{78, 1675500.5}, /* 10 Apr -0125 */
{79, 1672901.5}, /* 27 Feb -0132 */
{80, 1683474.5}, /* 07 Feb -0103 */
{81, 1713801.5}, /* 19 Feb -0020 */
{82, 1645349.5}, /* 21 Sep -0208 */
{83, 1649336.5}, /* 22 Aug -0197 */
{84, 1686249.5}, /* 13 Sep -0096 */
{85, 1683650.5}, /* 02 Aug -0103 */
{86, 1694222.5}, /* 13 Jul -0074 */
{87, 1731136.5}, /* 06 Aug 0027 */
{88, 1735122.5}, /* 05 Jul 0038 */
{89, 1745694.5}, /* 15 Jun 0067 */
{90, 1776022.5}, /* 27 Jun 0150 */
{91, 1786594.5}, /* 07 Jun 0179 */
{92, 1797166.5}, /* 17 May 0208 */
{93, 1827494.5}, /* 30 May 0291 */
{94, 1838066.5}, /* 09 May 0320 */
{95, 1848638.5}, /* 19 Apr 0349 */
{96, 1878966.5}, /* 01 May 0432 */
{97, 1882952.5}, /* 31 Mar 0443 */
{98, 1880354.5}, /* 18 Feb 0436 */
{99, 1923853.5}, /* 24 Mar 0555 */
{100, 1881741.5}, /* 06 Dec 0439 */
{101, 1852801.5}, /* 11 Sep 0360 */
{102, 1889715.5}, /* 05 Oct 0461 */
{103, 1893701.5}, /* 03 Sep 0472 */
{104, 1897688.5}, /* 04 Aug 0483 */
{105, 1928016.5}, /* 16 Aug 0566 */
{106, 1938588.5}, /* 27 Jul 0595 */
{107, 1942575.5}, /* 26 Jun 0606 */
{108, 1972903.5}, /* 08 Jul 0689 */
{109, 1990059.5}, /* 27 Jun 0736 */
{110, 1994046.5}, /* 28 May 0747 */
{111, 2024375.5}, /* 10 Jun 0830 */
{112, 2034946.5}, /* 20 May 0859 */
{113, 2045518.5}, /* 29 Apr 0888 */
{114, 2075847.5}, /* 13 May 0971 */
{115, 2086418.5}, /* 21 Apr 1000 */
{116, 2083820.5}, /* 11 Mar 0993 */
{117, 2120733.5}, /* 03 Apr 1094 */
{118, 2124719.5}, /* 02 Mar 1105 */
{119, 2062852.5}, /* 14 Oct 0935 */
{120, 2086596.5}, /* 16 Oct 1000 */
{121, 2103752.5}, /* 06 Oct 1047 */
{122, 2094568.5}, /* 14 Aug 1022 */
{123, 2118311.5}, /* 16 Aug 1087 */
{124, 2142054.5}, /* 17 Aug 1152 */
{125, 2146040.5}, /* 17 Jul 1163 */
{126, 2169783.5}, /* 18 Jul 1228 */
{127, 2186940.5}, /* 09 Jul 1275 */
{128, 2197512.5}, /* 18 Jun 1304 */
{129, 2214670.5}, /* 10 Jun 1351 */
{130, 2238412.5}, /* 10 Jun 1416 */
{131, 2242398.5}, /* 10 May 1427 */
{132, 2266142.5}, /* 12 May 1492 */
{133, 2289884.5}, /* 13 May 1557 */
{134, 2287285.5}, /* 01 Apr 1550 */
{135, 2311028.5}, /* 13 Apr 1615 */
{136, 2334770.5}, /* 13 Apr 1680 */
{137, 2292659.5}, /* 17 Dec 1564 */
{138, 2276890.5}, /* 15 Oct 1521 */
{139, 2326974.5}, /* 09 Dec 1658 */
{140, 2304619.5}, /* 25 Sep 1597 */
{141, 2308606.5}, /* 25 Aug 1608 */
{142, 2345520.5}, /* 19 Sep 1709 */
{143, 2349506.5}, /* 18 Aug 1720 */
{144, 2360078.5}, /* 29 Jul 1749 */
{145, 2390406.5}, /* 11 Aug 1832 */
{146, 2394392.5}, /* 11 Jul 1843 */
{147, 2411550.5}, /* 02 Jul 1890 */
{148, 2441878.5}, /* 15 Jul 1973 */
{149, 2445864.5}, /* 13 Jun 1984 */
{150, 2456437.5}, /* 25 May 2013 */
{151, 2486765.5}, /* 06 Jun 2096 */
{152, 2490751.5}, /* 07 May 2107 */
{153, 2501323.5}, /* 16 Apr 2136 */
{154, 2538236.5}, /* 10 May 2237 */
{155, 2529052.5}, /* 18 Mar 2212 */
{156, 2473771.5}, /* 08 Nov 2060 */
{157, 2563367.5}, /* 01 Mar 2306 */
{158, 2508085.5}, /* 21 Oct 2154 */
{159, 2505486.5}, /* 09 Sep 2147 */
{160, 2542400.5}, /* 03 Oct 2248 */
{161, 2546386.5}, /* 02 Sep 2259 */
{162, 2556958.5}, /* 12 Aug 2288 */
{163, 2587287.5}, /* 27 Aug 2371 */
{164, 2597858.5}, /* 05 Aug 2400 */
{165, 2601845.5}, /* 06 Jul 2411 */
{166, 2632173.5}, /* 18 Jul 2494 */
{167, 2649330.5}, /* 09 Jul 2541 */
{168, 2653317.5}, /* 08 Jun 2552 */
{169, 2683645.5}, /* 22 Jun 2635 */
{170, 2694217.5}, /* 01 Jun 2664 */
{171, 2698203.5}, /* 01 May 2675 */
{172, 2728532.5}, /* 15 May 2758 */
{173, 2739103.5}, /* 24 Apr 2787 */
{174, 2683822.5}, /* 16 Dec 2635 */
{175, 2740492.5}, /* 11 Feb 2791 */
{176, 2724722.5}, /* 09 Dec 2747 */
{177, 2708952.5}, /* 05 Oct 2704 */
{178, 2732695.5}, /* 07 Oct 2769 */
{179, 2749852.5}, /* 27 Sep 2816 */
{180, 2753839.5}, /* 28 Aug 2827 */
};
/* Computes geographic location and type of solar eclipse
* for a given tjd
* iflag: to indicate ephemeris to be used
* (SEFLG_JPLEPH, SEFLG_SWIEPH, SEFLG_MOSEPH)
*
* Algorithms for the central line is taken from Montenbruck, pp. 179ff.,
* with the exception, that we consider refraction for the maxima of
* partial and noncentral eclipses.
* Geographical positions are referred to sea level / the mean ellipsoid.
*
* Errors:
* - from uncertainty of JPL-ephemerides (0.01 arcsec):
* about 40 meters
* - from displacement of shadow points by atmospheric refraction:
* a few meters
* - from deviation of the geoid from the ellipsoid
* a few meters
* - from polar motion
* a few meters
* For geographical locations that are interesting for observation,
* the error is always < 100 m.
* However, if the sun is close to the horizon,
* all of these errors can grow up to a km or more.
*
* Function returns:
* -1 (ERR) on error (e.g. if swe_calc() for sun or moon fails)
* 0 if there is no solar eclipse at tjd
* SE_ECL_TOTAL
* SE_ECL_ANNULAR
* SE_ECL_TOTAL | SE_ECL_CENTRAL
* SE_ECL_TOTAL | SE_ECL_NONCENTRAL
* SE_ECL_ANNULAR | SE_ECL_CENTRAL
* SE_ECL_ANNULAR | SE_ECL_NONCENTRAL
* SE_ECL_PARTIAL
*
* geopos[0]: geographic longitude of central line
* geopos[1]: geographic latitude of central line
*
* not implemented so far:
*
* geopos[2]: geographic longitude of northern limit of umbra
* geopos[3]: geographic latitude of northern limit of umbra
* geopos[4]: geographic longitude of southern limit of umbra
* geopos[5]: geographic latitude of southern limit of umbra
* geopos[6]: geographic longitude of northern limit of penumbra
* geopos[7]: geographic latitude of northern limit of penumbra
* geopos[8]: geographic longitude of southern limit of penumbra
* geopos[9]: geographic latitude of southern limit of penumbra
*
* Attention: "northern" and "southern" limits of umbra do not
* necessarily correspond to the northernmost or southernmost
* geographic position, where the total, annular, or partial
* phase is visible at a given time.
* Imagine a situation in northern summer, when the sun illuminates
* the northern polar circle. The southernmost point of the core
* shadow may then touch the north pole, and therefore the
* northernmost point will be more in the south.
* Note also that with annular eclipses, the northern edge is
* usually geographically the southern one. With annular-total
* ones, the two lines cross, usually twice. The maximum is always
* total in such cases.
*
* attr[0] fraction of solar diameter covered by moon (magnitude)
* attr[1] ratio of lunar diameter to solar one
* attr[2] fraction of solar disc covered by moon (obscuration)
* attr[3] diameter of core shadow in km
* attr[4] azimuth of sun at tjd
* attr[5] true altitude of sun above horizon at tjd
* attr[6] apparent altitude of sun above horizon at tjd
* attr[7] angular distance of moon from sun in degrees
* declare as attr[20] at least !
*/
int32 CALL_CONV swe_sol_eclipse_where(
double tjd_ut,
int32 ifl,
double *geopos,
double *attr,
char *serr)
{
int32 retflag, retflag2;
double dcore[10];
ifl &= SEFLG_EPHMASK;
swi_set_tid_acc(tjd_ut, ifl, 0, serr);
if ((retflag = eclipse_where(tjd_ut, SE_SUN, NULL, ifl, geopos, dcore, serr)) < 0)
return retflag;
if ((retflag2 = eclipse_how(tjd_ut, SE_SUN, NULL, ifl, geopos[0], geopos[1], 0, attr, serr)) == ERR)
return retflag2;
attr[3] = dcore[0];
return retflag;
}
int32 CALL_CONV swe_lun_occult_where(
double tjd_ut,
int32 ipl,
char *starname,
int32 ifl,
double *geopos,
double *attr,
char *serr)
{
int32 retflag, retflag2;
double dcore[10];
if (ipl < 0) ipl = 0;
ifl &= SEFLG_EPHMASK;
swi_set_tid_acc(tjd_ut, ifl, 0, serr);
/* function calls for Pluto with asteroid number 134340
* are treated as calls for Pluto as main body SE_PLUTO */
if (ipl == SE_AST_OFFSET + 134340)
ipl = SE_PLUTO;
if ((retflag = eclipse_where(tjd_ut, ipl, starname, ifl, geopos, dcore, serr)) < 0)
return retflag;
if ((retflag2 = eclipse_how(tjd_ut, ipl, starname, ifl, geopos[0], geopos[1], 0, attr, serr)) == ERR)
return retflag2;
attr[3] = dcore[0];
return retflag;
}
/* Used by several swe_sol_eclipse_ functions.
* Like swe_sol_eclipse_where(), but instead of attr[0], it returns:
*
* dcore[0]: core shadow width in km
* dcore[2]: distance of shadow axis from geocenter r0
* dcore[3]: diameter of core shadow on fundamental plane d0
* dcore[4]: diameter of half-shadow on fundamental plane D0
*/
static int32 eclipse_where( double tjd_ut, int32 ipl, char *starname, int32 ifl, double *geopos, double *dcore,
char *serr)
{
int i;
int32 retc = 0, niter = 0;
double e[6], et[6], rm[6], rs[6], rmt[6], rst[6], xs[6], xst[6];
#if 0
double erm[6];
#endif
double x[6];
double lm[6], ls[6], lx[6];
double dsm, dsmt, d0, D0, s0, r0, d, s, dm;
double de = 6378140.0 / AUNIT;
double earthobl = 1 - EARTH_OBLATENESS;
double deltat, tjd, sidt;
double drad;
double sinf1, sinf2, cosf1, cosf2;
double rmoon = RMOON;
double dmoon = 2 * rmoon;
int32 iflag, iflag2;
/* double ecce = sqrt(2 * EARTH_OBLATENESS - EARTH_OBLATENESS * EARTH_OBLATENESS); */
AS_BOOL no_eclipse = FALSE;
struct epsilon *oe = &swed.oec;
for (i = 0; i < 10; i++)
dcore[i] = 0;
/* nutation need not be in lunar and solar positions,
* if mean sidereal time will be used */
iflag = SEFLG_SPEED | SEFLG_EQUATORIAL | ifl;
iflag2 = iflag | SEFLG_RADIANS;
iflag = iflag | SEFLG_XYZ;
deltat = swe_deltat_ex(tjd_ut, ifl, serr);
tjd = tjd_ut + deltat;
/* moon in cartesian coordinates */
if ((retc = swe_calc(tjd, SE_MOON, iflag, rm, serr)) == ERR)
return retc;
/* moon in polar coordinates */
if ((retc = swe_calc(tjd, SE_MOON, iflag2, lm, serr)) == ERR)
return retc;
/* sun in cartesian coordinates */
if ((retc = calc_planet_star(tjd, ipl, starname, iflag, rs, serr)) == ERR)
return retc;
/* sun in polar coordinates */
if ((retc = calc_planet_star(tjd, ipl, starname, iflag2, ls, serr)) == ERR)
return retc;
/* save sun position */
for (i = 0; i <= 2; i++)
rst[i] = rs[i];
/* save moon position */
for (i = 0; i <= 2; i++)
rmt[i] = rm[i];
if (iflag & SEFLG_NONUT)
sidt = swe_sidtime0(tjd_ut, oe->eps * RADTODEG, 0) * 15 * DEGTORAD;
else
sidt = swe_sidtime(tjd_ut) * 15 * DEGTORAD;
/*
* radius of planet disk in AU
*/
if (starname != NULL && *starname != '\0')
drad = 0;
else if (ipl < NDIAM)
drad = pla_diam[ipl] / 2 / AUNIT;
else if (ipl > SE_AST_OFFSET)
drad = swed.ast_diam / 2 * 1000 / AUNIT; /* km -> m -> AU */
else
drad = 0;
iter_where:
for (i = 0; i <= 2; i++) {
rs[i] = rst[i];
rm[i] = rmt[i];
}
/* Account for oblateness of earth:
* Instead of flattening the earth, we apply the
* correction to the z coordinate of the moon and
* the sun. This makes the calculation easier.
*/
for (i = 0; i <= 2; i++)
lx[i] = lm[i];
swi_polcart(lx, rm);
rm[2] /= earthobl;
/* distance of moon from geocenter */
dm = sqrt(square_sum(rm));
/* Account for oblateness of earth */
for (i = 0; i <= 2; i++)
lx[i] = ls[i];
swi_polcart(lx, rs);
rs[2] /= earthobl;
/* sun - moon vector */
for (i = 0; i <= 2; i++) {
e[i] = (rm[i] - rs[i]);
et[i] = (rmt[i] - rst[i]);
}
/* distance sun - moon */
dsm = sqrt(square_sum(e));
dsmt = sqrt(square_sum(et));
/* sun - moon unit vector */
for (i = 0; i <= 2; i++) {
e[i] /= dsm;
et[i] /= dsmt;
#if 0
erm[i] = rm[i] / dm;
#endif
}
sinf1 = ((drad - rmoon) / dsm);
cosf1 = sqrt(1 - sinf1 * sinf1);
sinf2 = ((drad + rmoon) / dsm);
cosf2 = sqrt(1 - sinf2 * sinf2);
/* distance of moon from fundamental plane */
s0 = -dot_prod(rm, e);
/* distance of shadow axis from geocenter */
r0 = sqrt(dm * dm - s0 * s0);
/* diameter of core shadow on fundamental plane */
d0 = (s0 / dsm * (drad * 2 - dmoon) - dmoon) / cosf1;
/* diameter of half-shadow on fundamental plane */
D0 = (s0 / dsm * (drad * 2 + dmoon) + dmoon) / cosf2;
dcore[2] = r0;
dcore[3] = d0;
dcore[4] = D0;
dcore[5] = cosf1;
dcore[6] = cosf2;
for (i = 2; i < 5; i++)
dcore[i] *= AUNIT / 1000.0;
/**************************
* central (total or annular) phase
**************************/
retc = 0;
if (de * cosf1 >= r0) {
retc |= SE_ECL_CENTRAL;
} else if (r0 <= de * cosf1 + fabs(d0) / 2) {
retc |= SE_ECL_NONCENTRAL;
} else if (r0 <= de * cosf2 + D0 / 2) {
retc |= (SE_ECL_PARTIAL | SE_ECL_NONCENTRAL);
} else {
if (serr != NULL)
sprintf(serr, "no solar eclipse at tjd = %f", tjd);
for (i = 0; i < 10; i++)
geopos[i] = 0;
*dcore = 0;
retc = 0;
d = 0;
no_eclipse = TRUE;
/*return retc;*/
}
/* distance of shadow point from fundamental plane */
d = s0 * s0 + de * de - dm * dm;
if (d > 0)
d = sqrt(d);
else
d = 0;
/* distance of moon from shadow point on earth */
s = s0 - d;
/* next: geographic position of eclipse center.
* if shadow axis does not touch the earth,
* place on earth with maximum occultation is computed.
*/
#if 0 /* the following stuff is meaningless for observations */
/*
* account for refraction at horizon
*/
if (d == 0) {
double ds, a, b;
/* distance of sun from geocenter */
ds = sqrt(square_sum(rs));
a = PI - acos(swi_dot_prod_unit(e, erm));
/* refraction at horizon + sun radius = about 0.83 degrees */
b = 34.4556 / 60.0 * DEGTORAD + asin(drad / ds);
# if 0
/* at edge of umbra and penumbra
* light rays are not parallel to shadow axis.
* for a short time close to contact of umbra and
* penumbra, an angle < 0.27 degrees would have
* to be subtracted from b;
*/
if (retc & SE_ECL_PARTIAL) {
d = d0;
sinf = sinf1;
} else {
d = D0;
sinf = sinf2;
}
c = (r0 - de) / d * 2 * sinf;
if (c > sinf1) {
b -= .....;
}
printf("%f %f %f", a * RADTODEG, b * RADTODEG, s);
printf(" %f\n", s);
# else
if (retc & SE_ECL_PARTIAL)
b -= asin(sinf2); /* maximum! */
else
b -= asin(sinf1);
# endif
s += tan(b) * cos(PI / 2 - a) * dm;
}
#endif
/* geographic position of eclipse center (maximum) */
for (i = 0; i <= 2; i++)
xs[i] = rm[i] + s * e[i];
/* we need geographic position with correct z, as well */
for (i = 0; i <= 2; i++)
xst[i] = xs[i];
xst[2] *= earthobl;
swi_cartpol(xst, xst);
if (niter <= 0) {
double cosfi = cos(xst[1]);
double sinfi = sin(xst[1]);
double eobl = EARTH_OBLATENESS;
double cc= 1 / sqrt(cosfi * cosfi + (1-eobl) * (1-eobl) * sinfi * sinfi);
double ss= (1-eobl) * (1-eobl) * cc;
earthobl = ss;
niter++;
goto iter_where;
}
swi_polcart(xst, xst);
/* to longitude and latitude */
swi_cartpol(xs, xs);
/* measure from sidereal time at greenwich */
xs[0] -= sidt;
xs[0] *= RADTODEG;
xs[1] *= RADTODEG;
xs[0] = swe_degnorm(xs[0]);
/* west is negative */
if (xs[0] > 180)
xs[0] -= 360;
geopos[0] = xs[0];
geopos[1] = xs[1];
/* diameter of core shadow:
* first, distance moon - place of eclipse on earth */
for (i = 0; i <= 2; i++)
x[i] = rmt[i] - xst[i];
s = sqrt(square_sum(x));
/* diameter of core shadow at place of maximum eclipse */
*dcore = (s / dsmt * ( drad * 2 - dmoon) - dmoon) * cosf1;
*dcore *= AUNIT / 1000.0;
/* diameter of penumbra at place of maximum eclipse */
dcore[1] = (s / dsmt * ( drad * 2 + dmoon) + dmoon) * cosf2;
dcore[1] *= AUNIT / 1000.0;
if (!(retc & SE_ECL_PARTIAL) && !no_eclipse) {
if (*dcore > 0) {
/*printf("annular\n");*/
retc |= SE_ECL_ANNULAR;
} else {
/*printf("total\n");*/
retc |= SE_ECL_TOTAL;
}
}
return retc;
}
static int32 calc_planet_star(double tjd_et, int32 ipl, char *starname, int32 iflag, double *x, char *serr)
{
int i;
int retc = OK;
if (starname == NULL || *starname == '\0') {
retc = swe_calc(tjd_et, ipl, iflag, x, serr);
} else {
if ((retc = swe_fixstar(starname, tjd_et, iflag, x, serr)) == OK) {
/* fixstars have the standard distance 1.
* in the occultation routines, this might lead to errors
* if interpreted as AU distance. To avoid this, we make it very high.
*/
if (iflag & SEFLG_XYZ) {
for (i = 0; i < 3; i++)
x[i] *= 100000000;
} else {
x[2] *= 100000000;
}
}
}
return retc;
}
/* Computes attributes of a solar eclipse for given tjd, geo. longitude,
* geo. latitude, and geo. height.
*
* retflag SE_ECL_TOTAL or SE_ECL_ANNULAR or SE_ECL_PARTIAL
* SE_ECL_NONCENTRAL
* if 0, no eclipse is visible at geogr. position.
*
* attr[0] fraction of solar diameter covered by moon;
* with total/annular eclipses, it results in magnitude acc. to IMCCE.
* attr[1] ratio of lunar diameter to solar one
* attr[2] fraction of solar disc covered by moon (obscuration)
* attr[3] diameter of core shadow in km
* attr[4] azimuth of sun at tjd
* attr[5] true altitude of sun above horizon at tjd
* attr[6] apparent altitude of sun above horizon at tjd
* attr[7] elongation of moon in degrees
* attr[8] magnitude acc. to NASA;
* = attr[0] for partial and attr[1] for annular and total eclipses
* attr[9] saros series number
* attr[10] saros series member number
* declare as attr[20] at least !
*
*/
int32 CALL_CONV swe_sol_eclipse_how(
double tjd_ut,
int32 ifl,
double *geopos,
double *attr,
char *serr)
{
int32 retflag, retflag2;
double dcore[10], ls[6], xaz[6];
double geopos2[20];
if (geopos[2] < SEI_ECL_GEOALT_MIN || geopos[2] > SEI_ECL_GEOALT_MAX) {
if (serr != NULL)
sprintf(serr, "location for eclipses must be between %.0f and %.0f m above sea", SEI_ECL_GEOALT_MIN, SEI_ECL_GEOALT_MAX);
return ERR;
}
ifl &= SEFLG_EPHMASK;
swi_set_tid_acc(tjd_ut, ifl, 0, serr);
if ((retflag = eclipse_how(tjd_ut, SE_SUN, NULL, ifl, geopos[0], geopos[1], geopos[2], attr, serr)) == ERR)
return retflag;
if ((retflag2 = eclipse_where(tjd_ut, SE_SUN, NULL, ifl, geopos2, dcore, serr)) == ERR)
return retflag2;
if (retflag)
retflag |= (retflag2 & (SE_ECL_CENTRAL | SE_ECL_NONCENTRAL));
attr[3] = dcore[0];
swe_set_topo(geopos[0], geopos[1], geopos[2]);
if (swe_calc_ut(tjd_ut, SE_SUN, ifl | SEFLG_TOPOCTR | SEFLG_EQUATORIAL, ls, serr) == ERR)
return ERR;
swe_azalt(tjd_ut, SE_EQU2HOR, geopos, 0, 10, ls, xaz);
attr[4] = xaz[0];
attr[5] = xaz[1];
attr[6] = xaz[2];
if (xaz[2] <= 0)
retflag = 0;
return retflag;
}
#define USE_AZ_NAV 0
static int32 eclipse_how( double tjd_ut, int32 ipl, char *starname, int32 ifl,
double geolon, double geolat, double geohgt,
double *attr, char *serr)
{
int i, j, k;
int32 retc = 0;
double te, d;
double xs[6], xm[6], ls[6], lm[6], x1[6], x2[6];
double rmoon, rsun, rsplusrm, rsminusrm;
double dctr;
double drad;
int32 iflag = SEFLG_EQUATORIAL | SEFLG_TOPOCTR | ifl;
int32 iflagcart = iflag | SEFLG_XYZ;
#if USE_AZ_NAV
double mdd, eps, sidt, armc;
#endif
double xh[6], hmin_appr;
double lsun, lmoon, lctr, lsunleft, a, b, sc1, sc2;
double geopos[3];
for (i = 0; i < 10; i++)
attr[i] = 0;
geopos[0] = geolon;
geopos[1] = geolat;
geopos[2] = geohgt;
te = tjd_ut + swe_deltat_ex(tjd_ut, ifl, serr);
swe_set_topo(geolon, geolat, geohgt);
if (calc_planet_star(te, ipl, starname, iflag, ls, serr) == ERR)
return ERR;
if (swe_calc(te, SE_MOON, iflag, lm, serr) == ERR)
return ERR;
if (calc_planet_star(te, ipl, starname, iflagcart, xs, serr) == ERR)
return ERR;
if (swe_calc(te, SE_MOON, iflagcart, xm, serr) == ERR)
return ERR;
/*
* radius of planet disk in AU
*/
if (starname != NULL && *starname != '\0')
drad = 0;
else if (ipl < NDIAM)
drad = pla_diam[ipl] / 2 / AUNIT;
else if (ipl > SE_AST_OFFSET)
drad = swed.ast_diam / 2 * 1000 / AUNIT; /* km -> m -> AU */
else
drad = 0;
/*
* azimuth and altitude of sun or planet
*/
#if USE_AZ_NAV /* old */
eps = swi_epsiln(te, iflag);
if (iflag & SEFLG_NONUT)
sidt = swe_sidtime0(tjd_ut, eps * RADTODEG, 0) * 15;
else
sidt = swe_sidtime(tjd_ut) * 15;
armc = sidt + geolon;
mdd = swe_degnorm(ls[0] - armc);
xh[0] = swe_degnorm(mdd - 90);
xh[1] = ls[1];
xh[2] = ls[2];
swe_cotrans(xh, xh, 90 - geolat); /* azimuth from east, counterclock, via north */
#else
swe_azalt(tjd_ut, SE_EQU2HOR, geopos, 0, 10, ls, xh); /* azimuth from south, clockwise, via west */
#endif
/* eclipse description */
rmoon = asin(RMOON / lm[2]) * RADTODEG;
rsun = asin(drad / ls[2]) * RADTODEG;
rsplusrm = rsun + rmoon;
rsminusrm = rsun - rmoon;
for (i = 0; i < 3; i++) {
x1[i] = xs[i] / ls[2];
x2[i] = xm[i] / lm[2];
}
dctr = acos(swi_dot_prod_unit(x1, x2)) * RADTODEG;
/*
* phase
*/
if (dctr < rsminusrm)
retc = SE_ECL_ANNULAR;
else if (dctr < fabs(rsminusrm))
retc = SE_ECL_TOTAL;
else if (dctr < rsplusrm)
retc = SE_ECL_PARTIAL;
else {
retc = 0;
if (serr != NULL)
sprintf(serr, "no solar eclipse at tjd = %f", tjd_ut);
}
/*
* ratio of diameter of moon to that of sun
*/
if (rsun > 0)
attr[1] = rmoon / rsun;
else
attr[1] = 0;
/*
* eclipse magnitude:
* fraction of solar diameter covered by moon
*/
lsun = asin(rsun / 2 * DEGTORAD) * 2;
lsunleft = (-dctr + rsun + rmoon);
if (lsun > 0) {
attr[0] = lsunleft / rsun / 2;
} else {
attr[0] = 100;
}
/*if (retc == SE_ECL_ANNULAR || retc == SE_ECL_TOTAL)
attr[0] = attr[1];*/
/*
* obscuration:
* fraction of solar disc obscured by moon
*/
lsun = rsun;
lmoon = rmoon;
lctr = dctr;
if (retc == 0 || lsun == 0) {
attr[2] = 100;
} else if (retc == SE_ECL_TOTAL || retc == SE_ECL_ANNULAR) {
attr[2] = lmoon * lmoon / lsun / lsun;
} else {
a = 2 * lctr * lmoon;
b = 2 * lctr * lsun;
if (a < 1e-9) {
attr[2] = lmoon * lmoon / lsun / lsun;
} else {
a = (lctr * lctr + lmoon * lmoon - lsun * lsun) / a;
if (a > 1) a = 1;
if (a < -1) a = -1;
b = (lctr * lctr + lsun * lsun - lmoon * lmoon) / b;
if (b > 1) b = 1;
if (b < -1) b = -1;
a = acos(a);
b = acos(b);
sc1 = a * lmoon * lmoon / 2;
sc2 = b * lsun * lsun / 2;
sc1 -= (cos(a) * sin(a)) * lmoon * lmoon / 2;
sc2 -= (cos(b) * sin(b)) * lsun * lsun / 2;
attr[2] = (sc1 + sc2) * 2 / PI / lsun / lsun;
}
}
attr[7] = dctr;
/* approximate minimum height for visibility, considering
* refraction and dip
* 34.4556': refraction at horizon, from Bennets formulae
* 1.75' / sqrt(geohgt): dip of horizon
* 0.37' / sqrt(geohgt): refraction between horizon and observer */
hmin_appr = -(34.4556 + (1.75 + 0.37) * sqrt(geohgt)) / 60;
if (xh[1] + rsun + fabs(hmin_appr) >= 0 && retc)
retc |= SE_ECL_VISIBLE; /* eclipse visible */
#if USE_AZ_NAV /* old */
attr[4] = swe_degnorm(90 - xh[0]); /* azimuth, from north, clockwise, via east */
#else
attr[4] = xh[0]; /* azimuth, from south, clockwise, via west */
#endif
attr[5] = xh[1]; /* height */
attr[6] = xh[2]; /* height */
if (ipl == SE_SUN && (starname == NULL || *starname == '\0')) {
/* magnitude of solar eclipse according to NASA */
attr[8] = attr[0]; /* fraction of diameter occulted */
if (retc & (SE_ECL_TOTAL | SE_ECL_ANNULAR))
attr[8] = attr[1]; /* ratio between diameters of sun and moon */
/* saros series and member */
for (i = 0; i < NSAROS_SOLAR; i++) {
d = (tjd_ut - saros_data_solar[i].tstart) / SAROS_CYCLE;
if (d < 0) continue;
j = (int) d;
if ((d - j) * SAROS_CYCLE < 2) {
attr[9] = (double) saros_data_solar[i].series_no;
attr[10] = (double) j + 1;
break;
}
k = j + 1;
if ((k - d) * SAROS_CYCLE < 2) {
attr[9] = (double) saros_data_solar[i].series_no;
attr[10] = (double) k + 1;
break;
}
}
if (i == NSAROS_SOLAR) {
attr[9] = attr[10] = -99999999;
}
}
return retc;
}
/* When is the next solar eclipse anywhere on earth?
*
* input parameters:
*
* tjd_start start time for search (UT)
* ifl ephemeris to be used (SEFLG_SWIEPH, etc.)
* ifltype eclipse type to be searched (SE_ECL_TOTAL, etc.)
* 0, if any type of eclipse is wanted
*
* return values:
*
* retflag SE_ECL_TOTAL or SE_ECL_ANNULAR or SE_ECL_PARTIAL
* or SE_ECL_ANNULAR_TOTAL
* SE_ECL_CENTRAL
* SE_ECL_NONCENTRAL
*
* tret[0] time of maximum eclipse
* tret[1] time, when eclipse takes place at local apparent noon
* tret[2] time of eclipse begin
* tret[3] time of eclipse end
* tret[4] time of totality begin
* tret[5] time of totality end
* tret[6] time of center line begin
* tret[7] time of center line end
* tret[8] time when annular-total eclipse becomes total
* not implemented so far
* tret[9] time when annular-total eclipse becomes annular again
* not implemented so far
* declare as tret[10] at least!
*
*/
int32 CALL_CONV swe_sol_eclipse_when_glob(double tjd_start, int32 ifl, int32 ifltype,
double *tret, int32 backward, char *serr)
{
int i, j, k, m, n, o, i1 = 0, i2 = 0;
int32 retflag = 0, retflag2 = 0;
double de = 6378.140, a;
double t, tt, tjd, tjds, dt, dtint, dta, dtb;
double T, T2, T3, T4, K, M, Mm;
double E, Ff;
double xs[6], xm[6], ls[6], lm[6];
double rmoon, rsun, dcore[10];
double dc[3], dctr;
double twohr = 2.0 / 24.0;
double tenmin = 10.0 / 24.0 / 60.0;
double dt1 = 0, dt2 = 0;
double geopos[20], attr[20];
double dtstart, dtdiv;
double xa[6], xb[6];
int direction = 1;
AS_BOOL dont_times = FALSE;
int32 iflag, iflagcart;
ifl &= SEFLG_EPHMASK;
swi_set_tid_acc(tjd_start, ifl, 0, serr);
iflag = SEFLG_EQUATORIAL | ifl;
iflagcart = iflag | SEFLG_XYZ;
if (ifltype == (SE_ECL_PARTIAL | SE_ECL_CENTRAL)) {
if (serr != NULL)
strcpy(serr, "central partial eclipses do not exist");
return ERR;
}
if (ifltype == (SE_ECL_ANNULAR_TOTAL | SE_ECL_NONCENTRAL)) {
if (serr != NULL)
strcpy(serr, "non-central hybrid (annular-total) eclipses do not exist");
return ERR;
}
if (ifltype == 0)
ifltype = SE_ECL_TOTAL | SE_ECL_ANNULAR | SE_ECL_PARTIAL
| SE_ECL_ANNULAR_TOTAL | SE_ECL_NONCENTRAL | SE_ECL_CENTRAL;
if (ifltype == SE_ECL_TOTAL || ifltype == SE_ECL_ANNULAR || ifltype == SE_ECL_ANNULAR_TOTAL)
ifltype |= (SE_ECL_NONCENTRAL | SE_ECL_CENTRAL);
if (ifltype == SE_ECL_PARTIAL)
ifltype |= SE_ECL_NONCENTRAL;
if (backward)
direction = -1;
K = (int) ((tjd_start - J2000) / 365.2425 * 12.3685);
K -= direction;
next_try:
retflag = 0;
dont_times = FALSE;
for (i = 0; i <= 9; i++)
tret[i] = 0;
T = K / 1236.85;
T2 = T * T; T3 = T2 * T; T4 = T3 * T;
Ff = swe_degnorm(160.7108 + 390.67050274 * K
- 0.0016341 * T2
- 0.00000227 * T3
+ 0.000000011 * T4);
if (Ff > 180)
Ff -= 180;
if (Ff > 21 && Ff < 159) { /* no eclipse possible */
K += direction;
goto next_try;
}
/* approximate time of geocentric maximum eclipse
* formula from Meeus, German, p. 381 */
tjd = 2451550.09765 + 29.530588853 * K
+ 0.0001337 * T2
- 0.000000150 * T3
+ 0.00000000073 * T4;
M = swe_degnorm(2.5534 + 29.10535669 * K
- 0.0000218 * T2
- 0.00000011 * T3);
Mm = swe_degnorm(201.5643 + 385.81693528 * K
+ 0.1017438 * T2
+ 0.00001239 * T3
+ 0.000000058 * T4);
E = 1 - 0.002516 * T - 0.0000074 * T2;
M *= DEGTORAD;
Mm *= DEGTORAD;
tjd = tjd - 0.4075 * sin(Mm)
+ 0.1721 * E * sin(M);
/*
* time of maximum eclipse (if eclipse) =
* minimum geocentric angle between sun and moon edges.
* After this time has been determined, check
* whether or not an eclipse is taking place with
* the functions eclipse_where() and _how().
*/
dtstart = 1;
if (tjd < 2000000 || tjd > 2500000)
dtstart = 5;
dtdiv = 4;
for (dt = dtstart;
dt > 0.0001;
dt /= dtdiv) {
for (i = 0, t = tjd - dt; i <= 2; i++, t += dt) {
if (swe_calc(t, SE_SUN, iflag, ls, serr) == ERR)
return ERR;
if (swe_calc(t, SE_MOON, iflag, lm, serr) == ERR)
return ERR;
if (swe_calc(t, SE_SUN, iflagcart, xs, serr) == ERR)
return ERR;
if (swe_calc(t, SE_MOON, iflagcart, xm, serr) == ERR)
return ERR;
for (m = 0; m < 3; m++) {
xa[m] = xs[m] / ls[2];
xb[m] = xm[m] / lm[2];
}
dc[i] = acos(swi_dot_prod_unit(xa, xb)) * RADTODEG;
rmoon = asin(RMOON / lm[2]) * RADTODEG;
rsun = asin(RSUN / ls[2]) * RADTODEG;
dc[i] -= (rmoon + rsun);
}
find_maximum(dc[0], dc[1], dc[2], dt, &dtint, &dctr);
tjd += dtint + dt;
}
tjds = tjd - swe_deltat_ex(tjd, ifl, serr);
tjds = tjd - swe_deltat_ex(tjds, ifl, serr);
tjds = tjd = tjd - swe_deltat_ex(tjds, ifl, serr);
if ((retflag = eclipse_where(tjd, SE_SUN, NULL, ifl, geopos, dcore, serr)) == ERR)
return retflag;
retflag2 = retflag;
/* in extreme cases _where() returns no eclipse, where there is
* actually a very small one, therefore call _how() with the
* coordinates returned by _where(): */
if ((retflag2 = eclipse_how(tjd, SE_SUN, NULL, ifl, geopos[0], geopos[1], 0, attr, serr)) == ERR)
return retflag2;
if (retflag2 == 0) {
K += direction;
goto next_try;
}
tret[0] = tjd;
if ((backward && tret[0] >= tjd_start - 0.0001)
|| (!backward && tret[0] <= tjd_start + 0.0001)) {
K += direction;
goto next_try;
}
/*
* eclipse type, SE_ECL_TOTAL, _ANNULAR, etc.
* SE_ECL_ANNULAR_TOTAL will be discovered later
*/
if ((retflag = eclipse_where(tjd, SE_SUN, NULL, ifl, geopos, dcore, serr)) == ERR)
return retflag;
if (retflag == 0) { /* can happen with extremely small percentage */
retflag = SE_ECL_PARTIAL | SE_ECL_NONCENTRAL;
tret[4] = tret[5] = tjd; /* fix this ???? */
dont_times = TRUE;
}
/*
* check whether or not eclipse type found is wanted
*/
/* non central eclipse is wanted: */
if (!(ifltype & SE_ECL_NONCENTRAL) && (retflag & SE_ECL_NONCENTRAL)) {
K += direction;
goto next_try;
}
/* central eclipse is wanted: */
if (!(ifltype & SE_ECL_CENTRAL) && (retflag & SE_ECL_CENTRAL)) {
K += direction;
goto next_try;
}
/* non annular eclipse is wanted: */
if (!(ifltype & SE_ECL_ANNULAR) && (retflag & SE_ECL_ANNULAR)) {
K += direction;
goto next_try;
}
/* non partial eclipse is wanted: */
if (!(ifltype & SE_ECL_PARTIAL) && (retflag & SE_ECL_PARTIAL)) {
K += direction;
goto next_try;
}
/* annular-total eclipse will be discovered later */
if (!(ifltype & (SE_ECL_TOTAL | SE_ECL_ANNULAR_TOTAL)) && (retflag & SE_ECL_TOTAL)) {
K += direction;
goto next_try;
}
if (dont_times)
goto end_search_global;
/*
* n = 0: times of eclipse begin and end
* n = 1: times of totality begin and end
* n = 2: times of center line begin and end
*/
if (retflag & SE_ECL_PARTIAL)
o = 0;
else if (retflag & SE_ECL_NONCENTRAL)
o = 1;
else
o = 2;
dta = twohr;
dtb = tenmin / 3.0;
for (n = 0; n <= o; n++) {
if (n == 0) {
/*dc[1] = dcore[3] / 2 + de - dcore[1];*/
i1 = 2; i2 = 3;
} else if (n == 1) {
if (retflag & SE_ECL_PARTIAL)
continue;
i1 = 4; i2 = 5;
} else if (n == 2) {
if (retflag & SE_ECL_NONCENTRAL)
continue;
i1 = 6; i2 = 7;
}
for (i = 0, t = tjd - dta; i <= 2; i += 1, t += dta) {
if ((retflag2 = eclipse_where(t, SE_SUN, NULL, ifl, geopos, dcore, serr)) == ERR)
return retflag2;
if (n == 0)
dc[i] = dcore[4] / 2 + de / dcore[5] - dcore[2];
else if (n == 1)
dc[i] = fabs(dcore[3]) / 2 + de / dcore[6] - dcore[2];
else if (n == 2)
dc[i] = de / dcore[6] - dcore[2];
}
find_zero(dc[0], dc[1], dc[2], dta, &dt1, &dt2);
tret[i1] = tjd + dt1 + dta;
tret[i2] = tjd + dt2 + dta;
for (m = 0, dt = dtb; m < 3; m++, dt /= 3) {
for (j = i1; j <= i2; j += (i2 - i1)) {
for (i = 0, t = tret[j] - dt; i < 2; i++, t += dt) {
if ((retflag2 = eclipse_where(t, SE_SUN, NULL, ifl, geopos, dcore, serr)) == ERR)
return retflag2;
if (n == 0)
dc[i] = dcore[4] / 2 + de / dcore[5] - dcore[2];
else if (n == 1)
dc[i] = fabs(dcore[3]) / 2 + de / dcore[6] - dcore[2];
else if (n == 2)
dc[i] = de / dcore[6] - dcore[2];
}
dt1 = dc[1] / ((dc[1] - dc[0]) / dt);
tret[j] -= dt1;
}
}
}
/*
* annular-total eclipses
*/
if (retflag & SE_ECL_TOTAL) {
if ((retflag2 = eclipse_where(tret[0], SE_SUN, NULL, ifl, geopos, dcore, serr)) == ERR)
return retflag2;
dc[0] = *dcore;
if ((retflag2 = eclipse_where(tret[4], SE_SUN, NULL, ifl, geopos, dcore, serr)) == ERR)
return retflag2;
dc[1] = *dcore;
if ((retflag2 = eclipse_where(tret[5], SE_SUN, NULL, ifl, geopos, dcore, serr)) == ERR)
return retflag2;
dc[2] = *dcore;
/* the maximum is always total, and there is either one or
* to times before and after, when the core shadow becomes
* zero and totality changes into annularity or vice versa.
*/
if (dc[0] * dc[1] < 0 || dc[0] * dc[2] < 0) {
retflag |= SE_ECL_ANNULAR_TOTAL;
retflag &= ~SE_ECL_TOTAL;
}
}
/* if eclipse is given but not wanted: */
if (!(ifltype & SE_ECL_TOTAL) && (retflag & SE_ECL_TOTAL)) {
K += direction;
goto next_try;
}
/* if annular_total eclipse is given but not wanted: */
if (!(ifltype & SE_ECL_ANNULAR_TOTAL) && (retflag & SE_ECL_ANNULAR_TOTAL)) {
K += direction;
goto next_try;
}
/*
* time of maximum eclipse at local apparent noon
*/
/* first, find out, if there is a solar transit
* between begin and end of eclipse */
k = 2;
for (i = 0; i < 2; i++) {
j = i + k;
tt = tret[j] + swe_deltat_ex(tret[j], ifl, serr);
if (swe_calc(tt, SE_SUN, iflag, ls, serr) == ERR)
return ERR;
if (swe_calc(tt, SE_MOON, iflag, lm, serr) == ERR)
return ERR;
dc[i] = swe_degnorm(ls[0] - lm[0]);
if (dc[i] > 180)
dc[i] -= 360;
}
if (dc[0] * dc[1] >= 0) /* no transit */
tret[1] = 0;
else {
tjd = tjds;
dt = 0.1;
dt1 = (tret[3] - tret[2]) / 2.0;
if (dt1 < dt)
dt = dt1 / 2.0;
for (j = 0;
dt > 0.01;
j++, dt /= 3) {
for (i = 0, t = tjd; i <= 1; i++, t -= dt) {
tt = t + swe_deltat_ex(t, ifl, serr);
if (swe_calc(tt, SE_SUN, iflag, ls, serr) == ERR)
return ERR;
if (swe_calc(tt, SE_MOON, iflag, lm, serr) == ERR)
return ERR;
dc[i] = swe_degnorm(ls[0] - lm[0]);
if (dc[i] > 180)
dc[i] -= 360;
if (dc[i] > 180)
dc[i] -= 360;
}
a = (dc[1] - dc[0]) / dt;
if (a < 1e-10)
break;
dt1 = dc[0] / a;
tjd += dt1;
}
tret[1] = tjd;
}
end_search_global:
return retflag;
/*
* the time of maximum occultation is practically identical
* with the time of maximum core shadow diameter.
*
* the time, when duration of totality is maximal,
* is not an interesting computation either. Near the maximum
* occulation, the time of totality can be the same by
* a second for hundreds of kilometers (for 10 minutes
* or more).
*
* for annular eclipses the maximum duration is close to the
* beginning and the end of the center lines, where is also
* the minimum of core shadow diameter.
*/
}
/* When is the next lunar occultation anywhere on earth?
* This function also finds solar eclipses, but is less efficient
* than swe_sol_eclipse_when_glob().
*
* input parameters:
*
* tjd_start start time for search (UT)
* ipl planet number of occulted body
* starname name of occulted star. Must be NULL or "", if a planetary
* occultation is to be calculated. For the use of this
* field, also see swe_fixstar().
* ifl ephemeris to be used (SEFLG_SWIEPH, etc.)
* ephemeris flag.
*
* ifltype eclipse type to be searched (SE_ECL_TOTAL, etc.)
* 0, if any type of eclipse is wanted
* this functionality also works with occultations
*
* backward if 1, causes search backward in time
*
* If you want to have only one conjunction
* of the moon with the body tested, add the following flag:
* backward |= SE_ECL_ONE_TRY. If this flag is not set,
* the function will search for an occultation until it
* finds one. For bodies with ecliptical latitudes > 5,
* the function may search successlessly until it reaches
* the end of the ephemeris.
* (Note: we do not add SE_ECL_ONE_TRY to ifl, because
* ifl may contain SEFLG_TOPOCTR (=SE_ECL_ONE_TRY) from
* the parameter iflag of swe_calc() etc. Although the
* topocentric flag is irrelevant here, it might cause
* confusion.)
*
* return values:
*
* retflag SE_ECL_TOTAL or SE_ECL_ANNULAR or SE_ECL_PARTIAL
* or SE_ECL_ANNULAR_TOTAL
* SE_ECL_CENTRAL
* SE_ECL_NONCENTRAL
*
* tret[0] time of maximum eclipse
* tret[1] time, when eclipse takes place at local apparent noon
* tret[2] time of eclipse begin
* tret[3] time of eclipse end
* tret[4] time of totality begin
* tret[5] time of totality end
* tret[6] time of center line begin
* tret[7] time of center line end
* tret[8] time when annular-total eclipse becomes total
* not implemented so far
* tret[9] time when annular-total eclipse becomes annular again
* not implemented so far
* declare as tret[10] at least!
*
*/
int32 CALL_CONV swe_lun_occult_when_glob(
double tjd_start, int32 ipl, char *starname, int32 ifl, int32 ifltype,
double *tret, int32 backward, char *serr)
{
int i, j, k, m, n, o, i1, i2;
int32 retflag = 0, retflag2 = 0;
double de = 6378.140, a;
double t, tt, tjd = 0, tjds, dt, dtint, dta, dtb;
double drad, dl;
double xs[6], xm[6], ls[6], lm[6];
double rmoon, rsun, dcore[10];
double dc[20], dctr;
double twohr = 2.0 / 24.0;
double tenmin = 10.0 / 24.0 / 60.0;
double dt1 = 0, dt2 = 0, dadd2 = 1;
double geopos[20];
double dtstart, dtdiv;
int direction = 1;
int32 ifltype2;
int32 iflag, iflagcart;
AS_BOOL dont_times = FALSE;
int32 one_try = backward & SE_ECL_ONE_TRY;
if (ipl < 0) ipl = 0;
/*if (backward & SEI_OCC_FAST)
dont_times = TRUE; */
/* function calls for Pluto with asteroid number 134340
* are treated as calls for Pluto as main body SE_PLUTO */
if (ipl == SE_AST_OFFSET + 134340)
ipl = SE_PLUTO;
ifl &= SEFLG_EPHMASK;
swi_set_tid_acc(tjd_start, ifl, 0, serr);
iflag = SEFLG_EQUATORIAL | ifl;
iflagcart = iflag | SEFLG_XYZ;
backward &= 1L;
/*
* initializations
*/
if (ifltype == (SE_ECL_PARTIAL | SE_ECL_CENTRAL)) {
if (serr != NULL)
strcpy(serr, "central partial eclipses do not exist");
return ERR;
}
if (ipl != SE_SUN) {
ifltype2 = (ifltype & ~(SE_ECL_NONCENTRAL | SE_ECL_CENTRAL));
if (ifltype2 == SE_ECL_ANNULAR || ifltype == SE_ECL_ANNULAR_TOTAL) {
if (serr != NULL)
sprintf(serr, "annular occulation do not exist for object %d %s\n", ipl, starname);
return ERR;
}
}
if (ipl != SE_SUN && (ifltype & (SE_ECL_ANNULAR | SE_ECL_ANNULAR_TOTAL)))
ifltype &= ~(SE_ECL_ANNULAR|SE_ECL_ANNULAR_TOTAL);
if (ifltype == 0) {
ifltype = SE_ECL_TOTAL | SE_ECL_PARTIAL | SE_ECL_NONCENTRAL | SE_ECL_CENTRAL;
if (ipl == SE_SUN)
ifltype |= (SE_ECL_ANNULAR | SE_ECL_ANNULAR_TOTAL);
}
if (ifltype & (SE_ECL_TOTAL | SE_ECL_ANNULAR | SE_ECL_ANNULAR_TOTAL))
ifltype |= (SE_ECL_NONCENTRAL | SE_ECL_CENTRAL);
if (ifltype & SE_ECL_PARTIAL)
ifltype |= SE_ECL_NONCENTRAL;
retflag = 0;
for (i = 0; i <= 9; i++)
tret[i] = 0;
if (backward)
direction = -1;
t = tjd_start;
tjd = t;
next_try:
if (calc_planet_star(t, ipl, starname, ifl, ls, serr) == ERR)
return ERR;
/* fixed stars with an ecliptic latitude > 7 or < -7 cannot have
* an occultation. Even lunar parallax andd proper motion of star
* will never allow it. */
if (fabs(ls[1]) > 7 && starname != NULL && *starname != '\0') {
if (serr != NULL)
sprintf(serr, "occultation never occurs: star %s has ecl. lat. %.1f", starname, ls[1]);
return ERR;
}
if (swe_calc(t, SE_MOON, ifl, lm, serr) == ERR)
return ERR;
dl = swe_degnorm(ls[0] - lm[0]);
if (direction < 0)
dl -= 360;
/* get rough conjunction in ecliptic longitude */
while (fabs(dl) > 0.1) {
t += dl / 13;
if (calc_planet_star(t, ipl, starname, ifl, ls, serr) == ERR)
return ERR;
if (swe_calc(t, SE_MOON, ifl, lm, serr) == ERR)
return ERR;
dl = swe_degnorm(ls[0] - lm[0]);
if (dl > 180) dl -= 360;
}
tjd = t;
/* difference in latitude too big for an occultation */
drad = fabs(ls[1] - lm[1]);
if (drad > 2) {
if (one_try) {
tret[0] = t + direction; /* return a date suitable for next try */
return 0;
}
t += direction * 20;
tjd = t;
goto next_try;
}
/*
* radius of planet disk in AU
*/
if (starname != NULL && *starname != '\0')
drad = 0;
else if (ipl < NDIAM)
drad = pla_diam[ipl] / 2 / AUNIT;
else if (ipl > SE_AST_OFFSET)
drad = swed.ast_diam / 2 * 1000 / AUNIT; /* km -> m -> AU */
else
drad = 0;
/*
* time of maximum eclipse (if eclipse) =
* minimum geocentric angle between sun and moon edges.
* After this time has been determined, check
* whether or not an eclipse is taking place with
* the functions eclipse_where() and _how().
*/
dtstart = dadd2; /* originally 1 */
dtdiv = 3;
for (dt = dtstart;
dt > 0.0001;
dt /= dtdiv) {
for (i = 0, t = tjd - dt; i <= 2; i++, t += dt) {
if (calc_planet_star(t, ipl, starname, iflag, ls, serr) == ERR)
return ERR;
if (swe_calc(t, SE_MOON, iflag, lm, serr) == ERR)
return ERR;
if (calc_planet_star(t, ipl, starname, iflagcart, xs, serr) == ERR)
return ERR;
if (swe_calc(t, SE_MOON, iflagcart, xm, serr) == ERR)
return ERR;
dc[i] = acos(swi_dot_prod_unit(xs, xm)) * RADTODEG;
rmoon = asin(RMOON / lm[2]) * RADTODEG;
rsun = asin(drad / ls[2]) * RADTODEG;
dc[i] -= (rmoon + rsun);
}
find_maximum(dc[0], dc[1], dc[2], dt, &dtint, &dctr);
tjd += dtint + dt;
}
tjd -= swe_deltat_ex(tjd, ifl, serr);
tjds = tjd;
if ((retflag = eclipse_where(tjd, ipl, starname, ifl, geopos, dcore, serr)) == ERR)
return retflag;
retflag2 = retflag;
/* in extreme cases _where() returns no eclipse, where there is
* actually a very small one, therefore call _how() with the
* coordinates returned by _where(): */
/* if ((retflag2 = eclipse_how(tjd, ipl, starname, ifl, geopos[0], geopos[1], 0, attr, serr)) == ERR)
return retflag2; */
if (retflag2 == 0) {
/* only one try! */
/* if (one_try && ((direction == 1 && tjd > tjd_start) || (direction == -1 && tjd < tjd_start))) {*/
if (one_try) {
tret[0] = tjd;
return 0;
}
/*t= tjd + direction * dadd;*/
t = tjd + direction * 20;
tjd = t;
goto next_try;
}
tret[0] = tjd;
/* should not happen anymore Version 2.01 */
if ((backward && tret[0] >= tjd_start - 0.0001)
|| (!backward && tret[0] <= tjd_start + 0.0001)) {
/*t= tjd + direction * dadd;*/
t = tjd + direction * 20;
tjd = t;
goto next_try;
}
/*
* eclipse type, SE_ECL_TOTAL, _ANNULAR, etc.
* SE_ECL_ANNULAR_TOTAL will be discovered later
*/
if ((retflag = eclipse_where(tjd, ipl, starname, ifl, geopos, dcore, serr)) == ERR)
return retflag;
if (retflag == 0) { /* can happen with extremely small percentage */
retflag = SE_ECL_PARTIAL | SE_ECL_NONCENTRAL;
tret[4] = tret[5] = tjd; /* fix this ???? */
dont_times = TRUE;
}
/*
* check whether or not eclipse type found is wanted
*/
/* non central eclipse is wanted: */
if (!(ifltype & SE_ECL_NONCENTRAL) && (retflag & SE_ECL_NONCENTRAL)) {
/*t= tjd + direction * dadd;*/
t = tjd + direction * 20;
if (one_try) {
tret[0] = tjd;
return 0;
}
tjd = t;
goto next_try;
}
/* central eclipse is wanted: */
if (!(ifltype & SE_ECL_CENTRAL) && (retflag & SE_ECL_CENTRAL)) {
/*t= tjd + direction * dadd;*/
t = tjd + direction * 20;
if (one_try) {
tret[0] = tjd;
return 0;
}
tjd = t;
goto next_try;
}
/* non annular eclipse is wanted: */
if (!(ifltype & SE_ECL_ANNULAR) && (retflag & SE_ECL_ANNULAR)) {
/*t= tjd + direction * dadd;*/
t = tjd + direction * 20;
if (one_try) {
tret[0] = tjd;
return 0;
}
tjd = t;
goto next_try;
}
/* non partial eclipse is wanted: */
if (!(ifltype & SE_ECL_PARTIAL) && (retflag & SE_ECL_PARTIAL)) {
/*t= tjd + direction * dadd;*/
t = tjd + direction * 20;
if (one_try) {
tret[0] = tjd;
return 0;
}
tjd = t;
goto next_try;
}
/* annular-total eclipse will be discovered later */
if (!(ifltype & (SE_ECL_TOTAL | SE_ECL_ANNULAR_TOTAL)) && (retflag & SE_ECL_TOTAL)) {
/*t= tjd + direction * dadd;*/
t = tjd + direction * 20;
if (one_try) {
tret[0] = tjd;
return 0;
}
tjd = t;
goto next_try;
}
if (dont_times)
goto end_search_global;
/*
* n = 0: times of eclipse begin and end
* n = 1: times of totality begin and end
* n = 2: times of center line begin and end
*/
if (retflag & SE_ECL_PARTIAL)
o = 0;
else if (retflag & SE_ECL_NONCENTRAL)
o = 1;
else
o = 2;
dta = twohr;
dtb = tenmin;
for (n = 0; n <= o; n++) {
if (n == 0) {
/*dc[1] = dcore[3] / 2 + de - dcore[1];*/
i1 = 2; i2 = 3;
} else if (n == 1) {
if (retflag & SE_ECL_PARTIAL)
continue;
i1 = 4; i2 = 5;
} else if (n == 2) {
if (retflag & SE_ECL_NONCENTRAL)
continue;
i1 = 6; i2 = 7;
}
for (i = 0, t = tjd - dta; i <= 2; i += 1, t += dta) {
if ((retflag2 = eclipse_where(t, ipl, starname, ifl, geopos, dcore, serr)) == ERR)
return retflag2;
if (n == 0)
dc[i] = dcore[4] / 2 + de / dcore[5] - dcore[2];
else if (n == 1)
dc[i] = fabs(dcore[3]) / 2 + de / dcore[6] - dcore[2];
else if (n == 2)
dc[i] = de / dcore[6] - dcore[2];
}
find_zero(dc[0], dc[1], dc[2], dta, &dt1, &dt2);
tret[i1] = tjd + dt1 + dta;
tret[i2] = tjd + dt2 + dta;
for (m = 0, dt = dtb; m < 3; m++, dt /= 3) {
for (j = i1; j <= i2; j += (i2 - i1)) {
for (i = 0, t = tret[j] - dt; i < 2; i++, t += dt) {
if ((retflag2 = eclipse_where(t, ipl, starname, ifl, geopos, dcore, serr)) == ERR)
return retflag2;
if (n == 0)
dc[i] = dcore[4] / 2 + de / dcore[5] - dcore[2];
else if (n == 1)
dc[i] = fabs(dcore[3]) / 2 + de / dcore[6] - dcore[2];
else if (n == 2)
dc[i] = de / dcore[6] - dcore[2];
}
dt1 = dc[1] / ((dc[1] - dc[0]) / dt);
tret[j] -= dt1;
}
}
}
/*
* annular-total eclipses
*/
if (retflag & SE_ECL_TOTAL) {
if ((retflag2 = eclipse_where(tret[0], ipl, starname, ifl, geopos, dcore, serr)) == ERR)
return retflag2;
dc[0] = *dcore;
if ((retflag2 = eclipse_where(tret[4], ipl, starname, ifl, geopos, dcore, serr)) == ERR)
return retflag2;
dc[1] = *dcore;
if ((retflag2 = eclipse_where(tret[5], ipl, starname, ifl, geopos, dcore, serr)) == ERR)
return retflag2;
dc[2] = *dcore;
/* the maximum is always total, and there is either one or
* to times before and after, when the core shadow becomes
* zero and totality changes into annularity or vice versa.
*/
if (dc[0] * dc[1] < 0 || dc[0] * dc[2] < 0) {
retflag |= SE_ECL_ANNULAR_TOTAL;
retflag &= ~SE_ECL_TOTAL;
}
}
/* if eclipse is given but not wanted: */
if (!(ifltype & SE_ECL_TOTAL) && (retflag & SE_ECL_TOTAL)) {
/*t= tjd + direction * dadd;*/
t = tjd + direction * 20;
if (one_try) {
tret[0] = tjd;
return 0;
}
tjd = t;
goto next_try;
}
/* if annular_total eclipse is given but not wanted: */
if (!(ifltype & SE_ECL_ANNULAR_TOTAL) && (retflag & SE_ECL_ANNULAR_TOTAL)) {
/*t= tjd + direction * dadd;*/
t = tjd + direction * 20;
if (one_try) {
tret[0] = tjd;
return 0;
}
tjd = t;
goto next_try;
}
/*
* time of maximum eclipse at local apparent noon
*/
/* first, find out, if there is a solar transit
* between begin and end of eclipse */
k = 2;
for (i = 0; i < 2; i++) {
j = i + k;
tt = tret[j] + swe_deltat_ex(tret[j], ifl, serr);
if (calc_planet_star(tt, ipl, starname, iflag, ls, serr) == ERR)
return ERR;
if (swe_calc(tt, SE_MOON, iflag, lm, serr) == ERR)
return ERR;
dc[i] = swe_degnorm(ls[0] - lm[0]);
if (dc[i] > 180)
dc[i] -= 360;
}
if (dc[0] * dc[1] >= 0) /* no transit */
tret[1] = 0;
else {
tjd = tjds;
dt = 0.1;
dt1 = (tret[3] - tret[2]) / 2.0;
if (dt1 < dt)
dt = dt1 / 2.0;
for (j = 0;
dt > 0.01;
j++, dt /= 3) {
for (i = 0, t = tjd; i <= 1; i++, t -= dt) {
tt = t + swe_deltat_ex(t, ifl, serr);
if (calc_planet_star(tt, ipl, starname, iflag, ls, serr) == ERR)
return ERR;
if (swe_calc(tt, SE_MOON, iflag, lm, serr) == ERR)
return ERR;
dc[i] = swe_degnorm(ls[0] - lm[0]);
if (dc[i] > 180)
dc[i] -= 360;
if (dc[i] > 180)
dc[i] -= 360;
}
a = (dc[1] - dc[0]) / dt;
if (a < 1e-10)
break;
dt1 = dc[0] / a;
tjd += dt1;
}
tret[1] = tjd;
}
end_search_global:
return retflag;
/*
* the time of maximum occultation is practically identical
* with the time of maximum core shadow diameter.
*
* the time, when duration of totality is maximal,
* is not an interesting computation either. Near the maximum
* occulation, the time of totality can be the same by
* a second for hundreds of kilometers (for 10 minutes
* or more).
*
* for annular eclipses the maximum duration is close to the
* beginning and the end of the center lines, where is also
* the minimum of core shadow diameter.
*/
}
/* When is the next solar eclipse at a given geographical position?
* Note the uncertainty of Delta T for the remote past and for
* the future.
*
* retflag SE_ECL_TOTAL or SE_ECL_ANNULAR or SE_ECL_PARTIAL
* SE_ECL_VISIBLE,
* SE_ECL_MAX_VISIBLE,
* SE_ECL_1ST_VISIBLE, SE_ECL_2ND_VISIBLE
* SE_ECL_3ST_VISIBLE, SE_ECL_4ND_VISIBLE
*
* tret[0] time of maximum eclipse
* tret[1] time of first contact
* tret[2] time of second contact
* tret[3] time of third contact
* tret[4] time of forth contact
* tret[5] time of sun rise between first and forth contact
* tret[6] time of sun set beween first and forth contact
*
* attr[0] fraction of solar diameter covered by moon;
* with total/annular eclipses, it results in magnitude acc. to IMCCE.
* attr[1] ratio of lunar diameter to solar one
* attr[2] fraction of solar disc covered by moon (obscuration)
* attr[3] diameter of core shadow in km
* attr[4] azimuth of sun at tjd
* attr[5] true altitude of sun above horizon at tjd
* attr[6] apparent altitude of sun above horizon at tjd
* attr[7] elongation of moon in degrees
* attr[8] magnitude acc. to NASA;
* = attr[0] for partial and attr[1] for annular and total eclipses
* attr[9] saros series number
* attr[10] saros series member number
* declare as attr[20] at least !
*/
int32 CALL_CONV swe_sol_eclipse_when_loc(double tjd_start, int32 ifl,
double *geopos, double *tret, double *attr, int32 backward, char *serr)
{
int32 retflag = 0, retflag2 = 0;
double geopos2[20], dcore[10];
if (geopos[2] < SEI_ECL_GEOALT_MIN || geopos[2] > SEI_ECL_GEOALT_MAX) {
if (serr != NULL)
sprintf(serr, "location for eclipses must be between %.0f and %.0f m above sea", SEI_ECL_GEOALT_MIN, SEI_ECL_GEOALT_MAX);
return ERR;
}
ifl &= SEFLG_EPHMASK;
swi_set_tid_acc(tjd_start, ifl, 0, serr);
if ((retflag = eclipse_when_loc(tjd_start, ifl, geopos, tret, attr, backward, serr)) <= 0)
return retflag;
/*
* diameter of core shadow
*/
if ((retflag2 = eclipse_where(tret[0], SE_SUN, NULL, ifl, geopos2, dcore, serr)) == ERR)
return retflag2;
retflag |= (retflag2 & SE_ECL_NONCENTRAL);
attr[3] = dcore[0];
return retflag;
}
/* When is the next solar eclipse at a given geographical position?
* Note the uncertainty of Delta T for the remote past and for
* the future.
*
* retflag SE_ECL_TOTAL or SE_ECL_ANNULAR or SE_ECL_PARTIAL
* SE_ECL_VISIBLE,
* SE_ECL_MAX_VISIBLE,
* SE_ECL_1ST_VISIBLE, SE_ECL_2ND_VISIBLE
* SE_ECL_3ST_VISIBLE, SE_ECL_4ND_VISIBLE
* SE_ECL_OCC_BEG_DAYLIGHT, SE_ECL_OCC_END_DAYLIGHT
* The latter two indicate that the beginning or end of the occultation takes
* place during the day. If Venus is occulted, it may be observable with the
* naked eye; if other objects, it may be observable with telescopes.
*
* int32 ipl planet number of occulted body
* char* starname name of occulted star. Must be NULL or "", if a planetary
* occultation is to be calculated. For the use of this
* field, also see swe_fixstar().
* int32 ifl ephemeris flag. If you want to have only one conjunction
* of the moon with the body tested, add the following flag:
* backward |= SE_ECL_ONE_TRY. If this flag is not set,
* the function will search for an occultation until it
* finds one. For bodies with ecliptical latitudes > 5,
* the function may search unsuccessfully until it reaches
* the end of the ephemeris.
*
* for all other parameters, see function swe_sol_eclipse_when_loc().
*/
int32 CALL_CONV swe_lun_occult_when_loc(double tjd_start, int32 ipl, char *starname, int32 ifl,
double *geopos, double *tret, double *attr, int32 backward, char *serr)
{
int32 retflag = 0, retflag2 = 0;
double geopos2[20], dcore[10];
/* function calls for Pluto with asteroid number 134340
* are treated as calls for Pluto as main body SE_PLUTO */
if (geopos[2] < SEI_ECL_GEOALT_MIN || geopos[2] > SEI_ECL_GEOALT_MAX) {
if (serr != NULL)
sprintf(serr, "location for occultations must be between %.0f and %.0f m above sea", SEI_ECL_GEOALT_MIN, SEI_ECL_GEOALT_MAX);
return ERR;
}
if (ipl < 0) ipl = 0;
if (ipl == SE_AST_OFFSET + 134340)
ipl = SE_PLUTO;
ifl &= SEFLG_EPHMASK;
swi_set_tid_acc(tjd_start, ifl, 0, serr);
if ((retflag = occult_when_loc(tjd_start, ipl, starname, ifl, geopos, tret, attr, backward, serr)) <= 0)
return retflag;
/*
* diameter of core shadow
*/
if ((retflag2 = eclipse_where(tret[0], ipl, starname, ifl, geopos2, dcore, serr)) == ERR)
return retflag2;
retflag |= (retflag2 & SE_ECL_NONCENTRAL);
attr[3] = dcore[0];
return retflag;
}
static int32 eclipse_when_loc(double tjd_start, int32 ifl, double *geopos, double *tret, double *attr, int32 backward, char *serr)
{
int i, j, k, m;
int32 retflag = 0, retc;
double t, tjd, dt, dtint, K, T, T2, T3, T4, F, M, Mm;
double tjdr, tjds;
double E, Ff, A1, Om;
double xs[6], xm[6], ls[6], lm[6], x1[6], x2[6], dm, ds;
double rmoon, rsun, rsplusrm, rsminusrm;
double dc[3], dctr, dctrmin;
double twomin = 2.0 / 24.0 / 60.0;
double tensec = 10.0 / 24.0 / 60.0 / 60.0;
double twohr = 2.0 / 24.0;
double tenmin = 10.0 / 24.0 / 60.0;
double dt1 = 0, dt2 = 0, dtdiv, dtstart;
int32 iflag = SEFLG_EQUATORIAL | SEFLG_TOPOCTR | ifl;
int32 iflagcart = iflag | SEFLG_XYZ;
swe_set_topo(geopos[0], geopos[1], geopos[2]);
K = (int) ((tjd_start - J2000) / 365.2425 * 12.3685);
if (backward)
K++;
else
K--;
next_try:
T = K / 1236.85;
T2 = T * T; T3 = T2 * T; T4 = T3 * T;
Ff = F = swe_degnorm(160.7108 + 390.67050274 * K
- 0.0016341 * T2
- 0.00000227 * T3
+ 0.000000011 * T4);
if (Ff > 180)
Ff -= 180;
if (Ff > 21 && Ff < 159) { /* no eclipse possible */
if (backward)
K--;
else
K++;
goto next_try;
}
/* approximate time of geocentric maximum eclipse.
* formula from Meeus, German, p. 381 */
tjd = 2451550.09765 + 29.530588853 * K
+ 0.0001337 * T2
- 0.000000150 * T3
+ 0.00000000073 * T4;
M = swe_degnorm(2.5534 + 29.10535669 * K
- 0.0000218 * T2
- 0.00000011 * T3);
Mm = swe_degnorm(201.5643 + 385.81693528 * K
+ 0.1017438 * T2
+ 0.00001239 * T3
+ 0.000000058 * T4);
Om = swe_degnorm(124.7746 - 1.56375580 * K
+ 0.0020691 * T2
+ 0.00000215 * T3);
E = 1 - 0.002516 * T - 0.0000074 * T2;
A1 = swe_degnorm(299.77 + 0.107408 * K - 0.009173 * T2);
M *= DEGTORAD;
Mm *= DEGTORAD;
F *= DEGTORAD;
Om *= DEGTORAD;
A1 *= DEGTORAD;
tjd = tjd - 0.4075 * sin(Mm)
+ 0.1721 * E * sin(M);
swe_set_topo(geopos[0], geopos[1], geopos[2]);
dtdiv = 2;
dtstart = 0.5;
if (tjd < 1900000 || tjd > 2500000) /* because above formula is not good (delta t?) */
dtstart = 2;
for (dt = dtstart;
dt > 0.00001;
dt /= dtdiv) {
if (dt < 0.1)
dtdiv = 3;
for (i = 0, t = tjd - dt; i <= 2; i++, t += dt) {
/* this takes some time, but is necessary to avoid
* missing an eclipse */
if (swe_calc(t, SE_SUN, iflagcart, xs, serr) == ERR)
return ERR;
if (swe_calc(t, SE_SUN, iflag, ls, serr) == ERR)
return ERR;
if (swe_calc(t, SE_MOON, iflagcart, xm, serr) == ERR)
return ERR;
if (swe_calc(t, SE_MOON, iflag, lm, serr) == ERR)
return ERR;
dm = sqrt(square_sum(xm));
ds = sqrt(square_sum(xs));
for (k = 0; k < 3; k++) {
x1[k] = xs[k] / ds /*ls[2]*/;
x2[k] = xm[k] / dm /*lm[2]*/;
}
dc[i] = acos(swi_dot_prod_unit(x1, x2)) * RADTODEG;
}
find_maximum(dc[0], dc[1], dc[2], dt, &dtint, &dctr);
tjd += dtint + dt;
}
if (swe_calc(tjd, SE_SUN, iflagcart, xs, serr) == ERR)
return ERR;
if (swe_calc(tjd, SE_SUN, iflag, ls, serr) == ERR)
return ERR;
if (swe_calc(tjd, SE_MOON, iflagcart, xm, serr) == ERR)
return ERR;
if (swe_calc(tjd, SE_MOON, iflag, lm, serr) == ERR)
return ERR;
dctr = acos(swi_dot_prod_unit(xs, xm)) * RADTODEG;
rmoon = asin(RMOON / lm[2]) * RADTODEG;
rsun = asin(RSUN / ls[2]) * RADTODEG;
rsplusrm = rsun + rmoon;
rsminusrm = rsun - rmoon;
if (dctr > rsplusrm) {
if (backward)
K--;
else
K++;
goto next_try;
}
tret[0] = tjd - swe_deltat_ex(tjd, ifl, serr);
tret[0] = tjd - swe_deltat_ex(tret[0], ifl, serr); /* these two lines are an iteration! */
if ((backward && tret[0] >= tjd_start - 0.0001)
|| (!backward && tret[0] <= tjd_start + 0.0001)) {
if (backward)
K--;
else
K++;
goto next_try;
}
if (dctr < rsminusrm)
retflag = SE_ECL_ANNULAR;
else if (dctr < fabs(rsminusrm))
retflag = SE_ECL_TOTAL;
else if (dctr <= rsplusrm)
retflag = SE_ECL_PARTIAL;
dctrmin = dctr;
/* contacts 2 and 3 */
if (dctr > fabs(rsminusrm)) /* partial, no 2nd and 3rd contact */
tret[2] = tret[3] = 0;
else {
dc[1] = fabs(rsminusrm) - dctrmin;
for (i = 0, t = tjd - twomin; i <= 2; i += 2, t = tjd + twomin) {
if (swe_calc(t, SE_SUN, iflagcart, xs, serr) == ERR)
return ERR;
if (swe_calc(t, SE_MOON, iflagcart, xm, serr) == ERR)
return ERR;
dm = sqrt(square_sum(xm));
ds = sqrt(square_sum(xs));
rmoon = asin(RMOON / dm) * RADTODEG;
rmoon *= 0.99916; /* gives better accuracy for 2nd/3rd contacts */
rsun = asin(RSUN / ds) * RADTODEG;
rsminusrm = rsun - rmoon;
for (k = 0; k < 3; k++) {
x1[k] = xs[k] / ds /*ls[2]*/;
x2[k] = xm[k] / dm /*lm[2]*/;
}
dctr = acos(swi_dot_prod_unit(x1, x2)) * RADTODEG;
dc[i] = fabs(rsminusrm) - dctr;
}
find_zero(dc[0], dc[1], dc[2], twomin, &dt1, &dt2);
tret[2] = tjd + dt1 + twomin;
tret[3] = tjd + dt2 + twomin;
for (m = 0, dt = tensec; m < 2; m++, dt /= 10) {
for (j = 2; j <= 3; j++) {
if (swe_calc(tret[j], SE_SUN, iflagcart | SEFLG_SPEED, xs, serr) == ERR)
return ERR;
if (swe_calc(tret[j], SE_MOON, iflagcart | SEFLG_SPEED, xm, serr) == ERR)
return ERR;
for (i = 0; i < 2; i++) {
if (i == 1) {
for(k = 0; k < 3; k++) {
xs[k] -= xs[k+3] * dt;
xm[k] -= xm[k+3] * dt;
}
}
dm = sqrt(square_sum(xm));
ds = sqrt(square_sum(xs));
rmoon = asin(RMOON / dm) * RADTODEG;
rmoon *= 0.99916; /* gives better accuracy for 2nd/3rd contacts */
rsun = asin(RSUN / ds) * RADTODEG;
rsminusrm = rsun - rmoon;
for (k = 0; k < 3; k++) {
x1[k] = xs[k] / ds /*ls[2]*/;
x2[k] = xm[k] / dm /*lm[2]*/;
}
dctr = acos(swi_dot_prod_unit(x1, x2)) * RADTODEG;
dc[i] = fabs(rsminusrm) - dctr;
}
dt1 = -dc[0] / ((dc[0] - dc[1]) / dt);
tret[j] += dt1;
}
}
tret[2] -= swe_deltat_ex(tret[2], ifl, serr);
tret[3] -= swe_deltat_ex(tret[3], ifl, serr);
}
/* contacts 1 and 4 */
dc[1] = rsplusrm - dctrmin;
for (i = 0, t = tjd - twohr; i <= 2; i += 2, t = tjd + twohr) {
if (swe_calc(t, SE_SUN, iflagcart, xs, serr) == ERR)
return ERR;
if (swe_calc(t, SE_MOON, iflagcart, xm, serr) == ERR)
return ERR;
dm = sqrt(square_sum(xm));
ds = sqrt(square_sum(xs));
rmoon = asin(RMOON / dm) * RADTODEG;
rsun = asin(RSUN / ds) * RADTODEG;
rsplusrm = rsun + rmoon;
for (k = 0; k < 3; k++) {
x1[k] = xs[k] / ds /*ls[2]*/;
x2[k] = xm[k] / dm /*lm[2]*/;
}
dctr = acos(swi_dot_prod_unit(x1, x2)) * RADTODEG;
dc[i] = rsplusrm - dctr;
}
find_zero(dc[0], dc[1], dc[2], twohr, &dt1, &dt2);
tret[1] = tjd + dt1 + twohr;
tret[4] = tjd + dt2 + twohr;
for (m = 0, dt = tenmin; m < 3; m++, dt /= 10) {
for (j = 1; j <= 4; j += 3) {
if (swe_calc(tret[j], SE_SUN, iflagcart | SEFLG_SPEED, xs, serr) == ERR)
return ERR;
if (swe_calc(tret[j], SE_MOON, iflagcart | SEFLG_SPEED, xm, serr) == ERR)
return ERR;
for (i = 0; i < 2; i++) {
if (i == 1) {
for(k = 0; k < 3; k++) {
xs[k] -= xs[k+3] * dt;
xm[k] -= xm[k+3] * dt;
}
}
dm = sqrt(square_sum(xm));
ds = sqrt(square_sum(xs));
rmoon = asin(RMOON / dm) * RADTODEG;
rsun = asin(RSUN / ds) * RADTODEG;
rsplusrm = rsun + rmoon;
for (k = 0; k < 3; k++) {
x1[k] = xs[k] / ds /*ls[2]*/;
x2[k] = xm[k] / dm /*lm[2]*/;
}
dctr = acos(swi_dot_prod_unit(x1, x2)) * RADTODEG;
dc[i] = fabs(rsplusrm) - dctr;
}
dt1 = -dc[0] / ((dc[0] - dc[1]) / dt);
tret[j] += dt1;
}
}
tret[1] -= swe_deltat_ex(tret[1], ifl, serr);
tret[4] -= swe_deltat_ex(tret[4], ifl, serr);
/*
* visibility of eclipse phases
*/
for (i = 4; i >= 0; i--) { /* attr for i = 0 must be kept !!! */
if (tret[i] == 0)
continue;
if (eclipse_how(tret[i], SE_SUN, NULL, ifl, geopos[0], geopos[1], geopos[2],
attr, serr) == ERR)
return ERR;
/*if (retflag2 & SE_ECL_VISIBLE) { could be wrong for 1st/4th contact */
if (attr[6] > 0) { /* this is safe, sun above horizon, using app. alt. */
retflag |= SE_ECL_VISIBLE;
switch(i) {
case 0: retflag |= SE_ECL_MAX_VISIBLE; break;
case 1: retflag |= SE_ECL_1ST_VISIBLE; break;
case 2: retflag |= SE_ECL_2ND_VISIBLE; break;
case 3: retflag |= SE_ECL_3RD_VISIBLE; break;
case 4: retflag |= SE_ECL_4TH_VISIBLE; break;
default: break;
}
}
}
#if 1
if (!(retflag & SE_ECL_VISIBLE)) {
if (backward)
K--;
else
K++;
goto next_try;
}
#endif
if ((retc = swe_rise_trans(tret[1] - 0.001, SE_SUN, NULL, iflag, SE_CALC_RISE|SE_BIT_DISC_BOTTOM, geopos, 0, 0, &tjdr, serr)) == ERR)
return ERR;
if (retc == -2) /* circumpolar sun */
return retflag;
if ((retc = swe_rise_trans(tret[1] - 0.001, SE_SUN, NULL, iflag, SE_CALC_SET|SE_BIT_DISC_BOTTOM, geopos, 0, 0, &tjds, serr)) == ERR)
return ERR;
if (retc == -2) /* circumpolar sun */
return retflag;
if (tjds < tret[1] || (tjds > tjdr && tjdr > tret[4])) {
if (backward)
K--;
else
K++;
goto next_try;
}
if (tjdr > tret[1] && tjdr < tret[4]) {
tret[5] = tjdr;
if (!(retflag & SE_ECL_MAX_VISIBLE)) {
tret[0] = tjdr;
if ((retc = eclipse_how(tret[5], SE_SUN, NULL, ifl, geopos[0], geopos[1], geopos[2], attr, serr)) == ERR)
return ERR;
retflag &= ~(SE_ECL_TOTAL|SE_ECL_ANNULAR|SE_ECL_PARTIAL);
retflag |= (retc & (SE_ECL_TOTAL|SE_ECL_ANNULAR|SE_ECL_PARTIAL));
}
}
if (tjds > tret[1] && tjds < tret[4]) {
tret[6] = tjds;
if (!(retflag & SE_ECL_MAX_VISIBLE)) {
tret[0] = tjds;
if ((retc = eclipse_how(tret[6], SE_SUN, NULL, ifl, geopos[0], geopos[1], geopos[2], attr, serr)) == ERR)
return ERR;
retflag &= ~(SE_ECL_TOTAL|SE_ECL_ANNULAR|SE_ECL_PARTIAL);
retflag |= (retc & (SE_ECL_TOTAL|SE_ECL_ANNULAR|SE_ECL_PARTIAL));
}
}
return retflag;
}
static int32 occult_when_loc(
double tjd_start, int32 ipl, char *starname,
int32 ifl, double *geopos, double *tret, double *attr,
int32 backward, char *serr)
{
int i, j, k, m;
int32 retflag = 0, retc;
double t, tjd, dt, dtint;
double tjdr, tjds;
double xs[6], xm[6], ls[6], lm[6], x1[6], x2[6], dm, ds;
double rmoon, rsun, rsplusrm, rsminusrm;
double dc[20], dctr, dctrmin;
double twomin = 2.0 / 24.0 / 60.0;
double tensec = 10.0 / 24.0 / 60.0 / 60.0;
double twohr = 2.0 / 24.0;
double tenmin = 10.0 / 24.0 / 60.0;
double dt1 = 0, dt2 = 0, dtdiv, dtstart;
double dadd2 = 1;
double drad, dl;
AS_BOOL is_partial = FALSE;
int32 iflag = SEFLG_TOPOCTR | ifl;
int32 iflaggeo = iflag & ~SEFLG_TOPOCTR;
int32 iflagcart = iflag | SEFLG_XYZ;
int direction = 1;
int32 one_try = backward & SE_ECL_ONE_TRY;
AS_BOOL stop_after_this = FALSE;
backward &= 1L;
retflag = 0;
swe_set_topo(geopos[0], geopos[1], geopos[2]);
for (i = 0; i <= 9; i++)
tret[i] = 0;
if (backward)
direction = -1;
t = tjd_start;
tjd = tjd_start;
next_try:
is_partial = FALSE;
if (calc_planet_star(t, ipl, starname, iflaggeo, ls, serr) == ERR)
return ERR;
/* fixed stars with an ecliptic latitude > 7 or < -7 cannot have
* an occultation. Even lunar parallax andd proper motion of star
* will never allow it. */
if (fabs(ls[1]) > 7 && starname != NULL && *starname != '\0') {
if (serr != NULL)
sprintf(serr, "occultation never occurs: star %s has ecl. lat. %.1f", starname, ls[1]);
return ERR;
}
if (swe_calc(t, SE_MOON, iflaggeo, lm, serr) == ERR)
return ERR;
dl = swe_degnorm(ls[0] - lm[0]);
if (direction < 0)
dl -= 360;
/* get rough conjunction in ecliptic longitude */
while (fabs(dl) > 0.1) {
t += dl / 13;
if (calc_planet_star(t, ipl, starname, iflaggeo, ls, serr) == ERR)
return ERR;
if (swe_calc(t, SE_MOON, iflaggeo, lm, serr) == ERR)
return ERR;
dl = swe_degnorm(ls[0] - lm[0]);
if (dl > 180) dl -= 360;
}
tjd = t;
/* difference in latitude too big for an occultation */
drad = fabs(ls[1] - lm[1]);
if (drad > 2) {
if (one_try) {
tret[0] = t + direction; /* return a date suitable for next try */
return 0;
}
t += direction * 20;
tjd = t;
goto next_try;
}
/*
* radius of planet disk in AU
*/
if (starname != NULL && *starname != '\0')
drad = 0;
else if (ipl < NDIAM)
drad = pla_diam[ipl] / 2 / AUNIT;
else if (ipl > SE_AST_OFFSET)
drad = swed.ast_diam / 2 * 1000 / AUNIT; /* km -> m -> AU */
else
drad = 0;
/* now find out, if there is an occultation at our geogr. location */
dtdiv = 3;
dtstart = dadd2; /* formerly 0.2 */
for (dt = dtstart;
dt > 0.00001;
dt /