http://www.k-wave.org/forum/topic/large-amplitudes-when-simulating-a-pressure-rise-from-a-gaussian-laser-pulse Support for the k-Wave MATLAB toolbox en-US Tue, 12 May 2026 23:07:56 +0000 http://bbpress.org/?v=1.0.2 q http://www.k-wave.org/forum/search.php http://www.k-wave.org/forum/topic/large-amplitudes-when-simulating-a-pressure-rise-from-a-gaussian-laser-pulse#post-5903 Wed, 26 Apr 2017 20:07:19 +0000 Bradley Treeby 5903@http://www.k-wave.org/forum/ <p>Glad you got it to work! </p> http://www.k-wave.org/forum/topic/large-amplitudes-when-simulating-a-pressure-rise-from-a-gaussian-laser-pulse#post-5897 Mon, 24 Apr 2017 20:31:09 +0000 samantp 5897@http://www.k-wave.org/forum/ <p>I actually found the problem. The time step was too short and that was causing the pressure to blow up. Sorry for the stupid question! </p> http://www.k-wave.org/forum/topic/large-amplitudes-when-simulating-a-pressure-rise-from-a-gaussian-laser-pulse#post-5894 Mon, 17 Apr 2017 20:17:10 +0000 samantp 5894@http://www.k-wave.org/forum/ <p>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?</p> <p>Thanks</p> <pre><code>% 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 (<a href="http://www.k-wave.org" rel="nofollow">http://www.k-wave.org</a>) % 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 &lt;http://www.gnu.org/licenses/&gt;. 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 = {&#39;p&#39;, &#39;p_final&#39;}; figure; % input arguments input_args = {&#39;DisplayMask&#39;, source.p_mask,&#39;PMLInside&#39;,false, &#39;DataCast&#39;, &#39;single&#39;, &#39;CartInterp&#39;, &#39;nearest&#39;}; % 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(&#39;Name&#39;,&#39;Pressure vs Time&#39;,&#39;NumberTitle&#39;,&#39;off&#39;) plot(kgrid.t_array,source.p); %Plot FFT of input pressure figure(&#39;Name&#39;,&#39;f vs amplitude&#39;,&#39;NumberTitle&#39;,&#39;off&#39;) [f, source.p_as] = spect(source.p, 1/dt); plot(f,source.p_as); xlabel(&#39;Frequency (Hz)&#39;) ylabel(&#39;Amplitude&#39;) % plot the simulated sensor data figure(&#39;Name&#39;,&#39;time vs sensor data&#39;,&#39;NumberTitle&#39;,&#39;off&#39;) plot(kgrid.t_array,sensor_data); % colormap(getColorMap); ylabel(&#39;Measured Pressure (Pa)&#39;); xlabel(&#39;Time (s)&#39;); % colorbar; figure(&#39;Name&#39;,&#39;f vs amplitude&#39;,&#39;NumberTitle&#39;,&#39;off&#39;) [f, sensor_data_as] = spect(sensor_data, 1/dt); plot(f,sensor_data_as); xlabel(&#39;Frequency (Hz)&#39;) ylabel(&#39;Amplitude&#39;)</code></pre>