diff --git a/gprMax/grid.py b/gprMax/grid.py index 964d10ad..a9b9ef4c 100644 --- a/gprMax/grid.py +++ b/gprMax/grid.py @@ -165,51 +165,54 @@ def dispersion_analysis(G): """ # Physical phase velocity error (percentage); grid sampling density; material with maximum permittivity; maximum significant frequency - results = {'deltavp': False, 'N': False, 'material': False, 'maxfreq': []} + results = {'deltavp': False, 'N': False, 'waveform': True, 'material': False, 'maxfreq': []} # Find maximum significant frequency - for waveform in G.waveforms: + if G.waveforms: + for waveform in G.waveforms: + if waveform.type == 'sine' or waveform.type == 'contsine': + results['maxfreq'].append(4 * waveform.freq) - if waveform.type == 'sine' or waveform.type == 'contsine': - results['maxfreq'].append(4 * waveform.freq) - - elif waveform.type == 'impulse': - pass - - else: - # User-defined waveform - if waveform.type == 'user': - waveformvalues = waveform.uservalues - - # Built-in waveform - else: - time = np.linspace(0, 1, G.iterations) - time *= (G.iterations * G.dt) - waveformvalues = np.zeros(len(time)) - timeiter = np.nditer(time, flags=['c_index']) - - while not timeiter.finished: - waveformvalues[timeiter.index] = waveform.calculate_value(timeiter[0], G.dt) - timeiter.iternext() - - # Ensure source waveform is not being overly truncated before attempting any FFT - if np.abs(waveformvalues[-1]) < np.abs(np.amax(waveformvalues)) / 100: - # Calculate magnitude of frequency spectra of waveform - power = 10 * np.log10(np.abs(np.fft.fft(waveformvalues))**2) - freqs = np.fft.fftfreq(power.size, d=G.dt) - - # Shift powers so that frequency with maximum power is at zero decibels - power -= np.amax(power) - - # Get frequency for max power - freqmaxpower = np.where(np.isclose(power[1::], np.amax(power[1::])))[0][0] - - # Set maximum frequency to a threshold drop from maximum power, ignoring DC value - freq = np.where((np.amax(power[freqmaxpower::]) - power[freqmaxpower::]) > G.highestfreqthres)[0][0] + 1 - results['maxfreq'].append(freqs[freq]) + elif waveform.type == 'impulse': + pass else: - results['waveformID'] = waveform.ID + # User-defined waveform + if waveform.type == 'user': + waveformvalues = waveform.uservalues + + # Built-in waveform + else: + time = np.linspace(0, 1, G.iterations) + time *= (G.iterations * G.dt) + waveformvalues = np.zeros(len(time)) + timeiter = np.nditer(time, flags=['c_index']) + + while not timeiter.finished: + waveformvalues[timeiter.index] = waveform.calculate_value(timeiter[0], G.dt) + timeiter.iternext() + + # Ensure source waveform is not being overly truncated before attempting any FFT + if np.abs(waveformvalues[-1]) < np.abs(np.amax(waveformvalues)) / 100: + # Calculate magnitude of frequency spectra of waveform + power = 10 * np.log10(np.abs(np.fft.fft(waveformvalues))**2) + freqs = np.fft.fftfreq(power.size, d=G.dt) + + # Shift powers so that frequency with maximum power is at zero decibels + power -= np.amax(power) + + # Get frequency for max power + freqmaxpower = np.where(np.isclose(power[1::], np.amax(power[1::])))[0][0] + + # Set maximum frequency to a threshold drop from maximum power, ignoring DC value + freq = np.where((np.amax(power[freqmaxpower::]) - power[freqmaxpower::]) > G.highestfreqthres)[0][0] + 1 + results['maxfreq'].append(freqs[freq]) + + # If waveform is truncated don't do any further analysis + else: + results['waveform'] = False + else: + results['waveform'] = False if results['maxfreq']: results['maxfreq'] = max(results['maxfreq']) diff --git a/gprMax/model_build_run.py b/gprMax/model_build_run.py index 88ff8124..d5694b90 100644 --- a/gprMax/model_build_run.py +++ b/gprMax/model_build_run.py @@ -172,8 +172,8 @@ def run_model(args, currentmodelrun, numbermodelruns, inputfile, usernamespace): # Check to see if numerical dispersion might be a problem results = dispersion_analysis(G) - if 'waveformID' in results: - print(Fore.RED + "\nWARNING: Numerical dispersion analysis not carried out as duration of waveform '{}' means it does not fit within specified time window and is therefore being truncated.".format(results['waveformID']) + Style.RESET_ALL) + if not results['waveform']: + print(Fore.RED + "\nWARNING: Numerical dispersion analysis not carried out as either no waveform detected or waveform does not fit within specified time window and is therefore being truncated." + Style.RESET_ALL) elif results['N'] < G.mingridsampling: raise GeneralError("Non-physical wave propagation: Material '{}' has wavelength sampled by {} cells, less than required minimum for physical wave propagation. Maximum significant frequency estimated as {:g}Hz".format(results['material'].ID, results['N'], results['maxfreq'])) elif results['deltavp'] and np.abs(results['deltavp']) > G.maxnumericaldisp: