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

# 1. NEC-Datei für volles 3D-Feld schreiben
# RP-Karte: 37 Theta-Schritte (0° bis 180° in 5°-Schritten), 73 Phi-Schritte (0° bis 360° in 5°-Schritten)
nec_content = """CM 10m Delta Loop 3D Pattern
CE
GW 1 11 0 0 5 -1.34 0 9 0.0007
GW 2 9 -1.34 0 9 1.34 0 9 0.0007
GW 3 11 1.34 0 9 0 0 5 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.nec", "w") as f:
    f.write(nec_content)

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

# 3. Daten für das 3D-Gitter parsen
theta_vals = []
phi_vals = []
gain_dict = {}

with open("pattern_3d.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)

# Matrix aufbauen
THETA, PHI = np.meshgrid(np.radians(theta_vals), np.radians(phi_vals))
R = np.zeros(THETA.shape)

for i, ph in enumerate(phi_vals):
    for j, th in enumerate(theta_vals):
        val = gain_dict.get((th, ph), -30.0)
        # Offset für Darstellung (Min-Wert auf 0 verschieben, damit der Radius positiv bleibt)
        R[i, j] = max(val + 30.0, 0.0)

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

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

# Farbskala basierend auf dem Radius (Gewinn)
norm = plt.Normalize(R.min(), R.max())
colors = plt.cm.jet(norm(R))

# 3D-Surface zeichnen
surf = ax.plot_surface(X, Y, Z, facecolors=colors, rstride=1, cstride=1, linewidth=0.2, antialiased=True, alpha=0.9)

# Blickwinkel einstellen (Elevation, Azimut)
ax.view_init(elev=25, azim=45)

# Achsenbeschriftungen & Titel
ax.set_title("10m Delta-Loop – 3D Strahlungsdiagramm", fontsize=14, fontweight='bold', pad=20)
ax.set_axis_off()  # Achsen ausblenden für eine cleane Ansicht der Antennenkeule

# Speichern
plt.savefig("delta_loop_3d.png", dpi=300, bbox_inches='tight')
print("3D-Diagramm erfolgreich gespeichert als 'delta_loop_3d.png'!")