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

# 1. Höhere Delta-Loop für flachere Strahlung (DX)
# Unterkante bei Z = 8m, Spitze bei Z = 12m (vorher 5m / 9m)
nec_content = """CM 10m Delta Loop 3D DX-Optimized (Higher)
CE
GW 1 11 0 0 8 -1.34 0 12 0.0007
GW 2 9 -1.34 0 12 1.34 0 12 0.0007
GW 3 11 1.34 0 12 0 0 8 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_dx.nec", "w") as f:
    f.write(nec_content)

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

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

with open("pattern_3d_dx.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. Optische Entzerrung: Lineare Leistungsskalierung (10^(dBi/10))
# Dadurch werden Unterschiede zwischen z.B. +8 dBi und +1 dBi optisch deutlich sichtbar!
R_lin = 10 ** (np.maximum(R_raw, -20.0) / 10.0)

# Umrechnung in Kartesische Koordinaten
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 (Höhe: 8-12m) – DX-3D-Abstrahlmuster", fontsize=14, fontweight='bold', pad=20)
ax.set_axis_off()

plt.savefig("delta_loop_3d_dx.png", dpi=300, bbox_inches='tight')
print("Neues DX-3D-Diagramm gespeichert als 'delta_loop_3d_dx.png'!")