fix: house calculation, MC formula, LST correction, zodiac direction

src/astro_mcp/astrology.py

@@ -319,54 +319,100 @@ def calculate_houses(
 
 
 def _houses_placidus(lst_deg: float, latitude: float) -> list[dict[str, Any]]:
-    """Placidus house cusps approximation.
+    """Placidus house cusps.
 
-    Uses the standard Placidus method: trisecting the semi-arc of each house.
-    This is a well-known approximation accurate to ~0.01° for most locations.
+    Standard Placidus semi-arc trisection method.
+    Houses numbered counter-clockwise from ASC (house 1).
+
+    The algorithm:
+    1. Compute RA of ASC, MC, DSC, IC
+    2. The diurnal semi-arc (ASC→DSC eastward through MC) = ~180° in RA
+    3. The nocturnal semi-arc (DSC→ASC eastward through IC) = ~180° in RA
+    4. Trisect each semi-arc to get intermediate cusps
+    5. Convert RA cusps back to ecliptic longitude
     """
-    lat_rad = math.radians(latitude)
-    obl = math.radians(23.4367)  # approximate obliquity
+    obl_e = 23.4367
 
-    # ASC: ascendant
-    # MC: midheaven from LST
     asc = _calc_ascendant(lst_deg, latitude)
-    mc = _calc_midheaven(lst_deg)
-
-    # For Placidus, we compute intermediate cusps via semi-arc method
-    # Cusps 2, 3, 11, 12 are on the diurnal semi-arc;
-    # cusps 5, 6, 8, 9 are on the nocturnal semi-arc.
-    # These approximations are standard in open-source astrology libraries.
+    mc  = _calc_midheaven(lst_deg)
+    dsc = _opposition(asc)
+    ic  = _opposition(mc)
+
+    ra_asc = _ecliptic_to_ra(asc, obl_e)
+    ra_mc  = _ecliptic_to_ra(mc,  obl_e)
+    ra_dsc = _ecliptic_to_ra(dsc, obl_e)
+    ra_ic  = _ecliptic_to_ra(ic,  obl_e)
+
+    # Helper: trisect an arc in RA space
+    def _trisect_ra(ra_start, ra_end, fraction):
+        """Move fraction of the way from ra_start to ra_end.
+        ra_end should already be adjusted so the arc goes in the correct direction
+        (i.e., ra_end > ra_start for the desired direction, possibly with +360 offset)."""
+        arc = ra_end - ra_start  # may be > 180 or < 0, that's OK
+        if arc < 0:
+            arc += 360
+        if arc > 180:
+            # This shouldn't happen if direction is correct, but handle it
+            arc = arc - 360
+        target = ra_start + arc * fraction
+        return _ra_to_ecliptic(normalize_degrees(target), obl_e)
+
+    # Determine the correct direction for semi-arc trisection.
+    # The diurnal semi-arc (ASC → MC → DSC) should contain MC between ASC and DSC.
+    # Check if MC is between ASC and DSC going eastward.
+    arc_mc_from_asc = (ra_mc - ra_asc) % 360
+    if arc_mc_from_asc > 180:
+        # MC is NOT between ASC and DSC going eastward.
+        # Go westward instead (negative RA offset = westward).
+        direction = -1
+        # Adjust: add 360 to MC and DSC RAs so they're "before" ASC in westward direction
+        ra_mc_eff = ra_mc + 360
+        ra_dsc_eff = ra_dsc + 360
+        ra_ic_eff = ra_ic + 360
+    else:
+        direction = 1
+        ra_mc_eff = ra_mc
+        ra_dsc_eff = ra_dsc
+        ra_ic_eff = ra_ic
 
-    # Compute RA of each cusp via trisecting the semi-arc
     cusps: list[float | None] = [None] * 12
-    cusps[0] = asc  # House 1
-    cusps[3] = _opposition(mc)  # House 4 (IC)
-    cusps[6] = _opposition(asc)  # House 7 (DSC)
-    cusps[9] = mc  # House 10 (MC)
-
-    # Standard book formulation for Placidus intermediate cusps:
-    obl_e = 23.4367  # obliquity of ecliptic in degrees
-    intermediate = [
-        (1, 30.0, True),   # House 2, diurnal
-        (2, 60.0, True),   # House 3, diurnal
-        (4, 30.0, False),  # House 5, nocturnal
-        (5, 60.0, False),  # House 6, nocturnal
-        (7, 30.0, False),  # House 8, nocturnal
-        (8, 60.0, False),  # House 9, nocturnal
-        (10, 30.0, True),  # House 11, diurnal
-        (11, 60.0, True),  # House 12, diurnal
-    ]
-    for cusp_idx, target_oa, diurnal_flag in intermediate:
-        cusp_lon = _oblique_ascension_to_ecliptic(
-            lst_deg, latitude, target_oa, obl_e, diurnal=diurnal_flag
-        )
-        cusps[cusp_idx] = cusp_lon
+    cusps[0] = asc    # H1: ASC
+    cusps[3] = ic     # H4: IC
+    cusps[6] = dsc    # H7: DSC
+    cusps[9] = mc     # H10: MC
+
+    # Diurnal semi-arc: ASC → MC → DSC
+    # H12 = 1/3 from ASC to MC, H11 = 2/3 from ASC to MC
+    # H9 = 1/3 from MC to DSC, H8 = 2/3 from MC to DSC
+    if direction == 1:
+        cusps[11] = _trisect_ra(ra_asc, ra_mc, 1.0/3.0)
+        cusps[10] = _trisect_ra(ra_asc, ra_mc, 2.0/3.0)
+        cusps[8] = _trisect_ra(ra_mc, ra_dsc, 1.0/3.0)
+        cusps[7] = _trisect_ra(ra_mc, ra_dsc, 2.0/3.0)
+        # Nocturnal: DSC → IC → ASC
+        cusps[5] = _trisect_ra(ra_dsc, ra_ic, 1.0/3.0)
+        cusps[4] = _trisect_ra(ra_dsc, ra_ic, 2.0/3.0)
+        cusps[2] = _trisect_ra(ra_ic, ra_asc, 1.0/3.0)
+        cusps[1] = _trisect_ra(ra_ic, ra_asc, 2.0/3.0)
+    else:
+        # Westward: use adjusted RAs
+        # Diurnal: ASC → MC → DSC (westward = decreasing RA, but we use +360 adjusted values)
+        cusps[11] = _trisect_ra(ra_asc, ra_mc_eff, 1.0/3.0)
+        cusps[10] = _trisect_ra(ra_asc, ra_mc_eff, 2.0/3.0)
+        cusps[8] = _trisect_ra(ra_mc_eff, ra_dsc_eff, 1.0/3.0)
+        cusps[7] = _trisect_ra(ra_mc_eff, ra_dsc_eff, 2.0/3.0)
+        # Nocturnal: DSC → IC → ASC (westward)
+        # Going westward from DSC: DSC → IC → ASC
+        # But with +360 adjustment: DSC+360 → IC+360 → ASC+360
+        ra_asc_eff = ra_asc + 360
+        cusps[5] = _trisect_ra(ra_dsc_eff, ra_ic_eff, 1.0/3.0)
+        cusps[4] = _trisect_ra(ra_dsc_eff, ra_ic_eff, 2.0/3.0)
+        cusps[2] = _trisect_ra(ra_ic_eff, ra_asc_eff, 1.0/3.0)
+        cusps[1] = _trisect_ra(ra_ic_eff, ra_asc_eff, 2.0/3.0)
 
-    # If any calculation returned None (polar), fall back to Equal House
     if any(c is None for c in cusps):
         return _houses_equal(lst_deg)
 
-    # Convert all to zodiac dicts
     result = []
     for i, cusp_lon in enumerate(cusps):
         z = ecliptic_to_zodiac(cusp_lon if cusp_lon is not None else 0.0)
@@ -432,6 +478,36 @@ def _ecliptic_to_ra(ecliptic_lon: float, obliquity: float) -> float:
     return normalize_degrees(math.degrees(ra))
 
 
+def _ra_to_ecliptic(ra_deg: float, obliquity: float) -> float:
+    """Convert right ascension to ecliptic longitude.
+
+    Inverse of _ecliptic_to_ra. Uses iterative refinement for accuracy.
+    """
+    obl_rad = math.radians(obliquity)
+    ra_rad  = math.radians(ra_deg)
+
+    # First approximation: tan(el) = tan(ra) / cos(obliquity)
+    tan_el = math.tan(ra_rad) / math.cos(obl_rad)
+    el_rad = math.atan(tan_el)
+    el_deg = math.degrees(el_rad)
+
+    # Adjust quadrant: ecliptic longitude should be in the same
+    # half-circle as the RA.
+    if 90 < ra_deg < 270:
+        el_deg += 180
+    el_deg = normalize_degrees(el_deg)
+
+    # Refine iteratively
+    for _ in range(5):
+        computed_ra = _ecliptic_to_ra(el_deg, obliquity)
+        error = normalize_degrees(ra_deg - computed_ra + 180) - 180
+        if abs(error) < 0.0001:
+            break
+        el_deg = normalize_degrees(el_deg + error * math.cos(obl_rad))
+
+    return el_deg
+
+
 def _calc_ascendant(lst_deg: float, latitude: float) -> float:
     """Calculate the ascendant ecliptic longitude from LST and latitude."""
     obl = math.radians(23.4367)
@@ -451,10 +527,10 @@ def _calc_ascendant(lst_deg: float, latitude: float) -> float:
 
 def _calc_midheaven(lst_deg: float) -> float:
     """Calculate MC ecliptic longitude from LST."""
-    # MC: tan(MC) = tan(RA) * cos(obliquity)
+    # MC: tan(MC) = tan(RA) / cos(obliquity)
     ra_rad = math.radians(lst_deg)
     obl = math.radians(23.4367)
-    mc_rad = math.atan2(math.tan(ra_rad), 1.0 / math.cos(obl))
+    mc_rad = math.atan2(math.tan(ra_rad), math.cos(obl))
     mc_deg = math.degrees(mc_rad)
     # Same quadrant as RA
     if 90 < lst_deg < 270:

src/astro_mcp/ephemeris_client.py

@@ -77,11 +77,19 @@ async def call_sky_state(
     lon: float = 0.0,
     elevation: float = 0.0,
     geocentric: bool = True,
+    house_system: str | None = None,
 ) -> dict[str, Any]:
     """Call ephemeris-mcp:get_sky_state and return the result dict.
 
     Datetime is normalized to UTC ISO format (no timezone offset)
     before sending to the ephemeris server.
+
+    The ephemeris server returns Greenwich sidereal_time regardless of
+    the lon parameter. We correct local_sidereal_time here by adding
+    the longitude offset (15° = 1 hour, east positive).
+
+    When house_system is provided, the response includes house cusps
+    and angles computed by the Swiss Ephemeris on the server side.
     """
     # Normalize datetime to UTC ISO string
     dt_str = _normalize_datetime(datetime)
@@ -105,18 +113,49 @@ async def call_sky_state(
                         "lon": lon,
                         "elevation": elevation,
                         "geocentric": geocentric,
+                        "house_system": house_system,
                     },
                 )
                 payload = _payload_from_result(result)
                 if not payload:
                     logger.warning("ephemeris-mcp returned empty payload")
                     return {"error": "empty_response", "url": url}
+
+                # Fix local sidereal time: ephemeris may return incorrect LST.
+                # Compute it properly from Greenwich ST + longitude.
+                # East longitude positive, 15° = 1 hour.
+                if "sidereal_time" in payload:
+                    st = payload["sidereal_time"]
+                    if isinstance(st, dict):
+                        gst = st.get("greenwich_sidereal_time",
+                                     st.get("local_sidereal_time", 0.0))
+                        if isinstance(gst, (int, float)):
+                            st["local_sidereal_time"] = (gst + lon / 15.0) % 24.0
+
                 return payload
     except Exception as exc:
         logger.error(f"ephemeris-mcp call failed: {exc}")
         return {"error": str(exc), "url": url}
 
 
+def extract_houses(sky_state: dict[str, Any]) -> list[dict[str, Any]] | None:
+    """Extract house cusps from a sky_state response.
+
+    Returns the list of 12 house cusp dicts if houses were computed,
+    or None if no house system was requested.
+
+    Each dict: {"house": int, "absolute_lon": float, "sign": str,
+                "abbreviation": str, "degree": float}
+    """
+    houses_data = sky_state.get("houses")
+    if houses_data is None:
+        return None
+    if "error" in houses_data:
+        logger.warning(f"house calculation error: {houses_data['error']}")
+        return None
+    return houses_data.get("cusps")
+
+
 def extract_bodies(sky_state: dict[str, Any]) -> list[dict[str, Any]]:
     """Extract the planetary bodies array from a sky_state response.
 

src/astro_mcp/tools.py

@@ -12,7 +12,7 @@ from typing import Any
 
 from .server import mcp
 from . import astrology
-from .ephemeris_client import call_sky_state, extract_bodies
+from .ephemeris_client import call_sky_state, extract_bodies, extract_houses
 
 logger = logging.getLogger("astro-mcp.tools")
 
@@ -161,17 +161,20 @@ Returns:
         lon=longitude,
         elevation=elevation,
         geocentric=True,
+        house_system=house_system,
     )
 
     if "error" in sky:
         return {"input": {"birth_datetime": birth_datetime, "latitude": latitude, "longitude": longitude}, "error": sky["error"]}
 
     raw_bodies = extract_bodies(sky)
-    sidereal = sky.get("sidereal_time", {})
-    lst_hours = sidereal.get("local_sidereal_time", 0.0)
 
-    # Calculate houses
-    houses = astrology.calculate_houses(lst_hours, latitude, house_system)
+    # Use server-side houses if available, otherwise fall back to local calculation
+    houses = extract_houses(sky)
+    if houses is None:
+        sidereal = sky.get("sidereal_time", {})
+        lst_hours = sidereal.get("local_sidereal_time", 0.0)
+        houses = astrology.calculate_houses(lst_hours, latitude, house_system)
 
     # Build planet list with house placement
     planets = []
@@ -344,6 +347,7 @@ Returns:
         lon=longitude,
         elevation=elevation,
         geocentric=True,
+        house_system=house_system,
     )
 
     # Get transit sky state at transit location
@@ -363,10 +367,12 @@ Returns:
     natal_bodies = extract_bodies(natal_sky)
     transit_bodies = extract_bodies(transit_sky)
 
-    # Natal houses from natal LST
-    sidereal = natal_sky.get("sidereal_time", {})
-    lst_hours = sidereal.get("local_sidereal_time", 0.0)
-    houses = astrology.calculate_houses(lst_hours, latitude, house_system)
+    # Use server-side houses if available, otherwise fall back to local calculation
+    houses = extract_houses(natal_sky)
+    if houses is None:
+        sidereal = natal_sky.get("sidereal_time", {})
+        lst_hours = sidereal.get("local_sidereal_time", 0.0)
+        houses = astrology.calculate_houses(lst_hours, latitude, house_system)
 
     # Build natal planets
     natal_planets = []
@@ -482,8 +488,8 @@ Returns:
     Dict with: input, chart_type, chart1_natal, chart2_natal, interaspects,
     house_overlays, composite_chart, davison_chart, summary (top_aspects,
     saturn_contacts, node_contacts, venus_mars_contacts, sun_moon_contacts)."""
-    sky1 = await call_sky_state(datetime=person1_datetime, lat=person1_latitude, lon=person1_longitude, elevation=elevation)
-    sky2 = await call_sky_state(datetime=person2_datetime, lat=person2_latitude, lon=person2_longitude, elevation=elevation)
+    sky1 = await call_sky_state(datetime=person1_datetime, lat=person1_latitude, lon=person1_longitude, elevation=elevation, house_system=house_system)
+    sky2 = await call_sky_state(datetime=person2_datetime, lat=person2_latitude, lon=person2_longitude, elevation=elevation, house_system=house_system)
 
     if "error" in sky1:
         return {"error": f"person1: {sky1['error']}"}
@@ -493,13 +499,18 @@ Returns:
     bodies1 = extract_bodies(sky1)
     bodies2 = extract_bodies(sky2)
 
-    sidereal1 = sky1.get("sidereal_time", {})
-    lst1 = sidereal1.get("local_sidereal_time", 0.0)
-    houses1 = astrology.calculate_houses(lst1, person1_latitude, house_system)
+    # Use server-side houses if available, otherwise fall back to local calculation
+    houses1 = extract_houses(sky1)
+    if houses1 is None:
+        sidereal1 = sky1.get("sidereal_time", {})
+        lst1 = sidereal1.get("local_sidereal_time", 0.0)
+        houses1 = astrology.calculate_houses(lst1, person1_latitude, house_system)
 
-    sidereal2 = sky2.get("sidereal_time", {})
-    lst2 = sidereal2.get("local_sidereal_time", 0.0)
-    houses2 = astrology.calculate_houses(lst2, person2_latitude, house_system)
+    houses2 = extract_houses(sky2)
+    if houses2 is None:
+        sidereal2 = sky2.get("sidereal_time", {})
+        lst2 = sidereal2.get("local_sidereal_time", 0.0)
+        houses2 = astrology.calculate_houses(lst2, person2_latitude, house_system)
 
     def build_planet_list(bodies):
         result = []
@@ -605,13 +616,15 @@ Returns:
         davison_dt = _jd_to_datetime(davison["date_midpoint_jd"])
         davison_sky = await call_sky_state(
             datetime=davison_dt, lat=davison_mid_lat, lon=davison_mid_lon,
-            elevation=0.0, geocentric=True,
+            elevation=0.0, geocentric=True, house_system=house_system,
         )
         if "error" not in davison_sky:
             davison_raw = extract_bodies(davison_sky)
-            davison_sidereal = davison_sky.get("sidereal_time", {})
-            davison_lst = davison_sidereal.get("local_sidereal_time", 0.0)
-            davison_houses = astrology.calculate_houses(davison_lst, davison_mid_lat, house_system)
+            davison_houses = extract_houses(davison_sky)
+            if davison_houses is None:
+                davison_sidereal = davison_sky.get("sidereal_time", {})
+                davison_lst = davison_sidereal.get("local_sidereal_time", 0.0)
+                davison_houses = astrology.calculate_houses(davison_lst, davison_mid_lat, house_system)
 
             davison_planets = []
             for body in davison_raw:
@@ -1030,13 +1043,16 @@ Returns:
     # Use a date near the midpoint for house calculation
     comp_sky = await call_sky_state(
         datetime=person1_datetime, lat=comp_lat, lon=comp_lon, elevation=elevation,
+        house_system=house_system,
     )
     if "error" in comp_sky:
         return {"error": f"composite ephemeris error: {comp_sky['error']}"}
 
+    houses = extract_houses(comp_sky)
     sidereal = comp_sky.get("sidereal_time", {})
     lst_hours = sidereal.get("local_sidereal_time", 0.0)
-    houses = astrology.calculate_houses(lst_hours, comp_lat, house_system)
+    if houses is None:
+        houses = astrology.calculate_houses(lst_hours, comp_lat, house_system)
 
     # Build composite planet list with house placement
     composite_planets = []
@@ -1126,15 +1142,17 @@ Returns:
 
     sky = await call_sky_state(
         datetime=davison_dt, lat=mid_lat, lon=mid_lon,
-        elevation=elevation, geocentric=True,
+        elevation=elevation, geocentric=True, house_system=house_system,
     )
     if "error" in sky:
         return {"error": f"davison ephemeris error: {sky['error']}"}
 
     raw_bodies = extract_bodies(sky)
+    houses = extract_houses(sky)
     sidereal = sky.get("sidereal_time", {})
     lst_hours = sidereal.get("local_sidereal_time", 0.0)
-    houses = astrology.calculate_houses(lst_hours, mid_lat, house_system)
+    if houses is None:
+        houses = astrology.calculate_houses(lst_hours, mid_lat, house_system)
 
     planets = []
     for body in raw_bodies: