From c9315541ec81ffd17e473c4b41c652e6a0cc549d Mon Sep 17 00:00:00 2001
From: Walter Simson <walter.simson@tum.de>
Date: Thu, 2 Nov 2023 12:28:57 -0700
Subject: [PATCH] Example of karray source (#193)

* WIP commit with two examples

* Working source example.

 - Pressure field is transposed. (might be related to open issue #173)
 - Needs to be cleaned up.

* sensor example running

* Change from open PR for 2D sim io

* Update branch to pass CI

* Pass Ruff CI

* Minor reformatting and adding future todos

* Add testing for linear_array_transducer setup

* Fix array as sensor example

* turn off simulation in tests

* how did that get in there?

* reorder sensor data to complete example

* pass ruff

* fix signatures and docs

* passes CI?

* make parse_executable_output static method

* clean up test_linear_array.py

* remove license header

* Use vectors

* Update checks to use Iterable type

* fix errors

* update variable names to clarify code

* clean up imports and magic numbers

* clean up commented code

* delete unused code

* clean up example file after review

* update linear array example indexing

* Clean up make cart circle test

* clean up ruff error

* clean up unused import for ruff
---
 examples/bmode_reconstruction_example.py      |   9 +-
 examples/example_array_as_a_sensor.py         | 149 +++++++++++++++
 examples/example_array_as_a_source.py         | 166 +++++++++++++++++
 .../example_at_linear_array_transducer.py     | 176 +++++++++---------
 kwave/executor.py                             |   3 +-
 kwave/kWaveSimulation.py                      |  65 +++----
 .../create_storage_variables.py               |   2 +-
 kwave/kgrid.py                                |   9 +
 kwave/kspaceFirstOrder2D.py                   |  12 +-
 kwave/utils/conversion.py                     |  16 +-
 kwave/utils/io.py                             |   2 +-
 kwave/utils/kwave_array.py                    |  24 ++-
 kwave/utils/signals.py                        |  84 ++++-----
 .../collect_at_linear_array_transducer.m      | 131 +++++++++++++
 ...ake_cart_bowl.m => collect_makeCartBowl.m} |   0
 .../collect_makeCartCircle.m                  |  25 +--
 ...ake_cart_disc.m => collect_makeCartDisc.m} |   0
 ...ake_cart_rect.m => collect_makeCartRect.m} |   0
 .../python_testers/makeCartCircle_test.py     |  27 ++-
 .../test_linear_array_transducer.py           |  63 +++++++
 20 files changed, 751 insertions(+), 212 deletions(-)
 create mode 100644 examples/example_array_as_a_sensor.py
 create mode 100644 examples/example_array_as_a_source.py
 create mode 100644 tests/matlab_test_data_collectors/matlab_collectors/collect_at_linear_array_transducer.m
 rename tests/matlab_test_data_collectors/matlab_collectors/{collect_make_cart_bowl.m => collect_makeCartBowl.m} (100%)
 rename tests/matlab_test_data_collectors/matlab_collectors/{collect_make_cart_disc.m => collect_makeCartDisc.m} (100%)
 rename tests/matlab_test_data_collectors/matlab_collectors/{collect_make_cart_rect.m => collect_makeCartRect.m} (100%)
 create mode 100644 tests/matlab_test_data_collectors/python_testers/test_linear_array_transducer.py

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/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/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..5cba1506 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:
 
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/io.py b/kwave/utils/io.py
index 8016043e..6a3a7a84 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():
diff --git a/kwave/utils/kwave_array.py b/kwave/utils/kwave_array.py
index 6a46db38..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
@@ -303,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"
 
@@ -436,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
 
@@ -642,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))
 
@@ -828,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)
@@ -860,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 f60354e1..21660c42 100644
--- a/kwave/utils/signals.py
+++ b/kwave/utils/signals.py
@@ -432,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:
@@ -607,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_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/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