Hi Dr.Treeby,
I have running a simple code shown below. I have set ample simulation time, and made sure the all the other parameters are met for stable simulations. But there seems to be an issue towards the end of my simulation. The reflected waves are not absorbed by the PML layer close to the sensor mask. I have tried increasing the PMLAlpha value but didn't see any difference. I'm afraid it might be affecting the results I am supposed to be getting. Any help would be appreciated!
Thanks, Susan
% create the computational grid
clear all;
pml_x_size = 20;
pml_y_size = 20;
Nx = 256 - 2*pml_x_size;
Ny = 96 - 2*pml_y_size;
dx = 0.2e-3; % grid point spacing in the x direction [m]
dy = 0.2e-3; % grid point spacing in the y direction [m]
kgrid = makeGrid(Nx, dx, Ny, dy);
gridx_pts =11;
c0 = 1540; % [m/s]
rho0 = 1000; % [kg/m^3]
medium.sound_speed_ref = c0;
medium.alpha_coeff = 0.75; % [dB/(MHz^y cm)]
medium.alpha_power = 1.5;
%medium.alpha_mode = 'no_absorption';
t_end = (Nx*dx)*2.4/c0; % [s]
kgrid.t_array = 0:15e-9:t_end;
source_strength = 1e6; % [Pa]
tone_burst_freq = 1.875e6; % [Hz]
tone_burst_cycles = 5;
input_signal = source_strength*toneBurst(1/kgrid.dt, tone_burst_freq, tone_burst_cycles);
source.p_mask = zeros(Nx,Ny);
num_elements = 32; % [grid points]
x_offset = 1 ; % [grid points]
start_index = Ny/2 - round(num_elements/2) + 1;
source.p_mask(x_offset, start_index:start_index + num_elements - 1) = 1;
source.p = input_signal;
sensor.mask = zeros(Nx,Ny);
sensor.mask(x_offset, start_index) = 1;
sensor.mask(x_offset,18) = 1;
%%
Ny_tot = Ny+(41);
Nx_tot = Nx;
sound_speed_map = c0*ones(Nx_tot, Ny_tot);
density_map = rho0*ones(Nx_tot, Ny_tot);
cwire = 3000;
scattering_rho0 = cwire/1.5;
% define the fishing wire point
radius = 0.4e-3; % [m]
xc = Nx_tot - 5;
yc = floor(Ny_tot/2);
%scattering_region = makeCircle(Nx_tot, Ny_tot, xc, yc,radius,'true');
scattering_region = makeDisc(Nx_tot, Ny_tot, xc, yc,1,'true'); % results
%are not clean if you use makeDisc
% assign scatterer the properties
sound_speed_map(scattering_region == 1) = cwire;
density_map(scattering_region == 1) = scattering_rho0;
%%
input_args = {'PMLInside', false, 'PMLSize', [pml_x_size, pml_y_size],'PlotSim',false,'PlotScale','auto'};
psny = 1;
medium.sound_speed = sound_speed_map(:,psny:psny+Ny-1);
medium.density = density_map(:,psny:psny+Ny-1);
sensor_data = kspaceFirstOrder2D(kgrid, medium, source, sensor,input_args{:});