import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import savgol_filter
import matplotlib

# グローバル設定
infile1 = "Eagle基板.TXT"
infile2 = "試料1_a-GeO2_eagle.TXT"
#infile2 = "試料2_a-GeO2_cAl2O3.TXT"

nsmooth = 15     # 平滑化点数（奇数）
norder_smooth = 2
xrange_analysis = [10.0, 50.0]

norder_bg = 3      # 多項式次数
bg_ranges = [(10, 11), (40, 50)]

# フォント設定（日本語対応）
matplotlib.rcParams['font.family'] = 'MS Gothic'


def read_xrd_data(filepath, xrange_analysis):
    x, y = [], []
    with open(filepath, 'r', encoding='utf-8') as f:
        for line in f:
            if line.strip():
                parts = line.strip().split('\t')
                if len(parts) == 2:
                    try:
                        x_val = float(parts[0])
                        y_val = float(parts[1])
                        if x_val < xrange_analysis[0] or x_val > xrange_analysis[1]:
                            continue

                        x.append(x_val)
                        y.append(y_val)
                    except ValueError:
                        continue

    return np.array(x), np.array(y)

def fit_background(x_bg, y_bg, x_cal, norder = 2):
    print()
    print("BG data:")
    for x, y in zip(x_bg, y_bg):
        print(f"   {x:6.3g}\t{y:6.3g}")

    coeffs = np.polyfit(x_bg, y_bg, deg=norder)

    return np.polyval(coeffs, x_cal)

def extract(x, y, ranges):
    """
    指定された範囲のデータを抽出して返す。
    x, y: 元のスペクトルデータ
    ranges: [(xmin1, xmax1), (xmin2, xmax2), ...] の形式
    戻り値: x_fit, y_fit（背景フィット用データ）
    """

    mask = np.zeros_like(x, dtype=bool)
    for r in ranges:
        mask |= (x >= r[0]) & (x <= r[1])

    return x[mask], y[mask]

def normalize_to_max(y, target_max):
    return y / np.max(y) * target_max

# データ読み込み
x1, y1_raw = read_xrd_data(infile1, xrange_analysis)
x2, y2_raw = read_xrd_data(infile2, xrange_analysis)

# 平滑化
y1 = savgol_filter(y1_raw, nsmooth, norder_smooth)
y2 = savgol_filter(y2_raw, nsmooth, norder_smooth)

# バックグラウンドフィット
x1fit, y1fit = extract(x1, y1, bg_ranges)
x2fit, y2fit = extract(x2, y2, bg_ranges)
bg1 = fit_background(x1fit, y1fit, x1, norder = norder_bg)
bg2 = fit_background(x2fit, y2fit, x2, norder = norder_bg)

# バックグラウンド除去
y1_corr = y1 - bg1
y2_corr = y2 - bg2
#print("x1=", len(x1), x1)
#print("y1=", len(y1_corr), y1_corr)

# xmaxの決定（infile1の補正後スペクトルで最大値の位置）
ymax1 = np.max(y1_corr)
xmax_index = np.argmax(y1_corr)
xmax = x1[xmax_index]
print(f"Max data in infile1 at ix={xmax_index} ({xmax}, {ymax1})")

# infile2のxmaxに最も近い点を探す
x2_index = np.argmin(np.abs(x2 - xmax))
y2_at_xmax = y2_corr[x2_index]
print(f"Data at ix={xmax_index} in infile2: i={x2_index}: ({x2[x2_index]}, {y2_at_xmax})")

# 規格化（infile2のxmax点の強度で両者をスケーリング）
y1_norm = y1_corr
y2_norm = y2_corr * ymax1 / y2_at_xmax

# グラフ①：元スペクトルとバックグラウンド
plt.figure(figsize=(10, 6))
plt.plot(x1, y1, label="Eagle基板（平滑化）", linestyle='-')
plt.plot(x1, bg1, label="Eagle基板 BG", linestyle='--')
plt.plot(x2, y2, label="試料1 a-GeO₂（平滑化）", linestyle='-')
plt.plot(x2, bg2, label="試料1 BG", linestyle='--')
plt.xlabel("2θ (deg)")
plt.ylabel("Intensity (a.u.)")
plt.title("XRDスペクトルとバックグラウンド（平滑化済）")
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.savefig('rawdata.png')
plt.pause(0.1)

# グラフ②：バックグラウンド除去・規格化後
plt.figure(figsize=(10, 6))
plt.plot(x1, y1_norm, label="Eagle基板（規格化）", linestyle='--')
plt.plot(x2, y2_norm, label="試料1 a-GeO₂（規格化）", linewidth=2)
plt.xlabel("2θ (deg)")
plt.ylabel("Intensity (a.u.)")
plt.title("XRDスペクトル比較（バックグラウンド除去・規格化済）")
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.savefig('bgsub.png')
plt.pause(0.1)

# 差スペクトル
plt.figure(figsize=(10, 6))
plt.plot(x1, y2_norm - y1_norm, label="差スペクトル", linestyle='--')
plt.xlabel("2θ (deg)")
plt.ylabel("delta I (a.u.)")
plt.title("差スペクトル")
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.savefig('difference.png')
plt.show()