Git Product home page Git Product logo

swiss-ephemeris's Introduction

swiss-ephemeris

build

Haskell bindings for the Swiss Ephemeris library.

See the tests in the spec folder for thorough example usage, but here's a simple "main" that demonstrates the current abilities, inspired by the sample program in the official library:

NOTE: this library is under very active development, as such, most releases in v1.x will probably show a fastly evolving API, which is reflected by the fact that new versions have been increasing the major version numbers (in PVP, unlike semver, the first two components of the version correspond to the major version.)

import SwissEphemeris

main :: IO
main = do 
  -- location of your ephemeris directory. We bundle a sample one in `swedist`.
  withEphemerides "./swedist/sweph_18" $ do
    let time = julianDay 1989 1 6 0.0
        place = GeographicPosition {geoLat = 14.0839053, geoLng = -87.2750137}
  
    -- locate all bodies between the Sun and Chiron
    forM_ [Sun .. Chiron] $ \planet -> do
      -- if no ephemerides data is available for the given planetary body, a `Left` value
      -- will be returned.
      coord <- calculateEclipticPosition time planet
      putStrLn $ show planet <> ": " <> show coord
   
    -- Calculate cusps for the given time and place, preferring the `Placidus` system.
    -- note that the underlying library may decide to use the `Porphyrius` system if it can't
    -- calculate cusps (happens for the Placidus and Koch systems in locations near the poles.)
    cusps <- calculateCusps Placidus time place
    putStrLn $ "Cusps: " <> show cusps

The above should print the ecliptic latitude and longitude (plus some velocities) for all planets, and then the cusps and other major angles (ascendant, mc, ARMC, alternative angles.)

There's withEphemerides to run calculations using a particular ephemerides directory and then close any used system resources, and withoutEphemerides to use the default ephemerides ("Moshier.")

To see actual results and more advanced usage, check out the tests. For some more advanced examples, see swetest.c and swemini.c in the csrc directory: they're the test/example programs provided by the original authors. You can also play around with the C library via the authors' test page.

Notes

All the code in the csrc folder comes directly from the latest official tarball (as of June 2021:) v2.10.01.

The swedist folder includes the original documentation from the tarball in PDF (see the doc) folder, and a copy of ephemeris data files.

For other formats of the original documentation, see: https://www.astro.com/ftp/swisseph/doc/

The authors also host HTML versions of the manuals. Two are provided, a general reference and a programming reference. Both are very useful to get acquainted with the functionality and implementation details.

Ephemerides files

As noted in the original documentation you can omit the setEphemerides call (or use setNoEphemerides, or the withoutEphemerides bracket function) and calculations will use a built-in analytical ephemeris ("Moshier") which:

provides "only" a precision of 0.1 arc seconds for the planets and 3" for the Moon. No asteroids will be available, and no barycentric option can be used.

Note that if you're interested in the asteroid Chiron (which is common in astrological practice these days,) you'll have to procure Ephemerides files and shouldn't use the default ephemerides.

For convenience, we bundle a few ephemerides files in this repository (see swedist) for the time range 1800 AD – 2399 AD. If you were born before that, or plan to code e.g. transits for after that (!) or you'd prefer even more precision, you can download more ephemerides files from the astro.com downloads page

I chose the bundled files due to this comment in the official docs:

If the [JPL] file is too big, then you can download the files sepl_18.se1 and semo_18.se1 from here: http://www.astro.com/ftp/swisseph/ephe/

For a full explanation of the files available, see the Description of the Ephemerides section of the original manual, also of interest is the comparison between the Swiss Ephemeris and the raw NASA JPL data.

Contributing

I've only made available the types and functions that are useful for my own, traditional, horoscope calculations. Feel free to add more! See the astro.com documentation for ideas.

Releasing

To generate docs (to proofread locally,) as well as the tarball for upload to hackage, run nix-build ./nix/release.nix.

swiss-ephemeris's People

Contributors

lfborjas avatar

Stargazers

 avatar  avatar  avatar  avatar  avatar  avatar

Watchers

 avatar  avatar

swiss-ephemeris's Issues

Consider adding equatorial coordinates, for declinations

Some explanation of declinations:

https://astro-charts.com/blog/2017/easy-way-calculate-parallels-contraparallels/

And https://www.astro.com/swisseph/swephprg.htm, the section 3.3. Options chosen by flag bits (long iflag), and the explanation in 3.4 Position and speed:

Ecliptic position Equatorial position (SEFLG_EQUATORIAL)
Longitude right ascension
Latitude declination
Distance in AU distance in AU
Speed in longitude (deg/day) speed in right ascension (deg/day)
Speed in latitude (deg/day) speed in declination (deg/day)
Speed in distance (AU/day) speed in distance (AU/day)

If you need rectangular coordinates (SEFLG_XYZ), swe_calc() returns x, y, z, dx, dy, dz in AU.

Once you have computed a planet, e.g., in ecliptic coordinates, its equatorial position or its rectangular coordinates are available, too. You can get them very cheaply (little CPU time used), calling again swe_calc() with the same parameters, but adding SEFLG_EQUATORIAL or SEFLG_XYZ to iflag, swe_calc() will not compute the body again, just return the data specified from internal storage.

And here's how to set the flags:

-- | Combine @CalcFlag@ as bitwise options. 
-- For example, can use it to obtain the equatorial,
-- not eliptical (default,) positions of planetary bodies.
calculationOptions :: [CalcFlag] -> CalcFlag
calculationOptions = CalcFlag . foldr ((.|.) . unCalcFlag) 0

(from an older version of this library, funnily enough.)

CI fails in ubuntu often, I suspect a segfault, have no idea how to diagnose

See the toils and travails in #6 -- many failures, some of those commits definitely deserved to fail, but even after I cleaned up the code considerably, I see "segfaults". Not sure if it's the C library doing something untoward, but it happens extremely often in Ubuntu and never in macOS (I develop in macOS, and the tests work a frustrating 99.999% of the time, though I've managed to see errors too.) I just don't know enough C, Haskell or CI to get to the bottom of this, and I spent an inordinate amount of hours the past two nights staring at this.

Add bracketed root finding for known small ranges of crossing

Some context:

Potential optimizations:

incorrect cusp

Hi,

am trying this version of swiss ehim code but am getting incorrect cusp values. i had an old version of swiss ephim 2.01-00.02 this version gave accurate cusp. could you please let me know what can cause cusp values change?

Fix warnings in C code

csrc/configurable_sweephe4.c: In function ‘eph4_posit’:

csrc/configurable_sweephe4.c:555:21: error:
     warning: ‘%s’ directive writing 5 bytes into a region of size between 1 and 256 [-Wformat-overflow=]
      555 |       sprintf(s, "%s%sM%d", ephe4d.ephe4path, EP4_FILE, -filenr);
          |                     ^~
    |
555 |       sprintf(s, "%s%sM%d", ephe4d.ephe4path, EP4_FILE, -filenr);
    |                     ^

csrc/configurable_sweephe4.c:555:18: error:
     note: directive argument in the range [1, 214749]
      555 |       sprintf(s, "%s%sM%d", ephe4d.ephe4path, EP4_FILE, -filenr);
          |                  ^~~~~~~~~
    |
555 |       sprintf(s, "%s%sM%d", ephe4d.ephe4path, EP4_FILE, -filenr);
    |                  ^
In file included from /usr/include/stdio.h:867,
                 from csrc/sweodef.h:169,
                 from csrc/swephexp.h:82,

                 from csrc/configurable_sweephe4.c:60:0: error: 

/usr/include/x86_64-linux-gnu/bits/stdio2.h:36:10: error:
     note: ‘__builtin___sprintf_chk’ output between 8 and 268 bytes into a destination of size 256
       36 |   return __builtin___sprintf_chk (__s, __USE_FORTIFY_LEVEL - 1,
    |                     ^

csrc/configurable_sweephe4.c:553:18: error:
     note: directive argument in the range [0, 214748]
      553 |       sprintf(s, "%s%s%d", ephe4d.ephe4path, EP4_FILE, filenr);
          |                  ^~~~~~~~
    |
553 |       sprintf(s, "%s%s%d", ephe4d.ephe4path, EP4_FILE, filenr);
    |                  ^
In file included from /usr/include/stdio.h:867,
                 from csrc/sweodef.h:169,
                 from csrc/swephexp.h:82,

                 from csrc/configurable_sweephe4.c:60:0: error: 

/usr/include/x86_64-linux-gnu/bits/stdio2.h:36:10: error:
     note: ‘__builtin___sprintf_chk’ output between 7 and 267 bytes into a destination of size 256
       36 |   return __builtin___sprintf_chk (__s, __USE_FORTIFY_LEVEL - 1,
          |          ^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
       37 |       __bos (__s), __fmt, __va_arg_pack ());
          |       ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
   |
36 |   return __builtin___sprintf_chk (__s, __USE_FORTIFY_LEVEL - 1,
   |          ^
csrc/configurable_sweephe4.c: In function ‘ephe4_write_file’:

csrc/configurable_sweephe4.c:990:61: error:
     warning: ‘%s’ directive output may be truncated writing up to 255 bytes into a region of size 235 [-Wformat-truncation=]
      990 |           snprintf(errtext, AS_MAXCH, "error in swe_calc(): %s", serr);
          |                                                             ^~   ~~~~
    |
990 |           snprintf(errtext, AS_MAXCH, "error in swe_calc(): %s", serr);
    |                                                             ^
In file included from /usr/include/stdio.h:867,
                 from csrc/sweodef.h:169,
                 from csrc/swephexp.h:82,

                 from csrc/configurable_sweephe4.c:60:0: error: 

/usr/include/x86_64-linux-gnu/bits/stdio2.h:67:10: error:
     note: ‘__builtin___snprintf_chk’ output between 22 and 277 bytes into a destination of size 256
       67 |   return __builtin___snprintf_chk (__s, __n, __USE_FORTIFY_LEVEL - 1,
          |          ^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
       68 |        __bos (__s), __fmt, __va_arg_pack ());
          |        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
   |
67 |   return __builtin___snprintf_chk (__s, __n, __USE_FORTIFY_LEVEL - 1,
   |          ^

Consider "parabola fit" for crossings, also, double-crossings

See:

/* y00, y11, y2 are values for -2*dx, -dx, 0.
* find zero points of parabola.
* return: 0 if none
* 1 if one zero in [-dx.. 0[
* 1 if both zeros in [-dx.. 0[
*/
static int find_zero(double y00, double y11, double y2, double dx,
double *dxret, double *dxret2)
{
double a, b, c, x1, x2;
c = y11;
b = (y2 - y00) / 2.0;
a = (y2 + y00) / 2.0 - c;
if (b * b - 4 * a * c < 0)
return 0;
if (fabs(a) < 1e-100) return 0;
x1 = (-b + sqrt(b * b - 4 * a * c)) / 2 / a;
x2 = (-b - sqrt(b * b - 4 * a * c)) / 2 / a;
// up to here the calcuation was made as if the x-values were -1, 0, 1.
// This is why below they are shifted by -1
if (x1 == x2) {
*dxret = (x1 - 1) * dx;
*dxret2 = (x1 - 1) * dx;
return 1;
}
if (x1 >=0 && x1 < 1 && x2 >= 0 && x2 < 1) {
if (x1 > x2) { // two zeroes, order return values
*dxret = (x2 - 1) * dx;
*dxret2 = (x1 - 1) * dx;
} else {
*dxret = (x1 - 1) * dx;
*dxret2 = (x2 - 1) * dx;
}
return 2;
}
if (x1 >=0 && x1 < 1) {
*dxret = (x1 - 1) * dx;
*dxret2 = (x2 - 1) * dx; // set this value just in case, should not be used.
return 1;
}
if (x2 >=0 && x2 < 1) {
*dxret = (x2 - 1) * dx;
*dxret2 = (x1 - 1) * dx;
return 1;
}
return 0; // should not happen!
}

A "parabola fit" is often mentioned in the swisseph forum, and I believe it's using the above. Unlike interpolation, it requires bracketing and some cajoling of values into it to account for "circle jumps," and out of it to actually find the equivalent date value; the gain being that it doesn't require as many ephemeris IO events as regular interpolation. This method is also able to detect double crossings, which we currently don't do (though the tools are there! See the comments for crossingBetween)

Handle failures to calculate cusps

I hadn't noticed that the C library can return -1 when unable to calculate cusps! See the java impl here for reference:
https://github.com/Kibo/AstroAPI/blob/7f2e884abea354aba6ff8c44b4db2aa110085b8f/src/main/java/cz/kibo/api/astrology/domain/Cusp.java#L97 (which references: https://github.com/Kibo/AstroAPI/blob/cfc3d1e5a4c1d8dd5cd9906b15cec855e9e6302e/src/main/java/swisseph/SweConst.java#L95)

It seems fairly unlikely to happen though, more for sunshine houses in polar regions, maybe?

int CALL_CONV swe_houses(double tjd_ut,
double geolat,
double geolon,
int hsys,
double *cusp,
double *ascmc)
{
int i, retc = 0;
double armc, eps, nutlo[2];
double tjde = tjd_ut + swe_deltat_ex(tjd_ut, -1, NULL);
eps = swi_epsiln(tjde, 0) * RADTODEG;
swi_nutation(tjde, 0, nutlo);
for (i = 0; i < 2; i++)
nutlo[i] *= RADTODEG;
armc = swe_degnorm(swe_sidtime0(tjd_ut, eps + nutlo[1], nutlo[0]) * 15 + geolon);
if (toupper(hsys) == 'I') { // compute sun declination for sunshine houses
int flags = SEFLG_SPEED| SEFLG_EQUATORIAL;
double xp[6];
int result = swe_calc_ut(tjd_ut, SE_SUN, flags, xp, NULL);
if (result < 0) {
// in case of failure, Porphyry houses
result = swe_houses_armc(armc, geolat, eps + nutlo[1], 'O', cusp, ascmc);
return ERR;
}
ascmc[9] = xp[1]; // declination in ascmc[9];
}
#ifdef TRACE
swi_open_trace(NULL);
if (swi_trace_count <= TRACE_COUNT_MAX) {
if (swi_fp_trace_c != NULL) {
fputs("\n/*SWE_HOUSES*/\n", swi_fp_trace_c);
fprintf(swi_fp_trace_c, "#if 0\n");
fprintf(swi_fp_trace_c, " tjd = %.9f;", tjd_ut);
fprintf(swi_fp_trace_c, " geolon = %.9f;", geolon);
fprintf(swi_fp_trace_c, " geolat = %.9f;", geolat);
fprintf(swi_fp_trace_c, " hsys = %d;\n", hsys);
fprintf(swi_fp_trace_c, " retc = swe_houses(tjd, geolat, geolon, hsys, cusp, ascmc);\n");
fprintf(swi_fp_trace_c, " /* swe_houses calls swe_houses_armc as follows: */\n");
fprintf(swi_fp_trace_c, "#endif\n");
fflush(swi_fp_trace_c);
}
}
#endif
retc = swe_houses_armc(armc, geolat, eps + nutlo[1], hsys, cusp, ascmc);
return retc;
}

Miscellaneous minor improvements

Library

  • Include documentation for the various types, in particular house systems.
  • Include support for the few celestial bodies past Chiron -- they're already in the C bindings, I just didn't have a need for the others! won't do, yet!
  • As mentioned in #1 , add tests for calculating cusps using other systems than Placidus. (done in #5)
  • Consider adding the scientific package to palliate precision issues in other build environments? (see note in README) Doubles are fine.
  • Consider adding a withEphemeridesPath :: a -> IO a function that will run a given calculation after setting an ephemerides path, and then close said path? Would be a bit more expensive than the authors intended, but safer on resources. (done in #5)
  • Make CI more deterministic: due to polar computations sometimes succeeding for Placidus and Koch, our property test of them failing seems to sometimes not hold! (this is why the build fails sometimes??) -- this was not why CI failed sometimes: it was a bona fide segfault caused by a memory leak. See notes in #9

Include namespace to work with precalculated ephemeris

The latest versions of swiss ephemeris include the following files, which allow for generation and consumption of optimized (for storage and reading) "precalculated" ephemeris. They're useful for transits and other use cases where several days of ephemeris may be examined at once:

I believe they need a little bit of tweaking to be able to generate and read the files from a non-hardcoded base path, in addition to the documentation/testing of the Haskell bindings. It would be interesting to benchmark them against calls to the actual ephemeris calculation functionality.

Introduction by the authors at: https://groups.io/g/swisseph/message/8864

Note that this is part of a thread about transit calculation, which seems to be the main use case.

Add dependency on `time` for UTCTime->JulianTime conversions

Edit

I believe any difference between functions is a misunderstanding. From my observations, calculating the julianTime yields a JulianTime in UTC, to convert to "terrestrial time", which accounts for the non-uniformity of Earth's rotation (more in the Delta T section of the official docs.)

The "Terrestrial Time" (TT) that the swiss ephemeris test page reports is then: julianTime + deltaTime. We should provide some conversions, namely UTCTime -> JulianTime, JulianTime -> TerrestrialTime (or EphemerisTime?), and UTCTime -> TerrestrialTime. The first of those would require either porting swe_utc_to_jd, or introducing a dependency on the time package, to replicate the below:

Original

I'm currently using the following in my application code, which glues Haskell's UTCTime with JulianTime via our very own julianDay:

-- | Convert between a UTC timestamp and the low-level JulianTime that SwissEphemeris requires.
utcToJulian :: UTCTime -> JulianTime
utcToJulian (UTCTime day time) = 
    julianDay (fromIntegral $ y) m d h
    where
        (y, m, d) = toGregorian day
        h         = 2.77778e-16 * (fromIntegral $ diffTimeToPicoseconds time)

julianToUTC :: JulianTime -> UTCTime
juliantoUTC jd = 
  UTCTime (fromGregorian y m d) (picosecondsToDiffTime $ round $ h * 3600 * 1e12)
  where
     (y, m, d, h) = gregorianDateTime jd

Which seems to yield pretty good results according to JPL's conversion tool (though a few seconds behind Swisseph!) for some test data:

import Data.Time
utc <- (parseTimeM True defaultTimeLocale "%Y-%-m-%-d %T" "1989-12-25 09:30:00" :: IO UTCTime)
utcToJulian tsUTC
2447885.89583365

-- and back
> julianToUTC $ JulianTime 2447885.89583365
1989-12-25 09:30:00.027371942996 UTC

JPL's data:

image

-- the advantage is that it's a pure function, the disadvantage is that it seems to be running afoul of some history-aware corrections that the swe_utc_to_jd function makes:

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)

The binding to this function would have to be in IO, and deal with some funky pointers; but it's clearly very well thought out and a useful utility to expose. One of the original authors says that the difference between the functions is only relevant in specific situations, so maybe the above method is okay for most Haskell users? 🤔

Change/improve signature of `calculateCusps`?

The house system is currently the last argument, perhaps it should be the first? That way people can say things like:

calculateCuspsP = calculateCusps Placidus

Additionally, we should add tests for the other popular systems (as per https://www.astro.com/swisseph/swephprg.htm#_Toc46406845).

I considered making the houses not be a record with 12 keys, but a simple list -- however, all popular house systems (except Gauquelin) have 12 houses, and each house has specific semantics. Perhaps a helper function that transform them to a list would be good, for those who only wish to iterate?

Make buildable in nix

It currently fails to build due to:

Configuring swiss-ephemeris-1.2.1.1...
CallStack (from HasCallStack):
  $, called at libraries/Cabal/Cabal/Distribution/Simple/Configure.hs:1024:20 in Cabal-3.2.1.0:Distribution.Simple.Configure
  configureFinalizedPackage, called at libraries/Cabal/Cabal/Distribution/Simple/Configure.hs:477:12 in Cabal-3.2.1.0:Distribution.Simple.Configure
  configure, called at libraries/Cabal/Cabal/Distribution/Simple.hs:625:20 in Cabal-3.2.1.0:Distribution.Simple
  confHook, called at libraries/Cabal/Cabal/Distribution/Simple/UserHooks.hs:65:5 in Cabal-3.2.1.0:Distribution.Simple.UserHooks
  configureAction, called at libraries/Cabal/Cabal/Distribution/Simple.hs:180:19 in Cabal-3.2.1.0:Distribution.Simple
  defaultMainHelper, called at libraries/Cabal/Cabal/Distribution/Simple.hs:116:27 in Cabal-3.2.1.0:Distribution.Simple
  defaultMain, called at Setup.hs:2:8 in main:Main
Setup: Encountered missing or private dependencies:
QuickCheck >=2.12 && <=2.14```

Recommend Projects

  • React photo React

    A declarative, efficient, and flexible JavaScript library for building user interfaces.

  • Vue.js photo Vue.js

    🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.

  • Typescript photo Typescript

    TypeScript is a superset of JavaScript that compiles to clean JavaScript output.

  • TensorFlow photo TensorFlow

    An Open Source Machine Learning Framework for Everyone

  • Django photo Django

    The Web framework for perfectionists with deadlines.

  • D3 photo D3

    Bring data to life with SVG, Canvas and HTML. 📊📈🎉

Recommend Topics

  • javascript

    JavaScript (JS) is a lightweight interpreted programming language with first-class functions.

  • web

    Some thing interesting about web. New door for the world.

  • server

    A server is a program made to process requests and deliver data to clients.

  • Machine learning

    Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.

  • Game

    Some thing interesting about game, make everyone happy.

Recommend Org

  • Facebook photo Facebook

    We are working to build community through open source technology. NB: members must have two-factor auth.

  • Microsoft photo Microsoft

    Open source projects and samples from Microsoft.

  • Google photo Google

    Google ❤️ Open Source for everyone.

  • D3 photo D3

    Data-Driven Documents codes.