exospective.py

#!/usr/bin/env python3
"""
exospective

Generates star charts showing the night sky from arbitrary positions in 3D space.
"""

import numpy as np
import pandas as pd
from skyfield.api import Star, load
from skyfield.data import hipparcos
from skyfield.positionlib import Barycentric
import matplotlib.pyplot as plt
from typing import Tuple, Optional, List
import argparse
import json
import os
import sys
from astroquery.simbad import Simbad
# Manual collision detection (adjustText removed for v1)


# Star name cache file
STAR_NAME_CACHE_FILE = '.star_name_cache.json'

# Load or initialize star name cache
def load_star_name_cache():
    """Load cached star names from disk."""
    if os.path.exists(STAR_NAME_CACHE_FILE):
        try:
            with open(STAR_NAME_CACHE_FILE, 'r') as f:
                return json.load(f)
        except (json.JSONDecodeError, IOError):
            return {}
    return {}

def save_star_name_cache(cache):
    """Save star name cache to disk."""
    try:
        with open(STAR_NAME_CACHE_FILE, 'w') as f:
            json.dump(cache, f, indent=2)
    except Exception as e:
        print(f"Warning: Could not save star name cache: {e}", file=sys.stderr)

def query_star_info_from_simbad(hip_id):
    """
    Query SIMBAD for the best name, spectral type, and temperature for a star.

    Parameters:
    -----------
    hip_id : int
        Hipparcos catalog number

    Returns:
    --------
    tuple : (name, spectral_type, teff) - spectral_type and teff may be None
    """
    try:
        custom_simbad = Simbad()
        custom_simbad.add_votable_fields('ids', 'sp_type', 'mesfe_h')

        result = custom_simbad.query_object(f"HIP {hip_id}")

        if result is None or len(result) == 0:
            return f"HIP {hip_id}", None, None

        # Get spectral type and temperature
        sp_type = result['sp_type'][0] if result['sp_type'][0] else None
        teff = None
        if 'mesfe_h.teff' in result.colnames:
            teff_val = result['mesfe_h.teff'][0]
            if teff_val and not np.ma.is_masked(teff_val):
                teff = float(teff_val)

        # Get all identifiers
        ids_str = result['ids'][0]
        if isinstance(ids_str, bytes):
            ids_str = ids_str.decode('utf-8')

        identifiers = [id.strip() for id in ids_str.split('|')]

        # Name priority hierarchy
        # 1. Check for NAME prefix first (proper names like "NAME Sirius")
        for id in identifiers:
            if id.startswith('NAME '):
                return id.replace('NAME ', ''), sp_type, teff

        # 2. Look for simple proper names (single word, capitalized, not a catalog designation)
        for id in identifiers:
            if id.startswith(('HIP ', 'HD ', 'HR ', 'SAO ', 'TYC ', 'BD', 'CD', 'CPD', 'GJ ', '2MASS', 'GCRV', 'GEN#', 'GSC', 'AG', 'PPM')) or '*' in id:
                continue
            if len(id.split()) == 1 and id[0].isupper():
                return id, sp_type, teff

        # 3. Look for Bayer designations (Greek letters)
        greek_letters = ['alf', 'bet', 'gam', 'del', 'eps', 'zet', 'eta', 'tet', 'iot', 'kap', 'lam', 'mu', 'nu', 'xi', 'omi', 'pi', 'rho', 'sig', 'tau', 'ups', 'phi', 'chi', 'psi', 'ome']
        for id in identifiers:
            for greek in greek_letters:
                if greek in id.lower() and len(id.split()) <= 3:
                    return id.replace('* ', '').strip(), sp_type, teff

        # 4. Look for Flamsteed designations (number + constellation)
        for id in identifiers:
            parts = id.split()
            if len(parts) == 2 and parts[0].isdigit():
                return id.replace('* ', '').strip(), sp_type, teff

        # 5. Look for HD catalog
        for id in identifiers:
            if id.startswith('HD '):
                return id, sp_type, teff

        # 6. Fallback to HIP
        return f"HIP {hip_id}", sp_type, teff

    except Exception as e:
        print(f"Warning: Could not query SIMBAD for HIP {hip_id}: {e}", file=sys.stderr)
        return f"HIP {hip_id}", None, None

def get_star_info(hip_id, cache):
    """
    Get star name, spectral type, and temperature using cache or SIMBAD.

    Parameters:
    -----------
    hip_id : int
        Hipparcos catalog number
    cache : dict
        Star info cache

    Returns:
    --------
    tuple : (name, spectral_type, teff)
    """
    hip_str = str(hip_id)

    if hip_str in cache:
        cached = cache[hip_str]
        # Handle old cache formats
        if isinstance(cached, str):
            return cached, None, None
        return cached.get('name', f"HIP {hip_id}"), cached.get('sp_type'), cached.get('teff')

    # Query SIMBAD
    name, sp_type, teff = query_star_info_from_simbad(hip_id)

    # Cache successful queries
    if not name.startswith('HIP '):
        cache[hip_str] = {'name': name, 'sp_type': sp_type, 'teff': teff}

    return name, sp_type, teff


# Common star names mapped to Hipparcos IDs
# Includes brightest stars and well-known named stars
def teff_to_rgb(teff):
    """
    Convert effective temperature to RGB color.
    Uses blackbody approximation for star colors.
    
    Parameters:
    -----------
    teff : float
        Effective temperature in Kelvin
    
    Returns:
    --------
    tuple : (r, g, b) values 0-1
    """
    # Approximate RGB from temperature (simplified blackbody)
    # Based on typical star color temperatures
    t = teff / 100.0
    
    # Red channel
    if t <= 66:
        r = 1.0
    else:
        r = 1.29293618 * ((t - 60) ** -0.1332047592)
        r = max(0, min(1, r))
    
    # Green channel
    if t <= 66:
        g = 0.390081579 * np.log(t) - 0.631841444
    else:
        g = 1.12989086 * ((t - 60) ** -0.0755148492)
    g = max(0, min(1, g))
    
    # Blue channel
    if t >= 66:
        b = 1.0
    elif t <= 19:
        b = 0.0
    else:
        b = 0.543206789 * np.log(t - 10) - 1.19625408
        b = max(0, min(1, b))
    
    return (r, g, b)


def spectral_to_teff(sp_type):
    """
    Estimate effective temperature from spectral type.
    
    Parameters:
    -----------
    sp_type : str
        Spectral type string (e.g., 'G2V', 'K1.5III', 'M3.5V')
    
    Returns:
    --------
    float : Estimated temperature in Kelvin, or None if unparseable
    """
    if not sp_type:
        return None
    
    # Base temperatures for each spectral class (main sequence)
    base_temps = {
        'O': 35000, 'B': 20000, 'A': 9000, 'F': 7000,
        'G': 5500, 'K': 4500, 'M': 3200, 'L': 1800, 'T': 1000
    }
    
    # Temperature range within each class (early to late subtype)
    temp_ranges = {
        'O': 15000, 'B': 10000, 'A': 2000, 'F': 1500,
        'G': 1000, 'K': 1000, 'M': 1400, 'L': 800, 'T': 700
    }
    
    # Parse spectral class (first letter)
    sp_class = sp_type[0].upper()
    if sp_class not in base_temps:
        return None
    
    base = base_temps[sp_class]
    range_t = temp_ranges[sp_class]
    
    # Try to parse subtype (0-9)
    subtype = 5.0  # default to middle
    try:
        # Find digits after the class letter
        i = 1
        num_str = ''
        while i < len(sp_type) and (sp_type[i].isdigit() or sp_type[i] == '.'):
            num_str += sp_type[i]
            i += 1
        if num_str:
            subtype = float(num_str)
            subtype = min(9.9, max(0, subtype))
    except:
        pass
    
    # Calculate temperature (higher subtype = cooler)
    teff = base - (subtype / 10.0) * range_t
    return teff


def get_star_color(sp_type, teff=None):
    """
    Get RGB color for a star based on temperature or spectral type.
    
    Parameters:
    -----------
    sp_type : str
        Spectral type string
    teff : float, optional
        Effective temperature in Kelvin (preferred if available)
    
    Returns:
    --------
    tuple : (r, g, b) values 0-1, or None for white default
    """
    # Use teff if available
    if teff and teff > 0:
        return teff_to_rgb(teff)
    
    # Fall back to spectral type
    estimated_teff = spectral_to_teff(sp_type)
    if estimated_teff:
        return teff_to_rgb(estimated_teff)
    
    return None  # Will use white


STAR_NAMES = {
    # Brightest stars (magnitude < 1.5)
    32349: 'Sirius',           # α CMa, mag -1.46
    30438: 'Canopus',          # α Car, mag -0.72
    71683: 'Rigil Kentaurus',  # α Cen A, mag -0.27
    69673: 'Arcturus',         # α Boo, mag -0.05
    91262: 'Vega',             # α Lyr, mag 0.03
    24608: 'Capella',          # α Aur, mag 0.08
    24436: 'Rigel',            # β Ori, mag 0.13
    37279: 'Procyon',          # α CMi, mag 0.34
    9884: 'Achernar',          # α Eri, mag 0.46
    21421: 'Betelgeuse',       # α Ori, mag 0.50 (variable)
    25336: 'Hadar',            # β Cen, mag 0.61
    80763: 'Altair',           # α Aql, mag 0.77
    68702: 'Acrux',            # α Cru, mag 0.77
    21589: 'Aldebaran',        # α Tau, mag 0.85
    113368: 'Antares',         # α Sco, mag 1.09
    7588: 'Pollux',            # β Gem, mag 1.14
    49669: 'Spica',            # α Vir, mag 1.04
    86032: 'Deneb',            # α Cyg, mag 1.25
    62434: 'Mimosa',           # β Cru, mag 1.30

    # Other well-known stars
    677: 'Alpheratz',          # α And
    3179: 'Schedar',           # α Cas
    5447: 'Mirach',            # β And
    8728: 'Diphda',            # β Cet
    11767: 'Polaris',          # α UMi, North Star
    15863: 'Menkar',           # α Cet
    25428: 'Bellatrix',        # γ Ori
    25930: 'Elnath',           # β Tau
    26311: 'Alnilam',          # ε Ori (Orion's Belt center)
    26727: 'Alnitak',          # ζ Ori (Orion's Belt)
    27989: 'Saiph',            # κ Ori
    30324: 'Miaplacidus',      # β Car
    31681: 'Adhara',           # ε CMa
    33579: 'Castor',           # α Gem
    34444: 'Shaula',           # λ Sco
    36850: 'Alhena',           # γ Gem
    41037: 'Mizar',            # ζ UMa
    44816: 'Dubhe',            # α UMa
    46390: 'Merak',            # β UMa
    49583: 'Alioth',           # ε UMa
    53910: 'Alkaid',           # η UMa
    54061: 'Menkent',          # θ Cen
    57632: 'Alphecca',         # α CrB
    59774: 'Antares',          # α Sco
    61084: 'Kochab',           # β UMi
    65474: 'Rasalhague',       # α Oph
    67301: 'Shaula',           # λ Sco
    72607: 'Eltanin',          # γ Dra
    76267: 'Rasalgethi',       # α Her
    77070: 'Vega',             # α Lyr
    84012: 'Sheliak',          # β Lyr
    85927: 'Albireo',          # β Cyg
    97649: 'Fomalhaut',        # α PsA
    102098: 'Deneb Kaitos',    # β Cet
    109268: 'Peacock',         # α Pav
    113881: 'Alnair',          # α Gru
}

# Constellation center positions (RA in degrees, Dec in degrees)
# Approximate centers for orientation purposes
CONSTELLATION_POSITIONS = {
    # Zodiac constellations
    'Aries': (32.0, 20.0),
    'Taurus': (65.0, 15.0),
    'Gemini': (105.0, 22.0),
    'Cancer': (130.0, 20.0),
    'Leo': (155.0, 15.0),
    'Virgo': (195.0, 0.0),
    'Libra': (230.0, -15.0),
    'Scorpius': (247.0, -25.0),
    'Sagittarius': (285.0, -25.0),
    'Capricornus': (315.0, -18.0),
    'Aquarius': (330.0, -10.0),
    'Pisces': (15.0, 10.0),

    # Major northern constellations
    'Ursa Major': (165.0, 55.0),
    'Ursa Minor': (225.0, 75.0),
    'Cassiopeia': (15.0, 62.0),
    'Cepheus': (315.0, 70.0),
    'Draco': (225.0, 65.0),
    'Cygnus': (305.0, 42.0),
    'Lyra': (280.0, 38.0),
    'Aquila': (295.0, 8.0),
    'Hercules': (258.0, 30.0),
    'Boötes': (210.0, 30.0),
    'Corona Borealis': (233.0, 30.0),
    'Andromeda': (20.0, 37.0),
    'Perseus': (50.0, 43.0),
    'Auriga': (80.0, 42.0),
    'Pegasus': (340.0, 20.0),

    # Orion region
    'Orion': (85.0, 0.0),
    'Canis Major': (105.0, -20.0),
    'Canis Minor': (115.0, 5.0),

    # Southern constellations
    'Centaurus': (200.0, -45.0),
    'Crux': (187.0, -60.0),
    'Carina': (140.0, -64.0),
    'Vela': (135.0, -47.0),
    'Puppis': (120.0, -30.0),
    'Hydra': (170.0, -15.0),
    'Corvus': (185.0, -18.0),
    'Crater': (175.0, -15.0),
    'Lupus': (225.0, -42.0),
    'Ara': (260.0, -55.0),
    'Triangulum Australe': (245.0, -65.0),
    'Pavo': (305.0, -65.0),
    'Grus': (335.0, -47.0),
    'Phoenix': (20.0, -48.0),
    'Eridanus': (55.0, -30.0),
    'Fornax': (45.0, -32.0),
    'Sculptor': (350.0, -32.0),
    'Tucana': (340.0, -65.0),
    'Indus': (315.0, -55.0),
    'Octans': (210.0, -85.0),
    'Columba': (85.0, -35.0),
    'Dorado': (75.0, -60.0),
    'Pictor': (85.0, -55.0),
}


def adjust_label_positions(ax, texts: List, constellation_texts: List = None,
                          text_magnitudes: List = None):
    """
    Manually adjust label positions to avoid overlaps.

    Star labels prefer to be above the star (initial offset already applied).
    Brighter stars (lower magnitude) are processed first and get priority for positions.
    Uses multiple position alternatives and collision detection.

    Parameters:
    -----------
    ax : matplotlib axis
        The axis containing the labels
    texts : list
        List of matplotlib Text objects (star labels and Sol)
    constellation_texts : list, optional
        List of matplotlib Text objects for constellation labels
    text_magnitudes : list, optional
        List of magnitudes corresponding to texts (for priority sorting)
    """
    if not texts and not constellation_texts:
        return

    # Combine texts and magnitudes for sorting
    # Brighter stars (lower magnitude) get processed first
    star_labels = []
    if texts and text_magnitudes:
        # Pair each text with its magnitude for sorting
        for text, mag in zip(texts, text_magnitudes):
            star_labels.append((text, mag))
        # Sort by magnitude (brighter/lower magnitude first)
        star_labels.sort(key=lambda x: x[1])
    elif texts:
        # No magnitude data, just use texts as-is
        star_labels = [(text, 999) for text in texts]

    # Add constellation labels (processed after stars)
    const_labels = []
    if constellation_texts:
        const_labels = [(text, 1000) for text in constellation_texts]

    # Combine: stars (sorted by brightness) then constellations
    all_labels = star_labels + const_labels

    if not all_labels:
        return

    # Get axis limits for calculating label dimensions
    xlim = ax.get_xlim()
    ylim = ax.get_ylim()
    x_range = xlim[1] - xlim[0]
    y_range = ylim[1] - ylim[0]

    # Base label dimensions (as fraction of axis range)
    base_label_width = x_range * 0.015
    base_label_height = y_range * 0.01

    # Store label positions and dimensions
    label_data = []

    # Position alternatives (prefer above, then other positions)
    # Order matters: above is first (preferred for stars)
    position_offsets = [
        (0, 1),      # Above (preferred)
        (1, 0),      # Right
        (0, -1),     # Below
        (-1, 0),     # Left
        (0.7, 0.7),  # Upper right
        (-0.7, 0.7), # Upper left
        (0.7, -0.7), # Lower right
        (-0.7, -0.7),# Lower left
        (0, 1.5),    # Far above
        (1.5, 0),    # Far right
        (0, -1.5),   # Far below
        (-1.5, 0),   # Far left
    ]

    def labels_collide(x1, y1, w1, h1, x2, y2, w2, h2):
        """Check if two labels overlap."""
        return not (x1 + w1/2 < x2 - w2/2 or
                   x1 - w1/2 > x2 + w2/2 or
                   y1 + h1/2 < y2 - h2/2 or
                   y1 - h1/2 > y2 + h2/2)

    # Process each label (in brightness order for stars, then constellations)
    for text, magnitude in all_labels:
        # Get font size and calculate scaled dimensions
        fontsize = text.get_fontsize()
        size_scale = fontsize / 8.0  # Scale relative to base fontsize
        label_width = base_label_width * size_scale
        label_height = base_label_height * size_scale

        # Star position (labels start here)
        star_x = text.get_position()[0]
        star_y = text.get_position()[1]

        # Direction-aware offset multiplier
        # Constellation labels need more space
        offset_multiplier = 1.3 if fontsize <= 7 else 1.0

        # Try each position alternative
        # Always offset labels from stars (start with "above" preference)
        best_position = None
        for dx, dy in position_offsets:
            # Use direction-aware offsets: width for horizontal, height for vertical
            # This prevents labels from being pushed too far in one direction
            offset_x = abs(dx) * label_width * offset_multiplier if dx != 0 else 0
            offset_y = abs(dy) * label_height * offset_multiplier if dy != 0 else 0
            # For diagonal moves, use both components
            test_x = star_x + (dx * label_width * offset_multiplier if dx != 0 else 0)
            test_y = star_y + (dy * label_height * offset_multiplier if dy != 0 else 0)

            # Check collision with all stored labels
            collides = False
            for stored in label_data:
                if labels_collide(test_x, test_y, label_width, label_height,
                                stored['x'], stored['y'], stored['w'], stored['h']):
                    collides = True
                    break

            if not collides:
                best_position = (test_x, test_y)
                break

        # If all positions collide, use the first alternative anyway
        if best_position is None:
            dx, dy = position_offsets[0]
            fallback_x = star_x + (dx * label_width * offset_multiplier if dx != 0 else 0)
            fallback_y = star_y + (dy * label_height * offset_multiplier if dy != 0 else 0)
            best_position = (fallback_x, fallback_y)

        # Update text position
        text.set_position(best_position)

        # Store this label's position and dimensions
        label_data.append({
            'x': best_position[0],
            'y': best_position[1],
            'w': label_width,
            'h': label_height
        })


class StarChartGenerator:
    """Main class for generating star charts from arbitrary observer positions."""

    def __init__(self, ra: float, dec: float, distance: float,
                 system_name: Optional[str] = None,
                 magnitude_cutoff: float = 6.0,
                 label_magnitude_threshold: float = 1.0,
                 show_constellations: bool = False,
                 label_constellations: bool = None,
                 grid_spacing: int = 15,
                 output_dir: str = ".",
                 show_grid: bool = True,
                 show_star_labels: bool = True,
                 show_constellation_labels: bool = True,
                 show_sol_label: bool = True,
                 highlight_sol: bool = True):
        """
        Initialize the star chart generator.

        Parameters:
        -----------
        ra : float
            Right Ascension in decimal degrees (0-360)
        dec : float
            Declination in decimal degrees (-90 to +90)
        distance : float
            Distance from Sol in parsecs
        system_name : str, optional
            Name for the system (default: coordinates string)
        magnitude_cutoff : float
            Faintest magnitude to display (default: 6.0)
        label_magnitude_threshold : float
            Brightest magnitude to label (default: 1.0)
        show_constellations : bool
            Whether to show constellation lines (default: False)
        label_constellations : bool, optional
            Whether to label constellations (auto-decides based on distance if None)
        grid_spacing : int
            Grid line spacing in degrees (default: 15)
        output_dir : str
            Output directory path (default: current directory)
        """
        self.ra = self._validate_ra(ra)
        self.dec = self._validate_dec(dec)
        self.distance = self._validate_distance(distance)
        self.system_name = system_name or f"RA {ra}° Dec {dec}° Distance {distance} pc"
        self.magnitude_cutoff = magnitude_cutoff
        self.label_magnitude_threshold = label_magnitude_threshold
        self.show_constellations = show_constellations

        # Auto-decide constellation labels based on distance (if enabled)
        if label_constellations is None:
            self.label_constellations = show_constellation_labels and (distance * 3.26156 < 40)
        else:
            self.label_constellations = label_constellations and show_constellation_labels

        self.grid_spacing = grid_spacing
        self.output_dir = output_dir
        self.show_grid = show_grid
        self.show_star_labels = show_star_labels
        self.show_sol_label = show_sol_label
        self.highlight_sol = highlight_sol

        # Will be populated during processing
        self.observer_position = None
        self.stars_data = None

        # Load star name cache
        self.star_name_cache = load_star_name_cache()

    @staticmethod
    def _validate_ra(ra: float) -> float:
        """Validate and normalize RA to 0-360 degrees."""
        if ra < 0 or ra > 360:
            raise ValueError(f"RA must be between 0 and 360 degrees, got {ra}")
        return ra

    @staticmethod
    def _validate_dec(dec: float) -> float:
        """Validate declination range."""
        if dec < -90 or dec > 90:
            raise ValueError(f"Dec must be between -90 and +90 degrees, got {dec}")
        return dec

    @staticmethod
    def _validate_distance(distance: float) -> float:
        """Validate distance is positive."""
        if distance <= 0:
            raise ValueError(f"Distance must be positive, got {distance}")
        return distance

    def ra_dec_distance_to_cartesian(self) -> np.ndarray:
        """
        Convert observer's RA/Dec/distance (spherical) to Cartesian coordinates.

        Returns:
        --------
        np.ndarray : [x, y, z] position in AU from Solar System barycenter
        """
        # Convert RA, Dec to radians
        ra_rad = np.radians(self.ra)
        dec_rad = np.radians(self.dec)

        # Convert parsecs to AU (1 parsec = 206265 AU)
        distance_au = self.distance * 206265.0

        # Spherical to Cartesian conversion
        # In equatorial coordinates:
        # x points to RA=0, Dec=0 (vernal equinox)
        # y points to RA=90, Dec=0
        # z points to Dec=90 (north celestial pole)
        x = distance_au * np.cos(dec_rad) * np.cos(ra_rad)
        y = distance_au * np.cos(dec_rad) * np.sin(ra_rad)
        z = distance_au * np.sin(dec_rad)

        return np.array([x, y, z])

    def load_hipparcos_catalog(self):
        """Load the Hipparcos star catalog using Skyfield."""
        print("Loading Hipparcos catalog...", file=sys.stderr)

        # Note: de421.bsp (JPL planetary ephemeris) was intentionally removed here.
        # It would be needed for time-dependent astrometric corrections (stellar
        # aberration, barycentric parallax refinement, date-specific sky rendering),
        # but exospective's current calculations work directly from catalog
        # RA/Dec/parallax and don't require ephemeris data. Re-add if adding
        # time-dependent or high-precision features in the future.

        # Load Hipparcos catalog
        with load.open(hipparcos.URL) as f:
            df = hipparcos.load_dataframe(f)

        print(f"Loaded {len(df)} stars from Hipparcos catalog", file=sys.stderr)
        return df

    def calculate_star_positions_from_observer(self, hipparcos_df: pd.DataFrame) -> pd.DataFrame:
        """
        Calculate positions and apparent magnitudes of stars as seen from observer position.

        Parameters:
        -----------
        hipparcos_df : pd.DataFrame
            Hipparcos catalog data

        Returns:
        --------
        pd.DataFrame : Processed star data with observer-relative positions and magnitudes
        """
        # Get observer position in Cartesian AU
        self.observer_position = self.ra_dec_distance_to_cartesian()

        print(f"Observer position: {self.observer_position / 206265.0} pc", file=sys.stderr)
        print(f"Processing {len(hipparcos_df)} stars...", file=sys.stderr)

        # Filter out stars without parallax data
        stars_with_parallax = hipparcos_df[hipparcos_df['parallax_mas'] > 0].copy()
        print(f"Stars with valid parallax: {len(stars_with_parallax)}", file=sys.stderr)

        # Calculate distance to each star from Earth (in parsecs)
        # parallax_mas is in milliarcseconds, so we use 1000.0 / parallax_mas to get parsecs
        stars_with_parallax['distance_from_earth_pc'] = 1000.0 / stars_with_parallax['parallax_mas']

        # Convert star positions to Cartesian coordinates (in AU)
        # Using the same coordinate system as observer position
        ra_rad = np.radians(stars_with_parallax['ra_degrees'])
        dec_rad = np.radians(stars_with_parallax['dec_degrees'])
        distance_au = stars_with_parallax['distance_from_earth_pc'] * 206265.0

        stars_with_parallax['x_earth'] = distance_au * np.cos(dec_rad) * np.cos(ra_rad)
        stars_with_parallax['y_earth'] = distance_au * np.cos(dec_rad) * np.sin(ra_rad)
        stars_with_parallax['z_earth'] = distance_au * np.sin(dec_rad)

        # Calculate position relative to observer
        stars_with_parallax['x_rel'] = stars_with_parallax['x_earth'] - self.observer_position[0]
        stars_with_parallax['y_rel'] = stars_with_parallax['y_earth'] - self.observer_position[1]
        stars_with_parallax['z_rel'] = stars_with_parallax['z_earth'] - self.observer_position[2]

        # Calculate distance from observer (in AU, then convert to parsecs)
        distance_from_observer_au = np.sqrt(
            stars_with_parallax['x_rel']**2 +
            stars_with_parallax['y_rel']**2 +
            stars_with_parallax['z_rel']**2
        )
        stars_with_parallax['distance_from_observer_pc'] = distance_from_observer_au / 206265.0

        # Recalculate apparent magnitude from observer position
        # Formula: m_new = m_earth + 5 * log10(d_new / d_earth)
        distance_ratio = (stars_with_parallax['distance_from_observer_pc'] /
                         stars_with_parallax['distance_from_earth_pc'])

        stars_with_parallax['magnitude_from_observer'] = (
            stars_with_parallax['magnitude'] + 5 * np.log10(distance_ratio)
        )

        # Filter out stars within 0.5 light-years (to exclude parent star system)
        # 0.5 ly = 0.1533 pc
        min_distance_pc = 0.5 / 3.26156
        stars_with_parallax = stars_with_parallax[
            stars_with_parallax['distance_from_observer_pc'] >= min_distance_pc
        ].copy()

        # Filter by magnitude cutoff
        visible_stars = stars_with_parallax[
            stars_with_parallax['magnitude_from_observer'] <= self.magnitude_cutoff
        ].copy()

        print(f"Stars visible from observer (mag <= {self.magnitude_cutoff}): {len(visible_stars)}", file=sys.stderr)

        # Add star names from STAR_NAMES dictionary
        visible_stars['name'] = visible_stars.index.map(STAR_NAMES)

        return visible_stars

    def add_sol(self, stars_df: pd.DataFrame) -> pd.DataFrame:
        """
        Add Sol (the Sun) to the star data.

        Parameters:
        -----------
        stars_df : pd.DataFrame
            Existing star data

        Returns:
        --------
        pd.DataFrame : Star data with Sol added
        """
        # Sol is at the origin in Earth-centered coordinates
        # Position relative to observer
        sol_x_rel = -self.observer_position[0]
        sol_y_rel = -self.observer_position[1]
        sol_z_rel = -self.observer_position[2]

        # Distance to Sol from observer
        distance_to_sol_au = np.sqrt(sol_x_rel**2 + sol_y_rel**2 + sol_z_rel**2)
        distance_to_sol_pc = distance_to_sol_au / 206265.0

        # Calculate Sol's apparent magnitude from observer
        # Absolute magnitude of Sol: M = 4.83
        # m = M + 5 * log10(d) - 5, where d is in parsecs
        if distance_to_sol_pc > 0:
            sol_magnitude = 4.83 + 5 * np.log10(distance_to_sol_pc) - 5
        else:
            sol_magnitude = -26.74  # As seen from Earth (observer at Sol)

        print(f"Sol distance: {distance_to_sol_pc:.2f} pc, magnitude: {sol_magnitude:.2f}", file=sys.stderr)

        # Create Sol entry
        sol_data = {
            'hip': 0,  # Special ID for Sol
            'magnitude': -26.74,  # From Earth
            'ra_degrees': 0.0,
            'dec_degrees': 0.0,
            'parallax': 1000000.0,  # Dummy value
            'distance_from_earth_pc': 0.0,
            'x_earth': 0.0,
            'y_earth': 0.0,
            'z_earth': 0.0,
            'x_rel': sol_x_rel,
            'y_rel': sol_y_rel,
            'z_rel': sol_z_rel,
            'distance_from_observer_pc': distance_to_sol_pc,
            'magnitude_from_observer': sol_magnitude,
            'name': 'Sol'
        }

        # Add to dataframe
        sol_df = pd.DataFrame([sol_data], index=[0])
        combined_df = pd.concat([stars_df, sol_df])

        return combined_df

    def transform_to_galactic_coordinates(self, stars_df: pd.DataFrame) -> pd.DataFrame:
        """
        Transform star positions from observer to galactic coordinates.

        Parameters:
        -----------
        stars_df : pd.DataFrame
            Star data with relative positions

        Returns:
        --------
        pd.DataFrame : Star data with galactic longitude and latitude
        """
        print("Transforming to galactic coordinates...", file=sys.stderr)

        # Calculate spherical coordinates from observer-relative Cartesian
        x = stars_df['x_rel'].values
        y = stars_df['y_rel'].values
        z = stars_df['z_rel'].values

        # Convert to spherical coordinates (RA, Dec from observer)
        r = np.sqrt(x**2 + y**2 + z**2)
        ra_from_obs = np.degrees(np.arctan2(y, x))
        ra_from_obs = np.where(ra_from_obs < 0, ra_from_obs + 360, ra_from_obs)
        dec_from_obs = np.degrees(np.arcsin(z / r))

        # Now convert from equatorial to galactic coordinates
        # Using standard transformation
        # Galactic north pole: RA = 192.859508°, Dec = 27.128336°
        # Galactic center: RA = 266.405°, Dec = -28.936°

        # Standard rotation matrices for equatorial to galactic conversion
        ra_rad = np.radians(ra_from_obs)
        dec_rad = np.radians(dec_from_obs)

        # Galactic pole in equatorial coordinates
        ra_gp = np.radians(192.859508)
        dec_gp = np.radians(27.128336)
        l_cp = np.radians(122.932)  # Galactic longitude of celestial pole

        # Conversion formulas
        sin_b = (np.sin(dec_rad) * np.sin(dec_gp) +
                np.cos(dec_rad) * np.cos(dec_gp) * np.cos(ra_rad - ra_gp))
        galactic_lat = np.degrees(np.arcsin(sin_b))

        y_gal = np.cos(dec_rad) * np.sin(ra_rad - ra_gp)
        x_gal = (np.sin(dec_rad) * np.cos(dec_gp) -
                np.cos(dec_rad) * np.sin(dec_gp) * np.cos(ra_rad - ra_gp))

        galactic_lon = np.degrees(np.arctan2(y_gal, x_gal)) + np.degrees(l_cp)
        galactic_lon = np.where(galactic_lon < 0, galactic_lon + 360, galactic_lon)
        galactic_lon = np.where(galactic_lon >= 360, galactic_lon - 360, galactic_lon)

        stars_df['galactic_lon'] = galactic_lon
        stars_df['galactic_lat'] = galactic_lat

        return stars_df

    def get_constellation_positions_from_observer(self):
        """
        Calculate where constellation region centers appear from observer's viewpoint.

        Returns:
        --------
        list of tuples: (constellation_name, galactic_lon, galactic_lat)
        """
        constellation_positions = []

        for const_name, (ra_earth, dec_earth) in CONSTELLATION_POSITIONS.items():
            # Convert constellation center position (as seen from Earth) to Cartesian
            # Treat it as a point at a large distance (arbitrary, just for direction)
            distance_au = 1e10  # Arbitrary large distance for direction

            ra_rad = np.radians(ra_earth)
            dec_rad = np.radians(dec_earth)

            # Position in Earth-centered coordinates
            x_earth = distance_au * np.cos(dec_rad) * np.cos(ra_rad)
            y_earth = distance_au * np.cos(dec_rad) * np.sin(ra_rad)
            z_earth = distance_au * np.sin(dec_rad)

            # Position relative to observer
            x_rel = x_earth - self.observer_position[0]
            y_rel = y_earth - self.observer_position[1]
            z_rel = z_earth - self.observer_position[2]

            # Convert to spherical (RA, Dec from observer)
            r = np.sqrt(x_rel**2 + y_rel**2 + z_rel**2)
            ra_from_obs = np.degrees(np.arctan2(y_rel, x_rel))
            if ra_from_obs < 0:
                ra_from_obs += 360
            dec_from_obs = np.degrees(np.arcsin(z_rel / r))

            # Convert to galactic coordinates
            ra_rad = np.radians(ra_from_obs)
            dec_rad = np.radians(dec_from_obs)

            ra_gp = np.radians(192.859508)
            dec_gp = np.radians(27.128336)
            l_cp = np.radians(122.932)

            sin_b = (np.sin(dec_rad) * np.sin(dec_gp) +
                    np.cos(dec_rad) * np.cos(dec_gp) * np.cos(ra_rad - ra_gp))
            galactic_lat = np.degrees(np.arcsin(sin_b))

            y_gal = np.cos(dec_rad) * np.sin(ra_rad - ra_gp)
            x_gal = (np.sin(dec_rad) * np.cos(dec_gp) -
                    np.cos(dec_rad) * np.sin(dec_gp) * np.cos(ra_rad - ra_gp))

            galactic_lon = np.degrees(np.arctan2(y_gal, x_gal)) + np.degrees(l_cp)
            if galactic_lon < 0:
                galactic_lon += 360
            if galactic_lon >= 360:
                galactic_lon -= 360

            constellation_positions.append((const_name, galactic_lon, galactic_lat))

        return constellation_positions

    def get_stars_to_label(self, stars_df, top_n=20):
        """
        Determine which stars should be labeled based on brightness from observer's perspective.
        Queries SIMBAD for proper names.

        Parameters:
        -----------
        stars_df : pd.DataFrame
            Star data with magnitude_from_observer
        top_n : int
            Number of brightest stars to label (default: 20)

        Returns:
        --------
        pd.DataFrame : Stars to label with display names and colors
        """
        # Exclude Sol from this logic (handled separately)
        non_sol_stars = stars_df[stars_df.get('name', pd.Series()) != 'Sol'].copy()

        # Get top N brightest stars from observer
        stars_to_label = non_sol_stars.nsmallest(top_n, 'magnitude_from_observer').copy()

        print(f"Querying SIMBAD for names of {len(stars_to_label)} bright stars...", file=sys.stderr)

        # Query SIMBAD for display names and colors
        display_names = []
        colors = []
        for hip_id in stars_to_label.index:
            name, sp_type, teff = get_star_info(hip_id, self.star_name_cache)
            display_names.append(name)
            color = get_star_color(sp_type, teff)
            colors.append(color if color else (1.0, 1.0, 1.0))

        stars_to_label['display_name'] = display_names
        stars_to_label['color'] = colors

        # Save cache after queries
        save_star_name_cache(self.star_name_cache)

        return stars_to_label

    def calculate_star_size(self, magnitudes):
        """
        Calculate star marker sizes based on apparent magnitude.
        Larger sizes for brighter (smaller magnitude) stars.

        Parameters:
        -----------
        magnitudes : array-like
            Apparent magnitudes

        Returns:
        --------
        array : Marker sizes in points
        """
        # Scale from magnitude to size
        # Magnitude 6 (faintest visible) -> smallest size
        # Magnitude -1.5 (brightest) -> largest size
        # Use exponential scaling for better visual effect

        # Normalize magnitudes to 0-1 range (inverted, so bright = 1)
        mag_range = self.magnitude_cutoff - (-2.0)  # Range from brightest possible to cutoff
        normalized = (self.magnitude_cutoff - magnitudes) / mag_range
        normalized = np.clip(normalized, 0, 1)

        # Exponential scaling for better visual distinction
        # Much smaller sizes with stronger exponential curve for dramatic difference
        size_min = 0.1
        size_max = 20.0
        sizes = size_min + (size_max - size_min) * (normalized ** 2.0)

        return sizes

    def plot_equirectangular(self, stars_df, filename):
        """
        Generate all-sky equirectangular (cylindrical) projection.

        Parameters:
        -----------
        stars_df : pd.DataFrame
            Star data with galactic coordinates
        filename : str
            Output filename
        """
        print(f"Generating equirectangular projection: {filename}", file=sys.stderr)

        fig, ax = plt.subplots(figsize=(16, 8), facecolor='black')
        ax.set_facecolor('black')

        # Get bright stars with colors first (need this for coloring)
        labeled_stars = self.get_stars_to_label(stars_df) if self.show_star_labels else pd.DataFrame()
        bright_hip_ids = set(labeled_stars.index) if len(labeled_stars) > 0 else set()
        bright_colors = {idx: row['color'] for idx, row in labeled_stars.iterrows()} if len(labeled_stars) > 0 else {}

        # Separate bright stars from dim stars for different coloring
        is_bright = stars_df.index.isin(bright_hip_ids)
        dim_stars = stars_df[~is_bright & (stars_df.get('name', pd.Series()) != 'Sol')]
        bright_stars_df = stars_df[is_bright]

        # Plot dim stars in white
        if len(dim_stars) > 0:
            lon = dim_stars['galactic_lon'].values
            lat = dim_stars['galactic_lat'].values
            mag = dim_stars['magnitude_from_observer'].values
            sizes = self.calculate_star_size(mag)
            ax.scatter(lon, lat, s=sizes, c='white', alpha=0.9, edgecolors='none')

        # Plot bright stars with colors
        if len(bright_stars_df) > 0:
            for idx, star in bright_stars_df.iterrows():
                color = bright_colors.get(idx, (1.0, 1.0, 1.0))
                size = self.calculate_star_size(np.array([star['magnitude_from_observer']]))[0]
                ax.scatter(star['galactic_lon'], star['galactic_lat'], 
                          s=size, c=[color], alpha=0.9, edgecolors='none')

        # Plot Sol
        sol_data = stars_df[stars_df.get('name', pd.Series()) == 'Sol']
        if len(sol_data) > 0:
            sol_lon = sol_data['galactic_lon'].values[0]
            sol_lat = sol_data['galactic_lat'].values[0]
            sol_mag = sol_data['magnitude_from_observer'].values[0]
            sol_size = self.calculate_star_size(np.array([sol_mag]))[0]

            # Sol color: yellow if highlighted, G2V star color otherwise
            sol_color = 'yellow' if self.highlight_sol else teff_to_rgb(5778)
            ax.scatter(sol_lon, sol_lat, s=sol_size, c=[sol_color],
                      alpha=1.0, edgecolors='none', zorder=10 if self.highlight_sol else 1)

        # Collect text labels for collision avoidance
        texts = []  # Star labels and Sol
        magnitudes = []  # Corresponding magnitudes for priority sorting
        constellation_texts = []  # Constellation labels

        # Add Sol label (no offset, it's special)
        if len(sol_data) > 0 and self.show_sol_label:
            label_color = 'yellow' if self.highlight_sol else 'white'
            sol_text = ax.text(sol_lon, sol_lat, 'Sol', color=label_color,
                              fontsize=10, ha='center', weight='bold')
            texts.append(sol_text)
            magnitudes.append(sol_mag)

        # Label bright stars (already fetched above for coloring)
        if self.show_star_labels and len(labeled_stars) > 0:
            for idx, star in labeled_stars.iterrows():
                star_text = ax.text(star['galactic_lon'], star['galactic_lat'],
                                   star['display_name'], color='white', fontsize=8,
                                   ha='center', alpha=0.8)
                texts.append(star_text)
                magnitudes.append(star['magnitude_from_observer'])

        # Add constellation region labels if close enough to Sol
        if self.label_constellations:
            const_positions = self.get_constellation_positions_from_observer()
            for const_name, const_lon, const_lat in const_positions:
                const_text = ax.text(const_lon, const_lat, const_name,
                                    color='gray', fontsize=7, alpha=0.5,
                                    ha='center', style='italic')
                constellation_texts.append(const_text)

        # Adjust text positions to avoid overlaps
        # Brighter stars (lower magnitude) get priority for position
        adjust_label_positions(ax, texts, constellation_texts, magnitudes)

        # Add grid
        if self.show_grid:
            ax.grid(True, color='gray', alpha=0.3, linewidth=0.5)
            ax.set_xticks(np.arange(0, 361, self.grid_spacing))
            ax.set_yticks(np.arange(-90, 91, self.grid_spacing))
            ax.set_xlabel('Galactic Longitude (degrees)', color='white', fontsize=12)
            ax.set_ylabel('Galactic Latitude (degrees)', color='white', fontsize=12)
            ax.tick_params(colors='white', labelsize=9)
        else:
            ax.set_xticks([])
            ax.set_yticks([])
            ax.spines['top'].set_visible(False)
            ax.spines['right'].set_visible(False)
            ax.spines['bottom'].set_visible(False)
            ax.spines['left'].set_visible(False)

        # Title
        ax.set_title(f'All-Sky View from {self.system_name}\nEquirectangular Projection',
                    color='white', fontsize=14, pad=20)

        # Set axis limits
        ax.set_xlim(0, 360)
        ax.set_ylim(-90, 90)

        plt.tight_layout()
        plt.savefig(filename, dpi=150, facecolor='black', edgecolor='none')
        plt.close()

        print(f"Saved: {filename}", file=sys.stderr)

    def plot_north_polar(self, stars_df, filename):
        """
        Generate north polar stereographic projection.

        Parameters:
        -----------
        stars_df : pd.DataFrame
            Star data with galactic coordinates
        filename : str
            Output filename
        """
        print(f"Generating north polar projection: {filename}", file=sys.stderr)

        # Filter to northern hemisphere (latitude > -30 for better coverage)
        north_stars = stars_df[stars_df['galactic_lat'] > -30].copy()

        fig, ax = plt.subplots(figsize=(10, 10), facecolor='black')
        ax.set_facecolor('black')

        # Get bright stars with colors
        labeled_stars = self.get_stars_to_label(stars_df) if self.show_star_labels else pd.DataFrame()
        bright_hip_ids = set(labeled_stars.index) if len(labeled_stars) > 0 else set()
        bright_colors = {idx: row['color'] for idx, row in labeled_stars.iterrows()} if len(labeled_stars) > 0 else {}

        # Separate bright from dim stars
        is_bright = north_stars.index.isin(bright_hip_ids)
        dim_stars = north_stars[~is_bright & (north_stars.get('name', pd.Series()) != 'Sol')]
        bright_stars_df = north_stars[is_bright]

        # Plot dim stars in white
        if len(dim_stars) > 0:
            lat = dim_stars['galactic_lat'].values
            lon = dim_stars['galactic_lon'].values
            mag = dim_stars['magnitude_from_observer'].values
            r = 2 * np.tan(np.radians(90 - lat) / 2)
            x = r * np.cos(np.radians(lon))
            y = r * np.sin(np.radians(lon))
            sizes = self.calculate_star_size(mag)
            ax.scatter(x, y, s=sizes, c='white', alpha=0.9, edgecolors='none')

        # Plot bright stars with colors
        for idx, star in bright_stars_df.iterrows():
            color = bright_colors.get(idx, (1.0, 1.0, 1.0))
            star_r = 2 * np.tan(np.radians(90 - star['galactic_lat']) / 2)
            star_x = star_r * np.cos(np.radians(star['galactic_lon']))
            star_y = star_r * np.sin(np.radians(star['galactic_lon']))
            size = self.calculate_star_size(np.array([star['magnitude_from_observer']]))[0]
            ax.scatter(star_x, star_y, s=size, c=[color], alpha=0.9, edgecolors='none')

        # Highlight Sol
        sol_data = north_stars[north_stars.get('name', pd.Series()) == 'Sol']
        if len(sol_data) > 0:
            sol_lat = sol_data['galactic_lat'].values[0]
            sol_lon = sol_data['galactic_lon'].values[0]
            sol_mag = sol_data['magnitude_from_observer'].values[0]

            sol_r = 2 * np.tan(np.radians(90 - sol_lat) / 2)
            sol_x = sol_r * np.cos(np.radians(sol_lon))
            sol_y = sol_r * np.sin(np.radians(sol_lon))
            sol_size = self.calculate_star_size(np.array([sol_mag]))[0]

            sol_color = 'yellow' if self.highlight_sol else teff_to_rgb(5778)
            ax.scatter(sol_x, sol_y, s=sol_size, c=[sol_color],
                      alpha=1.0, edgecolors='none', zorder=10 if self.highlight_sol else 1)

        # Collect text labels for collision avoidance
        texts = []  # Star labels and Sol
        magnitudes = []  # Corresponding magnitudes for priority sorting
        constellation_texts = []  # Constellation labels

        # Add Sol label (no offset, it's special)
        if len(sol_data) > 0 and self.show_sol_label:
            label_color = 'yellow' if self.highlight_sol else 'white'
            sol_text = ax.text(sol_x, sol_y, 'Sol', color=label_color,
                              fontsize=10, ha='center', weight='bold')
            texts.append(sol_text)
            magnitudes.append(sol_mag)

        # Label bright stars (already fetched above for coloring)
        if self.show_star_labels and len(labeled_stars) > 0:
            # Filter to stars in northern hemisphere for this view
            north_labeled = labeled_stars[labeled_stars['galactic_lat'] > -30]
            for idx, star in north_labeled.iterrows():
                star_r = 2 * np.tan(np.radians(90 - star['galactic_lat']) / 2)
                star_x = star_r * np.cos(np.radians(star['galactic_lon']))
                star_y = star_r * np.sin(np.radians(star['galactic_lon']))
                star_text = ax.text(star_x, star_y, star['display_name'], color='white',
                                   fontsize=8, ha='center', alpha=0.8)
                texts.append(star_text)
                magnitudes.append(star['magnitude_from_observer'])

        # Add constellation region labels if close enough to Sol
        if self.label_constellations:
            const_positions = self.get_constellation_positions_from_observer()
            for const_name, const_lon, const_lat in const_positions:
                # Only show constellations in northern hemisphere
                if const_lat > -30:
                    const_lat_rad = np.radians(const_lat)
                    const_lon_rad = np.radians(const_lon)
                    const_r = 2 * np.tan(np.radians(90 - const_lat) / 2)
                    const_x = const_r * np.cos(const_lon_rad)
                    const_y = const_r * np.sin(const_lon_rad)
                    const_text = ax.text(const_x, const_y, const_name,
                           color='gray', fontsize=6, alpha=0.4,
                           ha='center', style='italic')
                    constellation_texts.append(const_text)

        # Adjust text positions to avoid overlaps
        # Brighter stars (lower magnitude) get priority for position
        adjust_label_positions(ax, texts, constellation_texts, magnitudes)

        # Add radial grid lines (latitude circles)
        if self.show_grid:
            for lat_val in range(0, 91, self.grid_spacing):
                r_grid = 2 * np.tan(np.radians(90 - lat_val) / 2)
                circle = plt.Circle((0, 0), r_grid, fill=False, color='gray',
                                  alpha=0.3, linewidth=0.5)
                ax.add_patch(circle)

            # Add longitude lines
            for lon_val in range(0, 360, self.grid_spacing):
                lon_rad = np.radians(lon_val)
                r_max = 2 * np.tan(np.radians(90 + 30) / 2)  # Extend to -30° lat
                ax.plot([0, r_max * np.cos(lon_rad)],
                       [0, r_max * np.sin(lon_rad)],
                       color='gray', alpha=0.3, linewidth=0.5)

        # Title
        ax.set_title(f'North Polar View from {self.system_name}\nStereographic Projection',
                    color='white', fontsize=14, pad=20)

        # Equal aspect ratio
        ax.set_aspect('equal')
        max_r = 2 * np.tan(np.radians(90 + 30) / 2)
        ax.set_xlim(-max_r, max_r)
        ax.set_ylim(-max_r, max_r)

        # Hide axes
        ax.set_xticks([])
        ax.set_yticks([])
        ax.spines['top'].set_visible(False)
        ax.spines['right'].set_visible(False)
        ax.spines['bottom'].set_visible(False)
        ax.spines['left'].set_visible(False)

        plt.tight_layout()
        plt.savefig(filename, dpi=150, facecolor='black', edgecolor='none')
        plt.close()

        print(f"Saved: {filename}", file=sys.stderr)

    def plot_south_polar(self, stars_df, filename):
        """
        Generate south polar stereographic projection.

        Parameters:
        -----------
        stars_df : pd.DataFrame
            Star data with galactic coordinates
        filename : str
            Output filename
        """
        print(f"Generating south polar projection: {filename}", file=sys.stderr)

        # Filter to southern hemisphere (latitude < 30 for better coverage)
        south_stars = stars_df[stars_df['galactic_lat'] < 30].copy()

        fig, ax = plt.subplots(figsize=(10, 10), facecolor='black')
        ax.set_facecolor('black')

        # Get bright stars with colors
        labeled_stars = self.get_stars_to_label(stars_df) if self.show_star_labels else pd.DataFrame()
        bright_hip_ids = set(labeled_stars.index) if len(labeled_stars) > 0 else set()
        bright_colors = {idx: row['color'] for idx, row in labeled_stars.iterrows()} if len(labeled_stars) > 0 else {}

        # Separate bright from dim stars
        is_bright = south_stars.index.isin(bright_hip_ids)
        dim_stars = south_stars[~is_bright & (south_stars.get('name', pd.Series()) != 'Sol')]
        bright_stars_df = south_stars[is_bright]

        # Plot dim stars in white
        if len(dim_stars) > 0:
            lat = dim_stars['galactic_lat'].values
            lon = dim_stars['galactic_lon'].values
            mag = dim_stars['magnitude_from_observer'].values
            r = 2 * np.tan(np.radians(90 + lat) / 2)
            x = r * np.cos(np.radians(lon))
            y = r * np.sin(np.radians(lon))
            sizes = self.calculate_star_size(mag)
            ax.scatter(x, y, s=sizes, c='white', alpha=0.9, edgecolors='none')

        # Plot bright stars with colors
        for idx, star in bright_stars_df.iterrows():
            color = bright_colors.get(idx, (1.0, 1.0, 1.0))
            star_r = 2 * np.tan(np.radians(90 + star['galactic_lat']) / 2)
            star_x = star_r * np.cos(np.radians(star['galactic_lon']))
            star_y = star_r * np.sin(np.radians(star['galactic_lon']))
            size = self.calculate_star_size(np.array([star['magnitude_from_observer']]))[0]
            ax.scatter(star_x, star_y, s=size, c=[color], alpha=0.9, edgecolors='none')

        # Highlight Sol
        sol_data = south_stars[south_stars.get('name', pd.Series()) == 'Sol']
        if len(sol_data) > 0:
            sol_lat = sol_data['galactic_lat'].values[0]
            sol_lon = sol_data['galactic_lon'].values[0]
            sol_mag = sol_data['magnitude_from_observer'].values[0]

            sol_r = 2 * np.tan(np.radians(90 + sol_lat) / 2)
            sol_x = sol_r * np.cos(np.radians(sol_lon))
            sol_y = sol_r * np.sin(np.radians(sol_lon))
            sol_size = self.calculate_star_size(np.array([sol_mag]))[0]

            sol_color = 'yellow' if self.highlight_sol else teff_to_rgb(5778)
            ax.scatter(sol_x, sol_y, s=sol_size, c=[sol_color],
                      alpha=1.0, edgecolors='none', zorder=10 if self.highlight_sol else 1)

        # Collect text labels for collision avoidance
        texts = []  # Star labels and Sol
        magnitudes = []  # Corresponding magnitudes for priority sorting
        constellation_texts = []  # Constellation labels

        # Add Sol label (no offset, it's special)
        if len(sol_data) > 0 and self.show_sol_label:
            label_color = 'yellow' if self.highlight_sol else 'white'
            sol_text = ax.text(sol_x, sol_y, 'Sol', color=label_color,
                              fontsize=10, ha='center', weight='bold')
            texts.append(sol_text)
            magnitudes.append(sol_mag)

        # Label bright stars (already fetched above for coloring)
        if self.show_star_labels and len(labeled_stars) > 0:
            # Filter to stars in southern hemisphere for this view
            south_labeled = labeled_stars[labeled_stars['galactic_lat'] < 30]
            for idx, star in south_labeled.iterrows():
                star_r = 2 * np.tan(np.radians(90 + star['galactic_lat']) / 2)
                star_x = star_r * np.cos(np.radians(star['galactic_lon']))
                star_y = star_r * np.sin(np.radians(star['galactic_lon']))
                star_text = ax.text(star_x, star_y, star['display_name'], color='white',
                                   fontsize=8, ha='center', alpha=0.8)
                texts.append(star_text)
                magnitudes.append(star['magnitude_from_observer'])

        # Add constellation region labels if close enough to Sol
        if self.label_constellations:
            const_positions = self.get_constellation_positions_from_observer()
            for const_name, const_lon, const_lat in const_positions:
                # Only show constellations in southern hemisphere
                if const_lat < 30:
                    const_lat_rad = np.radians(const_lat)
                    const_lon_rad = np.radians(const_lon)
                    const_r = 2 * np.tan(np.radians(90 + const_lat) / 2)
                    const_x = const_r * np.cos(const_lon_rad)
                    const_y = const_r * np.sin(const_lon_rad)
                    const_text = ax.text(const_x, const_y, const_name,
                           color='gray', fontsize=6, alpha=0.4,
                           ha='center', style='italic')
                    constellation_texts.append(const_text)

        # Adjust text positions to avoid overlaps
        # Brighter stars (lower magnitude) get priority for position
        adjust_label_positions(ax, texts, constellation_texts, magnitudes)

        # Add radial grid lines (latitude circles)
        if self.show_grid:
            for lat_val in range(-90, 1, self.grid_spacing):
                r_grid = 2 * np.tan(np.radians(90 + lat_val) / 2)
                circle = plt.Circle((0, 0), r_grid, fill=False, color='gray',
                                  alpha=0.3, linewidth=0.5)
                ax.add_patch(circle)

            # Add longitude lines
            for lon_val in range(0, 360, self.grid_spacing):
                lon_rad = np.radians(lon_val)
                r_max = 2 * np.tan(np.radians(90 + 30) / 2)  # Extend to +30° lat
                ax.plot([0, r_max * np.cos(lon_rad)],
                       [0, r_max * np.sin(lon_rad)],
                       color='gray', alpha=0.3, linewidth=0.5)

        # Title
        ax.set_title(f'South Polar View from {self.system_name}\nStereographic Projection',
                    color='white', fontsize=14, pad=20)

        # Equal aspect ratio
        ax.set_aspect('equal')
        max_r = 2 * np.tan(np.radians(90 + 30) / 2)
        ax.set_xlim(-max_r, max_r)
        ax.set_ylim(-max_r, max_r)

        # Hide axes
        ax.set_xticks([])
        ax.set_yticks([])
        ax.spines['top'].set_visible(False)
        ax.spines['right'].set_visible(False)
        ax.spines['bottom'].set_visible(False)
        ax.spines['left'].set_visible(False)

        plt.tight_layout()
        plt.savefig(filename, dpi=150, facecolor='black', edgecolor='none')
        plt.close()

        print(f"Saved: {filename}", file=sys.stderr)

    def generate_all_charts(self):
        """
        Generate all three star chart projections and save to files.
        """
        if self.stars_data is None:
            raise RuntimeError("Must call process() before generating charts")

        import os

        # Create output directory if it doesn't exist
        os.makedirs(self.output_dir, exist_ok=True)

        # Generate safe filename base
        safe_name = "".join(c if c.isalnum() or c in (' ', '_') else '_'
                           for c in self.system_name)
        safe_name = safe_name.replace(' ', '_').lower()

        print(f"\n=== Generating Star Charts ===\n", file=sys.stderr)

        # Generate all three projections
        allsky_path = os.path.join(self.output_dir, f"{safe_name}_allsky.png")
        north_path = os.path.join(self.output_dir, f"{safe_name}_north_polar.png")
        south_path = os.path.join(self.output_dir, f"{safe_name}_south_polar.png")

        self.plot_equirectangular(self.stars_data, allsky_path)
        self.plot_north_polar(self.stars_data, north_path)
        self.plot_south_polar(self.stars_data, south_path)

        print(f"\n=== All Charts Generated ===\n", file=sys.stderr)
        return [allsky_path, north_path, south_path]

    def generate_table(self, n: int, format: str = 'pretty', include_names: bool = True,
                       sort_by: str = 'brightness', use_parsecs: bool = False) -> str:
        """
        Generate a table of stars from observer's position.

        Parameters:
        -----------
        n : int
            Number of stars to include
        format : str
            Output format: 'pretty', 'csv', or 'html'
        include_names : bool
            Whether to look up star names
        sort_by : str
            Sort order: 'brightness' or 'distance'
        use_parsecs : bool
            If True, show distances in parsecs; otherwise light-years

        Returns:
        --------
        str : Formatted table
        """
        if self.stars_data is None:
            raise RuntimeError("Must call process() before generating table")

        # Conversion factor
        PC_TO_LY = 3.26156

        # Exclude Sol, sort appropriately, take top N
        non_sol = self.stars_data[self.stars_data.get('name', pd.Series()) != 'Sol']
        sort_col = 'magnitude_from_observer' if sort_by == 'brightness' else 'distance_from_observer_pc'
        selected = non_sol.nsmallest(n, sort_col).copy()

        # Build output data
        # Distance column setup
        dist_col = 'Distance (pc)' if use_parsecs else 'Distance (ly)'

        rows = []
        for hip_id, star in selected.iterrows():
            if include_names:
                name, sp_type, _ = get_star_info(hip_id, self.star_name_cache)
            else:
                name, sp_type = f"HIP {hip_id}", None
            
            dist_pc = star['distance_from_observer_pc']
            dist_val = dist_pc if use_parsecs else dist_pc * PC_TO_LY
            
            rows.append({
                'Name': name,
                'HIP': hip_id,
                'Spectral': sp_type or '',
                'Magnitude': round(star['magnitude_from_observer'], 2),
                dist_col: round(dist_val, 2),
                'Gal. Lon': round(star['galactic_lon'], 2),
                'Gal. Lat': round(star['galactic_lat'], 2)
            })

        if include_names:
            save_star_name_cache(self.star_name_cache)

        df = pd.DataFrame(rows)

        if format == 'csv':
            return df.to_csv(index=False)
        elif format == 'html':
            return df.to_html(index=False)
        else:  # pretty
            return df.to_string(index=False)

    def process(self):
        """
        Main processing pipeline for Phase 1: Core Calculations.

        Returns:
        --------
        pd.DataFrame : Processed star data ready for plotting
        """
        print(f"\n=== Processing Star Chart for {self.system_name} ===\n", file=sys.stderr)

        # Load catalog
        hipparcos_df = self.load_hipparcos_catalog()

        # Calculate star positions from observer
        stars_df = self.calculate_star_positions_from_observer(hipparcos_df)

        # Add Sol
        stars_df = self.add_sol(stars_df)

        # Transform to galactic coordinates
        stars_df = self.transform_to_galactic_coordinates(stars_df)

        # Store for later use
        self.stars_data = stars_df

        print(f"\n=== Processing Complete ===", file=sys.stderr)
        print(f"Total visible stars: {len(stars_df)}", file=sys.stderr)
        print(f"Magnitude range: {stars_df['magnitude_from_observer'].min():.2f} to {stars_df['magnitude_from_observer'].max():.2f}", file=sys.stderr)

        return stars_df


def generate_star_charts(ra: float, dec: float, distance: float, **kwargs):
    """
    Main entry point for generating star charts.

    Parameters:
    -----------
    ra : float
        Right Ascension in decimal degrees (0-360)
    dec : float
        Declination in decimal degrees (-90 to +90)
    distance : float
        Distance from Sol in parsecs
    **kwargs :
        Additional optional parameters (see StarChartGenerator.__init__)
    """
    generator = StarChartGenerator(ra, dec, distance, **kwargs)
    stars_data = generator.process()

    # Generate all chart visualizations
    chart_paths = generator.generate_all_charts()

    return stars_data, chart_paths


def _spherical_to_cartesian_au(ra_deg: float, dec_deg: float, distance_pc: float) -> np.ndarray:
    """Convert RA/Dec/distance to Cartesian coordinates in AU from Sol."""
    ra_rad = np.radians(ra_deg)
    dec_rad = np.radians(dec_deg)
    distance_au = distance_pc * 206265.0
    x = distance_au * np.cos(dec_rad) * np.cos(ra_rad)
    y = distance_au * np.cos(dec_rad) * np.sin(ra_rad)
    z = distance_au * np.sin(dec_rad)
    return np.array([x, y, z])


def calculate_journey_encounters(
    start_ra: float, start_dec: float, start_distance_pc: float,
    end_ra: float, end_dec: float, end_distance_pc: float,
    threshold_pc: float = 1.0
) -> Tuple[pd.DataFrame, float]:
    """
    Find stars whose closest approach to the straight-line journey path is within
    threshold_pc parsecs.

    Parameters:
    -----------
    start_ra, start_dec, start_distance_pc : float
        Starting point in equatorial coordinates (distance 0 = Sol)
    end_ra, end_dec, end_distance_pc : float
        Ending point in equatorial coordinates
    threshold_pc : float
        Maximum closest-approach distance to include (default: 1.0 pc)

    Returns:
    --------
    tuple : (encounters DataFrame, journey_length_pc)
        DataFrame indexed by HIP ID (0 for Sol) with columns:
        approach_pc, journey_fraction, dist_along_journey_pc, magnitude_at_approach
    """
    PC_TO_AU = 206265.0
    PC_TO_LY = 3.26156
    EXCLUSION_PC = 0.5 / PC_TO_LY  # 0.5 ly exclusion zone around each endpoint

    start_pos = _spherical_to_cartesian_au(start_ra, start_dec, start_distance_pc)
    end_pos = _spherical_to_cartesian_au(end_ra, end_dec, end_distance_pc)

    seg = end_pos - start_pos
    seg_len_sq = float(np.dot(seg, seg))
    journey_pc = float(np.sqrt(max(seg_len_sq, 0.0))) / PC_TO_AU

    threshold_au = threshold_pc * PC_TO_AU
    exclusion_au = EXCLUSION_PC * PC_TO_AU

    print("Loading Hipparcos catalog...", file=sys.stderr)
    with load.open(hipparcos.URL) as f:
        hip_df = hipparcos.load_dataframe(f)
    print(f"Loaded {len(hip_df)} stars", file=sys.stderr)

    valid = hip_df[hip_df['parallax_mas'] > 0].copy()
    valid['distance_from_earth_pc'] = 1000.0 / valid['parallax_mas']

    ra_rad = np.radians(valid['ra_degrees'].values)
    dec_rad = np.radians(valid['dec_degrees'].values)
    dist_au_arr = valid['distance_from_earth_pc'].values * PC_TO_AU
    x_arr = dist_au_arr * np.cos(dec_rad) * np.cos(ra_rad)
    y_arr = dist_au_arr * np.cos(dec_rad) * np.sin(ra_rad)
    z_arr = dist_au_arr * np.sin(dec_rad)

    # Vectorized closest-point-on-segment calculation
    dx = x_arr - start_pos[0]
    dy = y_arr - start_pos[1]
    dz = z_arr - start_pos[2]

    if seg_len_sq > 0:
        t_raw = (dx * seg[0] + dy * seg[1] + dz * seg[2]) / seg_len_sq
        t = np.clip(t_raw, 0.0, 1.0)
    else:
        t_raw = np.zeros(len(valid))
        t = t_raw

    closest_x = start_pos[0] + t * seg[0]
    closest_y = start_pos[1] + t * seg[1]
    closest_z = start_pos[2] + t * seg[2]

    approach_au = np.sqrt(
        (x_arr - closest_x)**2 + (y_arr - closest_y)**2 + (z_arr - closest_z)**2
    )
    dist_from_start_sq = dx**2 + dy**2 + dz**2
    ex = x_arr - end_pos[0]
    ey = y_arr - end_pos[1]
    ez = z_arr - end_pos[2]
    dist_from_end_sq = ex**2 + ey**2 + ez**2

    # t_raw outside [0, 1] means the perpendicular foot falls behind the start or
    # past the destination — the ship would have to reverse or overshoot to divert,
    # so these are excluded regardless of how close they are to an endpoint.
    mask = (
        (approach_au <= threshold_au) &
        (t_raw >= 0.0) &
        (t_raw <= 1.0) &
        (dist_from_start_sq >= exclusion_au**2) &
        (dist_from_end_sq >= exclusion_au**2)
    )

    valid['_approach_au'] = approach_au
    valid['_t'] = t
    encounters = valid[mask].copy()
    encounters['approach_pc'] = encounters['_approach_au'] / PC_TO_AU
    encounters['journey_fraction'] = encounters['_t']
    encounters['dist_along_journey_pc'] = encounters['_t'] * journey_pc
    encounters = encounters[encounters['approach_pc'] > 1e-6].copy()
    encounters['magnitude_at_approach'] = (
        encounters['magnitude'] +
        5 * np.log10(encounters['approach_pc'] / encounters['distance_from_earth_pc'])
    )
    encounters = encounters.drop(columns=['_approach_au', '_t'])

    # Check whether Sol (at the origin) falls within the threshold
    sol_pos = np.array([0.0, 0.0, 0.0])
    if seg_len_sq > 0:
        t_sol_raw = float(np.dot(sol_pos - start_pos, seg) / seg_len_sq)
    else:
        t_sol_raw = 0.0
    t_sol = float(np.clip(t_sol_raw, 0.0, 1.0))
    sol_approach_au = float(np.linalg.norm(start_pos + t_sol * seg))
    dist_sol_start = float(np.linalg.norm(start_pos))
    dist_sol_end = float(np.linalg.norm(end_pos))
    if (0.0 <= t_sol_raw <= 1.0 and
            sol_approach_au <= threshold_au and
            dist_sol_start >= exclusion_au and dist_sol_end >= exclusion_au):
        sol_approach_pc = sol_approach_au / PC_TO_AU
        if sol_approach_pc > 1e-6:
            sol_mag = 4.83 + 5 * np.log10(sol_approach_pc) - 5
            sol_df = pd.DataFrame([{
                'approach_pc': sol_approach_pc,
                'journey_fraction': t_sol,
                'dist_along_journey_pc': t_sol * journey_pc,
                'magnitude_at_approach': sol_mag,
            }], index=pd.Index([0]))
            encounters = pd.concat([encounters, sol_df])

    print(f"Found {len(encounters)} stars within {threshold_pc:.2f} pc of the journey path", file=sys.stderr)
    return encounters, journey_pc


def generate_journey_table(
    encounters: pd.DataFrame,
    journey_pc: float,
    start_name: str,
    end_name: str,
    format: str = 'pretty',
    use_parsecs: bool = False,
    use_simbad: bool = True
) -> str:
    """
    Format journey encounter stars as a table.

    Parameters:
    -----------
    encounters : pd.DataFrame
        DataFrame from calculate_journey_encounters (indexed by HIP ID, 0 = Sol)
    journey_pc : float
        Total journey distance in parsecs
    start_name, end_name : str
        Display names for the departure and destination points
    format : str
        Output format: 'pretty', 'csv', or 'html'
    use_parsecs : bool
        If True, show distances in parsecs; otherwise light-years
    use_simbad : bool
        Whether to query SIMBAD for star names and spectral types

    Returns:
    --------
    str : Formatted table string
    """
    PC_TO_LY = 3.26156
    dist_unit = 'pc' if use_parsecs else 'ly'
    journey_display = journey_pc if use_parsecs else journey_pc * PC_TO_LY

    header = f"\nJourney: {start_name} → {end_name}  ({journey_display:.2f} {dist_unit})\n"

    if encounters.empty:
        return header + "No stars found within the specified approach distance.\n"

    star_name_cache = load_star_name_cache() if use_simbad else {}
    rows = []
    for hip_id, enc in encounters.iterrows():
        hip_id = int(hip_id)
        if hip_id == 0:
            name, sp_type = 'Sol', 'G2V'
        elif use_simbad:
            name, sp_type, _ = get_star_info(hip_id, star_name_cache)
        else:
            name, sp_type = f"HIP {hip_id}", None

        approach_pc = enc['approach_pc']
        approach_val = approach_pc if use_parsecs else approach_pc * PC_TO_LY
        dist_along_pc = enc['dist_along_journey_pc']
        dist_along_val = dist_along_pc if use_parsecs else dist_along_pc * PC_TO_LY

        rows.append({
            'Name': name,
            'HIP': '' if hip_id == 0 else hip_id,
            'Spectral': sp_type or '',
            f'Approach ({dist_unit})': round(approach_val, 3),
            f'From Start ({dist_unit})': round(dist_along_val, 2),
            'Journey %': round(enc['journey_fraction'] * 100, 1),
            'Mag at Approach': round(enc['magnitude_at_approach'], 2),
        })

    if use_simbad:
        save_star_name_cache(star_name_cache)

    df = pd.DataFrame(rows)
    df = df.sort_values(f'From Start ({dist_unit})').reset_index(drop=True)

    if format == 'csv':
        return df.to_csv(index=False)
    elif format == 'html':
        return df.to_html(index=False)
    else:
        return header + df.to_string(index=False)


def main():
    """Command-line interface."""
    parser = argparse.ArgumentParser(
        description='Generate star charts from arbitrary positions in 3D space'
    )
    parser.add_argument('--ra', type=float,
                       help='Right Ascension in decimal degrees (0-360)')
    parser.add_argument('--dec', type=float,
                       help='Declination in decimal degrees (-90 to +90)')
    parser.add_argument('--distance', type=float,
                       help='Distance from Sol in parsecs')
    parser.add_argument('--name', type=str, default=None,
                       help='System name for labeling')
    parser.add_argument('--magnitude-cutoff', type=float, default=6.0,
                       help='Faintest magnitude to display (default: 6.0)')
    parser.add_argument('--label-threshold', type=float, default=1.0,
                       help='Magnitude threshold for labeling stars (default: 1.0)')
    parser.add_argument('--constellation-lines', action='store_true',
                       help='Show constellation lines')
    parser.add_argument('--grid-spacing', type=int, default=15,
                       help='Grid spacing in degrees (default: 15)')
    parser.add_argument('--output-dir', type=str, default='.',
                       help='Output directory (default: current directory)')
    parser.add_argument('--clear-cache', action='store_true',
                       help='Clear the star name cache and exit')
    parser.add_argument('--no-grid', action='store_true',
                       help='Hide grid lines and axes')
    parser.add_argument('--no-star-labels', action='store_true',
                       help='Hide star name labels')
    parser.add_argument('--no-constellation-labels', action='store_true',
                       help='Hide constellation region labels')
    parser.add_argument('--no-labels', action='store_true',
                       help='Hide all labels (stars, constellations, Sol)')
    parser.add_argument('--no-sol-highlight', action='store_true',
                       help='Show Sol as a regular star instead of highlighted yellow')
    parser.add_argument('--table', type=int, metavar='N',
                       help='Generate table of N brightest stars (skips chart generation)')
    parser.add_argument('--format', choices=['pretty', 'csv', 'html'], default='pretty',
                       help='Table output format (default: pretty)')
    parser.add_argument('--sort-by', choices=['brightness', 'distance'], default='brightness',
                       help='Table sort order (default: brightness)')
    parser.add_argument('--no-simbad', action='store_true',
                       help='Skip SIMBAD lookups (names, spectral types) in table output')
    parser.add_argument('--parsecs', action='store_true',
                       help='Show distances in parsecs instead of light-years')

    # Journey mode arguments
    parser.add_argument('--journey', action='store_true',
                       help='Find stars near the straight-line path between two points in space')
    parser.add_argument('--from-ra', type=float,
                       help='Journey start: Right Ascension in degrees (default: Sol, 0)')
    parser.add_argument('--from-dec', type=float,
                       help='Journey start: Declination in degrees (default: Sol, 0)')
    parser.add_argument('--from-distance', type=float,
                       help='Journey start: distance from Sol in parsecs (default: 0 = Sol)')
    parser.add_argument('--from-name', type=str,
                       help='Display name for the departure point')
    parser.add_argument('--to-ra', type=float,
                       help='Journey end: Right Ascension in degrees')
    parser.add_argument('--to-dec', type=float,
                       help='Journey end: Declination in degrees')
    parser.add_argument('--to-distance', type=float,
                       help='Journey end: distance from Sol in parsecs')
    parser.add_argument('--to-name', type=str,
                       help='Display name for the destination point')
    parser.add_argument('--approach-distance', type=float, default=1.0,
                       help='Maximum closest-approach distance in parsecs (default: 1.0)')

    args = parser.parse_args()

    if args.clear_cache:
        if os.path.exists(STAR_NAME_CACHE_FILE):
            os.remove(STAR_NAME_CACHE_FILE)
            print(f"Cleared cache: {STAR_NAME_CACHE_FILE}")
        else:
            print("No cache file to clear.")
        return

    if args.journey:
        if args.to_ra is None or args.to_dec is None or args.to_distance is None:
            parser.error("--journey requires --to-ra, --to-dec, and --to-distance")

        from_ra = args.from_ra if args.from_ra is not None else 0.0
        from_dec = args.from_dec if args.from_dec is not None else 0.0
        from_distance = args.from_distance if args.from_distance is not None else 0.0

        if args.from_name:
            from_name = args.from_name
        elif args.from_distance is None:
            from_name = 'Sol'
        else:
            from_name = 'Start'
        to_name = args.to_name or 'Destination'

        encounters, journey_pc = calculate_journey_encounters(
            from_ra, from_dec, from_distance,
            args.to_ra, args.to_dec, args.to_distance,
            threshold_pc=args.approach_distance
        )
        print(generate_journey_table(
            encounters, journey_pc,
            from_name, to_name,
            format=args.format,
            use_parsecs=args.parsecs,
            use_simbad=not args.no_simbad
        ))
        return

    if args.ra is None or args.dec is None or args.distance is None:
        parser.error("--ra, --dec, and --distance are required unless using --clear-cache or --journey")

    if args.table:
        # Table-only mode
        # For distance sorting, use high cutoff to capture dim nearby stars (unless user specified)
        mag_cutoff = args.magnitude_cutoff
        if args.sort_by == 'distance' and args.magnitude_cutoff == 6.0:
            mag_cutoff = 15.0  # Capture dim nearby stars
        
        generator = StarChartGenerator(
            ra=args.ra,
            dec=args.dec,
            distance=args.distance,
            system_name=args.name,
            magnitude_cutoff=mag_cutoff
        )
        generator.process()
        print(generator.generate_table(args.table, args.format, not args.no_simbad, args.sort_by, args.parsecs))
    else:
        # Chart generation mode
        generate_star_charts(
            ra=args.ra,
            dec=args.dec,
            distance=args.distance,
            system_name=args.name,
            magnitude_cutoff=args.magnitude_cutoff,
            label_magnitude_threshold=args.label_threshold,
            show_constellations=args.constellation_lines,
            grid_spacing=args.grid_spacing,
            output_dir=args.output_dir,
            show_grid=not args.no_grid,
            show_star_labels=not (args.no_labels or args.no_star_labels),
            show_constellation_labels=not (args.no_labels or args.no_constellation_labels),
            show_sol_label=not args.no_labels,
            highlight_sol=not args.no_sol_highlight
        )


if __name__ == '__main__':
    main()