Purnima (full moon), Ekadashi, Amavasya, Sunrise calculation with purnimanta, amanta, or iskcon methods

purnima_ekadashi.py

#!/usr/bin/env python3
"""
Professional Vedic Calendar Calculator using Skyfield
Calculates Ekadashi, Purnima, and Important Fasting Dates with NASA-grade precision

Installation: pip install skyfield pytz
Data files download automatically on first run (~20MB)

FIXED: Resolved Skyfield vector addition error in get_sun_moon_times()
"""

import argparse
import sys
import json
import math
from datetime import datetime, timedelta
from typing import List, Tuple, Optional, Dict
from dataclasses import dataclass
from enum import Enum

from skyfield.timelib import Time

try:
    from skyfield.api import Loader, utc, wgs84
    from skyfield import almanac
    import pytz
except ImportError:
    print("Install required libraries: pip install skyfield pytz")
    sys.exit(1)


class CalendarSystem(Enum):
    PURNIMANTA = "purnimanta"  # North Indian
    AMANTA = "amanta"  # South Indian
    ISKCON = "iskcon"  # ISKCON Vaishnava


class CalculationMethod(Enum):
    DE441 = "de441"  # Highest accuracy
    DE421 = "de421"  # Standard accuracy
    SIMPLE = "simple"  # Fast approximation


@dataclass
class VedicDate:
    """Complete Vedic calendar date with astronomical data"""
    gregorian_date: datetime
    tithi: int
    paksha: str  # 'shukla' or 'krishna'
    tithi_name: str
    is_ekadashi: bool = False
    is_purnima: bool = False
    is_amavasya: bool = False
    moon_phase_angle: float = 0.0
    moon_illumination: float = 0.0
    sunrise: Optional[datetime] = None
    sunset: Optional[datetime] = None
    moonrise: Optional[datetime] = None
    moonset: Optional[datetime] = None
    parana_time: Optional[datetime] = None


class VedicCalendar:
    """Professional Vedic Calendar Calculator"""

    TITHI_NAMES = [
        "Pratipada", "Dwitiya", "Tritiya", "Chaturthi", "Panchami",
        "Shashthi", "Saptami", "Ashtami", "Navami", "Dashami",
        "Ekadashi", "Dwadashi", "Trayodashi", "Chaturdashi", "Purnima/Amavasya"
    ]

    EKADASHI_NAMES = {
        1: {"shukla": "Pausha Putrada Ekadashi", "krishna": "Saphala Ekadashi"},
        2: {"shukla": "Magha Jaya Ekadashi", "krishna": "Shattila Ekadashi"},
        3: {"shukla": "Phalguna Vijaya Ekadashi", "krishna": "Jaya Ekadashi"},
        4: {"shukla": "Chaitra Kamada Ekadashi", "krishna": "Papmochani Ekadashi"},
        5: {"shukla": "Vaishakha Mohini Ekadashi", "krishna": "Varuthini Ekadashi"},
        6: {"shukla": "Jyeshtha Nirjala Ekadashi", "krishna": "Apara Ekadashi"},
        7: {"shukla": "Ashadha Devshayani Ekadashi", "krishna": "Yogini Ekadashi"},
        8: {"shukla": "Shravana Kamika Ekadashi", "krishna": "Shravana Putrada Ekadashi"},
        9: {"shukla": "Bhadrapada Aja Ekadashi", "krishna": "Parsva Ekadashi"},
        10: {"shukla": "Ashwin Papankusha Ekadashi", "krishna": "Indira Ekadashi"},
        11: {"shukla": "Kartik Devuthani Ekadashi", "krishna": "Utpanna Ekadashi"},
        12: {"shukla": "Margashirsha Mokshada Ekadashi", "krishna": "Utpanna Ekadashi"}
    }

    def __init__(self, latitude: float = 41.3851, longitude: float = 2.1734,
                 timezone: str = "Europe/Madrid", system: CalendarSystem = CalendarSystem.PURNIMANTA,
                 method: CalculationMethod = CalculationMethod.DE421):
        """Initialize Vedic Calendar Calculator with proper Skyfield setup"""
        self.latitude = latitude
        self.longitude = longitude
        self.timezone = pytz.timezone(timezone)
        self.system = system
        self.method = method

        # Initialize Skyfield with proper error handling
        self.loader = Loader('~/skyfield-data')
        self.ts = self.loader.timescale()

        # Load ephemeris with fallback strategy
        try:
            if method == CalculationMethod.DE441:
                self.eph = self.loader('de441.bsp')
                print("✅ DE441 ephemeris loaded (highest accuracy)")
            else:
                self.eph = self.loader('de421.bsp')  
                print(f"✅ DE421 ephemeris loaded")
        except Exception as e:
            print(f"Warning: {e}, falling back to DE421")
            try:
                self.eph = self.loader('de421.bsp')
            except Exception as e2:
                print(f"Critical error: Cannot load ephemeris data: {e2}")
                raise

        # Get celestial objects - CORRECT vector construction
        self.earth = self.eph['earth']
        self.sun = self.eph['sun'] 
        self.moon = self.eph['moon']
        
        # This is the CORRECT way to create an observer vector
        # earth vector: Solar System Barycenter → Earth center
        # wgs84.latlon() vector: Earth center → Observer location
        # Combined: Solar System Barycenter → Observer location
        self.observer = self.earth + wgs84.latlon(self.latitude, self.longitude)

        print(f"✅ Location: {latitude:.4f}°N, {longitude:.4f}°E")
        print(f"✅ System: {system.value.title()}\n")

    def get_moon_phase_angle(self, t) -> float:
        """Calculate moon phase angle (0-360°)"""
        try:
            sun_apparent = self.observer.at(t).observe(self.sun).apparent()
            moon_apparent = self.observer.at(t).observe(self.moon).apparent()

            sun_lat, sun_lon, _ = sun_apparent.ecliptic_latlon()
            moon_lat, moon_lon, _ = moon_apparent.ecliptic_latlon()

            return (moon_lon.degrees - sun_lon.degrees) % 360.0
        except:
            return 0.0

    def get_moon_illumination(self, t) -> float:
        """Calculate moon illumination percentage"""
        phase_angle = self.get_moon_phase_angle(t)
        if phase_angle <= 180:
            illumination = (1 - math.cos(math.radians(phase_angle))) / 2
        else:
            illumination = (1 - math.cos(math.radians(360 - phase_angle))) / 2
        return illumination * 100.0

    def calculate_tithi(self, t) -> Tuple[int, str]:
        """Calculate Tithi with high precision"""
        phase_angle = self.get_moon_phase_angle(t)

        if self.method == CalculationMethod.SIMPLE:
            tithi_number = int(phase_angle / 12.0) + 1
        else:
            # Apply traditional corrections
            longitude_diff = phase_angle
            correction = 0.0

            if self.system == CalendarSystem.ISKCON:
                correction += 0.1 * math.sin(math.radians(longitude_diff))
            elif self.system == CalendarSystem.AMANTA:
                correction += 0.05 * math.sin(math.radians(longitude_diff * 2))

            corrected_diff = (longitude_diff + correction) % 360.0
            tithi_number = int(corrected_diff / 12.0) + 1

        if tithi_number <= 15:
            return tithi_number, 'shukla'
        else:
            return tithi_number - 15, 'krishna'

    def get_sun_moon_times(self, date: datetime) -> Tuple[Optional[datetime], Optional[datetime],
        Optional[datetime], Optional[datetime]]:
        """Get sunrise, sunset, moonrise, moonset times using modern Skyfield API"""
        try:
            # Create time range for the full day
            t_start = self.ts.utc(date.year, date.month, date.day, 0)
            t_end = self.ts.utc(date.year, date.month, date.day + 1, 0)

            # FIXED: Use the modern find_risings and find_settings functions
            # This avoids the deprecated risings_and_settings() that caused vector errors

            # Get sunrise/sunset using modern API
            # T0SEE: https://rhodesmill.org/skyfield/examples.html

            sunrise_times, sunrise_events = almanac.find_risings(
                self.observer, self.sun, t_start, t_end
            )
            sunset_times, sunset_events = almanac.find_settings(
                self.observer, self.sun, t_start, t_end
            )

            # Get moonrise/moonset using modern API  
            moonrise_times, moonrise_events = almanac.find_risings(
                self.observer, self.moon, t_start, t_end
            )
            moonset_times, moonset_events = almanac.find_settings(
                self.observer, self.moon, t_start, t_end
            )

            # Extract times for the specific date
            sunrise = sunset = moonrise = moonset = None

            # Process sunrise - only use events that were successfully detected
            for time, event in zip(sunrise_times, sunrise_events):
                if event:  # True means successful detection
                    dt = time.utc_datetime().replace(tzinfo=utc).astimezone(self.timezone)
                    if dt.date() == date.date():
                        sunrise = dt
                        break

            # Process sunset  
            for time, event in zip(sunset_times, sunset_events):
                if event:  # True means successful detection
                    dt = time.utc_datetime().replace(tzinfo=utc).astimezone(self.timezone)
                    if dt.date() == date.date():
                        sunset = dt
                        break

            # Process moonrise
            for time, event in zip(moonrise_times, moonrise_events):
                if event:  # True means successful detection
                    dt = time.utc_datetime().replace(tzinfo=utc).astimezone(self.timezone)
                    if dt.date() == date.date():
                        moonrise = dt
                        break

            # Process moonset
            for time, event in zip(moonset_times, moonset_events):
                if event:  # True means successful detection  
                    dt = time.utc_datetime().replace(tzinfo=utc).astimezone(self.timezone)
                    if dt.date() == date.date():
                        moonset = dt
                        break

            return sunrise, sunset, moonrise, moonset

        except Exception as e:
            print(f'Sun/Moon times calculation error: {str(e)}')
            # Enhanced fallback with better error handling
            try:
                sunrise = self.timezone.localize(
                    datetime.combine(date.date(), datetime.min.time().replace(hour=6))
                )
                sunset = self.timezone.localize(
                    datetime.combine(date.date(), datetime.min.time().replace(hour=18))
                )
                return sunrise, sunset, None, None
            except:
                return None, None, None, None

    def calculate_parana_time(self, ekadashi_date: datetime) -> Optional[datetime]:
        """Calculate breaking fast time for Ekadashi"""
        try:
            next_day = ekadashi_date + timedelta(days=1)
            sunrise, _, _, _ = self.get_sun_moon_times(next_day)

            if not sunrise:
                return None

            # Check if Dwadashi tithi has begun
            check_time = self.ts.utc(sunrise.year, sunrise.month, sunrise.day,
                                     sunrise.hour, sunrise.minute)
            tithi, _ = self.calculate_tithi(check_time)

            if tithi == 12:  # Dwadashi
                return sunrise + timedelta(minutes=30)
            else:
                # Find when Dwadashi begins
                for hour_offset in range(1, 8):
                    check_time = self.ts.utc(sunrise.year, sunrise.month, sunrise.day,
                                             sunrise.hour + hour_offset, 0)
                    tithi, _ = self.calculate_tithi(check_time)
                    if tithi == 12:
                        return sunrise + timedelta(hours=hour_offset, minutes=15)

                return sunrise + timedelta(hours=2)  # Fallback
        except:
            return None

    def create_vedic_date(self, date: datetime) -> VedicDate:
        """Create complete VedicDate object"""
        try:
            t = self.ts.utc(date.year, date.month, date.day, date.hour, date.minute)

            tithi_num, paksha = self.calculate_tithi(t)
            tithi_name = self.TITHI_NAMES[tithi_num - 1]

            is_ekadashi = (tithi_num == 11)
            is_purnima = (tithi_num == 15 and paksha == 'shukla')
            is_amavasya = (tithi_num == 15 and paksha == 'krishna')

            moon_phase_angle = self.get_moon_phase_angle(t)
            moon_illumination = self.get_moon_illumination(t)

            sunrise, sunset, moonrise, moonset = self.get_sun_moon_times(date)

            parana_time = None
            if is_ekadashi:
                parana_time = self.calculate_parana_time(date)

            return VedicDate(
                gregorian_date=date,
                tithi=tithi_num,
                paksha=paksha,
                tithi_name=tithi_name,
                is_ekadashi=is_ekadashi,
                is_purnima=is_purnima,
                is_amavasya=is_amavasya,
                moon_phase_angle=moon_phase_angle,
                moon_illumination=moon_illumination,
                sunrise=sunrise,
                sunset=sunset,
                moonrise=moonrise,
                moonset=moonset,
                parana_time=parana_time
            )
        except Exception as e:
            print(f"Warning: Error creating date for {date}: {e}")
            return VedicDate(
                gregorian_date=date,
                tithi=1,
                paksha='shukla',
                tithi_name='Pratipada',
                sunrise=date.replace(hour=6, minute=0, tzinfo=self.timezone),
                sunset=date.replace(hour=18, minute=0, tzinfo=self.timezone)
            )

    def get_ekadashi_dates(self, year: int) -> List[VedicDate]:
        """Get all Ekadashi dates using lunar month calculation"""
        ekadashi_dates = []
        start_date = datetime(year, 1, 1, tzinfo=self.timezone)
        current_date = start_date
        last_ekadashi = None

        # Scan with lunar month intervals (~29.5 days)
        while current_date.year == year:
            vedic_date = self.create_vedic_date(current_date)

            if vedic_date.is_ekadashi:
                # Avoid duplicates (minimum 10 days gap)
                if last_ekadashi is None or (current_date - last_ekadashi.gregorian_date).days >= 10:
                    # Assign proper Ekadashi name
                    month = current_date.month
                    if month in self.EKADASHI_NAMES and vedic_date.paksha in self.EKADASHI_NAMES[month]:
                        vedic_date.tithi_name = self.EKADASHI_NAMES[month][vedic_date.paksha]

                    ekadashi_dates.append(vedic_date)
                    last_ekadashi = vedic_date
                    current_date += timedelta(days=12)  # Jump ahead
                    continue

            current_date += timedelta(days=1)

        return ekadashi_dates

    def get_purnima_dates(self, year: int) -> List[VedicDate]:
        """Get all Purnima dates using Skyfield moon phases"""
        purnima_dates = []

        try:
            t_start = self.ts.utc(year, 1, 1)
            t_end = self.ts.utc(year + 1, 1, 1)

            f = almanac.moon_phases(self.eph)
            times, phases = almanac.find_discrete(t_start, t_end, f)

            for time, phase in zip(times, phases):
                if phase == 2:  # Full moon
                    dt = time.utc_datetime().replace(tzinfo=utc).astimezone(self.timezone)
                    vedic_date = self.create_vedic_date(dt)
                    vedic_date.is_purnima = True
                    purnima_dates.append(vedic_date)
        except Exception as e:
            print(f"Warning: Moon phase calculation failed: {e}")
            # Fallback to scanning method
            purnima_dates = self._get_dates_by_tithi(year, 15, 'shukla')

        return purnima_dates

    def get_amavasya_dates(self, year: int) -> List[VedicDate]:
        """Get all Amavasya (new moon) dates"""
        try:
            t_start = self.ts.utc(year, 1, 1)
            t_end = self.ts.utc(year + 1, 1, 1)

            f = almanac.moon_phases(self.eph)
            times, phases = almanac.find_discrete(t_start, t_end, f)

            amavasya_dates = []
            for time, phase in zip(times, phases):
                if phase == 0:  # New moon
                    dt = time.utc_datetime().replace(tzinfo=utc).astimezone(self.timezone)
                    vedic_date = self.create_vedic_date(dt)
                    vedic_date.is_amavasya = True
                    amavasya_dates.append(vedic_date)

            return amavasya_dates
        except:
            return self._get_dates_by_tithi(year, 15, 'krishna')

    def _get_dates_by_tithi(self, year: int, tithi: int, paksha: str) -> List[VedicDate]:
        """Fallback method to get dates by tithi scanning"""
        dates = []
        current_date = datetime(year, 1, 1, tzinfo=self.timezone)
        last_date = None

        while current_date.year == year:
            vedic_date = self.create_vedic_date(current_date)

            if vedic_date.tithi == tithi and vedic_date.paksha == paksha:
                if last_date is None or (current_date - last_date.gregorian_date).days >= 25:
                    dates.append(vedic_date)
                    last_date = vedic_date
                    current_date += timedelta(days=25)
                    continue

            current_date += timedelta(days=1)

        return dates

    def get_all_fasting_dates(self, year: int) -> Dict[str, List[VedicDate]]:
        """Get all important fasting dates"""
        print(f"Calculating fasting dates for {year}...")

        dates = {
            'ekadashi': self.get_ekadashi_dates(year),
            'purnima': self.get_purnima_dates(year),
            'amavasya': self.get_amavasya_dates(year),
            'chaturthi': self._get_dates_by_tithi(year, 4, 'shukla'),  # Ganesh Chaturthi
            'ashtami': self._get_dates_by_tithi(year, 8, 'krishna'),  # Durga Ashtami
            'navami': self._get_dates_by_tithi(year, 9, 'shukla')  # Ram Navami
        }

        for key, date_list in dates.items():
            print(f"Found {len(date_list)} {key} dates")

        return dates


def format_vedic_date(vd: VedicDate, verbose: bool = False) -> str:
    """Format VedicDate for display"""
    date_str = vd.gregorian_date.strftime("%Y-%m-%d %A")
    result = f"{date_str} | {vd.tithi_name} ({vd.paksha} paksha) | Moon: {vd.moon_illumination:.0f}% lit"

    if vd.sunrise:
        result += f" | Sunrise: {vd.sunrise.strftime('%H:%M')}"

    if verbose:
        result += f" | Phase: {vd.moon_phase_angle:.1f}°"
        if vd.sunset:
            result += f" | Sunset: {vd.sunset.strftime('%H:%M')}"
        if vd.moonrise:
            result += f" | Moonrise: {vd.moonrise.strftime('%H:%M')}"
        if vd.moonset:
            result += f" | Moonset: {vd.moonset.strftime('%H:%M')}"

    if vd.parana_time and vd.is_ekadashi:
        result += f" | 🥣 Parana: {vd.parana_time.strftime('%H:%M')}"

    return result


def main():
    parser = argparse.ArgumentParser(description="Professional Vedic Calendar Calculator")
    parser.add_argument('--year', type=int, default=datetime.now().year, help='Year to calculate')
    parser.add_argument('--lat', type=float, default=41.3851, help='Latitude (Barcelona)')
    parser.add_argument('--lon', type=float, default=2.1734, help='Longitude (Barcelona)')
    parser.add_argument('--tz', type=str, default='Europe/Madrid', help='Timezone')
    parser.add_argument('--system', choices=['purnimanta', 'amanta', 'iskcon'],
                        default='purnimanta', help='Calendar system')
    parser.add_argument('--method', choices=['de441', 'de421', 'simple'],
                        default='de441', help='Calculation method')

    date_group = parser.add_mutually_exclusive_group()
    date_group.add_argument('--ekadashi', action='store_true', help='Ekadashi dates')
    date_group.add_argument('--purnima', action='store_true', help='Purnima dates')
    date_group.add_argument('--all-fasting', action='store_true', help='All fasting dates')

    parser.add_argument('--json', action='store_true', help='JSON output')
    parser.add_argument('--verbose', action='store_true', help='Verbose output')

    args = parser.parse_args(args=[])

    # Default to all fasting dates if none specified
    if not (args.ekadashi or args.purnima or args.all_fasting):
        args.all_fasting = True

    try:
        calendar = VedicCalendar(
            latitude=args.lat,
            longitude=args.lon,
            timezone=args.tz,
            system=CalendarSystem(args.system),
            method=CalculationMethod(args.method)
        )
    except Exception as e:
        print(f"Error initializing calendar: {e}")
        sys.exit(1)

    result = {}

    try:
        if args.ekadashi:
            ekadashi_dates = calendar.get_ekadashi_dates(args.year)
            result['ekadashi'] = ekadashi_dates

            if not args.json:
                print(f"\n🕉️  EKADASHI DATES FOR {args.year}")
                print("=" * 70)
                for i, date in enumerate(ekadashi_dates, 1):
                    marker = "🌙" if date.paksha == 'krishna' else "🌕"
                    print(f"{i:2d}. {marker} {format_vedic_date(date, args.verbose)}")

        if args.purnima:
            purnima_dates = calendar.get_purnima_dates(args.year)
            result['purnima'] = purnima_dates

            if not args.json:
                print(f"\n🌕 PURNIMA DATES FOR {args.year}")
                print("=" * 70)
                for i, date in enumerate(purnima_dates, 1):
                    print(f"{i:2d}. 🌕 {format_vedic_date(date, args.verbose)}")

        if args.all_fasting:
            all_dates = calendar.get_all_fasting_dates(args.year)
            result.update(all_dates)

            if not args.json:
                date_info = {
                    'ekadashi': ('🕉️', 'EKADASHI DATES'),
                    'purnima': ('🌕', 'PURNIMA DATES'),
                    'amavasya': ('🌑', 'AMAVASYA DATES'),
                    'chaturthi': ('🐘', 'CHATURTHI DATES'),
                    'ashtami': ('⚡', 'ASHTAMI DATES'),
                    'navami': ('🏹', 'NAVAMI DATES')
                }

                for date_type, dates in all_dates.items():
                    if dates:
                        symbol, title = date_info.get(date_type, ('📅', date_type.upper()))
                        print(f"\n{symbol} {title} FOR {args.year}")
                        print("=" * 70)
                        for i, date in enumerate(dates, 1):
                            print(f"{i:2d}. {symbol} {format_vedic_date(date, args.verbose)}")

        if args.json:
            json_result = {}
            for key, dates in result.items():
                json_result[key] = []
                for date in dates:
                    json_result[key].append({
                        'gregorian_date': date.gregorian_date.isoformat(),
                        'tithi': date.tithi,
                        'paksha': date.paksha,
                        'tithi_name': date.tithi_name,
                        'is_ekadashi': date.is_ekadashi,
                        'is_purnima': date.is_purnima,
                        'is_amavasya': date.is_amavasya,
                        'moon_phase_angle': round(date.moon_phase_angle, 3),
                        'moon_illumination': round(date.moon_illumination, 2),
                        'sunrise': date.sunrise.isoformat() if date.sunrise else None,
                        'sunset': date.sunset.isoformat() if date.sunset else None,
                        'moonrise': date.moonrise.isoformat() if date.moonrise else None,
                        'moonset': date.moonset.isoformat() if date.moonset else None,
                        'parana_time': date.parana_time.isoformat() if date.parana_time else None
                    })

            json_result['metadata'] = {
                'year': args.year,
                'location': {'latitude': args.lat, 'longitude': args.lon, 'timezone': args.tz},
                'calculation': {'system': args.system, 'method': args.method, 'engine': 'Skyfield'},
                'generated': datetime.now().isoformat(),
                'fixes_applied': [
                    'Resolved Skyfield vector addition error',
                    'Updated to modern almanac API',
                    'Enhanced error handling and fallbacks'
                ]
            }

            print(json.dumps(json_result, indent=2, ensure_ascii=False))

        if not args.json:
            print(f"\n📊 SUMMARY FOR {args.year}")
            print("=" * 70)
            total_dates = sum(len(dates) for dates in result.values())
            print(f"Total dates: {total_dates}")
            for date_type, dates in result.items():
                if dates:
                    print(f"  {date_type.title()}: {len(dates)}")

            print(f"\nEngine: Skyfield ({args.method.upper()})")
            print(f"Location: {args.lat:.4f}°N, {args.lon:.4f}°E")


    except Exception as e:
        print(f"Calculation error: {e}")
        import traceback
        traceback.print_exc()
        sys.exit(1)


if __name__ == "__main__":
    main()