import os
import numpy as np
import pandas as pd
import shutil
from concurrent.futures import ProcessPoolExecutor
from tqdm import tqdm
import sys
import matplotlib.pyplot as plt
import csv
import fnmatch
import random
import copy
import concurrent.futures
import urllib.request
from concurrent import futures
import glob
import re

def extract_selected_log_values(filepath):
    targets = {
        "Final Time(sec)": None,
        "dt_sample": None,
        "tsPerSampling": None
    }

    with open(filepath, 'r') as f:
        lines = f.readlines()

    in_block = False
    for line in lines:
        if "==========================================" in line and "PERFORMING INFORMATION" in line:
            in_block = True
            continue
        if in_block and "==========================================" in line:
            break  # End of block

        for key in targets:
            if key in line:
                match = re.search(r'=\s*([\d.eE+-]+)', line)
                if match:
                    targets[key] = float(match.group(1))

    return {
        "FinalTime": targets["Final Time(sec)"],
        "dt_sample": targets["dt_sample"],
        "tsPerSampling": int(targets["tsPerSampling"])
    }
    

# Read Log File
log_file_path = "computeDGTD.log"  # Replace with your log file path
result = extract_selected_log_values(log_file_path)

FinalTime = result["FinalTime"]
dt_sample = result["dt_sample"]
tsPerSampling = result["tsPerSampling"]

print("FinalTime       =", FinalTime)
print("dt_sample       =", dt_sample)
print("tsPerSampling   =", tsPerSampling)

# === Constants and Paths ===
c0 = 3e8
current_dir = os.getcwd()

# ====================== USER DEFINED ========================== #
FileName = "PW"
sc_curFolder = current_dir+"/CURRENT_SC" 
inc_curFolder = current_dir+"/CURRENT_INC" 
sc_curFileName = sc_curFolder +"/Currents_"+FileName+"_" 
inc_curFileName = inc_curFolder +"/Currents_"+FileName+"_" 

DEST = current_dir+"/RCS"
# Time interval (Sampling interval dt)
dt = dt_sample
# Time step count
steps = tsPerSampling
# ====================== USER DEFINED ========================== #

# Sampling frequency f_s = 1 / dt
f_s = 1.0 / dt
# Nyquist frequency
f_nyq = f_s / 2.0

# Count .curJ files
nfiles = len(fnmatch.filter(os.listdir(sc_curFolder), '*.curJ'))
nfiles = int(nfiles)
print('Number of time domain current files:', nfiles)

# Number of FFT samples
NFFT = nfiles * 3







################################
####### Process .tri File ######
################################

def parse_mesh(filename):
    with open(filename, 'r') as f:
        lines = f.readlines()

    scale = float(lines[0].strip())
    num_nodes = int(lines[1].strip())

    # Read nodes
    nodes = np.array([list(map(float, line.strip().split())) for line in lines[2:2 + num_nodes]])

    num_triangles = int(lines[3 + num_nodes].strip())
    
    # Read triangles (1-based indexing in file, converting to 0-based)
    triangles = np.array([list(map(int, line.strip().split())) for line in lines[4 + num_nodes:3 + num_nodes + num_triangles]]) 

    print(num_triangles)
    print(triangles[0])
    print(triangles[-1])

    return scale, nodes, triangles
   

tri_filename = FileName + "_out.tri"
scale, nodes, triangles = parse_mesh(tri_filename)

Num_Nodes = len(nodes)
Num_tri = len(triangles)
print('Number of Triangles =', Num_tri)
print('Number of Nodes =', Num_Nodes)



def count_triangles_from_timedomain_curfile(filename):
    with open(filename, 'r') as f:
        lines = f.readlines()

    # Remove empty lines and strip
    values = [line.strip() for line in lines if line.strip() != '']

    total_values = len(values)
    if total_values % 9 != 0:
        print(f"Warning: File has {total_values} values which is not divisible by 9!")

    num_triangles = total_values // 9
    return num_triangles

sample_file = sc_curFileName + "00000_BC.curJ"
num_tri = count_triangles_from_timedomain_curfile(sample_file)
print('Number of Triangles tri file =', Num_tri)
print('Number of Triangles time domain current file =', num_tri)





#####################################
####### Process Incident Field ######
#####################################

# Folder containing your output
folder = "CURRENT_INC"
prefix = "Einc_field_" + FileName + "_"
extension = ".dat"

# Get sorted list of files
file_list = sorted(glob.glob(os.path.join(folder, f"{prefix}*{extension}")))

# Prepare time and field vectors
time_steps = []
Ex_vals = []
Ey_vals = []
Ez_vals = []

for idx, filepath in enumerate(file_list):
    with open(filepath, 'r') as f:
        line = f.readline().strip()
        if line:
            Ex, Ey, Ez = map(float, line.split())
            Ex_vals.append(Ex)
            Ey_vals.append(Ey)
            Ez_vals.append(Ez)
            time_steps.append(idx * dt)  # Or parse from filename if non-uniform



# Convert to numpy arrays
time = np.array(time_steps)
Ex_vals = np.array(Ex_vals)
Ey_vals = np.array(Ey_vals)
Ez_vals = np.array(Ez_vals)

# Plot
plt.figure(figsize=(10, 6))
plt.plot(time, Ex_vals, label='Ex')
plt.plot(time, Ey_vals, label='Ey')
plt.plot(time, Ez_vals, label='Ez')
plt.xlabel("Time (s)")
plt.ylabel("E-field (V/m)")
plt.title("Incident Electric Field at Selected Node Over Time")
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.savefig("Incident_Field.png")


# Compute max absolute values
max_Ex = np.max(np.abs(Ex_vals))
max_Ey = np.max(np.abs(Ey_vals))
max_Ez = np.max(np.abs(Ez_vals))

# Find the dominant component
if max_Ex >= max_Ey and max_Ex >= max_Ez:
    signal = Ex_vals
    dominant = "Ex"
elif max_Ey >= max_Ex and max_Ey >= max_Ez:
    signal = Ey_vals
    dominant = "Ey"
else:
    signal = Ez_vals
    dominant = "Ez"

print(f"Dominant component: {dominant}")
print(f"Max magnitude: {np.max(np.abs(signal))}")


# Frequency Resolution
df = f_s / NFFT

# time axis
time_axis = np.arange(nfiles) * dt

# frequency axis
freq_axis = np.arange(NFFT) * df / 1e9

incident_fft = np.fft.fft(signal , n=NFFT) / NFFT

# --- Step 2: Find indices where magnitude > 10% of peak ---
fft_inc_magnitude = np.abs(incident_fft)
threshold = 0.1 * np.max(fft_inc_magnitude)
indices_above_threshold = np.where(fft_inc_magnitude[0:NFFT//2] > threshold)[0]


# Extract frequency range where magnitude > 10% of peak
freq_selected = freq_axis[indices_above_threshold]
fft_selected = fft_inc_magnitude[indices_above_threshold]

# Find min and max frequency for box range
fmin, fmax = np.min(freq_selected), np.max(freq_selected)

# Find indices of fmin and fmax
idx_min = np.searchsorted(freq_axis, fmin)
idx_max = np.searchsorted(freq_axis, fmax)

# Confirm the values at those indices
actual_fmin = freq_axis[idx_min]
actual_fmax = freq_axis[idx_max]

num_points = idx_max - idx_min

print(f"Index range: {idx_min} to {idx_max}")
print(f"Frequency range: {actual_fmin:.3f} GHz to {actual_fmax:.3f} GHz")
print(f"Number of frequency points in range: {num_points}")


# --- Step 3: Plot with translucent box ---

# Create subplots: time domain (top), frequency domain (bottom)
fig, axs = plt.subplots(2, 1, figsize=(10, 8))

# Plot time-domain signal
axs[0].plot(time_axis, signal)
axs[0].set_title("Time-Domain Signal (Gauss Pulse)")
axs[0].set_xlabel("Time (s)")
axs[0].set_ylabel("Amplitude")
axs[0].grid(True)

# Plot frequency-domain magnitude
axs[1].plot(freq_axis, fft_inc_magnitude)
axs[1].axvspan(fmin, fmax, color='red', alpha=0.3, label=">10% of Peak")
axs[1].set_title("Magnitude Spectrum of Incident Electric Field")
axs[1].set_xlabel("Frequency (Hz)")
axs[1].set_ylabel("Magnitude")
axs[1].grid(True)
axs[1].set_xlim([0, f_nyq / 1e9])  # Adjust based on your center freq
plt.legend()
plt.tight_layout()
plt.savefig("Incident.png")


# Create DataFrame
df = pd.DataFrame({
    'Frequency (GHz)': freq_axis,
    'fft_inc Magnitude': fft_inc_magnitude,
})

# Save to Excel
excel_path = "./cur_inc_freqdomain_analytical.xlsx"
df.to_excel(excel_path, index=False)






# Multi-threading
import concurrent.futures
import urllib.request



# ----------------------------------------------- #
# Scattered Field 
# ----------------------------------------------- #
# Read curJ and curM in time domain (Incident Field)
from concurrent import futures


# Read curJ and curM files using multi-threading
def read_cur_files(i, currentFilename, ext, Num_tri):
    index = "%05d" % int(steps * i)
    pfile = open(currentFilename + index + "_BC." + ext, 'r')
    lines = pfile.readlines()
    
    cur = np.zeros([Num_tri, 3, 3])
    
    counter = 0
    for t in range(Num_tri):
        for n in range(3):
            cur[t, n, 0] = float(lines[counter])
            counter += 1
            cur[t, n, 1] = float(lines[counter])
            counter += 1
            cur[t, n, 2] = float(lines[counter])
            counter += 1                
        counter += 1
    pfile.close()     
    return cur



# curJ = n x H
# Time, Triangle, node, xyz component
curJ_sc = np.zeros([nfiles, Num_tri, 3, 3])

# curM = E x n
# Time, Triangle, node, xyz component
curM_sc = np.zeros([nfiles, Num_tri, 3, 3])

print("Reading curJ sc ...")
with concurrent.futures.ThreadPoolExecutor(max_workers=16) as executor:
    tickets = {executor.submit(read_cur_files, i, sc_curFileName, "curJ", Num_tri): i for i in range(nfiles)}
    for future in tqdm(concurrent.futures.as_completed(tickets), total=nfiles, desc="curJ files"):
        index = tickets[future]
        try:
            curJ_sc[index, :, :, :] = future.result()
        except Exception as exc:
            print(f"{index} generated an exception: {exc}")

print("Reading curM sc...")
with concurrent.futures.ThreadPoolExecutor(max_workers=16) as executor:
    tickets = {executor.submit(read_cur_files, i, sc_curFileName, "curM", Num_tri): i for i in range(nfiles)}
    for future in tqdm(concurrent.futures.as_completed(tickets), total=nfiles, desc="curM files"):
        index = tickets[future]
        try:
            curM_sc[index, :, :, :] = future.result()
        except Exception as exc:
            print(f"{index} generated an exception: {exc}")








    

def process_RCS(F_curJ, F_curM, FREQ, FREQ_STR, outfile, DEST, SURFACE_TRI_MESH, THETA_sc, PHI_sc):

    # Create frequency folder and copy geometry
    freq_dir = f"{current_dir}/freq/FREQ{FREQ_STR}"
    os.makedirs(freq_dir, exist_ok=True)
    print("Working Folder =",freq_dir)
    
    # Change directory into frequency folder
    os.chdir(freq_dir)
    
    # Save results for current and magnetic field components
    fid = open(outfile + ".curJ", 'w')
    for t in tqdm(range(Num_tri)):
        for n in range(3):
            for c in range(3):
                line = str(np.real(F_curJ[t, n, c])) + " " + str(np.imag(F_curJ[t, n, c])) + '\n'
                fid.write(line)
    fid.close()

    fid = open(outfile + ".curM", 'w')
    for t in tqdm(range(Num_tri)):
        for n in range(3):
            for c in range(3):
                line = str(np.real(F_curM[t, n, c])) + " " + str(np.imag(F_curM[t, n, c])) + '\n'
                fid.write(line)
    fid.close()


    # Create .region file
    REGIONFILE = outfile
    with open(REGIONFILE + ".region", "w") as filep:
        filep.write("1\n")
        filep.write("0 FEBI\n")
        filep.write(f"0 FEBI {outfile}\n")
        filep.write("Coupling\n")

    # Copy mesh to expected filename
    shutil.copy(f"{current_dir}/{SURFACE_TRI_MESH}_out.tri", f"{freq_dir}/{outfile}.tri")
    shutil.copy(f"{current_dir}/n2f_main", f"{freq_dir}/n2f_main")

    # Run field-to-far-field converter and translator
    os.system(f"./n2f_main {REGIONFILE} {int(FREQ)} 0 0 0 180 360 0 1")
    os.system(f"emsurftranslator -s {outfile}")

    print("\n\n==== Process Farfield Result ====\n")

    # Clean spacing in .cs file
    os.system(f"sed -i 's/  / /g' {REGIONFILE}.cs")

    # Read and split farfield RCS data
    data = pd.read_csv(f"./{REGIONFILE}.cs", delimiter=" ", skiprows=1, header=None)
    data1 = data.iloc[:361*10]
    data1.drop(columns=[0], inplace=True)
    data2 = data.iloc[361*10:]
    data2.drop(columns=[data2.columns[-1]], inplace=True)

    # Concatenate full pattern data
    result = np.concatenate((data1.to_numpy(), data2.to_numpy()))
    df = pd.DataFrame(result, columns=['Theta', 'Phi', 'RCS_V', 'RCS_H'])

    # Extract RCS at specific angles
    theta_df = df[df['Theta'] == THETA_sc]
    theta_phi_df = theta_df[theta_df['Phi'] == PHI_sc]
    print(theta_phi_df)
    RCS_V = theta_phi_df['RCS_V'].values[0]  # single value
    RCS_H = theta_phi_df['RCS_H'].values[0]  # single value
    
    # Change back to base directory
    os.chdir("..")
    os.chdir("..")

    return RCS_V, RCS_H    
    
    
    
    
    
freq_all = []
RCS_V_all = []
RCS_H_all = []
print("----- Compute RCS -----")

for freq_ind in tqdm(range(num_points)):
    freq_index = idx_min + freq_ind   
    SURFACE_TRI_MESH = FileName
    outfile = "out"
    FREQ = int(freq_axis[freq_index] * 1000)
    FREQ_STR = str(FREQ)
    THETA_sc = 180
    PHI_sc = 0
    
    print("---------------- Processing ", FREQ_STR, " MHz ---------------------")
    
    norm = 1.0 / fft_inc_magnitude[freq_index]
    curJ_FD = np.fft.fft(curJ_sc, n=NFFT, axis=0) / NFFT * norm
    curM_FD = np.fft.fft(curM_sc, n=NFFT, axis=0) / NFFT * norm        

    F_curJ_single = curJ_FD[freq_index, :, :, :].copy()
    F_curM_single = curM_FD[freq_index, :, :, :].copy()

    RCS_V,RCS_H = process_RCS(F_curJ_single, F_curM_single, FREQ, FREQ_STR, outfile, DEST, SURFACE_TRI_MESH, THETA_sc, PHI_sc)
    
    freq_all.append(FREQ/1000)
    RCS_V_all.append(RCS_V)
    RCS_H_all.append(RCS_H)
    


# Plotting
plt.figure(figsize=(8, 5))
plt.plot(freq_all, RCS_V_all, marker='o',label='RCS_V')
plt.plot(freq_all, RCS_H_all, marker='x',label='RCS_H')
plt.legend()
plt.grid(True)
plt.xlabel("Frequency (GHz)")
plt.ylabel("RCS (dBsm)")
plt.title("RCS vs Frequency")
plt.tight_layout()
plt.savefig("RCS_vs_Frequency_new.png", dpi=300)

# Save to Excel
df_rcs = pd.DataFrame({
    "Frequency (GHz)": freq_all,
    "RCS_V (dBsm)": RCS_V_all,
    "RCS_H (dBsm)": RCS_H_all
})
df_rcs.to_excel("RCS_vs_Frequency_new.xlsx", index=False)
print("Saved to RCS_vs_Frequency_new.xlsx")

# ----------------------------------------------- #