jesse-uxandGrok (agent loop)
a757af6b4e
10 张虚构盘逐体系对照。不一致的方法在标题下写明差在哪一类。本命 Tithi 顺行一致，逆行没有公开函数。Windows 上 1970 年以前的日出不再用 fromtimestamp。 自测：合同与研究章测试通过；200 张英文扫描约 1956 秒通过。未跑全量 pytest、快速门和前端。 Co-Authored-By: Grok (agent loop)
270 lines
10 KiB
Python
270 lines
10 KiB
Python
"""Return roots and elongation clocks on SiderealFrame.
|
|
|
|
The step and bisection follow the research return search. The longitude
|
|
call is ours. Windows stay at most 400 days.
|
|
"""
|
|
|
|
from __future__ import annotations
|
|
|
|
import bisect
|
|
from datetime import datetime, timedelta, timezone
|
|
|
|
import swisseph as swe
|
|
|
|
try:
|
|
from research_timing.frame import SiderealFrame
|
|
except ImportError: # pragma: no cover - package import
|
|
from scripts.research_timing.frame import SiderealFrame
|
|
|
|
try:
|
|
from muhurta import _find_next_boundary
|
|
except ImportError: # pragma: no cover - package import
|
|
from scripts.muhurta import _find_next_boundary
|
|
|
|
_STEP = 0.125
|
|
_BISECTIONS = 45
|
|
_TOLERANCE = 1e-6
|
|
_UNIX_EPOCH_JD = 2440587.5
|
|
|
|
|
|
def _utc_from_jd(jd: float) -> datetime:
|
|
"""UTC from a Julian day. Windows ``fromtimestamp`` rejects instants before 1970."""
|
|
unix = (float(jd) - _UNIX_EPOCH_JD) * 86400.0
|
|
seconds = int(unix)
|
|
micros = int(round((unix - seconds) * 1_000_000))
|
|
if micros >= 1_000_000:
|
|
seconds += 1
|
|
micros -= 1_000_000
|
|
elif micros < 0:
|
|
seconds -= 1
|
|
micros += 1_000_000
|
|
return datetime(1970, 1, 1, tzinfo=timezone.utc) + timedelta(seconds=seconds, microseconds=micros)
|
|
|
|
|
|
def find_returns(
|
|
frame: SiderealFrame,
|
|
technique: str,
|
|
birth_jd: float,
|
|
start_jd: float,
|
|
end_jd: float,
|
|
) -> dict:
|
|
if technique not in ("yoga_pravesha", "nakshatra_pravesha", "tithi_pravesha"):
|
|
raise ValueError(f"unsupported return technique: {technique}")
|
|
if not start_jd < end_jd or end_jd - start_jd > 400:
|
|
raise ValueError("search window must be positive and no longer than 400 days")
|
|
sun, moon = frame.sun_moon(birth_jd)
|
|
solar_sign = int(sun // 30)
|
|
target = _coordinate(technique, sun, moon)
|
|
|
|
def residual(jd: float) -> float:
|
|
left, right = frame.sun_moon(jd)
|
|
value = _coordinate(technique, left, right)
|
|
return (value - target + 180.0) % 360.0 - 180.0
|
|
|
|
rows: list[dict] = []
|
|
left = float(start_jd)
|
|
left_value = residual(left)
|
|
while left < end_jd:
|
|
right = min(left + _STEP, float(end_jd))
|
|
right_value = residual(right)
|
|
if left_value * right_value <= 0 and abs(right_value - left_value) < 180:
|
|
root = left if abs(left_value) < 1e-8 else _bisect(left, right, left_value, residual)
|
|
body_sun, body_moon = frame.sun_moon(root)
|
|
error = residual(root)
|
|
if (
|
|
start_jd <= root < end_jd
|
|
and abs(error) < _TOLERANCE
|
|
and int(body_sun // 30) == solar_sign
|
|
and all(abs(root - row["jd_ut"]) > _TOLERANCE for row in rows)
|
|
):
|
|
rows.append({
|
|
"jd_ut": root,
|
|
"sun_longitude": body_sun,
|
|
"moon_longitude": body_moon,
|
|
"residual_degrees": error,
|
|
})
|
|
left, left_value = right, right_value
|
|
return {
|
|
"technique": technique,
|
|
"target_degrees": target,
|
|
"natal_solar_sign_index": solar_sign,
|
|
"natal_sun": sun,
|
|
"natal_moon": moon,
|
|
"returns": rows,
|
|
"frame": "swiss_ephemeris_true_precession",
|
|
}
|
|
|
|
|
|
def _coordinate(technique: str, sun: float, moon: float) -> float:
|
|
if technique == "tithi_pravesha":
|
|
return (moon - sun) % 360.0
|
|
if technique == "yoga_pravesha":
|
|
return (sun + moon) % 360.0
|
|
return moon % 360.0
|
|
|
|
|
|
def _bisect(left: float, right: float, left_value: float, residual) -> float:
|
|
low, high, low_value = left, right, left_value
|
|
for _ in range(_BISECTIONS):
|
|
mid = (low + high) / 2.0
|
|
value = residual(mid)
|
|
if low_value * value <= 0:
|
|
high = mid
|
|
else:
|
|
low, low_value = mid, value
|
|
return (low + high) / 2.0
|
|
|
|
|
|
def jd_to_utc(jd: float) -> datetime:
|
|
year, month, day, hour = swe.revjul(float(jd), swe.GREG_CAL)
|
|
return datetime(int(year), int(month), int(day), tzinfo=timezone.utc) + timedelta(hours=float(hour))
|
|
|
|
|
|
def containing_tithi(frame: SiderealFrame, jd: float) -> dict:
|
|
"""Time fraction left in the 12-degree tithi that contains this instant."""
|
|
|
|
def phase(at: float) -> float:
|
|
sun, moon = frame.sun_moon(at)
|
|
return (moon - sun) % 360.0
|
|
|
|
angle = phase(jd)
|
|
cursor = jd - 1.6
|
|
previous = None
|
|
following = None
|
|
for _ in range(4):
|
|
boundary = _find_next_boundary(cursor, phase, 12.0, 1.6)
|
|
if boundary is None or boundary <= cursor:
|
|
raise RuntimeError("containing tithi boundaries unavailable")
|
|
if boundary > jd:
|
|
following = boundary
|
|
break
|
|
previous, cursor = boundary, boundary
|
|
if previous is None or following is None or not previous <= jd < following:
|
|
raise RuntimeError("containing tithi boundaries unavailable")
|
|
return {
|
|
"tithi_num": int(angle // 12) + 1,
|
|
"elongation": angle,
|
|
"remaining_time_fraction": (following - jd) / (following - previous),
|
|
"remaining_angular_fraction": 1.0 - (angle % 12.0) / 12.0,
|
|
"start_jd_ut": previous,
|
|
"end_jd_ut": following,
|
|
}
|
|
|
|
|
|
def measure_elongation(frame: SiderealFrame, start_jd: float, end_jd: float) -> float:
|
|
"""Unwrapped sidereal Moon-minus-Sun over at most 400 days."""
|
|
if not 0 < end_jd - start_jd <= 400:
|
|
raise ValueError("motion interval must span more than 0 and at most 400 days")
|
|
cursor = float(start_jd)
|
|
sun, moon = frame.sun_moon(cursor)
|
|
previous = (moon - sun) % 360.0
|
|
total = 0.0
|
|
while cursor < end_jd:
|
|
cursor = min(cursor + 1.0, float(end_jd))
|
|
sun, moon = frame.sun_moon(cursor)
|
|
current = (moon - sun) % 360.0
|
|
total += (current - previous + 180.0) % 360.0 - 180.0
|
|
previous = current
|
|
return total
|
|
|
|
|
|
class ElongationClock:
|
|
"""Map a 360-unit year onto accumulated sidereal elongation."""
|
|
|
|
def __init__(self, frame: SiderealFrame, return_jd: float, total_angle: float, known_end_jd: float | None = None) -> None:
|
|
if not 0 < total_angle <= 6000:
|
|
raise ValueError("positive total angle within the annual grid is required")
|
|
self.frame = frame
|
|
self.return_jd = float(return_jd)
|
|
self.total_angle = float(total_angle)
|
|
self._days = [self.return_jd + index for index in range(-100, 501)]
|
|
self._phases = [self._phase(day) for day in self._days]
|
|
self._angles = [0.0]
|
|
for left, right in zip(self._phases, self._phases[1:]):
|
|
step = (right - left + 180.0) % 360.0 - 180.0
|
|
if not 0 < step < 180:
|
|
raise ValueError("elongation clock requires forward daily motion")
|
|
self._angles.append(self._angles[-1] + step)
|
|
self._origin = self._angles[100]
|
|
self._cache = {0.0: self.return_jd}
|
|
if known_end_jd is not None:
|
|
if not self.return_jd < known_end_jd <= self.return_jd + 500:
|
|
raise ValueError("known endpoint must lie in the forward annual grid")
|
|
self._cache[360.0] = float(known_end_jd)
|
|
|
|
def _phase(self, jd: float) -> float:
|
|
sun, moon = self.frame.sun_moon(jd)
|
|
return (moon - sun) % 360.0
|
|
|
|
def solve(self, offset: float) -> float:
|
|
if abs(offset) < 1e-8:
|
|
offset = 0.0
|
|
if abs(offset - 360.0) < 1e-8:
|
|
offset = 360.0
|
|
cached = self._cache.get(offset)
|
|
if cached is not None:
|
|
return cached
|
|
target = self._origin + offset / 360.0 * self.total_angle
|
|
index = bisect.bisect_right(self._angles, target) - 1
|
|
if not 0 <= index < len(self._days) - 1:
|
|
raise ValueError("angular boundary outside the bounded search grid")
|
|
found = self._solve_day(index, target)
|
|
self._cache[offset] = found
|
|
return found
|
|
|
|
def _solve_day(self, index: int, target: float) -> float:
|
|
"""Day bracket, then a few slope steps. Residual stays within 1e-6 degree."""
|
|
left, right = self._days[index], self._days[index + 1]
|
|
left_angle = self._angles[index]
|
|
span = self._angles[index + 1] - left_angle
|
|
if span <= 0 or right <= left:
|
|
raise ValueError("elongation clock day is not forward")
|
|
origin_phase = self._phases[index]
|
|
|
|
def motion(jd: float) -> float:
|
|
return left_angle + (self._phase(jd) - origin_phase + 180.0) % 360.0 - 180.0
|
|
|
|
slope = span / (right - left)
|
|
guess = left + (target - left_angle) / span * (right - left)
|
|
for _ in range(6):
|
|
value = motion(guess)
|
|
error = value - target
|
|
if abs(error) <= 1e-6:
|
|
return guess
|
|
guess = min(right, max(left, guess - error / slope))
|
|
low, high = left, right
|
|
for _ in range(40):
|
|
mid = (low + high) / 2.0
|
|
error = motion(mid) - target
|
|
if abs(error) <= 1e-6 or high - low < 1e-9:
|
|
return mid
|
|
if error < 0:
|
|
low = mid
|
|
else:
|
|
high = mid
|
|
return (low + high) / 2.0
|
|
|
|
|
|
def hindu_solar_events(moment: datetime, latitude: float, longitude: float) -> list[dict]:
|
|
"""Sunrise and sunset from the Swiss ephemeris file. No Moshier and no true-position flag."""
|
|
if moment.utcoffset() is None:
|
|
raise ValueError("solar events need a timezone")
|
|
midnight = moment.astimezone(timezone.utc).replace(hour=0, minute=0, second=0, microsecond=0)
|
|
flags = swe.FLG_SWIEPH
|
|
events = []
|
|
for offset in range(-3, 4):
|
|
jd = (midnight + timedelta(days=offset)).timestamp() / 86400.0 + 2440587.5
|
|
for name, mode in (("sunrise", swe.CALC_RISE), ("sunset", swe.CALC_SET)):
|
|
status, times = swe.rise_trans(
|
|
jd,
|
|
swe.SUN,
|
|
mode | swe.BIT_HINDU_RISING,
|
|
(float(longitude), float(latitude), 0.0),
|
|
flags=flags,
|
|
)
|
|
if status == 0 and jd <= times[0] < jd + 1:
|
|
stamp = _utc_from_jd(times[0])
|
|
events.append({"event": name, "datetime_utc": stamp.isoformat()})
|
|
events.sort(key=lambda row: row["datetime_utc"])
|
|
return events
|