I have this script to simulate a calibration of a phased array I would like to add to the plots the phase calibration values that should be used in the beamformer 
import numpy as np
import os
import matplotlib.pyplot as plt


# === NEC PARSING ===
def parse_currents_from_outfile(filename, center_seg=6, num_elements=91):
    currents_x = np.zeros(num_elements, dtype=complex)
    currents_y = np.zeros(num_elements, dtype=complex)
    sequence1 = [132 * n - 126 for n in range(1, num_elements + 1)]
    sequence2 = [132 * n - 60 for n in range(1, num_elements + 1)]
    with open(filename, "r") as f:
        lines = f.readlines()
    idx_x = idx_y = 0
    reading = False
    for line in lines:
        if "CURRENTS AND LOCATION" in line:
            reading = True
            continue
        if reading:
            if line.strip() == "" or "----" in line:
                continue
            parts = line.split()
            if len(parts) < 8:
                continue
            try:
                seq_no = int(parts[0])
                real = float(parts[6])
                imag = float(parts[7])
            except ValueError:
                continue
            if seq_no in sequence1:
                currents_x[idx_x] = complex(real, imag)
                idx_x += 1
            if seq_no in sequence2:
                currents_y[idx_y] = complex(real, imag)
                idx_y += 1
    return currents_x, currents_y


# === LOAD NEC-BASED COUPLING MATRIX ===
def load_nec_coupling_matrix(out_dir, num_elements, center_seg=6):
    N = num_elements
    C = np.zeros((2 * N, 2 * N), dtype=complex)
    for i in range(1, N + 1):
        fx = os.path.join(out_dir, f"excitation_X_elem{i}.out")
        fy = os.path.join(out_dir, f"excitation_Y_elem{i}.out")
        cx, cy = parse_currents_from_outfile(fx, center_seg, N)
        C[:N, i - 1] = cx
        C[N:, i - 1] = cy
        cx, cy = parse_currents_from_outfile(fy, center_seg, N)
        C[:N, N + i - 1] = cx
        C[N:, N + i - 1] = cy
    return C


# === RX CHAIN ERRORS ===
def generate_channel_errors(N, phase_std_deg=10, amp_std_db=1, cross_pol_db=-25):
    """
    Generate G matrix: diagonal RX phase + amplitude errors, with optional cross-pol leakage

    Parameters:
        N             : number of elements
        phase_std_deg: standard deviation of per-port phase error
        amp_std_db   : standard deviation of per-port amplitude error
        cross_pol_db : X↔Y leakage level in dB (e.g., -25 dB = 0.056)

    Returns:
        G : complex 2N x 2N gain matrix
    """
    total_ports = 2 * N
    G = np.zeros((total_ports, total_ports), dtype=complex)

    # Main diagonal (independent per-port gain + phase)
    phase_errors = np.random.normal(0, phase_std_deg * np.pi / 180, total_ports)
    amp_errors = np.random.normal(1, 10 ** (amp_std_db / 20) - 1, total_ports)
    for i in range(total_ports):
        G[i, i] = amp_errors[i] * np.exp(1j * phase_errors[i])

    # Cross-pol leakage
    if cross_pol_db is not None:
        cross_mag = 10 ** (cross_pol_db / 20)
        for i in range(N):
            # H -> V and V -> H leakage between same antenna ports
            phase_offset1 = np.random.uniform(0, 2 * np.pi)
            phase_offset2 = np.random.uniform(0, 2 * np.pi)
            G[i + N, i] += cross_mag * np.exp(1j * phase_offset1)  # X to Y
            G[i, i + N] += cross_mag * np.exp(1j * phase_offset2)  # Y to X

    return G


# === SIMULATE RX DATA ===
def simulate_rx(C, G, tx_vec, snr_db=40, num_samples=100000):
    rx_true = C @ G @ tx_vec
    noise_power = 10 ** (-snr_db / 10)
    rx_samples = np.zeros((num_samples, len(tx_vec)), dtype=complex)
    for i in range(num_samples):
        noise = np.sqrt(noise_power / 2) * (
            np.random.randn(len(tx_vec)) + 1j * np.random.randn(len(tx_vec))
        )
        rx_samples[i, :] = rx_true + noise
    return np.mean(rx_samples, axis=0)


# === PHASE CALIBRATION USING COUPLING MATRIX ===
def calibrate(rx_vec, C, lambda_reg=1e-6):
    C_H = C.conj().T
    C_inv = np.linalg.inv(C_H @ C + lambda_reg * np.eye(C.shape[0])) @ C_H
    gain_est = C_inv @ rx_vec
    ref_phase = np.angle(gain_est[0])  # Anchor to antenna 0 (X-pol)
    return gain_est * np.exp(-1j * ref_phase)


# === PHASE ERROR METRIC ===
def residual_phase(estimated, truth):
    phase_diff = np.angle(estimated / truth, deg=True)
    return (phase_diff + 180) % 360 - 180  # Wrap to [-180, 180]


# === MAIN PROCEDURE ===
def main():
    # === Parameters ===
    num_elements = 91
    nec_dir = "nec_files"  # Update as needed
    snr_db = 40
    num_samples = 100000
    phase_std_deg = 10
    amp_std_db = 1
    lambda_reg = 1e-6
    cross_pol_db = -20  # ← Try -20, -30, -100, etc. NOTE the model will break with worser values than -20

    print("Loading C matrix from NEC...")
    C = load_nec_coupling_matrix(nec_dir, num_elements)
    print(f"Condition number of C: {np.linalg.cond(C):.2e}")

    print("Generating RX channel errors (cables + receivers)...")
    G = generate_channel_errors(num_elements, phase_std_deg, amp_std_db, cross_pol_db)

    # === TX-X (H-polarized dipole) ===
    tx_x = np.zeros(2 * num_elements, dtype=complex)
    tx_x[:num_elements] = 1.0  # Only excite X-pol chain

    # === TX-Y (V-polarized dipole) ===
    tx_y = np.zeros(2 * num_elements, dtype=complex)
    tx_y[num_elements:] = 1.0  # Only excite Y-pol chain

    print("Simulating RX data for TX-X...")
    rx_x = simulate_rx(C, G, tx_x, snr_db, num_samples)

    print("Simulating RX data for TX-Y...")
    rx_y = simulate_rx(C, G, tx_y, snr_db, num_samples)

    print("Applying mutual coupling correction...")
    gain_est_x = calibrate(rx_x, C, lambda_reg)
    gain_est_y = calibrate(rx_y, C, lambda_reg)

    # === Ground truth from system gain matrix ===
    true_gain = np.diag(G)
    true_x = true_gain[:num_elements]
    true_y = true_gain[num_elements:]

    # === Get residuals ===
    res_x = residual_phase(gain_est_x[:num_elements], true_x)
    res_y = residual_phase(gain_est_y[num_elements:], true_y)

    mean_x, std_x = np.mean(res_x), np.std(res_x)
    mean_y, std_y = np.mean(res_y), np.std(res_y)

    print(f"\n🧪 Calibration Results (SNR={snr_db} dB, Samples={num_samples}):")
    print(f"X-pol: Mean {mean_x:+.2f}°, Std {std_x:.2f}°")
    print(f"Y-pol: Mean {mean_y:+.2f}°, Std {std_y:.2f}°")

    # === Visualization ===
    plt.figure(figsize=(14, 6))

    plt.subplot(1, 2, 1)
    plt.plot(res_x, "bo-", label="X-pol residuals")
    plt.axhline(0, ls="--", color="gray")
    plt.title("X-pol Phase Residuals")
    plt.xlabel("Antenna index")
    plt.ylabel("Phase Error (deg)")
    plt.grid(True)

    plt.subplot(1, 2, 2)
    plt.plot(res_y, "ro-", label="Y-pol residuals")
    plt.axhline(0, ls="--", color="gray")
    plt.title("Y-pol Phase Residuals")
    plt.xlabel("Antenna index")
    plt.grid(True)

    plt.tight_layout()
    plt.show()


if __name__ == "__main__":
    main()