Hi, I've been writing a script to simulate the generation of pressure from a Photoacoustic point source. I have a source which is of the form of a derivative of a gaussian pulse with 7ps pulse duration However for some reason my sensor is detecting an enormous amplitude. My code is below. Can someone provide some insight into this?
Thanks
% Simulations In Three Dimensions Example
%
% This example provides a simple demonstration of using k-Wave for the
% simulation and detection of a time varying pressure source within a
% three-dimensional heterogeneous propagation medium. It builds on the
% Monopole Point Source In A Homogeneous Propagation Medium Example and
% Simulations In Three Dimensions examples.
%
% author: Bradley Treeby
% date: 20th January 2010
% last update: 25th August 2014
%
% This function is part of the k-Wave Toolbox (http://www.k-wave.org)
% Copyright (C) 2009-2014 Bradley Treeby and Ben Cox
% This file is part of k-Wave. k-Wave is free software: you can
% redistribute it and/or modify it under the terms of the GNU Lesser
% General Public License as published by the Free Software Foundation,
% either version 3 of the License, or (at your option) any later version.
%
% k-Wave is distributed in the hope that it will be useful, but WITHOUT ANY
% WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
% FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for
% more details.
%
% You should have received a copy of the GNU Lesser General Public License
% along with k-Wave. If not, see <http://www.gnu.org/licenses/>.
clear all;
close all;
% =========================================================================
% SIMULATION
% =========================================================================
% create the computational grid
Nx = 64; % number of grid points in the x direction
Ny = 64; % number of grid points in the y direction
Nz = 64; % number of grid points in the z direction
dx = 0.1e-9; % grid point spacing in the x direction [m]
dy = 0.1e-9; % grid point spacing in the y direction [m]
dz = 0.1e-9; % 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 = 1500*ones(Nx, Ny, Nz); % [m/s]
%medium.sound_speed(1:Nx/2, :, :) = 1800; % [m/s]
medium.density = 1000*ones(Nx, Ny, Nz); % [kg/m^3]
%medium.density(:, Ny/4:end, :) = 1200; % [kg/m^3]
medium.alpha_coeff = 0.0022; % [dB/(MHz^y cm)]
medium.alpha_power = 2;
dt=1e-13;
t_end=50e-12;
kgrid.t_array=[0:dt:t_end];
PhaseShift=30e-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_radius = 5; % [grid points]
source.p_mask = zeros(Nx, Ny, Nz);
source.p_mask(Nx/2, Ny/2, Nz/2) = 1;
%source.p_mask(Nx/2, Ny/2, Nz/2+5) = 1;
% define a time varying source
%source_freq = 2e11; % [Hz]
source_mag = 10; % [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
%plot(kgrid.t_array,source.p);
% 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(Nx/2,Ny/2,Nz/2+20)=1;
% define a series of Cartesian points to collect the data
% y = (-20:2:20)*dy; % [m]
% z = (-20:2:20)*dz; % [m]
% x = 20*dx*ones(size(z)); % [m]
% sensor.mask = [x; y; z];
% define the field parameters to record
%sensor.record = {'p', 'p_final'};
figure;
% input arguments
input_args = {'DisplayMask', source.p_mask,'PMLInside',false, 'DataCast', 'single', 'CartInterp', 'nearest'};
% run the simulation
sensor_data = kspaceFirstOrder3D(kgrid, medium, source, sensor, input_args{:});
% =========================================================================
% VISUALISATION
% =========================================================================
% view final pressure field slice by slice
%flyThrough(sensor_data.p_final);
% plot the position of the source and sensor
voxelPlot(double(source.p_mask | sensor.mask));
view(90, 0);
%Plot the input pressure
figure('Name','Pressure vs Time','NumberTitle','off')
plot(kgrid.t_array,source.p);
%Plot FFT of input pressure
figure('Name','f vs amplitude','NumberTitle','off')
[f, source.p_as] = spect(source.p, 1/dt);
plot(f,source.p_as);
xlabel('Frequency (Hz)')
ylabel('Amplitude')
% plot the simulated sensor data
figure('Name','time vs sensor data','NumberTitle','off')
plot(kgrid.t_array,sensor_data);
% colormap(getColorMap);
ylabel('Measured Pressure (Pa)');
xlabel('Time (s)');
% colorbar;
figure('Name','f vs amplitude','NumberTitle','off')
[f, sensor_data_as] = spect(sensor_data, 1/dt);
plot(f,sensor_data_as);
xlabel('Frequency (Hz)')
ylabel('Amplitude')