import numpy as np
import math

from scipy.signal import savgol_filter, find_peaks
from .constants import WAVE_MAP

def strip_ka2(x_data, y_data, ratio=0.5):
    """Rachinger法によるCuKα2成分の除去"""
    lambda1 = WAVE_MAP["CuKα1 (1.5406 Å)"]
    lambda2 = WAVE_MAP["CuKα2"] 
    
    y_stripped = np.copy(y_data)
    
    for i, two_theta in enumerate(x_data):
        sin_theta2 = math.sin(math.radians(two_theta / 2.0))
        sin_theta1 = (lambda1 / lambda2) * sin_theta2
        
        if sin_theta1 >= 1.0 or sin_theta1 <= -1.0:
            continue
            
        two_theta1 = math.degrees(math.asin(sin_theta1)) * 2.0
        
        if two_theta1 < x_data[0]:
            continue
            
        idx = np.searchsorted(x_data, two_theta1)
        if 0 < idx < len(x_data):
            x0, x1 = x_data[idx-1], x_data[idx]
            # ★修正: y_stripped ではなく y_data から補間（累積誤差防止）
            y0, y1 = y_data[idx-1], y_data[idx]
            
            y_interp = y0 + (y1 - y0) * (two_theta1 - x0) / (x1 - x0) if x1 != x0 else y0
            
            y_stripped[i] = y_data[i] - ratio * y_interp
            
    return np.maximum(y_stripped, 0)

def apply_smoothing(y_data, window_size, polyorder=3):
    window_size = int(window_size)

    if window_size < 3:
        return y_data

    # 奇数に揃える
    if window_size % 2 == 0:
        window_size += 1

    # polyorder より大きくなければならない
    if window_size <= polyorder:
        window_size = polyorder + 2
        if window_size % 2 == 0:
            window_size += 1

    # データ長より必ず小さくする
    if window_size >= len(y_data):
        window_size = len(y_data) - 1
        if window_size % 2 == 0:
            window_size -= 1
        if window_size <= polyorder:
            return y_data

    return savgol_filter(y_data, window_length=window_size, polyorder=polyorder, mode='interp')

def calc_ghost_2theta(primary_2theta, lambda_main, lambda_ghost):
    """ブラッグの法則に基づき、主波長での2θから副波長(ゴースト)の2θを計算する"""
    theta_rad = math.radians(primary_2theta / 2.0)
    sin_theta = math.sin(theta_rad)
    
    sin_theta_ghost = (lambda_ghost / lambda_main) * sin_theta
    
    if sin_theta_ghost > 1.0 or sin_theta_ghost < -1.0:
        return None
        
    theta_ghost_rad = math.asin(sin_theta_ghost)
    return math.degrees(theta_ghost_rad) * 2.0

def perform_peak_search(x_data, y_data, nsmooth=5, threshold=100.0):
    """
    XRDデータからピークを自動検索する。

    Args:
        x_data (np.ndarray): 2θ角度の配列
        y_data (np.ndarray): 強度の配列
        nsmooth (int): スムージング窓サイズ
        threshold (float): ピーク強度の最小しきい値

    Returns:
        list: [{"x": 角度, "intensity": 強度}, ...] の形式のリスト
    """
    if len(x_data) < 5 or len(y_data) < 5:
        return []

    y_smooth = apply_smoothing(y_data, window_size=int(nsmooth))

    # データ点数の0.5%を最小ピーク間距離とする
    min_distance = max(1, int(len(x_data) * 0.005))

    peak_indices, _ = find_peaks(
        y_smooth,
        height=threshold,
        distance=min_distance
    )

    return [
        {"x": float(x_data[i]), "intensity": float(y_data[i])}
        for i in peak_indices
    ]