diff --git a/.github/workflows/pytest.yml b/.github/workflows/pytest.yml
index d4ef3daa..0047694f 100644
--- a/.github/workflows/pytest.yml
+++ b/.github/workflows/pytest.yml
@@ -8,11 +8,11 @@ jobs:
     runs-on: ubuntu-latest
     steps:
       - name: checkout repository
-        uses: actions/checkout@v3
+        uses: actions/checkout@v4
         with:
           path: k-Wave-python
       - name: checkout reference k-Wave repo
-        uses: actions/checkout@v3
+        uses: actions/checkout@v4
         with:
           repository: ucl-bug/k-wave
           path: k-wave
@@ -50,7 +50,7 @@ jobs:
         python-version: [ "3.8", "3.9", "3.10", "3.11" ]
     runs-on: ${{matrix.os}}
     steps:
-      - uses: actions/checkout@v3
+      - uses: actions/checkout@v4
       - name: Download collectedValues
         uses: actions/download-artifact@v3
         with:
diff --git a/.github/workflows/test_example.yml b/.github/workflows/test_example.yml
index f8177b0a..c9b6c194 100644
--- a/.github/workflows/test_example.yml
+++ b/.github/workflows/test_example.yml
@@ -10,7 +10,7 @@ jobs:
         python-version: [ "3.8", "3.9", "3.10", "3.11"]
     runs-on: ${{matrix.os}}
     steps:
-      - uses: actions/checkout@v3
+      - uses: actions/checkout@v4
       - name: Set up Python ${{ matrix.python-version }}
         uses: actions/setup-python@v3
         with:
diff --git a/Dockerfile b/Dockerfile
index f23f9b5d..5a9d1027 100644
--- a/Dockerfile
+++ b/Dockerfile
@@ -4,7 +4,8 @@ RUN apt-get update
 RUN apt-get install -y gfortran libopenblas-dev liblapack-dev libgl1-mesa-glx libglib2.0-0 libgl1-mesa-glx libglib2.0-0 git
 RUN pip install --upgrade pip
 COPY pyproject.toml .
-COPY README.md .
+#COPY README.md .
+COPY docs/ docs
 COPY LICENSE .
 COPY kwave/ kwave
 RUN pip install '.[test]'
diff --git a/docs/README.md b/docs/README.md
index 5a13e691..ef9ec821 100644
--- a/docs/README.md
+++ b/docs/README.md
@@ -1,8 +1,8 @@
 # k-Wave-python
 [![Documentation Status](https://readthedocs.org/projects/k-wave-python/badge/?version=latest)](https://k-wave-python.readthedocs.io/en/latest/?badge=latest)
 
-This project is a Python implementation of most of the [MATLAB toolbox k-Wave](http://www.k-wave.org/) as well as an
-interface to the pre-compiled v1.3 of k-Wave simulation binaries which support NVIDIA sm 5.0 (Maxwell) to sm 9.0a (Hopper) GPUs.
+This project is a Python implementation of v1.4.0 of the [MATLAB toolbox k-Wave](http://www.k-wave.org/) as well as an
+interface to the pre-compiled v1.3 of k-Wave simulation binaries, which support NVIDIA sm 5.0 (Maxwell) to sm 9.0a (Hopper) GPUs.
 
 ## Mission
 
@@ -31,7 +31,7 @@ After installation, run the B-mode reconstruction example in the `examples` dire
 ```bash
 git clone https://github.com/waltsims/k-wave-python
 cd k-wave-python
-git checkout v0.2.1
+git checkout v0.3.0
 pip install '.[example]' 
 python3 examples/bmode_reconstruction_example.py
 ```
@@ -46,7 +46,7 @@ This example file steps through the process of:
 This example expects an NVIDIA GPU by default to simulate with k-Wave.
 
 To test the reconstruction on a machine with a GPU,
-set `RUN_SIMULATION` [on line 28 of `bmode_reconstruction_example.py`](https://github.com/waltsims/k-wave-python/blob/master/examples/bmode_reconstruction_example.py#L28)
+set `RUN_SIMULATION` [on line 30 of `bmode_reconstruction_example.py`](https://github.com/waltsims/k-wave-python/blob/master/examples/bmode_reconstruction_example.py#L30)
 to `True`, and the example will run without the pre-computed data.
 
 ## Development
@@ -68,4 +68,4 @@ The documentation for k-wave-python can be found [here](https://k-wave-python.re
 
 ## Contact
 
-e-mail [walter.simson@tum.de](mailto:walter.simson@tum.de).
+e-mail [wsimson@stanford.edu](mailto:wsimson@stanford.edu).
diff --git a/examples/bmode_reconstruction_example.py b/examples/bmode_reconstruction_example.py
index 181c8dcb..b52318b1 100644
--- a/examples/bmode_reconstruction_example.py
+++ b/examples/bmode_reconstruction_example.py
@@ -18,7 +18,8 @@
 from kwave.utils.dotdictionary import dotdict
 from kwave.utils.signals import tone_burst
 
-if __name__ == '__main__':
+
+def main():
     # pathname for the input and output files
     pathname = gettempdir()
     phantom_data_path = 'phantom_data.mat'
@@ -110,7 +111,7 @@
 
         # set the input settings
         input_filename = f'example_input_{scan_line_index}.h5'
-        input_file_full_path = os.path.join(pathname, input_filename)
+        input_file_full_path = os.path.join(pathname, input_filename) # noqa: F841
         # set the input settings
         simulation_options = SimulationOptions(
             pml_inside=False,
@@ -160,3 +161,7 @@
 
     logging.log(logging.INFO, "Beamforming channel data and reconstructing the image...")
     beamform(channel_data)
+
+
+if __name__ == '__main__':
+    main()
diff --git a/examples/example_array_as_a_sensor.py b/examples/example_array_as_a_sensor.py
new file mode 100644
index 00000000..9cd934d2
--- /dev/null
+++ b/examples/example_array_as_a_sensor.py
@@ -0,0 +1,149 @@
+import matplotlib.pyplot as plt
+import numpy as np
+
+from kwave.data import Vector
+from kwave.kgrid import kWaveGrid
+from kwave.kmedium import kWaveMedium
+from kwave.ksensor import kSensor
+from kwave.ksource import kSource
+from kwave.kspaceFirstOrder2D import kspace_first_order_2d_gpu
+from kwave.options.simulation_execution_options import SimulationExecutionOptions
+from kwave.options.simulation_options import SimulationOptions
+from kwave.utils.conversion import cart2grid
+from kwave.utils.kwave_array import kWaveArray
+from kwave.utils.mapgen import make_cart_circle, make_disc
+from kwave.utils.signals import reorder_binary_sensor_data
+
+
+def main():
+    # create empty array
+    karray = kWaveArray()
+
+    # define arc properties
+    radius = 100e-3  # [m]
+    diameter = 8e-3  # [m]
+    ring_radius = 50e-3  # [m]
+    num_elements = 20
+
+    # orient all elements towards the center of the grid
+    focus_pos = Vector([0, 0])  # [m]
+
+    element_pos = make_cart_circle(ring_radius, num_elements, focus_pos)
+
+    for idx in range(num_elements):
+        karray.add_arc_element(element_pos[:, idx], radius, diameter, focus_pos)
+
+    # grid properties
+    N = Vector([256, 256])
+    d = Vector([0.5e-3, 0.5e-3])
+    kgrid = kWaveGrid(N, d)
+
+    # medium properties
+    medium = kWaveMedium(sound_speed=1500)
+
+    # time array
+    kgrid.makeTime(medium.sound_speed)
+
+    source = kSource()
+    x_offset = 20  # [pixels]
+    # make a small disc in the top left of the domain
+    source.p0 = make_disc(N, Vector([N.x / 4 + x_offset, N.y / 4]), 4)
+    source.p0[99:119, 59:199] = 1
+    logical_p0 = source.p0.astype(bool)
+    sensor = kSensor()
+    sensor.mask = element_pos
+    simulation_options = SimulationOptions(
+        save_to_disk=True,
+        data_cast='single',
+    )
+
+    execution_options = SimulationExecutionOptions(is_gpu_simulation=True)
+    output = kspace_first_order_2d_gpu(kgrid, source, sensor, medium, simulation_options, execution_options)
+        # TODO (walter): This should be done by kspaceFirstOrder
+    _, _, reorder_index = cart2grid(kgrid, element_pos)
+    sensor_data_point = reorder_binary_sensor_data(output['p'].T, reorder_index=reorder_index)
+    
+    # assign binary mask from karray to the source mask
+    sensor.mask = karray.get_array_binary_mask(kgrid)
+
+    output = kspace_first_order_2d_gpu(kgrid, source, sensor, medium, simulation_options, execution_options)
+    sensor_data = output['p'].T
+    combined_sensor_data = karray.combine_sensor_data(kgrid, sensor_data)
+
+    # =========================================================================
+    # VISUALIZATION
+    # =========================================================================
+
+    # create pml mask (re-use default size of 20 grid points from simulation_options)
+    pml_size = simulation_options.pml_x_size  # 20 [grid_points] 
+    pml_mask = np.zeros((N.x, N.y), dtype=bool)
+    pml_mask[:pml_size, :] = 1
+    pml_mask[:, :pml_size] = 1
+    pml_mask[-pml_size:, :] = 1
+    pml_mask[:, -pml_size:] = 1
+
+    # Plot source, sensor, and pml masks
+
+    # Assign unique values to each mask
+    sensor_val = sensor.mask * 1
+    logical_p0_val = logical_p0 * 2
+    pml_mask_val = pml_mask * 3
+
+    # Combine masks
+    combined_mask = sensor_val + logical_p0_val + pml_mask_val
+    combined_mask = np.flipud(combined_mask) 
+
+    # Define custom colormap
+    colors = [
+        (1, 1, 1),  # White (Background)
+        (233/255, 131/255, 0/255),  # Orange (Sensor) 
+        (254/255, 221/255, 92/255),  # Yellow (Sources) 
+        (0.8, 0.8, 0.8),  # Light Grey (PML Mask)
+    ]
+    cmap = plt.matplotlib.colors.ListedColormap(colors)
+
+    fig, ax = plt.subplots()
+    c = ax.pcolormesh(combined_mask, cmap=cmap, shading='auto')
+    plt.axis('image')
+
+    # Define labels for the colorbar
+    labels = {
+        0: 'None',
+        1: 'Sensor',
+        2: 'Initial pressure p0',
+        3: 'PML Mask',
+    }
+
+    colorbar = plt.colorbar(c, ticks=list(labels.keys()), ax=ax)
+    colorbar.ax.set_yticklabels(list(labels.values()))
+
+    ax.set_xticks([])
+    ax.set_yticks([])
+
+    # Plot recorded sensor data
+    fig, (ax1, ax2)= plt.subplots(ncols=1, nrows=2)
+    ax1.imshow(sensor_data_point, aspect='auto')
+    ax1.set_xlabel(r'Time [$\mu$s]')
+    ax1.set_ylabel('Detector Number')
+    ax1.set_title('Cartesian point detectors')
+    
+    ax2.imshow(combined_sensor_data, aspect='auto')
+    ax2.set_xlabel(r'Time [$\mu$s]')
+    ax2.set_ylabel('Detector Number')
+    ax2.set_title('Arc detectors')
+
+    fig.subplots_adjust(hspace=0.5)
+
+    # Plot a trace from the recorded sensor data
+    fig = plt.figure()
+    plt.plot(kgrid.t_array.squeeze() * 1e6, sensor_data_point[0, :], label='Cartesian point detectors')
+    plt.plot(kgrid.t_array.squeeze() * 1e6, combined_sensor_data[0, :], label='Arc detectors')
+    plt.xlabel(r'Time [$\mu$s]')
+    plt.ylabel('Pressure [pa]')
+    plt.legend()
+
+    plt.show()
+
+
+if __name__ == '__main__':
+    main()
\ No newline at end of file
diff --git a/examples/example_array_as_a_source.py b/examples/example_array_as_a_source.py
new file mode 100644
index 00000000..89b75227
--- /dev/null
+++ b/examples/example_array_as_a_source.py
@@ -0,0 +1,166 @@
+import matplotlib.pyplot as plt
+import numpy as np
+from matplotlib import colors
+from matplotlib.animation import FuncAnimation
+
+from kwave.data import Vector
+from kwave.kgrid import kWaveGrid
+from kwave.kmedium import kWaveMedium
+from kwave.ksensor import kSensor
+from kwave.ksource import kSource
+from kwave.kspaceFirstOrder2D import kspace_first_order_2d_gpu
+from kwave.options.simulation_execution_options import SimulationExecutionOptions
+from kwave.options.simulation_options import SimulationOptions
+from kwave.utils.kwave_array import kWaveArray
+from kwave.utils.signals import tone_burst
+
+
+def main():
+    # =========================================================================
+    # DEFINE KWAVEARRAY
+    # =========================================================================
+
+    # create empty array
+    karray = kWaveArray()
+
+    # define arc properties
+    radius = 50e-3  # [m]
+    diameter = 30e-3  # [m]
+    focus_pos = [-20e-3, 0]  # [m]
+
+    # add arc-shaped element
+    elem_pos = [10e-3, -40e-3]  # [m]
+    karray.add_arc_element(elem_pos, radius, diameter, focus_pos)
+
+    # add arc-shaped element
+    elem_pos = [20e-3, 0]  # [m]
+    karray.add_arc_element(elem_pos, radius, diameter, focus_pos)
+
+    # add arc-shaped element
+    elem_pos = [10e-3, 40e-3]  # [m]
+    karray.add_arc_element(elem_pos, radius, diameter, focus_pos)
+
+    # move the array down 10 mm, and rotate by 10 degrees (this moves all the
+    # elements together)
+    karray.set_array_position([10e-3, 0], 10)
+
+    # =========================================================================
+    # DEFINE GRID PROPERTIES
+    # =========================================================================
+
+    # grid properties
+    Nx = 256
+    dx = 0.5e-3
+    Ny = 256
+    dy = 0.5e-3
+    kgrid = kWaveGrid(Vector([Nx, Ny]), Vector([dx, dy]))
+
+    # medium properties
+    medium = kWaveMedium(sound_speed=1500)
+
+    # time array
+    kgrid.makeTime(medium.sound_speed)
+
+    # =========================================================================
+    # SIMULATION
+    # =========================================================================
+
+    # assign binary mask from karray to the source mask
+    source_p_mask = karray.get_array_binary_mask(kgrid)
+
+    # set source signals, one for each physical array element
+    f1 = 100e3
+    f2 = 200e3
+    f3 = 500e3
+    sig1 = tone_burst(1 / kgrid.dt, f1, 3).squeeze()
+    sig2 = tone_burst(1 / kgrid.dt, f2, 5).squeeze()
+    sig3 = tone_burst(1 / kgrid.dt, f3, 5).squeeze()
+
+    # combine source signals into one array
+    source_signal = np.zeros((3, max(len(sig1), len(sig2))))
+    source_signal[0, :len(sig1)] = sig1
+    source_signal[1, :len(sig2)] = sig2
+    source_signal[2, :len(sig3)] = sig3
+
+    # get distributed source signals (this automatically returns a weighted
+    # source signal for each grid point that forms part of the source)
+    source_p = karray.get_distributed_source_signal(kgrid, source_signal)
+
+    simulation_options = SimulationOptions(
+        save_to_disk=True,
+        data_cast='single',
+    )
+    execution_options = SimulationExecutionOptions(is_gpu_simulation=True)
+    # run k-Wave simulation (no sensor is used for this example)
+    # TODO: I would say proper behaviour would be to return the entire pressure field if sensor is None
+    sensor = kSensor()
+    sensor.mask = np.ones((Nx, Ny), dtype=bool)
+
+    source = kSource()
+    source.p_mask = source_p_mask
+    source.p = source_p
+
+    p = kspace_first_order_2d_gpu(kgrid, source, sensor, medium, simulation_options, execution_options)
+
+    p_field = np.reshape(p['p'], (kgrid.Nt, Nx, Ny))
+    p_field = np.transpose(p_field, (0, 2, 1))
+    # =========================================================================
+    # VISUALIZATION
+    # =========================================================================
+
+    # Normalize frames based on the maximum value over all frames
+    max_value = np.max(p_field)
+    normalized_frames = p_field / max_value
+
+    # Create a custom colormap (replace 'viridis' with your preferred colormap)
+    cmap = plt.get_cmap('viridis')
+
+    # Create a figure and axis
+    fig, ax = plt.subplots()
+
+    # Create an empty image with the first normalized frame
+    image = ax.imshow(normalized_frames[0], cmap=cmap, norm=colors.Normalize(vmin=0, vmax=1))
+
+    # Function to update the image for each frame
+    def update(frame):
+        image.set_data(normalized_frames[frame])
+        ax.set_title(f'Frame {frame + 1}/{kgrid.Nt}')
+        return [image]
+
+    # Create the animation
+    ani = FuncAnimation(fig, update, frames=kgrid.Nt, interval=100)  # Adjust interval as needed (in milliseconds)
+
+    # Save the animation as a video file (e.g., MP4)
+    video_filename = 'output_video1.mp4'
+    ani.save('/tmp/' + video_filename, writer='ffmpeg', fps=30)  # Adjust FPS as needed
+
+    # Show the animation (optional)
+    plt.show()
+
+    # Show the animation (optional)
+    plt.show()
+
+    # create pml mask (default size in 2D is 20 grid points)
+    pml_size = 20
+    pml_mask = np.zeros((Nx, Ny), dtype=bool)
+    pml_mask[:pml_size, :] = 1
+    pml_mask[:, :pml_size] = 1
+    pml_mask[-pml_size:, :] = 1
+    pml_mask[:, -pml_size:] = 1
+
+    # plot source and pml masks
+    plt.figure()
+    plt.imshow(np.logical_not(np.squeeze(source.p_mask | pml_mask)), aspect='auto', cmap='gray')
+    plt.xlabel('x-position [m]')
+    plt.ylabel('y-position [m]')
+    plt.title('Source and PML Masks')
+    plt.show()
+
+    # overlay the physical source positions
+    plt.figure()
+    # TODO: missing karray.plot_array(show=True)
+    # karray.plot_array(show=True)
+
+
+if __name__ == "__main__":
+    main()
diff --git a/examples/example_at_linear_array_transducer.py b/examples/example_at_linear_array_transducer.py
index 20c06fbb..75e8129d 100644
--- a/examples/example_at_linear_array_transducer.py
+++ b/examples/example_at_linear_array_transducer.py
@@ -2,7 +2,9 @@
 import numpy as np
 
 import kwave.data
-from kwave import kWaveGrid, SimulationOptions, kWaveMedium
+from kwave.kWaveSimulation import SimulationOptions
+from kwave.kgrid import kWaveGrid
+from kwave.kmedium import kWaveMedium
 from kwave.ksensor import kSensor
 from kwave.ksource import kSource
 from kwave.kspaceFirstOrder3D import kspaceFirstOrder3DC
@@ -11,87 +13,91 @@
 from kwave.utils.plot import voxel_plot
 from kwave.utils.signals import tone_burst
 
-# DEFINE LITERALS
-model = 1
-c0 = 1500
-rho0 = 1000
-source_f0 = 1e6
-source_amp = 1e6
-source_cycles = 5
-source_focus = 20e-3
-element_num = 15
-element_width = 1e-3
-element_length = 10e-3
-element_pitch = 2e-3
-translation = kwave.data.Vector([5e-3, 0, 8e-3])
-rotation = kwave.data.Vector([0, 20, 0])
-grid_size_x = 40e-3
-grid_size_y = 20e-3
-grid_size_z = 40e-3
-ppw = 3
-t_end = 35e-6
-cfl = 0.5
-
-# GRID
-dx = c0 / (ppw * source_f0)
-Nx = round(grid_size_x / dx)
-Ny = round(grid_size_y / dx)
-Nz = round(grid_size_z / dx)
-kgrid = kWaveGrid([Nx, Ny, Nz], [dx, dx, dx])
-kgrid.makeTime(c0, cfl, t_end)
-
-# SOURCE
-if element_num % 2 != 0:
-    ids = np.arange(1, element_num + 1) - np.ceil(element_num / 2)
-else:
-    ids = np.arange(1, element_num + 1) - (element_num + 1) / 2
-
-time_delays = -(np.sqrt((ids * element_pitch) ** 2 + source_focus ** 2) - source_focus) / c0
-time_delays = time_delays - min(time_delays)
-
-source_sig = source_amp * tone_burst(1 / kgrid.dt, source_f0, source_cycles,
-                                     signal_offset=np.round(time_delays / kgrid.dt).astype(int))
-karray = kWaveArray(bli_tolerance=0.05, upsampling_rate=10)
-
-for ind in range(1, element_num + 1):
-    x_pos = 0 - (element_num * element_pitch / 2 - element_pitch / 2) + (ind - 1) * element_pitch
-    karray.add_rect_element([x_pos, 0, kgrid.z_vec[0]], element_width, element_length, rotation)
-
-karray.set_array_position(translation, rotation)
-source = kSource()
-source.p_mask = karray.get_array_binary_mask(kgrid)
-voxel_plot(np.single(source.p_mask))
-source.p = karray.get_distributed_source_signal(kgrid, source_sig)
-# MEDIUM
-
-medium = kWaveMedium(sound_speed=c0, density=rho0)
-
-# SENSOR
-sensor_mask = np.zeros((Nx, Ny, Nz))
-sensor_mask[:, Ny // 2, :] = 1
-sensor = kSensor(sensor_mask, record=['p_max'])
-
-# SIMULATION
-simulation_options = SimulationOptions(
-    pml_auto=True,
-    pml_inside=False,
-    save_to_disk=True,
-    data_cast='single',
-)
-
-execution_options = SimulationExecutionOptions(is_gpu_simulation=True)
-
-sensor_data = kspaceFirstOrder3DC(kgrid=kgrid, medium=medium, source=source, sensor=sensor,
-                                  simulation_options=simulation_options, execution_options=execution_options)
-
-p_max = np.reshape(sensor_data['p_max'], (Nx, Nz), order='F')
-
-# VISUALISATION
-plt.figure()
-plt.imshow(1e-6 * p_max, extent=[1e3 * kgrid.x_vec[0][0], 1e3 * kgrid.x_vec[-1][0], 1e3 * kgrid.z_vec[0][0],
-                                 1e3 * kgrid.z_vec[-1][0]], aspect='auto')
-plt.xlabel('z-position [mm]')
-plt.ylabel('x-position [mm]')
-plt.title('Pressure Field')
-plt.colorbar(label='[MPa]')
-plt.show()
+
+def main():
+    c0 = 1500
+    rho0 = 1000
+    source_f0 = 1e6
+    source_amp = 1e6
+    source_cycles = 5
+    source_focus = 20e-3
+    element_num = 15
+    element_width = 1e-3
+    element_length = 10e-3
+    element_pitch = 2e-3
+    translation = kwave.data.Vector([5e-3, 0, 8e-3])
+    rotation = kwave.data.Vector([0, 20, 0])
+    grid_size_x = 40e-3
+    grid_size_y = 20e-3
+    grid_size_z = 40e-3
+    ppw = 3
+    t_end = 35e-6
+    cfl = 0.5
+
+    # GRID
+    dx = c0 / (ppw * source_f0)
+    Nx = round(grid_size_x / dx)
+    Ny = round(grid_size_y / dx)
+    Nz = round(grid_size_z / dx)
+    kgrid = kWaveGrid([Nx, Ny, Nz], [dx, dx, dx])
+    kgrid.makeTime(c0, cfl, t_end)
+
+    # SOURCE
+    if element_num % 2 != 0:
+        ids = np.arange(1, element_num + 1) - np.ceil(element_num / 2)
+    else:
+        ids = np.arange(1, element_num + 1) - (element_num + 1) / 2
+
+    time_delays = -(np.sqrt((ids * element_pitch) ** 2 + source_focus ** 2) - source_focus) / c0
+    time_delays = time_delays - min(time_delays)
+
+    source_sig = source_amp * tone_burst(1 / kgrid.dt, source_f0, source_cycles,
+                                         signal_offset=np.round(time_delays / kgrid.dt).astype(int))
+    karray = kWaveArray(bli_tolerance=0.05, upsampling_rate=10)
+
+    for ind in range(element_num):
+        x_pos = 0 - (element_num * element_pitch / 2 - element_pitch / 2) + ind * element_pitch
+        karray.add_rect_element([x_pos, 0, kgrid.z_vec[0]], element_width, element_length, rotation)
+
+    karray.set_array_position(translation, rotation)
+    source = kSource()
+    source.p_mask = karray.get_array_binary_mask(kgrid)
+    voxel_plot(np.single(source.p_mask))
+    source.p = karray.get_distributed_source_signal(kgrid, source_sig)
+    # MEDIUM
+
+    medium = kWaveMedium(sound_speed=c0, density=rho0)
+
+    # SENSOR
+    sensor_mask = np.zeros((Nx, Ny, Nz))
+    sensor_mask[:, Ny // 2, :] = 1
+    sensor = kSensor(sensor_mask, record=['p_max'])
+
+    # SIMULATION
+    simulation_options = SimulationOptions(
+        pml_auto=True,
+        pml_inside=False,
+        save_to_disk=True,
+        data_cast='single',
+    )
+
+    execution_options = SimulationExecutionOptions(is_gpu_simulation=True)
+
+    sensor_data = kspaceFirstOrder3DC(kgrid=kgrid, medium=medium, source=source, sensor=sensor,
+                                      simulation_options=simulation_options, execution_options=execution_options)
+
+    p_max = np.reshape(sensor_data['p_max'], (Nx, Nz), order='F')
+
+    # VISUALISATION
+    plt.figure()
+    plt.imshow(1e-6 * p_max, extent=[1e3 * kgrid.x_vec[0][0], 1e3 * kgrid.x_vec[-1][0], 1e3 * kgrid.z_vec[0][0],
+                                     1e3 * kgrid.z_vec[-1][0]], aspect='auto')
+    plt.xlabel('z-position [mm]')
+    plt.ylabel('x-position [mm]')
+    plt.title('Pressure Field')
+    plt.colorbar(label='[MPa]')
+    plt.show()
+
+
+if __name__ == "__main__":
+    main()
diff --git a/kwave/__init__.py b/kwave/__init__.py
index 0ba35dc4..26b545b8 100644
--- a/kwave/__init__.py
+++ b/kwave/__init__.py
@@ -5,8 +5,9 @@
 from os import environ
 from pathlib import Path
 
-
-VERSION = '0.2.1'
+# Test installation with:
+# python3 -m pip install -i https://test.pypi.org/simple/ --extra-index-url=https://pypi.org/simple/ k-Wave-python==0.3.0
+VERSION = '0.3.0'
 # Set environment variable to binaries to get rid of user warning
 # This code is a crutch and should be removed when kspaceFirstOrder
 # is refactored
diff --git a/kwave/executor.py b/kwave/executor.py
index d4691c5c..dc7ab265 100644
--- a/kwave/executor.py
+++ b/kwave/executor.py
@@ -41,7 +41,8 @@ def run_simulation(self, input_filename: str, output_filename: str, options: str
 
         return sensor_data
 
-    def parse_executable_output(self, output_filename: str) -> dotdict:
+    @staticmethod
+    def parse_executable_output(output_filename: str) -> dotdict:
 
         # Load the simulation and pml sizes from the output file
         # with h5py.File(output_filename, 'r') as output_file:
diff --git a/kwave/kWaveSimulation.py b/kwave/kWaveSimulation.py
index 14d6baec..b777bf85 100644
--- a/kwave/kWaveSimulation.py
+++ b/kwave/kWaveSimulation.py
@@ -67,24 +67,24 @@ def __init__(self,
             # set time reversal flag
             self.userarg_time_rev = False
 
-        #: Whether sensor.mask should be re-ordered.
-        #: True if sensor.mask is Cartesian with nearest neighbour interpolation which is calculated using a binary mask
-        #: and thus must be re-ordered
-        self.reorder_data              = False
+            #: Whether sensor.mask should be re-ordered.
+            #: True if sensor.mask is Cartesian with nearest neighbour interpolation which is calculated using a binary mask
+            #: and thus must be re-ordered
+            self.reorder_data = False
 
-        #: Whether the sensor.mask is binary
-        self.binary_sensor_mask        = True
+            #: Whether the sensor.mask is binary
+            self.binary_sensor_mask = True
 
-        # check if the sensor mask is defined as a list of cuboid corners
-        if self.sensor.mask is not None and self.sensor.mask.shape[0] == (2 * self.kgrid.dim):
-            self.userarg_cuboid_corners = True
-        else:
-            self.userarg_cuboid_corners = False
+            # check if the sensor mask is defined as a list of cuboid corners
+            if self.sensor.mask is not None and self.sensor.mask.shape[0] == (2 * self.kgrid.dim):
+                self.userarg_cuboid_corners = True
+            else:
+                self.userarg_cuboid_corners = False
 
-        #: If tse sensor is an object of the kWaveTransducer class
-        self.transducer_sensor         = False
+            #: If tse sensor is an object of the kWaveTransducer class
+            self.transducer_sensor = False
 
-        self.record = Recorder()
+            self.record = Recorder()
 
         # transducer source flags
         #: transducer is object of kWaveTransducer class
@@ -502,8 +502,8 @@ def input_checking(self, calling_func_name) -> None:
         self.c_ref, self.c_ref_compression, self.c_ref_shear \
             = set_sound_speed_ref(self.medium, opt.simulation_type)
 
-        self.check_source(k_dim)
-        self.check_sensor(k_dim, self.kgrid.Nt)
+        self.check_source(k_dim, self.kgrid.Nt)
+        self.check_sensor(k_dim)
         self.check_kgrid_time()
         self.precision = self.select_precision(opt)
         self.check_input_combinations(opt, user_medium_density_input, k_dim, pml_size, self.kgrid.N)
@@ -606,12 +606,12 @@ def check_medium(medium, kgrid_k, simulation_type: SimulationType) -> bool:
             medium.check_fields(kgrid_k.shape)
         return user_medium_density_input
 
-    def check_source(self, kgrid_dim) -> None:
+    def check_sensor(self, kgrid_dim)-> None:
         """
-        Check the source properties for correctness and validity
+        Check the Sensor properties for correctness and validity
 
         Args:
-            kgrid_dim: kWaveGrid dimension
+            k_dim: kWaveGrid dimensionality
 
         Returns:
             None
@@ -631,17 +631,6 @@ def check_source(self, kgrid_dim) -> None:
             if not isinstance(self.sensor, NotATransducer):
                 if kgrid_dim == 2:
 
-                    # check field names, including the directivity inputs for the
-                    # regular 2D code, but not the axisymmetric code
-                    # TODO question to Walter: we do not need following checks anymore
-                    #  because we have kSensor that defines the structure, right?
-                    # if self.axisymmetric:
-                    #     check_field_names(self.sensor, *['mask', 'time_reversal_boundary_data', 'frequency_response',
-                    #                                        'record_mode', 'record', 'record_start_index'])
-                    # else:
-                    #     check_field_names(self.sensor, *['mask', 'directivity', 'time_reversal_boundary_data',
-                    #                                      'frequency_response', 'record_mode', 'record', 'record_start_index'])
-
                     # check for sensor directivity input and set flag
                     directivity = self.sensor.directivity
                     if directivity is not None and self.sensor.directivity.angle is not None:
@@ -663,15 +652,7 @@ def check_source(self, kgrid_dim) -> None:
                         directivity.set_unique_angles(self.sensor.mask)
                         directivity.set_wavenumbers(self.kgrid)
 
-                else:
-                    # TODO question to Walter: we do not need following checks anymore because
-                    #  we have kSensor that defines the structure, right?
-                    # check field names without directivity inputs (these are not supported in 1 or 3D)
-                    # check_field_names(self.sensor, *['mask', 'time_reversal_boundary_data', 'frequency_response',
-                    #                                    'record_mode', 'record', 'record_start_index'])
-                    pass
-
-                # check for time reversal inputs and set flgs
+                # check for time reversal inputs and set flags
                 if not self.options.simulation_type.is_elastic_simulation() and \
                         self.sensor.time_reversal_boundary_data is not None:
                     self.record.p = False
@@ -846,12 +827,12 @@ def check_source(self, kgrid_dim) -> None:
         if kgrid_dim == 2 and self.use_sensor and self.compute_directivity and self.time_rev:
             logging.log(logging.WARN, 'sensor directivity fields are not used for time reversal.')
 
-    def check_sensor(self, k_dim, k_Nt) -> None:
+    def check_source(self, k_dim, k_Nt) -> None:
         """
-        Check the Sensor properties for correctness and validity
+        Check the source properties for correctness and validity
 
         Args:
-            k_dim: kWaveGrid dimensionality
+            kgrid_dim: kWaveGrid dimension
             k_Nt: Number of time steps in kWaveGrid
 
         Returns:
diff --git a/kwave/kWaveSimulation_helper/create_storage_variables.py b/kwave/kWaveSimulation_helper/create_storage_variables.py
index 7cf8f844..c635d024 100644
--- a/kwave/kWaveSimulation_helper/create_storage_variables.py
+++ b/kwave/kWaveSimulation_helper/create_storage_variables.py
@@ -8,7 +8,7 @@
 from kwave.utils.dotdictionary import dotdict
 
 
-# Note from Farid: This function/file is very suspicios. I'm pretty sure that the implementation is not correct.
+# Note from Farid: This function/file is very suspicious. I'm pretty sure that the implementation is not correct.
 # Full test-coverage is required for bug-fixes!
 
 def create_storage_variables(
diff --git a/kwave/kWaveSimulation_helper/display_simulation_params.py b/kwave/kWaveSimulation_helper/display_simulation_params.py
index a35212ec..e0d50a54 100644
--- a/kwave/kWaveSimulation_helper/display_simulation_params.py
+++ b/kwave/kWaveSimulation_helper/display_simulation_params.py
@@ -75,7 +75,7 @@ def print_max_supported_freq(kgrid, c_min):
             logging.log(logging.INFO, '  maximum supported frequency: ', scale_SI(k_max_all * c_min / (2*np.pi))[0], 'Hz')
         else:
             logging.log(logging.INFO, '  maximum supported frequency: ', scale_SI(k_max.x * c_min / (2*np.pi))[0],
-                  'Hz by ', scale_SI(kgrid.ky_max * c_min / (2*np.pi))[0], 'Hz')
+                  'Hz by ', scale_SI(k_max.y * c_min / (2*np.pi))[0], 'Hz')
 
     elif kgrid.dim == 3:
         # display maximum supported frequency
diff --git a/kwave/kWaveSimulation_helper/save_to_disk_func.py b/kwave/kWaveSimulation_helper/save_to_disk_func.py
index 2fbd25e1..f4771034 100644
--- a/kwave/kWaveSimulation_helper/save_to_disk_func.py
+++ b/kwave/kWaveSimulation_helper/save_to_disk_func.py
@@ -16,7 +16,8 @@
 
 def save_to_disk_func(
         kgrid: kWaveGrid, medium: kWaveMedium, source,
-        opt: SimulationOptions, values: dotdict, flags: dotdict):
+        opt: SimulationOptions, auto_chunk: bool,
+        values: dotdict, flags: dotdict):
     # update command line status
     logging.log(logging.INFO, '  precomputation completed in ', scale_time(TicToc.toc()))
     TicToc.tic()
@@ -56,7 +57,8 @@ def save_to_disk_func(
     # =========================================================================
 
     remove_z_dimension(float_variables, kgrid.dim)
-    save_file(opt.input_filename, integer_variables, float_variables, opt.hdf_compression_level)
+    save_file(opt.input_filename, integer_variables, float_variables, opt.hdf_compression_level, 
+              auto_chunk=auto_chunk)
 
     # update command line status
     logging.log(logging.INFO, '  completed in ', scale_time(TicToc.toc()))
@@ -445,12 +447,12 @@ def enforce_filename_standards(filepath):
     return filepath, filename_ext
 
 
-def save_file(filepath, integer_variables, float_variables, hdf_compression_level):
+def save_file(filepath, integer_variables, float_variables, hdf_compression_level, auto_chunk):
     filepath, filename_ext = enforce_filename_standards(filepath)
 
     # save file
     if filename_ext == '.h5':
-        save_h5_file(filepath, integer_variables, float_variables, hdf_compression_level)
+        save_h5_file(filepath, integer_variables, float_variables, hdf_compression_level, auto_chunk)
 
     elif filename_ext == '.mat':
         save_mat_file(filepath, integer_variables, float_variables)
@@ -459,7 +461,7 @@ def save_file(filepath, integer_variables, float_variables, hdf_compression_leve
         raise NotImplementedError('unknown file extension for ''save_to_disk'' filename')
 
 
-def save_h5_file(filepath, integer_variables, float_variables, hdf_compression_level):
+def save_h5_file(filepath, integer_variables, float_variables, hdf_compression_level, auto_chunk):
     # ----------------
     # SAVE HDF5 FILE
     # ----------------
@@ -474,7 +476,7 @@ def save_h5_file(filepath, integer_variables, float_variables, hdf_compression_l
     for key, value in float_variables.items():
         # cast matrix to single precision
         value = np.array(value, dtype=np.float32)
-        write_matrix(filepath, value, key, hdf_compression_level)
+        write_matrix(filepath, value, key, hdf_compression_level, auto_chunk)
         del value
 
     # change all the index variables to be in 64-bit unsigned integers
@@ -482,7 +484,7 @@ def save_h5_file(filepath, integer_variables, float_variables, hdf_compression_l
     for key, value in integer_variables.items():
         # cast matrix to 64-bit unsigned integer
         value = np.array(value, dtype=np.uint64)
-        write_matrix(filepath, value, key, hdf_compression_level)
+        write_matrix(filepath, value, key, hdf_compression_level, auto_chunk)
         del value
 
     # set additional file attributes
diff --git a/kwave/kWaveSimulation_helper/scale_source_terms_func.py b/kwave/kWaveSimulation_helper/scale_source_terms_func.py
index 6423673b..d5bdf84a 100644
--- a/kwave/kWaveSimulation_helper/scale_source_terms_func.py
+++ b/kwave/kWaveSimulation_helper/scale_source_terms_func.py
@@ -6,7 +6,6 @@
 from kwave.kgrid import kWaveGrid
 from kwave.ksource import kSource
 from kwave.utils.dotdictionary import dotdict
-from kwave.utils.matlab import matlab_mask
 
 
 def scale_source_terms_func(
@@ -128,8 +127,11 @@ def scale_pressure_source_dirichlet(source_p, c0, N, p_source_pos_index):
     else:
         # compute the scale parameter seperately for each source
         # position based on the sound speed at that position
-        for p_index in range(source_p.size[0]):
-            source_p[p_index, :] = source_p[p_index, :] / (N * (c0[p_source_pos_index[p_index]] ** 2))
+        ind = range(source_p[:, 0].size)
+        mask = p_source_pos_index.flatten('F')[ind]
+        scale = 1.0 / (N * np.expand_dims(c0.ravel(order='F')[mask.ravel(order='F')] ** 2, axis=-1) )
+        source_p[ind, :] *= scale
+
     return source_p
 
 
@@ -182,9 +184,10 @@ def scale_pressure_source_uniform_grid(source_p, c0, N, dx, dt, p_source_pos_ind
     else:
         # compute the scale parameter seperately for each source
         # position based on the sound speed at that position
-        for p_index in range(source_p[:, 0].size):
-            source_p[p_index, :] = source_p[p_index, :] * \
-                                   (2 * dt / (N * matlab_mask(c0, p_source_pos_index.flatten('F')[p_index]) * dx))
+        ind = range(source_p[:, 0].size)
+        mask = p_source_pos_index.flatten('F')[ind]
+        scale = (2.0 * dt) / (N * np.expand_dims(c0.ravel(order='F')[mask.ravel(order='F')], axis=-1) * dx)
+        source_p[ind, :] *= scale
     return source_p
 
 
@@ -229,8 +232,8 @@ def scale_stress_source(source, c0, is_source_exists, is_p0_exists, source_val,
 
                 # compute the scale parameter seperately for each source
                 # position based on the sound speed at that position
-                for s_index in range(source_val.size[0]):
-                    source_val[s_index, :] = source_val[s_index, :] * (2 * dt * c0[s_source_pos_index[s_index]] / (N * dx))
+                s_index = range(source_val.size[0])
+                source_val[s_index, :] = source_val[s_index, :] * (2 * dt * c0[s_source_pos_index[s_index]] / (N * dx))
     return source_val
 
 
@@ -322,8 +325,8 @@ def scale_velocity_source(is_source, source_u_mode, source_val, c0, dt, u_source
     else:
         # compute the scale parameter seperately for each source position
         # based on the sound speed at that position
-        for u_index in range(source_val.size[0]):
-            source_val[u_index, :] = source_val[u_index, :] * (2 * c0(u_source_pos_index[u_index]) * dt / d_direction)
+        u_index = range(source_val.size[0])
+        source_val[u_index, :] = source_val[u_index, :] * (2 * c0[u_source_pos_index[u_index]] * dt / d_direction)
     return source_val
 
 
diff --git a/kwave/kgrid.py b/kwave/kgrid.py
index 9da1d550..de318769 100644
--- a/kwave/kgrid.py
+++ b/kwave/kgrid.py
@@ -93,10 +93,15 @@ def t_array(self):
         """
         time array [s]
         """
+        # TODO (walter): I would change this functionality to return a time array even if Nt or dt are not yet set
+        #  (e.g. if they are still 'auto')
+
         if self.Nt == 'auto' or self.dt == 'auto':
             return 'auto'
         else:
             t_array = np.arange(0, self.Nt) * self.dt
+            # TODO: adding this extra dimension seems unnecessary
+            # This leads to an extra squeeze when plotting e.g. in example "array as sensor" on lines 110 and 111
             return np.expand_dims(t_array, axis=0)
 
     @t_array.setter
@@ -395,6 +400,7 @@ def k_max_all(self):
     # functions that can only be accessed by class members
     ########################################
     @staticmethod
+    # TODO (walter): convert this name to snake case
     def makeDim(num_points, spacing):
         """
         Create the grid parameters for a single spatial direction
@@ -451,6 +457,7 @@ def highest_prime_factors(self, axisymmetric=None) -> np.ndarray:
                           largest_prime_factor(self.Nz)]
         return np.array(prime_facs)
 
+    # TODO (walter): convert this name to snake case
     def makeTime(self, c, cfl=CFL_DEFAULT, t_end=None):
         """
         Compute Nt and dt based on the cfl number and grid size, where
@@ -538,6 +545,7 @@ def kz_vec_dtt(self, dtt_type):
         return kz_vec_dtt, M
 
     @staticmethod
+    # TODO (walter): convert this name to snake case
     def makeDTTDim(Nx, dx, dtt_type):
         """
         Create the DTT grid parameters for a single spatial direction
@@ -597,6 +605,7 @@ def makeDTTDim(Nx, dx, dtt_type):
     ########################################
     # functions for non-uniform grids
     ########################################
+    # TODO (walter): convert this name to snake case
     def setNUGrid(self, dim, n_vec, dudn, n_vec_sg, dudn_sg):
         """
         Function to set non-uniform grid parameters in specified dimension
diff --git a/kwave/kspaceFirstOrder2D.py b/kwave/kspaceFirstOrder2D.py
index c77d35f5..3b151bc5 100644
--- a/kwave/kspaceFirstOrder2D.py
+++ b/kwave/kspaceFirstOrder2D.py
@@ -55,6 +55,12 @@ def kspace_first_order_2d_gpu(
     GPU binary.
     """
     execution_options.is_gpu_simulation = True  # force to GPU
+    assert isinstance(kgrid, kWaveGrid), 'kgrid must be a kWaveGrid object'
+    assert isinstance(medium, kWaveMedium), 'medium must be a kWaveMedium object'
+    assert isinstance(simulation_options, SimulationOptions), 'simulation_options must be a SimulationOptions object'
+    assert isinstance(execution_options,
+                      SimulationExecutionOptions), 'execution_options must be a SimulationExecutionOptions object'
+
     sensor_data = kspaceFirstOrder2DC(
         kgrid=kgrid,
         source=source,
@@ -113,7 +119,7 @@ def kspaceFirstOrder2DC(
     Args:
         kgrid: kWaveGrid instance
         source: kWaveSource instance
-        sensor: NotATransducer or kSensor instance
+        sensor: NotATransducer or kSensor instance or None
         medium: kWaveMedium instance
         simulation_options: SimulationOptions instance
         execution_options: SimulationExecutionOptions instance
@@ -137,7 +143,7 @@ def kspaceFirstOrder2DC(
 def kspaceFirstOrder2D(
         kgrid: kWaveGrid,
         source: kSource,
-        sensor: Union[NotATransducer, kSensor],
+        sensor: Union[NotATransducer, kSensor, None],
         medium: kWaveMedium,
         simulation_options: SimulationOptions,
         execution_options: SimulationExecutionOptions
@@ -281,7 +287,7 @@ def kspaceFirstOrder2D(
         kgrid: kWaveGrid instance
         medium: kWaveMedium instance
         source: kWaveSource instance
-        sensor: kWaveSensor instance
+        sensor: kWaveSensor instance or None
 
     Returns:
 
@@ -385,7 +391,7 @@ def kspaceFirstOrder2D(
         retract_size = [[options.pml_x_size, options.pml_y_size, options.pml_z_size]]
 
         # run subscript to save files to disk
-        save_to_disk_func(k_sim.kgrid, k_sim.medium, k_sim.source, k_sim.options,
+        save_to_disk_func(k_sim.kgrid, k_sim.medium, k_sim.source, k_sim.options, execution_options.auto_chunking,
                           dotdict({
                               'ddx_k_shift_pos': k_sim.ddx_k_shift_pos,
                               'ddx_k_shift_neg': k_sim.ddx_k_shift_neg,
diff --git a/kwave/kspaceFirstOrder3D.py b/kwave/kspaceFirstOrder3D.py
index b4a51a14..59144603 100644
--- a/kwave/kspaceFirstOrder3D.py
+++ b/kwave/kspaceFirstOrder3D.py
@@ -413,7 +413,7 @@ def kspaceFirstOrder3D(
         retract_size = [[options.pml_x_size, options.pml_y_size, options.pml_z_size]]
 
         # run subscript to save files to disk
-        save_to_disk_func(k_sim.kgrid, k_sim.medium, k_sim.source, k_sim.options,
+        save_to_disk_func(k_sim.kgrid, k_sim.medium, k_sim.source, k_sim.options, execution_options.auto_chunking,
                           dotdict({
                               'ddx_k_shift_pos': k_sim.ddx_k_shift_pos,
                               'ddx_k_shift_neg': k_sim.ddx_k_shift_neg,
diff --git a/kwave/kspaceFirstOrderAS.py b/kwave/kspaceFirstOrderAS.py
index b3f51d7b..3258b427 100644
--- a/kwave/kspaceFirstOrderAS.py
+++ b/kwave/kspaceFirstOrderAS.py
@@ -314,7 +314,7 @@ def kspaceFirstOrderAS(
         retract_size = [[options.pml_x_size, options.pml_y_size, options.pml_z_size]]
 
         # run subscript to save files to disk
-        save_to_disk_func(k_sim.kgrid, k_sim.medium, k_sim.source, k_sim.options,
+        save_to_disk_func(k_sim.kgrid, k_sim.medium, k_sim.source, k_sim.options, execution_options.auto_chunking,
                           dotdict({
                               'ddx_k_shift_pos': k_sim.ddx_k_shift_pos,
                               'ddx_k_shift_neg': k_sim.ddx_k_shift_neg,
diff --git a/kwave/options/simulation_execution_options.py b/kwave/options/simulation_execution_options.py
index 4b3bc833..4736665a 100644
--- a/kwave/options/simulation_execution_options.py
+++ b/kwave/options/simulation_execution_options.py
@@ -34,6 +34,9 @@ class SimulationExecutionOptions:
     system_call: Optional[str] = None
     verbose_level: int = 0
 
+    # determine whether chunking is handled automatically (the default), or manually
+    auto_chunking: Optional[bool] = True
+
     # show simulation log
     show_sim_log: bool = True
 
diff --git a/kwave/utils/conversion.py b/kwave/utils/conversion.py
index 015af5bc..1c645c4c 100644
--- a/kwave/utils/conversion.py
+++ b/kwave/utils/conversion.py
@@ -375,14 +375,14 @@ def tol_star(tolerance, kgrid, point, debug):
     if tol is None or tolerance != tol or len(subs0) != kgrid_dim:
         tol = tolerance
 
-        decay_subs = np.ceil(1 / (np.pi * tol))
+        decay_subs = int(np.ceil(1 / (np.pi * tol)))
 
         lin_ind = np.arange(-decay_subs, decay_subs + 1)
 
         if kgrid_dim == 1:
             is0 = lin_ind
         elif kgrid_dim == 2:
-            is0, js0 = np.meshgrid(lin_ind, lin_ind)
+            is0, js0 = np.meshgrid(lin_ind, lin_ind, indexing='ij')
         elif kgrid_dim == 3:
             is0, js0, ks0 = np.meshgrid(lin_ind, lin_ind, lin_ind, indexing='ij')
 
@@ -390,8 +390,8 @@ def tol_star(tolerance, kgrid, point, debug):
             subs0 = [is0]
         elif kgrid_dim == 2:
             instar = np.abs(is0 * js0) <= decay_subs
-            is0 = is0[instar]
-            js0 = js0[instar]
+            is0 = matlab_mask(is0, instar)
+            js0 = matlab_mask(js0, instar)
             subs0 = [is0, js0]
         elif kgrid_dim == 3:
             instar = np.abs(is0 * js0 * ks0) <= decay_subs
@@ -404,6 +404,7 @@ def tol_star(tolerance, kgrid, point, debug):
     js = subs0[1].copy() if kgrid_dim > 1 else []
     ks = subs0[2].copy() if kgrid_dim > 2 else []
 
+    # returns python index value (0 start) not MATLAB index (1 start)
     x_closest, x_closest_ind = find_closest(kgrid.x_vec, point[0])
 
     if np.abs(x_closest - point[0]) < ongrid_threshold:
@@ -428,6 +429,7 @@ def tol_star(tolerance, kgrid, point, debug):
             js = js[ks == 0]
             ks = ks[ks == 0]
 
+    # TODO: closest index is added to is_, +1 might not be correct
     is_ += x_closest_ind + 1
     if kgrid_dim > 1:
         js += y_closest_ind + 1
@@ -436,6 +438,7 @@ def tol_star(tolerance, kgrid, point, debug):
 
     if kgrid_dim == 1:
         inbounds = (1 <= is_) & (is_ <= kgrid.Nx)
+        # TODO: this should likely be matlabmask and indexes should be checked.
         is_ = is_[inbounds]
     elif kgrid_dim == 2:
         inbounds = (1 <= is_) & (is_ <= kgrid.Nx) & (1 <= js) & (js <= kgrid.Ny)
@@ -454,7 +457,10 @@ def tol_star(tolerance, kgrid, point, debug):
     elif kgrid_dim == 3:
         lin_ind = kgrid.Nx * kgrid.Ny * (ks - 1) + kgrid.Nx * (js - 1) + is_
 
-    return lin_ind, is_ - 1, js - 1, ks - 1  # -1 for mapping from Matlab indexing to Python indexing
+    # TODO: in current form linear indexing matches MATLAB, but might break in 0 indexed python
+    return lin_ind, np.array(is_) - 1, np.array(js) - 1, np.array(
+        ks) - 1  # -1 for mapping from Matlab indexing to Python indexing
+
 
 
 def find_closest(array, value):
diff --git a/kwave/utils/filters.py b/kwave/utils/filters.py
index 1498bea5..5f7f900b 100644
--- a/kwave/utils/filters.py
+++ b/kwave/utils/filters.py
@@ -44,12 +44,12 @@ def single_sided_correction(func_fft: np.ndarray, fft_len: int, dim: int) -> np.
 
         # even FFT length
         if dim == 0:
-            func_fft[1: -1, :, :, :] = func_fft[1: -1, :, :, :] * 2
+            func_fft[1: -1] = func_fft[1: -1] * 2
+        elif dim == 1:
+            func_fft[:, 1: -1] = func_fft[:, 1: -1] * 2
         elif dim == 2:
-            func_fft[:, 1: -1, :, :] = func_fft[:, 1: -1, :, :] * 2
+            func_fft[:, :, 1: -1] = func_fft[:, :, 1: -1] * 2
         elif dim == 3:
-            func_fft[:, :, 1: -1, :] = func_fft[:, :, 1: -1, :] * 2
-        elif dim == 4:
             func_fft[:, :, :, 1: -1] = func_fft[:, :, :, 1: -1] * 2
 
     return func_fft
diff --git a/kwave/utils/io.py b/kwave/utils/io.py
index 8016043e..da9f8237 100644
--- a/kwave/utils/io.py
+++ b/kwave/utils/io.py
@@ -9,10 +9,10 @@
 import h5py
 import numpy as np
 
+import kwave
 from .conversion import cast_to_type
 from .data import get_date_string
 from .dotdictionary import dotdict
-import kwave
 
 
 def get_h5_literals():
@@ -52,10 +52,12 @@ def get_h5_literals():
     return literals
 
 
-def write_matrix(filename, matrix: np.ndarray, matrix_name, compression_level=None):
+def write_matrix(filename, matrix: np.ndarray, matrix_name: str, compression_level:int =None, auto_chunk: bool =True):
     # get literals
     h5_literals = get_h5_literals()
 
+    assert isinstance(auto_chunk, bool), "auto_chunk must be a boolean."
+
     if compression_level is None:
         compression_level = h5_literals.HDF_COMPRESSION_LEVEL
 
@@ -181,8 +183,9 @@ def write_matrix(filename, matrix: np.ndarray, matrix_name, compression_level=No
     # allocate a holder for the new matrix within the file
     opts = {
         'dtype': data_type_matlab,
-        'chunks': tuple(chunk_size)
+        'chunks': auto_chunk if auto_chunk is True else tuple(chunk_size)
     }
+    
     if compression_level != 0:
         # use compression
         opts['compression'] = compression_level
diff --git a/kwave/utils/kwave_array.py b/kwave/utils/kwave_array.py
index 20bc3151..f79731c0 100644
--- a/kwave/utils/kwave_array.py
+++ b/kwave/utils/kwave_array.py
@@ -2,7 +2,7 @@
 import time
 from dataclasses import dataclass
 from math import ceil
-from typing import Optional
+from typing import Iterable, Optional
 
 import numpy as np
 from numpy import arcsin, pi, cos, size, array
@@ -26,6 +26,7 @@ class Element:
     active: bool
     measure: float
 
+    label: Optional[str] = None
     group_type: Optional[str] = None
     element_number: Optional[int] = None
 
@@ -36,10 +37,13 @@ class Element:
     radius_of_curvature: Optional[float] = None
     position: Optional[np.ndarray] = None
     focus_position: Optional[np.ndarray] = None
+    
+    # custom element
+    integration_points: Optional[np.ndarray] = None
 
     length: Optional[float] = None
     width: Optional[float] = None
-    orientation: Optional = None
+    orientation: Optional[np.ndarray] = None
 
     start_point: Optional[np.ndarray] = None
     end_point: Optional[np.ndarray] = None
@@ -220,7 +224,43 @@ def add_bowl_element(self, position, radius, diameter, focus_pos):
             active=True,
             measure=area
         ))
+        
+    def add_custom_element(self, integration_points, measure, element_dim, label):
+        
+        assert isinstance(integration_points, (np.ndarray)), "'integration_points' must be a numpy array"
+        assert isinstance(measure, (int, float)), "'measure' must be an integer or float"
+        assert isinstance(element_dim, (int)) and element_dim in [1, 2, 3], "'element_dim' must be an integer and either 1, 2 or 3"
+        assert isinstance(label, (str)), "'label' must be a string"
+        
+        # check the dimensionality of the integration points
+        input_dim = integration_points.shape[0]
+        if (input_dim < 1) or (input_dim > 3):
+            raise ValueError("Input integration_points must be a 1 x N (in 1D), 2 x N (in 2D), or 3 x N (in 3D) array.")
+
+        # check if this is the first element, and set the dimension
+        if self.number_elements == 0:
+            self.dim = input_dim
+        
+        # check that the element is being added to an array with the
+        # correct dimensions
+        if self.dim != input_dim:
+            raise ValueError(f"{input_dim}D custom element cannot be added to an array with {self.dim}D elements.")
 
+        self.number_elements += 1
+        
+        self.elements.append(
+            Element(
+                group_id=0,
+                type='custom',
+                dim=element_dim,
+                label=label,
+                integration_points=integration_points,
+                active=True,
+                measure=measure
+            )
+        )
+        
+        
     def add_rect_element(self, position, Lx, Ly, theta):
 
         assert isinstance(position, (list, tuple)), "'position' must be a list or tuple"
@@ -263,10 +303,10 @@ def add_rect_element(self, position, Lx, Ly, theta):
 
     def add_arc_element(self, position, radius, diameter, focus_pos):
 
-        assert isinstance(position, (list, tuple)), "'position' must be list or tuple"
+        assert isinstance(position, Iterable), "'position' must be list, tuple or Vector"
         assert isinstance(radius, (int, float)), "'radius' must be an integer or float"
         assert isinstance(diameter, (int, float)), "'diameter' must be an integer or float"
-        assert isinstance(focus_pos, (list, tuple)), "'focus_pos' must be list or tuple"
+        assert isinstance(focus_pos, Iterable), "'focus_pos' must be list, tuple or Vector"
         assert len(position) == 2, "'position' must have 2 elements"
         assert len(focus_pos) == 2, "'focus_pos' must have 2 elements"
 
@@ -396,7 +436,7 @@ def get_array_binary_mask(self, kgrid):
 
         for ind in range(self.number_elements):
             grid_weights = self.get_off_grid_points(kgrid, ind, True)
-            mask |= grid_weights
+            mask = np.bitwise_or(np.squeeze(mask), grid_weights)
 
         return mask
 
@@ -602,12 +642,15 @@ def get_distributed_source_signal(self, kgrid, source_signal):
         return distributed_source_signal
 
     def combine_sensor_data(self, kgrid, sensor_data):
+
         self.check_for_elements()
 
         mask = self.get_array_binary_mask(kgrid)
         mask_ind = matlab_find(mask).squeeze(axis=-1)
 
         Nt = np.shape(sensor_data)[1]
+        # TODO (Walter): this assertion does not work when "auto" is set
+        # assert kgrid.Nt == Nt, 'sensor_data must have the same number of time steps as kgrid'
 
         combined_sensor_data = np.zeros((self.number_elements, Nt))
 
@@ -788,14 +831,14 @@ def off_grid_points(kgrid, points,
         else:
             # create an array of neighbouring grid points for BLI evaluation
             if kgrid.dim == 1:
-                ind, is_ = tol_star(bli_tolerance, kgrid, point, debug)
+                ind, is_, _, _ = tol_star(bli_tolerance, kgrid, point, debug)
                 xs = x_vec[is_]
                 xyz = xs
             elif kgrid.dim == 2:
-                ind, is_, js = tol_star(bli_tolerance, kgrid, point, debug)
+                ind, is_, js, _ = tol_star(bli_tolerance, kgrid, point, debug)
                 xs = x_vec[is_]
                 ys = y_vec[js]
-                xyz = np.array([xs, ys]).T
+                xyz = np.array([xs, ys]).squeeze().T
             elif kgrid.dim == 3:
                 ind, is_, js, ks = tol_star(bli_tolerance, kgrid, point, debug)
                 xs = x_vec[is_.astype(int)].squeeze(axis=-1)
@@ -820,16 +863,17 @@ def off_grid_points(kgrid, points,
                         mask_t = sinc(pi_on_dxyz * (xyz - point.T))
                     else:
                         mask_t = sinc(pi_on_dxyz * (xyz - point.T))
-                mask_t = np.prod(mask_t, axis=1)
+                current_mask_t = np.prod(np.atleast_2d(mask_t), axis=1)
 
                 # apply scaling for non-uniform grid
                 if kgrid.nonuniform:
-                    mask_t = mask_t * BLIscale
+                    current_mask_t = mask_t * BLIscale
 
+                updated_mask_value = matlab_mask(mask, ind - 1).squeeze(axis=-1) + scale[point_ind] * current_mask_t
                 # add this contribution to the overall source mask
                 mask = matlab_assign(mask,
                                      ind - 1,
-                                     matlab_mask(mask, ind - 1).squeeze(axis=-1) + scale[point_ind] * mask_t)
+                                     updated_mask_value)
 
         # update the waitbar
         if display_wait_bar and (point_ind % wait_bar_update_freq == 0):
diff --git a/kwave/utils/signals.py b/kwave/utils/signals.py
index 3c1c34fe..21660c42 100644
--- a/kwave/utils/signals.py
+++ b/kwave/utils/signals.py
@@ -10,7 +10,7 @@
 from .data import scale_SI
 from .mapgen import ndgrid
 from .math import sinc, gaussian
-from .matlab import matlab_mask, unflatten_matlab_mask
+from .matlab import matlab_mask, unflatten_matlab_mask, rem
 from .matrix import broadcast_axis, num_dim
 
 
@@ -343,8 +343,10 @@ def tone_burst(sample_freq, signal_freq, num_cycles, envelope='Gaussian', plot_s
 
     # create the tone burst
     tone_length = num_cycles / signal_freq  # [s]
-    # We want to include the endpoint but only if it's divisible by the step-size
-    if tone_length % dt < 1e-18:
+    # We want to include the endpoint but only if it's divisible by the step-size. 
+    # Modulo operator is not stable, so multiple conditions included.
+    # if ( (tone_length % dt) < 1e-18 or (np.abs(tone_length % dt - dt) < 1e-18) ):
+    if (rem(tone_length, dt) < 1e-18):
         tone_t = np.linspace(0, tone_length, int(tone_length / dt) + 1)
     else:
         tone_t = np.arange(0, tone_length, dt)
@@ -430,6 +432,48 @@ def tone_burst(sample_freq, signal_freq, num_cycles, envelope='Gaussian', plot_s
     return signal
 
 
+def reorder_sensor_data(kgrid, sensor, sensor_data: np.ndarray) -> np.ndarray:
+    """
+    Reorders the sensor data based on the coordinates of the sensor points.
+
+    Args:
+        kgrid: The k-Wave grid object.
+        sensor: The k-Wave sensor object.
+        sensor_data: The sensor data to be reordered.
+
+    Returns:
+        np.ndarray of the reordered sensor data.
+
+    Raises:
+        ValueError: If the simulation is not 2D or the sensor is not defined as a binary mask.
+    """
+    # check simulation is 2D
+    if kgrid.dim != 2:
+        raise ValueError('The simulation must be 2D.')
+
+    # check sensor.mask is a binary mask
+    if sensor.mask.dtype != bool and set(np.unique(sensor.mask).tolist()) != {0, 1}:
+        raise ValueError('The sensor must be defined as a binary mask.')
+
+    # find the coordinates of the sensor points
+    x_sensor = matlab_mask(kgrid.x, sensor.mask == 1)
+    x_sensor = np.squeeze(x_sensor)
+    y_sensor = matlab_mask(kgrid.y, sensor.mask == 1)
+    y_sensor = np.squeeze(y_sensor)
+
+    # find the angle of each sensor point (from the centre)
+    angle = np.arctan2(-x_sensor, -y_sensor)
+    angle[angle < 0] = 2 * np.pi + angle[angle < 0]
+
+    # sort the sensor points in order of increasing angle
+    indices_new = np.argsort(angle, kind='stable')
+
+    # reorder the measure time series so that adjacent time series correspond
+    # to adjacent sensor points.
+    reordered_sensor_data = sensor_data[indices_new]
+    return reordered_sensor_data
+
+
 def reorder_binary_sensor_data(sensor_data: np.ndarray, reorder_index: np.ndarray):
     """
     Args:
@@ -605,48 +649,6 @@ def unmask_sensor_data(kgrid, sensor, sensor_data: np.ndarray) -> np.ndarray:
     return unmasked_sensor_data
 
 
-def reorder_sensor_data(kgrid, sensor, sensor_data: np.ndarray) -> np.ndarray:
-    """
-    Reorders the sensor data based on the coordinates of the sensor points.
-
-    Args:
-        kgrid: The k-Wave grid object.
-        sensor: The k-Wave sensor object.
-        sensor_data: The sensor data to be reordered.
-
-    Returns:
-        np.ndarray of the reordered sensor data.
-
-    Raises:
-        ValueError: If the simulation is not 2D or the sensor is not defined as a binary mask.
-    """
-    # check simulation is 2D
-    if kgrid.dim != 2:
-        raise ValueError('The simulation must be 2D.')
-
-    # check sensor.mask is a binary mask
-    if sensor.mask.dtype != bool and set(np.unique(sensor.mask).tolist()) != {0, 1}:
-        raise ValueError('The sensor must be defined as a binary mask.')
-
-    # find the coordinates of the sensor points
-    x_sensor = matlab_mask(kgrid.x, sensor.mask == 1)
-    x_sensor = np.squeeze(x_sensor)
-    y_sensor = matlab_mask(kgrid.y, sensor.mask == 1)
-    y_sensor = np.squeeze(y_sensor)
-
-    # find the angle of each sensor point (from the centre)
-    angle = np.arctan2(-x_sensor, -y_sensor)
-    angle[angle < 0] = 2 * np.pi + angle[angle < 0]
-
-    # sort the sensor points in order of increasing angle
-    indices_new = np.argsort(angle, kind='stable')
-
-    # reorder the measure time series so that adjacent time series correspond
-    # to adjacent sensor points.
-    reordered_sensor_data = sensor_data[indices_new]
-    return reordered_sensor_data
-
-
 def create_cw_signals(t_array: np.ndarray, freq: float, amp: np.ndarray, phase: np.ndarray,
                       ramp_length: int = 4) -> np.ndarray:
     """
diff --git a/tests/matlab_test_data_collectors/matlab_collectors/collect_at_linear_array_transducer.m b/tests/matlab_test_data_collectors/matlab_collectors/collect_at_linear_array_transducer.m
new file mode 100644
index 00000000..ddafd7cd
--- /dev/null
+++ b/tests/matlab_test_data_collectors/matlab_collectors/collect_at_linear_array_transducer.m
@@ -0,0 +1,131 @@
+% derived from linear_array_transducer script of k-Wave MATLAB
+
+% create data recorder
+recorder = utils.TestRecorder('collectedValues/linear_array_transducer.mat');
+
+% =========================================================================
+% DEFINE LITERALS
+% =========================================================================
+    
+% select which k-Wave code to run
+%   1: MATLAB CPU code
+%   2: MATLAB GPU code
+%   3: C++ code
+%   4: CUDA code
+model           = 4;
+
+% medium parameters
+c0              = 1500;     % sound speed [m/s]
+rho0            = 1000;     % density [kg/m^3]
+
+% source parameters
+source_f0       = 1e6;      % source frequency [Hz]
+source_amp      = 1e6;      % source pressure [Pa]
+source_cycles   = 5;        % number of toneburst cycles
+source_focus    = 20e-3;    % focal length [m]
+element_num     = 15;       % number of elements
+element_width   = 1e-3;     % width [m]
+element_length  = 10e-3;    % elevation height [m]
+element_pitch   = 2e-3;     % pitch [m]
+
+% transducer position
+translation     = [5e-3, 0, 8e-3];
+rotation        = [0, 20, 0];
+
+% grid parameters
+grid_size_x     = 40e-3;    % [m]
+grid_size_y     = 20e-3;    % [m]
+grid_size_z     = 40e-3;    % [m]
+
+% computational parameters
+ppw             = 3;        % number of points per wavelength
+t_end           = 35e-6;    % total compute time [s]
+cfl             = 0.5;      % CFL number
+
+% =========================================================================
+% RUN SIMULATION
+% =========================================================================
+
+% --------------------
+% GRID
+% --------------------
+
+% calculate the grid spacing based on the PPW and F0
+dx = c0 / (ppw * source_f0);   % [m]
+
+% compute the size of the grid
+Nx = roundEven(grid_size_x / dx);
+Ny = roundEven(grid_size_y / dx);
+Nz = roundEven(grid_size_z / dx);
+
+% create the computational grid
+kgrid = kWaveGrid(Nx, dx, Ny, dx, Nz, dx);
+
+% create the time array
+kgrid.makeTime(c0, cfl, t_end);
+
+recorder.recordObject('kgrid', kgrid);
+
+% --------------------
+% SOURCE
+% --------------------
+
+% set indices for each element
+if rem(element_num, 2)
+    ids = (1:element_num) - ceil(element_num/2);
+else
+    ids = (1:element_num) - (element_num + 1)/2;
+end
+
+% set time delays for each element to focus at source_focus
+time_delays = -(sqrt((ids .* element_pitch).^2 + source_focus.^2) - source_focus) ./ c0;
+time_delays = time_delays - min(time_delays);
+
+% create time varying source signals (one for each physical element)
+source_sig = source_amp .* toneBurst(1/kgrid.dt, source_f0, source_cycles, 'SignalOffset', round(time_delays / kgrid.dt));
+
+% create empty kWaveArray
+karray = kWaveArray('BLITolerance', 0.05, 'UpsamplingRate', 10);
+
+% add rectangular elements
+for ind = 1:element_num
+    
+    % set element y position
+    x_pos = 0 - (element_num * element_pitch / 2 - element_pitch / 2) + (ind - 1) * element_pitch;
+    
+    % add element (set rotation angle to match the global rotation angle)
+    karray.addRectElement([x_pos, 0, kgrid.z_vec(1)], element_width, element_length, rotation);
+    
+end
+
+% move the array
+karray.setArrayPosition(translation, rotation)
+
+recorder.recordObject('karray', karray);
+% assign binary mask
+source.p_mask = karray.getArrayBinaryMask(kgrid);
+
+% assign source signals
+source.p = karray.getDistributedSourceSignal(kgrid, source_sig);
+    
+% --------------------
+% MEDIUM
+% --------------------
+
+% assign medium properties
+medium.sound_speed = c0;
+medium.density = rho0;
+
+% --------------------
+% SENSOR
+% --------------------
+
+% set sensor mask to record central plane
+sensor.mask = zeros(Nx, Ny, Nz);
+sensor.mask(:, Ny/2, :) = 1;
+
+% record the pressure
+sensor.record = {'p_max'};
+
+recorder.recordObject('sensor', sensor);
+recorder.saveRecordsToDisk();
diff --git a/tests/matlab_test_data_collectors/matlab_collectors/collect_kWaveArray.m b/tests/matlab_test_data_collectors/matlab_collectors/collect_kWaveArray.m
index a073be2c..6e095717 100644
--- a/tests/matlab_test_data_collectors/matlab_collectors/collect_kWaveArray.m
+++ b/tests/matlab_test_data_collectors/matlab_collectors/collect_kWaveArray.m
@@ -38,6 +38,11 @@
 recorder.recordObject('kwave_array', kwave_array);
 recorder.increment();
 
+% Useful for testing addCustomElement in 3D
+kwave_array.addCustomElement(single([[1, 1, 1, 2, 2, 2, 3, 3, 3]; [1, 2, 3, 1, 2, 3, 1, 2, 3]; [0, 0, 0, 0, 0, 0, 0, 0, 0]]), 9, 2, 'custom_3d');
+recorder.recordObject('kwave_array', kwave_array);
+recorder.increment();
+
 % Useful for testing addRectElement in 3D
 kwave_array.addRectElement([12, -8, 0.3], 3, 4, [2, 4, 5]);
 recorder.recordObject('kwave_array', kwave_array);
@@ -61,6 +66,11 @@
 recorder.recordObject('kwave_array', kwave_array);
 recorder.increment();
 
+% Useful for testing addCustomElement in 2D
+kwave_array.addCustomElement(single([[1, 1, 1, 2, 2, 2, 3, 3, 3]; [1, 2, 3, 1, 2, 3, 1, 2, 3]]), 9, 1, 'custom_2d');
+recorder.recordObject('kwave_array', kwave_array);
+recorder.increment();
+
 % Useful for testing addRectElement in 2D
 kwave_array.addRectElement([12, -8], 3, 4, 2);
 recorder.recordObject('kwave_array', kwave_array);
diff --git a/tests/matlab_test_data_collectors/matlab_collectors/collect_make_cart_bowl.m b/tests/matlab_test_data_collectors/matlab_collectors/collect_makeCartBowl.m
similarity index 100%
rename from tests/matlab_test_data_collectors/matlab_collectors/collect_make_cart_bowl.m
rename to tests/matlab_test_data_collectors/matlab_collectors/collect_makeCartBowl.m
diff --git a/tests/matlab_test_data_collectors/matlab_collectors/collect_makeCartCircle.m b/tests/matlab_test_data_collectors/matlab_collectors/collect_makeCartCircle.m
index 2aac5558..882d044b 100644
--- a/tests/matlab_test_data_collectors/matlab_collectors/collect_makeCartCircle.m
+++ b/tests/matlab_test_data_collectors/matlab_collectors/collect_makeCartCircle.m
@@ -2,21 +2,24 @@
 all_params = {
     {5, 100, [1,1],  2*pi}, ...
     {1, 10, [5,5],  pi/4}, ...
-    {0.5, 10, [0,0],  2*pi}
+    {0.5, 10, [0,0],  2*pi}, ...
+    {0.05, 20, [0 0]}
     };  
-output_folder = 'collectedValues/makeCartCircle/';
+output_file = 'collectedValues/makeCartCircle.mat';
+recorder = utils.TestRecorder(output_file);
 
 for idx=1:length(all_params)
-        circle = makeCartCircle(all_params{idx}{1}, all_params{idx}{2}, all_params{idx}{3}, all_params{idx}{4});
         params = all_params{idx};
-        idx_padded = sprintf('%06d', idx - 1);
-
-        idx_padded = sprintf('%06d', idx - 1);
-        if ~exist(output_folder, 'dir')
-            mkdir(output_folder);
+        if length(all_params{idx}) < 4
+            circle = makeCartCircle(params{1}, params{2}, params{3});
+            params = {params{1}, params{2}, params{3}};
+        else
+            circle = makeCartCircle(params{1}, params{2}, params{3}, params{4});
+            params = {params{1}, params{2}, params{3}, params{4}};
         end
-        filename = [output_folder idx_padded '.mat'];
-        save(filename, 'params', 'circle');
+        recorder.recordVariable('params', params);
+        recorder.recordVariable('circle', circle);
+        recorder.increment()
 
 end
-disp('Done.')
+recorder.saveRecordsToDisk();
diff --git a/tests/matlab_test_data_collectors/matlab_collectors/collect_make_cart_disc.m b/tests/matlab_test_data_collectors/matlab_collectors/collect_makeCartDisc.m
similarity index 100%
rename from tests/matlab_test_data_collectors/matlab_collectors/collect_make_cart_disc.m
rename to tests/matlab_test_data_collectors/matlab_collectors/collect_makeCartDisc.m
diff --git a/tests/matlab_test_data_collectors/matlab_collectors/collect_make_cart_rect.m b/tests/matlab_test_data_collectors/matlab_collectors/collect_makeCartRect.m
similarity index 100%
rename from tests/matlab_test_data_collectors/matlab_collectors/collect_make_cart_rect.m
rename to tests/matlab_test_data_collectors/matlab_collectors/collect_makeCartRect.m
diff --git a/tests/matlab_test_data_collectors/matlab_collectors/collect_resize.m b/tests/matlab_test_data_collectors/matlab_collectors/collect_resize.m
index 73e5ae61..73eb3cd9 100644
--- a/tests/matlab_test_data_collectors/matlab_collectors/collect_resize.m
+++ b/tests/matlab_test_data_collectors/matlab_collectors/collect_resize.m
@@ -3,9 +3,7 @@
 output_file = 'collectedValues/resize.mat';
 recorder = utils.TestRecorder(output_file);
 
-% TODO: test also for nearest
-% interp_methods = {'nearest', 'linear'};
-interp_methods = {'linear'};
+interp_methods = {'nearest', 'linear'};
 
 % Generate a random 3D volume with dimensions Nx x Ny x Nz
 Nx = randi([10,20]);
@@ -30,4 +28,3 @@
 
 end
 recorder.saveRecordsToDisk(); 
-    
diff --git a/tests/matlab_test_data_collectors/python_testers/kWaveArray_test.py b/tests/matlab_test_data_collectors/python_testers/kWaveArray_test.py
index 2b79d09f..1e088153 100644
--- a/tests/matlab_test_data_collectors/python_testers/kWaveArray_test.py
+++ b/tests/matlab_test_data_collectors/python_testers/kWaveArray_test.py
@@ -2,6 +2,7 @@
 from pathlib import Path
 
 import numpy as np
+import pytest
 from kwave.kgrid import kWaveGrid
 
 from kwave.utils.kwave_array import kWaveArray
@@ -44,6 +45,22 @@ def test_kwave_array():
     check_kwave_array_equality(kwave_array, reader.expected_value_of('kwave_array'))
     reader.increment()
 
+    kwave_array.add_custom_element(
+        integration_points=np.array([[1, 1, 1, 2, 2, 2, 3, 3, 3], [1, 2, 3, 1, 2, 3, 1, 2, 3], [0, 0, 0, 0, 0, 0, 0, 0, 0]], dtype=np.float32),
+        measure=9,
+        element_dim=2, label='custom_3d'
+    )
+    check_kwave_array_equality(kwave_array, reader.expected_value_of('kwave_array'))
+
+    with pytest.raises(ValueError):
+        kwave_array.add_custom_element(
+            integration_points=np.array([[1, 1, 1, 2, 2, 2, 3, 3, 3], [1, 2, 3, 1, 2, 3, 1, 2, 3]], dtype=np.float32),
+            measure=9,
+            element_dim=2, label='custom_3d'
+        )
+
+    reader.increment()
+
     kwave_array.add_rect_element([12, -8, 0.3], 3, 4, [2, 4, 5])
     check_kwave_array_equality(kwave_array, reader.expected_value_of('kwave_array'))
     reader.increment()
@@ -61,6 +78,13 @@ def test_kwave_array():
     check_kwave_array_equality(kwave_array, reader.expected_value_of('kwave_array'))
     reader.increment()
 
+    kwave_array.add_custom_element(
+        np.array([[1, 1, 1, 2, 2, 2, 3, 3, 3], [1, 2, 3, 1, 2, 3, 1, 2, 3]], dtype=np.float32),
+        9, 1, label='custom_2d'
+    )
+    check_kwave_array_equality(kwave_array, reader.expected_value_of('kwave_array'))
+    reader.increment()
+
     # Useful for testing addRectElement in 2D
     kwave_array.add_rect_element([12, -8], 3, 4, 2)
     check_kwave_array_equality(kwave_array, reader.expected_value_of('kwave_array'))
diff --git a/tests/matlab_test_data_collectors/python_testers/makeCartCircle_test.py b/tests/matlab_test_data_collectors/python_testers/makeCartCircle_test.py
index e3a9af85..6b0782da 100644
--- a/tests/matlab_test_data_collectors/python_testers/makeCartCircle_test.py
+++ b/tests/matlab_test_data_collectors/python_testers/makeCartCircle_test.py
@@ -16,17 +16,26 @@ def test_makeCartCircle():
     for i in range(num_collected_values):
         filepath = os.path.join(collected_values_folder, f'{i:06d}.mat')
         recorded_data = loadmat(filepath)
-
-        radius, num_points, center, arc_angle = recorded_data['params'][0]
-        center = Vector(center[0])
-        num_points = num_points[0][0]
-        radius = radius[0][0]
-        arc_angle = arc_angle[0][0]
-        if np.isclose(arc_angle, 2 * np.pi):
-            arc_angle = 2 * np.pi
+        params = recorded_data['params'][0]
+
+        if len(params) == 4:
+            radius, num_points, center, arc_angle = params 
+            center = Vector(center[0])
+            num_points = num_points[0][0]
+            radius = radius[0][0]
+            arc_angle = arc_angle[0][0]
+            if np.isclose(arc_angle, 2 * np.pi):
+                arc_angle = 2 * np.pi
+
+            circle = make_cart_circle(radius, num_points, center, arc_angle)
+        else:
+            radius, num_points, center = params 
+            center = Vector(center[0])
+            num_points = num_points[0][0]
+            radius = radius[0][0]
+            circle = make_cart_circle(radius, num_points, center)
         expected_value = recorded_data['circle']
 
-        circle = make_cart_circle(radius, num_points, center, arc_angle)
 
         assert np.allclose(expected_value, circle)
 
diff --git a/tests/matlab_test_data_collectors/python_testers/test_linear_array_transducer.py b/tests/matlab_test_data_collectors/python_testers/test_linear_array_transducer.py
new file mode 100644
index 00000000..c3e8c97f
--- /dev/null
+++ b/tests/matlab_test_data_collectors/python_testers/test_linear_array_transducer.py
@@ -0,0 +1,63 @@
+from tests.matlab_test_data_collectors.python_testers.utils.record_reader import TestRecordReader
+from tests.matlab_test_data_collectors.python_testers.utils.check_equality import check_kwave_array_equality
+from tests.matlab_test_data_collectors.python_testers.utils.check_equality import check_kgrid_equality
+import numpy as np
+import os
+from pathlib import Path
+
+import kwave.data
+from kwave.kgrid import kWaveGrid
+from kwave.utils.kwave_array import kWaveArray
+
+
+def test_linear_array_transducer():
+    test_record_path = os.path.join(Path(__file__).parent, 'collectedValues/linear_array_transducer.mat')
+    reader = TestRecordReader(test_record_path)
+
+    c0 = 1500
+    source_f0 = 1e6
+    source_focus = 20e-3
+    element_num = 15
+    element_width = 1e-3
+    element_length = 10e-3
+    element_pitch = 2e-3
+    translation = kwave.data.Vector([5e-3, 0, 8e-3])
+    rotation = kwave.data.Vector([0, 20, 0])
+    grid_size_x = 40e-3
+    grid_size_y = 20e-3
+    grid_size_z = 40e-3
+    ppw = 3
+    t_end = 35e-6
+    cfl = 0.5
+    bli_tolerance = 0.05
+    upsampling_rate = 10
+
+    # GRID
+    dx = c0 / (ppw * source_f0)
+    Nx = round(grid_size_x / dx)
+    Ny = round(grid_size_y / dx)
+    Nz = round(grid_size_z / dx)
+    kgrid = kWaveGrid([Nx, Ny, Nz], [dx, dx, dx])
+    kgrid.makeTime(c0, cfl, t_end)
+
+    check_kgrid_equality(kgrid, reader.expected_value_of('kgrid'))
+    # SOURCE
+    if element_num % 2 != 0:
+        centering_offset = np.ceil(element_num / 2)
+    else:
+        centering_offset = (element_num + 1) / 2
+
+    positional_basis = np.arange(1, element_num + 1) - centering_offset
+    
+    time_delays = -(np.sqrt((positional_basis * element_pitch) ** 2 + source_focus ** 2) - source_focus) / c0
+    time_delays = time_delays - min(time_delays)
+
+    karray = kWaveArray(bli_tolerance=bli_tolerance, upsampling_rate=upsampling_rate)
+
+    for ind in range(1, element_num + 1):
+        x_pos = 0 - (element_num * element_pitch / 2 - element_pitch / 2) + (ind - 1) * element_pitch
+        karray.add_rect_element([x_pos, 0, kgrid.z_vec.flat[0]], element_width, element_length, rotation)
+
+    karray.set_array_position(translation, rotation)
+
+    check_kwave_array_equality(karray, reader.expected_value_of('karray'))
\ No newline at end of file
diff --git a/tests/matlab_test_data_collectors/python_testers/test_resize.py b/tests/matlab_test_data_collectors/python_testers/test_resize.py
index b95bbbae..422f5098 100644
--- a/tests/matlab_test_data_collectors/python_testers/test_resize.py
+++ b/tests/matlab_test_data_collectors/python_testers/test_resize.py
@@ -22,5 +22,6 @@ def test_resize():
 
         assert np.allclose(expected_resized_volume,
                            resized_volume), f"Results do not match for {i + 1} dimensional case."
+        reader.increment()
 
-    logging.log(logging.INFO, 'revolve2d(..) works as expected!')
+    logging.log(logging.INFO, 'resize(..) works as expected!')
diff --git a/tests/test_utils.py b/tests/test_utils.py
index 9f1aa81a..dd9266fd 100644
--- a/tests/test_utils.py
+++ b/tests/test_utils.py
@@ -74,12 +74,23 @@ def test_spect():
 
 
 def test_extract_amp_phase():
+    # odd length signal
     test_signal = tone_burst(sample_freq=10_000_000, signal_freq=2.5 * 1_000_000, num_cycles=2, envelope='Gaussian')
+    assert( np.shape(np.squeeze(test_signal))[0] % 2 != 0)
     a_t, b_t, c_t = extract_amp_phase(data=test_signal, Fs=10_000_000, source_freq=2.5 * 10 ** 6)
     a, b, c = 0.6547, -1.8035, 2.5926e06
     assert (abs(a_t - a) < 0.01).all()
     assert (abs(b_t - b) < 0.0001).all()
     assert (abs(c_t - c) < 100).all()
+    # even length signal
+    test_signal = tone_burst(sample_freq=18_000_000, signal_freq=6_000_000, num_cycles=5, envelope='Gaussian')
+    assert( np.shape(np.squeeze(test_signal))[0] % 2 == 0)
+    a_t, b_t, c_t = extract_amp_phase(data=test_signal, Fs=18_000_000, source_freq=6_000_000)
+    a, b, c = 0.6591, -1.5708, 6.000e06
+    assert (abs(a_t - a) < 0.01).all()
+    assert (abs(b_t - b) < 0.0001).all()
+    assert (abs(c_t - c) < 100).all()
+
 
 
 def test_apply_filter_lowpass():