k-Wave Toolbox |

## Diffraction Through A Slit Example

On this page…
Overview Defining the medium properties Running the simulation Diffraction through a double slit

### Overview

This example illustrates the diffraction of a plane acoustic wave through a slit. It builds on the Monopole Point Source In A Homogeneous Propagation Medium and Simulating Transducer Field Patterns examples.

- open the file in the MATLAB Editor
- run the file in MATLAB

Back to Top

### Defining the medium properties

The current version of k-Wave cannot explicitly enforce impedance boundary conditions. However, rigid boundaries can still be modelled by assigning a boundary with a large sound speed and density compared to the rest of the propagation medium. As this impedance difference is increased, the behaviour of the interface will approach that of a rigid boundary. However, increasing the sound speed of the boundary will similarly increase the computational requirements (the time step required for a stable simulation is dependent on the maximum sound speed within the medium).

Here, a diffraction slit is created by defining a thin layer with a significanly different sound speed and density to the background medium such that most of the incident wave is reflected. (Note, using `barrier_scale = 20`, some of the incident wave will still be transmitted through the barrier causing additional interference fringes. These can be reduced by increasing the barrier scale.)

    % define the ratio between the barrier and background c and rho
    barrier_scale = 20;

    % create the time array using the barrier sound speed
    t_end = 40e-6;                % [s]
    CFL = 0.5;                    % Courant-Friedrichs-Lewy number
    kgrid.t_array = makeTime(kgrid, c0*barrier_scale, CFL, t_end);

    % create a mask of a barrier with a slit
    slit_thickness = 2;           % [pixels]
    slit_width = 10;              % [pixels]
    slit_z_pos = Nz - Nz/4;
    slit_offset = Nx/2 - slit_width/2 - 1;
    slit_mask = zeros(Nz, Nx);
    slit_mask(slit_z_pos:slit_z_pos + slit_thickness, 1:1 + slit_offset) = 1;
    slit_mask(slit_z_pos:slit_z_pos + slit_thickness, end - slit_offset:end) = 1;

    % assign the slit to the properties of the propagation medium
    medium.sound_speed = c0*ones(Nz, Nx);
    medium.density = rho0*ones(Nz, Nx);
    medium.sound_speed(slit_mask == 1) = barrier_scale*c0;
    medium.density(slit_mask == 1) = barrier_scale*rho0;

Back to Top

### Running the simulation

In the first example (set `example_number = 1` within the example m-file), a sinusoidal plane wave source is created with a wavelength equal to the slit width. A visualisation of the barrier is also produced by assigning the `slit_mask` created above to the optional input parameter `'DisplayMask'`. The size and location of the perfectly matched layer are explicity defined (see Controlling The Absorbing Boundary Layer Example), the simulation is run in `single` precision to reduce the computation time (see Optimising k-Wave Performance Example), and the velocity data is also returned by setting the optional input parameter `'ReturnVelocity'` to `true` (see Recording The Particle Velocity Example).

    % set the input options and run the simulation
    input_args = {'ReturnVelocity', true, 'PMLInside', false, 'PMLSize', PML_size,...
        'DisplayMask', slit_mask, 'DataCast', 'single'};
    [sensor_data, field_data] = kspaceFirstOrder2D(kgrid, medium, source, sensor, input_args{:});

A visualisation of the pressure field is given below. The wavefield on the far side of the slit appears reasonably omni-directional.

In the second example (set `example_number = 2` within the example m-file), the slit size is increased and the wavelength is reduced to be a quarter of the slit width. A visualisation of the pressure field is given below. In this case, the directionality of the wavefield is increased, and interference fringes are visible.

Back to Top

### Diffraction through a double slit

In the third example (set `example_number = 3` within the example m-file), a double slit is used and the wavelength set equal to the slit width.

The particle velocity can also be returned by setting the optional input parameter `'ReturnVelocity'` to `true`.

Note, the pictures shown here were computed using `barrier_scale = 50` and the `scale = 2` within the example.

Back to Top

© 2009, 2010, 2011 Bradley Treeby and Ben Cox.