Dear Bradley:

Thank you for your help. I am trying to reconstruct a disc medium density model with two steps: 1. forward modeling (a. with disc density model, b. without disc density model, c. data is a-b), and 2. time reversal to reconstruct the medium density model. I tried to use a pulse source with the Mexican hat wavelet. I modified the 2D Time Reversal Reconstruction For A Line Sensor Example as below. I keep all the timing parameters same in forward and reversal modeling. However, the reconstructed image is not same as the disc density model (it seems reconstructed time is a little longer than it should be). I was wondering if you could please help with my questions?

My questions are as follows:

1. The time parameters in reconstruction are kept same as those in the forward modeling, could you please help with the possible reason why the reconstructed image is not the same as the medium disc2? How to restore the exact the original medium disc2 shape (e.g. do I need to control the timing)?

2. the reconstructed image 'density_recon' values have both positive and negative values,

why there are negative image values (the medium disc2 are postive values as medium density? Why the original disc values can't be restored?

"density_recon = kspaceFirstOrder2D(kgrid, medium, source, sensor, input_args{:}); % inverted medium parameters"

Thank you very much for your help in advance and best regards,

John

p.s. the Matlab script is as below:

--------------------------------------------------------------------------------

clear all;

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

% SET UP THE SIMULATION

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

% create the computational grid

PML_size = 20; % size of the PML in grid points

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

Ny = 256 - 2*PML_size; % 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);

% create the time array

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

% create initial pressure distribution using makeDisc

disc_magnitude = 5; % [au]

disc_x_pos = 10; % [grid points]

disc_y_pos = floor(Ny/2); % [grid points]

disc_radius = 5; % [grid points]

disc_1 = disc_magnitude*makeDisc(Nx, Ny, disc_x_pos, disc_y_pos, disc_radius);

% smooth the initial pressure distribution and restore the magnitude

% scale the initial pressure distribution

p0_magnitude = 1;

p0 = smooth(kgrid, disc_1, true);

% assign to the source structure

source_freq = 10;

source.p_mask = disc_1;

source.p = p0_magnitude .* (1 - 2 * pi^2 * source_freq^2 .* (kgrid.t_array).^2)...

.* exp(-pi^2 * source_freq^2 .* (kgrid.t_array).^2); % Ricker wavelet

% define the properties of the propagation medium

disc_x_pos = floor(Nx/2); % [grid points]

disc_y_pos = floor(Ny/2); % [grid points]

disc_radius = 5; % [grid points]

disc_2 = makeDisc(Nx, Ny, disc_x_pos, disc_y_pos, disc_radius);

p3 = smooth(kgrid, disc_2, true);

density_input = 1000*ones(Nx, Ny); % [kg/m^3]

medium.sound_speed = 1500 ; % [m/s]

medium.density = density_input + 1000 * disc_2; % [m/s]

% set the input arguements: force the PML to be outside the computational

% grid; switch off p0 smoothing within kspaceFirstOrder2D

input_args = {'PMLInside', false, 'PMLSize', PML_size, 'Smooth', false, 'PlotPML', false};

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

% DEFINE SENSOR INDICES FOR LATER USE

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

% define a four-sided box array and find the indices of the sensors

sensor.mask = zeros(Nx, Ny);

sensor.mask(1,:) = 1;

sensor.mask(:,1) = 1;

sensor.mask(end,:) = 1;

sensor.mask(:,end) = 1;

sensor_indices = find(sensor.mask==1);

% find the indices along two sides, with respect to sensor_indices

sensor.mask = zeros(Nx, Ny);

sensor.mask(end,:) = 1;

sensor.mask(:,end) = 1;

sensor_box_indices1 = find(sensor.mask(sensor_indices)==1);

% define an L-shaped array used for the measurements, find the sensor

% indices, and their indices with respect to sensor_indices

sensor.mask = zeros(Nx, Ny);

sensor.mask(1,:) = 1;

sensor.mask(:,1) = 1;

original_sensor_indices = find(sensor.mask == 1);

sensor_box_indices2 = find(sensor.mask(sensor_indices)==1);

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

% RUN THE SIMULATION AND RECORD THE DATA OVER THE L-SHAPED ARRAY

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

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

medium.density = density_input; % [m/s]

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

sensor_data = sensor_data_object - sensor_data_no_object;

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

% RECONSTRUCT AN IMAGE USING TIME REVERSAL AND LINE ARRAY DATA

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

% define a line array of sensors

sensor.mask = zeros(Nx, Ny);

sensor.mask(1, :) = 1;

% assign the time reversal data

sensor_data = sensor_data; % - sensor_data_no_object;

sensor.time_reversal_boundary_data = 1 * sensor_data((sensor.mask(original_sensor_indices) == 1),:);

% set the initial pressure to be zero

source.p = source.p * 0;

% run the time reversal reconstruction

density_recon = kspaceFirstOrder2D(kgrid, medium, source, sensor, input_args{:}); % inverted medium parameters

figure;

imagesc(disc_2)

title('medium density');

colorbar;

figure;

imagesc(sensor_data)

title('sensor_data');

colorbar;

figure;

imagesc(density_recon)

title('density_recon');

colorbar;