Asteroseismology

class Seismology:
    """
    Main asteroseismology analysis class for stellar oscillations.
    Provides comprehensive tools for seismic analysis and parameter estimation.
    """
    
    def __init__(self, periodogram):
        """
        Initialize seismology analysis.
        
        Parameters:
        - periodogram: Periodogram - Power spectrum of stellar oscillations
        """
    
    def estimate_numax(self, **kwargs):
        """
        Estimate frequency of maximum oscillation power (numax).
        
        Parameters:
        - method: str - Estimation method ('acf2d', 'gaussian_fit')
        - window_width: float - Frequency window for analysis
        
        Returns:
        SeismologyQuantity with numax estimate and uncertainty
        """
    
    def estimate_deltanu(self, **kwargs):
        """
        Estimate large frequency separation (deltanu).
        
        Parameters:
        - method: str - Estimation method ('acf2d', 'comb_fit')
        - numax: float - Frequency of maximum power (auto-estimated if None)
        
        Returns:
        SeismologyQuantity with deltanu estimate and uncertainty
        """
    
    def diagnose_numax(self, **kwargs):
        """Generate diagnostic plots for numax estimation"""
    
    def diagnose_deltanu(self, **kwargs):
        """Generate diagnostic plots for deltanu estimation"""

class SeismologyQuantity:
    """
    Specialized Astropy Quantity for seismology measurements.
    Stores value, uncertainty, and metadata for seismic parameters.
    """
    
    def __init__(self, value, uncertainty=None, **kwargs):
        """
        Parameters:
        - value: float - Parameter value  
        - uncertainty: float - Parameter uncertainty
        - method: str - Estimation method used
        - diagnostics: dict - Diagnostic information
        """

def estimate_numax_acf2d(periodogram, **kwargs):
    """
    Estimate frequency of maximum oscillation power using 2D autocorrelation.
    
    Parameters:
    - periodogram: Periodogram - Power spectrum 
    - frequency_step: float - Frequency resolution for analysis
    - window_width: float - Analysis window width
    
    Returns:
    SeismologyQuantity with numax estimate
    """

def diagnose_numax_acf2d(periodogram, numax=None, **kwargs):
    """
    Generate diagnostic plots for numax estimation.
    
    Parameters:
    - periodogram: Periodogram - Power spectrum
    - numax: float - Numax estimate to highlight
    """

def estimate_deltanu_acf2d(periodogram, numax=None, **kwargs):
    """
    Estimate large frequency separation using 2D autocorrelation.
    
    Parameters:
    - periodogram: Periodogram - Power spectrum
    - numax: float - Frequency of maximum power
    - deltanu_min: float - Minimum deltanu to search
    - deltanu_max: float - Maximum deltanu to search
    
    Returns:
    SeismologyQuantity with deltanu estimate
    """

def diagnose_deltanu_acf2d(periodogram, deltanu=None, numax=None, **kwargs):
    """
    Generate diagnostic plots for deltanu estimation.
    
    Parameters:
    - periodogram: Periodogram - Power spectrum
    - deltanu: float - Deltanu estimate to highlight  
    - numax: float - Numax estimate for context
    """

def estimate_radius(numax, deltanu, teff):
    """
    Estimate stellar radius from seismic scaling relations.
    
    Parameters:
    - numax: float or SeismologyQuantity - Frequency of maximum power (μHz)
    - deltanu: float or SeismologyQuantity - Large frequency separation (μHz)  
    - teff: float - Effective temperature (K)
    
    Returns:
    SeismologyQuantity with radius estimate in solar radii
    """

def estimate_mass(numax, deltanu, teff):
    """
    Estimate stellar mass from seismic scaling relations.
    
    Parameters:
    - numax: float or SeismologyQuantity - Frequency of maximum power (μHz)
    - deltanu: float or SeismologyQuantity - Large frequency separation (μHz)
    - teff: float - Effective temperature (K)
    
    Returns:  
    SeismologyQuantity with mass estimate in solar masses
    """

def estimate_logg(numax, deltanu, teff):
    """
    Estimate surface gravity from seismic scaling relations.
    
    Parameters:
    - numax: float or SeismologyQuantity - Frequency of maximum power (μHz)
    - deltanu: float or SeismologyQuantity - Large frequency separation (μHz)
    - teff: float - Effective temperature (K)
    
    Returns:
    SeismologyQuantity with log(g) estimate
    """

import lightkurve as lk
import numpy as np

# Download high-quality Kepler data for asteroseismology
search = lk.search_lightcurve('KIC 11904151', mission='Kepler')
lc_collection = search.download_all()

# Stitch quarters together for long baseline
lc = lc_collection.stitch()

# Prepare for asteroseismology (remove trends, normalize)
lc = lc.normalize().remove_outliers()

# Create high-resolution periodogram for seismology
pg = lc.to_periodogram(method='lombscargle',
                      minimum_frequency=1,   # μHz
                      maximum_frequency=300, # μHz  
                      samples_per_peak=10)

# Initialize seismology analysis
seismo = lk.Seismology(pg)

# Estimate seismic parameters
numax = seismo.estimate_numax()
deltanu = seismo.estimate_deltanu()

print(f"Numax: {numax.value:.1f} ± {numax.uncertainty:.1f} μHz")
print(f"Delta nu: {deltanu.value:.2f} ± {deltanu.uncertainty:.2f} μHz")

# Generate diagnostic plots for parameter estimates
seismo.diagnose_numax()
seismo.diagnose_deltanu()

# Or use standalone functions
lk.diagnose_numax_acf2d(pg, numax=numax.value)
lk.diagnose_deltanu_acf2d(pg, deltanu=deltanu.value, numax=numax.value)

# Stellar parameters from seismology (requires effective temperature)
teff = 5777  # Sun-like temperature in Kelvin

# Estimate fundamental stellar properties
radius = lk.estimate_radius(numax, deltanu, teff)
mass = lk.estimate_mass(numax, deltanu, teff)  
logg = lk.estimate_logg(numax, deltanu, teff)

print(f"Stellar radius: {radius.value:.2f} ± {radius.uncertainty:.2f} R_sun")
print(f"Stellar mass: {mass.value:.2f} ± {mass.uncertainty:.2f} M_sun")
print(f"Surface gravity: {logg.value:.2f} ± {logg.uncertainty:.2f} dex")

# Custom frequency analysis for detailed mode identification
# Focus on p-mode region around numax
freq_min = numax.value - 5 * deltanu.value
freq_max = numax.value + 5 * deltunu.value

pg_focused = lc.to_periodogram(
    minimum_frequency=freq_min,
    maximum_frequency=freq_max,
    samples_per_peak=20  # High resolution
)

# Look for individual oscillation modes
# Solar-like oscillations have regular frequency spacing
mode_frequencies = []
for n in range(-10, 11):  # Radial order range
    freq = numax.value + n * deltanu.value
    if freq_min <= freq <= freq_max:
        mode_frequencies.append(freq)

# Plot periodogram with predicted mode frequencies
import matplotlib.pyplot as plt
pg_focused.plot()
for freq in mode_frequencies:
    plt.axvline(freq, color='red', alpha=0.5, linestyle='--')
plt.xlabel('Frequency (μHz)')
plt.ylabel('Power')
plt.title('Predicted Solar-like Oscillation Modes')

# Analyze multiple stars for population study
targets = ['KIC 11904151', 'KIC 9139151', 'KIC 10068307']
seismic_results = []

for target in targets:
    try:
        # Download and process
        search = lk.search_lightcurve(target, mission='Kepler')
        lc = search.download_all().stitch().normalize()
        
        # Create periodogram
        pg = lc.to_periodogram(minimum_frequency=1, maximum_frequency=300)
        
        # Seismic analysis
        seismo = lk.Seismology(pg)
        numax = seismo.estimate_numax()
        deltanu = seismo.estimate_deltanu()
        
        # Store results
        seismic_results.append({
            'target': target,
            'numax': numax.value,
            'numax_err': numax.uncertainty,
            'deltanu': deltanu.value,
            'deltanu_err': deltanu.uncertainty
        })
        
        print(f"{target}: numax={numax.value:.1f}±{numax.uncertainty:.1f} μHz, "
              f"Δν={deltanu.value:.2f}±{deltanu.uncertainty:.2f} μHz")
              
    except Exception as e:
        print(f"Failed to analyze {target}: {e}")

# Convert to table for analysis
import pandas as pd
results_df = pd.DataFrame(seismic_results)
print(results_df)

# Use custom scaling relations or literature values
def custom_radius_estimate(numax, deltanu, teff, numax_sun=3100, deltanu_sun=135):
    """Custom radius scaling relation with different solar reference values"""
    radius_ratio = (numax / numax_sun) * (deltanu_sun / deltanu) * np.sqrt(teff / 5777)
    return radius_ratio  # in solar units

# Apply custom scaling
custom_radius = custom_radius_estimate(numax.value, deltanu.value, teff)
print(f"Custom radius estimate: {custom_radius:.2f} R_sun")

# Assess quality of seismic detection
def assess_seismic_quality(seismo_result):
    """Assess quality of seismic parameter estimation"""
    quality_metrics = {}
    
    # Signal-to-noise ratio
    if hasattr(seismo_result, 'snr'):
        quality_metrics['snr'] = seismo_result.snr
    
    # Relative uncertainty
    rel_uncertainty = seismo_result.uncertainty / seismo_result.value
    quality_metrics['relative_uncertainty'] = rel_uncertainty
    
    # Quality flag based on uncertainty
    if rel_uncertainty < 0.05:
        quality_metrics['quality'] = 'Excellent'
    elif rel_uncertainty < 0.10:
        quality_metrics['quality'] = 'Good'
    elif rel_uncertainty < 0.20:
        quality_metrics['quality'] = 'Fair'
    else:
        quality_metrics['quality'] = 'Poor'
    
    return quality_metrics

# Assess parameter quality
numax_quality = assess_seismic_quality(numax)
deltanu_quality = assess_seismic_quality(deltanu)

print(f"Numax quality: {numax_quality['quality']} "
      f"(σ/μ = {numax_quality['relative_uncertainty']:.1%})")
print(f"Deltanu quality: {deltanu_quality['quality']} "
      f"(σ/μ = {deltanu_quality['relative_uncertainty']:.1%})")