Hi,
I am trying to analyze the frequency attenuation of a photoacoustic point source as it propagates through the medium. Currently I have the source set up as a time varying function (a derivative of a Gaussian pulse of 7ps pulse duration)
I have a 3x3x1100 grid. I have the source at position (2, 2, 50), and sensors at points sensor.mask(2,2,50:1050) (i.e. there is a sensor at every point in the z direction starting at (2,2,50) and ending at (2,2,1050), The point spacing is dx=dy=dz=1e-9 %[m] in all directions.
Now, I run this simulation with a timestep of 1e-13 and t_end=2e-9 with sound speed 1480 in the medium and an alpha coefficient of 0.0022 and alpha power of 1.01.
I noticed some very strange effects here.
1) The recorded sensor signal is tripolar whereas the pressure vs time at the source is bipolar
source_mag = 2500; % [Pa] Max Amplitude
Pulsewidth=7e-12; %Pulse Width [s]
C=source_mag*Pulsewidth^2/(exp(-0.5)); %Constant adjusted so that peaks of
%pressure are equal to specified Source magnitude
source.p = (-C/Pulsewidth^3)*(kgrid.t_array-PhaseShift).*exp(-((kgrid.t_array-PhaseShift).^2)/(2*Pulsewidth^2));2) The FWHM of this signal is actually not consistently increasing as it gets recorded at further and further sensor points, and is displaying some sinusoidal behaviors.
3) The signal starts to have strange distortions introduced as it propagates through space.
Now I suspect that this strange behavior is due to my simulation paramaeters, but I have made the PML large in addition to decreasing the time step to ensure that the simulation is stable. nevertheless my output does not seem to be making sense. Can anyone provide some guidance as to why this could be happening? Are any of my parameters possible the problem? If anyone is interested the full source code is below.
Thank you in advance!
% create the computational grid
Nx = 3; % number of grid points in the x direction
Ny = 3; % number of grid points in the y direction
Nz = 1100; % number of grid points in the z direction
dx = 1e-9; % grid point spacing in the x direction [m]
dy = 1e-9; % grid point spacing in the y direction [m]
dz = 1e-10; % grid point spacing in the z direction [m]
kgrid = makeGrid(Nx, dx, Ny, dy, Nz, dz);
% define the properties of the propagation medium
medium.sound_speed = 1480*ones(Nx, Ny, Nz); % [m/s]
medium.density = 1000*ones(Nx, Ny, Nz); % [kg/m^3]
medium.alpha_coeff = 0.0022; % [dB/(MHz^y cm)]
medium.alpha_power = 1.01;
dt=1e-14;
t_end=2e-9;
kgrid.t_array=[0:dt:t_end];
PhaseShift=50e-12; %[s] Phase shift in time of pressure pulse
% create the time array
% define two point sources at (Nx/2, Ny/2, Nz/2) and Nx/2, Ny/2kgrid.t_array), Nz/2+5
source.p_mask = zeros(Nx, Ny, Nz);
source.p_mask(2, 2, 50) = 1;
source_mag = 2500; % [Pa] Max Amplitude
Pulsewidth=7e-12; %Pulse Width [s]
C=source_mag*Pulsewidth^2/(exp(-0.5)); %Constant adjusted so that peaks of pressure are equal to specified Source magnitude
source.p = (-C/Pulsewidth^3)*(kgrid.t_array-PhaseShift).*exp(-((kgrid.t_array-PhaseShift).^2)/(2*Pulsewidth^2)); %Derivation result as listed in .pdf
% filter the source to remove high frequencies not supported by the grid
source.p = filterTimeSeries(kgrid, medium, source.p);
sensor.mask=zeros(size(source.p_mask));
sensor.mask(2,2,50:1050)=1;
% define the field parameters to record
%sensor.record = {'p', 'p_final'};
% input arguments
input_args = {'DisplayMask', source.p_mask, 'DisplayMask',sensor.mask 'DataCast', 'single','PMLInside',false, 'PMLSize', [20, 20, 20],};
% run the simulation
sensor_data = kspaceFirstOrder3D(kgrid, medium, source, sensor, input_args{:});