From 72071473375c0a2123c94f63387ec760ffa2e1d4 Mon Sep 17 00:00:00 2001 From: Craig Warren Date: Tue, 24 Nov 2020 11:58:31 +0000 Subject: [PATCH] Tidied code formatting and improved comments. --- gprMax/sources.py | 155 +++++++++++++--------------------------------- 1 file changed, 43 insertions(+), 112 deletions(-) diff --git a/gprMax/sources.py b/gprMax/sources.py index 512f5ae7..fbbf32c9 100644 --- a/gprMax/sources.py +++ b/gprMax/sources.py @@ -97,21 +97,27 @@ class VoltageSource(Source): if self.polarisation == 'x': if self.resistance != 0: - Ex[i, j, k] -= updatecoeffsE[ID[G.IDlookup[componentID], i, j, k], 4] * self.waveformvaluesJ[iteration] * (1 / (self.resistance * G.dy * G.dz)) + Ex[i, j, k] -= (updatecoeffsE[ID[G.IDlookup[componentID], i, j, k], 4] + * self.waveformvaluesJ[iteration] + * (1 / (self.resistance * G.dy * G.dz))) else: - Ex[i, j, k] = -1 * self.waveformvaluesJ[iteration] / G.dx + Ex[i, j, k] = - self.waveformvaluesJ[iteration] / G.dx elif self.polarisation == 'y': if self.resistance != 0: - Ey[i, j, k] -= updatecoeffsE[ID[G.IDlookup[componentID], i, j, k], 4] * self.waveformvaluesJ[iteration] * (1 / (self.resistance * G.dx * G.dz)) + Ey[i, j, k] -= (updatecoeffsE[ID[G.IDlookup[componentID], i, j, k], 4] + * self.waveformvaluesJ[iteration] + * (1 / (self.resistance * G.dx * G.dz))) else: - Ey[i, j, k] = -1 * self.waveformvaluesJ[iteration] / G.dy + Ey[i, j, k] = - self.waveformvaluesJ[iteration] / G.dy elif self.polarisation == 'z': if self.resistance != 0: - Ez[i, j, k] -= updatecoeffsE[ID[G.IDlookup[componentID], i, j, k], 4] * self.waveformvaluesJ[iteration] * (1 / (self.resistance * G.dx * G.dy)) + Ez[i, j, k] -= (updatecoeffsE[ID[G.IDlookup[componentID], i, j, k], 4] + * self.waveformvaluesJ[iteration] + * (1 / (self.resistance * G.dx * G.dy))) else: - Ez[i, j, k] = -1 * self.waveformvaluesJ[iteration] / G.dz + Ez[i, j, k] = - self.waveformvaluesJ[iteration] / G.dz def create_material(self, G): """ @@ -173,13 +179,19 @@ class HertzianDipole(Source): componentID = 'E' + self.polarisation if self.polarisation == 'x': - Ex[i, j, k] -= updatecoeffsE[ID[G.IDlookup[componentID], i, j, k], 4] * self.waveformvaluesJ[iteration] * self.dl * (1 / (G.dx * G.dy * G.dz)) + Ex[i, j, k] -= (updatecoeffsE[ID[G.IDlookup[componentID], i, j, k], 4] + * self.waveformvaluesJ[iteration] + * self.dl * (1 / (G.dx * G.dy * G.dz))) elif self.polarisation == 'y': - Ey[i, j, k] -= updatecoeffsE[ID[G.IDlookup[componentID], i, j, k], 4] * self.waveformvaluesJ[iteration] * self.dl * (1 / (G.dx * G.dy * G.dz)) + Ey[i, j, k] -= (updatecoeffsE[ID[G.IDlookup[componentID], i, j, k], 4] + * self.waveformvaluesJ[iteration] + * self.dl * (1 / (G.dx * G.dy * G.dz))) elif self.polarisation == 'z': - Ez[i, j, k] -= updatecoeffsE[ID[G.IDlookup[componentID], i, j, k], 4] * self.waveformvaluesJ[iteration] * self.dl * (1 / (G.dx * G.dy * G.dz)) + Ez[i, j, k] -= (updatecoeffsE[ID[G.IDlookup[componentID], i, j, k], 4] + * self.waveformvaluesJ[iteration] + * self.dl * (1 / (G.dx * G.dy * G.dz))) class MagneticDipole(Source): @@ -206,17 +218,24 @@ class MagneticDipole(Source): componentID = 'H' + self.polarisation if self.polarisation == 'x': - Hx[i, j, k] -= updatecoeffsH[ID[G.IDlookup[componentID], i, j, k], 4] * self.waveformvaluesM[iteration] * (1 / (G.dx * G.dy * G.dz)) + Hx[i, j, k] -= (updatecoeffsH[ID[G.IDlookup[componentID], i, j, k], 4] + * self.waveformvaluesM[iteration] + * (1 / (G.dx * G.dy * G.dz))) elif self.polarisation == 'y': - Hy[i, j, k] -= updatecoeffsH[ID[G.IDlookup[componentID], i, j, k], 4] * self.waveformvaluesM[iteration] * (1 / (G.dx * G.dy * G.dz)) + Hy[i, j, k] -= (updatecoeffsH[ID[G.IDlookup[componentID], i, j, k], 4] + * self.waveformvaluesM[iteration] + * (1 / (G.dx * G.dy * G.dz))) elif self.polarisation == 'z': - Hz[i, j, k] -= updatecoeffsH[ID[G.IDlookup[componentID], i, j, k], 4] * self.waveformvaluesM[iteration] * (1 / (G.dx * G.dy * G.dz)) + Hz[i, j, k] -= (updatecoeffsH[ID[G.IDlookup[componentID], i, j, k], 4] + * self.waveformvaluesM[iteration] + * (1 / (G.dx * G.dy * G.dz))) def gpu_initialise_src_arrays(sources, G): - """Initialise arrays on GPU for source coordinates/polarisation, other source information, and source waveform values. + """Initialise arrays on GPU for source coordinates/polarisation, + other source information, and source waveform values. Args: sources (list): List of sources of one class, e.g. HertzianDipoles. @@ -264,7 +283,8 @@ def gpu_initialise_src_arrays(sources, G): class TransmissionLine(Source): """ A transmission line source is a one-dimensional transmission - line which is attached virtually to a grid cell. + line which is attached virtually to a grid cell. An example of this + type of model can be found in: https://doi.org/10.1109/8.277228 """ def __init__(self, G): @@ -304,7 +324,8 @@ class TransmissionLine(Source): def calculate_incident_V_I(self, G): """ Calculates the incident voltage and current with a long length - transmission line not connected to the main grid from: http://dx.doi.org/10.1002/mop.10415 + transmission line, initially not connected to the main grid. + This idea comes from: http://dx.doi.org/10.1002/mop.10415 Args: G (class): Grid class instance - holds essential parameters describing the model. @@ -327,7 +348,6 @@ class TransmissionLine(Source): """ h = (c * G.dt - self.dl) / (c * G.dt + self.dl) - self.voltage[0] = h * (self.voltage[1] - self.abcv0) + self.abcv1 self.abcv0 = self.voltage[0] self.abcv1 = self.voltage[1] @@ -341,7 +361,8 @@ class TransmissionLine(Source): """ # Update all the voltage values along the line - self.voltage[1:self.nl] -= self.resistance * (c * G.dt / self.dl) * (self.current[1:self.nl] - self.current[0:self.nl - 1]) + self.voltage[1:self.nl] -= (self.resistance * (c * G.dt / self.dl) + * (self.current[1:self.nl] - self.current[0:self.nl - 1])) # Update the voltage at the position of the one-way injector excitation self.voltage[self.srcpos] += (c * G.dt / self.dl) * self.waveformvaluesJ[iteration] @@ -358,10 +379,12 @@ class TransmissionLine(Source): """ # Update all the current values along the line - self.current[0:self.nl - 1] -= (1 / self.resistance) * (c * G.dt / self.dl) * (self.voltage[1:self.nl] - self.voltage[0:self.nl - 1]) + self.current[0:self.nl - 1] -= ((1 / self.resistance) * (c * G.dt / self.dl) + * (self.voltage[1:self.nl] - self.voltage[0:self.nl - 1])) # Update the current one cell before the position of the one-way injector excitation - self.current[self.srcpos - 1] += (1 / self.resistance) * (c * G.dt / self.dl) * self.waveformvaluesM[iteration] + self.current[self.srcpos - 1] += ((1 / self.resistance) * (c * G.dt / self.dl) + * self.waveformvaluesM[iteration]) def update_electric(self, iteration, updatecoeffsE, ID, Ex, Ey, Ez, G): """Updates electric field value in the main grid from voltage value in the transmission line. @@ -415,96 +438,4 @@ class TransmissionLine(Source): elif self.polarisation == 'z': self.current[self.antpos] = Iz(i, j, k, G.Hx, G.Hy, G.Hz, G) - self.update_current(iteration, G) - - -class PlaneWave(Source): - """A plane wave source. It uses a total-field/scattered-field (TF/SF) formulation.""" - - def __init__(self, G): - """ - Args: - G (class): Grid class instance - holds essential parameters describing the model. - """ - - super(Source, self).__init__() - - # Coordinates defining Huygen's surface - self.xs = 0 - self.xf = 0 - self.ys = 0 - self.yf = 0 - self.zs = 0 - self.zf = 0 - - # Spherical coordinates defining incident unit wavevector (k) - self.theta = 0 # 0 <= theta <= 180 - self.phi = 0 # 0 <= phi <= 360 - - # Angle that incident electric field makes with k cross z - self.psi = 0 # 0 <= psi <= 360 - - def calculate_origin(self, G): - """Calculate origin of TF/SF interface with incident wavefront.""" - - if self.theta >= 0 and self.theta <= 90: - if self.phi >= 0 and self.phi <= 90: - self.xcoordorigin = 0 - self.ycoordorigin = 0 - self.zcoordorigin = 0 - - elif self.phi > 90 and self.phi <= 180: - self.xcoordorigin = G.nx - self.ycoordorigin = 0 - self.zcoordorigin = 0 - - elif self.phi > 180 and self.phi <= 270: - self.xcoordorigin = G.nx - self.ycoordorigin = G.ny - self.zcoordorigin = 0 - - elif self.phi > 270 and self.phi <= 360: - self.xcoordorigin = 0 - self.ycoordorigin = G.ny - self.zcoordorigin = 0 - - elif self.theta > 90 and self.theta <= 180: - if self.phi >= 0 and self.phi <= 90: - self.xcoordorigin = 0 - self.ycoordorigin = 0 - self.zcoordorigin = G.nz - - elif self.phi > 90 and self.phi <= 180: - self.xcoordorigin = G.nx - self.ycoordorigin = 0 - self.zcoordorigin = G.nz - - elif self.phi > 180 and self.phi <= 270: - self.xcoordorigin = G.nx - self.ycoordorigin = G.ny - self.zcoordorigin = G.nz - - elif self.phi > 270 and self.phi <= 360: - self.xcoordorigin = 0 - self.ycoordorigin = G.ny - self.zcoordorigin = G.nz - - def calculate_vector_components(self): - """Calculate components of incident fields.""" - - self.theta = np.deg2rad(self.theta) - self.phi = np.deg2rad(self.phi) - self.psi = np.deg2rad(self.psi) - - # Components of incident unit wavevector - self.kx = np.sin(self.theta) * np.cos(self.phi) - self.ky = np.sin(self.theta) * np.sin(self.phi) - self.kz = np.cos(self.theta) - - # Components of incident field vectors - self.Exinc = np.cos(self.psi) * np.sin(self.phi) - np.sin(self.psi) * np.cos(self.theta) * np.cos(self.phi) - self.Eyinc = -np.cos(self.psi) * np.cos(self.phi) - np.sin(self.psi) * np.cos(self.theta) * np.sin(self.phi) - self.Ezinc = np.sin(self.psi) * np.sin(self.theta) - self.Hxinc = np.sin(self.psi) * np.sin(self.phi) + np.cos(self.psi) * np.cos(self.theta) * np.cos(self.phi) - self.Hyinc = -np.sin(self.psi) * np.cos(self.phi) + np.cos(self.psi) * np.cos(self.theta) * np.sin(self.phi) - self.Hzinc = -np.cos(self.psi) * np.sin(self.theta) + self.update_current(iteration, G) \ No newline at end of file