import subprocess
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

# 1. NEC-Datei für reale Aufbauhöhe aus dem Video
# Unterkante/Speisepunkt bei Z = 2m, Spitze bei Z = 6.5m
nec_content = """CM 10m Delta Loop 3D Video Setup (2m base)
CE
GW 1 11 0 0 2 -1.34 0 6.5 0.0007
GW 2 9 -1.34 0 6.5 1.34 0 6.5 0.0007
GW 3 11 1.34 0 6.5 0 0 2 0.0007
GE 1
GN 2 0 0 0 13 0.005
EK
EX 0 1 1 0 1 0
FR 0 1 0 0 28.5 0
RP 0 37 73 1000 0 0 5 5
EN
"""

with open("pattern_3d_2m.nec", "w") as f:
    f.write(nec_content)

# 2. Simulation ausführen
subprocess.run(["nec2c", "-i", "pattern_3d_2m.nec", "-o", "pattern_3d_2m.out"], stdout=subprocess.DEVNULL)

# 3. Daten parsen
theta_vals = []
phi_vals = []
gain_dict = {}

with open("pattern_3d_2m.out", "r") as f:
    lines = f.readlines()
    in_pattern = False
    for line in lines:
        if "RADIATION PATTERNS" in line:
            in_pattern = True
            continue
        
        if in_pattern:
            parts = line.strip().split()
            if len(parts) >= 8:
                try:
                    th = float(parts[0])
                    ph = float(parts[1])
                    g = float(parts[4])
                    
                    if th not in theta_vals: theta_vals.append(th)
                    if ph not in phi_vals: phi_vals.append(ph)
                    gain_dict[(th, ph)] = g
                except ValueError:
                    continue

theta_vals = sorted(theta_vals)
phi_vals = sorted(phi_vals)

THETA, PHI = np.meshgrid(np.radians(theta_vals), np.radians(phi_vals))
R_raw = np.zeros(THETA.shape)

for i, ph in enumerate(phi_vals):
    for j, th in enumerate(theta_vals):
        R_raw[i, j] = gain_dict.get((th, ph), -30.0)

# 4. Leistungsskalierung (10^(dBi/10)) für realistisches Formverhältnis
R_lin = 10 ** (np.maximum(R_raw, -20.0) / 10.0)

# Umrechnung in Kartesische Koordinaten (X, Y, Z)
X = R_lin * np.sin(THETA) * np.cos(PHI)
Y = R_lin * np.sin(THETA) * np.sin(PHI)
Z = R_lin * np.cos(THETA)

# 5. 3D-Plot erstellen
fig = plt.figure(figsize=(10, 8))
ax = fig.add_subplot(111, projection='3d')

# Farbskala basierend auf den echten dBi-Werten
norm = plt.Normalize(R_raw.min(), R_raw.max())
colors = plt.cm.jet(norm(R_raw))

surf = ax.plot_surface(X, Y, Z, facecolors=colors, rstride=1, cstride=1, linewidth=0.2, antialiased=True, alpha=0.9)

ax.view_init(elev=20, azim=45)
ax.set_title("10m Delta-Loop (Speisepunkt: 2m) – 3D Abstrahlmuster", fontsize=14, fontweight='bold', pad=20)
ax.set_axis_off()

plt.savefig("delta_loop_3d_2m.png", dpi=300, bbox_inches='tight')
print("Diagramm für 2m Aufbauhöhe gespeichert als 'delta_loop_3d_2m.png'!")