import numpy as np
import matplotlib
matplotlib.use('Agg')  # Non-interactive backend for saving plots
import matplotlib.pyplot as plt
from comtrade import Comtrade
import math


class DistanceRelay:
    """
    Algorytm zabezpieczenia odległościowego
    Implementacja charakterystyki wielokątnej (quadrilateral)
    """
    def __init__(self, Z_line_R=2.0, Z_line_X=8.0, line_angle=75.0):
        # Impedancja linii (obliczona z danych)
        self.Z_line_R = Z_line_R
        self.Z_line_X = Z_line_X
        self.Z_line_mag = np.sqrt(Z_line_R**2 + Z_line_X**2)
        self.line_angle = line_angle

        # === Nastawy stref jako % impedancji linii ===
        # Strefa 1 - 80% linii (natychmiastowa)
        self.Z1_R = self.Z_line_R * 0.8
        self.Z1_X = self.Z_line_X * 0.8
        self.tZ1 = 0       # Brak opóźnienia

        # Strefa 2 - 120% linii (koordynacja)
        self.Z2_R = self.Z_line_R * 1.2
        self.Z2_X = self.Z_line_X * 1.2
        self.tZ2 = 300     # 300ms

        # Strefa 3 - 200% linii (rezerwowa)
        self.Z3_R = self.Z_line_R * 2.0
        self.Z3_X = self.Z_line_X * 2.0
        self.tZ3 = 600     # 600ms

        # Kąt charakterystyki (na podstawie kąta linii)
        self.angle_r1 = line_angle

        # Minimalny prąd i napięcie (zabezpieczenie przed szumem)
        self.I_min = 0.5    # A
        self.U_min = 1.0    # V

        print(f"Nastawy stref (na podstawie Z linii = {self.Z_line_mag:.2f} Ohm, {line_angle:.1f} deg):")
        print(f"  Strefa 1: R={self.Z1_R:.2f} Ohm, X={self.Z1_X:.2f} Ohm (natychmiast)")
        print(f"  Strefa 2: R={self.Z2_R:.2f} Ohm, X={self.Z2_X:.2f} Ohm (300ms)")
        print(f"  Strefa 3: R={self.Z3_R:.2f} Ohm, X={self.Z3_X:.2f} Ohm (600ms)")

        # Stany wewnętrzne dla każdej fazy
        self.init_state()

    def init_state(self):
        """Inicjalizacja stanów dla każdej fazy"""
        # Timery dla każdej fazy i strefy
        self.timers = {
            'L1_Z1': 0, 'L1_Z2': 0, 'L1_Z3': 0,
            'L2_Z1': 0, 'L2_Z2': 0, 'L2_Z3': 0,
            'L3_Z1': 0, 'L3_Z2': 0, 'L3_Z3': 0,
        }
        # Flagi trip
        self.tripped = {'L1': False, 'L2': False, 'L3': False}

    def init_relay(self):
        print("Zabezpieczenie odległościowe zainicjalizowane")
        self.init_state()

    def in_polygon(self, R, X, reach_R, reach_X, angle_deg):
        """
        Sprawdza czy punkt (R, X) jest wewnątrz wielokąta
        Charakterystyka czworokątna (quadrilateral)
        """
        # Obrót punktu o -angle_deg aby wyprostować charakterystykę
        angle_rad = math.radians(angle_deg)
        cos_a = math.cos(-angle_rad)
        sin_a = math.sin(-angle_rad)

        R_rot = R * cos_a - X * sin_a
        X_rot = R * sin_a + X * cos_a

        # Sprawdź czy punkt jest wewnątrz prostokąta w układzie obróconym
        # R musi być dodatnie (kierunek forward)
        # X musi być w zakresie [-reach_X, reach_X]
        # R musi być mniejsze niż reach_R

        # Dodatkowo: nachylone linie R
        R_max = reach_R
        X_max = reach_X
        R_min = 0.1  # Minimalna wartość R (strefa aktywna)

        # Sprawdzenie podstawowych granic
        if R_rot < R_min or R_rot > R_max:
            return False
        if abs(X_rot) > X_max:
            return False

        # Sprawdzenie linii nachylonych (opcjonalnie)
        # Górna granica X
        X_upper = X_max * (1 - (R_rot / R_max) * math.tan(math.radians(90 - angle_deg + 10)))
        # Dolna granica X
        X_lower = -X_max * (1 - (R_rot / R_max) * math.tan(math.radians(90 - angle_deg + 10)))

        return True

    def check_direction(self, U1_zg_re, U1_zg_im, I1_zg_re, I1_zg_im):
        """
        Określenie kierunku na podstawie mocy
        P = Re(U * conj(I)) > 0 = forward
        """
        power = U1_zg_re * I1_zg_re + U1_zg_im * I1_zg_im
        return power > 0  # True = forward

    def step_relay(self, phase, u_re, u_im, i_re, i_im,
                   u0_re, u0_im, i0_re, i0_im,
                   u_zg_re, u_zg_im, i_zg_re, i_zg_im):
        """
        Krok algorytmu dla jednej fazy
        phase: 'L1', 'L2' lub 'L3'
        """
        # Oblicz moduł prądu
        i_mag = math.sqrt(i_re**2 + i_im**2)
        u_mag = math.sqrt(u_re**2 + u_im**2)

        # Sprawdź progi minimalne
        if i_mag < self.I_min or u_mag < self.U_min:
            return 0

        # Oblicz impedancję Z = U / I
        i_mag_sq = i_re**2 + i_im**2
        if i_mag_sq < 1e-9:
            return 0

        z_re = (u_re * i_re + u_im * i_im) / i_mag_sq
        z_x = (u_im * i_re - u_re * i_im) / i_mag_sq

        # Sprawdź kierunek (używamy składowej zgodnej)
        if i_zg_re is not None and i_zg_im is not None:
            forward = True  # Uproszczone - zakładamy forward
        else:
            forward = True

        # Jeśli już wyłączone, nie sprawdzaj dalej
        if self.tripped[phase]:
            return 1

        trip = 0

        if forward:
            # === Strefa 1 - natychmiastowa ===
            if self.in_polygon(z_re, z_x, self.Z1_R, self.Z1_X, self.angle_r1):
                self.tripped[phase] = True
                return 1

            # === Strefa 2 - opóźniona ===
            key_Z2 = f'{phase}_Z2'
            if self.in_polygon(z_re, z_x, self.Z2_R, self.Z2_X, self.angle_r1):
                if self.timers[key_Z2] < self.tZ2:
                    self.timers[key_Z2] += 1
                elif not self.tripped[phase]:
                    self.tripped[phase] = True
                    return 1
            else:
                self.timers[key_Z2] = 0

            # === Strefa 3 - rezerwowa ===
            key_Z3 = f'{phase}_Z3'
            if self.in_polygon(z_re, z_x, self.Z3_R, self.Z3_X, self.angle_r1):
                if self.timers[key_Z3] < self.tZ3:
                    self.timers[key_Z3] += 1
                elif not self.tripped[phase]:
                    self.tripped[phase] = True
                    return 1
            else:
                self.timers[key_Z3] = 0

        return 0

    def reset(self):
        """Reset stanów dla nowego uruchomienia"""
        self.init_state()


# Wczytanie pliku COMTRADE
import sys

# Obsluga argumentow wiersza polecen
if len(sys.argv) > 1:
    base_name = sys.argv[1]  # Nazwa pliku bez rozszerzenia
else:
    base_name = "zwarcie_testowe"  # Domyslna nazwa

cfg_file = f"{base_name}.cfg"
dat_file = f"{base_name}.dat"

print(f"Wczytywanie rejestracji: {cfg_file}")

rec = Comtrade()

# Probuj rozne kodowania dla plikow CFG
encodings = ['utf-8', 'cp1250', 'cp1252', 'latin-1', 'iso-8859-1']
loaded = False

for encoding in encodings:
    try:
        rec.load(cfg_file, dat_file, encoding=encoding)
        loaded = True
        break
    except (UnicodeDecodeError, Exception) as e:
        continue

if not loaded:
    print(f"BŁĄD: Nie można wczytać pliku {cfg_file}")
    sys.exit(1)

# Sprawdz czy mamy wystarczajaco kanalow analogowych
num_analog = len(rec.analog)
print(f"Liczba kanalow analogowych: {num_analog}")

if num_analog < 6:
    print(f"BŁĄD: Za mało kanałów analogowych (wymagane minimum 6)")
    sys.exit(1)

# Wyciagniecie danych (przyjmujemy indeksy kanałów analogowych)
t = np.array(rec.time)
Fs = rec.cfg.sample_rates[0][0] # Czestotliwosc probkowania
f_nom = 50.0 # Czestotliwosc sieci
N = int(Fs / f_nom) # Liczba probek na okres

if N < 1:
    N = 1  # Zabezpieczenie

print(f" Czestotliwosc probkowania: {Fs} Hz")
print(f" Liczba probek na okres: {N}")

# Zabezpieczenie przed index error
max_analog = min(6, num_analog)

# Pobierz tyle ile mamy
i1_raw = np.array(rec.analog[0])
i2_raw = np.array(rec.analog[1]) # Kanał 1: I_L2
i3_raw = np.array(rec.analog[2]) # Kanał 2: I_L3
u1_raw = np.array(rec.analog[3]) # Kanał 3: U_L1
u2_raw = np.array(rec.analog[4]) # Kanał 4: U_L2
u3_raw = np.array(rec.analog[5]) # Kanał 5: U_L3
# (musisz dopasować indeksy analog[] do swojego pliku .cfg)

# === ANALIZA DANYCH - wyznaczenie impedancji przed zwarciem ===
print(f" Czestotliwosc probkowania: {Fs} Hz")
print(f" Liczba probek na okres: {N}")
print(f" Liczba probek: {len(t)}")
print(f" Czas trwania: {t[-1]:.3f} s")

# 3. Definicja filtru Full-Cycle DFT (FCDFT)
def fcdft(samples):
    """Oblicza ortogonalne (Re, Im) dla okna N probek z wykorzystaniem transformaty Fouriera"""
    if len(samples) < N: return 0.0, 0.0
    k = np.arange(N)
    cos_wave = np.cos(2 * np.pi * k / N)
    sin_wave = np.sin(2 * np.pi * k / N)

    re = (2.0 / N) * np.sum(samples * cos_wave)
    im = -(2.0 / N) * np.sum(samples * sin_wave)
    return re, im

# Funkcja pomocnicza do obliczania impedancji
def calculate_impedance_from_raw(u_raw, i_raw, idx):
    """Oblicza impedancję dla próbki idx"""
    if idx < N:
        return 0, 0
    window_i = i_raw[idx-N:idx]
    window_u = u_raw[idx-N:idx]
    i_re, i_im = fcdft(window_i)
    u_re, u_im = fcdft(window_u)
    i_mag_sq = i_re**2 + i_im**2
    if i_mag_sq < 1e-9:
        return 0, 0
    z_re = (u_re * i_re + u_im * i_im) / i_mag_sq
    z_x = (u_im * i_re - u_re * i_im) / i_mag_sq
    return z_re, z_x

# Oblicz impedancję przed zwarciem (próbki 10-100)
pre_fault_start = 10
pre_fault_end = 100

# DFT dla przykładowych próbek
def calculate_impedance_from_raw(u_raw, i_raw, idx):
    """Oblicza impedancję dla próbki idx"""
    if idx < N:
        return 0, 0
    window_i = i_raw[idx-N:idx]
    window_u = u_raw[idx-N:idx]
    i_re, i_im = fcdft(window_i)
    u_re, u_im = fcdft(window_u)
    i_mag_sq = i_re**2 + i_im**2
    if i_mag_sq < 1e-9:
        return 0, 0
    z_re = (u_re * i_re + u_im * i_im) / i_mag_sq
    z_x = (u_im * i_re - u_re * i_im) / i_mag_sq
    return z_re, z_x

# Srednia impedancja przed zwarciem
z_r_list, z_x_list = [], []

# Poprawka: krok = N, ale N nie moze byc 0
step = max(N, 1)
for idx in range(pre_fault_start, min(pre_fault_end, len(t)), step):
    z_r, z_x = calculate_impedance_from_raw(u1_raw, i1_raw, idx)
    if z_r > 0:
        z_r_list.append(z_r)
        z_x_list.append(z_x)

if z_r_list:
    Z_line_R = np.median(z_r_list)
    Z_line_X = np.median(z_x_list)
    Z_line_mag = np.sqrt(Z_line_R**2 + Z_line_X**2)
    print(f"\nImpedancja linii (przed zwarciem):")
    print(f"  R = {Z_line_R:.2f} Ohm")
    print(f"  X = {Z_line_X:.2f} Ohm")
    print(f"  |Z| = {Z_line_mag:.2f} Ohm")

    # Kat linii
    line_angle = np.degrees(np.arctan2(Z_line_X, Z_line_R))
    print(f"  Kat = {line_angle:.1f} deg")
else:
    # Wartosci domyslne jesli nie mozna obliczyc
    Z_line_R = 2.0
    Z_line_X = 8.0
    Z_line_mag = np.sqrt(Z_line_R**2 + Z_line_X**2)
    line_angle = 75.0
    print("Nie mozna wyznaczyc impedancji linii, uzywam wartosci domyslnych")

# Utworzenie relay z automatycznie dobranymi nastawami
relay = DistanceRelay(Z_line_R=Z_line_R, Z_line_X=Z_line_X, line_angle=line_angle)

# Macierz operatora obrotu dla składowych symetrycznych
a = complex(-0.5, np.sqrt(3)/2)
a2 = complex(-0.5, -np.sqrt(3)/2)

# Symulacja "czasu rzeczywistego" próbka po próbce
relay.reset()
relay.init_relay()

trip_history_L1 = []
trip_history_L2 = []
trip_history_L3 = []

# Impedance phasors (R + jX)
z1_r_history = []
z1_x_history = []
z2_r_history = []
z2_x_history = []
z3_r_history = []
z3_x_history = []

def calculate_impedance(u_cpx, i_cpx):
    """Oblicza impedancję Z = U/I jako liczbę zespoloną"""
    i_mag_sq = i_cpx.real**2 + i_cpx.imag**2
    if i_mag_sq < 1e-9:  # Zabezpieczenie przed dzieleniem przez zero
        return 0.0, 0.0
    # Z = U / I = U * conj(I) / |I|^2
    z_cpx = u_cpx * complex(i_cpx.real, -i_cpx.imag) / i_mag_sq
    return z_cpx.real, z_cpx.imag

# Symulacja "czasu rzeczywistego" próbka po próbce
for i in range(N, len(t)):
    # Pobranie okna danych (historyczne N próbek aż do obecnej chwili i)
    window_i1 = i1_raw[i-N:i]
    window_i2 = i2_raw[i-N:i]
    window_i3 = i3_raw[i-N:i]
    window_u1 = u1_raw[i-N:i]
    window_u2 = u2_raw[i-N:i]
    window_u3 = u3_raw[i-N:i]

    # Filtracja DFT dla wszystkich faz
    i1_re, i1_im = fcdft(window_i1)
    i2_re, i2_im = fcdft(window_i2)
    i3_re, i3_im = fcdft(window_i3)
    u1_re, u1_im = fcdft(window_u1)
    u2_re, u2_im = fcdft(window_u2)
    u3_re, u3_im = fcdft(window_u3)

    # Tworzenie liczb zespolonych dla łatwiejszej matematyki symetrycznej
    I1_cpx = complex(i1_re, i1_im)
    I2_cpx = complex(i2_re, i2_im)
    I3_cpx = complex(i3_re, i3_im)
    U1_cpx = complex(u1_re, u1_im)
    U2_cpx = complex(u2_re, u2_im)
    U3_cpx = complex(u3_re, u3_im)

    # Obliczanie składowych symetrycznych
    I0_cpx = (I1_cpx + I2_cpx + I3_cpx) / 3.0
    I1zg_cpx = (I1_cpx + a * I2_cpx + a2 * I3_cpx) / 3.0
    I2pr_cpx = (I1_cpx + a2 * I2_cpx + a * I3_cpx) / 3.0
    U0_cpx = (U1_cpx + U2_cpx + U3_cpx) / 3.0
    U1zg_cpx = (U1_cpx + a * U2_cpx + a2 * U3_cpx) / 3.0
    U2pr_cpx = (U1_cpx + a2 * U2_cpx + a * U3_cpx) / 3.0

    # Przekazanie danych do algorytmu zabezpieczeniowego
    trip_l1 = relay.step_relay('L1',
        U1_cpx.real, U1_cpx.imag, I1_cpx.real, I1_cpx.imag,
        I0_cpx.real, I0_cpx.imag, U0_cpx.real, U0_cpx.imag,
        U1zg_cpx.real, U1zg_cpx.imag, I1zg_cpx.real, I1zg_cpx.imag
    )

    trip_l2 = relay.step_relay('L2',
        U2_cpx.real, U2_cpx.imag, I2_cpx.real, I2_cpx.imag,
        I0_cpx.real, I0_cpx.imag, U0_cpx.real, U0_cpx.imag,
        U1zg_cpx.real, U1zg_cpx.imag, I1zg_cpx.real, I1zg_cpx.imag
    )

    trip_l3 = relay.step_relay('L3',
        U3_cpx.real, U3_cpx.imag, I3_cpx.real, I3_cpx.imag,
        I0_cpx.real, I0_cpx.imag, U0_cpx.real, U0_cpx.imag,
        U1zg_cpx.real, U1zg_cpx.imag, I1zg_cpx.real, I1zg_cpx.imag
    )

    # Rejestracja wyniku
    trip_history_L1.append(trip_l1)
    trip_history_L2.append(trip_l2)
    trip_history_L3.append(trip_l3)

    # Obliczanie impedancji Z = U / I
    z1_r, z1_x = calculate_impedance(U1_cpx, I1_cpx)
    z2_r, z2_x = calculate_impedance(U2_cpx, I2_cpx)
    z3_r, z3_x = calculate_impedance(U3_cpx, I3_cpx)

    z1_r_history.append(z1_r)
    z1_x_history.append(z1_x)
    z2_r_history.append(z2_r)
    z2_x_history.append(z2_x)
    z3_r_history.append(z3_r)
    z3_x_history.append(z3_x)

# 5. Rysowanie wyników
plt.figure(figsize=(16, 12))

# Prądy faz L1, L2, L3
plt.subplot(3, 2, 1)
plt.plot(t[N:], i1_raw[N:], label='I_L1', color='blue')
plt.plot(t[N:], i2_raw[N:], label='I_L2', color='green')
plt.plot(t[N:], i3_raw[N:], label='I_L3', color='orange')
plt.title('Prady faz L1, L2, L3')
plt.xlabel('Czas [s]')
plt.ylabel('Prad [A]')
plt.legend()
plt.grid(True)

# Charakterystyka R-X
plt.subplot(3, 2, 2)
# Rysuj strefy
angle_rad = math.radians(relay.angle_r1)

# Strefa 1
r1_vals = np.linspace(0, relay.Z1_R, 50)
x1_upper = relay.Z1_X * (1 - r1_vals / relay.Z1_R * math.tan(math.radians(15)))
x1_lower = -relay.Z1_X * (1 - r1_vals / relay.Z1_R * math.tan(math.radians(15)))
plt.fill_between(r1_vals, x1_lower, x1_upper, alpha=0.2, color='green', label='Strefa 1')

# Strefa 2
r2_vals = np.linspace(0, relay.Z2_R, 50)
x2_upper = relay.Z2_X * (1 - r2_vals / relay.Z2_R * math.tan(math.radians(15)))
x2_lower = -relay.Z2_X * (1 - r2_vals / relay.Z2_R * math.tan(math.radians(15)))
plt.fill_between(r2_vals, x2_lower, x2_upper, alpha=0.15, color='yellow', label='Strefa 2')

# Strefa 3
r3_vals = np.linspace(0, relay.Z3_R, 50)
x3_upper = relay.Z3_X * (1 - r3_vals / relay.Z3_R * math.tan(math.radians(15)))
x3_lower = -relay.Z3_X * (1 - r3_vals / relay.Z3_R * math.tan(math.radians(15)))
plt.fill_between(r3_vals, x3_lower, x3_upper, alpha=0.1, color='red', label='Strefa 3')

# Linia impedancji linii
z_line = np.linspace(0, relay.Z_line_mag * 1.5, 50)
x_line = z_line * math.tan(angle_rad)
plt.plot(z_line, x_line, 'k--', linewidth=1, label='Linia Z')

# Trajektorie impedancji
# Filtruj tylko próbki podczas zwarcia (użyj indeksu gdzie napięcie spada)
fault_indices = np.where(np.array(z1_x_history) < 1.0)[0]
if len(fault_indices) > 0:
    plt.plot(np.array(z1_r_history)[fault_indices], np.array(z1_x_history)[fault_indices],
             'b-', linewidth=2, label='Z_L1 (zwarcie)')
plt.plot(z1_r_history, z1_x_history, 'b.', markersize=1, alpha=0.5, label='Trajektoria L1')
plt.plot(z2_r_history, z2_x_history, 'g.', markersize=1, alpha=0.5, label='Trajektoria L2')
plt.plot(z3_r_history, z3_x_history, 'o', color='orange', markersize=1, alpha=0.5, label='Trajektoria L3')

plt.xlim(-1, relay.Z3_R * 1.2)
plt.ylim(-relay.Z3_X * 1.2, relay.Z3_X * 1.2)
plt.xlabel('R [Ohm]')
plt.ylabel('X [Ohm]')
plt.title('Charakterystyka R-X zabezpieczenia odleglosciowego')
plt.legend(loc='upper right', fontsize=8)
plt.grid(True)
plt.axhline(y=0, color='k', linewidth=0.5)
plt.axvline(x=0, color='k', linewidth=0.5)

# Napięcia faz L1, L2, L3
plt.subplot(3, 2, 3)
plt.plot(t[N:], u1_raw[N:], label='U_L1', color='blue')
plt.plot(t[N:], u2_raw[N:], label='U_L2', color='green')
plt.plot(t[N:], u3_raw[N:], label='U_L3', color='orange')
plt.title('Napięcia faz L1, L2, L3')
plt.xlabel('Czas [s]')
plt.ylabel('Napięcie [V]')
plt.legend()
plt.grid(True)

# Impedancje faz L1, L2, L3
plt.subplot(3, 2, 4)
plt.plot(t[N:], z1_r_history, label='R_L1', color='blue', linestyle='-')
plt.plot(t[N:], z1_x_history, label='X_L1', color='blue', linestyle='--')
plt.plot(t[N:], z2_r_history, label='R_L2', color='green', linestyle='-')
plt.plot(t[N:], z2_x_history, label='X_L2', color='green', linestyle='--')
plt.plot(t[N:], z3_r_history, label='R_L3', color='orange', linestyle='-')
plt.plot(t[N:], z3_x_history, label='X_L3', color='orange', linestyle='--')
plt.title('Impedancje faz (R + jX)')
plt.xlabel('Czas [s]')
plt.ylabel('Impedancja [Ohm]')
plt.legend()
plt.grid(True)

# Wyjscia zabezpieczenia
plt.subplot(3, 2, 5)
max_val = max(max(i1_raw), max(i2_raw), max(i3_raw))
plt.plot(t[N:], i1_raw[N:], label='I_L1', color='blue', alpha=0.5)
plt.plot(t[N:], i2_raw[N:], label='I_L2', color='green', alpha=0.5)
plt.plot(t[N:], i3_raw[N:], label='I_L3', color='orange', alpha=0.5)
plt.plot(t[N:], np.array(trip_history_L1) * max_val, label='Trip L1-E', color='red', linewidth=2)
plt.plot(t[N:], np.array(trip_history_L2) * max_val, label='Trip L2-E', color='darkred', linewidth=2, linestyle='--')
plt.plot(t[N:], np.array(trip_history_L3) * max_val, label='Trip L3-E', color='darkorange', linewidth=2, linestyle=':')
plt.title('Wynik testu algorytmu ZDistA')
plt.xlabel('Czas [s]')
plt.ylabel('Wartosc')
plt.legend()
plt.grid(True)

# Modul impedancji
plt.subplot(3, 2, 6)
z1_mag = np.sqrt(np.array(z1_r_history)**2 + np.array(z1_x_history)**2)
z2_mag = np.sqrt(np.array(z2_r_history)**2 + np.array(z2_x_history)**2)
z3_mag = np.sqrt(np.array(z3_r_history)**2 + np.array(z3_x_history)**2)
plt.plot(t[N:], z1_mag, label='|Z_L1|', color='blue')
plt.plot(t[N:], z2_mag, label='|Z_L2|', color='green')
plt.plot(t[N:], z3_mag, label='|Z_L3|', color='orange')
z1_reach = np.sqrt(relay.Z1_R**2 + relay.Z1_X**2)
z2_reach = np.sqrt(relay.Z2_R**2 + relay.Z2_X**2)
z3_reach = np.sqrt(relay.Z3_R**2 + relay.Z3_X**2)
plt.axhline(y=z1_reach, color='green', linestyle='--', label='Z1 reach')
plt.axhline(y=z2_reach, color='orange', linestyle='--', label='Z2 reach')
plt.axhline(y=z3_reach, color='red', linestyle='--', label='Z3 reach')
plt.title('Modul impedancji |Z|')
plt.xlabel('Czas [s]')
plt.ylabel('|Z| [Ohm]')
plt.legend()
plt.grid(True)

plt.tight_layout()
plt.savefig('test_result.png')
print("Wynik zapisany do test_result.png")

# Generowanie pliku rezultat.md
def detect_fault_and_generate_report(t, trip_history, z_r_history, z_x_history,
                                     i1_raw, i2_raw, i3_raw,
                                     u1_raw, u2_raw, u3_raw,
                                     relay, base_name):
    """Wykrywa zwarcie i generuje raport"""

    # Znajdz indeks pierwszego trip
    fault_idx = None
    for i, trip in enumerate(trip_history):
        if trip > 0:
            fault_idx = i
            break

    # Przygotuj dane do raportu
    report = []
    report.append(f"# Wynik analizy zabezpieczenia odleglosciowego")
    report.append(f"")
    report.append(f"## Rejestracja: {base_name}")
    report.append(f"")
    report.append(f"## Parametry zabezpieczenia")
    report.append(f"- Impedancja linii: {relay.Z_line_mag:.2f} Ohm")
    report.append(f"- Kat linii: {relay.line_angle:.1f} deg")
    report.append(f"- Strefa 1: R={relay.Z1_R:.2f} Ohm, X={relay.Z1_X:.2f} Ohm (natychmiast)")
    report.append(f"- Strefa 2: R={relay.Z2_R:.2f} Ohm, X={relay.Z2_X:.2f} Ohm (300ms)")
    report.append(f"- Strefa 3: R={relay.Z3_R:.2f} Ohm, X={relay.Z3_X:.2f} Ohm (600ms)")
    report.append(f"")

    if fault_idx is not None:
        # Zwarcie wykryte
        fault_time = t[N + fault_idx]
        report.append(f"## Zwarcie wykryte")
        report.append(f"- Czas zwarcia: {fault_time*1000:.2f} ms")
        report.append(f"")

        # Ktora faza
        if trip_history[fault_idx] > 0:
            report.append(f"- Faza: L1 (trip)")
        report.append(f"")

        # Wartosci podczas zwarcia
        z_r_fault = z_r_history[fault_idx]
        z_x_fault = z_x_history[fault_idx]
        z_mag_fault = np.sqrt(z_r_fault**2 + z_x_fault**2)

        report.append(f"## Wartosci podczas zwarcia")
        report.append(f"- Impedancja Z: {z_mag_fault:.2f} Ohm")
        report.append(f"- R: {z_r_fault:.2f} Ohm")
        report.append(f"- X: {z_x_fault:.2f} Ohm")
        report.append(f"")

        # Prad i napiecie
        i1_fault = i1_raw[N + fault_idx]
        u1_fault = u1_raw[N + fault_idx]
        report.append(f"- Prad L1: {i1_fault:.2f} A")
        report.append(f"- Napiecie L1: {u1_fault:.2f} V")
    else:
        report.append(f"## Brak zwarcia")
        report.append(f"- Nie wykryto zadzialania zabezpieczenia")
        report.append(f"")

    # Podsumowanie
    any_trip = any(t > 0 for t in trip_history)
    report.append(f"## Podsumowanie")
    report.append(f"- Trip L1: {'TAK' if any(trip_history) else 'NIE'}")
    report.append(f"")

    # Zapisz raport
    report_content = "\n".join(report)
    with open('rezultat.md', 'w', encoding='utf-8') as f:
        f.write(report_content)

    print("Raport zapisany do rezultat.md")
    return fault_idx is not None

# Generuj raport
has_fault = detect_fault_and_generate_report(
    t,
    [trip_history_L1[i] + trip_history_L2[i] + trip_history_L3[i] for i in range(len(trip_history_L1))],
    z1_r_history, z1_x_history,
    i1_raw, i2_raw, i3_raw,
    u1_raw, u2_raw, u3_raw,
    relay, base_name
)

# Generowanie pliku rezultat.md
def generate_result_md(t, trip_history_L1, trip_history_L2, trip_history_L3,
                       z1_r_history, z1_x_history, z2_r_history, z2_x_history, z3_r_history, z3_x_history,
                       Z_line_R, Z_line_X, relay):
    """Generuje plik rezultat.md z informacjami o wykryciu zwarcia"""

    md_content = []
    md_content.append("# Wynik analizy zabezpieczenia odleglosciowego")
    md_content.append("")
    md_content.append(f"## Parametry zabezpieczenia")
    md_content.append(f"- Impedancja linii: R={Z_line_R:.2f} Ohm, X={Z_line_X:.2f} Ohm")
    md_content.append(f"- Kat linii: {relay.angle_r1:.1f} st.")
    md_content.append(f"- Strefa 1: R={relay.Z1_R:.2f} Ohm, X={relay.Z1_X:.2f} Ohm (natychmiast)")
    md_content.append(f"- Strefa 2: R={relay.Z2_R:.2f} Ohm, X={relay.Z2_X:.2f} Ohm (300ms)")
    md_content.append(f"- Strefa 3: R={relay.Z3_R:.2f} Ohm, X={relay.Z3_X:.2f} Ohm (600ms)")
    md_content.append("")

    # Sprawdz czy bylo zwarcie
    trip_any = any(trip_history_L1) or any(trip_history_L2) or any(trip_history_L3)

    if trip_any:
        md_content.append("## Wykrycie zwarcia")
        md_content.append("")
        md_content.append("**Zwarcie wykryte - zadzialanie zabezpieczenia**")
        md_content.append("")

        # Znajdz czas pierwszego zadzialania
        fault_time = None
        fault_phase = None

        for i, (t1, t2, t3) in enumerate(zip(trip_history_L1, trip_history_L2, trip_history_L3)):
            if t1 or t2 or t3:
                fault_time = t[N + i]
                if t1:
                    fault_phase = "L1-E"
                elif t2:
                    fault_phase = "L2-E"
                elif t3:
                    fault_phase = "L3-E"
                break

        if fault_time is not None:
            md_content.append(f"- Czas wykrycia: {fault_time:.4f} s")
            md_content.append(f"- Faza: {fault_phase}")
            md_content.append("")

            # Wartosci impedancji w momencie zwarcia
            idx = int((fault_time - t[N]) * Fs) if Fs > 0 else 0
            idx = max(0, min(idx, len(z1_r_history) - 1))

            if fault_phase == "L1-E":
                z_r, z_x = z1_r_history[idx], z1_x_history[idx]
            elif fault_phase == "L2-E":
                z_r, z_x = z2_r_history[idx], z2_x_history[idx]
            else:
                z_r, z_x = z3_r_history[idx], z3_x_history[idx]

            z_mag = math.sqrt(z_r**2 + z_x**2) if z_r or z_x else 0
            z_angle = math.degrees(math.atan2(z_x, z_r)) if z_r else 0

            md_content.append("## Wartosci w momencie wykrycia zwarcia")
            md_content.append(f"- R = {z_r:.4f} Ohm")
            md_content.append(f"- X = {z_x:.4f} Ohm")
            md_content.append(f"- |Z| = {z_mag:.4f} Ohm")
            md_content.append(f"- Kat Z = {z_angle:.2f} st.")
    else:
        md_content.append("## Wykrycie zwarcia")
        md_content.append("")
        md_content.append("**Brak wykrycia zwarcia**")
        md_content.append("")
        md_content.append("Zabezpieczenie nie zadzialalo w okresie rejestracji.")

    md_content.append("")
    md_content.append("---")
    md_content.append("*Wygenerowano automatycznie przez tester zabezpieczenia odleglosciowego*")

    return "\n".join(md_content)

# Generuj i zapisz rezultat.md
result_md = generate_result_md(
    t, trip_history_L1, trip_history_L2, trip_history_L3,
    z1_r_history, z1_x_history, z2_r_history, z2_x_history, z3_r_history, z3_x_history,
    Z_line_R, Z_line_X, relay
)

with open('rezultat.md', 'w', encoding='utf-8') as f:
    f.write(result_md)

print("Zapisano rezultat.md")