Changed importing of user-defined waveforms from excitation file to include a interpolation (see extra optional parameters for excitation file command in docs).

docs/source/input.rst

@@ -682,28 +682,29 @@ For example, to specify the normalised first derivate of a Gaussian waveform wit
 #excitation_file:
 -----------------
 
-Allows you to specify an ASCII file that contains columns of amplitude values that specify custom waveform shapes that can be used with sources in the model. The first row of each column must begin with a identifier string that will be used as the name of each waveform.
+Allows you to specify an ASCII file that contains columns of amplitude values that specify custom waveform shapes that can be used with sources in the model. The first row of each column must begin with a identifier string that will be used as the name of each waveform. Optionally, the first column of the file may contain a time vector of values (which must use the identifier ``time``) to interpolate the amplitude values of the waveform. If a time vector is not given, a vector of time values corresponding to the simulation time step and number of iterations will be used.
 
 If there are less amplitude values than the number of iterations that are going to be performed, the end of the sequence of amplitude values will be padded with zero values up to the number of iterations. If extra amplitude values are specified than needed then they are ignored. The syntax of the command is:
 
 .. code-block:: none
 
-    #excitation_file: file1
+    #excitation_file: file1 [str1 str2]
 
-``file1`` can be the name of the file containing the specified waveform in the same directory as the input file, or ``file`` can be the full path to the file containing the specified waveform (allowing you to specify any location).
+* ``file1`` can be the name of the file containing the specified waveform in the same directory as the input file, or ``file`` can be the full path to the file containing the specified waveform (allowing you to specify any location).
+* ``str1`` and ``str2`` are an optional parameter pair that allow values for ``kind`` and ``fill_value`` to be passed to the interpolation function (`scipy.interpolate.interp1d <https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.interp1d.html>`_). If they are not given the default values for the function will be used.
 
 For example, to specify the file ``my_waves.txt``, which contains two custom waveform shapes, use: ``#excitation_file: my_waves.txt``. The contents of the file ``my_waves.txt`` would take the form:
 
 .. code-block:: none
 
-    my_pulse1    my_pulse2
-    0           0
-    1.2e-6      0
-    1.3e-6      1e-1
-    5e-6        1.5e-1
-    ...         ...
-    ...         ...
-    ...         ...
+    time my_pulse1 my_pulse2
+    0 0 0
+    1.926e-12 1.2e-6 0
+    3.852e-12 1.3e-6 1.0e-1
+    5.778e-12 5.0e-6 1.5e-1
+    ...       ...    ...
+    ...       ...    ...
+    ...       ...    ...
 
 Then to use ``my_pulse1`` custom waveform shape with, for example, a z-polarised Hertzian dipole source:
 

gprMax/input_cmds_singleuse.py

@@ -18,11 +18,13 @@
 
 import os
 import decimal as d
+import inspect
 import sys
 
 from colorama import init, Fore, Style
 init()
 import numpy as np
+from scipy import interpolate
 
 from gprMax.constants import c
 from gprMax.constants import floattype
@@ -150,7 +152,7 @@ def process_singlecmds(singlecmds, G):
     if G.messages:
         print('Estimated memory (RAM) required: ~{}'.format(human_size(memestimate)))
 
-    # Time step CFL limit (use either 2D or 3D) and default PML thickness
+    # Time step CFL limit (either 2D or 3D); switch off appropriate PMLs for 2D
     if G.nx == 1:
         G.dt = 1 / (c * np.sqrt((1 / G.dy) * (1 / G.dy) + (1 / G.dz) * (1 / G.dz)))
         G.dimension = '2D'
@@ -216,7 +218,7 @@ def process_singlecmds(singlecmds, G):
     if singlecmds[cmd] is not None:
         tmp = singlecmds[cmd].split()
         if len(tmp) != 1 and len(tmp) != 6:
-            raise CmdInputError(cmd + ' requires either one or six parameters')
+            raise CmdInputError(cmd + ' requires either one or six parameter(s)')
         if len(tmp) == 1:
             for key in G.pmlthickness.keys():
                 G.pmlthickness[key] = int(tmp[0])
@@ -258,14 +260,26 @@ def process_singlecmds(singlecmds, G):
     cmd = '#excitation_file'
     if singlecmds[cmd] is not None:
         tmp = singlecmds[cmd].split()
-        if len(tmp) != 1:
-            raise CmdInputError(cmd + ' requires exactly one parameter')
+        if len(tmp) != 1 and len(tmp) != 3:
+            raise CmdInputError(cmd + ' requires either one or three parameter(s)')
         excitationfile = tmp[0]
 
+        # Optional parameters passed directly to scipy.interpolate.interp1d
+        kwargs = dict()
+        if len(tmp) > 1:
+            kwargs['kind'] = tmp[1]
+            kwargs['fill_value'] = tmp[2]
+        else:
+            args, varargs, keywords, defaults = inspect.getargspec(interpolate.interp1d)
+            kwargs = dict(zip(reversed(args), reversed(defaults)))
+
         # See if file exists at specified path and if not try input file directory
         if not os.path.isfile(excitationfile):
             excitationfile = os.path.abspath(os.path.join(G.inputdirectory, excitationfile))
 
+        if G.messages:
+            print('\nExcitation file: {}'.format(excitationfile))
+
         # Get waveform names
         with open(excitationfile, 'r') as f:
             waveformIDs = f.readline().split()
@@ -273,18 +287,37 @@ def process_singlecmds(singlecmds, G):
         # Read all waveform values into an array
         waveformvalues = np.loadtxt(excitationfile, skiprows=1, dtype=floattype)
 
+        # Time array (if specified) for interpolation, otherwise use simulation time
+        if waveformIDs[0].lower() == 'time':
+            waveformIDs = waveformIDs[1:]
+            waveformtime = waveformvalues[:, 0]
+            waveformvalues = waveformvalues[:, 1:]
+            timestr = 'user-defined time array'
+        else:
+            waveformtime = np.arange(0, G.timewindow + G.dt, G.dt)
+            timestr = 'simulation time array'
+
         for waveform in range(len(waveformIDs)):
             if any(x.ID == waveformIDs[waveform] for x in G.waveforms):
                 raise CmdInputError('Waveform with ID {} already exists'.format(waveformIDs[waveform]))
             w = Waveform()
             w.ID = waveformIDs[waveform]
             w.type = 'user'
-            if len(waveformvalues.shape) == 1:
-                w.uservalues = waveformvalues[:]
-            else:
-                w.uservalues = waveformvalues[:, waveform]
+
+            # Select correct column of waveform values depending on array shape
+            waveformvalues = waveformvalues[:] if len(waveformvalues.shape) == 1 else waveformvalues[:, waveform]
+
+            # Truncate waveform array if it is longer than time array
+            if len(waveformvalues) > len(waveformtime):
+                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)
+
+            # Interpolate waveform values
+            w.userfunc = interpolate.interp1d(waveformtime, waveformvalues, **kwargs)
 
             if G.messages:
-                print('User waveform {} created.'.format(w.ID))
+                print('User waveform {} created using {} and, if required, interpolation parameters (kind: {}, fill value: {}).'.format(w.ID, timestr, kwargs['kind'], kwargs['fill_value']))
 
             G.waveforms.append(w)

gprMax/waveforms.py

@@ -39,7 +39,7 @@ class Waveform(object):
         self.type = None
         self.amp = 1
         self.freq = None
-        self.uservalues = None
+        self.userfunc = None
         self.chi = 0
         self.zeta = 0
         self.delay = 0
@@ -66,7 +66,7 @@ class Waveform(object):
         """
 
         self.calculate_coefficients()
-        
+
         # Waveforms
         if self.type == 'gaussian':
             delay = time - self.chi
@@ -115,12 +115,7 @@ class Waveform(object):
                 ampvalue = 0
 
         elif self.type == 'user':
-            index = round_value(time / dt)
-            # Check to see if there are still user specified values and if not use zero
-            if index > len(self.uservalues) - 1:
-                ampvalue = 0
-            else:
-                ampvalue = self.uservalues[index]
+            ampvalue = self.userfunc(time)
 
         ampvalue *= self.amp