diff --git a/gprMax/fields_outputs.py b/gprMax/fields_outputs.py index 5ed0cdba..1cf8c445 100644 --- a/gprMax/fields_outputs.py +++ b/gprMax/fields_outputs.py @@ -67,7 +67,7 @@ __global__ void store_outputs(int NRX, int iteration, const int* __restrict__ rx // NRX: Total number of receivers in the model // rxs: Array to store field components for receivers - rows are field components; columns are iterations; pages are receivers // E, H: Access to field component arrays - + // Obtain the linear index corresponding to the current thread and use for each receiver int rx = blockIdx.x * blockDim.x + threadIdx.x; diff --git a/gprMax/fields_updates_gpu.py b/gprMax/fields_updates_gpu.py index 781915e0..e3d4af4f 100644 --- a/gprMax/fields_updates_gpu.py +++ b/gprMax/fields_updates_gpu.py @@ -29,7 +29,7 @@ kernels_template_fields = Template(""" #define INDEX4D_ID(p, i, j, k) (p)*($NX_ID)*($NY_ID)*($NZ_ID)+(i)*($NY_ID)*($NZ_ID)+(j)*($NZ_ID)+(k) #define INDEX4D_T(p, i, j, k) (p)*($NX_T)*($NY_T)*($NZ_T)+(i)*($NY_T)*($NZ_T)+(j)*($NZ_T)+(k) -// Material coefficients (read-only) in constant memory (64KB) +// Material coefficients (read-only) in constant memory (64KB)_ __device__ __constant__ $REAL updatecoeffsE[$N_updatecoeffsE]; __device__ __constant__ $REAL updatecoeffsH[$N_updatecoeffsH]; @@ -44,7 +44,7 @@ __global__ void update_e(int NX, int NY, int NZ, const unsigned int* __restrict_ // Args: // NX, NY, NZ: Number of cells of the model domain // ID, E, H: Access to ID and field component arrays - + // Obtain the linear index corresponding to the current thread int idx = blockIdx.x * blockDim.x + threadIdx.x; @@ -63,13 +63,13 @@ __global__ void update_e(int NX, int NY, int NZ, const unsigned int* __restrict_ int materialEx = ID[INDEX4D_ID(0,i_ID,j_ID,k_ID)]; Ex[INDEX3D_FIELDS(i,j,k)] = updatecoeffsE[INDEX2D_MAT(materialEx,0)] * Ex[INDEX3D_FIELDS(i,j,k)] + updatecoeffsE[INDEX2D_MAT(materialEx,2)] * (Hz[INDEX3D_FIELDS(i,j,k)] - Hz[INDEX3D_FIELDS(i,j-1,k)]) - updatecoeffsE[INDEX2D_MAT(materialEx,3)] * (Hy[INDEX3D_FIELDS(i,j,k)] - Hy[INDEX3D_FIELDS(i,j,k-1)]); } - + // Ey component if ((NX != 1 || NZ != 1) && i > 0 && i < NX && j >= 0 && j < NY && k > 0 && k < NZ) { int materialEy = ID[INDEX4D_ID(1,i_ID,j_ID,k_ID)]; Ey[INDEX3D_FIELDS(i,j,k)] = updatecoeffsE[INDEX2D_MAT(materialEy,0)] * Ey[INDEX3D_FIELDS(i,j,k)] + updatecoeffsE[INDEX2D_MAT(materialEy,3)] * (Hx[INDEX3D_FIELDS(i,j,k)] - Hx[INDEX3D_FIELDS(i,j,k-1)]) - updatecoeffsE[INDEX2D_MAT(materialEy,1)] * (Hz[INDEX3D_FIELDS(i,j,k)] - Hz[INDEX3D_FIELDS(i-1,j,k)]); } - + // Ez component if ((NX != 1 || NY != 1) && i > 0 && i < NX && j > 0 && j < NY && k >= 0 && k < NZ) { int materialEz = ID[INDEX4D_ID(2,i_ID,j_ID,k_ID)]; @@ -89,7 +89,7 @@ __global__ void update_e_dispersive_A(int NX, int NY, int NZ, int MAXPOLES, cons // NX, NY, NZ: Number of cells of the model domain // MAXPOLES: Maximum number of dispersive material poles present in model // updatedispersivecoeffs, T, ID, E, H: Access to update coefficients, dispersive, ID and field component arrays - + // Obtain the linear index corresponding to the current thread int idx = blockIdx.x * blockDim.x + threadIdx.x; @@ -118,7 +118,7 @@ __global__ void update_e_dispersive_A(int NX, int NY, int NZ, int MAXPOLES, cons } Ex[INDEX3D_FIELDS(i,j,k)] = updatecoeffsE[INDEX2D_MAT(materialEx,0)] * Ex[INDEX3D_FIELDS(i,j,k)] + updatecoeffsE[INDEX2D_MAT(materialEx,2)] * (Hz[INDEX3D_FIELDS(i,j,k)] - Hz[INDEX3D_FIELDS(i,j-1,k)]) - updatecoeffsE[INDEX2D_MAT(materialEx,3)] * (Hy[INDEX3D_FIELDS(i,j,k)] - Hy[INDEX3D_FIELDS(i,j,k-1)]) - updatecoeffsE[INDEX2D_MAT(materialEx,4)] * phi; } - + // Ey component if ((NX != 1 || NZ != 1) && i > 0 && i < NX && j >= 0 && j < NY && k > 0 && k < NZ) { int materialEy = ID[INDEX4D_ID(1,i_ID,j_ID,k_ID)]; @@ -129,7 +129,7 @@ __global__ void update_e_dispersive_A(int NX, int NY, int NZ, int MAXPOLES, cons } Ey[INDEX3D_FIELDS(i,j,k)] = updatecoeffsE[INDEX2D_MAT(materialEy,0)] * Ey[INDEX3D_FIELDS(i,j,k)] + updatecoeffsE[INDEX2D_MAT(materialEy,3)] * (Hx[INDEX3D_FIELDS(i,j,k)] - Hx[INDEX3D_FIELDS(i,j,k-1)]) - updatecoeffsE[INDEX2D_MAT(materialEy,1)] * (Hz[INDEX3D_FIELDS(i,j,k)] - Hz[INDEX3D_FIELDS(i-1,j,k)]) - updatecoeffsE[INDEX2D_MAT(materialEy,4)] * phi; } - + // Ez component if ((NX != 1 || NY != 1) && i > 0 && i < NX && j > 0 && j < NY && k >= 0 && k < NZ) { int materialEz = ID[INDEX4D_ID(2,i_ID,j_ID,k_ID)]; @@ -150,7 +150,7 @@ __global__ void update_e_dispersive_B(int NX, int NY, int NZ, int MAXPOLES, cons // NX, NY, NZ: Number of cells of the model domain // MAXPOLES: Maximum number of dispersive material poles present in model // updatedispersivecoeffs, T, ID, E, H: Access to update coefficients, dispersive, ID and field component arrays - + // Obtain the linear index corresponding to the current thread int idx = blockIdx.x * blockDim.x + threadIdx.x; @@ -176,7 +176,7 @@ __global__ void update_e_dispersive_B(int NX, int NY, int NZ, int MAXPOLES, cons Tx[INDEX4D_T(pole,i_T,j_T,k_T)] = Tx[INDEX4D_T(pole,i_T,j_T,k_T)] - updatecoeffsdispersive[INDEX2D_MATDISP(materialEx,2+(pole*3))] * Ex[INDEX3D_FIELDS(i,j,k)]; } } - + // Ey component if ((NX != 1 || NZ != 1) && i > 0 && i < NX && j >= 0 && j < NY && k > 0 && k < NZ) { int materialEy = ID[INDEX4D_ID(1,i_ID,j_ID,k_ID)]; @@ -184,7 +184,7 @@ __global__ void update_e_dispersive_B(int NX, int NY, int NZ, int MAXPOLES, cons Ty[INDEX4D_T(pole,i_T,j_T,k_T)] = Ty[INDEX4D_T(pole,i_T,j_T,k_T)] - updatecoeffsdispersive[INDEX2D_MATDISP(materialEy,2+(pole*3))] * Ey[INDEX3D_FIELDS(i,j,k)]; } } - + // Ez component if ((NX != 1 || NY != 1) && i > 0 && i < NX && j > 0 && j < NY && k >= 0 && k < NZ) { int materialEz = ID[INDEX4D_ID(2,i_ID,j_ID,k_ID)]; @@ -206,7 +206,7 @@ __global__ void update_h(int NX, int NY, int NZ, const unsigned int* __restrict_ // Args: // NX, NY, NZ: Number of cells of the model domain // ID, E, H: Access to ID and field component arrays - + // Obtain the linear index corresponding to the current thread int idx = blockIdx.x * blockDim.x + threadIdx.x; @@ -225,13 +225,13 @@ __global__ void update_h(int NX, int NY, int NZ, const unsigned int* __restrict_ int materialHx = ID[INDEX4D_ID(3,i_ID,j_ID,k_ID)]; Hx[INDEX3D_FIELDS(i,j,k)] = updatecoeffsH[INDEX2D_MAT(materialHx,0)] * Hx[INDEX3D_FIELDS(i,j,k)] - updatecoeffsH[INDEX2D_MAT(materialHx,2)] * (Ez[INDEX3D_FIELDS(i,j+1,k)] - Ez[INDEX3D_FIELDS(i,j,k)]) + updatecoeffsH[INDEX2D_MAT(materialHx,3)] * (Ey[INDEX3D_FIELDS(i,j,k+1)] - Ey[INDEX3D_FIELDS(i,j,k)]); } - + // Hy component if (NY != 1 && i >= 0 && i < NX && j > 0 && j < NY && k >= 0 && k < NZ) { int materialHy = ID[INDEX4D_ID(4,i_ID,j_ID,k_ID)]; Hy[INDEX3D_FIELDS(i,j,k)] = updatecoeffsH[INDEX2D_MAT(materialHy,0)] * Hy[INDEX3D_FIELDS(i,j,k)] - updatecoeffsH[INDEX2D_MAT(materialHy,3)] * (Ex[INDEX3D_FIELDS(i,j,k+1)] - Ex[INDEX3D_FIELDS(i,j,k)]) + updatecoeffsH[INDEX2D_MAT(materialHy,1)] * (Ez[INDEX3D_FIELDS(i+1,j,k)] - Ez[INDEX3D_FIELDS(i,j,k)]); } - + // Hz component if (NZ != 1 && i >= 0 && i < NX && j >= 0 && j < NY && k > 0 && k < NZ) { int materialHz = ID[INDEX4D_ID(5,i_ID,j_ID,k_ID)]; diff --git a/gprMax/fractals.py b/gprMax/fractals.py index 1c2945c3..19a2da2e 100644 --- a/gprMax/fractals.py +++ b/gprMax/fractals.py @@ -100,7 +100,7 @@ class FractalSurface(object): fractalmax = np.amax(self.fractalsurface) fractalrange = fractalmax - fractalmin self.fractalsurface = self.fractalsurface * ((self.fractalrange[1] - self.fractalrange[0]) / fractalrange) \ - + self.fractalrange[0] - ((self.fractalrange[1] - self.fractalrange[0]) / fractalrange) * fractalmin + + self.fractalrange[0] - ((self.fractalrange[1] - self.fractalrange[0]) / fractalrange) * fractalmin class FractalVolume(object): diff --git a/gprMax/geometry_outputs.py b/gprMax/geometry_outputs.py index 65760bcc..d152525b 100644 --- a/gprMax/geometry_outputs.py +++ b/gprMax/geometry_outputs.py @@ -218,25 +218,25 @@ class GeometryView(object): z_materials = np.zeros((n_z_lines), dtype=np.uint32) define_fine_geometry(self.nx, - self.ny, - self.nz, - self.xs, - self.xf, - self.ys, - self.yf, - self.zs, - self.zf, - G.dx, - G.dy, - G.dz, - G.ID, - points, - x_lines, - x_materials, - y_lines, - y_materials, - z_lines, - z_materials) + self.ny, + self.nz, + self.xs, + self.xf, + self.ys, + self.yf, + self.zs, + self.zf, + G.dx, + G.dy, + G.dz, + G.ID, + points, + x_lines, + x_materials, + y_lines, + y_materials, + z_lines, + z_materials) # Write point data datasize = points.nbytes diff --git a/gprMax/gprMax.py b/gprMax/gprMax.py index ccbae6df..1435c392 100644 --- a/gprMax/gprMax.py +++ b/gprMax/gprMax.py @@ -70,18 +70,18 @@ def main(): def api( - inputfile, - n=1, - task=None, - restart=None, - mpi=False, - gpu=None, - benchmark=False, - geometry_only=False, - geometry_fixed=False, - write_processed=False, - opt_taguchi=False - ): + inputfile, + n=1, + task=None, + restart=None, + mpi=False, + gpu=None, + benchmark=False, + geometry_only=False, + geometry_fixed=False, + write_processed=False, + opt_taguchi=False +): """If installed as a module this is the entry point.""" # Print gprMax logo, version, and licencing/copyright information @@ -165,7 +165,7 @@ def run_main(args): # Process for simulation with Taguchi optimisation # #################################################### elif args.opt_taguchi: - if args.mpi_worker: # Special case for MPI spawned workers - they do not need to enter the Taguchi optimisation mode + if args.mpi_worker: # Special case for MPI spawned workers - they do not need to enter the Taguchi optimisation mode run_mpi_sim(args, inputfile, usernamespace) else: from gprMax.optimisation_taguchi import run_opt_sim @@ -396,9 +396,9 @@ def run_mpi_sim(args, inputfile, usernamespace, optparams=None): # Connect to parent try: - comm = MPI.Comm.Get_parent() # get MPI communicator object + comm = MPI.Comm.Get_parent() # get MPI communicator object rank = comm.Get_rank() # rank of this process - except: + except ValueError: raise ValueError('Could not connect to parent') # Ask for work until stop sentinel diff --git a/gprMax/grid.py b/gprMax/grid.py index 84db0654..ee93fd18 100644 --- a/gprMax/grid.py +++ b/gprMax/grid.py @@ -50,10 +50,10 @@ class Grid(object): self.grid = grid def n_edges(self): - l = self.nx - m = self.ny - n = self.nz - e = (l * m * (n - 1)) + (m * n * (l - 1)) + (l * n * (m - 1)) + i = self.nx + j = self.ny + k = self.nz + e = (i * j * (k - 1)) + (j * k * (i - 1)) + (i * k * (j - 1)) return e def n_nodes(self): @@ -239,8 +239,8 @@ def dispersion_analysis(G): # Built-in waveform else: - # Time to analyse waveform - 4*pulse_width as using entire - # time window can result in demanding FFT + # Time to analyse waveform - 4*pulse_width as using entire + # time window can result in demanding FFT waveform.calculate_coefficients() iterations = round_value(4 * waveform.chi / G.dt) if iterations > G.iterations: @@ -261,7 +261,7 @@ def dispersion_analysis(G): try: freqthres = np.where(power[freqmaxpower:] < -G.highestfreqthres)[0][0] + freqmaxpower results['maxfreq'].append(freqs[freqthres]) - except: + except ValueError: results['error'] = 'unable to calculate maximum power from waveform, most likely due to undersampling.' # If waveform is truncated don't do any further analysis diff --git a/gprMax/input_cmd_funcs.py b/gprMax/input_cmd_funcs.py index c5c3e6e8..5b1078db 100644 --- a/gprMax/input_cmd_funcs.py +++ b/gprMax/input_cmd_funcs.py @@ -413,7 +413,7 @@ def cylinder(x1, y1, z1, x2, y2, z2, radius, material, averaging='', rotate90ori def cylindrical_sector(axis, ctr1, ctr2, t1, t2, radius, - startingangle, sweptangle, material, averaging=''): + startingangle, sweptangle, material, averaging=''): """Prints the gprMax #cylindrical_sector command. Args: @@ -471,7 +471,7 @@ def waveform(shape, amplitude, frequency, identifier): def hertzian_dipole(polarisation, f1, f2, f3, identifier, - t0=None, t_remove=None, dxdy=None, rotate90origin=()): + t0=None, t_remove=None, dxdy=None, rotate90origin=()): """Prints the #hertzian_dipole: polarisation, f1, f2, f3, identifier, [t0, t_remove] Args: @@ -545,7 +545,7 @@ def magnetic_dipole(polarisation, f1, f2, f3, identifier, def voltage_source(polarisation, f1, f2, f3, resistance, identifier, - t0=None, t_remove=None, dxdy=None, rotate90origin=()): + t0=None, t_remove=None, dxdy=None, rotate90origin=()): """Prints the #voltage_source: polarisation, f1, f2, f3, resistance, identifier, [t0, t_remove] Args: @@ -583,7 +583,7 @@ def voltage_source(polarisation, f1, f2, f3, resistance, identifier, def transmission_line(polarisation, f1, f2, f3, resistance, identifier, - t0=None, t_remove=None, dxdy=None, rotate90origin=()): + t0=None, t_remove=None, dxdy=None, rotate90origin=()): """Prints the #transmission_line: polarisation, f1, f2, f3, resistance, identifier, [t0, t_remove] Args: diff --git a/gprMax/input_cmds_geometry.py b/gprMax/input_cmds_geometry.py index 7eab34b8..9dc75d3c 100644 --- a/gprMax/input_cmds_geometry.py +++ b/gprMax/input_cmds_geometry.py @@ -877,10 +877,7 @@ def process_geometrycmds(geometry, G): volume.nbins = nbins volume.seed = seed volume.weighting = np.array([float(tmp[8]), float(tmp[9]), float(tmp[10])]) - try: - volume.averaging = averagefractalbox - except: - pass + volume.averaging = averagefractalbox if G.messages: if volume.averaging: diff --git a/gprMax/input_cmds_multiuse.py b/gprMax/input_cmds_multiuse.py index 3bb77a47..3ba653cc 100644 --- a/gprMax/input_cmds_multiuse.py +++ b/gprMax/input_cmds_multiuse.py @@ -499,7 +499,7 @@ def process_multicmds(multicmds, G): try: time = int(tmp[9]) # If real floating point value given - except: + except ValueError: time = float(tmp[9]) if time > 0: time = round_value((time / G.dt)) + 1 @@ -587,13 +587,8 @@ def process_multicmds(multicmds, G): material.type = 'debye' material.poles = poles material.averagable = False - # for pole in range(1, 2 * poles, 2): - # if float(tmp[pole]) > 0 and float(tmp[pole + 1]) > G.dt: - # material.deltaer.append(float(tmp[pole])) - # material.tau.append(float(tmp[pole + 1])) - # else: - # raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' requires positive values for the permittivity difference, and relaxation times that are greater than the time step for the model.') for pole in range(1, 2 * poles, 2): + # N.B Not checking if relaxation times are greater than time-step if float(tmp[pole]) > 0: material.deltaer.append(float(tmp[pole])) material.tau.append(float(tmp[pole + 1])) diff --git a/gprMax/input_cmds_singleuse.py b/gprMax/input_cmds_singleuse.py index 53bb3ee5..73ccab35 100644 --- a/gprMax/input_cmds_singleuse.py +++ b/gprMax/input_cmds_singleuse.py @@ -203,7 +203,7 @@ def process_singlecmds(singlecmds, G): G.timewindow = (tmp - 1) * G.dt G.iterations = tmp # If real floating point value given - except: + except ValueError: tmp = float(tmp) if tmp > 0: G.timewindow = tmp @@ -312,7 +312,7 @@ def process_singlecmds(singlecmds, G): waveformvalues = waveformvalues[:len(waveformtime)] # Zero-pad end of waveform array if it is shorter than time array elif len(waveformvalues) < len(waveformtime): - waveformvalues = np.lib.pad(waveformvalues, (0,len(waveformtime)-len(waveformvalues)), 'constant', constant_values=0) + waveformvalues = np.lib.pad(waveformvalues, (0, len(waveformtime) - len(waveformvalues)), 'constant', constant_values=0) # Interpolate waveform values w.userfunc = interpolate.interp1d(waveformtime, waveformvalues, **kwargs) diff --git a/gprMax/materials.py b/gprMax/materials.py index bb7687dd..6b26ad87 100644 --- a/gprMax/materials.py +++ b/gprMax/materials.py @@ -148,14 +148,14 @@ class Material(object): Args: freq (float): Frequency used to calculate complex relative permittivity. - + Returns: er (float): Complex relative permittivity. """ - - # This will be permittivity at infinite frequency if the material is dispersive + + # Permittivity at infinite frequency if the material is dispersive er = self.er - + if self.poles > 0: w = 2 * np.pi * freq er += self.se / (w * e0) @@ -166,10 +166,11 @@ class Material(object): for pole in range(self.poles): er += (self.deltaer[pole] * self.tau[pole]**2) / (self.tau[pole]**2 + 2j * w * self.alpha[pole] - w**2) elif 'drude' in self.type: + ersum = 0 for pole in range(self.poles): ersum += self.tau[pole]**2 / (w**2 - 1j * w * self.alpha[pole]) er -= ersum - + return er diff --git a/gprMax/receivers.py b/gprMax/receivers.py index 9299895b..a19f3362 100644 --- a/gprMax/receivers.py +++ b/gprMax/receivers.py @@ -44,7 +44,7 @@ class Rx(object): def gpu_initialise_rx_arrays(G): """Initialise arrays on GPU for receiver coordinates and to store field components for receivers. - + Args: G (class): Grid class instance - holds essential parameters describing the model. """ diff --git a/gprMax/source_updates_gpu.py b/gprMax/source_updates_gpu.py index 7b29aaf9..63292fae 100644 --- a/gprMax/source_updates_gpu.py +++ b/gprMax/source_updates_gpu.py @@ -52,7 +52,7 @@ __global__ void update_hertzian_dipole(int NHERTZDIPOLE, int iteration, $REAL dx int src = blockIdx.x * blockDim.x + threadIdx.x; if (src < NHERTZDIPOLE) { - + $REAL dl; int i, j, k, polarisation; @@ -104,7 +104,7 @@ __global__ void update_magnetic_dipole(int NMAGDIPOLE, int iteration, $REAL dx, int src = blockIdx.x * blockDim.x + threadIdx.x; if (src < NMAGDIPOLE) { - + int i, j, k, polarisation; i = srcinfo1[INDEX2D_SRCINFO(src,0)]; @@ -154,10 +154,10 @@ __global__ void update_voltage_source(int NVOLTSRC, int iteration, $REAL dx, $RE int src = blockIdx.x * blockDim.x + threadIdx.x; if (src < NVOLTSRC) { - + $REAL resistance; int i, j, k, polarisation; - + i = srcinfo1[INDEX2D_SRCINFO(src,0)]; j = srcinfo1[INDEX2D_SRCINFO(src,1)]; k = srcinfo1[INDEX2D_SRCINFO(src,2)]; diff --git a/gprMax/sources.py b/gprMax/sources.py index 5dd0179e..5881a018 100644 --- a/gprMax/sources.py +++ b/gprMax/sources.py @@ -217,17 +217,17 @@ class MagneticDipole(Source): def gpu_initialise_src_arrays(sources, G): """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. G (class): Grid class instance - holds essential parameters describing the model. - + Returns: srcinfo1_gpu (int): numpy array of source cell coordinates and polarisation information. srcinfo2_gpu (float): numpy array of other source information, e.g. length, resistance etc... srcwaves_gpu (float): numpy array of source waveform values. """ - + import pycuda.gpuarray as gpuarray srcinfo1 = np.zeros((len(sources), 4), dtype=np.int32) @@ -237,14 +237,14 @@ def gpu_initialise_src_arrays(sources, G): srcinfo1[i, 0] = src.xcoord srcinfo1[i, 1] = src.ycoord srcinfo1[i, 2] = src.zcoord - + if src.polarisation == 'x': srcinfo1[i, 3] = 0 elif src.polarisation == 'y': srcinfo1[i, 3] = 1 elif src.polarisation == 'z': srcinfo1[i, 3] = 2 - + if src.__class__.__name__ == 'HertzianDipole': srcinfo2[i] = src.dl srcwaves[i, :] = src.waveformvaluesJ @@ -253,7 +253,7 @@ def gpu_initialise_src_arrays(sources, G): srcwaves[i, :] = src.waveformvaluesJ elif src.__class__.__name__ == 'MagneticDipole': srcwaves[i, :] = src.waveformvaluesM - + srcinfo1_gpu = gpuarray.to_gpu(srcinfo1) srcinfo2_gpu = gpuarray.to_gpu(srcinfo2) srcwaves_gpu = gpuarray.to_gpu(srcwaves) diff --git a/gprMax/utilities.py b/gprMax/utilities.py index d26f68e4..23cd8294 100644 --- a/gprMax/utilities.py +++ b/gprMax/utilities.py @@ -152,7 +152,7 @@ def fft_power(waveform, dt): """ # Calculate magnitude of frequency spectra of waveform (ignore warning from taking a log of any zero values) - with np.errstate(divide='ignore'): # + with np.errstate(divide='ignore'): power = 10 * np.log10(np.abs(np.fft.fft(waveform))**2) # Replace any NaNs or Infs from zero division