python-tdd/pyswiss/apps/charts/calc.py

"""
Core ephemeris calculation logic — shared by views and management commands.
"""
from django.conf import settings as django_settings
import swisseph as swe


DEFAULT_HOUSE_SYSTEM = 'O'  # Porphyry

SIGNS = [
    'Aries', 'Taurus', 'Gemini', 'Cancer', 'Leo', 'Virgo',
    'Libra', 'Scorpio', 'Sagittarius', 'Capricorn', 'Aquarius', 'Pisces',
]

SIGN_ELEMENT = {
    'Aries': 'Fire',  'Leo': 'Fire',   'Sagittarius': 'Fire',
    'Taurus': 'Earth', 'Virgo': 'Earth', 'Capricorn': 'Earth',
    'Gemini': 'Air',  'Libra': 'Air',  'Aquarius': 'Air',
    'Cancer': 'Water', 'Scorpio': 'Water', 'Pisces': 'Water',
}

ASPECTS = [
    ('Conjunction',      0,   8.0),
    ('Semisextile',     30,   4.0),
    ('Semisquare',      45,   4.0),
    ('Sextile',         60,   6.0),
    ('Square',          90,   8.0),
    ('Trine',          120,   8.0),
    ('Sesquiquadrate', 135,   4.0),
    ('Quincunx',       150,   5.0),
    ('Opposition',     180,  10.0),
]

PLANET_CODES = {
    'Sun':     swe.SUN,
    'Moon':    swe.MOON,
    'Mercury': swe.MERCURY,
    'Venus':   swe.VENUS,
    'Mars':    swe.MARS,
    'Jupiter': swe.JUPITER,
    'Saturn':  swe.SATURN,
    'Uranus':  swe.URANUS,
    'Neptune': swe.NEPTUNE,
    'Pluto':   swe.PLUTO,
}


def set_ephe_path():
    ephe_path = getattr(django_settings, 'SWISSEPH_PATH', None)
    if ephe_path:
        swe.set_ephe_path(ephe_path)


def get_sign(lon):
    return SIGNS[int(lon // 30) % 12]


def get_julian_day(dt):
    return swe.julday(
        dt.year, dt.month, dt.day,
        dt.hour + dt.minute / 60 + dt.second / 3600,
    )


def get_planet_positions(jd):
    flag = swe.FLG_SWIEPH | swe.FLG_SPEED
    planets = {}
    for name, code in PLANET_CODES.items():
        pos, _ = swe.calc_ut(jd, code, flag)
        degree = pos[0]
        planets[name] = {
            'sign': get_sign(degree),
            'degree': degree,
            'speed': pos[3],
            'retrograde': pos[3] < 0,
        }
    return planets


# ── Reverse lookup — sign windows (Set the Game Clock narrowing) ─────────────
#
# Forward calculation answers "where is the planet at time T"; the Game Clock
# ritual needs the REVERSE — "when does planet P reside in sign S" — so each
# placement can narrow the real date window the ritual still resolves to.
#
# Scan strides (days) sit well under each planet's shortest typical sign
# residence. Retrograde re-dips briefer than the stride can be missed — that
# is CONSERVATIVE: a missed sliver only narrows the reachable window, it never
# admits a false moment (every returned window is forward-verified by its
# construction: sampled in-sign, edges bisected to the true crossing).

SIGN_SCAN_STRIDES = {
    'Moon': 0.5, 'Sun': 5.0, 'Mercury': 2.0, 'Venus': 2.0, 'Mars': 5.0,
    'Jupiter': 15.0, 'Saturn': 30.0, 'Uranus': 60.0, 'Neptune': 60.0,
    'Pluto': 60.0,
}
EDGE_PRECISION_DAYS = 1.0 / 24.0   # bisect sign crossings to the hour


def get_planet_sign(jd, planet):
    pos, _ = swe.calc_ut(jd, PLANET_CODES[planet], swe.FLG_SWIEPH)
    return get_sign(pos[0])


def jd_to_iso(jd):
    y, m, d, h = swe.revjul(jd)
    hh = int(h)
    mm = int(round((h - hh) * 60))
    if mm == 60:
        hh, mm = hh + 1, 0
    return f'{y:04d}-{m:02d}-{d:02d}T{hh:02d}:{mm:02d}:00Z'


def _bisect_edge(planet, sign, jd_match, jd_other):
    """Refine the sign crossing between a jd where `planet` IS in `sign` and
    one where it is not, to EDGE_PRECISION_DAYS. Returns the match-side jd."""
    while abs(jd_other - jd_match) > EDGE_PRECISION_DAYS:
        mid = (jd_match + jd_other) / 2
        if get_planet_sign(mid, planet) == sign:
            jd_match = mid
        else:
            jd_other = mid
    return jd_match


def sign_windows(planet, sign, windows):
    """The sub-windows of `windows` — a list of (jd_start, jd_end) — where
    `planet` resides in `sign`. Stride scan + bisection edge refine."""
    stride = SIGN_SCAN_STRIDES[planet]
    out = []
    for lo, hi in windows:
        if hi <= lo:
            continue
        samples = []
        jd = lo
        while jd < hi:
            samples.append(jd)
            jd += stride
        samples.append(hi)

        run_start = None
        prev = None
        for s in samples:
            in_sign = get_planet_sign(s, planet) == sign
            if in_sign and run_start is None:
                run_start = lo if prev is None else _bisect_edge(planet, sign, s, prev)
            elif not in_sign and run_start is not None:
                out.append((run_start, _bisect_edge(planet, sign, prev, s)))
                run_start = None
            prev = s
        if run_start is not None:
            out.append((run_start, hi))
    return out


def narrowed_windows(placements, windows):
    """Fold each placement's sign residences into the intersection. Slowest
    planets first (fewest fragments); each scan is restricted to the windows
    the prior placements already narrowed to, so the fast bodies only ever
    scan slivers."""
    for planet, sign in sorted(
            placements.items(), key=lambda kv: -SIGN_SCAN_STRIDES[kv[0]]):
        if not windows:
            break
        windows = sign_windows(planet, sign, windows)
    return windows


def present_signs(planet, windows):
    """The set of signs `planet` occupies at any point within `windows` —
    sampled at the scan stride plus both edges (presence needs no bisection)."""
    stride = SIGN_SCAN_STRIDES[planet]
    signs = set()
    for lo, hi in windows:
        jd = lo
        while jd < hi:
            signs.add(get_planet_sign(jd, planet))
            jd += stride
        signs.add(get_planet_sign(hi, planet))
    return signs


def get_element_counts(planets):
    sign_counts = {s: 0 for s in SIGNS}
    sign_planets = {s: [] for s in SIGNS}
    classic = {'Fire': [], 'Water': [], 'Earth': [], 'Air': []}

    for name, data in planets.items():
        sign = data['sign']
        el = SIGN_ELEMENT[sign]
        classic[el].append({'planet': name, 'sign': sign})
        sign_counts[sign] += 1
        sign_planets[sign].append({'planet': name, 'sign': sign})

    result = {
        el: {'count': len(contribs), 'contributors': contribs}
        for el, contribs in classic.items()
    }

    # Time: stellium — highest concentration in one sign, bonus = size - 1.
    # Collect all signs tied at the maximum.
    max_in_sign = max(sign_counts.values())
    stellia = [
        {'sign': s, 'planets': sign_planets[s]}
        for s in SIGNS
        if sign_counts[s] == max_in_sign and max_in_sign > 1
    ]
    result['Time'] = {
        'count': max_in_sign - 1,
        'stellia': stellia,
    }

    # Space: parade — longest consecutive run of occupied signs (circular),
    # bonus = run length - 1. Collect all runs tied at the maximum.
    index_set = {i for i, s in enumerate(SIGNS) if sign_counts[s] > 0}
    indices = sorted(index_set)
    max_seq = 0
    for start in range(len(indices)):
        seq_len = 1
        for offset in range(1, len(indices)):
            if (indices[start] + offset) % len(SIGNS) in index_set:
                seq_len += 1
            else:
                break
        max_seq = max(max_seq, seq_len)

    parades = []
    for start in range(len(indices)):
        run = []
        for offset in range(max_seq):
            idx = (indices[start] + offset) % len(SIGNS)
            if idx not in index_set:
                break
            run.append(idx)
        else:
            sign_run = [SIGNS[i] for i in run]
            parade_planets = [
                p for s in sign_run for p in sign_planets[s]
            ]
            parades.append({'signs': sign_run, 'planets': parade_planets})

    result['Space'] = {
        'count': max_seq - 1,
        'parades': parades,
    }

    return result


def calculate_aspects(planets):
    """Return a list of aspects between all planet pairs.

    Each entry: {planet1, planet2, type, angle (actual, rounded), orb (rounded)}.
    Only the first matching aspect type is reported per pair (aspects are
    well-separated enough that at most one can apply with standard orbs).
    """
    names = list(planets.keys())
    aspects = []
    for i, name1 in enumerate(names):
        for name2 in names[i + 1:]:
            deg1 = planets[name1]['degree']
            deg2 = planets[name2]['degree']
            angle = abs(deg1 - deg2)
            if angle > 180:
                angle = 360 - angle
            for aspect_name, target, max_orb in ASPECTS:
                orb = abs(angle - target)
                if orb <= max_orb:
                    s1 = abs(planets[name1].get('speed', 0))
                    s2 = abs(planets[name2].get('speed', 0))
                    applying = name1 if s1 >= s2 else name2
                    aspects.append({
                        'planet1': name1,
                        'planet2': name2,
                        'type': aspect_name,
                        'angle': round(angle, 2),
                        'orb': round(orb, 2),
                        'applying_planet': applying,
                    })
                    break
    return aspects