-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathswecl.cpp
6436 lines (6343 loc) · 221 KB
/
swecl.cpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
/* SWISSEPH
Ephemeris computations
Author: Dieter Koch
************************************************************/
/* Copyright (C) 1997 - 2021 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 Affero General Public License (AGPL)
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 AGPL 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 AGPL or a compatible license.
See https://www.gnu.org/licenses/agpl-3.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 "astrolog.h"
#ifdef SWISS
#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;};
// Saros cycle numbers of solar eclipses with date of initial
// eclipse. Table was derived from the table of the Nasa Eclipse Web Site:
// https://eclipse.gsfc.nasa.gov/SEsaros/SEsaros0-180.html
// Note, for eclipse dates =< 15 Feb -1604 and eclipse dates
// >= 2 Sep 2666, Saros cycle numbers cannot always be given.
#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 */
};
// Saros cycle numbers of lunar eclipses with date of initial
// eclipse. Table was derived from the table of the Nasa Eclipse Web Site:
// https://eclipse.gsfc.nasa.gov/LEsaros/LEsaroscat.html
// Note, for eclipse dates =< 29 April -1337 and eclipse dates
// >= 10 Aug 2892, Saros cycle numbers cannot always be given.
#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
* 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_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;
}
/*
double tjd_ut, time, Jul. day UT
int32 ipl, planet number
char* starname, star name, must be NULL or "" if not a star
int32 ifl, ephemeris flag
double *geopos, return array, 2 doubles, geo. long. and lat.
eastern longitude is positive,
western longitude is negative,
northern latitude is positive,
double *attr, return array, 20 doubles, see below
char *serr return error string
array attr:
attr[0] fraction of object's diameter covered by moon (magnitude)
attr[1] ratio of lunar diameter to object's diameter
attr[2] fraction of object's disc covered by moon (obscuration)
attr[3] diameter of core shadow in km
attr[4] azimuth of object at tjd
attr[5] true altitude of object above horizon at tjd
attr[6] apparent altitude of object above horizon at tjd
attr[7] angular distance of moon from object in degrees
*/
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 retc = OK;
if (starname == NULL || *starname == '\0') {
retc = swe_calc(tjd_et, ipl, iflag, x, serr);
} else {
retc = swe_fixstar(starname, tjd_et, iflag, x, serr);
}
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, i;
double dcore[10], ls[6], xaz[6];
double geopos2[20];
for (i = 0; i <= 10; i++)
attr[i] = 0;
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;
if (retflag == 0) {
for (i = 0; i <= 3; i++)
attr[i] = 0;
for (i = 8; i <= 10; i++)
attr[i] = 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)