Hi, Thank you very much for this code, its really helpfull!

I have concerns about using it in a specific calculation. The idea is to simulate the propagation of ultrasound at very low frequency (8kHz) inside a tube filled with water. There are two opposite sources on the side of the tube. Everything is supposed to be symetric here but the two sources behave very differently. How do you exlain the asymmetry in the calculation? Is there something wrong in my code?

'Nx = 101; % number of grid points in the x (row) direction

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

dx = 0.0005; % grid point spacing in the x direction [m]

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

kgrid = makeGrid(Nx, dx, Ny, dy);

% define the properties of the propagation medium

ray_int=18e-3; % [m]

medium.sound_speed_compression = 330*ones(Nx, Ny); % [m/s]

medium.sound_speed_compression(makeDisc(Nx,Ny,(Nx+1)/2,(Ny+1)/2,ray_int/dx+3)==1) = 5650;

medium.sound_speed_compression(makeDisc(Nx,Ny,(Nx+1)/2,(Ny+1)/2,ray_int/dx)==1) = 1480;

medium.sound_speed_shear = 0*ones(Nx, Ny); % [m/s]

medium.sound_speed_shear(makeDisc(Nx,Ny,(Nx+1)/2,(Ny+1)/2,ray_int/dx+3)==1) = 3100;

medium.sound_speed_shear(makeDisc(Nx,Ny,(Nx+1)/2,(Ny+1)/2,ray_int/dx)==1) = 0;

medium.density = 2*ones(Nx, Ny); % [kg/m^3]

medium.density(makeDisc(Nx,Ny,(Nx+1)/2,(Ny+1)/2,ray_int/dx+3)==1) = 8010;

medium.density(makeDisc(Nx,Ny,(Nx+1)/2,(Ny+1)/2,ray_int/dx)==1) = 1000;

% define a curved transducer element

source.u_mask = zeros(Nx,Ny);

source.u_mask((Nx+1)/2-20:(Nx+1)/2+20,(Nx+1)/2-ray_int/dx-2)=1;% create the time array

source.u_mask((Nx+1)/2-20:(Nx+1)/2+20,(Nx+1)/2+ray_int/dx+2)=1;% create the time array

% source.u_mask((Nx+1)/2-20:(Nx+1)/2+20,(Nx+1)/2-ray_int/dx-2)=1;% create the time array

% source.u_mask((Nx+1)/2-20:(Nx+1)/2+20,(Nx+1)/2+ray_int/dx+1)=1;% create the time array

% define a time varying sinusoidal source

source_freq1 = 8e3; % [Hz]

source_freq2 = 8e3; % [Hz]

source_mag = 0.03

t_end=1/source_freq2;

% create the time array

[kgrid.t_array, dt] = makeTime(kgrid, 16000, [], t_end);

signal_uy_plus = source_mag*sin(2*pi*source_freq1*kgrid.t_array);

signal_uy_moins = -source_mag*sin(2*pi*source_freq2*kgrid.t_array);

% filter the source to remove any high frequencies not supported by the grid

source.uy = [repmat(filterTimeSeries(kgrid, medium, signal_uy_plus),41,1);repmat(filterTimeSeries(kgrid, medium, signal_uy_moins),41,1)];

display_mask = source.u_mask;

% define a centered circular sensor

% create a sensor mask covering the entire computational domain using the

% opposing corners of a rectangle

Pos_x=(Nx+1)/2;

Pos_y=(Ny+1)/2;

sensor.mask = zeros(Nx,Ny);

sensor.mask (Pos_x, :)=1;

sensor.mask (:,Pos_y)=1;

sensor.record={'p'};

input_args = {'DisplayMask', display_mask, 'PMLInside', false, 'PlotPML', false,'PlotScale', [-100e3,100e3,-100e3,100e3],};

sensor_data = pstdElastic2D(kgrid, medium, source, sensor, input_args{:});'

I supposed the definition of the sources is not correct. By changing the location of the sources (commented lines), the results are slightly more symetric..

Subsidiary question: Would it be better with a coupling between sources and tube?

Thank you in advance.

Damien