Hi everybody, I'm new to k-wave and I am trying to learn by example -- in this case I'm interested in simulating a single-edge (knife-edge) diffraction simulation.

I started with the built-in example "example_tvsp_slit_diffraction.m" that Bradley Treeby wrote, and simply removed one of the two sides of the barrier. When I do so, the simulation seems to become divergent. I tried removing the "Datacast = single" option and it still doesn't work.

If I remove all the barriers, the simulation runs fine ... basically a homogeneous media.

Can somebody help me understand why the one-sided simulation doesn't run correctly, please? Thanks in advance!

% Diffraction Through A Slit Example Example

%

% 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.

%

% author: Bradley Treeby

% date: 1st November 2010

% last update: 20th June 2017

%

% This function is part of the k-Wave Toolbox (http://www.k-wave.org)

% Copyright (C) 2010-2017 Bradley Treeby% This file is part of k-Wave. k-Wave is free software: you can

% redistribute it and/or modify it under the terms of the GNU Lesser

% General Public License as published by the Free Software Foundation,

% either version 3 of the License, or (at your option) any later version.

%

% k-Wave is distributed in the hope that it will be useful, but WITHOUT ANY

% WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS

% FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for

% more details.

%

% You should have received a copy of the GNU Lesser General Public License

% along with k-Wave. If not, see <http://www.gnu.org/licenses/>.clearvars;

% =========================================================================

% SIMULATION

% =========================================================================% create the computational grid and define the bulk medium properties

scale = 1; % change to 2 to produce higher resolution images

PML_size = 10; % size of the perfectly matched layer [grid points]

Nx = scale * 128 - 2 * PML_size; % number of grid points in the x (row) direction

Ny = Nx; % number of grid points in the y (column) direction

dx = 50e-3 / Nx; % grid point spacing in the x direction [m]

dy = dx; % grid point spacing in the y direction [m]

c0 = 1500; % [m/s] speed of sound in water

rho0 = 1000; % [kg/m^3]

kgrid = kWaveGrid(Nx, dx, Ny, dy);% define the ratio between the barrier and background sound speed and

% density

%barrier_scale = 20;

barrier_scale = 20; % lower contrast values = better diffraction contrast%% define the barrier and the source wavelength depending on the example

% create a mask of a barrier with a slit

slit_thickness = scale * 2; % [grid points]

slit_width = scale * 50; % [grid points]

slit_x_pos = Nx - Nx/4; % [grid points]

slit_y_end = Ny/2; % [grid points]

slit_offset = Ny/2 - slit_width/2 - 1; % [grid points]

slit_mask = zeros(Nx, Ny);

%%% one sided -- crashes?

slit_mask(slit_x_pos:slit_x_pos + slit_thickness, 1:slit_y_end) = 1; % left side only

%%% two slits -- works

% slit_mask(slit_x_pos:slit_x_pos + slit_thickness, 1:slit_offset) = 1; % left side

% slit_mask(slit_x_pos:slit_x_pos + slit_thickness, end - 1:end) = 1; % right side% define the source wavelength to be a quarter of the slit size

source_wavelength = 0.25 * slit_width * dx; % [m]

source_freq = 500e3; % [Hz]

source_wavelength = c0/source_freq; % [m]

% source_wavelength = 0.0058;% assign the slit to the properties of the propagation medium

medium.sound_speed = c0 * ones(Nx, Ny);

medium.density = rho0 * ones(Nx, Ny);

medium.sound_speed(slit_mask == 1) = barrier_scale * c0;

medium.density(slit_mask == 1) = barrier_scale * rho0;% assign the reference sound speed to the background medium

medium.sound_speed_ref = c0;%% prepare for simulation

% find the time step at the stability limit

c_ref = medium.sound_speed_ref;

c_max = barrier_scale * c0;

k_max = max(kgrid.k(:));

dt_limit = 2 / (c_ref * k_max) * asin(c_ref / c_max);% create the time array, with the time step just below the stability limit

dt = 0.95 * dt_limit; % [s]

t_end = 40e-6; % [s]

kgrid.setTime(round(t_end / dt) + 1, dt);% create a source mask of a single line

source.p_mask = zeros(Nx, Ny);

source.p_mask(end, :) = 1;% create and filter the time varying sinusoidal source

source_mag = 2;

source_freq = c0 / source_wavelength;

source.p = source_mag * sin(2 * pi * source_freq * kgrid.t_array);

source.p = filterTimeSeries(kgrid, medium, source.p);% define the field parameters to record

sensor.mask = ones(Nx, Ny);

sensor.record = {'u_final', 'p_final'};%% run simulation and display results

% set the input options

% input_args = {'PMLInside', false, 'PMLSize', PML_size, 'PlotPML', false, ...

% 'DisplayMask', slit_mask, 'DataCast', 'single'};

input_args = {'PMLInside', false, 'PMLSize', PML_size, 'PlotPML', false, ...

'DisplayMask', slit_mask};% run the simulation

sensor_data = kspaceFirstOrder2D(kgrid, medium, source, sensor, input_args{:});% =========================================================================

% VISUALISATION

% =========================================================================% plot the final wave-field

figure;

mx = max(abs(sensor_data.p_final(:)));

sensor_data.p_final(slit_mask == 1) = mx;

imagesc(kgrid.y_vec * 1e3, kgrid.x_vec * 1e3, sensor_data.p_final, [-mx, mx]);

colormap(getColorMap);

ylabel('x-position [mm]');

xlabel('y-position [mm]');

axis image;

title('p');% plot the final wave-field

figure;

mx = max(abs(sensor_data.ux_final(:)));

sensor_data.ux_final(slit_mask == 1) = mx;

imagesc(kgrid.y_vec * 1e3, kgrid.x_vec * 1e3, sensor_data.ux_final, [-mx, mx]);

colormap(getColorMap);

ylabel('x-position [mm]');

xlabel('y-position [mm]');

axis image;

title('ux');% plot the final wave-field

figure;

mx = max(abs(sensor_data.uy_final(:)));

sensor_data.uy_final(slit_mask == 1) = mx;

imagesc(kgrid.y_vec * 1e3, kgrid.x_vec * 1e3, sensor_data.uy_final, [-mx, mx]);

colormap(getColorMap);

ylabel('x-position [mm]');

xlabel('y-position [mm]');

axis image;

title('uy');