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

# 输入实验数据 (电压和电流)
voltage = np.array([0.50, 1.00, 1.50, 2.00, 2.50, 3.00, 3.50, 4.00, 4.50, 5.00, 5.50, 6.00, 6.50, 7.00, 7.50, 8.00, 8.50, 9.00,
                    9.50, 10.00, 10.50, 11.00, 11.50, 12.00, 12.50, 13.00, 13.50, 14.00, 14.50, 15.00, 15.50, 16.00, 16.50,
                    17.00, 17.50, 18.00, 18.50, 19.00, 19.50, 20.00, 20.50, 21.00, 21.50, 22.00, 22.50, 23.00, 23.50, 24.00,
                    24.50, 25.00, 25.50, 26.00, 26.50, 27.00, 27.50, 28.00, 28.50, 29.00, 29.50, 30.00, 30.50, 31.00, 31.50,
                    32.00, 32.50, 33.00, 33.50, 34.00, 34.50, 35.00, 35.50, 36.00, 36.50, 37.00, 37.50, 38.00, 38.50, 39.00,
                    39.50, 40.00, 40.50, 41.00, 41.50, 42.00, 42.50, 43.00, 43.50, 44.00, 44.50, 45.00, 45.50, 46.00, 46.50,
                    47.00, 47.50, 48.00, 48.50, 49.00, 49.50, 50.00, 50.50, 51.00, 51.50, 52.00, 52.50, 53.00, 53.50, 54.00,
                    54.50, 55.00, 55.50, 56.00, 56.50, 57.00, 57.50, 58.00, 58.50, 59.00, 59.50, 60.00, 60.50, 61.00, 61.50,
                    62.00, 62.50, 63.00, 63.50, 64.00, 64.50, 65.00, 65.50, 66.00, 66.50, 67.00, 67.50, 68.00, 68.50, 69.00,
                    69.50, 70.00, 70.50, 71.00, 71.50, 72.00, 72.50, 73.00, 73.50, 74.00, 74.50, 75.00, 75.50, 76.00, 76.50,
                    77.00, 77.50, 78.00, 78.50])
current = np.array([0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.30, 2.30,
                    7.70, 13.90, 23.50, 31.30, 37.50, 42.10, 48.50, 53.20, 55.80, 59.80, 62.50, 64.50, 66.50, 67.80, 67.50,
                    66.30, 63.70, 60.10, 54.70, 48.30, 42.30, 34.80, 27.70, 21.80, 18.40, 17.80, 24.40, 35.80, 52.60, 68.90,
                    80.40, 92.70, 102.10, 107.50, 112.50, 114.10, 113.20, 109.10, 101.40, 92.00, 77.90, 62.00, 45.70, 32.50,
                    20.10, 11.70, 8.90, 14.30, 38.20, 56.20, 80.10, 101.00, 114.80, 129.30, 141.40, 148.70, 154.80, 157.00,
                    155.60, 150.40, 141.30, 130.10, 113.20, 93.30, 71.70, 53.60, 34.50, 19.70, 13.30, 18.50, 37.60, 57.90,
                    83.70, 106.80, 124.80, 145.00, 162.00, 173.40, 183.20, 190.30, 193.10, 192.50, 187.50, 178.10, 166.70,
                    149.60, 129.20, 110.30, 87.10, 64.90, 49.40, 39.20, 43.50, 56.50, 77.70, 100.50, 119.40, 141.80, 161.40,
                    176.60, 193.80, 208.10, 218.60, 224.10, 228.60, 228.60, 225.30, 217.50, 205.30, 192.20, 173.30, 152.20,
                    133.40, 112.60, 96.90, 91.00, 92.60, 101.40, 113.40, 131.20, 150.20, 166.70, 186.40, 205.80, 222.70,
                    235.70, 249.60, 260.60, 267.00, 271.70, 272.40, 269.90, 263.30, 253.00, 241.70, 225.50, 208.00, 193.90,
                    179.70, 170.30])

# 找到波峰和波谷
peaks, _ = find_peaks(current, distance=10)
valleys, _ = find_peaks(-current, distance=10)

# 平衡点（波峰和波谷中点）
mid_points = [(peaks[i] + valleys[i]) // 2 for i in range(min(len(peaks), len(valleys)))]

# 提取波峰、波谷和平衡点的 x 和 y 值
peak_voltage = voltage[peaks]
peak_current = current[peaks]

valley_voltage = voltage[valleys]
valley_current = current[valleys]

mid_voltage = voltage[mid_points]
mid_current = current[mid_points]

# 打印每个波峰的 xy 值
print("Peaks:")
for v, c in zip(peak_voltage, peak_current):
    print(f"Voltage: {v}, Current: {c}")

# 打印每个波谷的 xy 值
print("\nValleys:")
for v, c in zip(valley_voltage, valley_current):
    print(f"Voltage: {v}, Current: {c}")

# 打印每个中点的 xy 值
print("\nMidpoints:")
for v, c in zip(mid_voltage, mid_current):
    print(f"Voltage: {v}, Current: {c}")

# 绘图
plt.figure(figsize=(10, 6))
plt.plot(voltage, current, label="Original Data", color='blue', marker='o', markersize=4, linestyle='-')

# 标注波峰
plt.plot(peak_voltage, peak_current, "r^", label="Peaks")
# 标注波谷
plt.plot(valley_voltage, valley_current, "gv", label="Valleys")
# 标注平衡点
plt.plot(mid_voltage, mid_current, "ko", label="Midpoints")

plt.title("Voltage vs Current")
plt.xlabel("Voltage (V)")
plt.ylabel("Current (A)")
plt.legend()
plt.grid(True)
plt.savefig('frank_hertz_plot.png')