Initial import of version 2.10

This commit is contained in:
Gergely Polonkai 2021-04-27 08:06:35 +02:00
commit ded7b26f78
No known key found for this signature in database
GPG Key ID: 2D2885533B869ED4
109 changed files with 682348 additions and 0 deletions

7
Makefile Normal file
View File

@ -0,0 +1,7 @@
all:
cd src && ${MAKE}
install:
cd src && ${MAKE} install
cd doc && ${MAKE} install
cd data && ${MAKE} install

69
data/Makefile Normal file
View File

@ -0,0 +1,69 @@
SE1_FILES = \
seas_00.se1 \
seas_06.se1 \
seas_12.se1 \
seas_18.se1 \
seas_24.se1 \
seas_30.se1 \
seas_36.se1 \
seas_42.se1 \
seas_48.se1 \
seasm06.se1 \
seasm12.se1 \
seasm18.se1 \
seasm24.se1 \
seasm30.se1 \
seasm36.se1 \
seasm42.se1 \
seasm48.se1 \
seasm54.se1 \
semo_00.se1 \
semo_06.se1 \
semo_12.se1 \
semo_18.se1 \
semo_24.se1 \
semo_30.se1 \
semo_36.se1 \
semo_42.se1 \
semo_48.se1 \
semom06.se1 \
semom12.se1 \
semom18.se1 \
semom24.se1 \
semom30.se1 \
semom36.se1 \
semom42.se1 \
semom48.se1 \
semom54.se1 \
sepl_00.se1 \
sepl_06.se1 \
sepl_12.se1 \
sepl_18.se1 \
sepl_24.se1 \
sepl_30.se1 \
sepl_36.se1 \
sepl_42.se1 \
sepl_48.se1 \
seplm06.se1 \
seplm12.se1 \
seplm18.se1 \
seplm24.se1 \
seplm30.se1 \
seplm36.se1 \
seplm42.se1 \
seplm48.se1 \
seplm54.se1 \
$(NULL)
OTHER_FILES = \
seorbel.txt \
seasnam.txt \
sefstars.txt \
seleapsec.txt \
$(NULL)
all:
install:
mkdir -p $(DESTDIR)/usr/share/libswe
install -m 0644 $(SE1_FILES) $(OTHER_FILES) $(DESTDIR)/usr/share/libswe

BIN
data/seas_00.se1 Normal file

Binary file not shown.

BIN
data/seas_06.se1 Normal file

Binary file not shown.

BIN
data/seas_12.se1 Normal file

Binary file not shown.

BIN
data/seas_18.se1 Normal file

Binary file not shown.

BIN
data/seas_24.se1 Normal file

Binary file not shown.

BIN
data/seas_30.se1 Normal file

Binary file not shown.

BIN
data/seas_36.se1 Normal file

Binary file not shown.

BIN
data/seas_42.se1 Normal file

Binary file not shown.

BIN
data/seas_48.se1 Normal file

Binary file not shown.

BIN
data/seasm06.se1 Normal file

Binary file not shown.

BIN
data/seasm12.se1 Normal file

Binary file not shown.

BIN
data/seasm18.se1 Normal file

Binary file not shown.

BIN
data/seasm24.se1 Normal file

Binary file not shown.

BIN
data/seasm30.se1 Normal file

Binary file not shown.

BIN
data/seasm36.se1 Normal file

Binary file not shown.

BIN
data/seasm42.se1 Normal file

Binary file not shown.

BIN
data/seasm48.se1 Normal file

Binary file not shown.

BIN
data/seasm54.se1 Normal file

Binary file not shown.

547966
data/seasnam.txt Normal file

File diff suppressed because it is too large Load Diff

1589
data/sefstars.txt Normal file

File diff suppressed because it is too large Load Diff

6
data/seleapsec.txt Normal file
View File

@ -0,0 +1,6 @@
# This file contains the dates of leap seconds to be taken into account
# by the Swiss Ephemeris.
# For each new leap second add the date of its insertion in the format
# yyyymmdd, e.g. "20081231" for 31 december 2008.
# The leap second is inserted at the end of the day.
20081231

BIN
data/semo_00.se1 Normal file

Binary file not shown.

BIN
data/semo_06.se1 Normal file

Binary file not shown.

BIN
data/semo_12.se1 Normal file

Binary file not shown.

BIN
data/semo_18.se1 Normal file

Binary file not shown.

BIN
data/semo_24.se1 Normal file

Binary file not shown.

BIN
data/semo_30.se1 Normal file

Binary file not shown.

BIN
data/semo_36.se1 Normal file

Binary file not shown.

BIN
data/semo_42.se1 Normal file

Binary file not shown.

BIN
data/semo_48.se1 Normal file

Binary file not shown.

BIN
data/semom06.se1 Normal file

Binary file not shown.

BIN
data/semom12.se1 Normal file

Binary file not shown.

BIN
data/semom18.se1 Normal file

Binary file not shown.

BIN
data/semom24.se1 Normal file

Binary file not shown.

BIN
data/semom30.se1 Normal file

Binary file not shown.

BIN
data/semom36.se1 Normal file

Binary file not shown.

BIN
data/semom42.se1 Normal file

Binary file not shown.

BIN
data/semom48.se1 Normal file

Binary file not shown.

BIN
data/semom54.se1 Normal file

Binary file not shown.

87
data/seorbel.txt Normal file
View File

@ -0,0 +1,87 @@
# Orbital elements of ficticious planets
# 27 Jan. 2000
#
# This file is part of the Swiss Ephemeris, from Version 1.52 on.
#
# Warning! These planets do not exist!
#
# The user can add his or her own elements.
# 960 is the maximum number of ficticious planets.
#
# The elements order is as follows:
# 1. epoch of elements (Julian day)
# 2. equinox (Julian day or "J1900" or "B1950" or "J2000")
# 3. mean anomaly at epoch
# 4. semi-axis
# 5. eccentricity
# 6. argument of perihelion (ang. distance of perihelion from node)
# 7. ascending node
# 8. inclination
# 9. name of planet
#
# use '#' for comments
# to compute a body with swe_calc(), use planet number
# ipl = SE_FICT_OFFSET_1 + number_of_elements_set,
# e.g. number of Kronos is ipl = 39 + 4 = 43
#
# Witte/Sieggruen planets, refined by James Neely
J1900, J1900, 163.7409, 40.99837, 0.00460, 171.4333, 129.8325, 1.0833, Cupido # 1
J1900, J1900, 27.6496, 50.66744, 0.00245, 148.1796, 161.3339, 1.0500, Hades # 2
J1900, J1900, 165.1232, 59.21436, 0.00120, 299.0440, 0.0000, 0.0000, Zeus # 3
J1900, J1900, 169.0193, 64.81690, 0.00305, 208.8801, 0.0000, 0.0000, Kronos # 4
J1900, J1900, 138.0533, 70.29949, 0.00000, 0.0000, 0.0000, 0.0000, Apollon # 5
J1900, J1900, 351.3350, 73.62765, 0.00000, 0.0000, 0.0000, 0.0000, Admetos # 6
J1900, J1900, 55.8983, 77.25568, 0.00000, 0.0000, 0.0000, 0.0000, Vulcanus # 7
J1900, J1900, 165.5163, 83.66907, 0.00000, 0.0000, 0.0000, 0.0000, Poseidon # 8
#
# Isis-Transpluto; elements from "Die Sterne" 3/1952, p. 70ff.
# Strubell does not give an equinox. 1945 is taken in order to
# reproduce the as best as ASTRON ephemeris. (This is a strange
# choice, though.)
# The epoch according to Strubell is 1772.76.
# 1772 is a leap year!
# The fraction is counted from 1 Jan. 1772
2368547.66, 2431456.5, 0.0, 77.775, 0.3, 0.7, 0, 0, Isis-Transpluto # 9
# Nibiru, elements from Christian Woeltge, Hannover
1856113.380954, 1856113.380954, 0.0, 234.8921, 0.981092, 103.966, -44.567, 158.708, Nibiru # 10
# Harrington, elements from Astronomical Journal 96(4), Oct. 1988
2374696.5, J2000, 0.0, 101.2, 0.411, 208.5, 275.4, 32.4, Harrington # 11
# according to W.G. Hoyt, "Planets X and Pluto", Tucson 1980, p. 63
2395662.5, 2395662.5, 34.05, 36.15, 0.10761, 284.75, 0, 0, Leverrier (Neptune) # 12
2395662.5, 2395662.5, 24.28, 37.25, 0.12062, 299.11, 0, 0, Adams (Neptune) # 13
2425977.5, 2425977.5, 281, 43.0, 0.202, 204.9, 0, 0, Lowell (Pluto) # 14
2425977.5, 2425977.5, 48.95, 55.1, 0.31, 280.1, 100, 15, Pickering (Pluto) # 15
# intramercurian hypothetical Vulcan acc. to L.H. Weston
J1900,JDATE, 252.8987988 + 707550.7341 * T, 0.13744, 0.019, 322.212069+1670.056*T, 47.787931-1670.056*T, 7.5, Vulcan # 16
# Selena/White Moon
J2000,JDATE, 242.2205555 + 5143.5418158 * T, 0.05280098949, 0.0, 0.0, 0.0, 0.0, Selena/White Moon, geo # 17
# Hypothetical planet Proserpina, according to http://www.geocities.com/Hollywood/Academy/7519/proserpina.html
# J1900, 170.73 + 51.05 * T
J1900,JDATE, 170.73, 79.225630, 0, 0, 0, 0, Proserpina #18
# Waldemath's Second Earth Moon
# Elements were derived by D.Koch from Waldemaths original elements as given in
# David Walters' book on Vulcan. They differ from Solar Fire (Graham Dawsons)
# elements, which are based on the assumption that the "mean longitude" given
# by Waldemath is an observation (a true longitude)
# Neither Swisseph nor Solar fire elements agree with Delphine Jay's ephemeris,
# which is obviously wrong.
2414290.95827875,2414290.95827875, 70.3407215 + 109023.2634989 * T, 0.0068400705250028, 0.1587, 8.14049594 + 2393.47417444 * T, 136.24878256 - 1131.71719709 * T, 2.5, Waldemath, geo # 19
# Colin R. Nicholl's Comet, according to "The Great Christ Comet", p. 223.
1719500.7, J2000, 0.0, 1190.0, 0.9999, 9.47, 200.08, 178.3, Christ Comet # 20
# Planet 9, according to: Fienga & alii, Constraints on the location of a
# possible 9th planet, Astronomy & Astrophysics no. FiengaLaskar2016R5. The
# authors provide true anomaly 117.8. Mean anomaly derived by D. Koch
2457388.5,J2000,45.5272966,700,0.6,150,113,30,Planet_9 # 21
# Comet Halley 12 BCE, only good near 12 BCE;
# Yeomans/Kiang, Mon. Not. R. astr. Soc. (1981), 197, p. 643
1717323.34852, B1950, 0.0, 17.99261849, 0.9673664, 92.54399, 35.19064, 163.58392, Comet Halley 12 BCE # 22
# 2015 RR245
2457600.5,J2000,322.50413,81.2891975,0.5852663,261.41753,211.67680,7.57643,2015 RR245 # 23
#
# The following elements are for test only
# (Selena without T)
J2000,JDATE, 242.2205555, 0.05279142865925, 0.0, 0.0, 0.0, 0.0, Selena/White Moon, geo # 17
# (Selena with T, gives exactly the same position)
J2000,JDATE, 242.2205555 + 5143.5418158 * T, 0.05279142865925, 0.0, 0.0, 0.0, 0.0, Selena/White Moon with T Terms, geo # 17
J2000, JDATE, 174.794787 + 149472.5157715 * T, 0.38709831, 0.20563175 + 0.000020406 * T, 29.125226 + 0.3702885 * T, 48.330893 + 1.186189 * T, 7.004986 + 0.0018215 * T, Mercury elem. for equ. of date # 18
J2000, J2000, 174.794787 + 149472.5157715 * T, 0.38709831, 0.20563175 + 0.000020406 * T, 29.125226 + 0.2842872 * T, 48.330893 - 0.1254229 * T, 7.004986 - 0.0059516 * T, Mercury Test J2000 Elements# 18

BIN
data/sepl_00.se1 Normal file

Binary file not shown.

BIN
data/sepl_06.se1 Normal file

Binary file not shown.

BIN
data/sepl_12.se1 Normal file

Binary file not shown.

BIN
data/sepl_18.se1 Normal file

Binary file not shown.

BIN
data/sepl_24.se1 Normal file

Binary file not shown.

BIN
data/sepl_30.se1 Normal file

Binary file not shown.

BIN
data/sepl_36.se1 Normal file

Binary file not shown.

BIN
data/sepl_42.se1 Normal file

Binary file not shown.

BIN
data/sepl_48.se1 Normal file

Binary file not shown.

BIN
data/seplm06.se1 Normal file

Binary file not shown.

BIN
data/seplm12.se1 Normal file

Binary file not shown.

BIN
data/seplm18.se1 Normal file

Binary file not shown.

BIN
data/seplm24.se1 Normal file

Binary file not shown.

BIN
data/seplm30.se1 Normal file

Binary file not shown.

BIN
data/seplm36.se1 Normal file

Binary file not shown.

BIN
data/seplm42.se1 Normal file

Binary file not shown.

BIN
data/seplm48.se1 Normal file

Binary file not shown.

BIN
data/seplm54.se1 Normal file

Binary file not shown.

5
doc/Makefile Normal file
View File

@ -0,0 +1,5 @@
all:
install:
mkdir -p $(DESTDIR)/usr/share/doc/libswe
install -m 0644 sweph.gif swephin.gif swephprg.htm swephprg.pdf swisseph.htm swisseph.pdf $(DESTDIR)/usr/share/doc/libswe

BIN
doc/sweph.gif Normal file

Binary file not shown.

After

Width:  |  Height:  |  Size: 6.0 KiB

BIN
doc/swephin.gif Normal file

Binary file not shown.

After

Width:  |  Height:  |  Size: 9.3 KiB

36663
doc/swephprg.htm Normal file

File diff suppressed because it is too large Load Diff

BIN
doc/swephprg.pdf Normal file

Binary file not shown.

30561
doc/swisseph.htm Normal file

File diff suppressed because it is too large Load Diff

BIN
doc/swisseph.pdf Normal file

Binary file not shown.

71
perl/Changes Normal file
View File

@ -0,0 +1,71 @@
Revision history for Perl extension SwissEph.
First created on 3 April 2009, with Swiss Ephemeris 1.76.
Updated on 27 January 2011:
- Function swe_utc_time_zone() added for Swiss Ephemeris 1.77
- Delta t updated
- Fixed stars: algorithm in swe_fixstar() was improved, new input file sefstars.txt was created
Updated on 7 September 2012:
- Precession algorithm was updated
- Function swe_rise_trans_true_hor() added for Swiss Ephemeris 1.78
- several minor bug fixes (see documentation swephprg.pdf)
Updated on 23 April 2013, Swiss Ephemeris 1.79
- Improved precision in eclipse calculations: 2nd and 3rd contact with solar eclipses, penumbral and partial phases with lunar eclipses.
- Bug fix in function swe_sol_eclipse_when_loc().If the local maximum eclipse occurs at sunset or sunrise, tret[0] now gives the moment when the lower limb of the Sun touches the horizon. This was not correctly implemented in former versions
- Bug fix in Perl functions swe_house() etc. These functions had crashed with a segmention violation if called with the house parameter G.
- Bug fix in Perl function swe_utc_to_jd(), where gregflag had been read from the 4th instead of the 6th parameter.
- Bug fix in Perl functions to do with date conversion. The default mechanism for gregflag was buggy.
- For Hindu astrologers, some more ayanamshas were added that are related to Suryasiddhanta and Aryabhata and are of historical interest.
Updated on 31 August 2012
- Function swe_house_name() added for Swiss Ephemeris 1.80
Updated on 11 February 2014, for Swiss Ephemeris 2.00:
- Major release, Swisseph can handle JPL ephemeris DE431
- The following functions were added:
swe_lun_eclipse_when_loc() searches locally visible lunar eclipse
swe_lmt_to_lat() conversion of LMT to LAT
swe_lat_to_lmt() conversion of LAT to LMT
Updated on 18 March 2015, Swiss Ephemeris 2.01:
- a few constants added to SwissEph.pm
Updated on 5 August 2015, Swiss Ephemeris 2.02:
- The following functions were added:
swe_deltat_ex() provides Delta T adjusted to ephemeris chosen
swe_get_ayanamsa_ex() provides ayanamsa depending on ephemeris chosen
swe_get_ayanamsa_ex_ut() provides ayanamsa depending on ephemeris chosen
Updated on 14 August 2015, Swiss Ephemeris 2.02.01:
- see documentation swephprg.pdf
Updated on 16 Oct 2015, Swiss Ephemeris 2.03:
- see documentation swephprg.pdf
Updated on 21 Oct 2015, Swiss Ephemeris 2.04:
- see documentation swephprg.pdf
Updated on 26 May 2016, Swiss Ephemeris 2.05:
- see documentation swephprg.pdf
Updated on 12 Jan 2018, Swiss Ephemeris 2.07:
- see documentation swephprg.pdf
Updated on 19 Jun 2019, Swiss Ephemeris 2.08:
- a few defines (constants) were added in SwissEph.pm
Updated on 9 Aug 2019, Swiss Ephemeris 2.08:
- eclipse functions return saros series and eclipse numbers
Updated on 23 Jul 2020, Swiss Ephemeris 2.09:
- see documentation swephprg.pdf
Updated on 18 Aug 2020, Swiss Ephemeris 2.09.02:
- see documentation swephprg.pdf
Updated on 10 Dec 2020, Swiss Ephemeris 2.10:
- see documentation swephprg.pdf

85
perl/MANIFEST Normal file
View File

@ -0,0 +1,85 @@
Changes
Makefile.PL
MANIFEST
ppport.h
sweodef.h
swephexp.h
README
SwissEph.xs
t/SwissEph.t
lib/SwissEph.pm
blib/arch/auto/SwissEph/SwissEph.so
fix_selinux
/usr/local/lib64/swe/libswe.so.2.04
/usr/local/lib64/swe/libswe.so.2.05c
/usr/local/lib64/swe/libswe.so.2.05c
/usr/local/lib64/swe/libswe.so.2.05c
/usr/local/lib64/swe/libswe.so.2.05
/usr/local/lib64/swe/libswe.so.2.05.01
/usr/local/lib64/swe/libswe.so.2.05.02b02
/usr/local/lib64/swe/libswe.so.2.05.02b02
/usr/local/lib64/swe/libswe.so.2.05.02b03
/usr/local/lib64/swe/libswe.so.2.05.02b02
/usr/local/lib64/swe/libswe.so.2.05.02b04
/usr/local/lib64/swe/libswe.so.2.05.02b04
/usr/local/lib64/swe/libswe.so.2.05.02b05
/usr/local/lib64/swe/libswe.so.2.05.02b05
/usr/local/lib64/swe/libswe.so.2.05.02b05
/usr/local/lib64/swe/libswe.so.2.05.02b05
/usr/local/lib64/swe/libswe.so.2.05.02b05
/usr/local/lib64/swe/libswe.so.2.05.02b05
/usr/local/lib64/swe/libswe.so.2.05.02b05
/usr/local/lib64/swe/libswe.so.2.05.02b05
/usr/local/lib64/swe/libswe.so.2.05.02b05
/usr/local/lib64/swe/libswe.so.2.05.02b06
/usr/local/lib64/swe/libswe.so.2.06
/usr/local/lib64/swe/libswe.so.2.06.01b02
/usr/local/lib64/swe/libswe.so.2.06.01b03
/usr/local/lib64/swe/libswe.so.2.06.01b03
/usr/local/lib64/swe/libswe.so.2.06.01b03
/usr/local/lib64/swe/libswe.so.2.06.01b03
/usr/local/lib64/swe/libswe.so.2.07
/usr/local/lib64/swe/libswe.so.2.07
/usr/local/lib64/swe/libswe.so.2.07
/usr/local/lib64/swe/libswe.so.2.07
/usr/local/lib64/swe/libswe.so.2.07
/usr/local/lib64/swe/libswe.so.2.07
/usr/local/lib64/swe/libswe.so.2.07
/usr/local/lib64/swe/libswe.so.2.07
/usr/local/lib64/swe/libswe.so.2.07
/usr/local/lib64/swe/libswe.so.2.07
/usr/local/lib64/swe/libswe.so.2.07
/usr/local/lib64/swe/libswe.so.2.07
/usr/local/lib64/swe/libswe.so.2.07
/usr/local/lib64/swe/libswe.so.2.07
/usr/local/lib64/swe/libswe.so.2.07
/usr/local/lib64/swe/libswe.so.2.07
/usr/local/lib64/swe/libswe.so.2.07
/usr/local/lib64/swe/libswe.so.2.07
/usr/local/lib64/swe/libswe.so.2.07
/usr/local/lib64/swe/libswe.so.2.07
/usr/local/lib64/swe/libswe.so.2.07
/usr/local/lib64/swe/libswe.so.2.07
/usr/local/lib64/swe/libswe.so.2.07
/usr/local/lib64/swe/libswe.so.2.07
/usr/local/lib64/swe/libswe.so.2.07
/usr/local/lib64/swe/libswe.so.2.07.01
/usr/local/lib64/swe/libswe.so.2.05c
/usr/local/lib64/swe/libswe.so.2.07.02
/usr/local/lib64/swe/libswe.so.2.07.02
/usr/local/lib64/swe/libswe.so.2.08
/usr/local/lib64/swe/libswe.so.2.08
/usr/local/lib64/swe/libswe.so.2.08
/usr/local/lib64/swe/libswe.so.2.08
/usr/local/lib64/swe/libswe.so.2.08
/usr/local/lib64/swe/libswe.so.2.08
/usr/local/lib64/swe/libswe.so.2.08.00a
/usr/local/lib64/swe/libswe.so.2.09
/usr/local/lib64/swe/libswe.so.2.09
/usr/local/lib64/swe/libswe.so.2.09
/usr/local/lib64/swe/libswe.so.2.09.01
/usr/local/lib64/swe/libswe.so.2.09.02
/usr/local/lib64/swe/libswe.so.2.09.03
/usr/local/lib64/swe/libswe.so.2.10
META.yml Module YAML meta-data (added by MakeMaker)
META.json Module JSON meta-data (added by MakeMaker)

39
perl/META.json Normal file
View File

@ -0,0 +1,39 @@
{
"abstract" : "The Swiss Ephemeris made accessible for Perl \r",
"author" : [
"Dieter <dieter@astro.ch>"
],
"dynamic_config" : 1,
"generated_by" : "ExtUtils::MakeMaker version 6.68, CPAN::Meta::Converter version 2.120921",
"license" : [
"unknown"
],
"meta-spec" : {
"url" : "http://search.cpan.org/perldoc?CPAN::Meta::Spec",
"version" : "2"
},
"name" : "SwissEph",
"no_index" : {
"directory" : [
"t",
"inc"
]
},
"prereqs" : {
"build" : {
"requires" : {
"ExtUtils::MakeMaker" : "0"
}
},
"configure" : {
"requires" : {
"ExtUtils::MakeMaker" : "0"
}
},
"runtime" : {
"requires" : {}
}
},
"release_status" : "stable",
"version" : "2.10"
}

21
perl/META.yml Normal file
View File

@ -0,0 +1,21 @@
---
abstract: "The Swiss Ephemeris made accessible for Perl \r"
author:
- 'Dieter <dieter@astro.ch>'
build_requires:
ExtUtils::MakeMaker: 0
configure_requires:
ExtUtils::MakeMaker: 0
dynamic_config: 1
generated_by: 'ExtUtils::MakeMaker version 6.68, CPAN::Meta::Converter version 2.120921'
license: unknown
meta-spec:
url: http://module-build.sourceforge.net/META-spec-v1.4.html
version: 1.4
name: SwissEph
no_index:
directory:
- t
- inc
requires: {}
version: 2.10

21
perl/Makefile.PL Executable file
View File

@ -0,0 +1,21 @@
#!/usr/local/bin/perl
use 5.008007;
use ExtUtils::MakeMaker;
# See lib/ExtUtils/MakeMaker.pm for details of how to influence
# the contents of the Makefile that is written.
WriteMakefile(
NAME => 'SwissEph',
VERSION_FROM => 'lib/SwissEph.pm', # finds $VERSION
PREREQ_PM => {}, # e.g., Module::Name => 1.1
($] >= 5.005 ? ## Add these new keywords supported since 5.005
(ABSTRACT_FROM => 'lib/SwissEph.pm', # retrieve abstract from module
AUTHOR => 'Dieter <dieter@astro.ch>') : ()),
#LIBS => ['-lC:\sweph\bin\swedll32'], # e.g., '-lm'
LIBS => ['-L/usr/local/lib64/swe -lswe'], # e.g., '-lm'
#LIBS => ['-lswe'], # e.g., '-lm'
DEFINE => '', # e.g., '-DHAVE_SOMETHING'
INC => '-I.', # e.g., '-I. -I/usr/include/other'
# DISTVNAME => 'perl_swisseph', #
# Un-comment this if you add C files to link with later:
# OBJECT => '$(O_FILES)', # link all the C files too
);

112
perl/README Normal file
View File

@ -0,0 +1,112 @@
SwissEph version 2.056
=====================
The README is used to introduce the module and provide instructions on
how to install the module, any machine dependencies it may have (for
example C compilers and installed libraries) and any other information
that should be provided before the module is installed.
INSTALLATION
On Linux:
=========
Before installing this module, you have to install a shared library of
the Swiss Ephemeris functions.
Unpack http://www.astro.com/ftp/swisseph/swe_unix_src_2.06.tar.gz
or whatever the latest version is, or the one you want to use.
There will be a src directory in the unpacked files. Go there
and run
make libswe.so to build the shared library. Copy this file to
/usr/local/lib64/swe
Usually, this will also require that you create a file (as root)
/etc/ld.so.conf.d/swisseph.conf
with the content line
/usr/local/lib64/swe
Now run ldconfig(8) (as root) in order to generate the
symbolic links required.
Now you can install the Perl module for the Swiss Ephemeris.
Type the following:
perl Makefile.PL
make
make test
make install (requires root)
On Windows:
===========
Before installing this module, you have to install the following things
on your computer:
- the Swiss Ephemeris DLL swedll32.dll. If haven't done this yet,
download the Swiss Ephemeris from the download area at
www.astro.com/swisseph and unpack it.
- Visual C++ Express Edition, which can be downloaded for free from the
http://www.microsoft.com/express/download/.
After that you can install the Perl module. In the current directory
(where you read this README file), open the file Makefile.PL and
fix the LIBS parameter. It must contain the directory to the Swiss Ephemeris
DLL. After that run:
perl Makefile.PL
nmake
nmake test
nmake install
COPYRIGHT AND LICENCE
Copyright (C) 2017 Astrodienst, Zurich, Switzerland.
This library is free software; you can redistribute it and/or modify
it under the same terms as Perl itself, either Perl version 5.8.7 or,
at your option, any later version of Perl 5 you may have available.
-----------------------
Update 23-march-2016, example how to use PerlSwissEph-2.04.tar.gz
download PerlSwissEph-2.04.tar.gz
unpack it with
tar xzvf PerlSwissEph-2.04.tar.gz
it creates a directory SwissEph-2.04
cd SwissEph-2.04
perl Makefile.PL (to create Makefile)
there is a warning:
Warning: the following files are missing in your kit:
/usr/local/lib64/swe/libswe.so.2.04
The tarball contains ./usr/local/lib64/swe/libswe.so.2.04
you must copy or move this to /usr/local/lib64/swe/libswe.so.2.04 by hand.
If your architecture is not 64-bit Linux, you may have to create libswe.so.2.04
from the Swisseph C source distribution yourself, and then place it properly.
Then you must tell the Linux system how to find the dynamic library.
On Redhat RHEL 6 this goes like this:
as root, cd /etc/ld.so.conf.d
edit or create a file swisseph.conf with this line in it
/usr/local/lib64/swe
now run as root: ldconfig
afterwards, you should find in /usr/local/lib64/swe/
something like this
lrwxrwxrwx. 1 root root 14 Mar 23 17:44 libswe.so -> libswe.so.2.04
lrwxrwxrwx. 1 root root 14 Mar 23 17:43 libswe.so.1 -> libswe.so.2.04
-rwxr-xr-x. 1 root root 847686 Mar 23 17:10 libswe.so.2.04
Now, as normal user, go back to
cd SwissEph-2.04
make
make test
as root, (unless you can write in the install target directories)
make install

2833
perl/SwissEph.xs Normal file

File diff suppressed because it is too large Load Diff

5
perl/fix_selinux Executable file
View File

@ -0,0 +1,5 @@
/usr/sbin/semanage fcontext -a -t textrel_shlib_t /usr/lib/perl5/site_perl/5.8.8/i386-linux-thread-multi/auto/SwissEph/SwissEph.so
/sbin/restorecon -v /usr/lib/perl5/site_perl/5.8.8/i386-linux-thread-multi/auto/SwissEph/SwissEph.so
/usr/sbin/semanage fcontext -a -t textrel_shlib_t /usr/lib/libswe.so.1.76.00
/sbin/restorecon -v /usr/lib/libswe.so.1.76.00

1725
perl/lib/SwissEph.pm Normal file

File diff suppressed because it is too large Load Diff

4896
perl/ppport.h Normal file

File diff suppressed because it is too large Load Diff

341
perl/sweodef.h Normal file
View File

@ -0,0 +1,341 @@
/************************************************************
definitions and constants for all Swiss Ephemeris source files,
only required for compiling the libraries, not for the external
interface of the libraries.
The definitions are a subset of Astrodienst's ourdef.h content
and must be kept compatible. Everything not used in SwissEph
has been deleted.
Does auto-detection of MSDOS (TURBO_C or MS_C), HPUNIX, Linux.
Must be extended for more portability; there should be a way
to detect byte order and file system type.
************************************************************/
/* 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.
*/
#ifndef _OURDEF_INCLUDED /* ourdef.h is a superset of sweodef.h */
#ifndef _SWEODEF_INCLUDED /* allow multiple #includes */
#define _SWEODEF_INCLUDED
# define MY_TRUE 1 /* for use in other defines, before TRUE is defined */
# define MY_FALSE 0 /* for use in other defines, before TRUE is defined */
/* TLS support
*
* Sun Studio C/C++, IBM XL C/C++, GNU C and Intel C/C++ (Linux systems) -> __thread
* Borland, VC++ -> __declspec(thread)
*/
#if !defined(TLSOFF) && !defined( __APPLE__ ) && !defined(WIN32) && !defined(DOS32)
#if defined( __GNUC__ )
#define TLS __thread
#else
#define TLS __declspec(thread)
#endif
#else
#define TLS
#endif
#ifdef _WIN32 /* Microsoft VC 5.0 does not define MSDOS anymore */
# undef MSDOS
# define MSDOS MY_TRUE
#include <wtypes.h>
#include <objbase.h>
#include <wincon.h>
#include <winbase.h>
#include <io.h>
#include <windows.h>
# define sleep(x) Sleep((x) * 1000)
#endif
#ifdef _MSC_VER
# define MS_VC
#endif
#ifdef WIN32 /* Microsoft VC 5.0 does not define MSDOS anymore */
# define MSDOS MY_TRUE
#endif
#ifdef MSDOS /* already defined by some DOS compilers */
# undef MSDOS
# define MSDOS MY_TRUE
#endif
#ifdef __TURBOC__ /* defined by turboc */
# ifndef MSDOS
# define MSDOS MY_TRUE
# endif
# define TURBO_C
#endif
#ifdef __SC__ /* defined by Symantec C */
# ifndef MSDOS
# define MSDOS MY_TRUE
# endif
# define SYMANTEC_C
#endif
#ifdef __WATCOMC__ /* defined by WatcomC */
# ifndef MSDOS
# define MSDOS MY_TRUE
# endif
# define WATCOMC
#endif
#ifdef __MWERKS__ /* defined on Macintosh CodeWarrior */
# if macintosh && powerc
# define MACOS MY_TRUE /* let it undefined otherwise */
# define MSDOS MY_FALSE /* in case one above fired falsely */
# endif
#endif
#ifdef MSDOS
# define HPUNIX MY_FALSE
# define INTEL_BYTE_ORDER 1
# ifndef TURBO_C
# define MS_C /* assume Microsoft C compiler */
# endif
# define UNIX_FS MY_FALSE
#else
# ifdef MACOS
# define HPUNIX MY_FALSE
# define UNIX_FS MY_FALSE
# else
# define MSDOS MY_FALSE
# define HPUNIX MY_TRUE
# ifndef _HPUX_SOURCE
# define _HPUX_SOURCE
# endif
# define UNIX_FS MY_TRUE
# endif
#endif
#include <math.h>
#include <stdlib.h>
#ifndef FILE
# include <stdio.h>
#endif
#if HPUNIX
# include <unistd.h>
#endif
/*
* if we have 16-bit ints, we define INT_16; we will need %ld to printf an int32
* if we have 64-bit long, we define LONG_64
* If none is defined, we have int = long = 32 bit, and use %d to printf an int32
*/
#include <limits.h>
#if INT_MAX < 40000
# define INT_16
#else
# if LONG_MAX > INT_MAX
# define LONG_64
# endif
#endif
#ifdef BYTE_ORDER
#ifdef LITTLE_ENDIAN
# if BYTE_ORDER == LITTLE_ENDIAN
# define INTEL_BYTE_ORDER
# endif
#endif
#endif
#ifdef INT_16
typedef long int32;
typedef unsigned long uint32;
typedef int int16;
typedef double REAL8; /* real with at least 64 bit precision */
typedef long INT4; /* signed integer with at least 32 bit precision */
typedef unsigned long UINT4;
/* unsigned integer with at least 32 bit precision */
typedef int AS_BOOL;
typedef unsigned int UINT2; /* unsigned 16 bits */
# define ABS4 labs /* abs function for long */
#else
typedef int int32;
typedef long long int64;
typedef unsigned int uint32;
typedef short int16;
typedef double REAL8; /* real with at least 64 bit precision */
typedef int INT4; /* signed integer with at least 32 bit precision */
typedef unsigned int UINT4;
/* unsigned integer with at least 32 bit precision */
typedef int AS_BOOL;
typedef unsigned short UINT2; /* unsigned 16 bits */
# define ABS4 abs /* abs function for long */
#endif
#if MSDOS
# ifdef TURBO_C
# include <alloc.h> /* MSC needs malloc ! */
# else
# include <malloc.h>
# endif
# define SIGALRM SIGINT
#endif
#ifndef TRUE
# define TRUE 1
# define FALSE 0
#endif
#ifndef OK
# define OK (0)
# define ERR (-1)
#endif
/* hack because UCHAR is already used by mingw gcc */
#ifdef __GNUC__
#ifdef _WIN32
#define UCHAR SWE_UCHAR
#endif
#endif
typedef unsigned char UCHAR;
#define UCP (UCHAR*)
#define SCP (char*)
# define ODEGREE_STRING "°" /* degree as string, utf8 encoding */
#ifndef HUGE
# define HUGE 1.7E+308 /* biggest value for REAL8 */
#endif
#ifndef M_PI
# define M_PI 3.14159265358979323846
#endif
/* #define forward static obsolete */
#define AS_MAXCH 256 /* used for string declarations, allowing 255 char+\0 */
/*
#define DEGTORAD 0.0174532925199433
#define RADTODEG 57.2957795130823
*/
#define RADTODEG (180.0 / M_PI)
#define DEGTORAD (M_PI / 180.0)
typedef int32 centisec; /* centiseconds used for angles and times */
#define CS (centisec) /* use for casting */
#define CSEC centisec /* use for typing */
#define DEG 360000 /* degree expressed in centiseconds */
#define DEG7_30 (2700000) /* 7.5 degrees */
#define DEG15 (15 * DEG)
#define DEG24 (24 * DEG)
#define DEG30 (30 * DEG)
#define DEG60 (60 * DEG)
#define DEG90 (90 * DEG)
#define DEG120 (120 * DEG)
#define DEG150 (150 * DEG)
#define DEG180 (180 * DEG)
#define DEG270 (270 * DEG)
#define DEG360 (360 * DEG)
/* #define CSTORAD 4.84813681109536E-08 centisec to rad: pi / 180 /3600/100 */
/* #define RADTOCS 2.06264806247096E+07 rad to centisec 180*3600*100/pi */
#define CSTORAD (DEGTORAD / 360000.0)
#define RADTOCS (RADTODEG * 360000.0)
#define CS2DEG (1.0/360000.0) /* centisec to degree */
/* control strings for fopen() */
#if UNIX_FS
# define BFILE_R_ACCESS "r" /* open binary file for reading */
# define BFILE_RW_ACCESS "r+" /* open binary file for writing and reading */
# define BFILE_W_CREATE "w" /* create/open binary file for write*/
# define BFILE_A_ACCESS "a+" /* create/open binary file for append*/
# define FILE_R_ACCESS "r" /* open text file for reading */
# define FILE_RW_ACCESS "r+" /* open text file for writing and reading */
# define FILE_W_CREATE "w" /* create/open text file for write*/
# define FILE_A_ACCESS "a+" /* create/open text file for append*/
# define O_BINARY 0 /* for open(), not defined in Unix */
# define OPEN_MODE 0666 /* default file creation mode */
# define DIR_GLUE "/" /* glue string for directory/file */
# define PATH_SEPARATOR ";:" /* semicolon or colon may be used */
#else
# define BFILE_R_ACCESS "rb" /* open binary file for reading */
# define BFILE_RW_ACCESS "r+b" /* open binary file for writing and reading */
# define BFILE_W_CREATE "wb" /* create/open binary file for write*/
# define BFILE_A_ACCESS "a+b" /* create/open binary file for append*/
# define PATH_SEPARATOR ";" /* semicolon as PATH separator */
# define OPEN_MODE 0666 /* default file creation mode */
# ifdef MACOS
# define FILE_R_ACCESS "r" /* open text file for reading */
# define FILE_RW_ACCESS "r+" /* open text file for writing and reading */
# define FILE_W_CREATE "w" /* create/open text file for write*/
# define FILE_A_ACCESS "a+" /* create/open text file for append*/
# define DIR_GLUE ":" /* glue string for directory/file */
# else
# define FILE_R_ACCESS "rt" /* open text file for reading */
# define FILE_RW_ACCESS "r+t" /* open text file for writing and reading */
# define FILE_W_CREATE "wt" /* create/open text file for write*/
# define FILE_A_ACCESS "a+t" /* create/open text file for append*/
/* attention, all backslashes for msdos directry names must be written as \\,
because it is the C escape character */
# define DIR_GLUE "\\" /* glue string for directory/file */
# endif
#endif
#include <string.h>
#include <ctype.h>
#endif /* _SWEODEF_INCLUDED */
#endif /* _OURDEF_INCLUDED */

1015
perl/swephexp.h Normal file

File diff suppressed because it is too large Load Diff

567
perl/t/SwissEph.t Normal file
View File

@ -0,0 +1,567 @@
# Before `make install' is performed this script should be runnable with
# `make test'. After `make install' it should work as `perl SwissEph.t'
use Test::More tests => 248;
BEGIN { use_ok("SwissEph") };
#------------------------------------------------------------------------
# Sun for 2451544.5 ET : Longitude should be 279.858461
#------------------------------------------------------------------------
my $ref;
my $i;
my @xx;
my $serr="";
#------------------------------------------------------------------------
# Coordinate transformation
#------------------------------------------------------------------------
# ecliptic -> equator
my @rade = SwissEph::swe_cotrans([80, 5, 0], -23);
is(round_4($rade[0]), 78.7418, "swe_cotrans([80, 5, 0], -23) -> ra");
is(round_4($rade[1]), 27.6169, "swe_cotrans([80, 5, 0], -23) -> de");
# with speed
@rade = SwissEph::swe_cotrans_sp([80, 5, 0, 1, 0, 0], -23);
is(round_4($rade[0]), 78.7418, "swe_cotrans([80, 5, 0], -23) -> ra");
is(round_4($rade[1]), 27.6169, "swe_cotrans([80, 5, 0], -23) -> de");
is(round_4($rade[3]), 1.1210, "swe_cotrans([80, 5, 0], -23) -> ra_speed");
is(round_4($rade[4]), 0.0763, "swe_cotrans([80, 5, 0], -23) -> de_speed");
#------------------------------------------------------------------------
# Function Delta T - values may change, but the seconds for 2000-1-1 are sure
#------------------------------------------------------------------------
is( round_4(SwissEph::swe_deltat( 2451545 )* 86400), 63.8289, "swe_deltat(1 jan 2000)");
$ref = SwissEph::swe_deltat_ex(2451545, SwissEph::SEFLG_MOSEPH);
is(round_4($ref->{dt} *86400.0), 63.8289, "swe_deltat_ex(1 jan 2000)");
#------------------------------------------------------------------------
# Sidereal Time
#------------------------------------------------------------------------
is(round_4( SwissEph::swe_sidtime( 2451544.5 )), 6.6643, "swe_sidtime(1 jan 2000)");
is(round_4( SwissEph::swe_sidtime0( 2451544.5, 23.5, 0)), 6.6645, "swe_sidtime0(1 jan 2000,23.5,0)");
$ref = SwissEph::swe_time_equ(2415020.5);
is($ref->{retval}, 0, "swe_time_equ(1 jan 1900)->retval");
if (exists($ref->{serr})) {
print STDERR "swe_time_equ(1 jan 1900)->serr = $ref->{serr}\n";
}
is(round_4($ref->{time_equ}), -0.0023, "swe_time_equ(1 jan 1900)->time_equ");
$ref = SwissEph::swe_lat_to_lmt(-146780.0, 82.2);
if (exists($ref->{serr})) {
print STDERR "swe_time_equ(1 jan 1900)->serr = $ref->{serr}\n";
}
is(round_4($ref->{tjd_lmt}), -146779.9898, "swe_lat_to_lmt(-146780.0, 82.2)->tjd_lmt");
$ref = SwissEph::swe_lmt_to_lat(-146779.9898, 82.2);
if (exists($ref->{serr})) {
print STDERR "swe_time_equ(1 jan 1900)->serr = $ref->{serr}\n";
}
is(round_4($ref->{tjd_lat}), -146779.9998, "swe_lmt_to_lat(-146779.986742931, 82.2)->tjd_lat");
#------------------------------------------------------------------------
# Calendar conversion
#------------------------------------------------------------------------
is( SwissEph::swe_julday(2000,1,1,0.), 2451544.5, "swe_julday() 1 jan 2000");
is( SwissEph::swe_julday(2000,1,1,0.,1), 2451544.5, "swe_julday() 1 jan 2000, explicit greg flag");
is( SwissEph::swe_julday(2000,1,1,0.,0), 2451557.5, "swe_julday() 1 jan 2000, explicit jul flag");
$ref = SwissEph::swe_revjul(2436723.5888888889,1);
is($ref->{iyar}, 1959, "swe_revjul(2436723.5888888889), year");
is($ref->{imon}, 6, "swe_revjul(2436723.5888888889), month");
is($ref->{iday}, 4, "swe_revjul(2436723.5888888889), day");
is($ref->{ihou}, 2, "swe_revjul(2436723.5888888889), hour");
is($ref->{imin}, 8, "swe_revjul(2436723.5888888889), minute");
is($ref->{isec}, 0, "swe_revjul(2436723.5888888889), second");
is(round_4($ref->{dhou}), 2.1333, "swe_revjul(2436723.5888888889), hour/decimal");
#------------------------------------------------------------------------
# Wrapper for swe_houses
#------------------------------------------------------------------------
$ref = SwissEph::swe_houses( 2451544.5, 52., 7., "P");
is( round_4( $ref->{cusps}->[1] ), 191.8825, "swe_houses()->{1}");
is( round_4( $ref->{cusps}->[2] ), 217.0327, "swe_houses()->{2}");
is( round_4( $ref->{cusps}->[3] ), 248.4837, "swe_houses()->{3}");
is( round_4( $ref->{cusps}->[4] ), 285.6359, "swe_houses()->{4}");
is( round_4( $ref->{cusps}->[5] ), 321.0151, "swe_houses()->{5}");
is( round_4( $ref->{cusps}->[6] ), 349.6374, "swe_houses()->{6}");
is( round_4( $ref->{cusps}->[7] ), 11.8825, "swe_houses()->{7}");
is( round_4( $ref->{cusps}->[8] ), 37.0327, "swe_houses()->{8}");
is( round_4( $ref->{cusps}->[9] ), 68.4837, "swe_houses()->{9}");
is( round_4( $ref->{cusps}->[10] ), 105.6359, "swe_houses()->{10}");
is( round_4( $ref->{cusps}->[11] ), 141.0151, "swe_houses()->{11}");
is( round_4( $ref->{cusps}->[12] ), 169.6374, "swe_houses()->{12}");
is( round_4( $ref->{asc} ), 191.8825, "swe_houses()->{asc}");
is( round_4( $ref->{mc} ), 105.6359, "swe_houses()->{mc}");
is( round_4( $ref->{armc} ), 106.9642, "swe_houses()->{armc}");
is( round_4( $ref->{vertex} ), 27.2378, "swe_houses()->{vertex}");
is( round_4( $ref->{equasc} ), 198.3910, "swe_houses()->{equasc}");
is( round_4( $ref->{coasc1} ), 218.3740, "swe_houses()->{coasc1}");
is( round_4( $ref->{coasc2} ), 193.7952, "swe_houses()->{coasc2}");
is( round_4( $ref->{polasc} ), 38.3740, "swe_houses()->{polasc}");
$ref = SwissEph::swe_houses_armc( 106.96424825, 52., 23.4376796111111, "P");
is( round_4( $ref->{cusps}->[1] ), 191.8825, "swe_houses_armc()->{1}");
is( round_4( $ref->{cusps}->[2] ), 217.0327, "swe_houses_armc()->{2}");
is( round_4( $ref->{cusps}->[3] ), 248.4837, "swe_houses_armc()->{3}");
is( round_4( $ref->{cusps}->[4] ), 285.6359, "swe_houses_armc()->{4}");
is( round_4( $ref->{cusps}->[5] ), 321.0151, "swe_houses_armc()->{5}");
is( round_4( $ref->{cusps}->[6] ), 349.6374, "swe_houses_armc()->{6}");
is( round_4( $ref->{cusps}->[7] ), 11.8825, "swe_houses_armc()->{7}");
is( round_4( $ref->{cusps}->[8] ), 37.0327, "swe_houses_armc()->{8}");
is( round_4( $ref->{cusps}->[9] ), 68.4837, "swe_houses_armc()->{9}");
is( round_4( $ref->{cusps}->[10] ), 105.6359, "swe_houses_armc()->{10}");
is( round_4( $ref->{cusps}->[11] ), 141.0151, "swe_houses_armc()->{11}");
is( round_4( $ref->{cusps}->[12] ), 169.6374, "swe_houses_armc()->{12}");
is( round_4( $ref->{asc} ), 191.8825, "swe_houses_armc()->{asc}");
is( round_4( $ref->{mc} ), 105.6359, "swe_houses_armc()->{mc}");
is( round_4( $ref->{armc} ), 106.9642, "swe_houses_armc()->{armc}");
is( round_4( $ref->{vertex} ), 27.2378, "swe_houses_armc()->{vertex}");
is( round_4( $ref->{equasc} ), 198.3910, "swe_houses_armc()->{equasc}");
is( round_4( $ref->{coasc1} ), 218.3740, "swe_houses_armc()->{coasc1}");
is( round_4( $ref->{coasc2} ), 193.7952, "swe_houses_armc()->{coasc2}");
is( round_4( $ref->{polasc} ), 38.3740, "swe_houses_armc()->{polasc}");
$ref = SwissEph::swe_houses_ex( 2451544.5, 0, 52., 7., "P");
is( round_4( $ref->{cusps}->[1] ), 191.8825, "swe_houses_ex()->{1}");
is( round_4( $ref->{cusps}->[2] ), 217.0327, "swe_houses_ex()->{2}");
is( round_4( $ref->{cusps}->[3] ), 248.4837, "swe_houses_ex()->{3}");
is( round_4( $ref->{cusps}->[4] ), 285.6359, "swe_houses_ex()->{4}");
is( round_4( $ref->{cusps}->[5] ), 321.0151, "swe_houses_ex()->{5}");
is( round_4( $ref->{cusps}->[6] ), 349.6374, "swe_houses_ex()->{6}");
is( round_4( $ref->{cusps}->[7] ), 11.8825, "swe_houses_ex()->{7}");
is( round_4( $ref->{cusps}->[8] ), 37.0327, "swe_houses_ex()->{8}");
is( round_4( $ref->{cusps}->[9] ), 68.4837, "swe_houses_ex()->{9}");
is( round_4( $ref->{cusps}->[10] ), 105.6359, "swe_houses_ex()->{10}");
is( round_4( $ref->{cusps}->[11] ), 141.0151, "swe_houses_ex()->{11}");
is( round_4( $ref->{cusps}->[12] ), 169.6374, "swe_houses_ex()->{12}");
is( round_4( $ref->{asc} ), 191.8825, "swe_houses_ex()->{asc}");
is( round_4( $ref->{mc} ), 105.6359, "swe_houses_ex()->{mc}");
is( round_4( $ref->{armc} ), 106.9642, "swe_houses_ex()->{armc}");
is( round_4( $ref->{vertex} ), 27.2378, "swe_houses_ex()->{vertex}");
is( round_4( $ref->{equasc} ), 198.3910, "swe_houses_ex()->{equasc}");
is( round_4( $ref->{coasc1} ), 218.3740, "swe_houses_ex()->{coasc1}");
is( round_4( $ref->{coasc2} ), 193.7952, "swe_houses_ex()->{coasc2}");
is( round_4( $ref->{polasc} ), 38.3740, "swe_houses_ex()->{polasc}");
$ref = SwissEph::swe_house_pos(290,47,23.5,"P",72,0);
is($ref->{retval}, 0, "swe_house_pos()->retval");
if (exists($ref->{serr})) {
print STDERR "swe_house_pos()->serr = $ref->{serr}\n";
}
is($ref->{ihno}, 2, "swe_house_pos()->ihno");
is(round_4($ref->{dhpos}), 2.1459, "swe_house_pos()->dhpos");
is(round_4($ref->{dhpos_deg}), 34.3767, "swe_house_pos()->dhpos_deg");
$ref = SwissEph::swe_gauquelin_sector(2436723.5888888889,SwissEph::SE_MOON,"",0,0,[8.6,47.35,400],1013,15);
is($ref->{retval}, 0, "swe_gauquelin_sector()->retval");
if (exists($ref->{serr})) {
print STDERR "swe_gauquelin_sector(1 jan 1900)->serr = $ref->{serr}\n";
}
is(round_4($ref->{dsector}), 36.1724, "swe_gauquelin_sector()->dsector");
#------------------------------------------------------------------------
# Wrapper for swe_calc
#------------------------------------------------------------------------
$ref = SwissEph::swe_calc(2415020.5, 3, 260);
is($ref->{retval}, 260, "swe_calc(1900, Venus)->retval");
if (exists($ref->{serr})) {
print STDERR "swe_calc(1900, Venus)->serr = $ref->{serr}\n";
}
is(round_4($ref->{xx}->[0]), 306.3745, "swe_calc(1900, Venus)->xx[0]");
is(round_4($ref->{xx}->[1]), -1.6830, "swe_calc(1900, Venus)->xx[1]");
is(round_4($ref->{xx}->[2]), 1.4646, "swe_calc(1900, Venus)->xx[2]");
is(round_4($ref->{xx}->[3]), 1.2435, "swe_calc(1900, Venus)->xx[3]");
my $dt = SwissEph::swe_deltat(2415020.5);
$ref = SwissEph::swe_calc_ut(2415020.5-$dt, 3, 260);
is($ref->{retval}, 260, "swe_calc_ut(1900, Venus)->retval");
if (exists($ref->{serr})) {
print STDERR "swe_calc_ut(1900, Venus)->serr = $ref->{serr}\n";
}
is(round_4($ref->{xx}->[0]), 306.3745, "swe_calc_ut(1900, Venus)->xx[0]");
is(round_4($ref->{xx}->[1]), -1.6830, "swe_calc_ut(1900, Venus)->xx[1]");
is(round_4($ref->{xx}->[2]), 1.4646, "swe_calc_ut(1900, Venus)->xx[2]");
is(round_4($ref->{xx}->[3]), 1.2435, "swe_calc_ut(1900, Venus)->xx[3]");
#print STDERR "swe_calc_ut $ref->{xx}->[0] $ref->{xx}->[1]\n";
#------------------------------------------------------------------------
# Wrapper for swe_fixstar
#------------------------------------------------------------------------
$ref = SwissEph::swe_fixstar_mag("alcyone");
is($ref->{retval}, 0, "swe_fixstar_mag(Alcyone)->retval");
if (exists($ref->{serr})) {
print STDERR "swe_fixstar_mag(Alcyone)->serr = $ref->{serr}\n";
}
is($ref->{starname}, "Alcyone,etTau", "swe_fixstar_mag(Alcyone)->starname");
is(round_4($ref->{dmag}), 2.87, "swe_fixstar_mag(Alcyone)->dmag");
$ref = SwissEph::swe_fixstar("alcyone", 2415020.5, 4);
is($ref->{retval}, 4, "swe_fixstar(1900, Alcyone)->retval");
if (exists($ref->{serr})) {
print STDERR "swe_fixstar(1900, Alcyone)->serr = $ref->{serr}\n";
}
is($ref->{starname}, "Alcyone,etTau", "swe_fixstar(1900, Alcyone)->starname");
is(round_4($ref->{xx}->[0]), 58.6052, "swe_fixstar(1900, Alcyone)->xx[0]");
is(round_4($ref->{xx}->[2]), 25496153.13, "swe_fixstar(1900, Alcyone)->xx[2]");
$ref = SwissEph::swe_fixstar_ut("alcyone", 2415020.5-SwissEph::swe_deltat(2415020.5), 4);
is($ref->{retval}, 4, "swe_fixstar_ut(1900, Alcyone)->retval");
if (exists($ref->{serr})) {
print STDERR "swe_fixstar_ut(1900, Alcyone)->serr = $ref->{serr}\n";
}
is($ref->{starname}, "Alcyone,etTau", "swe_fixstar_ut(1900, Alcyone)->starname");
is(round_4($ref->{xx}->[0]), 58.6052, "swe_fixstar_ut(1900, Alcyone)->xx[0]");
is(round_4($ref->{xx}->[2]), 25496153.13, "swe_fixstar_ut(1900, Alcyone)->xx[2]");
$ref = SwissEph::swe_fixstar2("alcyone", 2415020.5, 4);
is($ref->{retval}, 4, "swe_fixstar2(1900, Alcyone)->retval");
if (exists($ref->{serr})) {
print STDERR "swe_fixstar2(1900, Alcyone)->serr = $ref->{serr}\n";
}
is($ref->{starname}, "Alcyone,etTau", "swe_fixstar2(1900, Alcyone)->starname");
is(round_4($ref->{xx}->[0]), 58.6052, "swe_fixstar2(1900, Alcyone)->xx[0]");
is(round_4($ref->{xx}->[2]), 25496153.13, "swe_fixstar2(1900, Alcyone)->xx[2]");
$ref = SwissEph::swe_fixstar2_ut("alcyone", 2415020.5-SwissEph::swe_deltat(2415020.5), 4);
is($ref->{retval}, 4, "swe_fixstar2_ut(1900, Alcyone)->retval");
if (exists($ref->{serr})) {
print STDERR "swe_fixstar2_ut(1900, Alcyone)->serr = $ref->{serr}\n";
}
is($ref->{starname}, "Alcyone,etTau", "swe_fixstar2_ut(1900, Alcyone)->starname");
is(round_4($ref->{xx}->[0]), 58.6052, "swe_fixstar2_ut(1900, Alcyone)->xx[0]");
is(round_4($ref->{xx}->[2]), 25496153.13, "swe_fixstar2_ut(1900, Alcyone)->xx[2]");
$ref = SwissEph::swe_fixstar2_mag("alcyone");
is($ref->{retval}, 0, "swe_fixstar2_mag(Alcyone)->retval");
if (exists($ref->{serr})) {
print STDERR "swe_fixstar2_mag(Alcyone)->serr = $ref->{serr}\n";
}
is($ref->{starname}, "Alcyone,etTau", "swe_fixstar2_mag(Alcyone)->starname");
is(round_4($ref->{dmag}), 2.87, "swe_fixstar2_mag(Alcyone)->dmag");
$ref = SwissEph::swe_pheno(2415020.5, 3, 4);
is($ref->{retval}, 4, "swe_pheno(1900, Venus)->retval");
if (exists($ref->{serr})) {
print STDERR "swe_pheno(1900, Venus)->serr = $ref->{serr}\n";
}
is(round_4($ref->{phase_angle}), 36.7449, "swe_pheno(1900, Venus)->phase_angle");
is(round_4($ref->{phase}), 0.9007, "swe_pheno(1900, Venus)->phase");
is(round_4($ref->{elongation}), 26.2712, "swe_pheno(1900, Venus)->elongation");
is(round_4($ref->{disc_diameter}), 0.0032, "swe_pheno(1900, Venus)->disc_diameter");
is(round_4($ref->{magnitude}), -3.9102, "swe_pheno(1900, Venus)->magnitude");
is(round_4($ref->{hor_parallax}), 0, "swe_pheno(1900, Venus)->hor_parallax");
is(round_4($ref->{attr}->[0]), 36.7449, "swe_pheno(1900, Venus)->attr[0]");
$ref = SwissEph::swe_pheno_ut(2415020.5 - SwissEph::swe_deltat(2415020.5), 3, 4);
is($ref->{retval}, 4, "swe_pheno_ut(1900, Venus)->retval");
if (exists($ref->{serr})) {
print STDERR "swe_pheno_ut(1900, Venus)->serr = $ref->{serr}\n";
}
is(round_4($ref->{phase_angle}), 36.7449, "swe_pheno_ut(1900, Venus)->phase_angle");
is(round_4($ref->{phase}), 0.9007, "swe_pheno_ut(1900, Venus)->phase");
is(round_4($ref->{elongation}), 26.2712, "swe_pheno_ut(1900, Venus)->elongation");
is(round_4($ref->{disc_diameter}), 0.0032, "swe_pheno_ut(1900, Venus)->disc_diameter");
is(round_4($ref->{magnitude}), -3.9102, "swe_pheno_ut(1900, Venus)->magnitude");
is(round_4($ref->{hor_parallax}), 0, "swe_pheno_ut(1900, Venus)->hor_parallax");
is(round_4($ref->{attr}->[0]), 36.7449, "swe_pheno_ut(1900, Venus)->attr[0]");
#------------------------------------------------------------------------
# ayanamsha
#------------------------------------------------------------------------
is(round_4(SwissEph::swe_get_ayanamsa(2451544.5)), 24.7403, "swe_get_ayanamsa(2000)");
is(round_4(SwissEph::swe_get_ayanamsa_ut(2451544.5-SwissEph::swe_deltat(2451544.5))), 24.7403, "swe_get_ayanamsa(2000)");
$ref = SwissEph::swe_get_ayanamsa_ex(2451544.5, SwissEph::SEFLG_MOSEPH);
is(round_4($ref->{daya}),24.7364, "swe_get_ayanamsa_ex(2000)");
$ref = SwissEph::swe_get_ayanamsa_ex_ut(2451544.5, SwissEph::SEFLG_MOSEPH);
is(round_4($ref->{daya}),24.7364, "swe_get_ayanamsa_ex(2000)");
is(SwissEph::swe_get_ayanamsa_name(3), "Raman", "swe_get_ayanamsa(Raman)");
is(SwissEph::swe_get_planet_name(SwissEph::SE_VENUS), "Venus", "swe_get_planet_name(Venus)");
#------------------------------------------------------------------------
# refraction, azimuth, and altitude, rise and transit
#------------------------------------------------------------------------
$ref = SwissEph::swe_refrac_extended(0.1,400,1013,15,0.02,1);
is($ref->{retval}, 0, "swe_refrac_extended()->retval");
if (exists($ref->{serr})) {
print STDERR "swe_refrac_extended()->serr = $ref->{serr}\n";
}
is(round_4($ref->{alt_true}), -0.4391, "swe_refrac_extended()->alt_true");
is(round_4($ref->{alt_apparent}), 0.1, "swe_refrac_extended()->alt_apparent");
is(round_4($ref->{refraction}), 0.5392, "swe_refrac_extended()->refraction");
is(round_4($ref->{dip}), -0.5238, "swe_refrac_extended()->dip");
is(round_4($ref->{dret}->[0]), -0.4391, "swe_refrac_extended()->dret[0]");
is(round_4(SwissEph::swe_refrac(0.6,1013,15,1)), 0.1441, "swe_refrac()");
@xx = SwissEph::swe_azalt(2436723.588888888,0,[8.55,47.35,0],1013,15,[252,0,1]);
is(round_4($xx[0]), 38.9512, "swe_azalt()->azimuth");
is(round_4($xx[1]), 10.8971, "swe_azalt()->alt_true");
is(round_4($xx[2]), 10.9771, "swe_azalt()->alt_apparent");
@xx = SwissEph::swe_azalt_rev(2436723.588888888,0,[8.55,47.35,0],\@xx);
is(round_4($xx[0]), 252, "swe_azalt_rev()->ecl_lon");
is(round_4($xx[1]), 0, "swe_azalt_rev()->ecl_lat");
$ref = SwissEph::swe_rise_trans(2415020.5,3,"", 4,0,[8.55,47.23,400],1013,15);
is($ref->{retval}, 0, "swe_rise_trans(1900, Venus)->retval");
if (exists($ref->{serr})) {
print STDERR "swe_rise_trans(1900, Venus)->serr = $ref->{serr}\n";
}
is(round_4($ref->{dret}), 2415020.8696, "swe_rise_trans(1900, Venus)->dret");
#------------------------------------------------------------------------
# eclipses and occultations
#------------------------------------------------------------------------
$ref = SwissEph::swe_lun_eclipse_how(2454517.643069,4,[278,0,0]);
is($ref->{retval}, 4, "swe_lun_eclipse_how(2008)->retval");
if (exists($ref->{serr})) {
print STDERR "swe_lun_eclipse_how(2008)->serr = $ref->{serr}\n";
}
is(round_4($ref->{mag_umbral}), 1.1059, "swe_lun_eclipse_how(2008)->mag_umbral");
is(round_4($ref->{mag_penumbral}), 2.1450, "swe_lun_eclipse_how(2008)->mag_penumbral");
is(round_4($ref->{attr}->[0]), 1.1059, "swe_lun_eclipse_how(2008)->attr[0]");
$ref = SwissEph::swe_sol_eclipse_how(2454503.663212,4,[-150.270493,-67.547072,0]);
is($ref->{retval}, 137, "swe_sol_eclipse_how(2008)->retval");
if (exists($ref->{serr})) {
print STDERR "swe_sol_eclipse_how(2008)->serr = $ref->{serr}\n";
}
is(round_4($ref->{disc_ratio}), 0.9658, "swe_sol_eclipse_how(2008)->disc_ratio");
is(round_4($ref->{fraction_diameter}), 0.9808, "swe_sol_eclipse_how(2008)->fraction_diameter");
is(round_4($ref->{fraction_disc}), 0.9327, "swe_sol_eclipse_how(2008)->fraction_disc");
is(round_4($ref->{core_shadow_km}), 123.5327, "swe_sol_eclipse_how(2008)->core_shadow_km");
is(round_4($ref->{sun_azimuth}), 88.5702, "swe_sol_eclipse_how(2008)->sun_azimuth");
is(round_4($ref->{sun_alt_true}), 16.2305, "swe_sol_eclipse_how(2008)->sun_alt_true");
is(round_4($ref->{separation_angle}), 0.0011, "swe_sol_eclipse_how(2008)->separation_angle");
$ref = SwissEph::swe_sol_eclipse_where(2454503.663212,4);
is($ref->{retval}, 9, "swe_sol_eclipse_where(2008)->retval");
if (exists($ref->{serr})) {
print STDERR "swe_sol_eclipse_where(2008)->serr = $ref->{serr}\n";
}
is(round_4($ref->{geopos}->[0]), -150.3452, "swe_sol_eclipse_where(2008)->geopos[0]");
is(round_4($ref->{geopos}->[1]), -67.5477, "swe_sol_eclipse_where(2008)->geopos[1]");
is(round_4($ref->{disc_ratio}), 0.9658, "swe_sol_eclipse_where(2008)->disc_ratio");
is(round_4($ref->{fraction_diameter}), 0.9809, "swe_sol_eclipse_where(2008)->fraction_diameter");
is(round_4($ref->{fraction_disc}), 0.9327, "swe_sol_eclipse_where(2008)->fraction_disc");
is(round_4($ref->{core_shadow_km}), 123.5327, "swe_sol_eclipse_where(2008)->core_shadow_km");
is(round_4($ref->{sun_azimuth}), 88.6393, "swe_sol_eclipse_where(2008)->sun_azimuth");
is(round_4($ref->{sun_alt_true}), 16.2591, "swe_sol_eclipse_where(2008)->sun_alt_true");
is(round_4($ref->{separation_angle}), 0.0011, "swe_sol_eclipse_where(2008)->separation_angle");
$ref = SwissEph::swe_lun_occult_where(2454531.296945,SwissEph::SE_VENUS,"",4);
is($ref->{retval}, 5, "swe_lun_occult_where(2008)->retval");
if (exists($ref->{serr})) {
print STDERR "swe_lun_occult_where(2008)->serr = $ref->{serr}\n";
}
is(round_4($ref->{geopos}->[0]), -132.4299, "swe_lun_occult_where(2008)->geopos[0]");
is(round_4($ref->{geopos}->[1]), -3.2154, "swe_lun_occult_where(2008)->geopos[1]");
is(round_4($ref->{disc_ratio}), 172.5202, "swe_lun_occult_where(2008)->disc_ratio");
is(round_4($ref->{fraction_diameter}), 1, "swe_lun_occult_where(2008)->fraction_diameter");
is(round_4($ref->{fraction_disc}), 1, "swe_lun_occult_where(2008)->fraction_disc");
is(round_4($ref->{core_shadow_km}), -3461.9130, "swe_lun_occult_where(2008)->core_shadow_km");
is(round_4($ref->{body_azimuth}), 336.2908, "swe_lun_occult_where(2008)->body_azimuth");
is(round_4($ref->{body_alt_true}), 76.8443, "swe_lun_occult_where(2008)->body_alt_true");
is(round_4($ref->{separation_angle}), 0.0000, "swe_lun_occult_where(2008)->separation_angle");
$ref = SwissEph::swe_lun_eclipse_when(2454466.500000,4,0,0);
is($ref->{retval}, 4, "swe_lun_eclipse_when(2008)->retval");
if (exists($ref->{serr})) {
print STDERR "swe_lun_eclipse_when(2008)->serr = $ref->{serr}\n";
}
is(round_4($ref->{ecl_maximum}), 2454517.6431, "swe_lun_eclipse_when(2008)->ecl_maximum");
is(round_4($ref->{ecl_partial_begin}), 2454517.5717, "swe_lun_eclipse_when(2008)->partial_begin");
is(round_4($ref->{ecl_partial_end}), 2454517.7144, "swe_lun_eclipse_when(2008)->partial_end");
is(round_4($ref->{ecl_total_begin}), 2454517.6258, "swe_lun_eclipse_when(2008)->ecl_total_begin");
is(round_4($ref->{ecl_total_end}), 2454517.6603, "swe_lun_eclipse_when(2008)->ecl_total_end");
is(round_4($ref->{ecl_penumbral_begin}), 2454517.5254, "swe_lun_eclipse_when(2008)->ecl_penumbral_begin");
is(round_4($ref->{ecl_penumbral_end}), 2454517.7609, "swe_lun_eclipse_when(2008)->ecl_penumbral_end");
$ref = SwissEph::swe_sol_eclipse_when_glob(2454466.500000,4,0,0);
is($ref->{retval}, 9, "swe_sol_eclipse_when_glob(2008)->retval");
if (exists($ref->{serr})) {
print STDERR "swe_sol_eclipse_when_glob(2008)->serr = $ref->{serr}\n";
}
is(round_4($ref->{ecl_maximum}), 2454503.6632, "swe_sol_eclipse_when_glob(2008)->ecl_maximum");
is(round_4($ref->{ecl_local_noon}), 2454503.6312, "swe_sol_eclipse_when_glob(2008)->ecl_local_noon");
is(round_4($ref->{ecl_begin}), 2454503.5686, "swe_sol_eclipse_when_glob(2008)->ecl_begin");
is(round_4($ref->{ecl_end}), 2454503.7582, "swe_sol_eclipse_when_glob(2008)->ecl_end");
is(round_4($ref->{ecl_total_begin}), 2454503.6388, "swe_sol_eclipse_when_glob(2008)->ecl_total_begin");
is(round_4($ref->{ecl_total_end}), 2454503.6879, "swe_sol_eclipse_when_glob(2008)->ecl_total_end");
is(round_4($ref->{ecl_central_begin}), 2454503.6417, "swe_sol_eclipse_when_glob(2008)->ecl_central_begin");
is(round_4($ref->{ecl_central_end}), 2454503.6850, "swe_sol_eclipse_when_glob(2008)->ecl_central_end");
$ref = SwissEph::swe_lun_occult_when_glob(2454466.500000,SwissEph::SE_VENUS,"",4,0,0);
is($ref->{retval}, 5, "swe_lun_occult_when_glob(2008, Venus)->retval");
if (exists($ref->{serr})) {
print STDERR "swe_lun_occult_when_glob(2008 Venus)->serr = $ref->{serr}\n";
}
is(round_4($ref->{occ_maximum}), 2454531.2969, "swe_lun_occult_when_glob(2008 Venus)->occ_maximum");
is(round_4($ref->{occ_local_noon}), 2454531.3051, "swe_lun_occult_when_glob(2008 Venus)->occ_local_noon");
is(round_4($ref->{occ_begin}), 2454531.1986, "swe_lun_occult_when_glob(2008 Venus)->occ_begin");
is(round_4($ref->{occ_end}), 2454531.3951, "swe_lun_occult_when_glob(2008 Venus)->occ_end");
is(round_4($ref->{occ_total_begin}), 2454531.1989, "swe_lun_occult_when_glob(2008 Venus)->occ_total_begin");
is(round_4($ref->{occ_total_end}), 2454531.3948, "swe_lun_occult_when_glob(2008 Venus)->occ_total_end");
is(round_4($ref->{occ_central_begin}), 2454531.2206, "swe_lun_occult_when_glob(2008 Venus)->occ_central_begin");
is(round_4($ref->{occ_central_end}), 2454531.3731, "swe_lun_occult_when_glob(2008 Venus)->occ_central_end");
$ref = SwissEph::swe_sol_eclipse_when_loc(2454466.500000,4,0,[8.55,47.35,400]);
is($ref->{retval}, 5008, "swe_sol_eclipse_when_loc(2008)->retval");
if (exists($ref->{serr})) {
print STDERR "swe_sol_eclipse_when_loc(2008)->serr = $ref->{serr}\n";
}
is(round_4($ref->{ecl_maximum}), 2454679.8966, "swe_sol_eclipse_when_loc(2008)->ecl_maximum");
is(round_4($ref->{t1st_contact}), 2454679.8708, "swe_sol_eclipse_when_loc(2008)->1st_contact");
is(round_4($ref->{t2nd_contact}), 0, "swe_sol_eclipse_when_loc(2008)->2nd_contact");
is(round_4($ref->{t3rd_contact}), 0, "swe_sol_eclipse_when_loc(2008)->3rd_contact");
is(round_4($ref->{t4th_contact}), 2454679.9230, "swe_sol_eclipse_when_loc(2008)->4th_contact");
is(round_4($ref->{disc_ratio}), 1.0449, "swe_sol_eclipse_when_loc(2008)->disc_ratio");
is(round_4($ref->{fraction_diameter}), 0.1203, "swe_sol_eclipse_when_loc(2008)->fraction_diameter");
is(round_4($ref->{fraction_disc}), 0.0497, "swe_sol_eclipse_when_loc(2008)->fraction_disc");
is(round_4($ref->{core_shadow_km}), -120.6315, "swe_sol_eclipse_when_loc(2008)->core_shadow_km");
is(round_4($ref->{sun_azimuth}), 309.5993, "swe_sol_eclipse_when_loc(2008)->sun_azimuth");
is(round_4($ref->{sun_alt_true}), 51.51, "swe_sol_eclipse_when_loc(2008)->sun_alt_true");
is(round_4($ref->{separation_angle}), 0.4739, "swe_sol_eclipse_when_loc(2008)->separation_angle");
$ref = SwissEph::swe_lun_occult_when_loc(2454466.500000,SwissEph::SE_VENUS,"",4,0,[8.55,47.35,400]);
is($ref->{retval}, 8070, "swe_lun_occult_when_loc(2008, Venus)->retval");
if (exists($ref->{serr})) {
print STDERR "swe_lun_occult_when_loc(2008 Venus)->serr = $ref->{serr}\n";
}
is(round_4($ref->{occ_maximum}), 2454802.1986, "swe_lun_occult_when_loc(2008 Venus)->occ_maximum");
is(round_4($ref->{t1st_contact}), 2454802.1698, "swe_lun_occult_when_loc(2008 Venus)->1st_contact");
is(round_4($ref->{t2nd_contact}), 2454802.1705, "swe_lun_occult_when_loc(2008 Venus)->2nd_contact");
is(round_4($ref->{t3rd_contact}), 2454802.2253, "swe_lun_occult_when_loc(2008 Venus)->3rd_contact");
is(round_4($ref->{t4th_contact}), 2454802.2259, "swe_lun_occult_when_loc(2008 Venus)->4th_contact");
is(round_4($ref->{disc_ratio}), 106.8891, "swe_lun_occult_when_loc(2008 Venus)->disc_ratio");
is(round_4($ref->{fraction_diameter}), 1, "swe_lun_occult_when_loc(2008 Venus)->fraction_diameter");
is(round_4($ref->{fraction_disc}), 1, "swe_lun_occult_when_loc(2008 Venus)->fraction_disc");
is(round_4($ref->{core_shadow_km}), -3452.9813, "swe_lun_occult_when_loc(2008 Venus)->core_shadow_km");
is(round_4($ref->{body_azimuth}), 33.3583, "swe_lun_occult_when_loc(2008 Venus)->body_azimuth");
is(round_4($ref->{body_alt_true}), 11.6903, "swe_lun_occult_when_loc(2008 Venus)->body_alt_true");
is(round_4($ref->{separation_angle}), 0.0995, "swe_lun_occult_when_loc(2008 Venus)->separation_angle");
#------------------------------------------------------------------------
# planetary nodes and apsides
#------------------------------------------------------------------------
$ref = SwissEph::swe_nod_aps(2436723.588888888, 9, 4, 0);
is(round_4($ref->{xnasc}->[0]), 108.997, "swe_nod_aps(pluto)->node_asc");
is(round_4($ref->{xndsc}->[0]), 290.9051, "swe_nod_aps(pluto)->node_dsc");
is(round_4($ref->{xperi}->[0]), 223.3911, "swe_nod_aps(pluto)->perhelion");
is(round_4($ref->{xaphe}->[0]), 44.9470, "swe_nod_aps(pluto)->aphelion");
#------------------------------------------------------------------------
# heliacal risings
#------------------------------------------------------------------------
$ref = SwissEph::swe_heliacal_ut(2454800,[8,47,400],[1000,15,15,0.17],[23,1,0,0,0,0],"venus",1,4);
is(round_4($ref->{topt}), 2454914.7067, "swe_heliacal_ut(2008 Venus)->topt");
#for (sort keys(%$ref)) {print STDERR "$_ = $ref->{$_}\n";}
#print STDERR "d=".(-146780.0 - $ref->{time_equ})."\n";
if (1) {
print STDERR "\nend of official test section **************************\n";
my $year = 0.9; my $hour = 0.7;
my @aa = (60.0, 0, 1, 1, 0, 0);
my $xa = -1; my $xb = -1;
my @geopos = (8.33, 47.35,400);
#my ($xa, $xb) = SwissEph::test1(\@aa);
#my $ref = SwissEph::swe_houses(2436723.588888888, 47.35, 8.72, "P");
#my $ref = SwissEph::swe_time_equ(2436723.588888888);
#my $ref = SwissEph::swe_revjul(2436723.58888);
#my $ref = SwissEph::swe_houses_armc($ref->{armc}, 47.35, 23.45, "P");
#my $ref = SwissEph::swe_houses_armc_ex2($ref->{armc}, 47.35, 23.45, "P");
#my $ref = SwissEph::swe_houses_ex(2436723.588888888, 0, 47.35, 8.72, "P");
#my $ref = SwissEph::swe_houses_ex2(2436723.588888888, 0, 47.35, 8.72, "P");
#my $ref = SwissEph::swe_house_pos($ref->{armc}, 47.35, 23.45, "P", 72.73, 0);
#my $ref = SwissEph::swe_calc(2436723.588888888, 8, 256);
#my $ref2 = SwissEph::swe_get_current_file_data(0);
#for (sort keys(%$ref2)) {print STDERR "$_ = $ref2->{$_}\n";}
#my $ref = SwissEph::swe_calc_ut(2436723.588888888, 8, 256);
#my $ref = SwissEph::swe_fixstar("aldeb", 2436723.588888888, 0);
#my $ref = SwissEph::swe_fixstar_ut("aldeb", 2436723.588888888, 0);
#my $ref = SwissEph::swe_fixstar_mag("sirius");
#my $ref = SwissEph::swe_gauquelin_sector(2436723.588888888,0,"sirius",0,0,\@geopos,0,0);
#my $ref = SwissEph::swe_lun_occult_when_loc(2436723.588888888,0,"",0,0,\@geopos);
#my $ref = SwissEph::swe_sol_eclipse_when_loc(2436723.588888888,0,0,\@geopos);
#my $ref = SwissEph::swe_lun_occult_when_glob(2436723.588888888,0,"",0,0,0);
#my $ref = SwissEph::swe_sol_eclipse_when_glob(2452845.588888888,0,0,0);
#my $ref = SwissEph::swe_lun_eclipse_when(2452845.588888888,0,0,0);
#my $ref = SwissEph::swe_lun_occult_where($ref->{ecl_maximum},0, "",0);
#my $ref = SwissEph::swe_sol_eclipse_how($ref->{ecl_maximum},0,\@geopos);
#my $ref = SwissEph::swe_lun_eclipse_how($ref->{ecl_maximum},0,\@geopos);
#my $ref = SwissEph::swe_pheno(2436723.588888888,3,0);
#my $ref = SwissEph::swe_pheno_ut(2436723.588888888,3,0);
#my $ref = SwissEph::swe_refrac_extended(0.1,400,1013,15,0.02,1);
#my $dref = SwissEph::swe_refrac(1,1013,15,0);
#my $ref = SwissEph::swe_nod_aps(2436723.588888888, 9, 0, 0);
#for (sort keys(%$ref)) {print STDERR "$_ = $ref->{$_}\n";}
#print STDERR "@{$ref->{xaphe}}\n";
#print STDERR "@{$ref->{xperi}}\n";
#print STDERR "@{$ref->{xnasc}}\n";
#print STDERR "@{$ref->{xndsc}}\n";
#@xin = (252, 0, 1);
#my @arf = SwissEph::swe_azalt(2436723.588888888,0,\@geopos,1013,15,\@xin);
#my @arf = SwissEph::swe_azalt_rev(2436723.588888888,0,\@geopos,\@arf);
#print STDERR "dref = @arf\n";
#my $ref = SwissEph::swe_revjul($ref->{ecl_maximum});
#for (sort keys(%$ref)) {print STDERR "$_ = $ref->{$_}\n";}
#for (@{$ref->{cusps_speed}}) {print STDERR "xx $_\n";}
#for (@{$ref->{tret}}) {print STDERR "xx $_\n";}
#for (@{$ref->{attr}}) {print STDERR "xx $_\n";}
#print STDERR "$retval\n@$cusp\n@$ascmc\n";
#my $shnam = SwissEph::swe_house_name('c');
#print STDERR "$shnam\n";
#$ref = SwissEph::swe_version();
#for (sort keys(%$ref)) {print STDERR "$_ = $ref->{$_}\n";}
#$ref = SwissEph::swe_heliacal_ut(2454800,[8,47,400],[1000,15,15,0.17],[23,1,0,0,0,0],"venus",1,260);
#for (sort keys(%$ref)) {print STDERR "$_ = $ref->{$_}\n";}
#print STDERR "dret: ".join(", ", @{$ref->{dret}})."\n";
#$ref = SwissEph::swe_vis_limit_mag(2454915.7040719,[8,47,400],[1000,15,15,0.17],[23,1,0,0,0,0],"venus",256);
#for (sort keys(%$ref)) {print STDERR "$_ = $ref->{$_}\n";}
#$ref = SwissEph::swe_utc_time_zone(2000,1,1,1,59,60.1,5.5);
#for (sort keys(%$ref)) {print STDERR "$_ = $ref->{$_}\n";}
#for (@{$ref->{dret}}) {print STDERR "$_\n";}
#$ref = SwissEph::swe_utc_to_jd(2008,12,31,23,59,57.3);
#for (sort keys(%$ref)) {print STDERR "$_ = $ref->{$_}\n";}
#my $tjd0 = $ref->{tjd_et};
#for (my $i = 0; $i < 10; $i++) {
# $ref = SwissEph::swe_jdet_to_utc($tjd0 + $i/86400);
# print STDERR "$ref->{iday}.$ref->{imon}.$ref->{iyar}, $ref->{ihou}:$ref->{imin}:$ref->{dsec}\n";
#}
#$ref = SwissEph::swe_get_orbital_elements(2436723.5, 3, 0);
#for (sort keys(%$ref)) {print STDERR "$_ = $ref->{$_}\n";}
#$ref = SwissEph::swe_orbit_max_min_true_distance(2436723.5, 3, 0);
#for (sort keys(%$ref)) {print STDERR "$_ = $ref->{$_}\n";}
#my $ref = SwissEph::swe_calc_pctr(2436723.588888888, 8, 5, 0);
#my $pp = $ref->{xx};
#for (sort keys(%$ref)) {print STDERR "$_ = $ref->{$_}\n";}
#print STDERR "@$pp\n";
#------------------------------------------------------------------------
# Obliquity
#------------------------------------------------------------------------
}
sub round6 {
return sprintf( "%.0f", 1000000 * shift);
}
sub round_4 {
my $a = shift;
$a = int(10000 * ($a + 0.00005));
return $a / 10000.0;
}

54
src/LICENSE Normal file
View File

@ -0,0 +1,54 @@
/* 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.
*/

73
src/Makefile Normal file
View File

@ -0,0 +1,73 @@
# $Header$
# this Makefile creates a SwissEph library and a swetest sample on 64-bit
# Redhat Enterprise Linux RHEL 6.
# The mode marked as 'Linux' should also work with the GNU C compiler
# gcc on other systems.
# If you modify this makefile for another compiler, please
# let us know. We would like to add as many variations as possible.
# If you get warnings and error messages from your compiler, please
# let us know. We like to fix the source code so that it compiles
# free of warnings.
# send email to the Swiss Ephemeris mailing list.
all: libswe.so swemini swetest
CFLAGS = -g -Wall -fPIC '-DSE_EPHE_PATH="/usr/share/libswe"' # for Linux and other gcc systems
OP=$(CFLAGS)
CC=cc #for Linux
# compilation rule for general cases
.o :
$(CC) $(OP) -o $@ $? -lm
.c.o:
$(CC) -c $(OP) $<
SWEOBJ = swedate.o swehouse.o swejpl.o swemmoon.o swemplan.o swepcalc.o sweph.o\
swepdate.o swephlib.o swecl.o swehel.o
swetest: swetest.o libswe.a
$(CC) $(OP) -o swetest swetest.o -L. -lswe -lm -ldl
swemini: swemini.o libswe.a
$(CC) $(OP) -o swemini swemini.o -L. -lswe -lm -ldl
# create an archive and a dynamic link libary fro SwissEph
# a user of this library will inlcude swephexp.h and link with -lswe
libswe.a: $(SWEOBJ)
ar r libswe.a $(SWEOBJ)
libswe.so: $(SWEOBJ)
$(CC) -shared -o libswe.so $(SWEOBJ)
INST_H_FILES = \
swephexp.h \
sweodef.h \
$(NULL)
install: libswe.so swemini
mkdir -p $(DESTDIR)/usr/bin
mkdir -p $(DESTDIR)/usr/lib64 $(DESTDIR)/usr/include
install swemini $(DESTDIR)/usr/bin/
install libswe.so $(DESTDIR)/usr/lib64/
install -m 0644 swephexp.h sweodef.h $(DESTDIR)/usr/include/
clean:
rm -f *.o swetest libswe*
###
swecl.o: swejpl.h sweodef.h swephexp.h swedll.h sweph.h swephlib.h
sweclips.o: sweodef.h swephexp.h swedll.h
swedate.o: swephexp.h sweodef.h swedll.h
swehel.o: swephexp.h sweodef.h swedll.h
swehouse.o: swephexp.h sweodef.h swedll.h swephlib.h swehouse.h
swejpl.o: swephexp.h sweodef.h swedll.h sweph.h swejpl.h
swemini.o: swephexp.h sweodef.h swedll.h
swemmoon.o: swephexp.h sweodef.h swedll.h sweph.h swephlib.h
swemplan.o: swephexp.h sweodef.h swedll.h sweph.h swephlib.h swemptab.h
swepcalc.o: swepcalc.h swephexp.h sweodef.h swedll.h
sweph.o: swejpl.h sweodef.h swephexp.h swedll.h sweph.h swephlib.h
swephlib.o: swephexp.h sweodef.h swedll.h sweph.h swephlib.h
swetest.o: swephexp.h sweodef.h swedll.h

680
src/ep4/sweephe4.c Normal file
View File

@ -0,0 +1,680 @@
/********************************************************************
ephe.c
access structures and functions for ephemeris file ep4_
*********************************************************************/
/* Copyright (C) 1997 - 2020 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 "swephexp.h"
# include "sweephe4.h"
# include <string.h>
# define INVALID_BASE 2000000000L
# define EPBS (2 * NDB) /* buffer size is 20 days */
# define EP_MIN_IX 2 /* load buffer when index below this */
# define EP_MAX_IX (EPBS - 4) /* load buffer when index above this */
FILE *ephfp = NULL;
const int qod[EP_NP] = {5,5,5,5,5,3,3,3,3,3,3,5,3,3,3};
static void inpolq_l(int n, int o, double p, centisec *x,
centisec *axu, centisec *adxu);
static int inpolq(int n, int o, double p, double *x,
double *axu, double *adxu);
static int ephe4_unpack(int jdl, int pflag, centisec lon[][EPBS], int i0,
char *errs);
static int ephe4_unpack_d(int jdl, int pflag, double lon[][EPBS], int i0,
char *errs);
static char *my_makepath(char *d, char *s);
# ifdef INTEL_BYTE_ORDER
/********************************************************************/
void shortreorder (UCHAR *p, int n)
/* p points to memory filled with 16-bit values; for
each of the values the seqeuence of the two bytes
has to be reversed, to translate HP-UX and VAX
ordering to MSDOS/Turboc ordering */
{
int i;
unsigned char c0;
for (i = 0; i < n; i += 2, p += 2) {
c0 = *p;
*p = *(p+1);
*(p + 1) = c0;
}
}
# endif
/****************************************************
read ephe file and return a pointer to normalized positions
for planets specified by pflag at julday jdl;
If the reading failes, NULL is returned and
an error text of max. 79 char in errtext, except when errtext = NULL.
If calc is used and succeeds, a message is put into errtext, otherwise
errtext is empty.
Attention: jd is an absolute Julian date
flag bits defined in sweephe4.h, onl
****************************************************/
centisec *ephread(double jd, int plalist, int flag, char *errtext)
{
static int jdbase = INVALID_BASE;
static int lastplalist = 0;
static centisec lon[EP_NP][EPBS]; /* buffer for 20 days unpacked ephe */
static centisec out[2 * EP_NP]; /* buffer for return longitude
and return speed */
int p, pf;
int ix, jdlong, iflagret;
centisec clp;
double jfract;
double x[6];
if (errtext != NULL)
*errtext = '\0';
if (plalist == 0)
plalist = EP_ALL_BITS; /* default: all logitudes, no speeds */
/*
* we must determine when to reload the lon buffer: if the contents do
* not allow immediate interpolation or if the plalist selection has
* changed since the last call.
*/
if ((plalist & lastplalist) != plalist) { /* new set is not contained in old */
jdbase = INVALID_BASE;
}
lastplalist = plalist;
jdlong = floor(jd - 0.5);
ix = jdlong - jdbase;
if (ix < EP_MIN_IX || ix >= EPBS) { /* must reload full buffer */
jdbase = ((jdlong - EP_MIN_IX) / NDB) * NDB; /* new base */
if (jdbase > jdlong - EP_MIN_IX) jdbase -= NDB; /* fix bug for neg. */
if (ephe4_unpack (jdbase, plalist, lon, 0, errtext) != OK)
goto err_exit;
if (ephe4_unpack (jdbase + NDB, plalist, lon, 0 + NDB, errtext) != OK)
goto err_exit;
ix = jdlong - jdbase;
} else if (ix > EP_MAX_IX) { /* must shift upper half down
and reload upper half of buffer */
jdbase += NDB; /* new base */
for (p = 0; p < EP_NP; p++)
memcpy (&lon[p][0], &lon[p][NDB], NDB * sizeof(centisec));
if (ephe4_unpack (jdbase + NDB, plalist, lon, 0 + NDB, errtext) != OK)
goto err_exit;
ix = jdlong - jdbase;
}
jfract = jd - 0.5 - jdlong;
/*
* we use the interpolator even for jfract = 0, because it delivers
* the speed term. The computation overhead is unimportant
* in any case.
*/
for (p = 0, pf = 1; p < EP_NP; p++, pf = pf << 1)
if ((plalist & pf) != 0) {
inpolq_l((int) ix, qod[p], jfract, &(lon[p][0]), &(out[p]), &clp);
if (p <= CHIRON) { /* normalize all except ecl and nut */
if (out[p] < 0)
out[p] += DEG360;
else if (out[p] >= DEG360)
out[p] -= DEG360;
}
#ifdef DEBUG
fprintf(stderr,"ephread p=%d, lon=%.3lf\n", p, out[p] * CS2DEG);
#endif
if (flag & EP_BIT_SPEED)
out[p+EP_NP] = clp;
}
return out;
err_exit:
jdbase = INVALID_BASE;
lastplalist = 0;
if ((flag & EP_BIT_MUST_USE_EPHE) == 0) { /* try using calc */
int sweflag = 0;
char serr[AS_MAXCH];
if (flag & EP_BIT_SPEED)
sweflag = SEFLG_SPEED;
if (errtext != NULL)
sprintf(errtext,"ephread failed for jd=%f; used swe_calc().", jd);
for (p = 0, pf = 1; p < CALC_N; p++, pf = pf << 1) {
if ((plalist & pf) != 0) {
if ((iflagret = swe_calc(jd, plac2swe(p), sweflag, x, serr)) != ERR) {
out[p] = d2l(x[0] * DEG);
if (flag & EP_BIT_SPEED)
out[p + EP_NP] = d2l(x[3] * DEG);
if (out[p] < 0)
out[p] += DEG360;
else if (out[p] >= DEG360)
out[p] -= DEG360;
} else {
swe_close();
if (errtext != NULL)
strcat(errtext," calc failed too.");
return NULL;
}
}
}
if ((iflagret = swe_calc(jd, SE_ECL_NUT, 0, x, serr)) == ERR) {
swe_close();
sprintf(errtext, "error in swe_calc() %s\n", serr);
return NULL;
}
out[EP_ECL_INDEX] = d2l(x[0] * DEG); /* true ecliptic */
out[EP_NUT_INDEX] = d2l(x[2] * DEG); /* nutation */
out[EP_ECL_INDEX + EP_NP] = 0;
out[EP_NUT_INDEX + EP_NP] = 0;
return out;
}
return NULL;
} /* ephread */
// same in double
double *dephread2(double jd, int plalist, int flag, char *errtext)
{
static int jdbase = INVALID_BASE;
static int lastplalist = 0;
static double lon[EP_NP][EPBS]; // buffer for 20 days unpacked ephe
static double out[2 * EP_NP]; // buffer for return longitude and return speed
int p, pf;
int ix, jdlong, iflagret;
double lp;
double jfract;
double x[6];
if (errtext != NULL)
*errtext = '\0';
if (plalist == 0)
plalist = EP_ALL_BITS; /* default: all logitudes, no speeds */
/*
* we must determine when to reload the lon buffer: if the contents do
* not allow immediate interpolation or if the plalist selection has
* changed since the last call.
*/
if ((plalist & lastplalist) != plalist) { /* new set is not contained in old */
jdbase = INVALID_BASE;
}
lastplalist = plalist;
jdlong = floor(jd - 0.5);
ix = jdlong - jdbase;
if (ix < EP_MIN_IX || ix >= EPBS) { /* must reload full buffer */
jdbase = ((jdlong - EP_MIN_IX) / NDB) * NDB; /* new base */
if (jdbase > jdlong - EP_MIN_IX) jdbase -= NDB; /* fix bug for neg. */
if (ephe4_unpack_d(jdbase, plalist, lon, 0, errtext) != OK)
goto err_exit;
if (ephe4_unpack_d(jdbase + NDB, plalist, lon, 0 + NDB, errtext) != OK)
goto err_exit;
ix = jdlong - jdbase;
} else if (ix > EP_MAX_IX) { /* must shift upper half down
and reload upper half of buffer */
jdbase += NDB; /* new base */
for (p = 0; p < EP_NP; p++)
memcpy(&lon[p][0], &lon[p][NDB], NDB * sizeof(double));
if (ephe4_unpack_d(jdbase + NDB, plalist, lon, 0 + NDB, errtext) != OK)
goto err_exit;
ix = jdlong - jdbase;
}
jfract = jd - 0.5 - jdlong;
/*
* we use the interpolator even for jfract = 0, because it delivers
* the speed term. The computation overhead is unimportant
* in any case.
*/
for (p = 0, pf = 1; p < EP_NP; p++, pf = pf << 1)
if ((plalist & pf) != 0) {
inpolq((int) ix, qod[p], jfract, &(lon[p][0]), &(out[p]), &lp);
if (p <= CHIRON) { /* normalize all except ecl and nut */
if (out[p] < 0)
out[p] += 360.0;
else if (out[p] >= 360.0)
out[p] -= 360.0;
}
#ifdef DEBUG
fprintf(stderr,"ephread p=%d, lon=%.3lf\n", p, out[p]);
#endif
if (flag & EP_BIT_SPEED)
out[p+EP_NP] = lp;
}
return out;
err_exit:
jdbase = INVALID_BASE;
lastplalist = 0;
if ((flag & EP_BIT_MUST_USE_EPHE) == 0) { /* try using calc */
int sweflag = 0;
char serr[AS_MAXCH];
if (flag & EP_BIT_SPEED)
sweflag = SEFLG_SPEED;
if (errtext != NULL)
sprintf(errtext,"ephread failed for jd=%f; used swe_calc().", jd);
for (p = 0, pf = 1; p < CALC_N; p++, pf = pf << 1) {
if ((plalist & pf) != 0) {
if ((iflagret = swe_calc(jd, plac2swe(p), sweflag, x, serr)) != ERR) {
out[p] = x[0];
if (flag & EP_BIT_SPEED)
out[p + EP_NP] = x[3];
if (out[p] < 0)
out[p] += 360.0;
else if (out[p] >= 360.0)
out[p] -= 360.0;
} else {
swe_close();
if (errtext != NULL)
strcat(errtext," calc failed too.");
return NULL;
}
}
}
if ((iflagret = swe_calc(jd, SE_ECL_NUT, 0, x, serr)) == ERR) {
swe_close();
sprintf(errtext, "error in swe_calc() %s\n", serr);
return NULL;
}
out[EP_ECL_INDEX] = x[0]; /* true ecliptic */
out[EP_NUT_INDEX] = x[2]; /* nutation */
out[EP_ECL_INDEX + EP_NP] = 0;
out[EP_NUT_INDEX + EP_NP] = 0;
return out;
}
return NULL;
}
/****************************************************
unpack an ephe file block specified by jlong
and the planets specified by pflag into
the array lon[p][EPBS], starting at index i0.
jdl is (long) floor(full julian date);
If the reading failes, ERR is returned and
an error text of max. 79 char in errs, except when errs = NULL.
****************************************************/
static int ephe4_unpack(int jdl, int plalist, centisec lon[][EPBS], int i0,char *errs)
{
int p, i, pf;
centisec l_ret, d_ret;
struct ep4 e;
if (eph4_posit (jdl, FALSE, errs) != OK)
return (ERR);
if (fread (&e, sizeof(struct ep4), 1, ephfp) != 1) {
if (errs != NULL)
sprintf (errs, "ephe4_unpack: fread for jd=%d failed", jdl);
return (ERR);
}
#ifdef INTEL_BYTE_ORDER
shortreorder((UCHAR *) &e, sizeof(struct ep4));
#endif
for (p = SUN, pf = 1; p <= CHIRON; p++, pf = pf << 1) {
if ((plalist & pf) == 0) continue;
l_ret = e.elo[p].p0m * 6000L + e.elo[p].p0s; /* csec */
d_ret = e.elo[p].pd1m * 6000L + e.elo[p].pd1s; /* csec */
lon[p][i0] = l_ret;
l_ret += d_ret;
if (l_ret < 0) {
lon[p][i0+1] = l_ret + DEG360;
} else if (l_ret >= DEG360) {
lon[p][i0+1] = l_ret - DEG360;
} else {
lon[p][i0+1] = l_ret;
}
for (i = 2; i < NDB; i++) {
if (p == MOON || p == MERCURY)
d_ret += e.elo[p].pd2[i-2] * 10L;
else
d_ret += e.elo[p].pd2[i-2];
l_ret += d_ret;
if (l_ret < 0) {
lon[p][i0+i] = l_ret + DEG360;
} else if (l_ret >= DEG360) {
lon[p][i0+i] = l_ret - DEG360;
} else {
lon[p][i0+i] = l_ret;
}
}
} /* for p */
if ( plalist & EP_ECL_BIT) { /* unpack ecl */
l_ret = e.ecl0m * 6000L + e.ecl0s;
lon[EP_ECL_INDEX][i0] = l_ret;
for (i = 1; i < NDB; i++)
lon[EP_ECL_INDEX][i0+i] = l_ret + e.ecld1[i-1];
}
if ( plalist & EP_NUT_BIT) { /* unpack nut */
for (i = 0; i < NDB; i++)
lon[EP_NUT_INDEX][i0+i] = e.nuts[i] ;
}
return OK;
} /* ephe4_unpack */
// same in double
static int ephe4_unpack_d(int jdl, int plalist, double lon[][EPBS], int i0,char *errs)
{
int p, i, pf;
double l_ret, d_ret;
struct ep4 e;
if (eph4_posit (jdl, FALSE, errs) != OK)
return (ERR);
if (fread (&e, sizeof(struct ep4), 1, ephfp) != 1) {
if (errs != NULL)
sprintf (errs, "ephe4_unpack: fread for jd=%d failed", jdl);
return (ERR);
}
#ifdef INTEL_BYTE_ORDER
shortreorder((UCHAR *) &e, sizeof(struct ep4));
#endif
for (p = SUN, pf = 1; p <= CHIRON; p++, pf = pf << 1) {
if ((plalist & pf) == 0) continue;
l_ret = (e.elo[p].p0m * 6000 + e.elo[p].p0s) * CS2DEG;
d_ret = (e.elo[p].pd1m * 6000 + e.elo[p].pd1s) * CS2DEG;
lon[p][i0] = l_ret;
l_ret += d_ret;
if (l_ret < 0) {
lon[p][i0+1] = l_ret + 360.0;
} else if (l_ret >= 360.0) {
lon[p][i0+1] = l_ret - 360.0;
} else {
lon[p][i0+1] = l_ret;
}
for (i = 2; i < NDB; i++) {
if (p == MOON || p == MERCURY)
d_ret += (e.elo[p].pd2[i-2] * 10 * CS2DEG);
else
d_ret += (e.elo[p].pd2[i-2] * CS2DEG);
l_ret += d_ret;
if (l_ret < 0) {
lon[p][i0+i] = l_ret + 360.0;
} else if (l_ret >= 360.0) {
lon[p][i0+i] = l_ret - 360.0;
} else {
lon[p][i0+i] = l_ret;
}
}
} /* for p */
if ( plalist & EP_ECL_BIT) { /* unpack ecl */
l_ret = (e.ecl0m * 6000L + e.ecl0s) * CS2DEG;
lon[EP_ECL_INDEX][i0] = l_ret;
for (i = 1; i < NDB; i++)
lon[EP_ECL_INDEX][i0+i] = l_ret + e.ecld1[i-1] * CS2DEG;
}
if ( plalist & EP_NUT_BIT) { /* unpack nut */
for (i = 0; i < NDB; i++)
lon[EP_NUT_INDEX][i0+i] = e.nuts[i] * CS2DEG ;
}
return OK;
}
/****************************************************
position ephe file at proper position for julian
date jd; if writeflag = TRUE (write mode), create file
if required. Return OK or ERR.
globals used: ephfp.
*****************************************************/
int eph4_posit (int jlong, AS_BOOL writeflag, char *errtext)
{
int filenr;
long posit;
static int open_filenr = -10000;
char fname[AS_MAXCH], s[80], *sp;
filenr = jlong / EP4_NDAYS;
if (jlong < 0 && filenr * EP4_NDAYS != jlong) filenr--;
posit = jlong - filenr * EP4_NDAYS;
posit = posit / NDB * sizeof(struct ep4);
if (open_filenr != filenr) {
if (ephfp != NULL) {
fclose(ephfp);
open_filenr = -10000;
}
if (filenr >= 0)
sprintf (s, "%s%s%d", EP4_PATH, EP4_FILE, filenr);
else
sprintf (s, "%s%sM%d", EP4_PATH, EP4_FILE, -filenr);
my_makepath(fname, s);
if (writeflag)
sp = BFILE_W_CREATE;
else
sp = BFILE_R_ACCESS;
ephfp = fopen (fname, sp);
if (ephfp == NULL) {
if (errtext != NULL) {
if (! writeflag) {
sprintf (errtext,"eph4_posit: file %s does not exist\n", fname);
} else {
sprintf (errtext,"eph4_posit: could not create file %s\n", fname);
}
}
return (ERR);
}
open_filenr = filenr;
}
if (fseek (ephfp, posit, 0) == 0 && ftell(ephfp) == posit) {
return (OK);
} else {
if (errtext != NULL)
sprintf (errtext,"eph4_posit: fseek(%ld) of file nr %d failed\n",
posit, open_filenr);
return (ERR);
}
} /* end eph4_posit */
/*****************************************************
quicker Everett interpolation, after Pottenger
version for long, 17.7.91 by Alois Treindl
*****************************************************/
static void inpolq_l(int n, int o, double p, centisec *x, centisec *axu, centisec *adxu)
/*
* interpolate between x[n] and x[n-1], at argument n+p
* o = order of interpolation, maximum 5
* p = argument in [0..1]
* x[] array of function values, x[n-2]..x[n+3] must exist
* axu pointer for storage of result
* adxu pointer for storage of dx/dt
*/
{
static double q,q2,q3,q4,q5,
p2,p3,p4,p5,
u,u0,u1,u2;
static double lastp = 9999;
double rl, rlp;
centisec dm2,dm1,d0,dp1,dp2,
d2m1,d20,d2p1,d2p2,
d30,d3p1,d3p2,
d4p1,d4p2;
centisec offset = 0;
if (lastp != p) { /* recompute the interpolator factors */
q=1.0-p;
q2 = q*q;
q3 = (q+1.0)*q*(q-1.0)/6.0; /* q - 1 over 3; u5 */
p2 = p*p;
p3 = (p+1.0)*p*(p-1.0)/6.0; /* p - 1 over 3; u8 */
u = (3.0*p2-1.0)/6.0;
u0 = (3.0*q2-1.0)/6.0;
q4 = q2*q2; /* f5 */
p4 = p2*p2; /* f4 */
u1 = (5.0*p4-15.0*p2+4.0)/120.0; /* u1 */
u2 = (5.0*q4-15.0*q2+4.0)/120.0; /* u2 */
q5 = q3*(q+2.0)*(q-2.0)/20.0; /* q - 2 over 5; u6 */
p5 = (p+2.0)*p3*(p-2.0)/20.0; /* p - 2 over 5; u9 */
lastp = p;
}
dm1 = x[n] - x[n-1];
if (dm1 >= DEG180)
dm1 -= DEG360;
else if (dm1 < -DEG180)
dm1 += DEG360;
d0 = x[n+1] - x[n];
if (d0 >= DEG180) {
d0 -= DEG360;
offset = DEG360;
} else if (d0 < -DEG180) {
d0 += DEG360;
offset = -DEG360;
}
dp1 = x[n+2] - x[n+1];
if (dp1 >= DEG180)
dp1 -= DEG360;
else if (dp1 < -DEG180)
dp1 += DEG360;
d20 = d0 - dm1; /* f8 */
d2p1 = dp1 - d0; /* f9 */
/*
* Everett interpolation 3rd order
*/
rl = q*(x[n] + offset) + q3*d20 + p*x[n+1] + p3*d2p1;
rlp = d0 + u*d2p1 - u0*d20;
if ( o > 3 ) { /* 5th order */
dm2 = x[n-1] - x[n-2];
if (dm2 >= DEG180)
dm2 -= DEG360;
else if (dm2 < -DEG180)
dm2 += DEG360;
dp2 = x[n+3] - x[n+2];
if (dp2 >= DEG180)
dp2 -= DEG360;
else if (dp2 < -DEG180)
dp2 += DEG360;
d2m1 = dm1 - dm2;
d2p2 = dp2 - dp1;
d30 = d20 - d2m1;
d3p1 = d2p1 - d20;
d3p2 = d2p2 - d2p1;
d4p1 = d3p1 - d30; /* f7 */
d4p2 = d3p2 - d3p1; /* f */
rl += p5*d4p2 + q5*d4p1;
rlp += u1*d4p2 - u2*d4p1;
}
*axu = d2l (rl);
*adxu = d2l (rlp);
} /* end inpolq_l() */
/*****************************************************
quicker Everett interpolation, after Pottenger
version for double 9 Jul 1988 by Alois Treindl
return OK, no error checking
Was used in Placalc to interpolate 80-day stored ephe for outer planets.
*****************************************************/
static int inpolq(int n, int o, double p, double *x, double *axu, double *adxu)
// n interpolate between x[n] and x[n-1], at argument n+p
// o order of interpolation, maximum 5
// p, argument , intervall [0..1]
// x[] array of function values, x[n-o]..x[n+o] must exist
// *axu pointer for storage of result
// *adxu pointer for storage of dx/dt
{
static double q,q2,q3,q4,q5,p2,p3,p4,p5, u,u0,u1,u2;
static double lastp = 9999.0;
double dm2,dm1,d0,dp1,dp2,
d2m1,d20,d2p1,d2p2,
d30,d3p1,d3p2,
d4p1,d4p2;
double offset = 0.0;
if (lastp != p) {
q=1.0-p;
q2 = q*q;
q3 = (q+1.0)*q*(q-1.0)/6.0; /* q - 1 over 3; u5 */
p2 = p*p;
p3 = (p+1.0)*p*(p-1.0)/6.0; /* p - 1 over 3; u8 */
u = (3.0*p2-1.0)/6.0;
u0 = (3.0*q2-1.0)/6.0;
q4 = q2*q2; /* f5 */
p4 = p2*p2; /* f4 */
u1 = (5.0*p4-15.0*p2+4.0)/120.0; /* u1 */
u2 = (5.0*q4-15.0*q2+4.0)/120.0; /* u2 */
q5 = q3*(q+2.0)*(q-2.0)/20.0; /* q - 2 over 5; u6 */
p5 = (p+2.0)*p3*(p-2.0)/20.0; /* p - 2 over 5; u9 */
lastp = p;
}
dm1 = x[n] - x[n-1];
if (dm1 > 180.0) dm1 -= 360.0;
if (dm1 < -180.0) dm1 += 360.0;
d0 = x[n+1] - x[n];
if (d0 > 180.0) {
d0 -= 360.0;
offset = 360.0;
}
if (d0 < -180.0) {
d0 += 360.0;
offset = -360.0;
}
dp1 = x[n+2] - x[n+1];
if (dp1 > 180.0) dp1 -= 360.0;
if (dp1 < -180.0) dp1 += 360.0;
d20 = d0 - dm1; /* f8 */
d2p1 = dp1 - d0; /* f9 */
/* Everett interpolation 3rd order */
*axu = q*(x[n] + offset) + q3*d20
+ p*x[n+1] + p3*d2p1;
*adxu = d0 + u*d2p1 - u0*d20;
if ( o > 3 ) { /* 5th order */
dm2 = x[n-1] - x[n-2];
if (dm2 > 180.0) dm2 -= 360.0;
if (dm2 < -180.0) dm2 += 360.0;
dp2 = x[n+3] - x[n+2];
if (dp2 > 180.0) dp2 -= 360.0;
if (dp2 < -180.0) dp2 += 360.0;
d2m1 = dm1 - dm2;
d2p2 = dp2 - dp1;
d30 = d20 - d2m1;
d3p1 = d2p1 - d20;
d3p2 = d2p2 - d2p1;
d4p1 = d3p1 - d30; /* f7 */
d4p2 = d3p2 - d3p1; /* f */
*axu += p5*d4p2 + q5*d4p1;
*adxu += u1*d4p2 - u2*d4p1;
}
return (OK);
} /* end inpolq() */
static char *my_makepath(char *d, char *s)
{
char *getenv();
if (*s == *DIR_GLUE || *s == '/' || strchr (s, ':') != NULL) {
strcpy (d, s); /* s is absolute path name */
}
# if MSDOS
while ((p = strchr(d, '/')) != NULL) *p = '\\';
# endif
return (d);
}

210
src/ep4/sweephe4.h Normal file
View File

@ -0,0 +1,210 @@
/*******************************************************
header file for structures and functions in module ephe.c
for reading and writing stored ephemerides in format ep4
The design of ephemeris type ep4:
In all ASYS and transit application of stored ephemerides
except Progressed Declinations Type 56 we need only the
longitudes of the planets or nodes.
The old EP3 format contains also latitudes, and uses ephemeris time.
Therefore we define a new ephemeris format, which is capable of
replacing EP3, when some ASYS programs are changed.
The ASYS programs requiring different data can receive them
by asking the calcserv module.
We therefore store now a daily ephemeris with only logitudes, ecl and nut.
The ephemeris is computed and stored for midnight ephemeris time, i.e.
for jd = 2400000.5, 2400001.5 etc.
In the ephemeris record for this date, only floor(jd) is kept.
In many cases universal time (UT) is desired, not ephemeris time.
Because computation with our current computers is very cheap for
everything except trigonometrci functions, we can afford to
build a good interpolation into the ephemeris package.
The user can request from ephread() ephemeris positions for
any (double) jd, not only for midnight ephemeris time.
Inside the package the same quick Everett 5th-order interpolator
is used as in placalc.
It delivers positions within 0.01" for all planets except moon, mercury
and true node. Moon and Mercury suffer, because their positions are
stored only with a precision of 0.1"; true node suffers because
it oscillates quickly with the fastest moon terms.
The maximum observed differences between placalc and ephread for 12.00 noon
are 0.25" for moon and true node and 0.1" for Mercury; in 80% of the days
the differences are less than 0.1". This is significantly better than
the implemented precision of the placalc package itself.
The Everett interpolator delivers without any additional cost also
the speed of the planets. This speed is very much better than the
speed derived for the inner planets from the mean orbits.
The returned data from ephread are in an array of centisec,
with ecl and nut behind the planets.
The default, pflag = 0, returns all.
The speeds are returned in the second half of the array;
the speed is always there, even when the speed bit has not been set.
***********************************************************/
/* Copyright (C) 1997 - 2020 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.
*/
# ifndef _EPHE_INCLUDED
# define _EPHE_INCLUDED
#ifdef __cplusplus
extern "C" {
#endif
#include "astrolib.h"
# define EP4_BLOCKSIZE sizeof(struct ep4)
# if HPUNIX
# if 0
# define EP4_PATH "/users/dieter/"
# else
# define EP4_PATH "/users/ephe/"
# endif
# else
# define EP4_PATH "ephe\\"
# endif
# define EP4_FILE "sep4_" /* packed ephemeris */
# define EP4_NDAYS 10000L /* days per EP4_ file */
# define NDB 10L /* 10 days per block */
/*
* bits for plalist in ephread():
* the planet flag bits SUN .. CHIRON, ECl, NUT can be set individually.
* plalist = 0 is equivalent to all planets and ecl,nut
* EP_ALL_PLANETS sets all planet bits SUN .. CHIRON
* EP_ALL_BITS sets all bits.
*/
# define EP_NP (CHIRON + 3) /* total number of factors in ep4 */
/* sun .. chiron, ecl, nut */
# define EP_ALL_PLANETS ((1 << (CHIRON + 1)) - 1) /* bits 0..12 set */
# define EP_CALC_N (CHIRON+1) /* 13 planets, SUN .. CHIRON */
# define EP_ECL_INDEX (CHIRON + 1) /* index for ecliptic centisec */
# define EP_NUT_INDEX (CHIRON + 2) /* index for nutation centisec */
# define EP_ECL_BIT (1 << EP_ECL_INDEX)
# define EP_NUT_BIT (1 << EP_NUT_INDEX)
# define EP_ALL_BITS (EP_ALL_PLANETS|EP_ECL_BIT|EP_NUT_BIT)
// bits for flag in ephread(), values come from placalc compatibility
# define EP_BIT_SPEED 16 // must get speed
# define EP_BIT_MUST_USE_EPHE 256
struct elon { /* longitudes for 10 days */
short p0m; /* longitude at day 0, minutes */
short p0s; /* 0.01" */
short pd1m; /* delta of days 1, 0.01" */
short pd1s; /* 0.01" */
short pd2[NDB-2]; /* second differences, day 2 .. 9,
0.1" moon, mercury, 0.01" others*/
};
/*
* ep4 is the new ephemeris format for longitude only, ephemeris time
*/
struct ep4 {
short j_10000; /* floor(jd - 0.5) / 10000L; */
short j_rest; /* (jd - 0.5 ) - 10000L * j_10000
j_rest is always positive;
jd = j_10000 * 10000L + j_rest + 0.5 */
short ecl0m; /* true ecliptic day 0, min ; */
short ecl0s; /* 0.01" */
short ecld1[NDB-1]; /* first differences 0.01", day 1..9 */
short nuts[NDB]; /* nutation in 0.01", day 0..9 */
struct elon elo[CHIRON +1]; /* longitude sun...chiron */
};
/******************************************
globals exported by module ephe.c
********************************************/
extern FILE *ephfp;
/******************************************
functions exported by module ephe.c
********************************************/
extern centisec *ephread(double jd, int plalist, int flag, char *errtext);
/*
* This is the only function normally used by applications.
* ATTENTION: jd is an absolute Julian date, whereas calc() and deltat()
* require Astrodienst-relative Julian dates.
* plalist can consist of individual planet bit flags, to indicate that
* only these planets are wanted.
* plalist = 0 returns all planets, and ecl and nut.
* Because the computation is so fast, it is recommended to use pflag = 0.
* flag recognizes only the bits CALC_BIT_SPEED and CALC_MUST_USE_EPHE.
* If CALC_BIT_SPEED is set, of the planets are returned in the result array
* after all longitudes (speeds for ecl and nut are always set to zero).
* If CALC_BIT_MUST_USE_EPHE is NOT set, calc() will be used if reading
* the ephemeris fails.
* cp is the returned pointer, the longitude of the planet p is
* at cp[p] and the speed is at cp[p + EP_NP].
* The returned longitudes are always normalized into [0..DEG360[,
* except for nut, which is small and close to zero, negative or positive.
*/
extern double *dephread2(double jd, int plalist, int flag, char *errtext);
extern int eph4_posit (int jlong, AS_BOOL writeflag, char *errtext);
#ifdef __cplusplus
}
#endif
# endif /* _EPHE_INCLUDED */

310
src/ep4/swephgen4.c Normal file
View File

@ -0,0 +1,310 @@
/********************************************************************
sweephgen4.c
Create ephemeris file type 4 ep4_
options: -fYYY (start) file number, required option
-nNN number of files to be created, default 1
-v verbose: print differences (default: no)
-t test by reading
File format:
1000 blocks of xxx bytes
File names: ep4_243, ep4_244
corresponding to the absolute julian day number
*********************************************************************/
/* Copyright (C) 1997 - 2020 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 "swephexp.h"
# include "sweephe4.h"
# include "ourfiles.h"
# include "astrolib.h"
# define EPHR_NPL (CHIRON + 1)
char *arg0;
int32 max_dd[EP_CALC_N]; /* remember maximum of second dfifferences */
double max_err[EP_CALC_N]; /* remember maximum error */
AS_BOOL verbose = FALSE;
char errtext[AS_MAXCH];
int split(w, m, min, sec)
int32 w; /* position in seconds/m */
int m; /* factor for seconds */
short *min, /* storage for degrees and minutes */
*sec; /* storage for seconds * m */
{
if (w >= 0) {
*sec = w % (60 * m);
*min = w / (60 * m);
} else {
*sec = -(-w % (60 * m));
*min = -(-w / (60 * m));
}
return OK;
}
/*************************************************************
Pack positions of 10 days and write to file
The longitude is packed with second differences in such a way,
that the accumulating rounding erros do not exceed half of
the last stored digit, i.e. 0.05" moon, 0.005" other planets
**************************************************************/
int eph4_pack (int32 jd, double (*l)[NDB], double ecliptic[],
double nutation[])
{
int i, p,ps;
int32 d1, d2, dd, d_ret, w0, w_ret;
double err;
struct ep4 e;
e.j_10000 = jd / 10000.0;
e.j_rest = jd - 10000.0 * e.j_10000;
w0 = d2l( ecliptic[0] * DEG);
split( w0, 100, &e.ecl0m, &e.ecl0s );
for (i = 1; i < NDB; i++)
e.ecld1[i-1] = d2l(ecliptic[i] * DEG - w0);
for (i = 0; i < NDB; i++)
e.nuts[i] = d2l( nutation[i] * DEG ); /* int32 casted into short */
for (p = SUN; p <= CHIRON ; p++) {
ps = p;
w0 = d2l( l[ps][0] * DEG);
d1 = d2l( l[ps][1] * DEG - w0);
if (d1 >= DEG180)
d1 -= DEG360;
else if (d1 <= -DEG180)
d1 += DEG360;
split(w0, 100, &e.elo[p].p0m, &e.elo[p].p0s);
split(d1, 100, &e.elo[p].pd1m, &e.elo[p].pd1s);
d_ret = d1; /* recalculated diff */
w_ret = w0 + d_ret; /* recalculated position */
for (i = 2; i < NDB; i++) {
d2 = d2l( l[ps][i] * DEG - w_ret);
if (d2 >= DEG180)
d2 -= DEG360;
else if (d2 <= -DEG180)
d2 += DEG360;
dd = d2 - d_ret; /* second difference */
if (p == MOON || p == MERCURY)
dd = d2l(dd / 10.0); /* moon only 0.1" */
if (verbose && abs(dd) > abs(max_dd[ps]))
max_dd[ps] = dd;
e.elo[p].pd2[i-2] = dd;
if (p == MOON || p == MERCURY)
d_ret += e.elo[p].pd2[i-2] * 10L;
else
d_ret += e.elo[p].pd2[i-2];
w_ret += d_ret;
if (verbose) {
err = swe_difdeg2n(w_ret/360000.0, l[ps][i]); /* error */
if (fabs(err) > fabs(max_err[ps]))
max_err[ps] = err;
}
}
} /* for p */
#ifdef INTEL_BYTE_ORDER
shortreorder((UCHAR *) &e, sizeof(struct ep4));
#endif
fwrite (&e, sizeof(struct ep4), 1, ephfp);
return (OK);
}
/*************************************/
char *degstr (t)
double t;
{
static char a[20]; /* must survive call */
double min, sec;
int ideg, imin;
char sign = ' ';
if ( t < 0) sign = '-';
t = fabs (t);
ideg = (int) floor (t);
min = ( t - ideg ) * 60.0;
imin = (int) floor(min);
sec = ( min - imin ) * 60.0;
sprintf (a, "%c%3d %2d'%5.2f\"", sign, ideg, imin, sec);
return (a);
} /* degstr */
/********************************************************/
int eph_test()
{
char cal;
int p, jday, jmon, jyear;
double al, jd;
centisec *cp;
while (TRUE) {
printf ("date ?");
if (scanf ("%d%d%d", &jday,&jmon,&jyear) < 1) exit(1);
if (jyear < 1600)
cal = 'j';
else
cal = 'g';
swe_date_conversion (jyear, jmon, jday, 0, cal, &jd);
if ((cp = ephread(jd, 0,0, errtext)) == NULL) {
fprintf (stderr,"%s: %s", arg0, errtext);
exit (1);
}
printf ("ephgen test d=%12.1f dmy %d.%d.%d", jd, jday, jmon, jyear);
if (cal == 'g')
printf (" greg");
else
printf (" julian");
printf ("\n\tecliptic %s ", degstr(cp[EP_ECL_INDEX]*CS2DEG));
printf ("nutation %s\n", degstr(cp[EP_NUT_INDEX] * CS2DEG));
for (p = 0; p <= CHIRON; p++) {
al = cp[p] * CS2DEG;
printf ("%2d%18s\n", p, degstr(al));
}
}
} /* end ephtest */
int main(int argc, char **argv)
{
int day, i, n, p;
char serr[AS_MAXCH];
double l[EPHR_NPL][NDB], ecliptic[NDB], nutation[NDB];
double jd0, jd;
double x[6];
int32 jlong;
int file;
int nfiles = 1;
int fnr = -10000;
int32 iflagret;
arg0 = argv[0];
for (i = 1; i < argc; i++) {
if (strncmp(argv[i], "-f", 2) == 0) {
fnr = atoi (argv[i] + 2);
if (fnr < -20 || fnr > 300) {
printf("file number out of range -20 ... 300");
exit (1);
}
}
if (strncmp(argv[i], "-n", 2) == 0) {
nfiles = atoi (argv[i] + 2);
}
if (strncmp(argv[i], "-v", 2) == 0) {
verbose = TRUE;
}
if (strncmp(argv[i], "-t", 2) == 0) {
eph_test();
exit(0);
}
}
if (fnr == -10000) {
fprintf(stderr,"missing file number -fNNN\n");
exit(1);
}
for (file = fnr; file < fnr + nfiles; file++) {
if (file > fnr) printf ("\n");
printf ("file = %d\n", file);
jd0 = EP4_NDAYS * file + 0.5;
jlong = floor(jd0);
if (eph4_posit (jlong, TRUE, errtext) != OK) {
fprintf (stderr,"%s: %s", arg0, errtext);
exit(1);
}
for (n = 0; n < EP4_NDAYS; n += NDB, jd0 += NDB) {
if (n % 500 == 0) {
if ( n > 0 && verbose) {
printf ("\ndd");
for (p = 0; p < 11; p++) {
printf("%6d ",max_dd[p]);
max_dd[p] = 0;
}
printf("\ner");
for (p = 0; p < 11; p++) {
printf("%6.3f ",max_err[p] * 3600);
max_err[p] = 0;
}
}
printf ("\n%d ", n);
} else {
printf (".");
}
fflush( stdout );
for (day = 0; day < NDB; day++) { /* compute positions for 10 days */
jd = jd0 + day;
for (p = SUN; p <= EP_CALC_N; p++) {
if ((iflagret = swe_calc(jd, plac2swe(p), 0, x, serr)) == ERR) {
swe_close();
printf("error in swe_calc() %s\n", serr);
exit (1);
}
l[p][day] = x[0];
}
if ((iflagret = swe_calc(jd, SE_ECL_NUT, 0, x, serr)) == ERR) {
swe_close();
printf("error in swe_calc() %s\n", serr);
exit (1);
}
ecliptic[day] = x[0];
nutation[day] = x[2];
}
jlong = floor(jd0);
eph4_pack (jlong, l, ecliptic, nutation);
}
putchar('\n');
fclose (ephfp);
ephfp = NULL;
} /* for file */
swe_close();
return(0);
} /* end main */

6437
src/swecl.c Normal file

File diff suppressed because it is too large Load Diff

589
src/swedate.c Normal file
View File

@ -0,0 +1,589 @@
/*********************************************************
version 15-feb-89 16:30
swe_date_conversion()
swe_revjul()
swe_julday()
************************************************************/
/* 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.
*/
/*
swe_date_conversion():
This function converts some date+time input {y,m,d,uttime}
into the Julian day number tjd.
The function checks that the input is a legal combination
of dates; for illegal dates like 32 January 1993 it returns ERR
but still converts the date correctly, i.e. like 1 Feb 1993.
The function is usually used to convert user input of birth data
into the Julian day number. Illegal dates should be notified to the user.
Be aware that we always use astronomical year numbering for the years
before Christ, not the historical year numbering.
Astronomical years are done with negative numbers, historical
years with indicators BC or BCE (before common era).
Year 0 (astronomical) = 1 BC historical.
year -1 (astronomical) = 2 BC
etc.
Many users of Astro programs do not know about this difference.
Return: OK or ERR (for illegal date)
*********************************************************/
# include "swephexp.h"
# include "sweph.h"
static TLS AS_BOOL init_leapseconds_done = FALSE;
int CALL_CONV swe_date_conversion(int y,
int m,
int d, /* day, month, year */
double uttime, /* UT in hours (decimal) */
char c, /* calendar g[regorian]|j[ulian] */
double *tjd)
{
int rday, rmon, ryear;
double rut, jd;
int gregflag = SE_JUL_CAL;
if (c == 'g')
gregflag = SE_GREG_CAL;
rut = uttime; /* hours UT */
jd = swe_julday(y, m, d, rut, gregflag);
swe_revjul(jd, gregflag, &ryear, &rmon, &rday, &rut);
*tjd = jd;
if (rmon == m && rday == d && ryear == y) {
return OK;
} else {
return ERR;
}
}
/*************** swe_julday ********************************************
* This function returns the absolute Julian day number (JD)
* for a given calendar date.
* The arguments are a calendar date: day, month, year as integers,
* hour as double with decimal fraction.
* If gregflag = SE_GREG_CAL (1), Gregorian calendar is assumed,
* if gregflag = SE_JUL_CAL (0),Julian calendar is assumed.
The Julian day number is a system of numbering all days continously
within the time range of known human history. It should be familiar
to every astrological or astronomical programmer. The time variable
in astronomical theories is usually expressed in Julian days or
Julian centuries (36525 days per century) relative to some start day;
the start day is called 'the epoch'.
The Julian day number is a double representing the number of
days since JD = 0.0 on 1 Jan -4712, 12:00 noon (in the Julian calendar).
Midnight has always a JD with fraction .5, because traditionally
the astronomical day started at noon. This was practical because
then there was no change of date during a night at the telescope.
From this comes also the fact the noon ephemerides were printed
before midnight ephemerides were introduced early in the 20th century.
NOTE: The Julian day number must not be confused with the Julian
calendar system.
Be aware the we always use astronomical year numbering for the years
before Christ, not the historical year numbering.
Astronomical years are done with negative numbers, historical
years with indicators BC or BCE (before common era).
Year 0 (astronomical) = 1 BC
year -1 (astronomical) = 2 BC
etc.
Original author: Marc Pottenger, Los Angeles.
with bug fix for year < -4711 15-aug-88 by Alois Treindl
References: Oliver Montenbruck, Grundlagen der Ephemeridenrechnung,
Verlag Sterne und Weltraum (1987), p.49 ff
related functions: swe_revjul() reverse Julian day number: compute the
calendar date from a given JD
swe_date_conversion() includes test for legal date values
and notifies errors like 32 January.
****************************************************************/
double CALL_CONV swe_julday(int year, int month, int day, double hour, int gregflag)
{
double jd;
double u,u0,u1,u2;
u = year;
if (month < 3) u -=1;
u0 = u + 4712.0;
u1 = month + 1.0;
if (u1 < 4) u1 += 12.0;
jd = floor(u0*365.25)
+ floor(30.6*u1+0.000001)
+ day + hour/24.0 - 63.5;
if (gregflag == SE_GREG_CAL) {
u2 = floor(fabs(u) / 100) - floor(fabs(u) / 400);
if (u < 0.0) u2 = -u2;
jd = jd - u2 + 2;
if ((u < 0.0) && (u/100 == floor(u/100)) && (u/400 != floor(u/400)))
jd -=1;
}
return jd;
}
/*** swe_revjul ******************************************************
swe_revjul() is the inverse function to swe_julday(), see the description
there.
Arguments are julian day number, calendar flag (0=julian, 1=gregorian)
return values are the calendar day, month, year and the hour of
the day with decimal fraction (0 .. 23.999999).
Be aware the we use astronomical year numbering for the years
before Christ, not the historical year numbering.
Astronomical years are done with negative numbers, historical
years with indicators BC or BCE (before common era).
Year 0 (astronomical) = 1 BC historical year
year -1 (astronomical) = 2 BC historical year
year -234 (astronomical) = 235 BC historical year
etc.
Original author Mark Pottenger, Los Angeles.
with bug fix for year < -4711 16-aug-88 Alois Treindl
*************************************************************************/
void CALL_CONV swe_revjul (double jd, int gregflag,
int *jyear, int *jmon, int *jday, double *jut)
{
double u0,u1,u2,u3,u4;
u0 = jd + 32082.5;
if (gregflag == SE_GREG_CAL) {
u1 = u0 + floor (u0/36525.0) - floor (u0/146100.0) - 38.0;
if (jd >= 1830691.5) u1 +=1;
u0 = u0 + floor (u1/36525.0) - floor (u1/146100.0) - 38.0;
}
u2 = floor (u0 + 123.0);
u3 = floor ( (u2 - 122.2) / 365.25);
u4 = floor ( (u2 - floor (365.25 * u3) ) / 30.6001);
*jmon = (int) (u4 - 1.0);
if (*jmon > 12) *jmon -= 12;
*jday = (int) (u2 - floor (365.25 * u3) - floor (30.6001 * u4));
*jyear = (int) (u3 + floor ( (u4 - 2.0) / 12.0) - 4800);
*jut = (jd - floor (jd + 0.5) + 0.5) * 24.0;
}
/* transform local time to UTC or UTC to local time
*
* input
* iyear ... dsec date and time
* d_timezone timezone offset
* output
* iyear_out ... dsec_out
*
* For time zones east of Greenwich, d_timezone is positive.
* For time zones west of Greenwich, d_timezone is negative.
*
* For conversion from local time to utc, use +d_timezone.
* For conversion from utc to local time, use -d_timezone.
*/
void CALL_CONV swe_utc_time_zone(
int32 iyear, int32 imonth, int32 iday,
int32 ihour, int32 imin, double dsec,
double d_timezone,
int32 *iyear_out, int32 *imonth_out, int32 *iday_out,
int32 *ihour_out, int32 *imin_out, double *dsec_out
)
{
double tjd, d;
AS_BOOL have_leapsec = FALSE;
double dhour;
if (dsec >= 60.0) {
have_leapsec = TRUE;
dsec -= 1.0;
}
dhour = ((double) ihour) + ((double) imin) / 60.0 + dsec / 3600.0;
tjd = swe_julday(iyear, imonth, iday, 0, SE_GREG_CAL);
dhour -= d_timezone;
if (dhour < 0.0) {
tjd -= 1.0;
dhour += 24.0;
}
if (dhour >= 24.0) {
tjd += 1.0;
dhour -= 24.0;
}
swe_revjul(tjd + 0.001, SE_GREG_CAL, iyear_out, imonth_out, iday_out, &d);
*ihour_out = (int) dhour;
d = (dhour - (double) *ihour_out) * 60;
*imin_out = (int) d;
*dsec_out = (d - (double) *imin_out) * 60;
if (have_leapsec)
*dsec_out += 1.0;
}
/*
* functions for the handling of UTC
*/
/* Leap seconds were inserted at the end of the following days:*/
#define NLEAP_SECONDS 27 // ignoring end mark '0'
#define NLEAP_SECONDS_SPACE 100
static TLS int leap_seconds[NLEAP_SECONDS_SPACE] = {
19720630,
19721231,
19731231,
19741231,
19751231,
19761231,
19771231,
19781231,
19791231,
19810630,
19820630,
19830630,
19850630,
19871231,
19891231,
19901231,
19920630,
19930630,
19940630,
19951231,
19970630,
19981231,
20051231,
20081231,
20120630,
20150630,
20161231,
0 /* keep this 0 as end mark */
};
#define J1972 2441317.5
#define NLEAP_INIT 10
/* Read additional leap second dates from external file, if given.
*/
static int init_leapsec(void)
{
FILE *fp;
int ndat, ndat_last;
int tabsiz = 0;
int i;
char s[AS_MAXCH];
char *sp;
if (!init_leapseconds_done) {
init_leapseconds_done = TRUE;
tabsiz = NLEAP_SECONDS;
ndat_last = leap_seconds[NLEAP_SECONDS - 1];
/* no error message if file is missing */
if ((fp = swi_fopen(-1, "seleapsec.txt", swed.ephepath, NULL)) == NULL)
return NLEAP_SECONDS;
while(fgets(s, AS_MAXCH, fp) != NULL) {
sp = s;
while (*sp == ' ' || *sp == '\t') sp++;
sp++;
if (*sp == '#' || *sp == '\n')
continue;
ndat = atoi(s);
if (ndat <= ndat_last)
continue;
/* table space is limited. no error msg, if exceeded */
if (tabsiz >= NLEAP_SECONDS_SPACE)
return tabsiz;
leap_seconds[tabsiz] = ndat;
tabsiz++;
}
if (tabsiz > NLEAP_SECONDS) leap_seconds[tabsiz] = 0; /* end mark */
fclose(fp);
return tabsiz;
}
/* find table size */
tabsiz = 0;
for (i = 0; i < NLEAP_SECONDS_SPACE; i++) {
if (leap_seconds[i] == 0)
break;
else
tabsiz++;
}
return tabsiz;
}
/*
* Input: Clock time UTC, year, month, day, hour, minute, second (decimal).
* gregflag Calendar flag
* serr error string
* Output: An array of doubles:
* dret[0] = Julian day number TT (ET)
* dret[1] = Julian day number UT1
*
* Function returns OK or Error.
*
* - Before 1972, swe_utc_to_jd() treats its input time as UT1.
* Note: UTC was introduced in 1961. From 1961 - 1971, the length of the
* UTC second was regularly changed, so that UTC remained very close to UT1.
* - From 1972 on, input time is treated as UTC.
* - If delta_t - nleap - 32.184 > 1, the input time is treated as UT1.
* Note: Like this we avoid errors greater than 1 second in case that
* the leap seconds table (or the Swiss Ephemeris version) is not updated
* for a long time.
*/
int32 CALL_CONV swe_utc_to_jd(int32 iyear, int32 imonth, int32 iday, int32 ihour, int32 imin, double dsec, int32 gregflag, double *dret, char *serr)
{
double tjd_ut1, tjd_et, tjd_et_1972, dhour, d;
int iyear2, imonth2, iday2;
int i, j, ndat, nleap, tabsiz_nleap;
/*
* error handling: invalid iyear etc.
*/
tjd_ut1 = swe_julday(iyear, imonth, iday, 0, gregflag);
swe_revjul(tjd_ut1, gregflag, &iyear2, &imonth2, &iday2, &d);
if (iyear != iyear2 || imonth != imonth2 || iday != iday2) {
if (serr != NULL)
sprintf(serr, "invalid date: year = %d, month = %d, day = %d", iyear, imonth, iday);
return ERR;
}
if (ihour < 0 || ihour > 23
|| imin < 0 || imin > 59
|| dsec < 0 || dsec >= 61
|| (dsec >= 60 && (imin < 59 || ihour < 23 || tjd_ut1 < J1972))) {
if (serr != NULL)
sprintf(serr, "invalid time: %d:%d:%.2f", ihour, imin, dsec);
return ERR;
}
dhour = (double) ihour + ((double) imin) / 60.0 + dsec / 3600.0;
/*
* before 1972, we treat input date as UT1
*/
if (tjd_ut1 < J1972) {
dret[1] = swe_julday(iyear, imonth, iday, dhour, gregflag);
dret[0] = dret[1] + swe_deltat_ex(dret[1], -1, NULL);
return OK;
}
/*
* if gregflag = Julian calendar, convert to gregorian calendar
*/
if (gregflag == SE_JUL_CAL) {
gregflag = SE_GREG_CAL;
swe_revjul(tjd_ut1, gregflag, &iyear, &imonth, &iday, &d);
}
/*
* number of leap seconds since 1972:
*/
tabsiz_nleap = init_leapsec();
nleap = NLEAP_INIT; /* initial difference between UTC and TAI in 1972 */
ndat = iyear * 10000 + imonth * 100 + iday;
for (i = 0; i < tabsiz_nleap; i++) {
if (ndat <= leap_seconds[i])
break;
nleap++;
}
/*
* For input dates > today:
* If leap seconds table is not up to date, we'd better interpret the
* input time as UT1, not as UTC. How do we find out?
* Check, if delta_t - nleap - 32.184 > 0.9
*/
d = swe_deltat_ex(tjd_ut1, -1, NULL) * 86400.0;
if (d - (double) nleap - 32.184 >= 1.0) {
dret[1] = tjd_ut1 + dhour / 24.0;
dret[0] = dret[1] + swe_deltat_ex(dret[1], -1, NULL);
return OK;
}
/*
* if input second is 60: is it a valid leap second ?
*/
if (dsec >= 60) {
j = 0;
for (i = 0; i < tabsiz_nleap; i++) {
if (ndat == leap_seconds[i]) {
j = 1;
break;
}
}
if (j != 1) {
if (serr != NULL)
sprintf(serr, "invalid time (no leap second!): %d:%d:%.2f", ihour, imin, dsec);
return ERR;
}
}
/*
* convert UTC to ET and UT1
*/
/* the number of days between input date and 1 jan 1972: */
d = tjd_ut1 - J1972;
/* SI time since 1972, ignoring leap seconds: */
d += (double) ihour / 24.0 + (double) imin / 1440.0 + dsec / 86400.0;
/* ET (TT) */
tjd_et_1972 = J1972 + (32.184 + NLEAP_INIT) / 86400.0;
tjd_et = tjd_et_1972 + d + ((double) (nleap - NLEAP_INIT)) / 86400.0;
d = swe_deltat_ex(tjd_et, -1, NULL);
tjd_ut1 = tjd_et - swe_deltat_ex(tjd_et - d, -1, NULL);
tjd_ut1 = tjd_et - swe_deltat_ex(tjd_ut1, -1, NULL);
dret[0] = tjd_et;
dret[1] = tjd_ut1;
return OK;
}
/*
* Input: tjd_et Julian day number, terrestrial time (ephemeris time).
* gregfalg Calendar flag
* Output: UTC year, month, day, hour, minute, second (decimal).
*
* - Before 1 jan 1972 UTC, output UT1.
* Note: UTC was introduced in 1961. From 1961 - 1971, the length of the
* UTC second was regularly changed, so that UTC remained very close to UT1.
* - From 1972 on, output is UTC.
* - If delta_t - nleap - 32.184 > 1, the output is UT1.
* Note: Like this we avoid errors greater than 1 second in case that
* the leap seconds table (or the Swiss Ephemeris version) has not been
* updated for a long time.
*/
void CALL_CONV swe_jdet_to_utc(double tjd_et, int32 gregflag, int32 *iyear, int32 *imonth, int32 *iday, int32 *ihour, int32 *imin, double *dsec)
{
int i;
int second_60 = 0;
int iyear2, imonth2, iday2, nleap, ndat, tabsiz_nleap;
double d, tjd, tjd_et_1972, tjd_ut, dret[10];
/*
* if tjd_et is before 1 jan 1972 UTC, return UT1
*/
tjd_et_1972 = J1972 + (32.184 + NLEAP_INIT) / 86400.0;
d = swe_deltat_ex(tjd_et, -1, NULL);
tjd_ut = tjd_et - swe_deltat_ex(tjd_et - d, -1, NULL);
tjd_ut = tjd_et - swe_deltat_ex(tjd_ut, -1, NULL);
if (tjd_et < tjd_et_1972) {
swe_revjul(tjd_ut, gregflag, iyear, imonth, iday, &d);
*ihour = (int32) d;
d -= (double) *ihour;
d *= 60;
*imin = (int32) d;
*dsec = (d - (double) *imin) * 60.0;
return;
}
/*
* minimum number of leap seconds since 1972; we may be missing one leap
* second
*/
tabsiz_nleap = init_leapsec();
swe_revjul(tjd_ut-1, SE_GREG_CAL, &iyear2, &imonth2, &iday2, &d);
ndat = iyear2 * 10000 + imonth2 * 100 + iday2;
nleap = 0;
for (i = 0; i < tabsiz_nleap; i++) {
if (ndat <= leap_seconds[i])
break;
nleap++;
}
/* date of potentially missing leapsecond */
if (nleap < tabsiz_nleap) {
i = leap_seconds[nleap];
iyear2 = i / 10000;
imonth2 = (i % 10000) / 100;;
iday2 = i % 100;
tjd = swe_julday(iyear2, imonth2, iday2, 0, SE_GREG_CAL);
swe_revjul(tjd+1, SE_GREG_CAL, &iyear2, &imonth2, &iday2, &d);
swe_utc_to_jd(iyear2,imonth2,iday2, 0, 0, 0, SE_GREG_CAL, dret, NULL);
d = tjd_et - dret[0];
if (d >= 0) {
nleap++;
} else if (d < 0 && d > -1.0/86400.0) {
second_60 = 1;
}
}
/*
* UTC, still unsure about one leap second
*/
tjd = J1972 + (tjd_et - tjd_et_1972) - ((double) nleap + second_60) / 86400.0;
swe_revjul(tjd, SE_GREG_CAL, iyear, imonth, iday, &d);
*ihour = (int32) d;
d -= (double) *ihour;
d *= 60;
*imin = (int32) d;
*dsec = (d - (double) *imin) * 60.0 + second_60;
/*
* For input dates > today:
* If leap seconds table is not up to date, we'd better interpret the
* input time as UT1, not as UTC. How do we find out?
* Check, if delta_t - nleap - 32.184 > 0.9
*/
d = swe_deltat_ex(tjd_et, -1, NULL);
d = swe_deltat_ex(tjd_et - d, -1, NULL);
if (d * 86400.0 - (double) (nleap + NLEAP_INIT) - 32.184 >= 1.0) {
swe_revjul(tjd_et - d, SE_GREG_CAL, iyear, imonth, iday, &d);
*ihour = (int32) d;
d -= (double) *ihour;
d *= 60;
*imin = (int32) d;
*dsec = (d - (double) *imin) * 60.0;
}
if (gregflag == SE_JUL_CAL) {
tjd = swe_julday(*iyear, *imonth, *iday, 0, SE_GREG_CAL);
swe_revjul(tjd, gregflag, iyear, imonth, iday, &d);
}
}
/*
* Input: tjd_ut Julian day number, universal time (UT1).
* gregfalg Calendar flag
* Output: UTC year, month, day, hour, minute, second (decimal).
*
* - Before 1 jan 1972 UTC, output UT1.
* Note: UTC was introduced in 1961. From 1961 - 1971, the length of the
* UTC second was regularly changed, so that UTC remained very close to UT1.
* - From 1972 on, output is UTC.
* - If delta_t - nleap - 32.184 > 1, the output is UT1.
* Note: Like this we avoid errors greater than 1 second in case that
* the leap seconds table (or the Swiss Ephemeris version) has not been
* updated for a long time.
*/
void CALL_CONV swe_jdut1_to_utc(double tjd_ut, int32 gregflag, int32 *iyear, int32 *imonth, int32 *iday, int32 *ihour, int32 *imin, double *dsec)
{
double tjd_et = tjd_ut + swe_deltat_ex(tjd_ut, -1, NULL);
swe_jdet_to_utc(tjd_et, gregflag, iyear, imonth, iday, ihour, imin, dsec);
}

81
src/swedate.h Normal file
View File

@ -0,0 +1,81 @@
/*********************************************************
version 15-feb-89 16:30
*********************************************************/
/* 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.
*/
#ifdef __cplusplus
extern "C" {
#endif
#ifndef _SWEDLL_H
extern EXP32 int swe_date_conversion (
int y , int m , int d , /* year, month, day */
double utime, /* universal time in hours (decimal) */
char c, /* calendar g[regorian]|j[ulian]|a[stro = greg] */
double *tgmt);
extern EXP32 double *swe_julday(
int year, int month, int day, double hour,
int gregflag);
extern EXP32 void swe_revjul (
double jd,
int gregflag,
int *jyear, int *jmon, int *jday, double *jut);
#endif
#ifdef __cplusplus
} /* extern C */
#endif

388
src/swedll.h Normal file
View File

@ -0,0 +1,388 @@
/* SWISSEPH
*
* Windows DLL interface imports for the Astrodienst SWISSEPH package
*
**************************************************************/
/* 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.
*/
/* $Id: swedll.h,v 1.75 2009/04/08 07:19:08 dieter Exp $ */
#ifdef __cplusplus
extern "C" {
#endif
#ifndef _SWEDLL_H
#define _SWEDLL_H
#ifndef _SWEPHEXP_INCLUDED
#include "swephexp.h"
#endif
# ifdef __cplusplus
#define DllImport extern "C" __declspec( dllimport )
# else
#define DllImport __declspec( dllimport )
# endif
/* DLL defines
Define UNDECO_DLL for un-decorated dll
verify compiler option __cdecl for un-decorated and __stdcall for decorated */
/*#define UNDECO_DLL */
#if defined (PASCAL) || defined(__stdcall)
#if defined UNDECO_DLL
#define CALL_CONV_IMP __cdecl
#else
#define CALL_CONV_IMP __stdcall
#endif
#else
#define CALL_CONV_IMP
#endif
DllImport int32 CALL_CONV_IMP swe_heliacal_ut(double JDNDaysUTStart, double *geopos, double *datm, double *dobs, char *ObjectName, int32 TypeEvent, int32 iflag, double *dret, char *serr);
DllImport int32 CALL_CONV_IMP swe_heliacal_pheno_ut(double JDNDaysUT, double *geopos, double *datm, double *dobs, char *ObjectName, int32 TypeEvent, int32 helflag, double *darr, char *serr);
DllImport int32 CALL_CONV_IMP swe_vis_limit_mag(double tjdut, double *geopos, double *datm, double *dobs, char *ObjectName, int32 helflag, double *dret, char *serr);
/* the following are secret, for Victor Reijs' */
DllImport int32 CALL_CONV_IMP swe_heliacal_angle(double tjdut, double *dgeo, double *datm, double *dobs, int32 helflag, double mag, double azi_obj, double azi_sun, double azi_moon, double alt_moon, double *dret, char *serr);
DllImport int32 CALL_CONV_IMP swe_topo_arcus_visionis(double tjdut, double *dgeo, double *datm, double *dobs, int32 helflag, double mag, double azi_obj, double alt_obj, double azi_sun, double azi_moon, double alt_moon, double *dret, char *serr);
DllImport double CALL_CONV_IMP swe_degnorm(double deg);
DllImport char * CALL_CONV_IMP swe_version(char *);
DllImport char * CALL_CONV_IMP swe_get_library_path(char *);
DllImport int32 CALL_CONV_IMP swe_calc(
double tjd, int ipl, int32 iflag,
double *xx,
char *serr);
DllImport int32 CALL_CONV_IMP swe_calc_pctr(
double tjd, int32 ipl, int32 iplctr, int32 iflag,
double *xxret,
char *serr);
DllImport int32 CALL_CONV_IMP swe_calc_ut(
double tjd_ut, int32 ipl, int32 iflag,
double *xx,
char *serr);
DllImport int32 CALL_CONV_IMP swe_fixstar(
char *star, double tjd, int32 iflag,
double *xx,
char *serr);
DllImport int32 CALL_CONV_IMP swe_fixstar_ut(
char *star, double tjd_ut, int32 iflag,
double *xx,
char *serr);
DllImport int32 CALL_CONV_IMP swe_fixstar_mag(
char *star, double *xx, char *serr);
DllImport int32 CALL_CONV_IMP swe_fixstar2(
char *star, double tjd, int32 iflag,
double *xx,
char *serr);
DllImport int32 CALL_CONV_IMP swe_fixstar2_ut(
char *star, double tjd_ut, int32 iflag,
double *xx,
char *serr);
DllImport int32 CALL_CONV_IMP swe_fixstar2_mag(
char *star, double *xx, char *serr);
DllImport double CALL_CONV_IMP swe_sidtime0(double tjd_ut, double ecl, double nut);
DllImport double CALL_CONV_IMP swe_sidtime(double tjd_ut);
DllImport double CALL_CONV_IMP swe_deltat_ex(double tjd, int32 iflag, char *serr);
DllImport double CALL_CONV_IMP swe_deltat(double tjd);
DllImport int CALL_CONV_IMP swe_houses(
double tjd_ut, double geolat, double geolon, int hsys,
double *hcusps, double *ascmc);
DllImport int CALL_CONV_IMP swe_houses_ex(
double tjd_ut, int32 iflag, double geolat, double geolon, int hsys,
double *hcusps, double *ascmc);
DllImport int CALL_CONV_IMP swe_houses_ex2(
double tjd_ut, int32 iflag, double geolat, double geolon, int hsys,
double *hcusps, double *ascmc, double *cusp_speed, double *ascmc_speed, char *serr);
DllImport int CALL_CONV_IMP swe_houses_armc(
double armc, double geolat, double eps, int hsys,
double *hcusps, double *ascmc);
DllImport int CALL_CONV_IMP swe_houses_armc_ex2(
double armc, double geolat, double eps, int hsys,
double *hcusps, double *ascmc, double *cusp_speed, double *ascmc_speed, char *serr);
DllImport double CALL_CONV_IMP swe_house_pos(
double armc, double geolon, double eps, int hsys, double *xpin, char *serr);
DllImport char * CALL_CONV_IMP swe_house_name(int hsys);
DllImport int32 CALL_CONV_IMP swe_gauquelin_sector(
double t_ut, int32 ipl, char *starname, int32 iflag, int32 imeth, double *geopos, double atpress, double attemp, double *dgsect, char *serr);
DllImport void CALL_CONV_IMP swe_set_sid_mode(
int32 sid_mode, double t0, double ayan_t0);
DllImport int32 CALL_CONV_IMP swe_get_ayanamsa_ex(double tjd_et, int32 iflag, double *daya, char *serr);
DllImport int32 CALL_CONV_IMP swe_get_ayanamsa_ex_ut(double tjd_ut, int32 iflag, double *daya, char *serr);
DllImport double CALL_CONV_IMP swe_get_ayanamsa(double tjd_et);
DllImport double CALL_CONV_IMP swe_get_ayanamsa_ut(double tjd_ut);
DllImport char * CALL_CONV_IMP swe_get_ayanamsa_name(int32 isidmode);
DllImport char * CALL_CONV_IMP swe_get_current_file_data(int ifno, double *tfstart, double *tfend, int *denum);
DllImport int CALL_CONV_IMP swe_date_conversion(
int y , int m , int d , /* year, month, day */
double utime, /* universal time in hours (decimal) */
char c, /* calendar g[regorian]|j[ulian]|a[stro = greg] */
double *tjd);
DllImport double CALL_CONV_IMP swe_julday(
int year, int mon, int mday,
double hour,
int gregflag);
DllImport void CALL_CONV_IMP swe_revjul(
double jd, int gregflag,
int *year, int *mon, int *mday,
double *hour);
DllImport void CALL_CONV_IMP swe_utc_time_zone(
int32 iyear, int32 imonth, int32 iday,
int32 ihour, int32 imin, double dsec,
double d_timezone,
int32 *iyear_out, int32 *imonth_out, int32 *iday_out,
int32 *ihour_out, int32 *imin_out, double *dsec_out);
DllImport int32 CALL_CONV_IMP swe_utc_to_jd(
int32 iyear, int32 imonth, int32 iday,
int32 ihour, int32 imin, double dsec,
int32 gregflag, double *dret, char *serr);
DllImport void CALL_CONV_IMP swe_jdet_to_utc(
double tjd_et, int32 gregflag,
int32 *iyear, int32 *imonth, int32 *iday,
int32 *ihour, int32 *imin, double *dsec);
DllImport void CALL_CONV_IMP swe_jdut1_to_utc(
double tjd_ut, int32 gregflag,
int32 *iyear, int32 *imonth, int32 *iday,
int32 *ihour, int32 *imin, double *dsec);
DllImport int CALL_CONV_IMP swe_time_equ(
double tjd, double *e, char *serr);
DllImport int CALL_CONV_IMP swe_lmt_to_lat(double tjd_lmt, double geolon, double *tjd_lat, char *serr);
DllImport int CALL_CONV_IMP swe_lat_to_lmt(double tjd_lat, double geolon, double *tjd_lmt, char *serr);
DllImport double CALL_CONV_IMP swe_get_tid_acc(void);
DllImport void CALL_CONV_IMP swe_set_tid_acc(double tidacc);
DllImport void CALL_CONV_IMP swe_set_delta_t_userdef(double dt);
DllImport void CALL_CONV_IMP swe_set_ephe_path(char *path);
DllImport void CALL_CONV_IMP swe_set_jpl_file(char *fname);
DllImport void CALL_CONV_IMP swe_close(void);
DllImport char * CALL_CONV_IMP swe_get_planet_name(int ipl, char *spname);
DllImport void CALL_CONV_IMP swe_cotrans(double *xpo, double *xpn, double eps);
DllImport void CALL_CONV_IMP swe_cotrans_sp(double *xpo, double *xpn, double eps);
DllImport void CALL_CONV_IMP swe_set_topo(double geolon, double geolat, double height);
DllImport void CALL_CONV_IMP swe_set_astro_models(char *samod, int32 iflag);
DllImport void CALL_CONV_IMP swe_get_astro_models(char *samod, char *sdet, int32 iflag);
/****************************
* from swecl.c
****************************/
/* computes geographic location and attributes of solar
* eclipse at a given tjd */
DllImport int32 CALL_CONV_IMP swe_sol_eclipse_where(double tjd, int32 ifl, double *geopos, double *attr, char *serr);
DllImport int32 CALL_CONV_IMP swe_lun_occult_where(double tjd, int32 ipl, char *starname, int32 ifl, double *geopos, double *attr, char *serr);
/* computes attributes of a solar eclipse for given tjd, geolon, geolat */
DllImport int32 CALL_CONV_IMP swe_sol_eclipse_how(double tjd, int32 ifl, double *geopos, double *attr, char *serr);
/* finds time of next local eclipse */
DllImport int32 CALL_CONV_IMP swe_sol_eclipse_when_loc(double tjd_start, int32 ifl, double *geopos, double *tret, double *attr, int32 backward, char *serr);
DllImport int32 CALL_CONV_IMP swe_lun_occult_when_loc(double tjd_start, int32 ipl, char *starname, int32 ifl, double *geopos, double *tret, double *attr, int32 backward, char *serr);
/* finds time of next eclipse globally */
DllImport int32 CALL_CONV_IMP swe_sol_eclipse_when_glob(double tjd_start, int32 ifl, int32 ifltype, double *tret, int32 backward, char *serr);
/* finds time of next occultation globally */
DllImport int32 CALL_CONV_IMP swe_lun_occult_when_glob(double tjd_start, int32 ipl, char *starname, int32 ifl, int32 ifltype, double *tret, int32 backward, char *serr);
/* computes attributes of a lunar eclipse for given tjd */
DllImport int32 CALL_CONV_IMP swe_lun_eclipse_how(
double tjd_ut,
int32 ifl,
double *geopos,
double *attr,
char *serr);
DllImport int32 CALL_CONV_IMP swe_lun_eclipse_when(double tjd_start, int32 ifl, int32 ifltype, double *tret, int32 backward, char *serr);
DllImport int32 CALL_CONV_IMP swe_lun_eclipse_when_loc(double tjd_start, int32 ifl, double *geopos, double *tret, double *attr, int32 backward, char *serr);
/* planetary phenomena */
DllImport int32 CALL_CONV_IMP swe_pheno(double tjd, int32 ipl, int32 iflag, double *attr, char *serr);
DllImport int32 CALL_CONV_IMP swe_pheno_ut(double tjd_ut, int32 ipl, int32 iflag, double *attr, char *serr);
DllImport double CALL_CONV_IMP swe_refrac(double inalt, double atpress, double attemp, int32 calc_flag);
DllImport double CALL_CONV_IMP swe_refrac_extended(double inalt, double geoalt, double atpress, double attemp, double lapse_rate, int32 calc_flag, double *dret);
DllImport void CALL_CONV_IMP swe_set_lapse_rate(double lapse_rate);
DllImport void CALL_CONV_IMP swe_azalt(
double tjd_ut,
int32 calc_flag,
double *geopos,
double atpress,
double attemp,
double *xin,
double *xaz);
DllImport void CALL_CONV_IMP swe_azalt_rev(
double tjd_ut,
int32 calc_flag,
double *geopos,
double *xin,
double *xout);
DllImport int32 CALL_CONV_IMP swe_rise_trans(
double tjd_ut, int32 ipl, char *starname,
int32 epheflag, int32 rsmi,
double *geopos,
double atpress, double attemp,
double *tret,
char *serr);
DllImport int32 CALL_CONV_IMP swe_rise_trans_true_hor(
double tjd_ut, int32 ipl, char *starname,
int32 epheflag, int32 rsmi,
double *geopos,
double atpress, double attemp,
double horhgt,
double *tret,
char *serr);
DllImport int32 CALL_CONV_IMP swe_nod_aps(double tjd_et, int32 ipl, int32 iflag,
int32 method,
double *xnasc, double *xndsc,
double *xperi, double *xaphe,
char *serr);
DllImport int32 CALL_CONV_IMP swe_nod_aps_ut(double tjd_ut, int32 ipl, int32 iflag,
int32 method,
double *xnasc, double *xndsc,
double *xperi, double *xaphe,
char *serr);
DllImport int32 CALL_CONV_IMP swe_get_orbital_elements(double tjd_et, int32 ipl, int32 iflag, double *dret, char *serr);
DllImport int32 CALL_CONV_IMP swe_orbit_max_min_true_distance(double tjd_et, int32 ipl, int32 iflag, double *dmax, double *dmin, double *dtrue, char *serr);
/*******************************************************
* other functions from swephlib.c;
* they are not needed for Swiss Ephemeris,
* but may be useful to former Placalc users.
********************************************************/
/* normalize argument into interval [0..DEG360] */
DllImport centisec CALL_CONV_IMP swe_csnorm(centisec p);
/* distance in centisecs p1 - p2 normalized to [0..360[ */
DllImport centisec CALL_CONV_IMP swe_difcsn (centisec p1, centisec p2);
DllImport double CALL_CONV_IMP swe_difdegn (double p1, double p2);
/* distance in centisecs p1 - p2 normalized to [-180..180[ */
DllImport centisec CALL_CONV_IMP swe_difcs2n(centisec p1, centisec p2);
DllImport double CALL_CONV_IMP swe_difdeg2n(double p1, double p2);
DllImport double CALL_CONV_IMP swe_difdeg2n(double p1, double p2);
DllImport double CALL_CONV_IMP swe_difrad2n(double p1, double p2);
DllImport double CALL_CONV_IMP swe_rad_midp(double x1, double x0);
DllImport double CALL_CONV_IMP swe_deg_midp(double x1, double x0);
/* round second, but at 29.5959 always down */
DllImport centisec CALL_CONV_IMP swe_csroundsec(centisec x);
/* double to int32 with rounding, no overflow check */
DllImport int32 CALL_CONV_IMP swe_d2l(double x);
DllImport void CALL_CONV_IMP swe_split_deg(double ddeg, int32 roundflag, int32 *ideg, int32 *imin, int32 *isec, double *dsecfr, int32 *isgn);
/* monday = 0, ... sunday = 6 */
DllImport int CALL_CONV_IMP swe_day_of_week(double jd);
DllImport char * CALL_CONV_IMP swe_cs2timestr(CSEC t, int sep, AS_BOOL suppressZero, char *a);
DllImport char * CALL_CONV_IMP swe_cs2lonlatstr(CSEC t, char pchar, char mchar, char *s);
DllImport char * CALL_CONV_IMP swe_cs2degstr(CSEC t, char *a);
DllImport void CALL_CONV_IMP swe_set_interpolate_nut(AS_BOOL do_interpolate);
#endif /* !_SWEDLL_H */
#ifdef __cplusplus
} /* extern C */
#endif

3518
src/swehel.c Normal file

File diff suppressed because it is too large Load Diff

3047
src/swehouse.c Normal file

File diff suppressed because it is too large Load Diff

98
src/swehouse.h Normal file
View File

@ -0,0 +1,98 @@
/*******************************************************
module swehouse.h
house and (simple) aspect calculation
*******************************************************/
/* 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.
*/
struct houses {
double cusp[37];
double cusp_speed[37];
double ac;
double ac_speed; // speed of ac
double mc;
double mc_speed; // speed of mc
double armc_speed; // speed of armc
double vertex;
double vertex_speed; // speed of vertex
double equasc;
double equasc_speed; // speed
double coasc1;
double coasc1_speed; // speed
double coasc2;
double coasc2_speed; // speed
double polasc;
double polasc_speed; // speed
double sundec; // declination of Sun for Sunshine houses
AS_BOOL do_speed;
AS_BOOL do_hspeed;
AS_BOOL do_interpol;
char serr[AS_MAXCH];
};
#define HOUSES struct houses
#define VERY_SMALL 1E-10
#define degtocs(x) (d2l((x) * DEG))
#define cstodeg(x) (double)((x) * CS2DEG)
#define sind(x) sin((x) * DEGTORAD)
#define cosd(x) cos((x) * DEGTORAD)
#define tand(x) tan((x) * DEGTORAD)
#define asind(x) (asin(x) * RADTODEG)
#define acosd(x) (acos(x) * RADTODEG)
#define atand(x) (atan(x) * RADTODEG)
#define atan2d(y, x) (atan2(y, x) * RADTODEG)

952
src/swejpl.c Normal file
View File

@ -0,0 +1,952 @@
/*
|
| Subroutines for reading JPL ephemerides.
| derived from testeph.f as contained in DE403 distribution July 1995.
| works with DE200, DE102, DE403, DE404, DE405, DE406, DE431
| (attention, these ephemerides do not have exactly the same reference frame)
Authors: Dieter Koch and Alois Treindl, Astrodienst Zurich
************************************************************/
/* 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.
*/
#if MSDOS
#else
#define _FILE_OFFSET_BITS 64
#endif
#include <string.h>
#include "swephexp.h"
#include "sweph.h"
#include "swejpl.h"
#if MSDOS
typedef __int64 off_t64;
#define FSEEK _fseeki64
#define FTELL _ftelli64
#else
typedef off_t off_t64;
#define FSEEK fseeko
#define FTELL ftello
#endif
#define DEBUG_DO_SHOW FALSE
/*
* local globals
*/
struct jpl_save {
char *jplfname;
char *jplfpath;
FILE *jplfptr;
short do_reorder;
double eh_cval[400];
double eh_ss[3], eh_au, eh_emrat;
int32 eh_denum, eh_ncon, eh_ipt[39];
char ch_cnam[6*400];
double pv[78];
double pvsun[6];
double buf[1500];
double pc[18], vc[18], ac[18], jc[18];
short do_km;
};
static TLS struct jpl_save *js;
static int state (double et, int32 *list, int do_bary,
double *pv, double *pvsun, double *nut, char *serr);
static int interp(double *buf, double t, double intv, int32 ncfin,
int32 ncmin, int32 nain, int32 ifl, double *pv);
static int32 fsizer(char *serr);
static void reorder(char *x, int size, int number);
static int read_const_jpl(double *ss, char *serr);
/* information about eh_ipt[] and buf[]
DE200 DE102 DE403
3 3 ipt[0] 3 body 0 (mercury) starts at buf[2]
12 15 ipt[1] 14 body 0, ncf = coefficients per component
4 2 ipt[2] 4 na = nintervals, tot 14*4*3=168
147 93 ipt[3] 171 body 1 (venus) starts at buf[170]
12 15 ipt[4] 10 ncf = coefficients per component
1 1 ipt[5] 2 total 10*2*3=60
183 138 ipt[6] 231 body 2 (earth) starts at buf[230]
15 15 ipt[7] 13 ncf = coefficients per component
2 2 ipt[8] 2 total 13*2*3=78
273 228 ipt[9] 309 body 3 (mars) starts at buf[308]
10 10 ipt[10] 11 ncf = coefficients per component
1 1 ipt[11] 1 total 11*1*3=33
303 258 ipt[12] 342 body 4 (jupiter) at buf[341]
9 9 ipt[13] 8 total 8 * 1 * 3 = 24
1 1 ipt[14] 1
330 285 ipt[15] 366 body 5 (saturn) at buf[365]
8 8 ipt[16] 7 total 7 * 1 * 3 = 21
1 1 ipt[17] 1
354 309 ipt[18] 387 body 6 (uranus) at buf[386]
8 8 ipt[19] 6 total 6 * 1 * 3 = 18
1 1 ipt[20] 1
378 333 ipt[21] 405 body 7 (neptune) at buf[404]
6 6 ipt[22] 6 total 18
1 1 ipt[23] 1
396 351 ipt[24] 423 body 8 (pluto) at buf[422]
6 6 ipt[25] 6 total 18
1 1 ipt[26] 1
414 369 ipt[27] 441 body 9 (moon) at buf[440]
12 15 ipt[28] 13 total 13 * 8 * 3 = 312
8 8 ipt[29] 8
702 729 ipt[30] 753 SBARY SUN, starts at buf[752]
15 15 ipt[31] 11 SBARY SUN, ncf = coeff per component
1 1 ipt[32] 2 total 11*2*3=66
747 774 ipt[33] 819 nutations, starts at buf[818]
10 0 ipt[34] 10 total 10 * 4 * 2 = 80
4 0 ipt[35] 4 (nutation only two coordinates)
0 0 ipt[36] 899 librations, start at buf[898]
0 0 ipt[37] 10 total 10 * 4 * 3 = 120
0 0 ipt[38] 4
last element of buf[1017]
buf[0] contains start jd and buf[1] end jd of segment;
each segment is 32 days in de403, 64 days in DE102, 32 days in DE200
Length of blocks: DE406 = 1456*4=5824 bytes = 728 double
DE405 = 2036*4=8144 bytes = 1018 double
DE404 = 1456*4=5824 bytes = 728 double
DE403 = 2036*4=8144 bytes = 1018 double
DE200 = 1652*4=6608 bytes = 826 double
DE102 = 1546*4=6184 bytes = 773 double
each DE102 record has 53*8=424 fill bytes so that
the records have the same length as DE200.
*/
/*
* This subroutine opens the file jplfname, with a phony record length,
* reads the first record, and uses the info to compute ksize,
* the number of single precision words in a record.
* RETURN: ksize (record size of ephemeris data)
* jplfptr is opened on return.
* note 26-aug-2008: now record size is computed by fsizer(), not
* set to a fixed value depending as in previous releases. The caller of
* fsizer() will verify by data comparison whether it computed correctly.
*/
static int32 fsizer(char *serr)
{
/* Local variables */
int32 ncon;
double emrat;
int32 numde;
double au, ss[3];
int i, kmx, khi, nd;
int32 ksize, lpt[3];
char ttl[6*14*3];
size_t nrd; /* unused, removes compile warnings */
if ((js->jplfptr = swi_fopen(SEI_FILE_PLANET, js->jplfname, js->jplfpath, serr)) == NULL) {
return NOT_AVAILABLE;
}
/* ttl = ephemeris title, e.g.
* "JPL Planetary Ephemeris DE404/LE404
* Start Epoch: JED= 625296.5-3001 DEC 21 00:00:00
* Final Epoch: JED= 2817168.5 3001 JAN 17 00:00:00c */
nrd = fread((void *) &ttl[0], 1, 252, js->jplfptr);
if (nrd != 252) return NOT_AVAILABLE;
/* cnam = names of constants */
nrd = fread((void *) js->ch_cnam, 1, 6*400, js->jplfptr);
if (nrd != 6*400) return NOT_AVAILABLE;
/* ss[0] = start epoch of ephemeris
* ss[1] = end epoch
* ss[2] = segment size in days */
nrd = fread((void *) &ss[0], sizeof(double), 3, js->jplfptr);
if (nrd != 3) return NOT_AVAILABLE;
/* reorder ? */
if (ss[2] < 1 || ss[2] > 200)
js->do_reorder = TRUE;
else
js->do_reorder = 0;
for (i = 0; i < 3; i++)
js->eh_ss[i] = ss[i];
if (js->do_reorder)
reorder((char *) &js->eh_ss[0], sizeof(double), 3);
/* plausibility test of these constants. Start and end date must be
* between -20000 and +20000, segment size >= 1 and <= 200 */
if (js->eh_ss[0] < -5583942 || js->eh_ss[1] > 9025909 || js->eh_ss[2] < 1 || js->eh_ss[2] > 200) {
if (serr != NULL) {
strcpy(serr, "alleged ephemeris file has invalid format.");
if (strlen(serr) + strlen(js->jplfname) + 3 < AS_MAXCH) {
sprintf(serr, "alleged ephemeris file (%s) has invalid format.", js->jplfname);
}
}
return(NOT_AVAILABLE);
}
/* ncon = number of constants */
nrd = fread((void *) &ncon, sizeof(int32), 1, js->jplfptr);
if (nrd != 1) return NOT_AVAILABLE;
if (js->do_reorder)
reorder((char *) &ncon, sizeof(int32), 1);
/* au = astronomical unit */
nrd = fread((void *) &au, sizeof(double), 1, js->jplfptr);
if (nrd != 1) return NOT_AVAILABLE;
if (js->do_reorder)
reorder((char *) &au, sizeof(double), 1);
/* emrat = earth moon mass ratio */
nrd = fread((void *) &emrat, sizeof(double), 1, js->jplfptr);
if (nrd != 1) return NOT_AVAILABLE;
if (js->do_reorder)
reorder((char *) &emrat, sizeof(double), 1);
/* ipt[i+0]: coefficients of planet i start at buf[ipt[i+0]-1]
* ipt[i+1]: number of coefficients (interpolation order - 1)
* ipt[i+2]: number of intervals in segment */
nrd = fread((void *) &js->eh_ipt[0], sizeof(int32), 36, js->jplfptr);
if (nrd != 36) return NOT_AVAILABLE;
if (js->do_reorder)
reorder((char *) &js->eh_ipt[0], sizeof(int32), 36);
/* numde = number of jpl ephemeris "404" with de404 */
nrd = fread((void *) &numde, sizeof(int32), 1, js->jplfptr);
if (nrd != 1) return NOT_AVAILABLE;
if (js->do_reorder)
reorder((char *) &numde, sizeof(int32), 1);
/* read librations */
nrd = fread(&lpt[0], sizeof(int32), 3, js->jplfptr);
if (nrd != 3) return NOT_AVAILABLE;
if (js->do_reorder)
reorder((char *) &lpt[0], sizeof(int32), 3);
/* fill librations into eh_ipt[36]..[38] */
for (i = 0; i < 3; ++i)
js->eh_ipt[i + 36] = lpt[i];
rewind(js->jplfptr);
/* find the number of ephemeris coefficients from the pointers */
/* re-activated this code on 26-aug-2008 */
kmx = 0;
khi = 0;
for (i = 0; i < 13; i++) {
if (js->eh_ipt[i * 3] > kmx) {
kmx = js->eh_ipt[i * 3];
khi = i + 1;
}
}
if (khi == 12)
nd = 2;
else
nd = 3;
ksize = (js->eh_ipt[khi * 3 - 3] + nd * js->eh_ipt[khi * 3 - 2] * js->eh_ipt[khi * 3 - 1] - 1L) * 2L;
/*
* de102 files give wrong ksize, because they contain 424 empty bytes
* per record. Fixed by hand!
*/
if (ksize == 1546)
ksize = 1652;
#if 0 /* we prefer to compute ksize to be comaptible
with new DE releases */
switch (numde) {
case 403:
case 405:
case 410:
case 413:
case 414:
case 418:
case 421:
ksize = 2036;
break;
case 404:
case 406:
ksize = 1456;
break;
case 200:
ksize = 1652;
break;
case 102:
ksize = 1652; /* de102 is filled with blanks to length of de200 */
break;
default:
if (serr != NULL)
sprintf(serr,"unknown numde value %d;", numde);
return ERR;
}
#endif
if (ksize < 1000 || ksize > 5000) {
if (serr != NULL)
sprintf(serr, "JPL ephemeris file does not provide valid ksize (%d)", ksize);/**/
return NOT_AVAILABLE;
}
return ksize;
}
/*
* This subroutine reads the jpl planetary ephemeris
* and gives the position and velocity of the point 'ntarg'
* with respect to 'ncent'.
* calling sequence parameters:
* et = d.p. julian ephemeris date at which interpolation
* is wanted.
* ** note the entry dpleph for a doubly-dimensioned time **
* the reason for this option is discussed in the
* subroutine state
* ntarg = integer number of 'target' point.
* ncent = integer number of center point.
* the numbering convention for 'ntarg' and 'ncent' is:
* 0 = mercury 7 = neptune
* 1 = venus 8 = pluto
* 2 = earth 9 = moon
* 3 = mars 10 = sun
* 4 = jupiter 11 = solar-system barycenter
* 5 = saturn 12 = earth-moon barycenter
* 6 = uranus 13 = nutations (longitude and obliq)
* 14 = librations, if on eph file
* (if nutations are wanted, set ntarg = 13. for librations,
* set ntarg = 14. set ncent=0.)
* rrd = output 6-word d.p. array containing position and velocity
* of point 'ntarg' relative to 'ncent'. the units are au and
* au/day. for librations the units are radians and radians
* per day. in the case of nutations the first four words of
* rrd will be set to nutations and rates, having units of
* radians and radians/day.
* The option is available to have the units in km and km/sec.
* For this, set do_km=TRUE (default FALSE).
*/
int swi_pleph(double et, int ntarg, int ncent, double *rrd, char *serr)
{
int i, retc;
int32 list[12];
double *pv = js->pv;
double *pvsun = js->pvsun;
for (i = 0; i < 6; ++i)
rrd[i] = 0.0;
if (ntarg == ncent)
return 0;
for (i = 0; i < 12; ++i)
list[i] = 0;
/* check for nutation call */
if (ntarg == J_NUT) {
if (js->eh_ipt[34] > 0) {
list[10] = 2;
return(state(et, list, FALSE, pv, pvsun, rrd, serr));
} else {
if (serr != NULL)
sprintf(serr,"No nutations on the JPL ephemeris file;");
return (NOT_AVAILABLE);
}
}
if (ntarg == J_LIB) {
if (js->eh_ipt[37] > 0) {
list[11] = 2;
if ((retc = state(et, list, FALSE, pv, pvsun, rrd, serr)) != OK)
return (retc);
for (i = 0; i < 6; ++i)
rrd[i] = pv[i + 60];
return 0;
} else {
if (serr != NULL)
sprintf(serr,"No librations on the ephemeris file;");
return (NOT_AVAILABLE);
}
}
/* set up proper entries in 'list' array for state call */
if (ntarg < J_SUN)
list[ntarg] = 2;
if (ntarg == J_MOON) /* Mooon needs Earth */
list[J_EARTH] = 2;
if (ntarg == J_EARTH) /* Earth needs Moon */
list[J_MOON] = 2;
if (ntarg == J_EMB) /* EMB needs Earth */
list[J_EARTH] = 2;
if (ncent < J_SUN)
list[ncent] = 2;
if (ncent == J_MOON) /* Mooon needs Earth */
list[J_EARTH] = 2;
if (ncent == J_EARTH) /* Earth needs Moon */
list[J_MOON] = 2;
if (ncent == J_EMB) /* EMB needs Earth */
list[J_EARTH] = 2;
if ((retc = state(et, list, TRUE, pv, pvsun, rrd, serr)) != OK)
return (retc);
if (ntarg == J_SUN || ncent == J_SUN) {
for (i = 0; i < 6; ++i)
pv[i + 6*J_SUN] = pvsun[i];
}
if (ntarg == J_SBARY || ncent == J_SBARY) {
for (i = 0; i < 6; ++i) {
pv[i + 6*J_SBARY] = 0.;
}
}
if (ntarg == J_EMB || ncent == J_EMB) {
for (i = 0; i < 6; ++i)
pv[i + 6*J_EMB] = pv[i + 6*J_EARTH];
}
if ((ntarg==J_EARTH && ncent==J_MOON) || (ntarg == J_MOON && ncent==J_EARTH)){
for (i = 0; i < 6; ++i)
pv[i + 6*J_EARTH] = 0.;
} else {
if (list[J_EARTH] == 2) {
for (i = 0; i < 6; ++i)
pv[i + 6*J_EARTH] -= pv[i + 6*J_MOON] / (js->eh_emrat + 1.);
}
if (list[J_MOON] == 2) {
for (i = 0; i < 6; ++i) {
pv[i + 6*J_MOON] += pv[i + 6*J_EARTH];
}
}
}
for (i = 0; i < 6; ++i)
rrd[i] = pv[i + ntarg * 6] - pv[i + ncent * 6];
return OK;
}
/*
* This subroutine differentiates and interpolates a
* set of chebyshev coefficients to give pos, vel, acc, and jerk
* calling sequence parameters:
* input:
* buf 1st location of array of d.p. chebyshev coefficients of position
* t is dp fractional time in interval covered by
* coefficients at which interpolation is wanted, 0 <= t <= 1
* intv is dp length of whole interval in input time units.
* ncf number of coefficients per component
* ncm number of components per set of coefficients
* na number of sets of coefficients in full array
* (i.e., number of sub-intervals in full interval)
* ifl int flag: =1 for positions only
* =2 for pos and vel
* =3 for pos, vel, and acc
* =4 for pos, vel, acc, and jerk
* output:
* pv d.p. interpolated quantities requested.
* assumed dimension is pv(ncm,fl).
*/
static int interp(double *buf, double t, double intv, int32 ncfin,
int32 ncmin, int32 nain, int32 ifl, double *pv)
{
/* Initialized data */
static TLS int np, nv;
static TLS int nac;
static TLS int njk;
static TLS double twot = 0.;
double *pc = js->pc;
double *vc = js->vc;
double *ac = js->ac;
double *jc = js->jc;
int ncf = (int) ncfin;
int ncm = (int) ncmin;
int na = (int) nain;
/* Local variables */
double temp;
int i, j, ni;
double tc;
double dt1, bma;
double bma2, bma3;
/*
| get correct sub-interval number for this set of coefficients and then
| get normalized chebyshev time within that subinterval.
*/
if (t >= 0)
dt1 = floor(t);
else
dt1 = -floor(-t);
temp = na * t;
ni = (int) (temp - dt1);
/* tc is the normalized chebyshev time (-1 <= tc <= 1) */
tc = (fmod(temp, 1.0) + dt1) * 2. - 1.;
/*
* check to see whether chebyshev time has changed,
* and compute new polynomial values if it has.
* (the element pc(2) is the value of t1(tc) and hence
* contains the value of tc on the previous call.)
*/
if (tc != pc[1]) {
np = 2;
nv = 3;
nac = 4;
njk = 5;
pc[1] = tc;
twot = tc + tc;
}
/*
* be sure that at least 'ncf' polynomials have been evaluated
* and are stored in the array 'pc'.
*/
if (np < ncf) {
for (i = np; i < ncf; ++i)
pc[i] = twot * pc[i - 1] - pc[i - 2];
np = ncf;
}
/* interpolate to get position for each component */
for (i = 0; i < ncm; ++i) {
pv[i] = 0.;
for (j = ncf-1; j >= 0; --j)
pv[i] += pc[j] * buf[j + (i + ni * ncm) * ncf];
}
if (ifl <= 1)
return 0;
/*
* if velocity interpolation is wanted, be sure enough
* derivative polynomials have been generated and stored.
*/
bma = (na + na) / intv;
vc[2] = twot + twot;
if (nv < ncf) {
for (i = nv; i < ncf; ++i)
vc[i] = twot * vc[i - 1] + pc[i - 1] + pc[i - 1] - vc[i - 2];
nv = ncf;
}
/* interpolate to get velocity for each component */
for (i = 0; i < ncm; ++i) {
pv[i + ncm] = 0.;
for (j = ncf-1; j >= 1; --j)
pv[i + ncm] += vc[j] * buf[j + (i + ni * ncm) * ncf];
pv[i + ncm] *= bma;
}
if (ifl == 2)
return 0;
/* check acceleration polynomial values, and */
/* re-do if necessary */
bma2 = bma * bma;
ac[3] = pc[1] * 24.;
if (nac < ncf) {
nac = ncf;
for (i = nac; i < ncf; ++i)
ac[i] = twot * ac[i - 1] + vc[i - 1] * 4. - ac[i - 2];
}
/* get acceleration for each component */
for (i = 0; i < ncm; ++i) {
pv[i + ncm * 2] = 0.;
for (j = ncf-1; j >= 2; --j)
pv[i + ncm * 2] += ac[j] * buf[j + (i + ni * ncm) * ncf];
pv[i + ncm * 2] *= bma2;
}
if (ifl == 3)
return 0;
/* check jerk polynomial values, and */
/* re-do if necessary */
bma3 = bma * bma2;
jc[4] = pc[1] * 192.;
if (njk < ncf) {
njk = ncf;
for (i = njk; i < ncf; ++i)
jc[i] = twot * jc[i - 1] + ac[i - 1] * 6. - jc[i - 2];
}
/* get jerk for each component */
for (i = 0; i < ncm; ++i) {
pv[i + ncm * 3] = 0.;
for (j = ncf-1; j >= 3; --j)
pv[i + ncm * 3] += jc[j] * buf[j + (i + ni * ncm) * ncf];
pv[i + ncm * 3] *= bma3;
}
return 0;
}
/*
| ********** state ********************
| this subroutine reads and interpolates the jpl planetary ephemeris file
| calling sequence parameters:
| input:
| et dp julian ephemeris epoch at which interpolation is wanted.
| list 12-word integer array specifying what interpolation
| is wanted for each of the bodies on the file.
| list(i)=0, no interpolation for body i
| =1, position only
| =2, position and velocity
| the designation of the astronomical bodies by i is:
| i = 0: mercury
| = 1: venus
| = 2: earth-moon barycenter, NOT earth!
| = 3: mars
| = 4: jupiter
| = 5: saturn
| = 6: uranus
| = 7: neptune
| = 8: pluto
| = 9: geocentric moon
| =10: nutations in longitude and obliquity
| =11: lunar librations (if on file)
| If called with list = NULL, only the header records are read and
| stored in the global areas.
| do_bary short, if true, barycentric, if false, heliocentric.
| only the 9 planets 0..8 are affected by it.
| output:
| pv dp 6 x 11 array that will contain requested interpolated
| quantities. the body specified by list(i) will have its
| state in the array starting at pv(1,i). (on any given
| call, only those words in 'pv' which are affected by the
| first 10 'list' entries (and by list(11) if librations are
| on the file) are set. the rest of the 'pv' array
| is untouched.) the order of components starting in
| pv is: x,y,z,dx,dy,dz.
| all output vectors are referenced to the earth mean
| equator and equinox of epoch. the moon state is always
| geocentric; the other nine states are either heliocentric
| or solar-system barycentric, depending on the setting of
| common flags (see below).
| lunar librations, if on file, are put into pv(k,10) if
| list(11) is 1 or 2.
| pvsun dp 6-word array containing the barycentric position and
| velocity of the sun.
| nut dp 4-word array that will contain nutations and rates,
| depending on the setting of list(10). the order of
| quantities in nut is:
| d psi (nutation in longitude)
| d epsilon (nutation in obliquity)
| d psi dot
| d epsilon dot
| globals used:
| do_km logical flag defining physical units of the output states.
| TRUE = return km and km/sec, FALSE = return au and au/day
| default value = FALSE (km determines time unit
| for nutations and librations. angle unit is always radians.)
*/
static int state(double et, int32 *list, int do_bary,
double *pv, double *pvsun, double *nut, char *serr)
{
int i, j, k;
int32 nseg;
off_t64 flen, nb;
double *buf = js->buf;
double aufac, s, t, intv, ts[4];
int32 nrecl, ksize;
int32 nr;
double et_mn, et_fr;
int32 *ipt = js->eh_ipt;
char ch_ttl[252];
static TLS int32 irecsz;
static TLS int32 nrl, lpt[3], ncoeffs;
size_t nrd; /* unused, removes compile warnings */
if (js->jplfptr == NULL) {
ksize = fsizer(serr); /* the number of single precision words in a record */
nrecl = 4;
if (ksize == NOT_AVAILABLE)
return NOT_AVAILABLE;
irecsz = nrecl * ksize; /* record size in bytes */
ncoeffs = ksize / 2; /* # of coefficients, doubles */
/* ttl = ephemeris title, e.g.
* "JPL Planetary Ephemeris DE404/LE404
* Start Epoch: JED= 625296.5-3001 DEC 21 00:00:00
* Final Epoch: JED= 2817168.5 3001 JAN 17 00:00:00c */
nrd = fread((void *) ch_ttl, 1, 252, js->jplfptr);
if (nrd != 252) return NOT_AVAILABLE;
/* cnam = names of constants */
nrd = fread((void *) js->ch_cnam, 1, 2400, js->jplfptr);
if (nrd != 2400) return NOT_AVAILABLE;
/* ss[0] = start epoch of ephemeris
* ss[1] = end epoch
* ss[2] = segment size in days */
nrd = fread((void *) &js->eh_ss[0], sizeof(double), 3, js->jplfptr);
if (nrd != 3) return NOT_AVAILABLE;
if (js->do_reorder)
reorder((char *) &js->eh_ss[0], sizeof(double), 3);
/* ncon = number of constants */
nrd = fread((void *) &js->eh_ncon, sizeof(int32), 1, js->jplfptr);
if (nrd != 1) return NOT_AVAILABLE;
if (js->do_reorder)
reorder((char *) &js->eh_ncon, sizeof(int32), 1);
/* au = astronomical unit */
nrd = fread((void *) &js->eh_au, sizeof(double), 1, js->jplfptr);
if (nrd != 1) return NOT_AVAILABLE;
if (js->do_reorder)
reorder((char *) &js->eh_au, sizeof(double), 1);
/* emrat = earth moon mass ratio */
nrd = fread((void *) &js->eh_emrat, sizeof(double), 1, js->jplfptr);
if (nrd != 1) return NOT_AVAILABLE;
if (js->do_reorder)
reorder((char *) &js->eh_emrat, sizeof(double), 1);
/* ipt[i+0]: coefficients of planet i start at buf[ipt[i+0]-1]
* ipt[i+1]: number of coefficients (interpolation order - 1)
* ipt[i+2]: number of intervals in segment */
nrd = fread((void *) &ipt[0], sizeof(int32), 36, js->jplfptr);
if (nrd != 36) return NOT_AVAILABLE;
if (js->do_reorder)
reorder((char *) &ipt[0], sizeof(int32), 36);
/* numde = number of jpl ephemeris "404" with de404 */
nrd = fread((void *) &js->eh_denum, sizeof(int32), 1, js->jplfptr);
if (nrd != 1) return NOT_AVAILABLE;
if (js->do_reorder)
reorder((char *) &js->eh_denum, sizeof(int32), 1);
nrd = fread((void *) &lpt[0], sizeof(int32), 3, js->jplfptr);
if (nrd != 3) return NOT_AVAILABLE;
if (js->do_reorder)
reorder((char *) &lpt[0], sizeof(int32), 3);
/* cval[]: other constants in next record */
FSEEK(js->jplfptr, (off_t64) (1L * irecsz), 0);
nrd = fread((void *) &js->eh_cval[0], sizeof(double), 400, js->jplfptr);
if (nrd != 400) return NOT_AVAILABLE;
if (js->do_reorder)
reorder((char *) &js->eh_cval[0], sizeof(double), 400);
/* new 26-aug-2008: verify correct block size */
for (i = 0; i < 3; ++i)
ipt[i + 36] = lpt[i];
nrl = 0;
/* is file length correct? */
/* file length */
FSEEK(js->jplfptr, (off_t64) 0L, SEEK_END);
flen = FTELL(js->jplfptr);
/* # of segments in file */
nseg = (int32) ((js->eh_ss[1] - js->eh_ss[0]) / js->eh_ss[2]);
/* sum of all cheby coeffs of all planets and segments */
for(i = 0, nb = 0; i < 13; i++) {
k = 3;
if (i == 11)
k = 2;
nb += (ipt[i*3+1] * ipt[i*3+2]) * k * nseg;
}
/* add start and end epochs of segments */
nb += 2 * nseg;
/* doubles to bytes */
nb *= 8;
/* add size of header and constants section */
nb += 2 * ksize * nrecl;
if (flen != nb
/* some of our files are one record too long */
&& flen - nb != ksize * nrecl
) {
if (serr != NULL) {
sprintf(serr, "JPL ephemeris file is mutilated; length = %d instead of %d.", (unsigned int) flen, (unsigned int) nb);
if (strlen(serr) + strlen(js->jplfname) < AS_MAXCH - 1) {
sprintf(serr, "JPL ephemeris file %s is mutilated; length = %d instead of %d.", js->jplfname, (unsigned int) flen, (unsigned int) nb);
}
}
return(NOT_AVAILABLE);
}
/* check if start and end dates in segments are the same as in
* file header */
FSEEK(js->jplfptr, (off_t64) (2L * irecsz), 0);
nrd = fread((void *) &ts[0], sizeof(double), 2, js->jplfptr);
if (nrd != 2) return NOT_AVAILABLE;
if (js->do_reorder)
reorder((char *) &ts[0], sizeof(double), 2);
FSEEK(js->jplfptr, (off_t64) ((nseg + 2 - 1) * ((off_t64) irecsz)), 0);
nrd = fread((void *) &ts[2], sizeof(double), 2, js->jplfptr);
if (nrd != 2) return NOT_AVAILABLE;
if (js->do_reorder)
reorder((char *) &ts[2], sizeof(double), 2);
if (ts[0] != js->eh_ss[0] || ts[3] != js->eh_ss[1]) {
if (serr != NULL)
sprintf(serr, "JPL ephemeris file is corrupt; start/end date check failed. %.1f != %.1f || %.1f != %.1f", ts[0],js->eh_ss[0],ts[3],js->eh_ss[1]);
return NOT_AVAILABLE;
}
}
if (list == NULL)
return 0;
s = et - .5;
et_mn = floor(s);
et_fr = s - et_mn; /* fraction of days since previous midnight */
et_mn += .5; /* midnight before epoch */
/* error return for epoch out of range */
if (et < js->eh_ss[0] || et > js->eh_ss[1]) {
if (serr != NULL)
sprintf(serr,"jd %f outside JPL eph. range %.2f .. %.2f;", et, js->eh_ss[0], js->eh_ss[1]);
return BEYOND_EPH_LIMITS;
}
/* calculate record # and relative time in interval */
nr = (int32) ((et_mn - js->eh_ss[0]) / js->eh_ss[2]) + 2;
if (et_mn == js->eh_ss[1])
--nr; /* end point of ephemeris, use last record */
t = (et_mn - ((nr - 2) * js->eh_ss[2] + js->eh_ss[0]) + et_fr) / js->eh_ss[2];
/* read correct record if not in core */
if (nr != nrl) {
nrl = nr;
if (FSEEK(js->jplfptr, (off_t64) (nr * ((off_t64) irecsz)), 0) != 0) {
if (serr != NULL)
sprintf(serr, "Read error in JPL eph. at %f\n", et);
return NOT_AVAILABLE;
}
for (k = 1; k <= ncoeffs; ++k) {
if ( fread((void *) &buf[k - 1], sizeof(double), 1, js->jplfptr) != 1) {
if (serr != NULL)
sprintf(serr, "Read error in JPL eph. at %f\n", et);
return NOT_AVAILABLE;
}
if (js->do_reorder)
reorder((char *) &buf[k-1], sizeof(double), 1);
}
}
if (js->do_km) {
intv = js->eh_ss[2] * 86400.;
aufac = 1.;
} else {
intv = js->eh_ss[2];
aufac = 1. / js->eh_au;
}
/* interpolate ssbary sun */
interp(&buf[(int) ipt[30] - 1], t, intv, ipt[31], 3L, ipt[32], 2L, pvsun);
for (i = 0; i < 6; ++i) {
pvsun[i] *= aufac;
}
/* check and interpolate whichever bodies are requested */
for (i = 0; i < 10; ++i) {
if (list[i] > 0) {
interp(&buf[(int) ipt[i * 3] - 1], t, intv, ipt[i * 3 + 1], 3L,
ipt[i * 3 + 2], list[i], &pv[i * 6]);
for (j = 0; j < 6; ++j) {
if (i < 9 && ! do_bary) {
pv[j + i * 6] = pv[j + i * 6] * aufac - pvsun[j];
} else {
pv[j + i * 6] *= aufac;
}
}
}
}
/* do nutations if requested (and if on file) */
if (list[10] > 0 && ipt[34] > 0) {
interp(&buf[(int) ipt[33] - 1], t, intv, ipt[34], 2L, ipt[35],
list[10], nut);
}
/* get librations if requested (and if on file) */
if (list[11] > 0 && ipt[37] > 0) {
interp(&buf[(int) ipt[36] - 1], t, intv, ipt[37], 3L, ipt[38], list[1],
&pv[60]);
}
return OK;
}
/*
* this entry obtains the constants from the ephemeris file
* call state to initialize the ephemeris and read in the constants
*/
static int read_const_jpl(double *ss, char *serr)
{
int i, retc;
retc = state(0.0, NULL, FALSE, NULL, NULL, NULL, serr);
if (retc != OK)
return (retc);
for (i = 0; i < 3; i++)
ss[i] = js->eh_ss[i];
#if DEBUG_DO_SHOW
{
static const char *bname[] = {
"Mercury", "Venus", "EMB", "Mars", "Jupiter", "Saturn",
"Uranus", "Neptune", "Pluto", "Moon", "SunBary", "Nut", "Libr"};
int j, k;
int32 nb, nc;
printf(" JPL TEST-EPHEMERIS program. Version October 1995.\n");
for (i = 0; i < 13; i++) {
j = i * 3;
k = 3;
if (i == 11) k = 2;
nb = js->eh_ipt[j+1] * js->eh_ipt[j+2] * k;
nc = (int32) (nb * 36525L / js->eh_ss[2] * 8L);
printf("%s\t%d\tipt[%d]\t%3ld %2ld %2ld,\t",
bname[i], i, j, js->eh_ipt[j], js->eh_ipt[j+1], js->eh_ipt[j+2]);
printf("%3ld double, bytes per century = %6ld\n", nb, nc);
fflush(stdout);
}
printf("%16.2f %16.2f %16.2f\n", js->eh_ss[0], js->eh_ss[1], js->eh_ss[2]);
for (i = 0; i < js->eh_ncon; ++i)
printf("%.6s\t%24.16f\n", js->ch_cnam + i * 6, js->eh_cval[i]);
fflush(stdout);
}
#endif
return OK;
}
static void reorder(char *x, int size, int number)
{
int i, j;
char s[8];
char *sp1 = x;
char *sp2 = &s[0];
for (i = 0; i < number; i++) {
for (j = 0; j < size; j++)
*(sp2 + j) = *(sp1 + size - j - 1);
for (j = 0; j < size; j++)
*(sp1 + j) = *(sp2 + j);
sp1 += size;
}
}
void swi_close_jpl_file(void)
{
if (js != NULL) {
if (js->jplfptr != NULL)
fclose(js->jplfptr);
if (js->jplfname != NULL)
FREE((void *) js->jplfname);
if (js->jplfpath != NULL)
FREE((void *) js->jplfpath);
FREE((void *) js);
js = NULL;
}
}
int swi_open_jpl_file(double *ss, char *fname, char *fpath, char *serr)
{
int retc = OK;
/* if open, return */
if (js != NULL && js->jplfptr != NULL)
return OK;
if ((js = (struct jpl_save *) CALLOC(1, sizeof(struct jpl_save))) == NULL
|| (js->jplfname = (char *) MALLOC(strlen(fname)+1)) == NULL
|| (js->jplfpath = (char *) MALLOC(strlen(fpath)+1)) == NULL
) {
if (serr != NULL)
strcpy(serr, "error in malloc() with JPL ephemeris.");
return ERR;
}
strcpy(js->jplfname, fname);
strcpy(js->jplfpath, fpath);
retc = read_const_jpl(ss, serr);
if (retc != OK)
swi_close_jpl_file();
else {
/* intializations for function interpol() */
js->pc[0] = 1;
js->pc[1] = 2;
js->vc[1] = 1;
js->ac[2] = 4;
js->jc[3] = 24;
}
return retc;
}
int32 swi_get_jpl_denum()
{
return js->eh_denum;
}

103
src/swejpl.h Normal file
View File

@ -0,0 +1,103 @@
/*
|
| Subroutines for reading JPL ephemerides.
| derived from testeph.f as contained in DE403 distribution July 1995.
| works with DE200, DE102, DE403, DE404, DE405, DE406, DE431
| (attention, these ephemerides do not have exactly the same reference frame)
Authors: Dieter Koch and Alois Treindl, Astrodienst Zurich
**************************************************************/
/* 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 "sweodef.h"
#define J_MERCURY 0 /* jpl body indices, modified by Alois */
#define J_VENUS 1 /* now they start at 0 and not at 1 */
#define J_EARTH 2
#define J_MARS 3
#define J_JUPITER 4
#define J_SATURN 5
#define J_URANUS 6
#define J_NEPTUNE 7
#define J_PLUTO 8
#define J_MOON 9
#define J_SUN 10
#define J_SBARY 11
#define J_EMB 12
#define J_NUT 13
#define J_LIB 14
/*
* compute position and speed at time et, for body ntarg with center
* ncent. rrd must be double[6] to contain the return vectors.
* ntarg can be all of the above, ncent all except J_NUT and J_LIB.
* Librations and Nutations are not affected by ncent.
*/
extern int swi_pleph(double et, int ntarg, int ncent, double *rrd, char *serr);
/*
* read the ephemeris constants. ss[0..2] returns start, end and granule size.
* If do_show is TRUE, a list of constants is printed to stdout.
*/
extern void swi_close_jpl_file(void);
extern int swi_open_jpl_file(double *ss, char *fname, char *fpath, char *serr);
extern int32 swi_get_jpl_denum(void);
extern void swi_IERS_FK5(double *xin, double *xout, int dir);

131
src/swemini.c Normal file
View File

@ -0,0 +1,131 @@
/*
swemini.c A minimal program to test the Swiss Ephemeris.
Input: a date (in gregorian calendar, sequence day.month.year)
Output: Planet positions at midnight Universal time, ecliptic coordinates,
geocentric apparent positions relative to true equinox of date, as
usual in western astrology.
Authors: Dieter Koch and Alois Treindl, Astrodienst Zurich
**************************************************************/
/* 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 "swephexp.h" /* this includes "sweodef.h" */
int main()
{
char sdate[AS_MAXCH], snam[40], serr[AS_MAXCH];
int jday = 1, jmon = 1, jyear = 2000;
double jut = 0.0;
double tjd, te, x2[6];
int32 iflag, iflgret;
int p;
swe_set_ephe_path(NULL);
iflag = SEFLG_SPEED;
while (TRUE) {
printf("\nDate (d.m.y) ?");
/*gets(sdate);*/
if( !fgets(sdate, sizeof(sdate)-1, stdin) ) return OK;
/*
* stop if a period . is entered
*/
if (*sdate == '.')
return OK;
if (sscanf (sdate, "%d%*c%d%*c%d", &jday,&jmon,&jyear) < 1) exit(1);
/*
* we have day, month and year and convert to Julian day number
*/
tjd = swe_julday(jyear,jmon,jday,jut,SE_GREG_CAL);
/*
* compute Ephemeris time from Universal time by adding delta_t
*/
te = tjd + swe_deltat(tjd);
printf("date: %02d.%02d.%d at 0:00 Universal time\n", jday, jmon, jyear);
printf("planet \tlongitude\tlatitude\tdistance\tspeed long.\n");
/*
* a loop over all planets
*/
for (p = SE_SUN; p <= SE_CHIRON; p++) {
if (p == SE_EARTH) continue;
/*
* do the coordinate calculation for this planet p
*/
iflgret = swe_calc(te, p, iflag, x2, serr);
/*
* if there is a problem, a negative value is returned and an
* errpr message is in serr.
*/
if (iflgret < 0)
printf("error: %s\n", serr);
else if (iflgret != iflag)
printf("warning: iflgret != iflag. %s\n", serr);
/*
* get the name of the planet p
*/
swe_get_planet_name(p, snam);
/*
* print the coordinates
*/
printf("%10s\t%11.7f\t%10.7f\t%10.7f\t%10.7f\n",
snam, x2[0], x2[1], x2[2], x2[3]);
}
}
return OK;
}

1930
src/swemmoon.c Normal file

File diff suppressed because it is too large Load Diff

967
src/swemplan.c Normal file
View File

@ -0,0 +1,967 @@
/* SWISSEPH
Moshier planet routines
modified for SWISSEPH by 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 <string.h>
#include "swephexp.h"
#include "sweph.h"
#include "swephlib.h"
#include "swemptab.h"
#define TIMESCALE 3652500.0
#define mods3600(x) ((x) - 1.296e6 * floor ((x)/1.296e6))
#define FICT_GEO 1
#define KGAUSS_GEO 0.0000298122353216 /* Earth only */
/* #define KGAUSS_GEO 0.00002999502129737 Earth + Moon */
static void embofs_mosh(double J, double *xemb);
static int check_t_terms(double t, char *sinp, double *doutp);
static int read_elements_file(int32 ipl, double tjd,
double *tjd0, double *tequ,
double *mano, double *sema, double *ecce,
double *parg, double *node, double *incl,
char *pname, int32 *fict_ifl, char *serr);
static const int pnoint2msh[] = {2, 2, 0, 1, 3, 4, 5, 6, 7, 8, };
/* From Simon et al (1994) */
static const double freqs[] =
{
/* Arc sec per 10000 Julian years. */
53810162868.8982,
21066413643.3548,
12959774228.3429,
6890507749.3988,
1092566037.7991,
439960985.5372,
154248119.3933,
78655032.0744,
52272245.1795
};
static const double phases[] =
{
/* Arc sec. */
252.25090552 * 3600.,
181.97980085 * 3600.,
100.46645683 * 3600.,
355.43299958 * 3600.,
34.35151874 * 3600.,
50.07744430 * 3600.,
314.05500511 * 3600.,
304.34866548 * 3600.,
860492.1546,
};
static const struct plantbl *planets[] =
{
&mer404,
&ven404,
&ear404,
&mar404,
&jup404,
&sat404,
&ura404,
&nep404,
&plu404
};
static TLS double ss[9][24];
static TLS double cc[9][24];
static void sscc (int k, double arg, int n);
int swi_moshplan2 (double J, int iplm, double *pobj)
{
int i, j, k, m, k1, ip, np, nt;
signed char *p;
double *pl, *pb, *pr;
double su, cu, sv, cv, T;
double t, sl, sb, sr;
const struct plantbl *plan = planets[iplm];
T = (J - J2000) / TIMESCALE;
/* Calculate sin( i*MM ), etc. for needed multiple angles. */
for (i = 0; i < 9; i++)
{
if ((j = plan->max_harmonic[i]) > 0)
{
sr = (mods3600 (freqs[i] * T) + phases[i]) * STR;
sscc (i, sr, j);
}
}
/* Point to start of table of arguments. */
p = plan->arg_tbl;
/* Point to tabulated cosine and sine amplitudes. */
pl = plan->lon_tbl;
pb = plan->lat_tbl;
pr = plan->rad_tbl;
sl = 0.0;
sb = 0.0;
sr = 0.0;
for (;;)
{
/* argument of sine and cosine */
/* Number of periodic arguments. */
np = *p++;
if (np < 0)
break;
if (np == 0)
{ /* It is a polynomial term. */
nt = *p++;
/* Longitude polynomial. */
cu = *pl++;
for (ip = 0; ip < nt; ip++)
{
cu = cu * T + *pl++;
}
sl += mods3600 (cu);
/* Latitude polynomial. */
cu = *pb++;
for (ip = 0; ip < nt; ip++)
{
cu = cu * T + *pb++;
}
sb += cu;
/* Radius polynomial. */
cu = *pr++;
for (ip = 0; ip < nt; ip++)
{
cu = cu * T + *pr++;
}
sr += cu;
continue;
}
k1 = 0;
cv = 0.0;
sv = 0.0;
for (ip = 0; ip < np; ip++)
{
/* What harmonic. */
j = *p++;
/* Which planet. */
m = *p++ - 1;
if (j)
{
k = j;
if (j < 0)
k = -k;
k -= 1;
su = ss[m][k]; /* sin(k*angle) */
if (j < 0)
su = -su;
cu = cc[m][k];
if (k1 == 0)
{ /* set first angle */
sv = su;
cv = cu;
k1 = 1;
}
else
{ /* combine angles */
t = su * cv + cu * sv;
cv = cu * cv - su * sv;
sv = t;
}
}
}
/* Highest power of T. */
nt = *p++;
/* Longitude. */
cu = *pl++;
su = *pl++;
for (ip = 0; ip < nt; ip++)
{
cu = cu * T + *pl++;
su = su * T + *pl++;
}
sl += cu * cv + su * sv;
/* Latitiude. */
cu = *pb++;
su = *pb++;
for (ip = 0; ip < nt; ip++)
{
cu = cu * T + *pb++;
su = su * T + *pb++;
}
sb += cu * cv + su * sv;
/* Radius. */
cu = *pr++;
su = *pr++;
for (ip = 0; ip < nt; ip++)
{
cu = cu * T + *pr++;
su = su * T + *pr++;
}
sr += cu * cv + su * sv;
}
pobj[0] = STR * sl;
pobj[1] = STR * sb;
pobj[2] = STR * plan->distance * sr + plan->distance;
return OK;
}
/* Moshier ephemeris.
* computes heliocentric cartesian equatorial coordinates of
* equinox 2000
* for earth and a planet
* tjd julian day
* ipli internal SWEPH planet number
* xp array of 6 doubles for planet's position and speed
* xe earth's
* serr error string
*/
int swi_moshplan(double tjd, int ipli, AS_BOOL do_save, double *xpret, double *xeret, char *serr)
{
int i;
int do_earth = FALSE;
double dx[3], x2[3], xxe[6], xxp[6];
double *xp, *xe;
double dt;
char s[AS_MAXCH];
int iplm = pnoint2msh[ipli];
struct plan_data *pdp = &swed.pldat[ipli];
struct plan_data *pedp = &swed.pldat[SEI_EARTH];
double seps2000 = swed.oec2000.seps;
double ceps2000 = swed.oec2000.ceps;
if (do_save) {
xp = pdp->x;
xe = pedp->x;
} else {
xp = xxp;
xe = xxe;
}
if (do_save || ipli == SEI_EARTH || xeret != NULL)
do_earth = TRUE;
/* tjd beyond ephemeris limits, give some margin for spped at edge */
if (tjd < MOSHPLEPH_START - 0.3 || tjd > MOSHPLEPH_END + 0.3) {
if (serr != NULL) {
sprintf(s, "jd %f outside Moshier planet range %.2f .. %.2f ",
tjd, MOSHPLEPH_START, MOSHPLEPH_END);
if (strlen(serr) + strlen(s) < AS_MAXCH)
strcat(serr, s);
}
return(ERR);
}
/* earth, for geocentric position */
if (do_earth) {
if (tjd == pedp->teval
&& pedp->iephe == SEFLG_MOSEPH) {
xe = pedp->x;
} else {
/* emb */
swi_moshplan2(tjd, pnoint2msh[SEI_EMB], xe); /* emb hel. ecl. 2000 polar */
swi_polcart(xe, xe); /* to cartesian */
swi_coortrf2(xe, xe, -seps2000, ceps2000);/* and equator 2000 */
embofs_mosh(tjd, xe); /* emb -> earth */
if (do_save) {
pedp->teval = tjd;
pedp->xflgs = -1;
pedp->iephe = SEFLG_MOSEPH;
}
/* one more position for speed. */
swi_moshplan2(tjd - PLAN_SPEED_INTV, pnoint2msh[SEI_EMB], x2);
swi_polcart(x2, x2);
swi_coortrf2(x2, x2, -seps2000, ceps2000);
embofs_mosh(tjd - PLAN_SPEED_INTV, x2);/**/
for (i = 0; i <= 2; i++)
dx[i] = (xe[i] - x2[i]) / PLAN_SPEED_INTV;
/* store speed */
for (i = 0; i <= 2; i++) {
xe[i+3] = dx[i];
}
}
if (xeret != NULL)
for (i = 0; i <= 5; i++)
xeret[i] = xe[i];
}
/* earth is the planet wanted */
if (ipli == SEI_EARTH) {
xp = xe;
} else {
/* other planet */
/* if planet has already been computed, return */
if (tjd == pdp->teval && pdp->iephe == SEFLG_MOSEPH) {
xp = pdp->x;
} else {
swi_moshplan2(tjd, iplm, xp);
swi_polcart(xp, xp);
swi_coortrf2(xp, xp, -seps2000, ceps2000);
if (do_save) {
pdp->teval = tjd;/**/
pdp->xflgs = -1;
pdp->iephe = SEFLG_MOSEPH;
}
/* one more position for speed.
* the following dt gives good speed for light-time correction
*/
#if 0
for (i = 0; i <= 2; i++)
dx[i] = xp[i] - pedp->x[i];
dt = LIGHTTIME_AUNIT * sqrt(square_sum(dx));
#endif
dt = PLAN_SPEED_INTV;
swi_moshplan2(tjd - dt, iplm, x2);
swi_polcart(x2, x2);
swi_coortrf2(x2, x2, -seps2000, ceps2000);
for (i = 0; i <= 2; i++)
dx[i] = (xp[i] - x2[i]) / dt;
/* store speed */
for (i = 0; i <= 2; i++) {
xp[i+3] = dx[i];
}
}
if (xpret != NULL)
for (i = 0; i <= 5; i++)
xpret[i] = xp[i];
}
return(OK);
}
/* Prepare lookup table of sin and cos ( i*Lj )
* for required multiple angles
*/
static void sscc (int k, double arg, int n)
{
double cu, su, cv, sv, s;
int i;
su = sin (arg);
cu = cos (arg);
ss[k][0] = su; /* sin(L) */
cc[k][0] = cu; /* cos(L) */
sv = 2.0 * su * cu;
cv = cu * cu - su * su;
ss[k][1] = sv; /* sin(2L) */
cc[k][1] = cv;
for (i = 2; i < n; i++)
{
s = su * cv + cu * sv;
cv = cu * cv - su * sv;
sv = s;
ss[k][i] = sv; /* sin( i+1 L ) */
cc[k][i] = cv;
}
}
/* Adjust position from Earth-Moon barycenter to Earth
*
* J = Julian day number
* xemb = rectangular equatorial coordinates of Earth
*/
static void embofs_mosh(double tjd, double *xemb)
{
double T, M, a, L, B, p;
double smp, cmp, s2mp, c2mp, s2d, c2d, sf, cf;
double s2f, sx, cx, xyz[6];
double seps = swed.oec.seps;
double ceps = swed.oec.ceps;
int i;
/* Short series for position of the Moon
*/
T = (tjd-J1900)/36525.0;
/* Mean anomaly of moon (MP) */
a = swe_degnorm(((1.44e-5*T + 0.009192)*T + 477198.8491)*T + 296.104608);
a *= DEGTORAD;
smp = sin(a);
cmp = cos(a);
s2mp = 2.0*smp*cmp; /* sin(2MP) */
c2mp = cmp*cmp - smp*smp; /* cos(2MP) */
/* Mean elongation of moon (D) */
a = swe_degnorm(((1.9e-6*T - 0.001436)*T + 445267.1142)*T + 350.737486);
a = 2.0 * DEGTORAD * a;
s2d = sin(a);
c2d = cos(a);
/* Mean distance of moon from its ascending node (F) */
a = swe_degnorm((( -3.e-7*T - 0.003211)*T + 483202.0251)*T + 11.250889);
a *= DEGTORAD;
sf = sin(a);
cf = cos(a);
s2f = 2.0*sf*cf; /* sin(2F) */
sx = s2d*cmp - c2d*smp; /* sin(2D - MP) */
cx = c2d*cmp + s2d*smp; /* cos(2D - MP) */
/* Mean longitude of moon (LP) */
L = ((1.9e-6*T - 0.001133)*T + 481267.8831)*T + 270.434164;
/* Mean anomaly of sun (M) */
M = swe_degnorm((( -3.3e-6*T - 1.50e-4)*T + 35999.0498)*T + 358.475833);
/* Ecliptic longitude of the moon */
L = L
+ 6.288750*smp
+ 1.274018*sx
+ 0.658309*s2d
+ 0.213616*s2mp
- 0.185596*sin( DEGTORAD * M )
- 0.114336*s2f;
/* Ecliptic latitude of the moon */
a = smp*cf;
sx = cmp*sf;
B = 5.128189*sf
+ 0.280606*(a+sx) /* sin(MP+F) */
+ 0.277693*(a-sx) /* sin(MP-F) */
+ 0.173238*(s2d*cf - c2d*sf); /* sin(2D-F) */
B *= DEGTORAD;
/* Parallax of the moon */
p = 0.950724
+0.051818*cmp
+0.009531*cx
+0.007843*c2d
+0.002824*c2mp;
p *= DEGTORAD;
/* Elongation of Moon from Sun
*/
L = swe_degnorm(L);
L *= DEGTORAD;
/* Distance in au */
a = 4.263523e-5/sin(p);
/* Convert to rectangular ecliptic coordinates */
xyz[0] = L;
xyz[1] = B;
xyz[2] = a;
swi_polcart(xyz, xyz);
/* Convert to equatorial */
swi_coortrf2(xyz, xyz, -seps, ceps);
/* Precess to equinox of J2000.0 */
swi_precess(xyz, tjd, 0, J_TO_J2000);/**/
/* now emb -> earth */
for (i = 0; i <= 2; i++)
xemb[i] -= xyz[i] / (EARTH_MOON_MRAT + 1.0);
}
/* orbital elements of planets that are computed from osculating elements
* epoch
* equinox
* mean anomaly,
* semi axis,
* eccentricity,
* argument of perihelion,
* ascending node
* inclination
*/
#define SE_NEELY /* use James Neely's revised elements
* of Uranian planets*/
static const char *plan_fict_nam[SE_NFICT_ELEM] =
{"Cupido", "Hades", "Zeus", "Kronos",
"Apollon", "Admetos", "Vulkanus", "Poseidon",
"Isis-Transpluto", "Nibiru", "Harrington",
"Leverrier", "Adams",
"Lowell", "Pickering",};
char *swi_get_fict_name(int32 ipl, char *snam)
{
if (read_elements_file(ipl, 0, NULL, NULL,
NULL, NULL, NULL, NULL, NULL, NULL,
snam, NULL, NULL) == ERR)
strcpy(snam, "name not found");
return snam;
}
static const double plan_oscu_elem[SE_NFICT_ELEM][8] = {
#ifdef SE_NEELY
{J1900, J1900, 163.7409, 40.99837, 0.00460, 171.4333, 129.8325, 1.0833},/* Cupido Neely */
{J1900, J1900, 27.6496, 50.66744, 0.00245, 148.1796, 161.3339, 1.0500},/* Hades Neely */
{J1900, J1900, 165.1232, 59.21436, 0.00120, 299.0440, 0.0000, 0.0000},/* Zeus Neely */
{J1900, J1900, 169.0193, 64.81960, 0.00305, 208.8801, 0.0000, 0.0000},/* Kronos Neely */
{J1900, J1900, 138.0533, 70.29949, 0.00000, 0.0000, 0.0000, 0.0000},/* Apollon Neely */
{J1900, J1900, 351.3350, 73.62765, 0.00000, 0.0000, 0.0000, 0.0000},/* Admetos Neely */
{J1900, J1900, 55.8983, 77.25568, 0.00000, 0.0000, 0.0000, 0.0000},/* Vulcanus Neely */
{J1900, J1900, 165.5163, 83.66907, 0.00000, 0.0000, 0.0000, 0.0000},/* Poseidon Neely */
#else
{J1900, J1900, 104.5959, 40.99837, 0, 0, 0, 0}, /* Cupido */
{J1900, J1900, 337.4517, 50.667443, 0, 0, 0, 0}, /* Hades */
{J1900, J1900, 104.0904, 59.214362, 0, 0, 0, 0}, /* Zeus */
{J1900, J1900, 17.7346, 64.816896, 0, 0, 0, 0}, /* Kronos */
{J1900, J1900, 138.0354, 70.361652, 0, 0, 0, 0}, /* Apollon */
{J1900, J1900, -8.678, 73.736476, 0, 0, 0, 0}, /* Admetos */
{J1900, J1900, 55.9826, 77.445895, 0, 0, 0, 0}, /* Vulkanus */
{J1900, J1900, 165.3595, 83.493733, 0, 0, 0, 0}, /* Poseidon */
#endif
/* Isis-Transpluto; elements from "Die Sterne" 3/1952, p. 70ff.
* Strubell does not give an equinox. 1945 is taken to best reproduce
* ASTRON ephemeris. (This is a strange choice, though.)
* The epoch is 1772.76. The year is understood to have 366 days.
* The fraction is counted from 1 Jan. 1772 */
{2368547.66, 2431456.5, 0.0, 77.775, 0.3, 0.7, 0, 0},
/* Nibiru, elements from Christian Woeltge, Hannover */
{1856113.380954, 1856113.380954, 0.0, 234.8921, 0.981092, 103.966, -44.567, 158.708},
/* Harrington, elements from Astronomical Journal 96(4), Oct. 1988 */
{2374696.5, J2000, 0.0, 101.2, 0.411, 208.5, 275.4, 32.4},
/* Leverrier's Neptune,
according to W.G. Hoyt, "Planets X and Pluto", Tucson 1980, p. 63 */
{2395662.5, 2395662.5, 34.05, 36.15, 0.10761, 284.75, 0, 0},
/* Adam's Neptune */
{2395662.5, 2395662.5, 24.28, 37.25, 0.12062, 299.11, 0, 0},
/* Lowell's Pluto */
{2425977.5, 2425977.5, 281, 43.0, 0.202, 204.9, 0, 0},
/* Pickering's Pluto */
{2425977.5, 2425977.5, 48.95, 55.1, 0.31, 280.1, 100, 15}, /**/
#if 0 /* Ceres JPL 1600, without perturbations from other minor planets,
* from following initial elements:
* 2450600.5 2000 0 1 164.7073602 73.0340746 80.5995101
* 10.5840296 0.07652422 0.0 2.770176095 */
{2305447.5, J2000, 0.5874558977449977e+02, 0.2766536058742327e+01,
0.7870946565779195e-01, 0.5809199028919189e+02,
0.8650119410725021e+02, 0.1066835622280712e+02},
/* Chiron, Bowell database 18-mar-1997 */
{2450500.5, J2000, 7.258191, 13.67387471, 0.38174778, 339.558345, 209.379239, 6.933360}, /**/
#endif
};
/* computes a planet from osculating elements *
* tjd julian day
* ipl body number
* ipli body number in planetary data structure
* iflag flags
*/
int swi_osc_el_plan(double tjd, double *xp, int ipl, int ipli, double *xearth, double *xsun, char *serr)
{
double pqr[9], x[6];
double eps, K, fac, rho, cose, sine;
double alpha, beta, zeta, sigma, M2, Msgn, M_180_or_0;
double tjd0, tequ, mano, sema, ecce, parg, node, incl, dmot;
double cosnode, sinnode, cosincl, sinincl, cosparg, sinparg;
double M, E;
struct plan_data *pedp = &swed.pldat[SEI_EARTH];
struct plan_data *pdp = &swed.pldat[ipli];
int32 fict_ifl = 0;
int i;
/* orbital elements, either from file or, if file not found,
* from above built-in set
*/
if (read_elements_file(ipl, tjd, &tjd0, &tequ,
&mano, &sema, &ecce, &parg, &node, &incl,
NULL, &fict_ifl, serr) == ERR)
return ERR;
dmot = 0.9856076686 * DEGTORAD / sema / sqrt(sema); /* daily motion */
if (fict_ifl & FICT_GEO)
dmot /= sqrt(SUN_EARTH_MRAT);
cosnode = cos(node);
sinnode = sin(node);
cosincl = cos(incl);
sinincl = sin(incl);
cosparg = cos(parg);
sinparg = sin(parg);
/* Gaussian vector */
pqr[0] = cosparg * cosnode - sinparg * cosincl * sinnode;
pqr[1] = -sinparg * cosnode - cosparg * cosincl * sinnode;
pqr[2] = sinincl * sinnode;
pqr[3] = cosparg * sinnode + sinparg * cosincl * cosnode;
pqr[4] = -sinparg * sinnode + cosparg * cosincl * cosnode;
pqr[5] = -sinincl * cosnode;
pqr[6] = sinparg * sinincl;
pqr[7] = cosparg * sinincl;
pqr[8] = cosincl;
/* Kepler problem */
E = M = swi_mod2PI(mano + (tjd - tjd0) * dmot); /* mean anomaly of date */
/* better E for very high eccentricity and small M */
if (ecce > 0.975) {
M2 = M * RADTODEG;
if (M2 > 150 && M2 < 210) {
M2 -= 180;
M_180_or_0 = 180;
} else
M_180_or_0 = 0;
if (M2 > 330)
M2 -= 360;
if (M2 < 0) {
M2 = -M2;
Msgn = -1;
} else
Msgn = 1;
if (M2 < 30) {
M2 *= DEGTORAD;
alpha = (1 - ecce) / (4 * ecce + 0.5);
beta = M2 / (8 * ecce + 1);
zeta = pow(beta + sqrt(beta * beta + alpha * alpha), 1/3);
sigma = zeta - alpha / 2;
sigma = sigma - 0.078 * sigma * sigma * sigma * sigma * sigma / (1 + ecce);
E = Msgn * (M2 + ecce * (3 * sigma - 4 * sigma * sigma * sigma))
+ M_180_or_0;
}
}
E = swi_kepler(E, M, ecce);
/* position and speed, referred to orbital plane */
if (fict_ifl & FICT_GEO)
K = KGAUSS_GEO / sqrt(sema);
else
K = KGAUSS / sqrt(sema);
cose = cos(E);
sine = sin(E);
fac = sqrt((1 - ecce) * (1 + ecce));
rho = 1 - ecce * cose;
x[0] = sema * (cose - ecce);
x[1] = sema * fac * sine;
x[3] = -K * sine / rho;
x[4] = K * fac * cose / rho;
/* transformation to ecliptic */
xp[0] = pqr[0] * x[0] + pqr[1] * x[1];
xp[1] = pqr[3] * x[0] + pqr[4] * x[1];
xp[2] = pqr[6] * x[0] + pqr[7] * x[1];
xp[3] = pqr[0] * x[3] + pqr[1] * x[4];
xp[4] = pqr[3] * x[3] + pqr[4] * x[4];
xp[5] = pqr[6] * x[3] + pqr[7] * x[4];
/* transformation to equator */
eps = swi_epsiln(tequ, 0);
swi_coortrf(xp, xp, -eps);
swi_coortrf(xp+3, xp+3, -eps);
/* precess to J2000 */
if (tequ != J2000) {
swi_precess(xp, tequ, 0, J_TO_J2000);
swi_precess(xp+3, tequ, 0, J_TO_J2000);
}
/* to solar system barycentre */
if (fict_ifl & FICT_GEO) {
for (i = 0; i <= 5; i++) {
xp[i] += xearth[i];
}
} else {
for (i = 0; i <= 5; i++) {
xp[i] += xsun[i];
}
}
if (pdp->x == xp) {
pdp->teval = tjd; /* for precession! */
pdp->iephe = pedp->iephe;
}
return OK;
}
#if 1
/* note: input parameter tjd is required for T terms in elements */
static int read_elements_file(int32 ipl, double tjd,
double *tjd0, double *tequ,
double *mano, double *sema, double *ecce,
double *parg, double *node, double *incl,
char *pname, int32 *fict_ifl, char *serr)
{
int i, iline, iplan, retc, ncpos;
FILE *fp = NULL;
char s[AS_MAXCH], *sp;
char *cpos[20], serri[AS_MAXCH];
AS_BOOL elem_found = FALSE;
double tt = 0;
/* -1, because file information is not saved, file is always closed */
if ((fp = swi_fopen(-1, SE_FICTFILE, swed.ephepath, serr)) == NULL) {
/* file does not exist, use built-in bodies */
if (ipl >= SE_NFICT_ELEM) {
if (serr != NULL)
sprintf(serr, "error no elements for fictitious body no %7.0f", (double) ipl);
return ERR;
}
if (tjd0 != NULL)
*tjd0 = plan_oscu_elem[ipl][0]; /* epoch */
if (tequ != NULL)
*tequ = plan_oscu_elem[ipl][1]; /* equinox */
if (mano != NULL)
*mano = plan_oscu_elem[ipl][2] * DEGTORAD; /* mean anomaly */
if (sema != NULL)
*sema = plan_oscu_elem[ipl][3]; /* semi-axis */
if (ecce != NULL)
*ecce = plan_oscu_elem[ipl][4]; /* eccentricity */
if (parg != NULL)
*parg = plan_oscu_elem[ipl][5] * DEGTORAD; /* arg. of peri. */
if (node != NULL)
*node = plan_oscu_elem[ipl][6] * DEGTORAD; /* asc. node */
if (incl != NULL)
*incl = plan_oscu_elem[ipl][7] * DEGTORAD; /* inclination */
if (pname != NULL)
strcpy(pname, plan_fict_nam[ipl]);
return OK;
}
/*
* find elements in file
*/
iline = 0;
iplan = -1;
while (fgets(s, AS_MAXCH, fp) != NULL) {
iline++;
sp = s;
while(*sp == ' ' || *sp == '\t')
sp++;
swi_strcpy(s, sp);
if (*s == '#')
continue;
if (*s == '\r')
continue;
if (*s == '\n')
continue;
if (*s == '\0')
continue;
if ((sp = strchr(s, '#')) != NULL)
*sp = '\0';
ncpos = swi_cutstr(s, ",", cpos, 20);
sprintf(serri, "error in file %s, line %7.0f:", SE_FICTFILE, (double) iline);
if (ncpos < 9) {
if (serr != NULL) {
sprintf(serr, "%s nine elements required", serri);
}
goto return_err;
}
iplan++;
if (iplan != ipl)
continue;
elem_found = TRUE;
/* epoch of elements */
if (tjd0 != NULL) {
sp = cpos[0];
for (i = 0; i < 5; i++)
sp[i] = tolower(sp[i]);
if (strncmp(sp, "j2000", 5) == OK)
*tjd0 = J2000;
else if (strncmp(sp, "b1950", 5) == OK)
*tjd0 = B1950;
else if (strncmp(sp, "j1900", 5) == OK)
*tjd0 = J1900;
else if (*sp == 'j' || *sp == 'b') {
if (serr != NULL) {
sprintf(serr, "%s invalid epoch", serri);
}
goto return_err;
} else
*tjd0 = atof(sp);
tt = tjd - *tjd0;
}
/* equinox */
if (tequ != NULL) {
sp = cpos[1];
while(*sp == ' ' || *sp == '\t')
sp++;
for (i = 0; i < 5; i++)
sp[i] = tolower(sp[i]);
if (strncmp(sp, "j2000", 5) == OK)
*tequ = J2000;
else if (strncmp(sp, "b1950", 5) == OK)
*tequ = B1950;
else if (strncmp(sp, "j1900", 5) == OK)
*tequ = J1900;
else if (strncmp(sp, "jdate", 5) == OK)
*tequ = tjd;
else if (*sp == 'j' || *sp == 'b') {
if (serr != NULL) {
sprintf(serr, "%s invalid equinox", serri);
}
goto return_err;
} else
*tequ = atof(sp);
}
/* mean anomaly t0 */
if (mano != NULL) {
retc = check_t_terms(tt, cpos[2], mano);
*mano = swe_degnorm(*mano);
if (retc == ERR) {
if (serr != NULL) {
sprintf(serr, "%s mean anomaly value invalid", serri);
}
goto return_err;
}
/* if mean anomaly has t terms (which happens with fictitious
* planet Vulcan), we set
* epoch = tjd, so that no motion will be added anymore
* equinox = tjd */
if (retc == 1) {
*tjd0 = tjd;
}
*mano *= DEGTORAD;
}
/* semi-axis */
if (sema != NULL) {
retc = check_t_terms(tt, cpos[3], sema);
if (*sema <= 0 || retc == ERR) {
if (serr != NULL) {
sprintf(serr, "%s semi-axis value invalid", serri);
}
goto return_err;
}
}
/* eccentricity */
if (ecce != NULL) {
retc = check_t_terms(tt, cpos[4], ecce);
if (*ecce >= 1 || *ecce < 0 || retc == ERR) {
if (serr != NULL) {
sprintf(serr, "%s eccentricity invalid (no parabolic or hyperbolic orbits allowed)", serri);
}
goto return_err;
}
}
/* perihelion argument */
if (parg != NULL) {
retc = check_t_terms(tt, cpos[5], parg);
*parg = swe_degnorm(*parg);
if (retc == ERR) {
if (serr != NULL) {
sprintf(serr, "%s perihelion argument value invalid", serri);
}
goto return_err;
}
*parg *= DEGTORAD;
}
/* node */
if (node != NULL) {
retc = check_t_terms(tt, cpos[6], node);
*node = swe_degnorm(*node);
if (retc == ERR) {
if (serr != NULL) {
sprintf(serr, "%s node value invalid", serri);
}
goto return_err;
}
*node *= DEGTORAD;
}
/* inclination */
if (incl != NULL) {
retc = check_t_terms(tt, cpos[7], incl);
*incl = swe_degnorm(*incl);
if (retc == ERR) {
if (serr != NULL) {
sprintf(serr, "%s inclination value invalid", serri);
}
goto return_err;
}
*incl *= DEGTORAD;
}
/* planet name */
if (pname != NULL) {
sp = cpos[8];
while(*sp == ' ' || *sp == '\t')
sp++;
swi_right_trim(sp);
strcpy(pname, sp);
}
/* geocentric */
if (fict_ifl != NULL && ncpos > 9) {
for (sp = cpos[9]; *sp != '\0'; sp++)
*sp = tolower(*sp);
if (strstr(cpos[9], "geo") != NULL)
*fict_ifl |= FICT_GEO;
}
break;
}
if (!elem_found) {
if (serr != NULL) {
sprintf(serr, "%s elements for planet %7.0f not found", serri, (double) ipl);
}
goto return_err;
}
fclose(fp);
return OK;
return_err:
fclose(fp);
return ERR;
}
#endif
static int check_t_terms(double t, char *sinp, double *doutp)
{
int i, isgn = 1, z;
int retc = 0;
char *sp;
double tt[5], fac;
tt[0] = t / 36525;
tt[1] = tt[0];
tt[2] = tt[1] * tt[1];
tt[3] = tt[2] * tt[1];
tt[4] = tt[3] * tt[1];
if ((sp = strpbrk(sinp, "+-")) != NULL)
retc = 1; /* with additional terms */
sp = sinp;
*doutp = 0;
fac = 1;
z = 0;
while (1) {
while(*sp != '\0' && strchr(" \t", *sp) != NULL)
sp++;
if (strchr("+-", *sp) || *sp == '\0') {
if (z > 0)
*doutp += fac;
isgn = 1;
if (*sp == '-')
isgn = -1;
fac = 1 * isgn;
if (*sp == '\0')
return retc;
sp++;
} else {
while(*sp != '\0' && strchr("* \t", *sp) != NULL)
sp++;
if (*sp != '\0' && strchr("tT", *sp) != NULL) {
/* a T */
sp++;
if (*sp != '\0' && strchr("+-", *sp))
fac *= tt[0];
else if ((i = atoi(sp)) <= 4 && i >= 0)
fac *= tt[i];
} else {
/* a number */
if (atof(sp) != 0 || *sp == '0')
fac *= atof(sp);
}
while (*sp != '\0' && strchr("0123456789.", *sp))
sp++;
}
z++;
}
return retc; /* there have been additional terms */
}

10640
src/swemptab.h Normal file

File diff suppressed because it is too large Load Diff

2819
src/swenut2000a.h Normal file

File diff suppressed because it is too large Load Diff

341
src/sweodef.h Normal file
View File

@ -0,0 +1,341 @@
/************************************************************
definitions and constants for all Swiss Ephemeris source files,
only required for compiling the libraries, not for the external
interface of the libraries.
The definitions are a subset of Astrodienst's ourdef.h content
and must be kept compatible. Everything not used in SwissEph
has been deleted.
Does auto-detection of MSDOS (TURBO_C or MS_C), HPUNIX, Linux.
Must be extended for more portability; there should be a way
to detect byte order and file system type.
************************************************************/
/* 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.
*/
#ifndef _OURDEF_INCLUDED /* ourdef.h is a superset of sweodef.h */
#ifndef _SWEODEF_INCLUDED /* allow multiple #includes */
#define _SWEODEF_INCLUDED
# define MY_TRUE 1 /* for use in other defines, before TRUE is defined */
# define MY_FALSE 0 /* for use in other defines, before TRUE is defined */
/* TLS support
*
* Sun Studio C/C++, IBM XL C/C++, GNU C and Intel C/C++ (Linux systems) -> __thread
* Borland, VC++ -> __declspec(thread)
*/
#if !defined(TLSOFF) && !defined( __APPLE__ ) && !defined(WIN32) && !defined(DOS32)
#if defined( __GNUC__ )
#define TLS __thread
#else
#define TLS __declspec(thread)
#endif
#else
#define TLS
#endif
#ifdef _WIN32 /* Microsoft VC 5.0 does not define MSDOS anymore */
# undef MSDOS
# define MSDOS MY_TRUE
#include <wtypes.h>
#include <objbase.h>
#include <wincon.h>
#include <winbase.h>
#include <io.h>
#include <windows.h>
# define sleep(x) Sleep((x) * 1000)
#endif
#ifdef _MSC_VER
# define MS_VC
#endif
#ifdef WIN32 /* Microsoft VC 5.0 does not define MSDOS anymore */
# define MSDOS MY_TRUE
#endif
#ifdef MSDOS /* already defined by some DOS compilers */
# undef MSDOS
# define MSDOS MY_TRUE
#endif
#ifdef __TURBOC__ /* defined by turboc */
# ifndef MSDOS
# define MSDOS MY_TRUE
# endif
# define TURBO_C
#endif
#ifdef __SC__ /* defined by Symantec C */
# ifndef MSDOS
# define MSDOS MY_TRUE
# endif
# define SYMANTEC_C
#endif
#ifdef __WATCOMC__ /* defined by WatcomC */
# ifndef MSDOS
# define MSDOS MY_TRUE
# endif
# define WATCOMC
#endif
#ifdef __MWERKS__ /* defined on Macintosh CodeWarrior */
# if macintosh && powerc
# define MACOS MY_TRUE /* let it undefined otherwise */
# define MSDOS MY_FALSE /* in case one above fired falsely */
# endif
#endif
#ifdef MSDOS
# define HPUNIX MY_FALSE
# define INTEL_BYTE_ORDER 1
# ifndef TURBO_C
# define MS_C /* assume Microsoft C compiler */
# endif
# define UNIX_FS MY_FALSE
#else
# ifdef MACOS
# define HPUNIX MY_FALSE
# define UNIX_FS MY_FALSE
# else
# define MSDOS MY_FALSE
# define HPUNIX MY_TRUE
# ifndef _HPUX_SOURCE
# define _HPUX_SOURCE
# endif
# define UNIX_FS MY_TRUE
# endif
#endif
#include <math.h>
#include <stdlib.h>
#ifndef FILE
# include <stdio.h>
#endif
#if HPUNIX
# include <unistd.h>
#endif
/*
* if we have 16-bit ints, we define INT_16; we will need %ld to printf an int32
* if we have 64-bit long, we define LONG_64
* If none is defined, we have int = long = 32 bit, and use %d to printf an int32
*/
#include <limits.h>
#if INT_MAX < 40000
# define INT_16
#else
# if LONG_MAX > INT_MAX
# define LONG_64
# endif
#endif
#ifdef BYTE_ORDER
#ifdef LITTLE_ENDIAN
# if BYTE_ORDER == LITTLE_ENDIAN
# define INTEL_BYTE_ORDER
# endif
#endif
#endif
#ifdef INT_16
typedef long int32;
typedef unsigned long uint32;
typedef int int16;
typedef double REAL8; /* real with at least 64 bit precision */
typedef long INT4; /* signed integer with at least 32 bit precision */
typedef unsigned long UINT4;
/* unsigned integer with at least 32 bit precision */
typedef int AS_BOOL;
typedef unsigned int UINT2; /* unsigned 16 bits */
# define ABS4 labs /* abs function for long */
#else
typedef int int32;
typedef long long int64;
typedef unsigned int uint32;
typedef short int16;
typedef double REAL8; /* real with at least 64 bit precision */
typedef int INT4; /* signed integer with at least 32 bit precision */
typedef unsigned int UINT4;
/* unsigned integer with at least 32 bit precision */
typedef int AS_BOOL;
typedef unsigned short UINT2; /* unsigned 16 bits */
# define ABS4 abs /* abs function for long */
#endif
#if MSDOS
# ifdef TURBO_C
# include <alloc.h> /* MSC needs malloc ! */
# else
# include <malloc.h>
# endif
# define SIGALRM SIGINT
#endif
#ifndef TRUE
# define TRUE 1
# define FALSE 0
#endif
#ifndef OK
# define OK (0)
# define ERR (-1)
#endif
/* hack because UCHAR is already used by mingw gcc */
#ifdef __GNUC__
#ifdef _WIN32
#define UCHAR SWE_UCHAR
#endif
#endif
typedef unsigned char UCHAR;
#define UCP (UCHAR*)
#define SCP (char*)
# define ODEGREE_STRING "°" /* degree as string, utf8 encoding */
#ifndef HUGE
# define HUGE 1.7E+308 /* biggest value for REAL8 */
#endif
#ifndef M_PI
# define M_PI 3.14159265358979323846
#endif
/* #define forward static obsolete */
#define AS_MAXCH 256 /* used for string declarations, allowing 255 char+\0 */
/*
#define DEGTORAD 0.0174532925199433
#define RADTODEG 57.2957795130823
*/
#define RADTODEG (180.0 / M_PI)
#define DEGTORAD (M_PI / 180.0)
typedef int32 centisec; /* centiseconds used for angles and times */
#define CS (centisec) /* use for casting */
#define CSEC centisec /* use for typing */
#define DEG 360000 /* degree expressed in centiseconds */
#define DEG7_30 (2700000) /* 7.5 degrees */
#define DEG15 (15 * DEG)
#define DEG24 (24 * DEG)
#define DEG30 (30 * DEG)
#define DEG60 (60 * DEG)
#define DEG90 (90 * DEG)
#define DEG120 (120 * DEG)
#define DEG150 (150 * DEG)
#define DEG180 (180 * DEG)
#define DEG270 (270 * DEG)
#define DEG360 (360 * DEG)
/* #define CSTORAD 4.84813681109536E-08 centisec to rad: pi / 180 /3600/100 */
/* #define RADTOCS 2.06264806247096E+07 rad to centisec 180*3600*100/pi */
#define CSTORAD (DEGTORAD / 360000.0)
#define RADTOCS (RADTODEG * 360000.0)
#define CS2DEG (1.0/360000.0) /* centisec to degree */
/* control strings for fopen() */
#if UNIX_FS
# define BFILE_R_ACCESS "r" /* open binary file for reading */
# define BFILE_RW_ACCESS "r+" /* open binary file for writing and reading */
# define BFILE_W_CREATE "w" /* create/open binary file for write*/
# define BFILE_A_ACCESS "a+" /* create/open binary file for append*/
# define FILE_R_ACCESS "r" /* open text file for reading */
# define FILE_RW_ACCESS "r+" /* open text file for writing and reading */
# define FILE_W_CREATE "w" /* create/open text file for write*/
# define FILE_A_ACCESS "a+" /* create/open text file for append*/
# define O_BINARY 0 /* for open(), not defined in Unix */
# define OPEN_MODE 0666 /* default file creation mode */
# define DIR_GLUE "/" /* glue string for directory/file */
# define PATH_SEPARATOR ";:" /* semicolon or colon may be used */
#else
# define BFILE_R_ACCESS "rb" /* open binary file for reading */
# define BFILE_RW_ACCESS "r+b" /* open binary file for writing and reading */
# define BFILE_W_CREATE "wb" /* create/open binary file for write*/
# define BFILE_A_ACCESS "a+b" /* create/open binary file for append*/
# define PATH_SEPARATOR ";" /* semicolon as PATH separator */
# define OPEN_MODE 0666 /* default file creation mode */
# ifdef MACOS
# define FILE_R_ACCESS "r" /* open text file for reading */
# define FILE_RW_ACCESS "r+" /* open text file for writing and reading */
# define FILE_W_CREATE "w" /* create/open text file for write*/
# define FILE_A_ACCESS "a+" /* create/open text file for append*/
# define DIR_GLUE ":" /* glue string for directory/file */
# else
# define FILE_R_ACCESS "rt" /* open text file for reading */
# define FILE_RW_ACCESS "r+t" /* open text file for writing and reading */
# define FILE_W_CREATE "wt" /* create/open text file for write*/
# define FILE_A_ACCESS "a+t" /* create/open text file for append*/
/* attention, all backslashes for msdos directry names must be written as \\,
because it is the C escape character */
# define DIR_GLUE "\\" /* glue string for directory/file */
# endif
#endif
#include <string.h>
#include <ctype.h>
#endif /* _SWEODEF_INCLUDED */
#endif /* _OURDEF_INCLUDED */

Some files were not shown because too many files have changed in this diff Show More