test_single_planet.py

#!/usr/bin/env python3
"""
Single Planet Exomoon Detection Test Script
==========================================
Test the exomoon detection pipeline on a single planet before running
on the full dataset.

Usage:
    python test_single_planet.py "WASP-10 b"
    python test_single_planet.py --tic 123456789
    python test_single_planet.py --name "WASP-10" --period 3.09 --radius 3.15
"""

import os
import sys
import numpy as np
import pandas as pd
import argparse
import warnings
warnings.filterwarnings('ignore')

from lightkurve import search_lightcurve, LightCurveCollection
from astropy.timeseries import BoxLeastSquares, LombScargle
from astropy import units as u
from astroquery.nasa_exoplanet_archive import NasaExoplanetArchive

# Import helper functions from comprehensive script
# (We'll include them here for standalone use)

def validate_parameters(P, dur, t, flux):
    """Validate computed parameters are physically reasonable."""
    if np.isnan(P) or P <= 0 or P > 1000:
        raise ValueError(f"Invalid period: {P}")
    if np.isnan(dur) or dur <= 0 or dur > P:
        raise ValueError(f"Invalid duration: {dur} (P={P})")
    if len(t) < 10:
        raise ValueError(f"Insufficient data points: {len(t)}")
    if np.all(np.isnan(flux)) or np.nanstd(flux) == 0:
        raise ValueError("Invalid flux data")
    return True


def ingress_egress_area_per_transit(t, resid, P, T0, dur):
    """Compute ingress and egress areas for each transit."""
    nmin = int(np.floor((t.min() - T0) / P))
    nmax = int(np.ceil((t.max() - T0) / P))
    transits = T0 + np.arange(nmin, nmax+1) * P
    Tn_list, area_ing_list, area_egr_list = [], [], []
    
    for tm in transits:
        ing_mask = (t >= tm - 1.5*dur) & (t < tm - 0.5*dur)
        egr_mask = (t >  tm + 0.5*dur) & (t <= tm + 1.5*dur)
        if ing_mask.sum() < 1 or egr_mask.sum() < 1:
            continue
        dt_ing = np.median(np.diff(t[ing_mask])) if ing_mask.sum() > 1 else np.nanmedian(np.diff(t))
        dt_egr = np.median(np.diff(t[egr_mask])) if egr_mask.sum() > 1 else np.nanmedian(np.diff(t))
        area_ing = np.nansum(-resid[ing_mask]) * dt_ing
        area_egr = np.nansum(-resid[egr_mask]) * dt_egr
        Tn_list.append(tm)
        area_ing_list.append(area_ing)
        area_egr_list.append(area_egr)
    
    if not Tn_list:
        raise ValueError("No transits with sufficient ingress/egress data")
    return np.array(Tn_list), np.array(area_ing_list), np.array(area_egr_list)


def compute_skew_series(Tn, area_ing, area_egr):
    """Compute skew time series."""
    denom = area_ing + area_egr
    valid = (denom != 0) & np.isfinite(denom) & np.isfinite(area_ing) & np.isfinite(area_egr)
    if valid.sum() < 3:
        raise ValueError("Insufficient points for skew series")
    return Tn[valid], (area_ing[valid] - area_egr[valid]) / denom[valid]


def detect_skew_period(Tn, S_n):
    """Detect periodic signal in skew time series."""
    if len(Tn) < 3:
        return np.nan, np.nan, np.nan
    
    time_span = Tn.max() - Tn.min()
    if time_span <= 0:
        return np.nan, np.nan, np.nan
    
    ls = LombScargle(Tn, S_n)
    freqs, power = ls.autopower(
        minimum_frequency=1/time_span,
        maximum_frequency=0.5
    )
    if len(power) == 0:
        return np.nan, np.nan, np.nan
    
    idx = np.nanargmax(power)
    best_period = 1/freqs[idx]
    best_power = float(power[idx])
    
    # Estimate FAP
    N_indep = len(Tn)
    fap = N_indep * np.exp(-best_power)
    
    return best_period, best_power, fap


def extract_transit_times(t, flux, flux_err, P, T0, dur, window_factor=2.0):
    """Extract individual transit mid-times."""
    n_min = int(np.floor((t.min() - T0) / P))
    n_max = int(np.ceil((t.max() - T0) / P))
    
    transit_times = []
    transit_errors = []
    
    for n in range(n_min, n_max + 1):
        expected_time = T0 + n * P
        window = window_factor * dur
        
        mask = np.abs(t - expected_time) < window
        if mask.sum() < 5:
            continue
            
        t_transit = t[mask]
        flux_transit = flux[mask]
        err_transit = flux_err[mask] if flux_err is not None else np.ones_like(flux_transit) * np.nanstd(flux[~mask])
        
        try:
            bls = BoxLeastSquares(t_transit, flux_transit)
            periods = np.array([P])
            durations = np.array([dur])
            
            result = bls.power(periods, durations)
            if len(result.transit_time) > 0:
                transit_time = float(result.transit_time[0])
                model = bls.model(t_transit, P, dur, transit_time)
                residuals = flux_transit - model
                time_uncertainty = dur / (2 * np.sqrt(np.sum(1 / err_transit**2)))
                
                transit_times.append(transit_time)
                transit_errors.append(time_uncertainty)
        except:
            continue
    
    return np.array(transit_times), np.array(transit_errors)


def fit_linear_ephemeris(transit_times, transit_errors=None):
    """Fit linear ephemeris to transit times."""
    if len(transit_times) < 2:
        raise ValueError("Need at least 2 transits")
    
    if transit_errors is None:
        transit_errors = np.ones_like(transit_times) * np.nanstd(np.diff(transit_times))
    
    n = np.arange(len(transit_times))
    weights = 1.0 / transit_errors**2
    
    sum_w = np.sum(weights)
    sum_wn = np.sum(weights * n)
    sum_wt = np.sum(weights * transit_times)
    sum_wn2 = np.sum(weights * n**2)
    sum_wnt = np.sum(weights * n * transit_times)
    
    det = sum_w * sum_wn2 - sum_wn**2
    if det == 0:
        raise ValueError("Singular matrix in ephemeris fit")
    
    P = (sum_w * sum_wnt - sum_wn * sum_wt) / det
    T0 = (sum_wn2 * sum_wt - sum_wn * sum_wnt) / det
    
    residuals = transit_times - (T0 + n * P)
    chi2_val = np.sum(weights * residuals**2)
    dof = len(transit_times) - 2
    sigma = np.sqrt(chi2_val / dof) if dof > 0 else 1.0
    
    P_err = sigma * np.sqrt(sum_w / det)
    T0_err = sigma * np.sqrt(sum_wn2 / det)
    
    return P, T0, P_err, T0_err


def detect_ttv_period(transit_times, oc_residuals, oc_errors=None, 
                     min_period=None, max_period=None):
    """Search for periodic TTV signal."""
    if len(transit_times) < 3:
        return np.nan, np.nan, np.nan, np.nan, np.nan
    
    time_span = transit_times.max() - transit_times.min()
    if min_period is None:
        min_period = 0.1 * np.median(np.diff(transit_times))
    if max_period is None:
        max_period = min(10 * np.median(np.diff(transit_times)), time_span / 2)
    
    if min_period >= max_period:
        return np.nan, np.nan, np.nan, np.nan, np.nan
    
    try:
        if oc_errors is not None:
            ls = LombScargle(transit_times, oc_residuals, dy=oc_errors)
        else:
            ls = LombScargle(transit_times, oc_residuals)
        
        freqs = np.linspace(1/max_period, 1/min_period, 10000)
        power = ls.power(freqs)
        periods = 1.0 / freqs
        
        idx_max = np.argmax(power)
        best_period = periods[idx_max]
        best_power = float(power[idx_max])
        
        N = len(transit_times)
        fap = N * np.exp(-best_power)
        
        return periods, power, best_period, best_power, fap
    except:
        return np.nan, np.nan, np.nan, np.nan, np.nan


def compute_shoulder_snr(t, resid, P, T0, dur):
    """Compute shoulder signal-to-noise ratios."""
    phase = ((t - T0 + 0.5*P) % P) / P - 0.5
    window = dur / P
    
    lead_mask = (phase >= -1.5*window) & (phase < -0.5*window)
    trail_mask = (phase > 0.5*window) & (phase <= 1.5*window)
    out_mask = np.abs(phase) > 2*window
    
    if out_mask.sum() < 10:
        return np.nan, np.nan
    
    sigma = np.nanstd(resid[out_mask])
    if sigma == 0 or np.isnan(sigma) or not np.isfinite(sigma):
        return np.nan, np.nan
    
    # Check if we have enough points in each region
    if lead_mask.sum() < 3:
        S_lead = np.nan
    else:
        lead_mean = np.nanmean(resid[lead_mask])
        if np.isfinite(lead_mean):
            S_lead = float(lead_mean / sigma)
        else:
            S_lead = np.nan
    
    if trail_mask.sum() < 3:
        S_trail = np.nan
    else:
        trail_mean = np.nanmean(resid[trail_mask])
        if np.isfinite(trail_mean):
            S_trail = float(trail_mean / sigma)
        else:
            S_trail = np.nan
    
    return S_lead, S_trail


def compute_variability_metrics(t, resid, P, T0, dur):
    """Compute variability in ingress, egress, and baseline regions."""
    nmin = int(np.floor((t.min() - T0) / P))
    nmax = int(np.ceil((t.max() - T0) / P))
    transits = T0 + np.arange(nmin, nmax+1) * P
    
    sig_ing, sig_egr, sig_base = [], [], []
    
    for tm in transits:
        mask = (t >= tm - 2*dur) & (t <= tm + 2*dur)
        if mask.sum() < 10:
            continue
        
        t_loc, r_loc = t[mask], resid[mask]
        
        # Remove any NaN or infinite values
        valid = np.isfinite(r_loc)
        if valid.sum() < 5:
            continue
        t_loc = t_loc[valid]
        r_loc = r_loc[valid]
        
        ing_mask = (t_loc >= tm - 1.5*dur) & (t_loc < tm - 0.5*dur)
        egr_mask = (t_loc > tm + 0.5*dur) & (t_loc <= tm + 1.5*dur)
        base_mask = np.abs(t_loc - tm) > 2*dur
        
        if ing_mask.sum() >= 5:
            sig_val = np.nanstd(r_loc[ing_mask])
            if np.isfinite(sig_val) and sig_val > 0:
                sig_ing.append(sig_val)
        if egr_mask.sum() >= 5:
            sig_val = np.nanstd(r_loc[egr_mask])
            if np.isfinite(sig_val) and sig_val > 0:
                sig_egr.append(sig_val)
        if base_mask.sum() >= 5:
            sig_val = np.nanstd(r_loc[base_mask])
            if np.isfinite(sig_val) and sig_val > 0:
                sig_base.append(sig_val)
    
    if not (sig_ing and sig_egr and sig_base):
        return np.nan, np.nan, np.nan, np.nan
    
    sigma_ing = np.nanmedian(sig_ing)
    sigma_egr = np.nanmedian(sig_egr)
    sigma_base = np.nanmedian(sig_base)
    
    if not np.isfinite(sigma_base) or sigma_base == 0:
        V_ratio = np.nan
    else:
        max_sig = max(sigma_ing if np.isfinite(sigma_ing) else 0,
                     sigma_egr if np.isfinite(sigma_egr) else 0)
        V_ratio = max_sig / sigma_base if max_sig > 0 else np.nan
    
    return sigma_ing, sigma_egr, sigma_base, V_ratio


def lookup_planet_parameters(planet_name):
    """
    Look up planet parameters from NASA Exoplanet Archive.
    
    Parameters
    ----------
    planet_name : str
        Planet name (e.g., "WASP-10 b")
        
    Returns
    -------
    params : dict
        Dictionary with planet parameters, or None if not found
    """
    try:
        # Query NASA Exoplanet Archive
        tbl = NasaExoplanetArchive.query_criteria(
            table="ps",
            select="pl_name,pl_orbper,pl_rade,hostname,tic_id",
            where=f"pl_name='{planet_name}'"
        )
        
        if len(tbl) == 0:
            # Try searching by hostname if planet name format doesn't match
            star_name = planet_name.split()[0]  # e.g., "WASP-10" from "WASP-10 b"
            tbl = NasaExoplanetArchive.query_criteria(
                table="ps",
                select="pl_name,pl_orbper,pl_rade,hostname,tic_id",
                where=f"hostname LIKE '%{star_name}%'"
            )
        
        if len(tbl) > 0:
            df = tbl.to_pandas()
            # Take first result (or best match)
            row = df.iloc[0]
            
            params = {
                'planet_name': str(row.get('pl_name', planet_name)),
                'P': float(row.get('pl_orbper', np.nan)) if pd.notna(row.get('pl_orbper')) else np.nan,
                'Rp': float(row.get('pl_rade', np.nan)) if pd.notna(row.get('pl_rade')) else np.nan,
                'hostname': str(row.get('hostname', '')),
                'tic_id': str(row.get('tic_id', '')) if pd.notna(row.get('tic_id')) else None
            }
            
            return params
        else:
            return None
    except Exception as e:
        print(f"  ⚠ Could not query exoplanet archive: {e}")
        return None


def process_single_planet(planet_name, tic_id=None, P_expected=None, Rp=None, lookup_params=True):
    """
    Process a single planet with all detection methods.
    
    Parameters
    ----------
    planet_name : str
        Planet name (e.g., "WASP-10 b")
    tic_id : str or int, optional
        TIC ID if known
    P_expected : float, optional
        Expected orbital period in days
    Rp : float, optional
        Planet radius in Earth radii
        
    Returns
    -------
    result : dict
        Comprehensive results dictionary
    """
    result = {
        'planet_name': planet_name,
        'tic_id': tic_id if tic_id else 'Unknown',
        'P_expected': P_expected if P_expected else np.nan,
        'Rp': Rp if Rp else np.nan,
        'Status': 'Error',
        'Reason': '',
        'num_transits_observed': 0,
        'data_quality_score': 0.0,
    }
    
    # Initialize all method results
    for prefix in ['skew', 'ttv', 'shoulder', 'variability']:
        for key in ['P_refined', 'T0', 'dur', 'power', 'fap', 'amplitude', 'snr']:
            result[f'{prefix}_{key}'] = np.nan
    
    try:
        print(f"\n{'='*60}")
        print(f"Processing: {planet_name}")
        print(f"{'='*60}")
        
        # Look up planet parameters from database if not provided
        if lookup_params and (P_expected is None or tic_id is None or Rp is None):
            print("Looking up planet parameters from NASA Exoplanet Archive...")
            params = lookup_planet_parameters(planet_name)
            
            if params:
                if tic_id is None and params.get('tic_id'):
                    tic_id = params['tic_id']
                    print(f"  Found TIC ID: {tic_id}")
                if P_expected is None and not np.isnan(params.get('P', np.nan)):
                    P_expected = params['P']
                    print(f"  Found orbital period: {P_expected:.4f} days")
                if Rp is None and not np.isnan(params.get('Rp', np.nan)):
                    Rp = params['Rp']
                    print(f"  Found planet radius: {Rp:.2f} Earth radii")
                if params.get('hostname'):
                    print(f"  Host star: {params['hostname']}")
            else:
                print("  ⚠ Planet not found in database - proceeding without prior parameters")
        
        # Download light curve
        if tic_id:
            # Use TIC ID if provided
            tic_str = str(tic_id).strip()
            tic_num_str = tic_str.replace('TIC', '').strip()
            try:
                tic_num = int(tic_num_str)
                search_str = f"TIC {tic_num}"
            except ValueError:
                search_str = tic_str
            print(f"Searching using TIC ID: {search_str}")
        else:
            # Use planet name
            search_str = planet_name.split()[0]  # Use star name (e.g., "WASP-10" from "WASP-10 b")
            print(f"Searching using star name: {search_str}")
        
        search_result = search_lightcurve(search_str, mission='TESS', author='SPOC', 
                                         cadence='short')
        
        if len(search_result) == 0:
            raise ValueError(f'No TESS lightcurve found for {search_str}')
        
        print(f"Found {len(search_result)} TESS sector(s)")
        lcf = search_result.download_all()
        
        if not lcf or len(lcf) == 0:
            raise ValueError(f'No TESS lightcurve files downloaded for {search_str}')
        
        # Process light curve
        print("Processing light curve...")
        lc = (LightCurveCollection(lcf)
              .stitch()
              .flatten(window_length=401, break_tolerance=0.5)
              .normalize())
        
        t = lc.time.value
        flux = lc.flux.value
        flux_err = lc.flux_err.value if hasattr(lc, 'flux_err') else None
        
        print(f"  Data points: {len(t)}")
        print(f"  Time span: {t.max() - t.min():.1f} days")
        
        if len(t) < 50:
            raise ValueError(f'Too few data points: {len(t)}')
        
        # Refine transit parameters with BLS
        print("Searching for transits...")
        bls = BoxLeastSquares(t, flux)
        
        # Strategy: If we have expected period, search around it with appropriate duration
        # Otherwise, do a wide search with adaptive duration ranges
        if P_expected and not np.isnan(P_expected) and P_expected > 0:
            print(f"  Expected period: {P_expected:.4f} d")
            
            # Primary search: Around expected period with appropriate duration
            periods_narrow = np.linspace(0.8*P_expected, 1.2*P_expected, 5000)
            # Duration for this period range: 0.01 to 0.2 of period
            durations_narrow = np.linspace(0.01*P_expected, 0.2*P_expected, 50)
            print(f"  Primary search: P={0.8*P_expected:.2f}-{1.2*P_expected:.2f} d, dur={0.01*P_expected:.3f}-{0.2*P_expected:.3f} d")
            
            bls_power_narrow = bls.power(periods_narrow, durations_narrow)
            best_power_narrow = np.nanmax(bls_power_narrow.power)
            idx_narrow = np.nanargmax(bls_power_narrow.power)
            i_period_n, i_dur_n = np.unravel_index(idx_narrow, (len(periods_narrow), len(durations_narrow)))
            P_narrow = float(periods_narrow[i_period_n])
            power_narrow = float(bls_power_narrow.power[idx_narrow])
            
            # Secondary search: Wide range for comparison (to catch aliases)
            periods_wide = np.linspace(0.5, 50, 10000)
            # Use shorter durations for wide search (catches short-period planets)
            durations_wide = np.linspace(0.005, 0.1, 50)
            print(f"  Secondary search: P=0.5-50 d, dur=0.005-0.1 d (for comparison)")
            
            bls_power_wide = bls.power(periods_wide, durations_wide)
            best_power_wide = np.nanmax(bls_power_wide.power)
            idx_wide = np.nanargmax(bls_power_wide.power)
            i_period_w, i_dur_w = np.unravel_index(idx_wide, (len(periods_wide), len(durations_wide)))
            P_wide = float(periods_wide[i_period_w])
            power_wide = float(bls_power_wide.power[idx_wide])
            
            # Choose result: Prefer narrow search if power is comparable (within 20%)
            # This prioritizes the expected period unless wide search is much better
            if power_narrow >= 0.8 * power_wide:
                # Use narrow search result
                P_ref = P_narrow
                dur = float(durations_narrow[i_dur_n])
                T0 = float(bls_power_narrow.transit_time[idx_narrow])
                power_used = power_narrow
                print(f"  ✓ Using primary search result (power={power_narrow:.2f} vs {power_wide:.2f})")
            else:
                # Wide search is significantly better - use it but warn
                P_ref = P_wide
                dur = float(durations_wide[i_dur_w])
                T0 = float(bls_power_wide.transit_time[idx_wide])
                power_used = power_wide
                period_diff = abs(P_wide - P_expected) / P_expected
                print(f"  ⚠ Using secondary search result (power={power_wide:.2f} vs {power_narrow:.2f})")
                print(f"  ⚠ Found period ({P_wide:.4f} d) differs from expected by {period_diff*100:.1f}%")
                print(f"  ⚠ This may indicate a period alias or different signal")
        else:
            # No expected period - do wide search
            periods = np.linspace(0.5, 50, 10000)
            durations = np.linspace(0.005, 0.1, 50)
            print(f"  Searching period range: 0.5 - 50 days (no prior period)")
            print(f"  Duration range: 0.005 - 0.1 days")
            
            bls_power = bls.power(periods, durations)
            idx = np.nanargmax(bls_power.power)
            i_period, i_dur = np.unravel_index(idx, (len(periods), len(durations)))
            
            P_ref = float(periods[i_period])
            dur = float(durations[i_dur])
            T0 = float(bls_power.transit_time[idx])
            power_used = float(bls_power.power[idx])
        
        print(f"  Found transit: P={P_ref:.4f} d, T0={T0:.2f}, dur={dur:.4f} d, power={power_used:.2f}")
        
        # Validate parameters
        validate_parameters(P_ref, dur, t, flux)
        
        # Create transit model and compute residuals
        model = bls.model(t, P_ref, dur, T0)
        resid = flux - model
        
        # Compute data quality metrics
        num_transits = int((t.max() - t.min()) / P_ref)
        snr_per_transit = np.nanstd(model) / np.nanstd(resid) if np.nanstd(resid) > 0 else 0
        data_completeness = 1.0 - np.sum(np.isnan(flux)) / len(flux)
        
        result['num_transits_observed'] = num_transits
        result['data_quality_score'] = min(snr_per_transit * data_completeness / 10.0, 1.0)
        
        print(f"  Transits observed: {num_transits}")
        print(f"  Data quality score: {result['data_quality_score']:.3f}")
        
        # ====================================================================
        # 1. SKEW DETECTION
        # ====================================================================
        print("\n[1/4] Running Skew Detection...")
        try:
            Tn, Ai, Ae = ingress_egress_area_per_transit(t, resid, P_ref, T0, dur)
            Tn_ok, S_n = compute_skew_series(Tn, Ai, Ae)
            P_skew, skew_power, skew_fap = detect_skew_period(Tn_ok, S_n)
            
            result['skew_P_refined'] = P_ref
            result['skew_T0'] = T0
            result['skew_dur'] = dur
            result['skew_P_skew'] = P_skew
            result['skew_power'] = skew_power
            result['skew_fap'] = skew_fap
            result['skew_amplitude'] = float(np.nanstd(S_n))
            result['skew_num_transits'] = len(Tn_ok)
            
            print(f"  ✓ Skew period: {P_skew:.4f} d (power={skew_power:.2f}, FAP={skew_fap:.4e})")
        except Exception as e:
            print(f"  ✗ Skew detection failed: {e}")
            result['skew_Status'] = f'Error: {str(e)}'
        
        # ====================================================================
        # 2. TTV DETECTION
        # ====================================================================
        print("\n[2/4] Running TTV Detection...")
        try:
            transit_times, transit_errors = extract_transit_times(
                t, flux, flux_err, P_ref, T0, dur
            )
            
            if len(transit_times) >= 3:
                P_linear, T0_linear, P_err, T0_err = fit_linear_ephemeris(
                    transit_times, transit_errors
                )
                
                # Compute O-C residuals
                n = np.round((transit_times - T0_linear) / P_linear).astype(int)
                expected_times = T0_linear + n * P_linear
                oc_residuals = transit_times - expected_times
                oc_errors = np.sqrt(transit_errors**2 + T0_err**2)
                
                # Search for periodic TTV
                _, _, P_ttv, ttv_power, ttv_fap = detect_ttv_period(
                    transit_times, oc_residuals, oc_errors,
                    min_period=0.1*P_ref,
                    max_period=10*P_ref
                )
                
                ttv_amplitude = float(np.nanstd(oc_residuals))  # in days
                ttv_amplitude_minutes = ttv_amplitude * 24 * 60  # convert to minutes
                
                result['ttv_P_refined'] = P_ref
                result['ttv_T0'] = T0_linear
                result['ttv_dur'] = dur
                result['ttv_P_ttv'] = P_ttv
                result['ttv_power'] = ttv_power
                result['ttv_fap'] = ttv_fap
                result['ttv_amplitude'] = ttv_amplitude
                result['ttv_amplitude_minutes'] = ttv_amplitude_minutes
                result['ttv_num_transits'] = len(transit_times)
                
                print(f"  ✓ TTV period: {P_ttv:.4f} d (power={ttv_power:.2f}, FAP={ttv_fap:.4e})")
                if ttv_amplitude_minutes < 1000:
                    print(f"  ✓ TTV amplitude: {ttv_amplitude_minutes:.2f} minutes ({ttv_amplitude*24:.3f} hours)")
                else:
                    print(f"  ⚠ TTV amplitude: {ttv_amplitude_minutes:.2f} minutes ({ttv_amplitude:.3f} days) - very large, may indicate poor fit")
            else:
                print(f"  ✗ Insufficient transits for TTV: {len(transit_times)} (need ≥3)")
                result['ttv_Status'] = f'Insufficient transits: {len(transit_times)}'
        except Exception as e:
            print(f"  ✗ TTV detection failed: {e}")
            result['ttv_Status'] = f'Error: {str(e)}'
        
        # ====================================================================
        # 3. SHOULDER DETECTION
        # ====================================================================
        print("\n[3/4] Running Shoulder Detection...")
        try:
            S_lead, S_trail = compute_shoulder_snr(t, resid, P_ref, T0, dur)
            
            result['shoulder_P_refined'] = P_ref
            result['shoulder_T0'] = T0
            result['shoulder_dur'] = dur
            result['shoulder_S_lead'] = S_lead
            result['shoulder_S_trail'] = S_trail
            result['shoulder_snr'] = max(abs(S_lead) if not np.isnan(S_lead) else 0,
                                        abs(S_trail) if not np.isnan(S_trail) else 0)
            
            if np.isfinite(S_lead):
                print(f"  ✓ Lead shoulder SNR: {S_lead:.2f}")
            else:
                print(f"  ⚠ Lead shoulder SNR: NaN (insufficient data)")
            if np.isfinite(S_trail):
                print(f"  ✓ Trail shoulder SNR: {S_trail:.2f}")
            else:
                print(f"  ⚠ Trail shoulder SNR: NaN (insufficient data)")
        except Exception as e:
            print(f"  ✗ Shoulder detection failed: {e}")
            result['shoulder_Status'] = f'Error: {str(e)}'
        
        # ====================================================================
        # 4. VARIABILITY DETECTION
        # ====================================================================
        print("\n[4/4] Running Variability Detection...")
        try:
            sigma_ing, sigma_egr, sigma_base, V_ratio = compute_variability_metrics(
                t, resid, P_ref, T0, dur
            )
            
            result['variability_P_refined'] = P_ref
            result['variability_T0'] = T0
            result['variability_dur'] = dur
            result['variability_sigma_ing'] = sigma_ing
            result['variability_sigma_egr'] = sigma_egr
            result['variability_sigma_base'] = sigma_base
            result['variability_V_ratio'] = V_ratio
            
            if np.isfinite(V_ratio):
                print(f"  ✓ Variability ratio: {V_ratio:.3f}")
                print(f"    (σ_ing={sigma_ing:.6f}, σ_egr={sigma_egr:.6f}, σ_base={sigma_base:.6f})")
            else:
                print(f"  ⚠ Variability ratio: NaN (insufficient data in transit regions)")
        except Exception as e:
            print(f"  ✗ Variability detection failed: {e}")
            result['variability_Status'] = f'Error: {str(e)}'
        
        # ====================================================================
        # COMPUTE COMBINED SCORE
        # ====================================================================
        score = 0.0
        
        # Skew contribution (0-0.3)
        if not np.isnan(result.get('skew_power', np.nan)):
            skew_score = min(result['skew_power'] / 10.0, 1.0)
            if not np.isnan(result.get('skew_fap', np.nan)) and result['skew_fap'] > 0:
                fap_penalty = min(-np.log10(result['skew_fap'] + 1e-10) / 3.0, 1.0)
                skew_score *= fap_penalty
            score += 0.25 * skew_score
        
        # TTV contribution (0-0.4) - highest weight
        if not np.isnan(result.get('ttv_power', np.nan)):
            ttv_score = min(result['ttv_power'] / 10.0, 1.0)
            if not np.isnan(result.get('ttv_fap', np.nan)) and result['ttv_fap'] > 0:
                fap_penalty = min(-np.log10(result['ttv_fap'] + 1e-10) / 3.0, 1.0)
                ttv_score *= fap_penalty
            score += 0.40 * ttv_score
        
        # Shoulder contribution (0-0.15)
        if not np.isnan(result.get('shoulder_snr', np.nan)):
            shoulder_score = min(result['shoulder_snr'] / 5.0, 1.0)
            score += 0.15 * shoulder_score
        
        # Variability contribution (0-0.15)
        if not np.isnan(result.get('variability_V_ratio', np.nan)):
            var_score = min(result['variability_V_ratio'] / 3.0, 1.0)
            score += 0.15 * var_score
        
        # Data quality bonus (0-0.05)
        score += 0.05 * result['data_quality_score']
        
        # Clamp score to [0, 1]
        result['combined_exomoon_score'] = max(0.0, min(score, 1.0))
        result['Status'] = 'OK'
        
        print(f"\n{'='*60}")
        print(f"RESULTS SUMMARY")
        print(f"{'='*60}")
        print(f"Status: {result['Status']}")
        print(f"Combined Exomoon Score: {result['combined_exomoon_score']:.3f} / 1.0")
        print(f"  - Skew: {result.get('skew_power', np.nan):.2f} (FAP={result.get('skew_fap', np.nan):.4e})")
        print(f"  - TTV: {result.get('ttv_power', np.nan):.2f} (FAP={result.get('ttv_fap', np.nan):.4e})")
        print(f"  - Shoulder SNR: {result.get('shoulder_snr', np.nan):.2f}")
        print(f"  - Variability Ratio: {result.get('variability_V_ratio', np.nan):.3f}")
        print(f"{'='*60}\n")
        
    except Exception as e:
        result['Status'] = 'Error'
        result['Reason'] = str(e)[:200]
        print(f"\n✗ ERROR: {result['Reason']}\n")
    
    return result


def main():
    parser = argparse.ArgumentParser(
        description='Test exomoon detection on a single planet',
        formatter_class=argparse.RawDescriptionHelpFormatter,
        epilog="""
Examples:
  python test_single_planet.py "WASP-10 b"
  python test_single_planet.py --name "WASP-10" --period 3.09
  python test_single_planet.py --tic 123456789
  python test_single_planet.py "WASP-10 b" --period 3.09 --radius 3.15
        """
    )
    parser.add_argument('planet', nargs='?', type=str,
                       help='Planet name (e.g., "WASP-10 b")')
    parser.add_argument('--name', '--planet', type=str,
                       help='Planet or star name')
    parser.add_argument('--tic', type=str,
                       help='TIC ID (e.g., "123456789" or "TIC 123456789")')
    parser.add_argument('--period', '--P', type=float,
                       help='Orbital period in days')
    parser.add_argument('--radius', '--Rp', type=float,
                       help='Planet radius in Earth radii')
    parser.add_argument('--no-lookup', action='store_true',
                       help='Disable automatic lookup from exoplanet database')
    parser.add_argument('--output', '-o', type=str,
                       help='Output CSV file (optional)')
    
    args = parser.parse_args()
    
    # Determine planet name
    planet_name = args.planet or args.name
    if not planet_name and not args.tic:
        parser.print_help()
        print("\nError: Must provide either planet name or TIC ID")
        sys.exit(1)
    
    if not planet_name:
        planet_name = f"TIC {args.tic}"
    
    # Process the planet
    result = process_single_planet(
        planet_name=planet_name,
        tic_id=args.tic,
        P_expected=args.period,
        Rp=args.radius,
        lookup_params=not args.no_lookup
    )
    
    # Save to CSV if requested
    if args.output:
        df_result = pd.DataFrame([result])
        df_result.to_csv(args.output, index=False)
        print(f"Results saved to {args.output}")
    
    # Exit with error code if failed
    if result['Status'] == 'Error':
        sys.exit(1)
    else:
        sys.exit(0)


if __name__ == '__main__':
    main()