You can not select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
513 lines
14 KiB
513 lines
14 KiB
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")
|
|
|
|
# ----------------------------------------------- #
|
|
|
|
|
|
|
|
|
|
|