lfborjas / swiss-ephemeris Goto Github PK

7.13. swe_pheno_ut() and swe_pheno()

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?

Some context:

• https://groups.io/g/swisseph/topic/27493138#3009
• https://groups.io/g/swisseph/message/8036 + https://groups.io/g/swisseph/message/8037
• https://groups.io/g/swisseph/message/8861

Potential optimizations:

• Instead of calling swe_calc on each iteration, call it only four times: three points to draw a parabola, and one final time after pure approximation on that parabola is done. See the find_zero function, and usage for e.g. zodiac crossings (note that the code that does the > 180 check is now immortalized in the difdeg2n function.
• Consider a more sophisticated implementation of brent-dekker, from scipy: https://github.com/scipy/scipy/blob/39bf11b96f771dcecf332977fb2c7843a9fd55f2/scipy/optimize/Zeros/brentq.c#L37

See: https://groups.io/g/swisseph/topic/27496438#5564 here with syntax highlighting

Introduction in: https://groups.io/g/swisseph/message/5560

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```

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.

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:

• swephgen4.c
• sweephe4.h
• sweephe4.c

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.

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?

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

https://groups.io/g/swisseph/topic/new_release_2_10_02/84663960?p=,,,20,0,0,0::recentpostdate%2Fsticky,,,20,2,0,84663960

https://www.astro.com/swisseph/swephprg.htm#_Toc78973568

See:

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)

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:

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.)

• Assess the precision (or loss thereof) when converting between double types: https://ro-che.info/articles/2019-05-14-convert-cdouble-to-double
• Could ameliorate this by using scientific to return values, and let downstream consumers convert to other floats: https://hackage.haskell.org/package/scientific-0.3.7.0/docs/Data-Scientific.html#v:fromFloatDigits

See:
https://groups.io/g/swisseph/topic/swiss_ephemeris_version_2_10/78852296?p=,,,20,0,0,0::recentpostdate%2Fsticky,,,20,2,0,78852296

and

https://www.astro.com/swisseph/swephprg.htm?lang=g#_Toc58481317

The current behavior includes second fractions, meaning that people who wish to not display second fractions have to do rounding up themselves!

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:

-- 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:

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? 🤔

See: https://groups.io/g/swisseph/message/7781 -- which is a follow up of: https://groups.io/g/swisseph/topic/64346793#7779

And: https://gist.github.com/lfborjas/ccde203851906f4b81305f401d03787c

See: lfborjas/timezone-detect#4 -- one-off users should be able to run a computation that needs the files to be open and later closed, so we should provide a withEphePath bracket for those use cases.

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?

I believe it's because the swedist folder, at least the ephemeris, is not included. Test it by unzipping and trying to run tests.