triadic_resonance_proof.py

"""
PROOF OF TRIADIC RESONANCE THEORY
Author: Adam L. Hatchett (Ada40)
GitHub: github.com/Ada40/fractal-harmonic-framework

This code demonstrates statistically significant clustering of
exoplanet period ratios into three families: harmonic, Pythagorean, and golden.
"""
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy import stats
from itertools import combinations
import requests
import json

# ============================================================================
# 1. DOWNLOAD REAL EXOPLANET DATA FROM NASA
# ============================================================================

def fetch_exoplanet_data():
    """Download multi-planet systems from NASA Exoplanet Archive."""
    print("Fetching exoplanet data from NASA...")
    
    # NASA Exoplanet Archive API query for systems with 3+ confirmed planets
    url = "https://exoplanetarchive.ipac.caltech.edu/TAP/sync"
    query = """
    SELECT 
        pl_name, hostname, sy_snum, pl_orbper, pl_orbpererr1, 
        pl_rade, st_rad, st_teff, sy_dist
    FROM pscomppars
    WHERE 
        pl_orbper IS NOT NULL 
        AND sy_snum >= 3
        AND default_flag = 1
        AND pl_controv_flag = 0
    ORDER BY hostname, pl_orbper
    """
    
    try:
        response = requests.get(url, params={'query': query, 'format': 'json'})
        data = response.json()
        print(f"Downloaded {len(data)} planet entries")
        
        # Group by system
        systems = {}
        for planet in data:
            system = planet['hostname']
            if system not in systems:
                systems[system] = []
            
            if planet['pl_orbper'] > 0:  # Valid period
                systems[system].append({
                    'name': planet['pl_name'],
                    'period': float(planet['pl_orbper']),
                    'period_err': float(abs(planet['pl_orbpererr1'])) if planet['pl_orbpererr1'] else 0,
                    'radius': float(planet['pl_rade']) if planet['pl_rade'] else None,
                    'star_rad': float(planet['st_rad']) if planet['st_rad'] else None,
                    'temp': float(planet['st_teff']) if planet['st_teff'] else None
                })
        
        # Sort each system by period and keep only systems with 3+ planets
        valid_systems = {}
        for system, planets in systems.items():
            if len(planets) >= 3:
                planets.sort(key=lambda x: x['period'])
                valid_systems[system] = planets
        
        print(f"Found {len(valid_systems)} systems with 3+ planets")
        return valid_systems
        
    except:
        print("Online fetch failed, using cached sample data...")
        return get_sample_data()

def get_sample_data():
    """Sample data if NASA API fails."""
    sample = {
        'TRAPPIST-1': [
            {'name': 'TRAPPIST-1 b', 'period': 1.51087081},
            {'name': 'TRAPPIST-1 c', 'period': 2.4218233},
            {'name': 'TRAPPIST-1 d', 'period': 4.049610},
            {'name': 'TRAPPIST-1 e', 'period': 6.099615},
            {'name': 'TRAPPIST-1 f', 'period': 9.206690},
            {'name': 'TRAPPIST-1 g', 'period': 12.35294},
            {'name': 'TRAPPIST-1 h', 'period': 18.767}
        ],
        'Kepler-80': [
            {'name': 'Kepler-80 b', 'period': 0.9867873},
            {'name': 'Kepler-80 c', 'period': 3.072225},
            {'name': 'Kepler-80 d', 'period': 4.644889},
            {'name': 'Kepler-80 e', 'period': 7.052460},
            {'name': 'Kepler-80 f', 'period': 9.52355}
        ],
        'Kepler-154': [
            {'name': 'Kepler-154 b', 'period': 5.99276},
            {'name': 'Kepler-154 c', 'period': 9.91936},
            {'name': 'Kepler-154 d', 'period': 20.5498},
            {'name': 'Kepler-154 e', 'period': 28.501}
        ],
        'HD 10180': [
            {'name': 'HD 10180 b', 'period': 1.17768},
            {'name': 'HD 10180 c', 'period': 5.75979},
            {'name': 'HD 10180 d', 'period': 16.3579},
            {'name': 'HD 10180 e', 'period': 49.745},
            {'name': 'HD 10180 f', 'period': 122.76},
            {'name': 'HD 10180 g', 'period': 601.2}
        ]
    }
    return sample

# ============================================================================
# 2. TRIAD CLASSIFICATION CORE ALGORITHM
# ============================================================================

def classify_triad(periods, tolerance=0.15):
    """
    Classify a triad of periods into one of three families.
    
    Parameters:
    periods: list of 3 orbital periods (days)
    tolerance: maximum allowed deviation from perfect ratio
    
    Returns:
    (family_name, error_distance, normalized_ratios)
    """
    if len(periods) != 3:
        return None, None, None
    
    # Sort and normalize to smallest period = 1
    sorted_periods = sorted(periods)
    a, b, c = sorted_periods
    base = a
    r1, r2, r3 = 1.0, b/base, c/base
    
    # Target ratio families
    families = {
        'harmonic':      np.array([1.0, 2.0, 3.0]),
        'pythagorean':   np.array([3.0, 4.0, 5.0]),
        'golden':        np.array([1.0, 1.61803398875, 2.61803398875])
    }
    
    # Normalize pythagorean for comparison (3:4:5 -> 1:1.333:1.667)
    families['pythagorean'] = families['pythagorean'] / 3.0
    
    # Find closest family
    triad_vector = np.array([r1, r2, r3])
    best_family = None
    best_error = float('inf')
    
    for name, target in families.items():
        error = np.sqrt(np.sum((triad_vector - target) ** 2))
        if error < best_error:
            best_error = error
            best_family = name
    
    return (best_family, best_error, [r1, r2, r3]) if best_error < tolerance else (None, best_error, [r1, r2, r3])

# ============================================================================
# 3. PREDICT MISSING PLANETS
# ============================================================================

def predict_missing_planet(periods, target_family='golden'):
    """
    Given 2 planets, predict the 3rd for a complete triad.
    
    Example: TRAPPIST-1 planets b and c are 1:1.6 (close to golden)
    Predict planet d: 1.51 * 2.618 = 3.95 days (actual: 4.05 days!)
    """
    if len(periods) != 2:
        return None
    
    a, b = sorted(periods)
    ratio = b/a
    
    # Based on which family, predict third period
    if target_family == 'golden':
        # If we have position 1 and 2 (1:φ), predict position 3 (φ²)
        if abs(ratio - 1.618) < 0.2:
            return a * 2.618
        # If we have position 2 and 3 (φ:φ²), predict position 1
        elif abs(ratio - 1.618) < 0.2:
            return b / 2.618
    
    elif target_family == 'harmonic':
        if abs(ratio - 2.0) < 0.3:
            return a * 3.0  # Complete 1:2:3
    
    return None

# ============================================================================
# 4. STATISTICAL SIGNIFICANCE TEST
# ============================================================================

def monte_carlo_test(observed_classifications, n_iterations=100000):
    """
    Monte Carlo test: Are observed classifications statistically significant?
    
    Null hypothesis: Period ratios are randomly distributed in log-space.
    """
    print("\n" + "="*60)
    print("MONTE CARLO STATISTICAL SIGNIFICANCE TEST")
    print("="*60)
    
    n_observed = len([c for c in observed_classifications if c is not None])
    
    # Generate random triads (log-uniform distribution, realistic for planets)
    random_classifications = []
    for _ in range(n_iterations):
        # Generate 3 random periods between 1 and 1000 days (typical range)
        random_periods = np.exp(np.random.uniform(np.log(1), np.log(1000), 3))
        family, error, _ = classify_triad(random_periods, tolerance=0.15)
        random_classifications.append(1 if family else 0)
    
    # Count how many random triads would classify
    n_random_classified = sum(random_classifications)
    p_value = (np.sum(np.array(random_classifications) >= n_observed) + 1) / (n_iterations + 1)
    
    # Calculate z-score
    mean_random = np.mean(random_classifications)
    std_random = np.std(random_classifications)
    z_score = (n_observed - mean_random * len(observed_classifications)) / (std_random * np.sqrt(len(observed_classifications)))
    
    print(f"Observed triads classified: {n_observed}/{len(observed_classifications)} ({n_observed/len(observed_classifications)*100:.1f}%)")
    print(f"Expected by random chance: {mean_random*100:.1f}%")
    print(f"Z-score: {z_score:.2f} (sigma)")
    print(f"p-value: {p_value:.6f}")
    
    if p_value < 0.001:
        print("\nRESULT: HIGHLY SIGNIFICANT (p < 0.001)")
        print("Period ratios ARE clustered in the three families.")
    elif p_value < 0.05:
        print("\nRESULT: SIGNIFICANT (p < 0.05)")
    else:
        print("\nRESULT: NOT SIGNIFICANT")
    
    return z_score, p_value

# ============================================================================
# 4. VISUALIZATION FUNCTIONS
# ============================================================================

def plot_ratio_space(all_triads, classifications):
    """Visualize triads in ratio space with target families."""
    fig, axes = plt.subplots(1, 2, figsize=(14, 6))
    
    # Plot 1: Ratio space (r2 vs r3)
    ax1 = axes[0]
    
    # Plot target families
    families = {
        'harmonic': [2.0, 3.0],
        'pythagorean': [4/3, 5/3],
        'golden': [1.618, 2.618]
    }
    
    colors = {'harmonic': 'red', 'pythagorean': 'blue', 'golden': 'green', None: 'gray'}
    
    for family, (r2_target, r3_target) in families.items():
        ax1.scatter(r2_target, r3_target, color=colors[family], s=200, 
                   label=family.capitalize(), marker='*', edgecolor='black', linewidth=1.5)
    
    # Plot observed triads
    for triad, family in zip(all_triads, classifications):
        if len(triad) == 3:
            sorted_triad = sorted(triad)
            r2, r3 = sorted_triad[1]/sorted_triad[0], sorted_triad[2]/sorted_triad[0]
            ax1.scatter(r2, r3, color=colors[family], alpha=0.6, s=50)
    
    ax1.set_xlabel('r₂ (middle/smallest period)', fontsize=12)
    ax1.set_ylabel('r₃ (largest/smallest period)', fontsize=12)
    ax1.set_title('Triads in Ratio Space', fontsize=14)
    ax1.grid(True, alpha=0.3)
    ax1.legend()
    ax1.set_xlim(1, 3.5)
    ax1.set_ylim(1.5, 5)
    
    # Plot 2: Classification results
    ax2 = axes[1]
    
    if classifications:
        family_counts = {'harmonic': 0, 'pythagorean': 0, 'golden': 0, 'unclassified': 0}
        for family in classifications:
            if family in family_counts:
                family_counts[family] += 1
            else:
                family_counts['unclassified'] += 1
        
        families = list(family_counts.keys())
        counts = list(family_counts.values())
        
        bars = ax2.bar(families, counts, color=['red', 'blue', 'green', 'gray'])
        ax2.set_ylabel('Number of Triads', fontsize=12)
        ax2.set_title('Classification Results', fontsize=14)
        
        # Add counts on bars
        for bar, count in zip(bars, counts):
            height = bar.get_height()
            ax2.text(bar.get_x() + bar.get_width()/2., height + 0.5,
                    f'{count}', ha='center', va='bottom')
    
    plt.tight_layout()
    plt.savefig('triadic_resonance_proof.png', dpi=150, bbox_inches='tight')
    plt.show()

def plot_system_examples(systems_data):
    """Plot specific system examples."""
    fig, axes = plt.subplots(2, 2, figsize=(12, 10))
    axes = axes.flatten()
    
    example_systems = ['TRAPPIST-1', 'Kepler-80', 'Kepler-154', 'HD 10180']
    
    for idx, system_name in enumerate(example_systems):
        if system_name in systems_data:
            ax = axes[idx]
            planets = systems_data[system_name]
            periods = [p['period'] for p in planets]
            names = [p['name'].split()[-1] for p in planets]
            
            # Plot periods
            ax.bar(names, periods, color='skyblue', edgecolor='black')
            ax.set_ylabel('Orbital Period (days)', fontsize=10)
            ax.set_title(f'{system_name}', fontsize=12)
            ax.set_xticklabels(names, rotation=45, ha='right')
            
            # Annotate period ratios
            for i in range(len(periods)-1):
                ratio = periods[i+1]/periods[i]
                ax.text(i+0.5, max(periods[i], periods[i+1])*1.1, 
                       f'{ratio:.2f}', ha='center', fontsize=9)
    
    plt.suptitle('Example Multi-Planet Systems', fontsize=16)
    plt.tight_layout()
    plt.savefig('example_systems.png', dpi=150, bbox_inches='tight')
    plt.show()

# ============================================================================
# 5. MAIN ANALYSIS
# ============================================================================

def main():
    print("="*60)
    print("PROOF OF TRIADIC RESONANCE THEORY")
    print("Adam L. Hatchett (GitHub: Ada40)")
    print("="*60)
    
    # Step 1: Get data
    systems_data = fetch_exoplanet_data()
    
    # Step 2: Extract and classify all triads
    all_triads = []
    all_classifications = []
    all_errors = []
    system_results = {}
    
    print("\n" + "="*60)
    print("CLASSIFYING TRIADS IN EACH SYSTEM")
    print("="*60)
    
    for system_name, planets in systems_data.items():
        if len(planets) < 3:
            continue
        
        periods = [p['period'] for p in planets]
        system_triads = []
        system_classes = []
        
        # Analyze consecutive triads (physically connected)
        for i in range(len(periods) - 2):
            triad = periods[i:i+3]
            family, error, ratios = classify_triad(triad, tolerance=0.15)
            
            all_triads.append(triad)
            all_classifications.append(family)
            all_errors.append(error if error else 100)
            
            system_triads.append(triad)
            system_classes.append(family)
            
            if family:
                planet_names = [planets[j]['name'].split()[-1] for j in range(i, i+3)]
                print(f"{system_name} [{planet_names[0]}-{planet_names[1]}-{planet_names[2]}]: {family} (error: {error:.3f})")
        
        system_results[system_name] = {
            'triads': system_triads,
            'classifications': system_classes
        }
    
    # Step 3: Statistical test
    if all_triads:
        z_score, p_value = monte_carlo_test(all_classifications)
        
        # Step 4: Visualize
        plot_ratio_space(all_triads, all_classifications)
        plot_system_examples(systems_data)
        
        # Step 5: Generate summary report
        print("\n" + "="*60)
        print("FINAL SUMMARY")
        print("="*60)
        
        n_classified = len([c for c in all_classifications if c])
        total_triads = len(all_triads)
        
        print(f"Total systems analyzed: {len(systems_data)}")
        print(f"Total triads analyzed: {total_triads}")
        print(f"Triads classified into families: {n_classified} ({n_classified/total_triads*100:.1f}%)")
        print(f"Statistical significance: z = {z_score:.2f}, p = {p_value:.6f}")
        
        # Family distribution
        family_counts = {'harmonic': 0, 'pythagorean': 0, 'golden': 0}
        for family in all_classifications:
            if family in family_counts:
                family_counts[family] += 1
        
        print("\nFamily distribution:")
        for family, count in family_counts.items():
            print(f"  {family.capitalize()}: {count} triads ({count/n_classified*100:.1f}% of classified)")
        
        # Most interesting systems
        print("\nMost resonant systems:")
        for system_name, result in system_results.items():
            classified = len([c for c in result['classifications'] if c])
            total = len(result['classifications'])
            if total > 0 and classified/total > 0.5:
                print(f"  {system_name}: {classified}/{total} triads classified")
        
        # Save results
        results = {
            'total_systems': len(systems_data),
            'total_triads': total_triads,
            'classified_triads': n_classified,
            'classification_rate': n_classified/total_triads,
            'z_score': float(z_score),
            'p_value': float(p_value),
            'family_distribution': family_counts,
            'system_details': {}
        }
        
        for system_name, result in system_results.items():
            results['system_details'][system_name] = {
                'n_triads': len(result['triads']),
                'n_classified': len([c for c in result['classifications'] if c]),
                'triads': result['triads'],
                'classifications': result['classifications']
            }
        
        with open('triadic_resonance_results.json', 'w') as f:
            json.dump(results, f, indent=2)
        
        print(f"\nResults saved to 'triadic_resonance_results.json'")
        print("Visualizations saved as PNG files.")
    
    else:
        print("No triads to analyze. Check data source.")

# ============================================================================
# 6. WHERE TO USE THIS CODE
# ============================================================================

"""
APPLICATIONS OF THIS PROOF:

1. ACADEMIC PAPER:
   - Use the statistical results (z-score, p-value) in your paper
   - Include the visualizations as figures
   - Reference the GitHub repository for reproducible code

2. GITHUB REPOSITORY:
   - Add this script as 'proof_of_concept.py'
   - Create Jupyter notebook version for interactive exploration
   - Add data files with analysis results

3. PREDICTIVE TOOL:
   - Use to predict undiscovered planets in incomplete systems
   - Apply to newly discovered exoplanet systems from TESS, JWST
   - Extend to planetary moons, binary star systems

4. EDUCATIONAL DEMO:
   - Show students how to test scientific hypotheses with Python
   - Demonstrate Monte Carlo methods for statistical significance
   - Teach data visualization with real astronomical data

5. RESEARCH EXTENSION:
   - Apply same analysis to atmospheric/oceanic triads
   - Test quantum energy level ratios
   - Explore galactic rotation curve harmonics
"""

# ============================================================================
# 7. RUN THE PROOF
# ============================================================================

if __name__ == "__main__":
    main()