你已经派生过 gprMax
镜像自地址
https://gitee.com/sunhf/gprMax.git
已同步 2025-08-07 04:56:51 +08:00
Updated for input impedance and admittance calculations.
这个提交包含在：
Craig Warren
@@ -19,16 +19,15 @@
|
||||
import os, argparse
|
||||
import h5py
|
||||
import numpy as np
|
||||
np.seterr(divide='ignore', invalid='ignore')
|
||||
import matplotlib.pyplot as plt
|
||||
import matplotlib.gridspec as gridspec
|
||||
|
||||
from gprMax.exceptions import CmdInputError
|
||||
from gprMax.constants import complextype
|
||||
|
||||
"""Plots antenna parameters (s11 parameter and input impedance) from an output file containing a transmission line source."""
|
||||
"""Plots antenna parameters (s11 parameter and input impedance and admittance) from an output file containing a transmission line source."""
|
||||
|
||||
# Parse command line arguments
|
||||
parser = argparse.ArgumentParser(description='Plots antenna parameters (s11 parameter and input impedance) from an output file containing a transmission line source.', usage='cd gprMax; python -m tools.plot_antenna_params outputfile')
|
||||
parser = argparse.ArgumentParser(description='Plots antenna parameters (s11 parameter and input impedance and admittance) from an output file containing a transmission line source.', usage='cd gprMax; python -m tools.plot_antenna_params outputfile')
|
||||
parser.add_argument('outputfile', help='name of output file including path')
|
||||
parser.add_argument('-tln', default=1, type=int, help='transmission line number')
|
||||
args = parser.parse_args()
|
||||
@@ -41,116 +40,250 @@ iterations = f.attrs['Iterations']
|
||||
time = np.arange(0, dt * iterations, dt)
|
||||
time = time / 1e-9
|
||||
|
||||
# Read/calculate voltages and currents
|
||||
path = '/tls/tl' + str(args.tln) + '/'
|
||||
Vinc = f[path + 'Vinc'][:]
|
||||
Vscat = f[path + 'Vscat'][:]
|
||||
Iscat = f[path + 'Iscat'][:]
|
||||
Iinc = f[path + 'Iinc'][:]
|
||||
Vtotal = f[path +'Vtotal'][:]
|
||||
Itotal = f[path +'Itotal'][:]
|
||||
Vref = Vtotal - Vinc
|
||||
Iref = Itotal - Iinc
|
||||
|
||||
# Calculate magnitude of frequency spectra
|
||||
Vincp = np.abs(np.fft.fft(Vinc))**2
|
||||
# Calculate magnitude of FFTs of voltages and currents
|
||||
Vincp = np.abs(np.fft.fft(Vinc))
|
||||
freqs = np.fft.fftfreq(Vincp.size, d=dt)
|
||||
Vscatp = np.abs(np.fft.fft(Vscat))**2
|
||||
s11 = Vscatp / Vincp
|
||||
zin = np.zeros(iterations, dtype=np.complex)
|
||||
zin = np.abs(np.fft.fft(Vscat)) / np.abs(np.fft.fft(Iscat))
|
||||
delaycorrection = np.zeros(Vincp.size, dtype=complextype)
|
||||
delaycorrection = np.exp(-1j * np.pi * freqs * dt)
|
||||
Iincp = np.abs(np.fft.fft(Iinc))
|
||||
Vrefp = np.abs(np.fft.fft(Vref))
|
||||
Irefp = np.abs(np.fft.fft(Iref))
|
||||
Vtotalp = np.abs(np.fft.fft(Vtotal))
|
||||
Itotalp = np.abs(np.fft.fft(Itotal))
|
||||
|
||||
# Calculate s11
|
||||
s11 = Vrefp / Vincp
|
||||
|
||||
# Calculate input impedance
|
||||
zin = np.zeros(iterations, dtype=complextype)
|
||||
zin = (np.fft.fft(Vtotal) * delaycorrection) / np.fft.fft(Itotal)
|
||||
|
||||
# Calculate input admittance
|
||||
yin = np.zeros(iterations, dtype=complextype)
|
||||
yin = np.fft.fft(Itotal) / (np.fft.fft(Vtotal) * delaycorrection)
|
||||
|
||||
# Convert to decibels
|
||||
Vincp = 10 * np.log10(Vincp)
|
||||
Vscatp = 10 * np.log10(Vscatp)
|
||||
s11 = 10 * np.log10(s11)
|
||||
Vincp = 20 * np.log10(Vincp)
|
||||
Iincp = 20 * np.log10(Iincp)
|
||||
Vrefp = 20 * np.log10(Vrefp)
|
||||
Irefp = 20 * np.log10(Irefp)
|
||||
Vtotalp = 20 * np.log10(Vtotalp)
|
||||
Itotalp = 20 * np.log10(Itotalp)
|
||||
s11 = 20 * np.log10(s11)
|
||||
|
||||
# Set plotting range to -60dB from maximum power
|
||||
pltrange = np.where((np.amax(Vincp) - Vincp) > 60)[0][0] + 1
|
||||
pltrange = np.s_[0:pltrange]
|
||||
# Set plotting range
|
||||
# To a certain drop from maximum power
|
||||
#pltrange = np.where((np.amax(Vincp) - Vincp) > 30)[0][0] + 1
|
||||
# To a maximum frequency
|
||||
pltrange = np.where(freqs > 2e9)[0][0]
|
||||
pltrange = np.s_[1:pltrange]
|
||||
|
||||
# Figure 1
|
||||
# Plot incident voltage
|
||||
fig1, ax = plt.subplots(num='Transmission line parameters', figsize=(20, 10), facecolor='w', edgecolor='w')
|
||||
gs1 = gridspec.GridSpec(2, 2, hspace=0.3)
|
||||
ax1 = plt.subplot(gs1[0, 0])
|
||||
ax1.plot(time, Vinc, 'r', lw=2, label='Vinc')
|
||||
ax1.set_title('Incident voltage')
|
||||
ax1.set_xlabel('Time [ns]')
|
||||
ax1.set_ylabel('Voltage [V]')
|
||||
ax1.set_xlim([0, np.amax(time)])
|
||||
ax1.grid()
|
||||
fig1, ax = plt.subplots(num='Transmission line parameters', figsize=(20, 12), facecolor='w', edgecolor='w')
|
||||
gs1 = gridspec.GridSpec(4, 2, hspace=0.7)
|
||||
ax = plt.subplot(gs1[0, 0])
|
||||
ax.plot(time, Vinc, 'r', lw=2, label='Vinc')
|
||||
ax.set_title('Incident voltage')
|
||||
ax.set_xlabel('Time [ns]')
|
||||
ax.set_ylabel('Voltage [V]')
|
||||
ax.set_xlim([0, np.amax(time)])
|
||||
ax.grid()
|
||||
|
||||
# Plot frequency spectra of incident voltage
|
||||
ax2 = plt.subplot(gs1[0, 1])
|
||||
markerline, stemlines, baseline = ax2.stem(freqs[pltrange]/1e9, Vincp[pltrange], '-.')
|
||||
ax = plt.subplot(gs1[0, 1])
|
||||
markerline, stemlines, baseline = ax.stem(freqs[pltrange]/1e9, Vincp[pltrange], '-.')
|
||||
plt.setp(baseline, 'linewidth', 0)
|
||||
plt.setp(stemlines, 'color', 'r')
|
||||
plt.setp(markerline, 'markerfacecolor', 'r', 'markeredgecolor', 'r')
|
||||
ax2.plot(freqs[pltrange]/1e9, Vincp[pltrange], 'r', lw=2)
|
||||
ax2.set_title('Incident voltage')
|
||||
ax2.set_xlabel('Frequency [GHz]')
|
||||
ax2.set_ylabel('Power [dB]')
|
||||
ax2.grid()
|
||||
ax.plot(freqs[pltrange]/1e9, Vincp[pltrange], 'r', lw=2)
|
||||
ax.set_title('Incident voltage')
|
||||
ax.set_xlabel('Frequency [GHz]')
|
||||
ax.set_ylabel('Power [dB]')
|
||||
ax.grid()
|
||||
|
||||
# Plot scattered (field) voltage
|
||||
ax3 = plt.subplot(gs1[1, 0])
|
||||
ax3.plot(time, Vscat, 'r', lw=2, label='Vscat')
|
||||
ax3.set_title('Reflected voltage')
|
||||
ax3.set_xlabel('Time [ns]')
|
||||
ax3.set_ylabel('Voltage [V]')
|
||||
ax3.set_xlim([0, np.amax(time)])
|
||||
ax3.grid()
|
||||
# Plot incident current
|
||||
ax = plt.subplot(gs1[1, 0])
|
||||
ax.plot(time, Iinc, 'b', lw=2, label='Vinc')
|
||||
ax.set_title('Incident current')
|
||||
ax.set_xlabel('Time [ns]')
|
||||
ax.set_ylabel('Current [A]')
|
||||
ax.set_xlim([0, np.amax(time)])
|
||||
ax.grid()
|
||||
|
||||
# Plot frequency spectra of scattered voltage
|
||||
ax4 = plt.subplot(gs1[1, 1])
|
||||
markerline, stemlines, baseline = ax4.stem(freqs[pltrange]/1e9, Vscatp[pltrange], '-.')
|
||||
# Plot frequency spectra of incident current
|
||||
ax = plt.subplot(gs1[1, 1])
|
||||
markerline, stemlines, baseline = ax.stem(freqs[pltrange]/1e9, Iincp[pltrange], '-.')
|
||||
plt.setp(baseline, 'linewidth', 0)
|
||||
plt.setp(stemlines, 'color', 'b')
|
||||
plt.setp(markerline, 'markerfacecolor', 'b', 'markeredgecolor', 'b')
|
||||
ax.plot(freqs[pltrange]/1e9, Iincp[pltrange], 'b', lw=2)
|
||||
ax.set_title('Incident current')
|
||||
ax.set_xlabel('Frequency [GHz]')
|
||||
ax.set_ylabel('Power [dB]')
|
||||
ax.grid()
|
||||
|
||||
# Plot total voltage
|
||||
ax = plt.subplot(gs1[2, 0])
|
||||
ax.plot(time, Vtotal, 'r', lw=2, label='Vinc')
|
||||
ax.set_title('Total (incident + reflected) voltage')
|
||||
ax.set_xlabel('Time [ns]')
|
||||
ax.set_ylabel('Voltage [V]')
|
||||
ax.set_xlim([0, np.amax(time)])
|
||||
ax.grid()
|
||||
|
||||
# Plot frequency spectra of total voltage
|
||||
ax = plt.subplot(gs1[2, 1])
|
||||
markerline, stemlines, baseline = ax.stem(freqs[pltrange]/1e9, Vtotalp[pltrange], '-.')
|
||||
plt.setp(baseline, 'linewidth', 0)
|
||||
plt.setp(stemlines, 'color', 'r')
|
||||
plt.setp(markerline, 'markerfacecolor', 'r', 'markeredgecolor', 'r')
|
||||
ax4.plot(freqs[pltrange]/1e9, Vscatp[pltrange], 'r', lw=2)
|
||||
ax4.set_title('Reflected voltage')
|
||||
ax4.set_xlabel('Frequency [GHz]')
|
||||
ax4.set_ylabel('Power [dB]')
|
||||
ax4.grid()
|
||||
ax.plot(freqs[pltrange]/1e9, Vtotalp[pltrange], 'r', lw=2)
|
||||
ax.set_title('Total (incident + reflected) voltage')
|
||||
ax.set_xlabel('Frequency [GHz]')
|
||||
ax.set_ylabel('Power [dB]')
|
||||
ax.grid()
|
||||
|
||||
# Plot total current
|
||||
ax = plt.subplot(gs1[3, 0])
|
||||
ax.plot(time, Itotal, 'b', lw=2, label='Vinc')
|
||||
ax.set_title('Total (incident + reflected) current')
|
||||
ax.set_xlabel('Time [ns]')
|
||||
ax.set_ylabel('Current [A]')
|
||||
ax.set_xlim([0, np.amax(time)])
|
||||
ax.grid()
|
||||
|
||||
# Plot frequency spectra of reflected current
|
||||
ax = plt.subplot(gs1[3, 1])
|
||||
markerline, stemlines, baseline = ax.stem(freqs[pltrange]/1e9, Itotalp[pltrange], '-.')
|
||||
plt.setp(baseline, 'linewidth', 0)
|
||||
plt.setp(stemlines, 'color', 'b')
|
||||
plt.setp(markerline, 'markerfacecolor', 'b', 'markeredgecolor', 'b')
|
||||
ax.plot(freqs[pltrange]/1e9, Itotalp[pltrange], 'b', lw=2)
|
||||
ax.set_title('Total (incident + reflected) current')
|
||||
ax.set_xlabel('Frequency [GHz]')
|
||||
ax.set_ylabel('Power [dB]')
|
||||
ax.grid()
|
||||
|
||||
## Plot reflected (reflected) voltage
|
||||
#ax = plt.subplot(gs1[4, 0])
|
||||
#ax.plot(time, Vref, 'r', lw=2, label='Vref')
|
||||
#ax.set_title('Reflected voltage')
|
||||
#ax.set_xlabel('Time [ns]')
|
||||
#ax.set_ylabel('Voltage [V]')
|
||||
#ax.set_xlim([0, np.amax(time)])
|
||||
#ax.grid()
|
||||
#
|
||||
## Plot frequency spectra of reflected voltage
|
||||
#ax = plt.subplot(gs1[4, 1])
|
||||
#markerline, stemlines, baseline = ax.stem(freqs[pltrange]/1e9, Vrefp[pltrange], '-.')
|
||||
#plt.setp(baseline, 'linewidth', 0)
|
||||
#plt.setp(stemlines, 'color', 'r')
|
||||
#plt.setp(markerline, 'markerfacecolor', 'r', 'markeredgecolor', 'r')
|
||||
#ax.plot(freqs[pltrange]/1e9, Vrefp[pltrange], 'r', lw=2)
|
||||
#ax.set_title('Reflected voltage')
|
||||
#ax.set_xlabel('Frequency [GHz]')
|
||||
#ax.set_ylabel('Power [dB]')
|
||||
#ax.grid()
|
||||
#
|
||||
## Plot reflected (reflected) current
|
||||
#ax = plt.subplot(gs1[5, 0])
|
||||
#ax.plot(time, Iref, 'b', lw=2, label='Iref')
|
||||
#ax.set_title('Reflected current')
|
||||
#ax.set_xlabel('Time [ns]')
|
||||
#ax.set_ylabel('Current [A]')
|
||||
#ax.set_xlim([0, np.amax(time)])
|
||||
#ax.grid()
|
||||
#
|
||||
## Plot frequency spectra of reflected current
|
||||
#ax = plt.subplot(gs1[5, 1])
|
||||
#markerline, stemlines, baseline = ax.stem(freqs[pltrange]/1e9, Irefp[pltrange], '-.')
|
||||
#plt.setp(baseline, 'linewidth', 0)
|
||||
#plt.setp(stemlines, 'color', 'b')
|
||||
#plt.setp(markerline, 'markerfacecolor', 'b', 'markeredgecolor', 'b')
|
||||
#ax.plot(freqs[pltrange]/1e9, Irefp[pltrange], 'b', lw=2)
|
||||
#ax.set_title('Reflected current')
|
||||
#ax.set_xlabel('Frequency [GHz]')
|
||||
#ax.set_ylabel('Power [dB]')
|
||||
#ax.grid()
|
||||
|
||||
# Figure 2
|
||||
# Plot frequency spectra of s11
|
||||
fig2, ax = plt.subplots(num='Antenna parameters', figsize=(20, 10), facecolor='w', edgecolor='w')
|
||||
gs2 = gridspec.GridSpec(2, 2, hspace=0.3)
|
||||
ax5 = plt.subplot(gs2[0, 0])
|
||||
markerline, stemlines, baseline = ax5.stem(freqs[pltrange]/1e9, s11[pltrange], '-.')
|
||||
fig2, ax = plt.subplots(num='Antenna parameters', figsize=(20, 12), facecolor='w', edgecolor='w')
|
||||
gs2 = gridspec.GridSpec(3, 2, hspace=0.5)
|
||||
ax = plt.subplot(gs2[0, 0])
|
||||
markerline, stemlines, baseline = ax.stem(freqs[pltrange]/1e9, s11[pltrange], '-.')
|
||||
plt.setp(baseline, 'linewidth', 0)
|
||||
plt.setp(stemlines, 'color', 'r')
|
||||
plt.setp(markerline, 'markerfacecolor', 'r', 'markeredgecolor', 'r')
|
||||
ax5.plot(freqs[pltrange]/1e9, s11[pltrange], 'r', lw=2)
|
||||
ax5.set_title('s11 parameter')
|
||||
ax5.set_xlabel('Frequency [GHz]')
|
||||
ax5.set_ylabel('Power [dB]')
|
||||
ax5.grid()
|
||||
plt.setp(stemlines, 'color', 'g')
|
||||
plt.setp(markerline, 'markerfacecolor', 'g', 'markeredgecolor', 'g')
|
||||
ax.plot(freqs[pltrange]/1e9, s11[pltrange], 'g', lw=2)
|
||||
ax.set_title('s11')
|
||||
ax.set_xlabel('Frequency [GHz]')
|
||||
ax.set_ylabel('Power [dB]')
|
||||
ax.grid()
|
||||
|
||||
# Plot input resistance (real part of impedance)
|
||||
ax6 = plt.subplot(gs2[1, 0])
|
||||
markerline, stemlines, baseline = ax6.stem(freqs[pltrange]/1e9, zin[pltrange].real, '-.')
|
||||
ax = plt.subplot(gs2[1, 0])
|
||||
markerline, stemlines, baseline = ax.stem(freqs[pltrange]/1e9, np.abs(zin[pltrange]), '-.')
|
||||
plt.setp(baseline, 'linewidth', 0)
|
||||
plt.setp(stemlines, 'color', 'r')
|
||||
plt.setp(markerline, 'markerfacecolor', 'r', 'markeredgecolor', 'r')
|
||||
ax6.plot(freqs[pltrange]/1e9, zin[pltrange].real, 'r', lw=2)
|
||||
ax6.set_title('Input impedance')
|
||||
ax6.set_xlabel('Frequency [GHz]')
|
||||
ax6.set_ylabel('Resistance [Ohms]')
|
||||
ax6.set_ylim(bottom=0)
|
||||
ax6.grid()
|
||||
plt.setp(stemlines, 'color', 'g')
|
||||
plt.setp(markerline, 'markerfacecolor', 'g', 'markeredgecolor', 'g')
|
||||
ax.plot(freqs[pltrange]/1e9, np.abs(zin[pltrange]), 'g', lw=2)
|
||||
ax.set_title('Input impedance (resistive)')
|
||||
ax.set_xlabel('Frequency [GHz]')
|
||||
ax.set_ylabel('Resistance [Ohms]')
|
||||
ax.set_ylim([0, 2000])
|
||||
ax.grid()
|
||||
|
||||
# Plot input reactance (imaginery part of impedance)
|
||||
ax7 = plt.subplot(gs2[1, 1])
|
||||
markerline, stemlines, baseline = ax7.stem(freqs[pltrange]/1e9, zin[pltrange].imag, '-.')
|
||||
ax = plt.subplot(gs2[1, 1])
|
||||
markerline, stemlines, baseline = ax.stem(freqs[pltrange]/1e9, zin[pltrange].imag, '-.')
|
||||
plt.setp(baseline, 'linewidth', 0)
|
||||
plt.setp(stemlines, 'color', 'r')
|
||||
plt.setp(markerline, 'markerfacecolor', 'r', 'markeredgecolor', 'r')
|
||||
ax7.plot(freqs[pltrange]/1e9, zin[pltrange].imag, 'r', lw=2)
|
||||
ax7.set_title('Input impedance')
|
||||
ax7.set_xlabel('Frequency [GHz]')
|
||||
ax7.set_ylabel('Reactance [Ohms]')
|
||||
ax7.set_ylim(bottom=0)
|
||||
ax7.grid()
|
||||
plt.setp(stemlines, 'color', 'g')
|
||||
plt.setp(markerline, 'markerfacecolor', 'g', 'markeredgecolor', 'g')
|
||||
ax.plot(freqs[pltrange]/1e9, zin[pltrange].imag, 'g', lw=2)
|
||||
ax.set_title('Input impedance (reactive)')
|
||||
ax.set_xlabel('Frequency [GHz]')
|
||||
ax.set_ylabel('Reactance [Ohms]')
|
||||
ax.set_ylim(-2000, 2000)
|
||||
ax.grid()
|
||||
|
||||
# Plot input admittance (magnitude)
|
||||
ax = plt.subplot(gs2[2, 0])
|
||||
markerline, stemlines, baseline = ax.stem(freqs[pltrange]/1e9, np.abs(yin[pltrange]), '-.')
|
||||
plt.setp(baseline, 'linewidth', 0)
|
||||
plt.setp(stemlines, 'color', 'g')
|
||||
plt.setp(markerline, 'markerfacecolor', 'g', 'markeredgecolor', 'g')
|
||||
ax.plot(freqs[pltrange]/1e9, np.abs(yin[pltrange]), 'g', lw=2)
|
||||
ax.set_title('Input admittance (magnitude)')
|
||||
ax.set_xlabel('Frequency [GHz]')
|
||||
ax.set_ylabel('Admittance [Siemens]')
|
||||
ax.grid()
|
||||
|
||||
# Plot input admittance (phase)
|
||||
ax = plt.subplot(gs2[2, 1])
|
||||
markerline, stemlines, baseline = ax.stem(freqs[pltrange]/1e9, np.angle(yin[pltrange], deg=True), '-.')
|
||||
plt.setp(baseline, 'linewidth', 0)
|
||||
plt.setp(stemlines, 'color', 'g')
|
||||
plt.setp(markerline, 'markerfacecolor', 'g', 'markeredgecolor', 'g')
|
||||
ax.plot(freqs[pltrange]/1e9, np.angle(yin[pltrange], deg=True), 'g', lw=2)
|
||||
ax.set_title('Input admittance (phase)')
|
||||
ax.set_xlabel('Frequency [GHz]')
|
||||
ax.set_ylabel('Phase [degrees]')
|
||||
ax.grid()
|
||||
|
||||
plt.show()
|
||||
fig1.savefig(os.path.splitext(os.path.abspath(file))[0] + '_tl_params.png', dpi=150, format='png', bbox_inches='tight', pad_inches=0.1)
|
||||
fig2.savefig(os.path.splitext(os.path.abspath(file))[0] + '_ant_params.png', dpi=150, format='png', bbox_inches='tight', pad_inches=0.1)
|
||||
#fig1.savefig(os.path.splitext(os.path.abspath(file))[0] + '_tl_params.png', dpi=150, format='png', bbox_inches='tight', pad_inches=0.1)
|
||||
#fig2.savefig(os.path.splitext(os.path.abspath(file))[0] + '_ant_params.png', dpi=150, format='png', bbox_inches='tight', pad_inches=0.1)
|
||||
#fig1.savefig(os.path.splitext(os.path.abspath(file))[0] + '_tl_params.pdf', dpi=None, format='pdf', bbox_inches='tight', pad_inches=0.1)
|
||||
#fig2.savefig(os.path.splitext(os.path.abspath(file))[0] + '_ant_params.pdf', dpi=None, format='pdf', bbox_inches='tight', pad_inches=0.1)
|
||||
f.close()
|
||||
在新工单中引用
屏蔽一个用户