horizonEvents function

List<HorizonEvent> horizonEvents(
  1. Ephemeris ephemeris,
  2. int target,
  3. double latitude,
  4. double longitude,
  5. double startJdTT,
  6. double endJdTT, {
  7. double? horizonAltitude,
})

Risings, settings and meridian transits of body target for a place at latitude and geographic longitude (degrees, east positive) in startJdTT … endJdTT, in time order.

The body rises or sets when its geocentric centre is at horizonAltitude degrees. By default that is the almanac convention: sunHorizonAltitude for the Sun, starHorizonAltitude for planets, and for the Moon 0.7275·π − 34′, where π is its horizontal parallax (J. Meeus, Astronomical Algorithms, chapter 15). Refraction is a standard value; real refraction near the horizon varies by several arcminutes, i.e. by tens of seconds or more in time.

Near the poles a body may stay above or below the horizon: then the list has only transits.

Implementation

List<HorizonEvent> horizonEvents(
  Ephemeris ephemeris,
  int target,
  double latitude,
  double longitude,
  double startJdTT,
  double endJdTT, {
  double? horizonAltitude,
}) {
  double h0(EquatorialPosition p) =>
      horizonAltitude ??
      switch (target) {
        Naif.sun => sunHorizonAltitude,
        Naif.moon =>
          0.7275 * math.asin(earthEquatorialRadius / p.distance) / _radian +
              starHorizonAltitude,
        _ => starHorizonAltitude,
      };
  double hourAngleAt(double t) =>
      _hourAngle(equatorialOfDate(ephemeris, target, t), longitude, t);
  double aboveHorizon(double t) {
    final p = equatorialOfDate(ephemeris, target, t);
    return _altitude(p, _hourAngle(p, longitude, t), latitude) - h0(p);
  }

  // Culminations: hour angle 0° (upper) and 180° (lower). Between them
  // the altitude changes in one direction, so each piece has at most one
  // rising or setting.
  final transits = longitudeCrossings(
    hourAngleAt,
    startJdTT,
    endJdTT,
    period: 180,
    step: 0.25,
  );
  final events = <HorizonEvent>[
    for (final c in transits)
      HorizonEvent(
        c.jdTT,
        c.longitude < 90
            ? HorizonEventKind.upperTransit
            : HorizonEventKind.lowerTransit,
      ),
  ];
  final cuts = [startJdTT, for (final c in transits) c.jdTT, endJdTT];
  for (var i = 0; i + 1 < cuts.length; i++) {
    final a = cuts[i], b = cuts[i + 1];
    final fa = aboveHorizon(a), fb = aboveHorizon(b);
    if (fa < 0 == fb < 0) continue;
    events.add(
      HorizonEvent(
        _bisect(aboveHorizon, a, b, fa),
        fa < 0 ? HorizonEventKind.rise : HorizonEventKind.set,
      ),
    );
  }
  return events..sort((x, y) => x.jdTT.compareTo(y.jdTT));
}