import numpy as np
import math
from core.constants import WAVE_MAP

def estimate_lattice_from_peaks(positions_2theta, hkls_dicts, lambda_A):
        """
        ピーク位置(2θ)とミラー指数(hkl)から、最小二乗法を用いて一般結晶系の格子定数(a,b,c,α,β,γ)を逆算します。
        十分な独立したピーク(最低6つ)がない場合や、特異行列で計算が収束しない場合は None を返します。
        """
        if len(positions_2theta) < 6 or not hkls_dicts:
            return None
            
        A = []
        Q = []
        for th_deg, hkl in zip(positions_2theta, hkls_dicts):
            th = math.radians(th_deg / 2.0)
            if th <= 0 or th >= math.pi/2: continue
            
            # d^* = 2 * sin(θ) / λ
            d_star = (2.0 * math.sin(th)) / lambda_A
            q = d_star ** 2
            
            h = hkl.get('h', 0)
            k = hkl.get('k', 0)
            l = hkl.get('l', 0)
            
            A.append([h**2, k**2, l**2, 2*h*k, 2*k*l, 2*h*l])
            Q.append(q)
            
        if len(A) < 6:
            return None
            
        A = np.array(A)
        Q = np.array(Q)
        
        try:
            # 最小二乗法で A * x = Q を解く
            res = np.linalg.lstsq(A, Q, rcond=None)
            x = res[0]
            
            astar2, bstar2, cstar2, ab_cos_gastar, bc_cos_alstar, ca_cos_bestar = x
            
            # 物理的にあり得ない解（負の平方根）が出た場合は収束失敗とみなす
            if astar2 <= 0 or bstar2 <= 0 or cstar2 <= 0:
                print("デバッグ: 負の平方根が発生したため、格子定数の推定を中止しました。")
                return None
                
            astar = math.sqrt(astar2)
            bstar = math.sqrt(bstar2)
            cstar = math.sqrt(cstar2)
            
            cos_gastar = max(-1.0, min(1.0, ab_cos_gastar / (astar * bstar)))
            cos_alstar = max(-1.0, min(1.0, bc_cos_alstar / (bstar * cstar)))
            cos_bestar = max(-1.0, min(1.0, ca_cos_bestar / (cstar * astar)))
            
            sin_gastar = math.sqrt(1 - cos_gastar**2)
            sin_alstar = math.sqrt(1 - cos_alstar**2)
            sin_bestar = math.sqrt(1 - cos_bestar**2)
            
            # 逆格子の体積 V*
            vol_term = 1 - cos_alstar**2 - cos_bestar**2 - cos_gastar**2 + 2*cos_alstar*cos_bestar*cos_gastar
            if vol_term <= 0: return None
            Vstar = astar * bstar * cstar * math.sqrt(vol_term)
            
            if Vstar <= 0: return None
            V = 1.0 / Vstar
            
            # 実格子定数の計算
            a = bstar * cstar * sin_alstar * V
            b = cstar * astar * sin_bestar * V
            c = astar * bstar * sin_gastar * V
            
            cos_al = max(-1.0, min(1.0, (cos_bestar * cos_gastar - cos_alstar) / (sin_bestar * sin_gastar)))
            cos_be = max(-1.0, min(1.0, (cos_gastar * cos_alstar - cos_bestar) / (sin_gastar * sin_alstar)))
            cos_ga = max(-1.0, min(1.0, (cos_alstar * cos_bestar - cos_gastar) / (sin_alstar * sin_bestar)))
            
            alpha = math.degrees(math.acos(cos_al))
            beta = math.degrees(math.acos(cos_be))
            gamma = math.degrees(math.acos(cos_ga))
            
            # 極端な値（計算崩壊）を弾く
            if any(val > 1000 or val < 0.1 for val in [a, b, c]):
                print(f"デバッグ: 格子定数の推定が異常値(a={a:.1f}, b={b:.1f}, c={c:.1f})となったため破棄しました。")
                return None
                
            # ★★★ 成功時のデバッグ出力を追加 ★★★
            print(f"[格子定数推定 成功] 利用ピーク数: {len(Q)}")
            print(f"  a = {a:.4f} Å, b = {b:.4f} Å, c = {c:.4f} Å")
            print(f"  α = {alpha:.2f}°, β = {beta:.2f}°, γ = {gamma:.2f}°")
            
            return {"a": a, "b": b, "c": c, "alpha": alpha, "beta": beta, "gamma": gamma}
            
        except Exception as e:
            print(f"デバッグ: 格子定数の推定計算が収束しませんでした。詳細: {e}")
            return None
        
def _gixrd_oop_best(h, k, l,
                        a=None, b=None, c=None,
                        alpha=90.0, beta=90.0, gamma=90.0,
                        omega_deg=0.5, tth_min=3.0, tth_max=100.0,
                        lambda_A=1.5418,
                        du=0.01, window=0.02,
                        out_h=1, out_k=0, out_l=0):
        """
        Grazing Incidence (GI) X線回折の Out-of-plane 配置において、
        指定された反射 (h,k,l) がエワルド球を横切る際の「最近接距離」と「その時の条件」を計算します。
        
        【概要】
        任意の結晶系 (a,b,c, α,β,γ) に対応し、実格子ベクトルから逆格子ベクトル G を構築します。
        面外配向軸 N に基づき、G を面内成分 g_par と面外成分 g_per に分解し、
        入射角 ω と散乱角 2θ の走査範囲内で、回折ベクトル Q と G の距離が
        最小となる出射角の正弦 u = sin(θ_f) を探索します。
        
        戻り値: (d_par, d_per, g_par, g_per, best_u) または計算不可の場合は None
        """
        try:
            if any(val is None or val <= 0 for val in [a, b, c, lambda_A]):
                return None

            h = float(h); k = float(k); l = float(l)
            a = float(a); b = float(b); c = float(c)
            out_h = float(out_h); out_k = float(out_k); out_l = float(out_l)
            
            al = math.radians(float(alpha))
            be = math.radians(float(beta))
            ga = math.radians(float(gamma))

            # 実格子ベクトルの構築 (a軸をX軸, b軸をXY平面に置く標準配置)
            ax, ay, az = a, 0.0, 0.0
            bx, by, bz = b * math.cos(ga), b * math.sin(ga), 0.0
            cx = c * math.cos(be)
            cy = c * (math.cos(al) - math.cos(be)*math.cos(ga)) / math.sin(ga)
            cz = math.sqrt(max(0.0, c**2 - cx**2 - cy**2))
            
            v_a = np.array([ax, ay, az])
            v_b = np.array([bx, by, bz])
            v_c = np.array([cx, cy, cz])
            
            V = np.dot(v_a, np.cross(v_b, v_c))
            if V <= 0: return None
            
            # 逆格子ベクトルの構築 (GIXRDの散乱ベクトル定義に合わせ 2π を含む)
            star_a = 2.0 * math.pi * np.cross(v_b, v_c) / V
            star_b = 2.0 * math.pi * np.cross(v_c, v_a) / V
            star_c = 2.0 * math.pi * np.cross(v_a, v_b) / V

            def gvec(hh, kk, ll):
                return hh * star_a + kk * star_b + ll * star_c

            g = gvec(h, k, l)
            n = gvec(out_h, out_k, out_l)
            n_norm = float(np.linalg.norm(n))
            if n_norm < 1e-12: return None
            n_hat = n / n_norm

            g_par = abs(float(np.dot(g, n_hat)))
            g2 = float(np.dot(g, g))
            g_per = math.sqrt(max(0.0, g2 - g_par*g_par))

            k0 = 2.0 * math.pi / lambda_A
            om = math.radians(omega_deg)
            s = math.sin(om)
            co = math.cos(om)

            tf_min = math.radians(max(0.0, tth_min - omega_deg))
            tf_max = math.radians(min(179.9, tth_max - omega_deg))
            u_min = max(-1.0, min(1.0, math.sin(tf_min)))
            u_max = max(-1.0, min(1.0, math.sin(tf_max)))
            if u_min > u_max: u_min, u_max = u_max, u_min

            u0 = (g_par / k0) - s
            def clamp_u(u): return min(u_max, max(u_min, u))
            uc = clamp_u(u0)

            nstep = int(math.ceil(max(0.0, float(window)) / max(1e-5, float(du))))
            candidates = {uc, u_min, u_max}
            for i in range(1, nstep + 1):
                candidates.add(clamp_u(uc + i*du))
                candidates.add(clamp_u(uc - i*du))

            best = float("inf")
            best_tuple = None

            for u in candidates:
                Q_par = k0 * (s + u)
                root = math.sqrt(max(0.0, 1.0 - u*u))
                Q_per = k0 * abs(root - co)
                d_par = (g_par - Q_par)
                d_per = (g_per - Q_per)
                d = math.sqrt(d_par*d_par + d_per*d_per)
                if d < best:
                    best = d
                    best_tuple = (d_par, d_per, g_par, g_per, u)

            return best_tuple
        except Exception:
            return None
        
def gixrd_oop_dist(h, k, l, out_h=1, out_k=0, out_l=0, **kwargs):
        """
        GI-out-of-plane のエワルド球と逆格子点との最近接距離 [Å^-1] を計算します。
        
        引数:
            out_h, out_k, out_l (float): 面外配向している結晶面のミラー指数（デフォルトは 1, 0, 0）
        """
        best = _gixrd_oop_best(h, k, l, out_h=out_h, out_k=out_k, out_l=out_l, **kwargs)
        if best is None:
            return float("inf")
        d_par, d_per, *_ = best
        return float(math.sqrt(d_par*d_par + d_per*d_per))

def _gixrd_oop_ndist(h, k, l,
                         fwhm_out_deg=0.05, fwhm_in_deg=0.50,
                         out_h=1, out_k=0, out_l=0,
                         **kwargs):
        """
        GI-out-of-plane配置における、ピークの「異方的な広がり」を考慮した正規化距離(σ単位)を計算します。
        
        引数:
            fwhm_out_deg (float): 面外方向(Out-of-plane)のロッキングカーブ半値幅 [deg]
            fwhm_in_deg (float): 面内方向(In-plane)の広がり半値幅 [deg]
            out_h, out_k, out_l (float): 面外配向している結晶面のミラー指数（デフォルトは 1, 0, 0）
        """
        best = _gixrd_oop_best(h, k, l, out_h=out_h, out_k=out_k, out_l=out_l, **kwargs)
        if best is None:
            return float("inf")

        d_par, d_per, g_par, g_per, _u = best

        # FWHM[deg] -> σ[rad] への変換 (ガウス関数近似: FWHM ≈ 2.355 * σ)
        sig_out = math.radians(float(fwhm_out_deg)) / 2.355
        sig_in  = math.radians(float(fwhm_in_deg))  / 2.355

        eps = 1e-12
        sigma_par = math.sqrt((g_per*sig_out)**2 + eps)
        sigma_per = math.sqrt((g_par*sig_out)**2 + (g_per*sig_in)**2 + eps)

        dn = math.sqrt((d_par/sigma_par)**2 + (d_per/sigma_per)**2)
        return float(dn)

def gixrd_oop_vis_fwhm(h, k, l,
                            nsigma=2.0,
                            fwhm_out_deg=0.05, fwhm_in_deg=0.50,
                            out_h=1, out_k=0, out_l=0,
                            **kwargs):
        """
        GI-out-of-plane配置において、ピークが走査線にかかって「観測可能(Visible)」かどうかを判定します。
        
        引数:
            nsigma (float): 観測可能と判定する閾値（デフォルトは 2.0σ 以内）
            out_h, out_k, out_l (float): 面外配向している結晶面のミラー指数（デフォルトは 1, 0, 0）
        """
        return bool(_gixrd_oop_ndist(
            h, k, l,
            fwhm_out_deg=fwhm_out_deg,
            fwhm_in_deg=fwhm_in_deg,
            out_h=out_h, out_k=out_k, out_l=out_l,
            **kwargs
        ) <= float(nsigma))