From 3bf30d069e862082b2aa6310ce40ecb43580311c Mon Sep 17 00:00:00 2001 From: Craig Warren Date: Thu, 26 Jan 2017 12:15:04 +0000 Subject: [PATCH] Corrected issues with dispersion analysis for high values of dielectric, and for pulses with low power at low frequencies (affected finding highest frequency). --- gprMax/grid.py | 16 +++++++++++++--- 1 file changed, 13 insertions(+), 3 deletions(-) diff --git a/gprMax/grid.py b/gprMax/grid.py index 85a021f0..17064531 100644 --- a/gprMax/grid.py +++ b/gprMax/grid.py @@ -21,6 +21,7 @@ from collections import OrderedDict from colorama import init, Fore, Style init() import numpy as np +np.seterr(invalid='raise') from gprMax.constants import c, floattype, complextype from gprMax.materials import Material @@ -193,9 +194,12 @@ def dispersion_analysis(G): # Shift powers so that frequency with maximum power is at zero decibels power -= np.amax(power) + + # Get frequency for max power + freqmaxpower = np.where(power[1::] == np.amax(power[1::]))[0][0] - # Set maximum frequency to -60dB from maximum power, ignoring DC value - freq = np.where((np.amax(power[1::]) - power[1::]) > G.highestfreqthres)[0][0] + 1 + # 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 maxfreqs.append(freqs[freq]) else: @@ -240,8 +244,14 @@ def dispersion_analysis(G): # Grid sampling density results['N'] = minwavelength / delta + # If there is an invalid inverse sine value set to maximum, i.e. pi/2 + try: + tmp = np.arcsin((1 / S) * np.sin((np.pi * S) / results['N'])) + except FloatingPointError: + tmp = np.pi / 2 + # Numerical phase velocity - vp = np.pi / (results['N'] * np.arcsin((1 / S) * np.sin((np.pi * S) / results['N']))) + vp = np.pi / (results['N'] * tmp) # Physical phase velocity error (percentage) results['deltavp'] = (((vp * c) - c) / c) * 100