I'm new to this toolbox and I'm working through some of the simulatons. I tried copy pasting the examples and the ones I tried worked until I tried to run the Single Slit Diffusion. Could anyone look at my code and tell me what I'm doing wrong?
=====================================================================
%%%%% load homogenous Propogation Medium
% define the properties used in the simulation
p0 = 7*10^6; % source pressure [Pa]
c0 = 1500; % sound speed [m/s]
rho0 = 1000; % density [kg/m^3]
alpha_0 = 0.25; % absorption coefficient [dB/(MHz^2 cm)]
sigma = 1; % shock parameter
source_freq = 1e6; % frequency [Hz]
points_per_wavelength = 50; % number of grid points per wavelength at f0
wavelength_separation = 15; % separation between the source and detector
PML_size = 64; % PML size
CFL = 0.25; %
% create the computational grid
Nx = 128; % number of grid points in the x (row) direction
Ny = 128; % number of grid points in the y (column) direction
dx = 0.1e-3; % grid point spacing in the x direction [m]
dy = 0.1e-3; % grid point spacing in the y direction [m]
kgrid = makeGrid(Nx, dx, Ny, dy);
% define the properties of the propagation medium
medium.sound_speed = 1500; % [m/s]
medium.alpha_coeff = 0.75; % [dB/(MHz^y cm)]
medium.alpha_power = 1.5;
% define the ratio between the barrier and background sound speed and density
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; % [grid points]
slit_width = 10; % [grid points]
slit_x_pos = Nx - Nx/4; % [grid points]
slit_offset = Ny/2 - slit_width/2 - 1; % [grid points]
slit_mask = zeros(Nx, Ny);
slit_mask(slit_x_pos:slit_x_pos + slit_thickness, 1:1 + slit_offset) = 1;
slit_mask(slit_x_pos:slit_x_pos + slit_thickness, end - slit_offset:end) = 1;
% 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;
% define a centered circular sensor
sensor_radius = 4e-3; % [m]
num_sensor_points = 50;
sensor.mask = makeCartCircle(sensor_radius, num_sensor_points);
% create time array
t_end = 3e-5;
kgrid.t_array = makeTime(kgrid, medium.sound_speed, [], t_end);
% define four sensor points centered about source.p0
sensor_radius = 40; % [grid points]
sensor.mask = zeros(Nx, Ny);
sensor.mask(Nx/2 + sensor_radius, Ny/2) = 1;
sensor.mask(Nx/2 - sensor_radius, Ny/2) = 1;
sensor.mask(Nx/2, Ny/2 + sensor_radius) = 1;
sensor.mask(Nx/2, Ny/2 - sensor_radius) = 1;
% define a single source point
source.p_mask = zeros(Nx, Ny);
source.p_mask(end - Nx/4, Ny/2) = 1;
% define a time varying sinusoidal source
source_freq = 0.25e6; % [Hz]
source_mag = 2; % [Pa]
source.p = source_mag*sin(2*pi*source_freq*kgrid.t_array);
% filter the source to remove high frequencies not supported by the grid
source.p = filterTimeSeries(kgrid, medium, source.p);
% define a curved transducer element
source.p_mask = makeCircle(Nx, Ny, 61, 61, 60, pi/2);
% define a time varying sinusoidal source
source_freq = 0.25e6; % [Hz]
source_mag = 0.5; % [Pa]
source.p = source_mag*sin(2*pi*source_freq*kgrid.t_array);
% filter the source to remove any high frequencies not supported by the grid
source.p = filterTimeSeries(kgrid, medium, source.p);
% create a sensor mask covering the entire computational domain using the
% opposing corners of a rectangle
sensor.mask = [1, 1, Nx, Ny].';
% set the record mode to capture the final wave-field and the statistics at
% each sensor point
sensor.record = {'p_final', 'p_max', 'p_rms'};
% set the input options
input_args = {'PMLInside', false, 'PMLSize', PML_size, 'PlotPML', false, ...
'DisplayMask', slit_mask, 'DataCast', 'single'};
% run the simulation
sensor_data = kspaceFirstOrder2D(kgrid, medium, source, sensor, input_args{:});
================================================================
Thank you!
Pauline